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
packages or download required scripts. © mathischeap.com

"""

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 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 = 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 * 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}.')
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()