# Source code for L2_error_dual

"""
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 quadrature import Gauss
from projection_dual import ReconstructionDual
from projection import Reconstruction

[docs]class L2ErrorDual(object):
"""A wrapper of all functions for computing :math:L^2-error. in
the dual spaces
:math:\\widetilde{\\text{NP}}_{N}(\\Omega),
:math:\\widetilde{\\text{EP}}_{N-1}(\\Omega),
:math:\\widetilde{\\text{FP}}_{N-1}(\\Omega) and
:math:\\widetilde{\\text{VP}}_{N-1}(\\Omega).

: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_dual import ReductionDual
>>> 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 = ReductionDual(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 = L2ErrorDual(bf, ct)
>>> L2e.NP(loc_dofs_N, pressure) # doctest: +ELLIPSIS
0.0228947239...
>>> L2e.EP(loc_dofs_E, (velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
0.3612213119...
>>> L2e.FP(loc_dofs_F, (velocity_x, velocity_y, velocity_z)) # doctest: +ELLIPSIS
0.4178056796...
>>> L2e.VP(loc_dofs_V, source) # doctest: +ELLIPSIS
2.8203143473...

"""
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
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.xi_et_sg = xi_et_sg
self.reconstruction = ReconstructionDual(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
\\widetilde{\\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^1(\\Omega) and
:math:\\boldsymbol{\\omega}^h=\\pi\\left(
\\boldsymbol{\\omega}\\right)\\in
\\widetilde{\\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^1(\\Omega) and
:math:\\boldsymbol{u}^h=\\pi\\left(
\\boldsymbol{u}\\right)\\in
\\widetilde{\\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
\\widetilde{\\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 *
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}.')

def _FP_diff_(self, loc_dofs, primal_dofs, exact_solution):
primal_reconstruction = Reconstruction(self.bf, self.ct)
xyz, v = self.reconstruction.FP(loc_dofs, *self.quad_nodes,
ravel=True)
___, V = primal_reconstruction.FP(primal_dofs, *self.quad_nodes,
ravel=True)
detJ = self.ct.Jacobian(*self.xi_et_sg)
LEIntermediate = ((v[0]-V[0]) - exact_solution[0](*xyz))**2 \
+ ((v[1]-V[1]) - exact_solution[1](*xyz))**2 \
+ ((v[2]-V[2]) - exact_solution[2](*xyz))**2
localError = np.sum(LEIntermediate * detJ * self.quad_weights)
return localError**0.5

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_dual import ReductionDual
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)
for N in Ns:
print(f"{N} of {Ns}...")
bf = MimeticBasisPolynomials(
'Lobatto-' + str(N), 'Lobatto-' + str(N), 'Lobatto-' + str(N))

rdd = ReductionDual(bf, ct)
loc_dofs_N = rdd.NP(pressure)
loc_dofs_E = rdd.EP((velocity_x, velocity_y, velocity_z))
loc_dofs_F = rdd.FP((velocity_x, velocity_y, velocity_z))
loc_dofs_V = rdd.VP(source)
L2ed = L2ErrorDual(bf, ct)
eN = L2ed.NP(loc_dofs_N, pressure)
eE = L2ed.EP(loc_dofs_E, (velocity_x, velocity_y, velocity_z))
eF = L2ed.FP(loc_dofs_F, (velocity_x, velocity_y, velocity_z))
eV = L2ed.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('dual(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('dual(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('dual(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('dual(VP) L2-error; exponential convergence')
plt.xlabel('N')
plt.ylabel(r'$\log(L2-error)$')

plt.show()