"""
Compute the :math:`L^2`-error between an infinite dimensional variable
and its finite dimensional (discrete) projection in the mimetic
polynomial space.
⭕ To access the source code, click on the [source] button at the right
side or click on
:download:`[L2_error.py]</contents/LIBRARY/ptc/mathischeap_ptc/L2_error.py>`.
Dependence may exist. In case of error, check import and install required
packages or download required scripts. © mathischeap.com
"""
import numpy as np
from quadrature import Gauss
from projection import Reconstruction
[docs]
class L2Error(object):
"""A wrapper of all functions for computing :math:`L^2`-error.
:param bf: The basis functions in :math:`\\Omega_{\\mathrm{ref}}`.
:type bf: MimeticBasisPolynomials
:param ct: The coordinate transformation representing the
mapping :math:`\\Phi`,
.. math::
\\Phi: \\Omega_{\\mathrm{ref}}\\to\\Omega
:type ct: CoordinateTransformation
:param quad_degree: (default: ``None``) The degree used for the
numerical integral. It should be a list or tuple of three
positive integers. If it is ``None``, a suitable degree will be
obtained from ``bf``.
:type quad_degree: list, tuple
:example:
>>> import numpy as np
>>> from numpy import sin, cos, pi
>>> from coordinate_transformation import CoordinateTransformation
>>> from coordinate_transformation import Phi, d_Phi
>>> from mimetic_basis_polynomials import MimeticBasisPolynomials
>>> from projection import Reduction
>>> ct = CoordinateTransformation(Phi, d_Phi)
>>> bf = MimeticBasisPolynomials('Lobatto-5','Lobatto-5','Lobatto-5')
>>> def pressure(x,y,z):
... return sin(np.pi*x) * sin(pi*y) * sin(pi*z)
>>> def velocity_x(x,y,z):
... return pi * cos(pi*x) * sin(pi*y) * sin(pi*z)
>>> def velocity_y(x,y,z):
... return pi * sin(pi*x) * cos(pi*y) * sin(pi*z)
>>> def velocity_z(x,y,z):
... return pi * sin(pi*x) * sin(pi*y) * cos(pi*z)
>>> def source(x,y,z):
... return -3 * pi**2 * sin(np.pi*x) * sin(pi*y) * sin(pi*z)
>>> rd = Reduction(bf, ct)
>>> loc_dofs_N = rd.NP(pressure)
>>> loc_dofs_E = rd.EP((velocity_x, velocity_y, velocity_z))
>>> loc_dofs_F = rd.FP((velocity_x, velocity_y, velocity_z))
>>> loc_dofs_V = rd.VP(source)
>>> L2e = L2Error(bf, ct)
>>> L2e.NP(loc_dofs_N, pressure) # doctest: +ELLIPSIS
0.0228947...
>>> L2e.EP(loc_dofs_E, (velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
0.3612213...
>>> L2e.FP(loc_dofs_F, (velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
0.4178056...
>>> L2e.VP(loc_dofs_V, source) # doctest: +ELLIPSIS
2.8203143...
"""
def __init__(self, bf, ct, quad_degree=None):
assert bf.__class__.__name__ == 'MimeticBasisPolynomials'
assert ct.__class__.__name__ == 'CoordinateTransformation'
self.bf = bf
self.ct = ct
if quad_degree is None:
quad_degree = [bf.degree[i]+5 for i in range(3)]
else:
pass
self.quad_degree = quad_degree
qn0, qw0 = Gauss(quad_degree[0])
qn1, qw1 = Gauss(quad_degree[1])
qn2, qw2 = Gauss(quad_degree[2])
quad_nodes, QW = [qn0, qn1, qn2], [qw0, qw1, qw2]
xi, et, sg = np.meshgrid(*quad_nodes, indexing='ij')
xi_et_sg = [xi.ravel('F'), et.ravel('F'), sg.ravel('F')]
quad_weights = np.kron(QW[2], np.kron(QW[1], QW[0]))
self.quad_nodes = quad_nodes
self.quad_weights = quad_weights
self.xi_et_sg = xi_et_sg
self.reconstruction = Reconstruction(bf, ct)
[docs]
def NP(self, loc_dofs, exact_solution):
"""Let :math:`\\psi\in H^1(\\Omega)` and
:math:`\\psi^h=\\pi\\left(\\psi\\right)\\in
\\text{NP}_{N}(\\Omega)`, we compute
:math:`\\left\|\\psi^h-\\psi\\right\|_{L^2\\text{-error}}`.
:param loc_dofs: The local dofs of :math:`\\psi^h`.
:param exact_solution: The scalar/function :math:`\\psi`.
:returns: A float representing the :math:`L^2`-error.
"""
xyz, v = self.reconstruction.NP(loc_dofs, *self.quad_nodes,
ravel=True)
detJ = self.ct.Jacobian(*self.xi_et_sg)
LEIntermediate = (v - exact_solution(*xyz))**2
localError = np.sum(LEIntermediate * detJ * self.quad_weights)
return localError**0.5
[docs]
def EP(self, loc_dofs, exact_solution):
"""Let :math:`\\boldsymbol{\\omega}\in H(\\mathrm{curl};\\Omega)` and
:math:`\\boldsymbol{\\omega}^h=\\pi\\left(
\\boldsymbol{\\omega}\\right)\\in
\\text{EP}_{N-1}(\\Omega)`, we compute
:math:`\\left\|\\boldsymbol{\\omega}^h-\\boldsymbol{\\omega}
\\right\|_{L^2\\text{-error}}`.
:param loc_dofs: The local dofs of
:math:`\\boldsymbol{\\omega}^h`.
:param exact_solution: A tuple of three functions representing
the three components of the vector
:math:`\\boldsymbol{\\omega}`.
:returns: A float representing the :math:`L^2`-error.
"""
xyz, v = self.reconstruction.EP(loc_dofs, *self.quad_nodes,
ravel=True)
detJ = self.ct.Jacobian(*self.xi_et_sg)
LEIntermediate = (v[0] - exact_solution[0](*xyz))**2 \
+ (v[1] - exact_solution[1](*xyz))**2 \
+ (v[2] - exact_solution[2](*xyz))**2
localError = np.sum(LEIntermediate * detJ * self.quad_weights)
return localError**0.5
[docs]
def FP(self, loc_dofs, exact_solution):
"""Let :math:`\\boldsymbol{\\omega}\in H(\\mathrm{div};\\Omega)` and
:math:`\\boldsymbol{u}^h=\\pi\\left(
\\boldsymbol{u}\\right)\\in
\\text{FP}_{N-1}(u)`, we compute
:math:`\\left\|\\boldsymbol{u}^h-\\boldsymbol{u}
\\right\|_{L^2\\text{-error}}`.
:param loc_dofs: The local dofs of
:math:`\\boldsymbol{u}^h`.
:param exact_solution: A tuple of three functions representing
the three components of the vector
:math:`\\boldsymbol{u}`.
:returns: A float representing the :math:`L^2`-error.
"""
xyz, v = self.reconstruction.FP(loc_dofs, *self.quad_nodes,
ravel=True)
detJ = self.ct.Jacobian(*self.xi_et_sg)
LEIntermediate = (v[0] - exact_solution[0](*xyz))**2 \
+ (v[1] - exact_solution[1](*xyz))**2 \
+ (v[2] - exact_solution[2](*xyz))**2
localError = np.sum(LEIntermediate * detJ * self.quad_weights)
return localError**0.5
[docs]
def VP(self, loc_dofs, exact_solution, n=2):
"""Let :math:`f\in L^2(\\Omega)` and
:math:`f^h=\\pi\\left(f\\right)\\in
\\text{VP}_{N-1}(\\Omega)`, we compute
:math:`\\left\|f^h-f\\right\|_{L^n\\text{-error}}`.
:param loc_dofs: The local dofs of :math:`f^h`.
:param exact_solution: The scalar/function :math:`f`.
:param n: (default: :math:`2`) We compute :math:`L^{n}`-error.
:returns: A float representing the :math:`L^n`-error.
"""
if n == 2:
xyz, v = self.reconstruction.VP(loc_dofs, *self.quad_nodes,
ravel=True)
detJ = self.ct.Jacobian(*self.xi_et_sg)
LEIntermediate = (v - exact_solution(*xyz))**n
localError = np.sum(LEIntermediate * detJ *
self.quad_weights)
return localError**(1/n)
elif n == 'infty':
qn0, _ = Gauss(self.quad_degree[0]+3)
qn1, _ = Gauss(self.quad_degree[1]+3)
qn2, _ = Gauss(self.quad_degree[2]+3)
quad_nodes = [qn0, qn1, qn2]
xyz, v = self.reconstruction.VP(loc_dofs, *quad_nodes,
ravel=True)
diff = np.abs(v - exact_solution(*xyz))
return np.max(diff)
else:
raise NotImplementedError(f'not implemented for n={n}.')
if __name__ == '__main__':
import doctest
doctest.testmod()
from numpy import sin, cos, pi
from coordinate_transformation import CoordinateTransformation
from coordinate_transformation import Phi, d_Phi
from mimetic_basis_polynomials import MimeticBasisPolynomials
from projection import Reduction
ct = CoordinateTransformation(Phi, d_Phi)
def pressure(x,y,z):
return sin(np.pi*x) * sin(pi*y) * sin(pi*z)
def velocity_x(x,y,z):
return pi * cos(pi*x) * sin(pi*y) * sin(pi*z)
def velocity_y(x,y,z):
return pi * sin(pi*x) * cos(pi*y) * sin(pi*z)
def velocity_z(x,y,z):
return pi * sin(pi*x) * sin(pi*y) * cos(pi*z)
def source(x,y,z):
return -3 * pi**2 * sin(np.pi*x) * sin(pi*y) * sin(pi*z)
EN, EE, EF, EV = list(), list(), list(), list()
Ns = (3, 4, 5, 6, 7, 8, 9, 10)
for N in Ns:
bf = MimeticBasisPolynomials('Lobatto-' + str(N),
'Lobatto-' + str(N),
'Lobatto-' + str(N))
rd = Reduction(bf, ct)
loc_dofs_N = rd.NP(pressure)
loc_dofs_E = rd.EP((velocity_x, velocity_y, velocity_z))
loc_dofs_F = rd.FP((velocity_x, velocity_y, velocity_z))
loc_dofs_V = rd.VP(source)
L2e = L2Error(bf, ct)
eN = L2e.NP(loc_dofs_N, pressure)
eE = L2e.EP(loc_dofs_E, (velocity_x, velocity_y, velocity_z))
eF = L2e.FP(loc_dofs_F, (velocity_x, velocity_y, velocity_z))
eV = L2e.VP(loc_dofs_V, source)
EN.append(np.log10(eN))
EE.append(np.log10(eE))
EF.append(np.log10(eF))
EV.append(np.log10(eV))
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(12, 12))
plt.subplot2grid((2, 2), (0, 0))
plt.plot(Ns, EN)
plt.title('NP L2-error; exponential convergence')
plt.xlabel('N')
plt.ylabel(r'$\log(L2-error)$')
plt.subplot2grid((2, 2), (0, 1))
plt.plot(Ns, EE)
plt.title('EP L2-error; exponential convergence')
plt.xlabel('N')
plt.ylabel(r'$\log(L2-error)$')
plt.subplot2grid((2, 2), (1, 0))
plt.plot(Ns, EF)
plt.title('FP L2-error; exponential convergence')
plt.xlabel('N')
plt.ylabel(r'$\log(L2-error)$')
plt.subplot2grid((2, 2), (1, 1))
plt.plot(Ns, EV)
plt.title('VP L2-error; exponential convergence')
plt.xlabel('N')
plt.ylabel(r'$\log(L2-error)$')
plt.show()
