ChEn-5310: Computational Continuum Transport Phenomena Fall 2021 UMass Lowell; Prof. V. F. de Almeida 16Oct21
$ \newcommand{\Amtrx}{\boldsymbol{\mathsf{A}}} \newcommand{\Bmtrx}{\boldsymbol{\mathsf{B}}} \newcommand{\Mmtrx}{\boldsymbol{\mathsf{M}}} \newcommand{\Imtrx}{\boldsymbol{\mathsf{I}}} \newcommand{\Pmtrx}{\boldsymbol{\mathsf{P}}} \newcommand{\Lmtrx}{\boldsymbol{\mathsf{L}}} \newcommand{\Umtrx}{\boldsymbol{\mathsf{U}}} \newcommand{\Smtrx}{\boldsymbol{\mathsf{S}}} \newcommand{\xvec}{\boldsymbol{\mathsf{x}}} \newcommand{\avec}{\boldsymbol{\mathsf{a}}} \newcommand{\bvec}{\boldsymbol{\mathsf{b}}} \newcommand{\cvec}{\boldsymbol{\mathsf{c}}} \newcommand{\rvec}{\boldsymbol{\mathsf{r}}} \newcommand{\fvec}{\boldsymbol{\mathsf{f}}} \newcommand{\mvec}{\boldsymbol{\mathsf{m}}} \newcommand{\gvec}{\boldsymbol{\mathsf{g}}} \newcommand{\zerovec}{\boldsymbol{\mathsf{0}}} \newcommand{\norm}[1]{\bigl\lVert{#1}\bigr\rVert} \newcommand{\transpose}[1]{{#1}^\top} \DeclareMathOperator{\rank}{rank} \newcommand{\Reals}{\mathbb{R}} \newcommand{\thetavec}{\boldsymbol{\theta}} \newcommand{\Ecal}{\mathcal{E}} $
The following sections describe what is referred to in the literature as the one-dimensional Peclet problem with Dirichlet boundary conditions. This is a classical boundary-value problem of mathematics.
Solve the Peclet model problem. Find $u:[a,b]\subset\Reals\rightarrow\Reals$ such that:
\begin{align*} c(x)\,u'(x) &= -\bigl(-D(x)\, u'\bigr)'(x) + s(x)\,u(x) + f(x) \quad\quad \forall \quad\quad x\in\ ]a,b[, \\ u(a) &= u_a, \\ u(b) &= u_b. \end{align*}This problem is linear and has an analytical solution for given data: convective speed, $c(x)$, diffusion coefficient, $D(x)$, source, $s(x)$ slope, source bias, $f(x)$. The diffusion flux associated to the quantity $u$, is denoted $q := -D(x)\,u'$, and it is often of interest as a derived quantity.
The values of the dependent variable are given on the two end points of the domain. This is called essential boundary conditions or Dirichlet boundary conditions. If the values are equal to zero, the boundary condition is referred to as homogeneous.
Give the scaffolding below, build a manufactured solution by computing the corresponding source function.
'''Domain'''
x_a = 0
x_b = 25
'''Parameters and data'''
conv_speed = 0.1
diff_coeff = 0.1
source_slope_value = -1e-2
u_a_0 = 3.5 # initial value
u_b_0 = 3.5 # initial value
'''Generate the scaffolding of the manufactured solution'''
from engy_5310.toolkit import TargetFunction
import numpy as np
shape_pts_y = [u_a_0,10,-7,-5,9,-8,28,-23,12,1,10,15,18,u_b_0]
shape_pts_x = [x for x in np.linspace(x_a, x_b, len(shape_pts_y))]
shape_pts = list(zip(shape_pts_x, shape_pts_y))
f = TargetFunction(shape_pts, type='linear')
f.plot(n_plot_pts=300, y_label=r'u_{frame}(x)', show_shape_pts=True, title='Target Function')
'''Build the basis functions list'''
import math
wavelength = x_b - x_a
kappa = 2*math.pi/wavelength
N = 10 # number of pairs of sine/cosine
try:
from engy_5310.toolkit import FourierBasis
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
fb = FourierBasis(kappa, N, x_a, x_b)
fb.plot()
'''Build the Gram matrix'''
gram_mtrx = fb.gram_matrix(constrain_end_points=True)
#np.set_printoptions(precision=1, threshold=1000, edgeitems=20, linewidth=200)
#print(gram_mtrx)
if gram_mtrx.shape[0] > gram_mtrx.shape[1]:
print('G is overdetermined.')
elif gram_mtrx.shape[0] < gram_mtrx.shape[1]:
print('G is underdetermined.')
else:
print('G is determined.')
try:
from engy_5310.toolkit import matrix_rank
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
if matrix_rank(gram_mtrx) == min(gram_mtrx.shape):
print('G is full rank.')
else:
print('G is rank deficient.')
G is determined. G is full rank.
'''Build load vector'''
b_vec = np.zeros(len(fb.basis_func_list), dtype=np.float64)
for i, phi_i in enumerate(fb.basis_func_list):
(b_vec[i], _) = fb.inner_product(f, phi_i, epsabs=1e-3, epsrel=1.e-3, limit=1500) \
+ f(x_a)*phi_i(x_a) + f(x_b)*phi_i(x_b)
'''Compute optimal coefficient vector'''
beta_vec = np.linalg.solve(gram_mtrx, b_vec)
'''Build the best approximant function'''
def u_manufactured(x):
return fb.evaluation_matrix(x)@beta_vec
def u_prime_manufactured(x):
return fb.evaluation_matrix(x, derivative=True)@beta_vec
def diff_flux_x_manufactured(x):
return -diff_coeff*u_prime_manufactured(x)
'''Plot comparison of '''
n_pts = 300
x_pts = np.linspace(x_a, x_b, n_pts)
u = u_manufactured(x_pts)
f.plot(g=(x_pts,u), n_plot_pts=n_pts,
y_label=r'$u_\mathrm{mfac}(x)$',
title=r'Manufactured Function $u_\mathrm{mfac}$ ($m\rightarrow\infty$, n='+str(len(fb.basis_func_list))+')',
f_line_style='dashed', g_line_style='continuous',
f_legend=None, g_legend=r'$u_\mathrm{mfac}$')
print('u_mfac(x_a) = ', u_manufactured(x_a), '; u_a_0 = ', u_a_0)
print('u_mfac(x_b) = ', u_manufactured(x_b), '; u_b_0 = ', u_b_0)
u_mfac(x_a) = [3.84494] ; u_a_0 = 3.5 u_mfac(x_b) = [3.84494] ; u_b_0 = 3.5
print("u'_mfac(x_a) = ", u_prime_manufactured(x_a))
print("u'_mfac(x_b) = ", u_prime_manufactured(x_b))
u'_mfac(x_a) = [-2.6005] u'_mfac(x_b) = [-2.6005]
'''Update Boundary Condition Values'''
u_a = u_manufactured(x_a)
u_b = u_manufactured(x_b)
'''Build the best approximant function'''
def u_prime2_manufactured(x):
return fb.evaluation_matrix(x, derivative_2=True)@beta_vec
'''Build the best approximant function'''
def source_bias(x):
return conv_speed * u_prime_manufactured(x) - diff_coeff * u_prime2_manufactured(x) - source_slope_value * u_manufactured(x)
'''Source function'''
n_pts = 500
x = np.linspace(x_a, x_b, n_pts)
f_values = source_bias(x)
s_u_values = source_slope_value * u_manufactured(x)
import matplotlib.pyplot as plt
#%matplotlib inline
#plt.style.use('dark_background')
(fig, ax1) = plt.subplots(1, figsize=(14, 5))
ax1.plot(x, f_values, 'r-', label='Source Bias function')
ax1.set_xlabel(r'$x$ [cm]', fontsize=18)
ax1.set_ylabel(r'$f(x)$', fontsize=18, color='red')
ax1.tick_params(axis='y', labelcolor='red', labelsize=16)
ax1.tick_params(axis='x', labelsize=16)
ax1.legend(loc='upper left', fontsize=12)
ax1.grid(True)
ax2 = ax1.twinx()
ax2.plot(x, s_u_values,'-', color='blue', label='Source linear term')
ax2.set_ylabel(r"$S(x)\, u_\mathrm{mfac}(x)$", fontsize=16, color='blue')
ax2.tick_params(axis='y', labelcolor='blue', labelsize=16)
ax2.legend(loc='upper right', fontsize=12)
ax3 = ax1.twinx()
ax3.plot(x, s_u_values+f_values,'-', color='black', label='Source')
ax3.set_ylabel(r"$S(x)\, u_\mathrm{mfac}(x) + f(x)$", fontsize=16, color='black')
ax3.tick_params(axis='y', labelcolor='black', labelsize=16)
ax3.legend(loc='upper center', fontsize=12)
ax3.spines["right"].set_position(("axes", 1.125))
plt.title(r'Manufactured Source', fontsize=20)
plt.show()
'''Source bias evaluated at the boundaries'''
np.set_printoptions(precision=5)
print('f(x_a) = ',source_bias(x_a), ' f(x_b) = ',source_bias(x_b))
f(x_a) = [-1.54527] f(x_b) = [-1.54527]
The Galerkin weak formulation is as follows. Find $u = u_0 + w \in H^1\!\bigl([a,b]\bigr)$ so that
\begin{align*} \int\limits_a^b c\, u'(x)\, v(x)\,dx - \int\limits_a^b -D\, u'(x)\,v'(x)\,dx - \int\limits_a^b S\,v(x)\,dx &= 0 \quad\quad \forall \quad\quad v \in H^1_0\!\bigl([a,b]\bigr), \end{align*}or \begin{align*} \bigl(c\, u', v\bigr) - \bigl(-D\, u',v'\bigr) - \bigl(S,v\bigr) &= 0 \quad\quad \forall \quad\quad v \in H^1_0\!\bigl([a,b]\bigr), \end{align*}
where $H^1\!\bigl([a,b]\bigr) := \bigl\{ u:[a,b]\subset\Reals\rightarrow \Reals \mid \int_a^b u'^2 + u^2\,dx < \infty\bigr\}$ and $H^1_0\!\bigl([a,b]\bigr) := \bigl\{ v \mid v \in H^1(a,b), v(a) = 0, v(b) =0 \bigr\}$. Both function sets as just defined are Hilbert spaces. The function $v$ is called a test function. Because $v$ and $u$ are sought in very similar sets of functions, this weak form is called Galerkin's weak form.
The coefficients $\cvec_0 := \{c_i\mid i=1,\ldots,N\}$ solve
\begin{equation*} \overset{(N\times N)}{\Amtrx}\,\overset{(N\times 1)}{\cvec_0} = \overset{(N\times 1)}\bvec , \end{equation*}where:
This formulation uses basis functions that satisfy homogeneous boundary conditions, however the linear algebraic problem for the optimum coefficients accounts for the inhomogeneous boundary condition data through the lift function $w$. Therefore the solution is
\begin{equation*} u_N = u_0 + w. \end{equation*}Since $w\in V(a,b)$, construct as follows $w = \sum\limits_{i=1}^N\,\alpha_i\,\phi_i$ such that $w(a) = u_a$ and $w(b) = u_b$.
'''Build the lift basis functions'''
degree = 3
try:
from engy_5310.toolkit import LagrangeFEMBasis
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
lift_basis = LagrangeFEMBasis(degree=degree, n_elem=3, x_partition=[x_a, x_b])
lift_basis.plot(n_plot_pts=300)
'''Build the coefficients of the lift function'''
#alpha_vec = np.zeros(len(lift_basis.basis_func_list))
alpha_vec = np.random.random(len(lift_basis.basis_func_list))
alpha_vec[0] = u_a
if degree == 2:
alpha_vec[-2] = u_b
elif degree == 1:
alpha_vec[-1] = u_b
elif degree == 3:
alpha_vec[-3] = u_b
'''Build the lift function'''
def w_lift(x):
return lift_basis.evaluation_matrix(x)@alpha_vec
'''Build the lift function derivative'''
def w_lift_prime(x):
return lift_basis.evaluation_matrix(x, derivative=True)@alpha_vec
'''Test boundary values'''
print('w(a) = ',w_lift(x_a), ' w(b) = ', w_lift(x_b))
w(a) = [3.84494] w(b) = [3.84494]
'''Rayleigh Ritz Method with Lagrange Basis Functions'''
n_pts = 500
x = np.linspace(x_a, x_b, n_pts)
u_values = w_lift(x)
import matplotlib.pyplot as plt
#%matplotlib inline
#plt.style.use('dark_background')
plt.figure(1, figsize=(14, 5))
plt.plot(x, u_values, 'r-', label='Lift function')
plt.title(r'Rayleigh-Ritz Method with FE Lagrange Basis Functions (n='+str(len(lift_basis.basis_func_list))+')', fontsize=20)
plt.ylabel(r'$w(x)$', fontsize=18)
plt.xlabel(r'$x$', fontsize=18)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.legend(loc='best',fontsize=12)
plt.grid(True)
plt.show()
'''Build the basis functions'''
degree = 3
n_elem=[4,10,4]
degree = [2,3,2]
#n_elem = [3,8,3]
try:
from engy_5310.toolkit import LagrangeFEMBasis
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
femlb = LagrangeFEMBasis(degree=degree, n_elem=n_elem, x_partition=[x_a, x_a+2, x_b-2, x_b], bc_x_min='dirichlet', bc_x_max='dirichlet')
femlb.plot(n_plot_pts=800)
'''Mesh artifact'''
for (p_id, elem_gnode_ids) in enumerate(femlb.local_to_global_node_id_map):
print('Partition %i'%p_id, ' Degree %i'%(femlb.degree_partition[p_id]))
for (e_id, gnode_ids) in enumerate(elem_gnode_ids):
printout = '\t element #%i:'%e_id
for (lnode_id, gnode_id) in enumerate(gnode_ids):
printout += ' x_%i = %+5.2e'%(gnode_id, femlb.patches[p_id][e_id][lnode_id])
print(printout)
Partition 0 Degree 2 element #0: x_-1 = +0.00e+00 x_2 = +5.00e-01 x_1 = +2.50e-01 element #1: x_2 = +5.00e-01 x_4 = +1.00e+00 x_3 = +7.50e-01 element #2: x_4 = +1.00e+00 x_6 = +1.50e+00 x_5 = +1.25e+00 element #3: x_6 = +1.50e+00 x_8 = +2.00e+00 x_7 = +1.75e+00 Partition 1 Degree 3 element #0: x_8 = +2.00e+00 x_11 = +4.10e+00 x_9 = +2.70e+00 x_10 = +3.40e+00 element #1: x_11 = +4.10e+00 x_14 = +6.20e+00 x_12 = +4.80e+00 x_13 = +5.50e+00 element #2: x_14 = +6.20e+00 x_17 = +8.30e+00 x_15 = +6.90e+00 x_16 = +7.60e+00 element #3: x_17 = +8.30e+00 x_20 = +1.04e+01 x_18 = +9.00e+00 x_19 = +9.70e+00 element #4: x_20 = +1.04e+01 x_23 = +1.25e+01 x_21 = +1.11e+01 x_22 = +1.18e+01 element #5: x_23 = +1.25e+01 x_26 = +1.46e+01 x_24 = +1.32e+01 x_25 = +1.39e+01 element #6: x_26 = +1.46e+01 x_29 = +1.67e+01 x_27 = +1.53e+01 x_28 = +1.60e+01 element #7: x_29 = +1.67e+01 x_32 = +1.88e+01 x_30 = +1.74e+01 x_31 = +1.81e+01 element #8: x_32 = +1.88e+01 x_35 = +2.09e+01 x_33 = +1.95e+01 x_34 = +2.02e+01 element #9: x_35 = +2.09e+01 x_38 = +2.30e+01 x_36 = +2.16e+01 x_37 = +2.23e+01 Partition 2 Degree 2 element #0: x_38 = +2.30e+01 x_40 = +2.35e+01 x_39 = +2.32e+01 element #1: x_40 = +2.35e+01 x_42 = +2.40e+01 x_41 = +2.38e+01 element #2: x_42 = +2.40e+01 x_44 = +2.45e+01 x_43 = +2.42e+01 element #3: x_44 = +2.45e+01 x_-1 = +2.50e+01 x_45 = +2.48e+01
Full assembly of $\Amtrx$
'''Assemble A matrix directly (slow)'''
import numpy as np
def get_a_matrix(femlb):
n = len(femlb.basis_func_list)
a_mtrx = np.zeros((n, n), dtype=np.float64)
for i,phi_i in enumerate(femlb.basis_func_list):
for j,phi_prime_j in enumerate(femlb.basis_func_prime_list):
(a_ij, _) = femlb.inner_product(phi_prime_j, phi_i)
a_mtrx[i,j] = a_ij * conv_speed
for i,phi_prime_i in enumerate(femlb.basis_func_prime_list):
for j,phi_prime_j in enumerate(femlb.basis_func_prime_list):
(a_ij, _) = femlb.inner_product(phi_prime_j, phi_prime_i)
a_mtrx[i,j] += a_ij * diff_coeff
for i,phi_i in enumerate(femlb.basis_func_list):
for j,phi_j in enumerate(femlb.basis_func_list):
(a_ij, _) = femlb.inner_product(phi_j, phi_i)
a_mtrx[i,j] -= a_ij * source_slope_value
return a_mtrx
'''Assemble A matrix directly (slow)'''
a1_mtrx = get_a_matrix(femlb)
try:
from engy_5310.toolkit import matrix_rank
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
if matrix_rank(a1_mtrx) == min(a1_mtrx.shape):
print('A is full rank.')
else:
print('A is rank deficient.')
A is full rank.
np.set_printoptions(precision=2, threshold=200, edgeitems=6, linewidth=170)
print(a1_mtrx)
a1_mtrx.shape
try:
from engy_5310.toolkit import plot_matrix
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
plot_matrix(a1_mtrx, plt.get_cmap('copper'))
[[ 0.93 -0.6 0.05 -0.47 0. 0. ... 0. 0. 0. 0. 0. 0. ] [-0.47 1.07 0. 0. 0. 0. ... 0. 0. 0. 0. 0. 0. ] [ 0.08 0. 0.93 -0.6 0.05 -0.47 ... 0. 0. 0. 0. 0. 0. ] [-0.6 0. -0.47 1.07 0. 0. ... 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.08 0. 0.93 -0.6 ... 0. 0. 0. 0. 0. 0. ] [ 0. 0. -0.6 0. -0.47 1.07 ... 0. 0. 0. 0. 0. 0. ] ... [ 0. 0. 0. 0. 0. 0. ... 1.07 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. 0.93 -0.6 0.05 -0.47 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. -0.47 1.07 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. 0.08 0. 0.93 -0.6 -0.47] [ 0. 0. 0. 0. 0. 0. ... 0. -0.6 0. -0.47 1.07 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. 0. 0. -0.6 0. 1.07]] matrix shape = (45, 45)
Element-by-element assembly of $\Amtrx$
'''Assemble A matrix element-by-element'''
from scipy.integrate import quad
def assemble_a_matrix(fembf):
n_dof = fembf.n_dof # total number of degrees of freedom
a_mtrx = np.zeros((n_dof,n_dof), dtype=np.float64)
for (p_id, degree) in enumerate(fembf.degree_partition):
# number of local nodes
n = len(fembf.parent_basis_func_list[degree])
local_a_conv_mtrx = np.zeros((n, n), dtype=np.float64)
local_a_diff_mtrx = np.zeros((n, n), dtype=np.float64)
local_a_source_mtrx = np.zeros((n, n), dtype=np.float64)
for I in range(n):
parent_basis_func_prime_I = fembf.parent_basis_func_prime_list[degree][I]
parent_basis_func_I = fembf.parent_basis_func_list[degree][I]
for J in range(n):
integrand = lambda zeta: conv_speed * fembf.parent_basis_func_prime_list[degree][J](zeta) * \
parent_basis_func_I(zeta)
(local_a_conv_mtrx[I,J], _) = quad(integrand, -1, 1)
integrand = lambda zeta: diff_coeff * fembf.parent_basis_func_prime_list[degree][J](zeta) * \
parent_basis_func_prime_I(zeta)
(local_a_diff_mtrx[I,J], _) = quad(integrand, -1, 1)
integrand = lambda zeta: - source_slope_value * fembf.parent_basis_func_list[degree][J](zeta) * \
parent_basis_func_I(zeta)
(local_a_source_mtrx[I,J], _) = quad(integrand, -1, 1)
for (e, dof_ids) in enumerate(fembf.dof_ids[p_id]):
patch_nodes_x = fembf.patches[p_id][e]
h_e = patch_nodes_x[1] - patch_nodes_x[0]
parent_mapping_jacobian = fembf.parent_mapping_prime(h_e)
for (I, i) in enumerate(dof_ids):
if i < 0: continue
for (J, j) in enumerate(dof_ids):
if j < 0: continue
a_mtrx[i,j] += local_a_conv_mtrx[I,J] + \
local_a_diff_mtrx[I,J] / parent_mapping_jacobian + \
local_a_source_mtrx[I,J] * parent_mapping_jacobian
return a_mtrx
'''Assemble A matrix element-by-element'''
a2_mtrx = assemble_a_matrix(femlb)
try:
from engy_5310.toolkit import matrix_rank
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
if matrix_rank(a2_mtrx) == min(a2_mtrx.shape):
print('A is full rank.')
else:
print('A is rank deficient.')
A is full rank.
print(a2_mtrx)
a2_mtrx.shape
try:
from engy_5310.toolkit import plot_matrix
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
plot_matrix(a2_mtrx, plt.get_cmap('copper'))
[[ 0.93 -0.6 0.05 -0.47 0. 0. ... 0. 0. 0. 0. 0. 0. ] [-0.47 1.07 0. 0. 0. 0. ... 0. 0. 0. 0. 0. 0. ] [ 0.08 0. 0.93 -0.6 0.05 -0.47 ... 0. 0. 0. 0. 0. 0. ] [-0.6 0. -0.47 1.07 0. 0. ... 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0.08 0. 0.93 -0.6 ... 0. 0. 0. 0. 0. 0. ] [ 0. 0. -0.6 0. -0.47 1.07 ... 0. 0. 0. 0. 0. 0. ] ... [ 0. 0. 0. 0. 0. 0. ... 1.07 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. 0.93 -0.6 0.05 -0.47 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. -0.47 1.07 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. 0.08 0. 0.93 -0.6 -0.47] [ 0. 0. 0. 0. 0. 0. ... 0. -0.6 0. -0.47 1.07 0. ] [ 0. 0. 0. 0. 0. 0. ... 0. 0. 0. -0.6 0. 1.07]] matrix shape = (45, 45)
print('||A1-A2|| = ',np.linalg.norm(a1_mtrx-a2_mtrx))
||A1-A2|| = 7.651892437423125e-15
Full assembly of $\bvec$
'''Assemble b vector directly (slow)'''
def get_b_vector(femlb):
b_vec = np.zeros(len(femlb.basis_func_list), dtype=np.float64)
for i, phi_i in enumerate(femlb.basis_func_list):
(b_vec[i], _) = femlb.inner_product(source_bias, phi_i)
c_x_w_lift_prime = lambda x: conv_speed * w_lift_prime(x)
(term1, _) = femlb.inner_product(c_x_w_lift_prime, phi_i)
b_vec[i] -= term1
d_x_w_lift_prime = lambda x: diff_coeff * w_lift_prime(x)
(term2, _) = femlb.inner_product(d_x_w_lift_prime, femlb.basis_func_prime_list[i])
b_vec[i] -= term2
s_x_w = lambda x: source_slope_value * w_lift(x)
(term3, _) = femlb.inner_product(s_x_w, phi_i)
b_vec[i] += term3
return b_vec
b1_vec = get_b_vector(femlb)
print(b1_vec)
plot_matrix(b1_vec, color_map=plt.get_cmap('copper'))
[-0.06 -0.29 0.17 0.12 0.3 0.49 0.29 0.53 -0.38 -0.23 -1.35 0.02 -0.19 -0.12 0.17 0.73 1.3 0.64 -1.62 -1.53 -1.45 4.06 1.4 1.36 -5.14 -0.51 -0.57 3.24 -0.29 0.32 -0.95 0.53 0.22 0.51 0.36 0.26 0.25 0.52 0.05 0.19 -0.11 -0.06 -0.27 -0.39 -0.59] matrix shape = (45, 1)
Element-by-element assembly of $\bvec$
'''Assemble b vector element by element'''
from scipy.integrate import quad
def assemble_b_vector(fembf):
n_dof = fembf.n_dof # total number of degrees of freedom
b_vec = np.zeros(n_dof, dtype=np.float64)
for (p_id, degree) in enumerate(fembf.degree_partition):
# number of local nodes
n = len(fembf.parent_basis_func_list[degree])
local_b_bias_vec = np.zeros(n, dtype=np.float64)
local_b_conv_vec = np.zeros(n, dtype=np.float64)
local_b_diff_vec = np.zeros(n, dtype=np.float64)
local_b_source_vec = np.zeros(n, dtype=np.float64)
for (e, dof_ids) in enumerate(fembf.dof_ids[p_id]):
patch_nodes_x = fembf.patches[p_id][e]
h_e = patch_nodes_x[1] - patch_nodes_x[0]
x_e_bar = (patch_nodes_x[0] + patch_nodes_x[1])/2.0
parent_mapping_jacobian = fembf.parent_mapping_prime(h_e)
source_bias_parent = lambda zeta, x_e_bar=x_e_bar, h_e=h_e: \
source_bias(fembf.parent_mapping(zeta, x_e_bar, h_e))
w_lift_prime_parent = lambda zeta, x_e_bar=x_e_bar, h_e=h_e: \
w_lift_prime(fembf.parent_mapping(zeta, x_e_bar, h_e)) * h_e/2.0
w_lift_parent = lambda zeta, x_e_bar=x_e_bar, h_e=h_e: \
w_lift(fembf.parent_mapping(zeta, x_e_bar, h_e))
for I in range(n):
parent_basis_func_I = fembf.parent_basis_func_list[degree][I]
parent_basis_func_prime_I = fembf.parent_basis_func_prime_list[degree][I]
integrand = lambda zeta: source_bias_parent(zeta) * parent_basis_func_I(zeta)
(local_b_bias_vec[I], _) = quad(integrand, -1, 1)
integrand = lambda zeta: - conv_speed * w_lift_prime_parent(zeta) * parent_basis_func_I(zeta)
(local_b_conv_vec[I], _) = quad(integrand, -1, 1)
integrand = lambda zeta: - diff_coeff * w_lift_prime_parent(zeta) * parent_basis_func_prime_I(zeta)
(local_b_diff_vec[I], _) = quad(integrand, -1, 1)
integrand = lambda zeta: source_slope_value * w_lift_parent(zeta) * parent_basis_func_I(zeta)
(local_b_source_vec[I], _) = quad(integrand, -1, 1)
for (I, i) in enumerate(dof_ids):
if i < 0: continue
b_vec[i] += local_b_bias_vec[I] * parent_mapping_jacobian + \
local_b_conv_vec[I] + \
local_b_diff_vec[I] / parent_mapping_jacobian + \
local_b_source_vec[I] * parent_mapping_jacobian
return b_vec
b2_vec = assemble_b_vector(femlb)
print(b2_vec)
plot_matrix(b2_vec, color_map=plt.get_cmap('copper'))
[-0.05831 -0.28887 0.1731 0.11626 0.29566 0.48974 0.27979 0.52237 -0.39621 -0.24096 -1.37212 0.0061 -0.20593 -0.143 0.13639 0.71391 1.28863 0.651 -1.57949 -1.5009 -1.45943 4.06766 1.39698 1.35002 -5.15114 -0.52665 -0.55869 3.22209 -0.29379 0.32012 -0.91789 0.52827 0.2042 0.49767 0.33162 0.24681 0.21638 0.49727 0.04896 0.18329 -0.11258 -0.05686 -0.2693 -0.38944 -0.58126] matrix shape = (45, 1)
print('||b1-b2|| = ',np.linalg.norm(b1_vec-b2_vec))
||b1-b2|| = 3.0664784075787275e-08
'''Compute optimal coefficient vector'''
#c_star_vec = np.linalg.solve(a1_mtrx, b1_vec)
c_star_vec = np.linalg.solve(a2_mtrx, b2_vec)
'''Build the best approximation function in V_N'''
def u_star(x, femlb, c_star_vec):
u_0 = femlb.evaluation_matrix(x)@c_star_vec
w = lift_basis.evaluation_matrix(x)@alpha_vec
return u_0 + w
'''Build the best approximation function derivative in V_N'''
def u_prime_star(x, femlb, c_star_vec):
u_0 = femlb.evaluation_matrix(x, derivative=True)@c_star_vec
w = lift_basis.evaluation_matrix(x, derivative=True)@alpha_vec
return u_0 + w
def diff_flux_x_star(x, femlb, c_star_vec):
return -diff_coeff*u_prime_star(x, femlb, c_star_vec)
'''Rayleigh Ritz Method with FEM Lagrange Basis Functions'''
def plot_solution(cases=[(femlb, c_star_vec)], bf_type='lagrange', plot_flux=True):
import matplotlib.pyplot as plt
#%matplotlib inline
#plt.style.use('dark_background')
(fig, ax1) = plt.subplots(1, figsize=(14, 5))
if plot_flux:
ax2 = ax1.twinx()
n_pts = 500
x = np.linspace(x_a, x_b, n_pts)
ax1.plot(x, u_manufactured(x), 'k--', label=r'$u$: exact')
if plot_flux:
ax2.plot(x, diff_flux_x_manufactured(x), 'k-.', label=r"$-k\,u'$: exact")
for (case, (femlb, c_star_vec)) in enumerate(cases):
u_values = u_star(x, femlb, c_star_vec)
diff_flux_x_values = diff_flux_x_star(x, femlb, c_star_vec)
ax1.plot(x, u_values, 'r-', label=r'$u_{N_0}$: N = '+str(len(femlb.basis_func_list)))
ax1.set_xlabel(r'$x$ [cm]', fontsize=18)
ax1.set_ylabel(r'$u_{N}$ [g/cc]', color='red', fontsize=18)
ax1.tick_params(axis='y', labelcolor='red', labelsize=16)
ax1.tick_params(axis='x', labelsize=16)
ax1.legend(loc='upper left', fontsize=12)
ax1.grid(True)
if plot_flux:
ax2.plot(x, diff_flux_x_values,'-', color='blue', label=r'$q_x$: N = '+str(len(femlb.basis_func_list)))
ax2.set_ylabel(r"$q_x$ [g/cm$^2$-s]", fontsize=16, color='blue')
ax2.tick_params(axis='y', labelcolor='blue', labelsize=16)
ax2.legend(loc='upper right', fontsize=12)
if bf_type == 'lagrange':
plt.title(r'Galerkin Method with FE Lagrange Basis Functions (degree='+str(degree)+' dim='+str(len(femlb.basis_func_list))+')', fontsize=20)
else:
plt.title(r'Galerkin Method with FE Hermite Basis Functions (degree='+str(degree)+' dim='+str(len(femlb.basis_func_list))+')', fontsize=20)
plt.show()
plot_solution([(femlb, c_star_vec)], bf_type='lagrange', plot_flux=True)
'''L2 norm relative error'''
from scipy.integrate import quad
import math
u_mfac2_integral, error = quad(lambda x: u_manufactured(x)**2, x_a, x_b, limit=1000)
print('||u_mfac||_L2 =', math.sqrt(u_mfac2_integral))
#print('quadrature error =', error)
print('')
integrand = lambda x, bf=femlb, c=c_star_vec: (u_manufactured(x)-u_star(x, bf, c))**2
(u_mfac_minus_u_star_2_integral, _) = quad(integrand, x_a, x_b, limit=1000)
print('||u_mfac-u||_L2 =', math.sqrt(u_mfac_minus_u_star_2_integral))
#print('quadrature error =', error)
print('')
print('Relative error [%%] = %3.3e'%(math.sqrt(u_mfac_minus_u_star_2_integral)/math.sqrt(u_mfac2_integral)*100))
||u_mfac||_L2 = 51.19407408360285 ||u_mfac-u||_L2 = 1.4852399726924892 Relative error [%] = 2.901e+00
'''Nodal errors 2-norm'''
error_2norm = np.linalg.norm(u_manufactured(femlb.gnodes_x) - u_star(femlb.gnodes_x, femlb, c_star_vec))
u_mfac_2norm = np.linalg.norm(u_manufactured(femlb.gnodes_x))
print('||u-u_mfac||_2/||u_mfac||_2 [%%] = %3.3e'%(error_2norm/u_mfac_2norm*100))
||u-u_mfac||_2/||u_mfac||_2 [%] = 1.205e+00
'''Rayleigh-Ritz solution evaluated at the boundaries'''
np.set_printoptions(precision=5)
print('u_N(a) = ',u_star(x_a, femlb, c_star_vec), ' u_N(b) = ',u_star(x_b, femlb, c_star_vec))
u_N(a) = [3.84494] u_N(b) = [3.84494]
'''L2 norm relative error'''
from scipy.integrate import quad
import math
u_prime_mfac2_integral, error = quad(lambda x: u_prime_manufactured(x)**2, x_a, x_b, limit=1000)
print("||u'_mfac||_L2 =", math.sqrt(u_prime_mfac2_integral))
#print('quadrature error =', error)
print('')
integrand = lambda x, bf=femlb, c=c_star_vec: (u_prime_manufactured(x)-u_prime_star(x, bf, c))**2
(u_prime_mfac_minus_u_prime_star_2_integral, _) = quad(integrand, x_a, x_b, limit=1000)
print("||u'_mfac-u'||_L2 =", math.sqrt(u_prime_mfac_minus_u_prime_star_2_integral))
#print('quadrature error =', error)
print('')
print('Relative error [%%] = %3.3e'%(math.sqrt(u_prime_mfac_minus_u_prime_star_2_integral)/math.sqrt(u_prime_mfac2_integral)*100))
||u'_mfac||_L2 = 52.57069490716376 ||u'_mfac-u'||_L2 = 6.829303983137099 Relative error [%] = 1.299e+01
'''Rayleigh-Ritz solution evaluated at the boundaries'''
np.set_printoptions(precision=5)
print("u'_N(a) = ", u_prime_star(x_a, femlb, c_star_vec), " u'_N(b) = ", u_prime_star(x_b, femlb, c_star_vec))
print("u'_mfac(a) = ", u_prime_manufactured(x_a), " u'_mfac(b) = ", u_prime_manufactured(x_b))
u'_N(a) = [-2.5043] u'_N(b) = [-2.8115] u'_mfac(a) = [-2.6005] u'_mfac(b) = [-2.6005]
#'''Hermite finite element basis functions'''
try:
from engy_5310.toolkit import HermiteFEMBasis
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
degree = 3
femhb = HermiteFEMBasis(degree=degree, n_elem=[8,8], x_partition=[x_a, 12, x_b], bc_x_min='dirichlet', bc_x_max='dirichlet')
a_mtrx = get_a_matrix(femhb)
b_vec = get_b_vector(femhb)
c_star_vec = np.linalg.solve(a_mtrx, b_vec)
plot_solution([(femhb, c_star_vec)], bf_type='hermite')
'''L2 norm relative error'''
from scipy.integrate import quad
import math
u_mfac2_integral, error = quad(lambda x: u_manufactured(x)**2, x_a, x_b, limit=1000)
print('||u_mfac||_L2 =', math.sqrt(u_mfac2_integral))
#print('quadrature error =', error)
print('')
integrand = lambda x, bf=femhb, c=c_star_vec: (u_manufactured(x)-u_star(x, bf, c))**2
(u_mfac_minus_u_star_2_integral, _) = quad(integrand, x_a, x_b, epsabs=1e-6, epsrel=1e-6, limit=1000)
print('||u_mfac-u||_L2 =', math.sqrt(u_mfac_minus_u_star_2_integral))
#print('quadrature error =', error)
print('')
print('Relative error [%%] = %3.3e'%(math.sqrt(u_mfac_minus_u_star_2_integral)/math.sqrt(u_mfac2_integral)*100))
||u_mfac||_L2 = 51.19407408360285 ||u_mfac-u||_L2 = 0.4209340958526319 Relative error [%] = 8.222e-01
'''L2 norm relative error'''
from scipy.integrate import quad
import math
u_prime_mfac2_integral, error = quad(lambda x: u_prime_manufactured(x)**2, x_a, x_b, limit=1000)
print("||u'_mfac||_L2 =", math.sqrt(u_prime_mfac2_integral))
#print('quadrature error =', error)
print('')
integrand = lambda x, bf=femhb, c=c_star_vec: (u_prime_manufactured(x)-u_prime_star(x, bf, c))**2
(u_prime_mfac_minus_u_prime_star_2_integral, _) = quad(integrand, x_a, x_b, limit=1000)
print("||u'_mfac-u'||_L2 =", math.sqrt(u_prime_mfac_minus_u_prime_star_2_integral))
#print('quadrature error =', error)
print('')
print('Relative error [%%] = %3.3e'%(math.sqrt(u_prime_mfac_minus_u_prime_star_2_integral)/math.sqrt(u_prime_mfac2_integral)*100))
||u'_mfac||_L2 = 52.57069490716376 ||u'_mfac-u'||_L2 = 2.4576785296875676 Relative error [%] = 4.675e+00
'''Mesh convergence test'''
import math
def mesh_convergence(type, degree, n_elem, x_partition, fembf, c_star_vec):
assert type in ['lagrange', 'hermite']
fembf_0 = fembf
c_star_vec_0 = np.copy(c_star_vec)
u_0 = u_star(fembf_0.gnodes_x, fembf_0, c_star_vec_0)
cases = [(fembf_0, c_star_vec_0)] # mesh refinement cases
degree = degree
h = list()
for elem_patches in femlb.patches:
for patch in elem_patches:
h.append(patch[1]-patch[0])
h = np.array(h).mean()
integrand = lambda x, bf=fembf, c=c_star_vec: (u_manufactured(x) - u_star(x, bf, c))**2
(u_exact_minus_u_2_integral, _) = quad(integrand, x_a, x_b, limit=1000)
error_0 = math.sqrt(u_exact_minus_u_2_integral)
power_law_data = list()
for i in range(1,15):
n_elem = [j+2 for j in n_elem]
print(n_elem)
if type == 'lagrange':
fembf = LagrangeFEMBasis(degree=degree, n_elem=n_elem, x_partition=x_partition, bc_x_min='dirichlet', bc_x_max='dirichlet')
if type == 'hermite':
fembf = HermiteFEMBasis(degree=degree, n_elem=n_elem, x_partition=x_partition, bc_x_min='dirichlet', bc_x_max='dirichlet')
a_mtrx = assemble_a_matrix(fembf)
b_vec = assemble_b_vector(fembf)
c_star_vec = np.linalg.solve(a_mtrx, b_vec)
h = list()
for elem_patches in femlb.patches:
for patch in elem_patches:
h.append(patch[1]-patch[0])
h = np.array(h).mean()
integrand = lambda x, bf=fembf, c=c_star_vec: (u_manufactured(x) - u_star(x, bf, c))**2
(u_exact_minus_u_2_integral, _) = quad(integrand, x_a, x_b, limit=1000)
error = math.sqrt(u_exact_minus_u_2_integral)
(u_mfac2_integral, _) = quad(u_manufactured, x_a, x_b, limit=1000)
relative_error = math.sqrt(error)/math.sqrt(u_mfac2_integral)*100
if h <= (x_b-x_a)/20: # discard coarse meshes case
power_law_data.append((math.log(h), math.log(error)))
p = (math.log(error) - math.log(error_0))/(math.log(h)-math.log(h_0))
print('# elem = %i; relative ||u_mfac - u_%i||_2 [%%] = %3.3e; h [cm] = %3.2f; p [] = %3.2f'%(sum(n_elem),i,relative_error,h,p))
cases.append((fembf, c_star_vec))
error_0 = error
h_0 = h
if relative_error < 1:
break
print('done.')
return (cases, power_law_data)
'''Build the basis functions'''
degree = 2
n_elem = [2]
x_partition = [x_a, x_b]
try:
from engy_5310.toolkit import LagrangeFEMBasis
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
femlb = LagrangeFEMBasis(degree=degree, n_elem=n_elem, x_partition=x_partition, bc_x_min='dirichlet', bc_x_max='dirichlet')
femlb.plot(n_plot_pts=800)
a_mtrx = assemble_a_matrix(femlb)
b_vec = assemble_b_vector(femlb)
c_star_vec = np.linalg.solve(a_mtrx, b_vec)
(cases, power_law_data) = mesh_convergence('lagrange', degree, n_elem, x_partition, femlb, c_star_vec)
plot_solution(cases)
plot_solution([cases[-1]])
f = TargetFunction(list(reversed(power_law_data)), type='stepwise')
#f.plot(show_shape_pts=True)
try:
from engy_5310.toolkit import LegendreBasis
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
x_tilde_vec = np.array([x[0] for x in reversed(power_law_data)])
lb = LegendreBasis(degree=1, x_min=x_tilde_vec[0], x_max=x_tilde_vec[-1])
a_tilde_mtrx = lb.evaluation_matrix(x_tilde_vec)
f_tilde_vec = np.array([x[1] for x in reversed(power_law_data)])
c_star_vec = np.linalg.solve(a_tilde_mtrx.transpose()@a_tilde_mtrx, a_tilde_mtrx.transpose()@f_tilde_vec)
def g_best_vec_func(x_pts):
return lb.evaluation_matrix(x_pts)@c_star_vec
'''Plot comparison of f and g_best_vec'''
n_pts = 300
x_pts = np.linspace(x_tilde_vec[0], x_tilde_vec[-1], n_pts)
g_best_vec_2 = g_best_vec_func(x_pts)
f.plot(g=(x_pts,g_best_vec_2), n_plot_pts=n_pts, g_line_style='plain', show_shape_pts=True, f_line_style='dashed',
title=r'Compare $f$ and $g$ Functions ($m\rightarrow\infty$, n='+str(len(lb.basis_func_list))+')',
y_label=r'$\ln$ Error', x_label=r'$\ln h$')
print('Power exponent of error estimate = %2.3f'%c_star_vec[1])
power_law_data