# Source code for L2_error

"""
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

"""

import numpy as np
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.

: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...

"""
assert bf.__class__.__name__ == 'MimeticBasisPolynomials'
assert ct.__class__.__name__ == 'CoordinateTransformation'
self.bf = bf
self.ct = ct
quad_degree = [bf.degree[i]+5 for i in range(3)]
else:
pass
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')]
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.
"""
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.
"""
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.
"""
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:
ravel=True)
detJ = self.ct.Jacobian(*self.xi_et_sg)
LEIntermediate = (v - exact_solution(*xyz))**n
localError = np.sum(LEIntermediate * detJ *
return localError**(1/n)
elif n == 'infty':
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()