ChEn-5310: Computational Continuum Transport Phenomena Fall 2021 UMass Lowell; Prof. V. F. de Almeida 29Sep21
$ \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 Poisson problem with Dirichlet boundary conditions. This is a classical boundary-value problem of mathematics.
Solve the Poisson model problem. Find $u:[a,b]\subset\Reals\rightarrow\Reals$ such that:
\begin{align*} -\bigl(-D(x)\, u'\bigr)'(x) + S(x)\,u(x) + f(x) &= 0 \quad\quad \forall \quad\quad x\in\ ]a,b[, \\ u(a) &= u_a, \\ q_n(b) &= h\,\bigl(u(b)-u_\text{ref}\bigr). \end{align*}This problem is linear and has an analytical solution for given data: 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 normal diffusive flux at $x=b$ is $q_n(b) = -D\,u'(b) = h\bigl(u(b)-u_\text{ref}\bigr)$
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.
Find $u^*_N \in V_N(a,b)\subset V(a,b) = \bigl\{ u:[a,b]\subset\Reals\rightarrow\Reals \bigr\}$ such that it minimizes the constrained Poisson energy norm:
\begin{equation*} \min\limits_{u_0\, \in\, V} \norm{u_0+w-u_N}^2_{\Ecal}, \end{equation*}where $V_N(a,b) := \bigl\{ u_N = \sum\limits_{i=1}^N\,c_i\,\phi_i \mid u_N(a) = 0, (u'_N,u'_N) < \infty \bigr\}$, $\{\phi_i\mid i=1\ldots,N\}$ is a basis of $V_N(a,b)$, and any $w \in V$ with $w(a) = u_a$ is called the lift function.
The optimum coefficients $\cvec^* := \{c_i\mid i=1,\ldots,N\}$ solve
\begin{equation*} \overset{(N\times N)}{\Amtrx}\,\overset{(N\times 1)}{\cvec^*} = \overset{(N\times 1)}\bvec , \end{equation*}where:
This formulation uses basis functions that satisfy the left homogeneous boundary condition, however the linear algebraic problem for the optimum coefficients accounts for the inhomogeneous boundary condition data through the lift function $w$.
'''Domain'''
x_a = 0
x_b = 25
'''Parameters and data'''
diff_coeff = 0.1
source_bias_value = 1e-2
source_slope_value = 1.87e-2
u_a = 2
h = 1
u_ref = 4.8125
'''Generate the source bias'''
try:
from engy_5310.toolkit import TargetFunction
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
shape_pts = [(x_a,source_bias_value), (x_b,source_bias_value)]
#shape_pts = [(x_a,source_bias_value), ((x_a+x_b)/2,1.3*source_bias_value), (x_b,source_bias_value)]
shape_pts = [(x_a,source_bias_value), ((x_b+x_a)/2,-2*source_bias_value/4), (x_b,source_bias_value)]
#shape_pts = [(x_a,source_bias_value), ((x_b+x_a)/8,2*source_bias_value), ((x_b+x_a)/4,source_bias_value), (3*(x_b+x_a)/4,2*source_bias_value), (x_b,source_bias_value)]
print('# shape pts = ', len(shape_pts))
source_bias = TargetFunction(shape_pts, type='stepwise')
source_bias = TargetFunction(shape_pts, type='linear')
#source_bias = TargetFunction(shape_pts, type='legendre')
source_bias.plot(n_plot_pts=200, show_shape_pts=True, title='Source Bias')
# shape pts = 3
'''Generate the source slope'''
try:
from engy_5310.toolkit import TargetFunction
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
shape_pts = [(x_a,source_slope_value), (x_b,source_slope_value)]
#shape_pts = [(x_a,source_slope_value), ((x_b+x_a)/2,source_slope_value/4), (x_b,source_slope_value)]
shape_pts = [(x_a,source_slope_value), ((x_b+x_a)/8,4*source_slope_value), ((x_b+x_a)/4,source_slope_value/3), (3*(x_b+x_a)/4,2*source_slope_value), (x_b,source_slope_value)]
print('# shape pts = ', len(shape_pts))
source_slope = TargetFunction(shape_pts, type='legendre')
#source_slope = TargetFunction(shape_pts, type='linear')
source_slope.plot(n_plot_pts=200, show_shape_pts=True, title='Source Slope', y_label=r'$S(x)$')
# shape pts = 5
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$.
'''Build the basis functions'''
degree = 1
n_elem = 1
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=n_elem, x_partition=[x_a, x_b])
lift_basis.plot(n_plot_pts=300)
'''Build the coefficients of the lift function'''
import numpy as np
alpha_vec = np.zeros(len(lift_basis.basis_func_list))
alpha_vec[0] = u_a
alpha_vec[-1] = u_a/2 # any value
'''Build the lift function'''
def w_lift(x):
a_mtrx = lift_basis.evaluation_matrix(x)
return a_mtrx@alpha_vec
'''Build the lift function derivative'''
def w_lift_prime(x):
a_mtrx = lift_basis.evaluation_matrix(x, derivative=True)
return a_mtrx@alpha_vec
'''Test boundary values'''
print('w(a) = ',w_lift(x_a), ' w(b) = ',w_lift(x_b))
w(a) = [2.] w(b) = [1.]
'''Rayleigh Ritz Method with FEM 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=(18, 5))
plt.plot(x, u_values, 'r-', label='Lift function')
plt.title(r'Rayleigh-Ritz Method with Finite Element 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 = [2,2]
n_elem = [10,2]
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_b-1.5, x_b], bc_x_min='dirichlet')
femlb.plot(n_plot_pts=300, n_legend_cols=16)
'''Build the matrix of coefficients of the linear system (using global basis functions)'''
# Too slow (keep for information)
import numpy as np
n = len(femlb.basis_func_list)
a_mtrx = np.zeros((n, n), dtype=np.float64)
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):
s_x_phi_j = lambda x: source_slope(x) * phi_j(x)
(a_ij, _) = femlb.inner_product(s_x_phi_j, phi_i)
a_mtrx[i,j] -= a_ij
for i,phi_i in enumerate(femlb.basis_func_list):
for j,phi_j in enumerate(femlb.basis_func_list):
a_mtrx[i,j] += h * phi_j(x_b) * phi_i(x_b)
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(a_mtrx) == min(a_mtrx.shape):
print('A is full rank.')
else:
print('A is rank deficient.')
A is full rank.
'''Build load vector (using global basis functions)'''
# Too slow (keep for information)
b_vec = np.zeros(n, dtype=np.float64)
for i, phi_i in enumerate(femlb.basis_func_list):
(b_vec[i], _) = femlb.inner_product(source_bias, phi_i)
(term1, _) = femlb.inner_product(w_lift_prime, femlb.basis_func_prime_list[i])
b_vec[i] -= diff_coeff * term1
s_x_w = lambda x: source_slope(x) * w_lift(x)
(term2, _) = femlb.inner_product(s_x_w, phi_i)
b_vec[i] += term2
b_vec[i] -= h * (w_lift(x_b) - u_ref) * phi_i(x_b)
'''Compute optimal coefficient vector'''
c_star_vec = np.linalg.solve(a_mtrx, b_vec)
'''Build the best approximation function in V_N'''
def u_star(x):
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 in V_N'''
def u_prime_star(x):
u_prime_0 = femlb.evaluation_matrix(x, derivative=True)@c_star_vec
w_prime = lift_basis.evaluation_matrix(x, derivative=True)@alpha_vec
return u_prime_0 + w_prime
'''Rayleigh Ritz Method with Fourier Basis Functions'''
n_pts = 500
x = np.linspace(x_a, x_b, n_pts)
u_values = u_star(x)
import matplotlib.pyplot as plt
#%matplotlib inline
#plt.style.use('dark_background')
plt.figure(1, figsize=(18, 5))
plt.plot(x, u_values, 'r-', label='Rayleigh-Ritz Solution w/ Dirichlet/Robin BC')
plt.title(r'Rayleigh-Ritz Method with Finite Element Lagrange Basis Functions (n='+str(len(femlb.basis_func_list))+')', fontsize=20)
plt.ylabel(r'$u_N(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()
'''Rayleigh-Ritz solution evaluated at the boundaries'''
np.set_printoptions(precision=5)
print('u^*_N(a) = ',u_star(x_a), 'u^*_N(b) = ', u_star(x_b), " u^*'_N(b) = ",u_prime_star(x_b))
u^*_N(a) = [2.] u^*_N(b) = [4.29112] u^*'_N(b) = [5.23683]
'''Flux at boundary'''
print('q_nb computed = ', -diff_coeff*u_prime_star(x_b))
print('h (u(b)-u_ref) = ',h*(u_star(x_b)-u_ref))
print('flux error [%]= ',(-diff_coeff*u_prime_star(x_b) - h*(u_star(x_b)-u_ref))/(h*(u_star(x_b)-u_ref))*100)
q_nb computed = [-0.52368] h (u(b)-u_ref) = [-0.52138] flux error [%]= [0.44102]
The following sections describe what is referred to in the literature as the one-dimensional Poisson problem with Dirichlet boundary conditions. This is a classical boundary-value problem of mathematics.
Solve the Poisson model problem. Find $u:[a,b]\subset\Reals\rightarrow\Reals$ such that:
\begin{align*} -\bigl(-D(x)\, u'\bigr)'(x) + S(x)\,u(x) + f(x) &= 0 \quad\quad \forall \quad\quad x\in\ ]a,b[, \\ q_n(a) &= h\,\bigl(u(a)-u_\text{ref}\bigr) \\ u(b) &= u_b. \end{align*}This problem is linear and has an analytical solution for given data: 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 normal diffusive flux at $x=a$ is $q_n(a) = -\bigl(-D\,u'(a)\bigr) = h\bigl(u(a)-u_\text{ref}\bigr)$.
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.
Find $u^*_N \in V_N(a,b)\subset V(a,b) = \bigl\{ u:[a,b]\subset\Reals\rightarrow\Reals \bigr\}$ such that it minimizes the constrained Poisson energy norm:
\begin{equation*} \min\limits_{u_0\, \in\, V} \norm{u_0+w-u_N}^2_{\Ecal}, \end{equation*}where $V_N(a,b) := \bigl\{ u_N = \sum\limits_{i=1}^N\,c_i\,\phi_i \mid u_N(a) = 0, (u'_N,u'_N) < \infty \bigr\}$, $\{\phi_i\mid i=1\ldots,N\}$ is a basis of $V_N(a,b)$, and any $w \in V$ with $w(a) = u_a$ is called the lift function.
The optimum coefficients $\cvec^* := \{c_i\mid i=1,\ldots,N\}$ solve
\begin{equation*} \overset{(N\times N)}{\Amtrx}\,\overset{(N\times 1)}{\cvec^*} = \overset{(N\times 1)}\bvec , \end{equation*}where:
This formulation uses basis functions that satisfy the left homogeneous boundary condition, however the linear algebraic problem for the optimum coefficients accounts for the inhomogeneous boundary condition data through the lift function $w$.
'''Domain'''
x_a = 0
x_b = 25
'''Parameters and data'''
diff_coeff = 0.1
source_bias_value = 1e-2
source_slope_value = 1.87e-2
u_b = 2
h = 1
u_ref = 4.8125
'''Generate the source bias'''
try:
from engy_5310.toolkit import TargetFunction
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
shape_pts = [(x_a,source_bias_value), (x_b,source_bias_value)]
#shape_pts = [(x_a,source_bias_value), ((x_a+x_b)/2,1.3*source_bias_value), (x_b,source_bias_value)]
shape_pts = [(x_a,source_bias_value), ((x_b+x_a)/2,-2*source_bias_value/4), (x_b,source_bias_value)]
#shape_pts = [(x_a,source_bias_value), ((x_b+x_a)/8,2*source_bias_value), ((x_b+x_a)/4,source_bias_value), (3*(x_b+x_a)/4,2*source_bias_value), (x_b,source_bias_value)]
print('# shape pts = ', len(shape_pts))
source_bias = TargetFunction(shape_pts, type='stepwise')
source_bias = TargetFunction(shape_pts, type='linear')
#source_bias = TargetFunction(shape_pts, type='legendre')
source_bias.plot(n_plot_pts=200, show_shape_pts=True, title='Source Bias')
# shape pts = 3
'''Generate the source slope'''
try:
from engy_5310.toolkit import TargetFunction
except ModuleNotFoundError:
assert False, 'You need to provide your own code here. Bailing out.'
shape_pts = [(x_a,source_slope_value), (x_b,source_slope_value)]
#shape_pts = [(x_a,source_slope_value), ((x_b+x_a)/2,source_slope_value/4), (x_b,source_slope_value)]
shape_pts = [(x_a,source_slope_value), ((x_b+x_a)/8,4*source_slope_value), ((x_b+x_a)/4,source_slope_value/3), (3*(x_b+x_a)/4,2*source_slope_value), (x_b,source_slope_value)]
print('# shape pts = ', len(shape_pts))
source_slope = TargetFunction(shape_pts, type='legendre')
#source_slope = TargetFunction(shape_pts, type='linear')
source_slope.plot(n_plot_pts=200, show_shape_pts=True, title='Source Slope', y_label=r'$S(x)$')
# shape pts = 5
Since $w\in V(a,b)$, construct as follows $w = \sum\limits_{i=1}^N\,\alpha_i\,\phi_i$ such that $w(b) = u_b$.
'''Build the basis functions'''
degree = 1
n_elem = 1
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=n_elem, x_partition=[x_a, x_b])
lift_basis.plot(n_plot_pts=300)
'''Build the coefficients of the lift function'''
import numpy as np
alpha_vec = np.zeros(len(lift_basis.basis_func_list))
alpha_vec[0] = u_b/2 # any value
alpha_vec[-1] = u_b
'''Build the lift function'''
def w_lift(x):
a_mtrx = lift_basis.evaluation_matrix(x)
return a_mtrx@alpha_vec
'''Build the lift function derivative'''
def w_lift_prime(x):
a_mtrx = lift_basis.evaluation_matrix(x, derivative=True)
return a_mtrx@alpha_vec
'''Test boundary values'''
print('w(a) = ',w_lift(x_a), ' w(b) = ',w_lift(x_b))
w(a) = [1.] w(b) = [2.]
'''Rayleigh Ritz Method with FEM 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=(18, 5))
plt.plot(x, u_values, 'r-', label='Lift function')
plt.title(r'Rayleigh-Ritz Method with Finite Element 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 = [2,2]
n_elem = [3,10]
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+1, x_b], bc_x_max='dirichlet')
femlb.plot(n_plot_pts=300, n_legend_cols=16)
'''Build the matrix of coefficients of the linear system (using global basis functions)'''
# Too slow (keep for information)
import numpy as np
n = len(femlb.basis_func_list)
a_mtrx = np.zeros((n, n), dtype=np.float64)
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):
s_x_phi_j = lambda x: source_slope(x) * phi_j(x)
(a_ij, _) = femlb.inner_product(s_x_phi_j, phi_i)
a_mtrx[i,j] -= a_ij
for i,phi_i in enumerate(femlb.basis_func_list):
for j,phi_j in enumerate(femlb.basis_func_list):
a_mtrx[i,j] += h * phi_j(x_a) * phi_i(x_a)
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(a_mtrx) == min(a_mtrx.shape):
print('A is full rank.')
else:
print('A is rank deficient.')
A is full rank.
'''Build load vector (using global basis functions)'''
# Too slow (keep for information)
b_vec = np.zeros(n, dtype=np.float64)
for i, phi_i in enumerate(femlb.basis_func_list):
(b_vec[i], _) = femlb.inner_product(source_bias, phi_i)
(term1, _) = femlb.inner_product(w_lift_prime, femlb.basis_func_prime_list[i])
b_vec[i] -= diff_coeff * term1
s_x_w = lambda x: source_slope(x) * w_lift(x)
(term2, _) = femlb.inner_product(s_x_w, phi_i)
b_vec[i] += term2
b_vec[i] -= h * (w_lift(x_a) - u_ref) * phi_i(x_a)
'''Compute optimal coefficient vector'''
c_star_vec = np.linalg.solve(a_mtrx, b_vec)
'''Build the best approximation function in V_N'''
def u_star(x):
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 in V_N'''
def u_prime_star(x):
u_prime_0 = femlb.evaluation_matrix(x, derivative=True)@c_star_vec
w_prime = lift_basis.evaluation_matrix(x, derivative=True)@alpha_vec
return u_prime_0 + w_prime
'''Rayleigh Ritz Method with Fourier Basis Functions'''
n_pts = 500
x = np.linspace(x_a, x_b, n_pts)
u_values = u_star(x)
import matplotlib.pyplot as plt
#%matplotlib inline
#plt.style.use('dark_background')
plt.figure(1, figsize=(18, 5))
plt.plot(x, u_values, 'r-', label='Rayleigh-Ritz Solution w/ Dirichlet/Robin BC')
plt.title(r'Rayleigh-Ritz Method with Finite Element Lagrange Basis Functions (n='+str(len(femlb.basis_func_list))+')', fontsize=20)
plt.ylabel(r'$u_N(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()
'''Rayleigh-Ritz solution evaluated at the boundaries'''
np.set_printoptions(precision=5)
print('u^*_N(a) = ', u_star(x_a), 'u^*_N(b) = ', u_star(x_b), " u^*'_N(a) = ", u_prime_star(x_a))
u^*_N(a) = [5.00917] u^*_N(b) = [2.] u^*'_N(a) = [1.99502]
'''Flux at boundary'''
print('q_na computed = ', -(-diff_coeff*u_prime_star(x_a)))
print('h (u(a)-u_ref) = ', h*(u_star(x_a)-u_ref))
print('flux error [%]= ', (-(-diff_coeff*u_prime_star(x_a)) - h*(u_star(x_a)-u_ref))/(h*(u_star(x_a)-u_ref))*100)
q_na computed = [0.1995] h (u(a)-u_ref) = [0.19667] flux error [%]= [1.44137]