diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index 597baa731..ddf01893b 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -563,6 +563,12 @@ def discretize(a, *args, **kwargs): if isinstance(a, sym_BasicForm): if isinstance(a, (sym_Norm, sym_SemiNorm)): + if hasattr(a._space, 'codomain_type') and a._space.codomain_type == 'complex': + print('WARNING (56436574):', + 'discretization of norms is not correctly tested yet for complex fields: see https://github.com/campospinto/psydac_dev/issues/12.') + if isinstance(a, sym_SemiNorm): + print('WARNING (87676545):', + 'discretization of semi-norms is not correctly tested yet: see https://github.com/campospinto/psydac_dev/issues/12.') kernel_expr = TerminalExpr(a, domain) if not mapping is None: kernel_expr = tuple(LogicalExpr(i, domain) for i in kernel_expr) diff --git a/psydac/api/tests/test_2d_complex.py b/psydac/api/tests/test_2d_complex.py index 57e6a37db..495a8ac4f 100644 --- a/psydac/api/tests/test_2d_complex.py +++ b/psydac/api/tests/test_2d_complex.py @@ -1,11 +1,14 @@ # -*- coding: UTF-8 -*- import os +from collections import OrderedDict + from mpi4py import MPI import pytest import numpy as np from sympy import pi, cos, sin, symbols, conjugate, exp from sympy import Tuple, Matrix +from sympy import lambdify from sympde.calculus import grad, dot, cross, curl from sympde.calculus import minus, plus @@ -224,7 +227,7 @@ def run_helmholtz_2d(solution, kappa, e_w_0, dx_e_w_0, domain, ncells=None, degr equation = find(u, forall=v, lhs=a(u,v), rhs=l(v)) l2norm = Norm(error, domain, kind='l2') - h1norm = SemiNorm(error, domain, kind='h1') + h1norm = SemiNorm(error, domain, kind='h1') # [MCP 19.07.2024] This is probably wrong: see https://github.com/campospinto/psydac_dev/issues/12 #+++++++++++++++++++++++++++++++ # 2. Discretization @@ -244,7 +247,7 @@ def run_helmholtz_2d(solution, kappa, e_w_0, dx_e_w_0, domain, ncells=None, degr l2_error = l2norm_h.assemble(u=uh) h1_error = h1norm_h.assemble(u=uh) - return l2_error, h1_error + return l2_error, h1_error, uh #============================================================================== def run_maxwell_2d(uex, f, alpha, domain, *, ncells=None, degree=None, filename=None, comm=None): @@ -400,7 +403,7 @@ def test_complex_poisson_2d_multipatch_mapping(): assert abs(h1_error - expected_h1_error) < 1e-7 -def test_complex_helmholtz_2d(): +def test_complex_helmholtz_2d(plot_sol=False): # This test solves the homogeneous Helmhotz equation with impedance BC. # In particular, we impose an incoming wave from the left and absorbing boundary conditions at the right. # Along y, periodic boundary conditions are considered. @@ -412,13 +415,50 @@ def test_complex_helmholtz_2d(): e_w_0 = sin(kappa * y) # value of incoming wave at x=0, forall y dx_e_w_0 = 1j*kappa*sin(kappa * y) # derivative wrt. x of incoming wave at x=0, forall y - l2_error, h1_error = run_helmholtz_2d(solution, kappa, e_w_0, dx_e_w_0, domain, ncells=[2**2,2**2], degree=[2,2]) + l2_error, h1_error, uh = run_helmholtz_2d(solution, kappa, e_w_0, dx_e_w_0, domain, ncells=[2**2,2**2], degree=[2,2]) expected_l2_error = 0.01540947560953227 - expected_h1_error = 0.19040207344639598 + # commented because semi-norms of complex-fields are not working properly for now: see https://github.com/campospinto/psydac_dev/issues/12 + expected_h1_error = 'NOT IMPLEMENTED' # + + print(f'errors: l2 = {l2_error}, h1 = {h1_error} (THIS IS PROBABLY WRONG)') + print(f'expected errors: l2 = {expected_l2_error}, h1 = {expected_h1_error}') + + if plot_sol: + from psydac.feec.multipatch.plotting_utilities import get_plotting_grid, get_grid_vals + from psydac.feec.multipatch.plotting_utilities import get_patch_knots_gridlines, my_small_plot + from psydac.feec.pull_push import pull_2d_h1 + + Id_mapping = IdentityMapping('M', 2) + mappings = OrderedDict([(domain, Id_mapping)]) + mappings_list = [m for m in mappings.values()] + call_mappings_list = [m.get_callable_mapping() for m in mappings_list] + + uh = [uh] # single-patch cast as multi-patch solution + + u = lambdify(domain.coordinates, solution) + u_log = [pull_2d_h1(u, f) for f in call_mappings_list] + + etas, xx, yy = get_plotting_grid(mappings, N=20) + grid_vals_h1 = lambda v: get_grid_vals(v, etas, mappings_list, space_kind='h1') + + uh_vals = grid_vals_h1(uh) + u_vals = grid_vals_h1(u_log) + + u_err = [(u1 - u2) for u1, u2 in zip(u_vals, uh_vals)] + + my_small_plot( + title=r'approximation of solution $u$', + vals=[u_vals, uh_vals, u_err], + titles=[r'$u^{ex}(x,y)$', r'$u^h(x,y)$', r'$|(u^{ex}-u^h)(x,y)|$'], + xx=xx, + yy=yy, + gridlines_x1=None, + gridlines_x2=None, + ) assert( abs(l2_error - expected_l2_error) < 1.e-7) - assert( abs(h1_error - expected_h1_error) < 1.e-7) + # assert( abs(h1_error - expected_h1_error) < 1.e-7) def test_maxwell_2d_2_patch_dirichlet_2(): # This test solve the maxwell problem with non-homogeneous dirichlet condition with penalization on the border of the exact solution @@ -487,59 +527,67 @@ def teardown_function(): if __name__ == '__main__': - from collections import OrderedDict - - from sympy import lambdify - - from psydac.feec.multipatch.plotting_utilities import get_plotting_grid, get_grid_vals - from psydac.feec.multipatch.plotting_utilities import get_patch_knots_gridlines, my_small_plot - from psydac.api.tests.build_domain import build_pretzel - from psydac.feec.pull_push import pull_2d_hcurl + # we do some visual verifications... - domain = build_pretzel() - x,y = domain.coordinates + verify = 'helmholtz' - omega = 1.5 - alpha = -omega**2 - Eex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y)) - f = Tuple(alpha*sin(pi*y) - pi**2*sin(pi*y)*cos(pi*x) + pi**2*sin(pi*y), - alpha*sin(pi*x)*cos(pi*y) + pi**2*sin(pi*x)*cos(pi*y)) + if verify == 'helmholtz': + test_complex_helmholtz_2d(plot_sol=True) - l2_error, Eh = run_maxwell_2d(Eex, f, alpha, domain, ncells=[2**2, 2**2], degree=[2,2]) - - mappings = OrderedDict([(P.logical_domain, P.mapping) for P in domain.interior]) - mappings_list = list(mappings.values()) - mappings_list = [mapping.get_callable_mapping() for mapping in mappings_list] - - Eex_x = lambdify(domain.coordinates, Eex[0]) - Eex_y = lambdify(domain.coordinates, Eex[1]) - Eex_log = [pull_2d_hcurl([Eex_x,Eex_y], f) for f in mappings_list] - - etas, xx, yy = get_plotting_grid(mappings, N=20) - grid_vals_hcurl = lambda v: get_grid_vals(v, etas, mappings_list, space_kind='hcurl') - - Eh_x_vals, Eh_y_vals = grid_vals_hcurl(Eh) - E_x_vals, E_y_vals = grid_vals_hcurl(Eex_log) - - E_x_err = [(u1 - u2) for u1, u2 in zip(E_x_vals, Eh_x_vals)] - E_y_err = [(u1 - u2) for u1, u2 in zip(E_y_vals, Eh_y_vals)] - - my_small_plot( - title=r'approximation of solution $u$, $x$ component', - vals=[E_x_vals, Eh_x_vals, E_x_err], - titles=[r'$u^{ex}_x(x,y)$', r'$u^h_x(x,y)$', r'$|(u^{ex}-u^h)_x(x,y)|$'], - xx=xx, - yy=yy, - gridlines_x1=None, - gridlines_x2=None, - ) - - my_small_plot( - title=r'approximation of solution $u$, $y$ component', - vals=[E_y_vals, Eh_y_vals, E_y_err], - titles=[r'$u^{ex}_y(x,y)$', r'$u^h_y(x,y)$', r'$|(u^{ex}-u^h)_y(x,y)|$'], - xx=xx, - yy=yy, - gridlines_x1=None, - gridlines_x2=None, - ) + else: + from collections import OrderedDict + + from sympy import lambdify + + from psydac.feec.multipatch.plotting_utilities import get_plotting_grid, get_grid_vals + from psydac.feec.multipatch.plotting_utilities import get_patch_knots_gridlines, my_small_plot + from psydac.api.tests.build_domain import build_pretzel + from psydac.feec.pull_push import pull_2d_hcurl + + domain = build_pretzel() + x,y = domain.coordinates + + omega = 1.5 + alpha = -omega**2 + Eex = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y)) + f = Tuple(alpha*sin(pi*y) - pi**2*sin(pi*y)*cos(pi*x) + pi**2*sin(pi*y), + alpha*sin(pi*x)*cos(pi*y) + pi**2*sin(pi*x)*cos(pi*y)) + + l2_error, Eh = run_maxwell_2d(Eex, f, alpha, domain, ncells=[2**2, 2**2], degree=[2,2]) + + mappings = OrderedDict([(P.logical_domain, P.mapping) for P in domain.interior]) + mappings_list = list(mappings.values()) + call_mappings_list = [m.get_callable_mapping() for m in mappings_list] + + Eex_x = lambdify(domain.coordinates, Eex[0]) + Eex_y = lambdify(domain.coordinates, Eex[1]) + Eex_log = [pull_2d_hcurl([Eex_x,Eex_y], f) for f in call_mappings_list] + + etas, xx, yy = get_plotting_grid(mappings, N=20) + grid_vals_hcurl = lambda v: get_grid_vals(v, etas, mappings_list, space_kind='hcurl') + + Eh_x_vals, Eh_y_vals = grid_vals_hcurl(Eh) + E_x_vals, E_y_vals = grid_vals_hcurl(Eex_log) + + E_x_err = [(u1 - u2) for u1, u2 in zip(E_x_vals, Eh_x_vals)] + E_y_err = [(u1 - u2) for u1, u2 in zip(E_y_vals, Eh_y_vals)] + + my_small_plot( + title=r'approximation of solution $u$, $x$ component', + vals=[E_x_vals, Eh_x_vals, E_x_err], + titles=[r'$u^{ex}_x(x,y)$', r'$u^h_x(x,y)$', r'$|(u^{ex}-u^h)_x(x,y)|$'], + xx=xx, + yy=yy, + gridlines_x1=None, + gridlines_x2=None, + ) + + my_small_plot( + title=r'approximation of solution $u$, $y$ component', + vals=[E_y_vals, Eh_y_vals, E_y_err], + titles=[r'$u^{ex}_y(x,y)$', r'$u^h_y(x,y)$', r'$|(u^{ex}-u^h)_y(x,y)|$'], + xx=xx, + yy=yy, + gridlines_x1=None, + gridlines_x2=None, + ) diff --git a/psydac/api/tests/test_assembly.py b/psydac/api/tests/test_assembly.py index 853935199..310629c8c 100644 --- a/psydac/api/tests/test_assembly.py +++ b/psydac/api/tests/test_assembly.py @@ -1,14 +1,13 @@ import pytest import numpy as np from mpi4py import MPI -from sympy import pi, sin, cos, tan, atan, atan2, exp, sinh, cosh, tanh, atanh, Tuple, I - - +from sympy import Domain, pi, sin, cos, tan, atan, atan2, exp, sinh, cosh, tanh, atanh, Tuple, I, sqrt +from sympy import lambdify from sympde.topology import Line, Square from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace -from sympde.topology import element_of, Derham +from sympde.topology import element_of, Derham, elements_of from sympde.core import Constant -from sympde.expr import LinearForm, BilinearForm, Functional, Norm +from sympde.expr import LinearForm, BilinearForm, Functional, Norm, SemiNorm from sympde.expr import integral from sympde.calculus import Inner @@ -530,8 +529,272 @@ def test_assembly_no_synchr_args(backend): assert( abs(inte_lin) < 1.e-12) assert( abs(inte_norm) < 1.e-12) +# def norm_projected(expr, Vh, nquads): +# from psydac.feec.global_projectors import Projector_L2 +# P = Projector_L2(Vh, nquads=nquads) +# solution_call = lambdify(Domain.coordinates, solution) +# sol_h = P(solution_call) +# sol_c = sol_h.coeffs +# a = BilinearForm((u, v), integral(domain, u*v)) +# a_h = discretize(a, domain_h, [Vh,Vh], backend=None) +# M = a_h.assemble() +# l2_norm_psydac = np.sqrt(sol_c.dot(M.dot(sol_c))) + + + +def test_sympde_norm(backend): + + kwargs = {'backend': PSYDAC_BACKENDS[backend]} if backend else {} + + domain = Square('domain', bounds1=(0, 1), bounds2=(0, 1)) + + U = ScalarFunctionSpace('U', domain, kind='h1') + U.codomain_type='real' + V = ScalarFunctionSpace('V', domain, kind='h1') + V.codomain_type='complex' + u1, u2 = elements_of(U, names='u1, u2') + v1, v2 = elements_of(V, names='v1, v2') + + x, y = domain.coordinates + kappa = 2*pi + rho = sin(kappa * y) # real + dx_rho = 0 * y + dy_rho = kappa * cos(kappa * y) + phi = exp(1j * kappa * x) * rho # complex + # other possible expressions? also fail: + # phi = exp(I * kappa * x) * rho + # phi = (cos(kappa * x) + 1j * sin(kappa * x)) * rho + # phi = (cos(kappa * x) + 1.j * sin(kappa * x)) * rho + dx_phi = 1j * kappa * phi + dy_phi = exp(1j * kappa * x) * dy_rho + + # sympde L2 norms + rho_l2_sym = Norm(rho, domain, kind='l2') + phi_l2_sym = Norm(phi, domain, kind='l2') + # sympde H1 semi-norms + rho_h1s_sym = SemiNorm(rho, domain, kind='h1') + phi_h1s_sym = SemiNorm(phi, domain, kind='h1') + + # Q: can we evaluate these norms directly (in sympde)? + + # aux + rho_l2_usym = Norm(rho - u1, domain, kind='l2') + rho_l2_vsym = Norm(rho - v1, domain, kind='l2') + phi_l2_vsym = Norm(phi - v1, domain, kind='l2') + rho_h1s_usym = SemiNorm(rho - u1, domain, kind='h1') + rho_h1s_vsym = SemiNorm(rho - v1, domain, kind='h1') + phi_h1s_vsym = SemiNorm(phi - v1, domain, kind='h1') + + # aux 2 + l2_usym = Norm(u1, domain, kind='l2') + l2_vsym = Norm(v1, domain, kind='l2') + h1s_usym = SemiNorm(u1, domain, kind='h1') + h1s_vsym = SemiNorm(v1, domain, kind='h1') + + # exact norms + rho_l2_ex = sqrt((1-sin(2*kappa)/(2*kappa))/2) + rho_h1s_ex = kappa * sqrt(1-rho_l2_ex**2) # todo + rho_h1_ex = sqrt(rho_l2_ex**2 + rho_h1s_ex**2) + + phi_l2_ex = rho_l2_ex + phi_h1s_ex = kappa + phi_h1_ex = sqrt(phi_l2_ex**2 + phi_h1s_ex**2) + + # discretize norms + ncells = [8, 8] + degree = [4, 4] + domain_h = discretize(domain, ncells=ncells, periodic=[False, False]) + Uh = discretize(U, domain_h, degree=degree) + Vh = discretize(V, domain_h, degree=degree) + + # commented because currently not working. TODO: fix this + # rho_l2_h = discretize(rho_l2_sym, domain_h, Uh, **kwargs) # todo: also try in Vh + # phi_l2_h = discretize(phi_l2_sym, domain_h, Vh, **kwargs) + # rho_l2 = rho_l2_h.assemble() + # phi_l2 = phi_l2_h.assemble() + rho_l2 = 'NOT IMPLEMENTED (todo)' + phi_l2 = 'NOT IMPLEMENTED (todo)' + + # aux + rho_l2_usym_h = discretize(rho_l2_usym, domain_h, Uh, **kwargs) + zero_h = FemField(Uh) + rho_l2_u0 = rho_l2_usym_h.assemble(u1=zero_h) + + rho_l2_vsym_h = discretize(rho_l2_vsym, domain_h, Vh, **kwargs) + zero_h = FemField(Vh) + rho_l2_v0 = rho_l2_vsym_h.assemble(v1=zero_h) + + phi_l2_vsym_h = discretize(phi_l2_vsym, domain_h, Vh, **kwargs) + zero_h = FemField(Vh) + phi_l2_v0 = phi_l2_vsym_h.assemble(v1=zero_h) + + rho_h1s_usym_h = discretize(rho_h1s_usym, domain_h, Uh, **kwargs) + zero_h = FemField(Uh) + rho_h1s_u0 = rho_h1s_usym_h.assemble(u1=zero_h) + + rho_h1s_vsym_h = discretize(rho_h1s_vsym, domain_h, Vh, **kwargs) + zero_h = FemField(Vh) + rho_h1s_v0 = rho_h1s_vsym_h.assemble(v1=zero_h) + + phi_h1s_vsym_h = discretize(phi_h1s_vsym, domain_h, Vh, **kwargs) + zero_h = FemField(Vh) + phi_h1s_v0 = phi_h1s_vsym_h.assemble(v1=zero_h) + + # discretize norms through L2 projection + m_U = BilinearForm((u1, v1), integral(domain, u1*v1)) + m_V = BilinearForm((u2, v2), integral(domain, u2*v2)) + mh_U = discretize(m_U, domain_h, [Uh,Uh], **kwargs) + mh_V = discretize(m_V, domain_h, [Vh,Vh], **kwargs) + M_U = mh_U.assemble() + M_V = mh_V.assemble() + + from psydac.feec.global_projectors import Projector_L2 + P_U = Projector_L2(Uh, nquads=[2*d+1 for d in degree]) + P_V = Projector_L2(Vh, nquads=[2*d+1 for d in degree]) + + def proj(fun_sym, space = 'U'): + fun = lambdify(domain.coordinates, fun_sym) + if space == 'U': + return P_U(fun) + elif space == 'V': + return P_V(fun) + else: + raise ValueError('space must be either "U" or "V"') + + rho_h = proj(rho, space = 'U') + dx_rho_h = proj(dx_rho, space = 'U') + dy_rho_h = proj(dy_rho, space = 'U') + vrho_h = proj(rho, space = 'V') + dx_vrho_h = proj(dx_rho, space = 'V') + dy_vrho_h = proj(dy_rho, space = 'V') + phi_h = proj(phi, space = 'V') + dx_phi_h = proj(dx_phi, space = 'V') + dy_phi_h = proj(dy_phi, space = 'V') + + rho_c = rho_h.coeffs + dx_rho_c = dx_rho_h.coeffs + dy_rho_c = dy_rho_h.coeffs + vrho_c = vrho_h.coeffs + dx_vrho_c = dx_vrho_h.coeffs + dy_vrho_c = dy_vrho_h.coeffs + phi_c = phi_h.coeffs + dx_phi_c = dx_phi_h.coeffs + dy_phi_c = dy_phi_h.coeffs + + rho_h_l2 = np.sqrt(rho_c.dot(M_U.dot(rho_c))) + rho_h_h1s = np.sqrt(dx_rho_c.dot(M_U.dot(dx_rho_c)) + dy_rho_c.dot(M_U.dot(dy_rho_c))) + + vrho_h_l2 = np.sqrt(vrho_c.dot(M_V.dot(vrho_c))) + vrho_h_h1s = np.sqrt(dx_vrho_c.dot(M_V.dot(dx_vrho_c)) + dy_vrho_c.dot(M_V.dot(dy_vrho_c))) + + phi_h_l2 = np.sqrt(phi_c.dot(M_V.dot(phi_c))) + phi_h_h1s = np.sqrt(dx_phi_c.dot(M_V.dot(dx_phi_c)) + dy_phi_c.dot(M_V.dot(dy_phi_c))) + + # aux 2 + l2_usym_h = discretize(l2_usym, domain_h, Uh, **kwargs) # todo: try with V + h1s_usym_h = discretize(h1s_usym, domain_h, Uh, **kwargs) + l2_vsym_h = discretize(l2_vsym, domain_h, Vh, **kwargs) + h1s_vsym_h = discretize(h1s_vsym, domain_h, Vh, **kwargs) + + l2_urho = l2_usym_h.assemble(u1=rho_h) + h1s_urho = h1s_usym_h.assemble(u1=rho_h) + + l2_vrho = l2_vsym_h.assemble(v1=vrho_h) + h1s_vrho = h1s_vsym_h.assemble(v1=vrho_h) + + l2_vphi = l2_vsym_h.assemble(v1=phi_h) + h1s_vphi = h1s_vsym_h.assemble(v1=phi_h) + + # aux 3 + phi_I, phi_R = phi.as_real_imag() + phi_I_l2_sym = Norm(phi_I - u1, domain, kind='l2') + phi_I_l2_sym_h = discretize(phi_I_l2_sym, domain_h, Uh, **kwargs) + zero_h = FemField(Uh) + phi_I_l2_u0 = phi_I_l2_sym_h.assemble(u1=zero_h) + phi_R_l2_sym = Norm(phi_R - u1, domain, kind='l2') + phi_R_l2_sym_h = discretize(phi_R_l2_sym, domain_h, Uh, **kwargs) + zero_h = FemField(Uh) + phi_R_l2_u0 = phi_R_l2_sym_h.assemble(u1=zero_h) + phi_IR_l2_u0 = np.sqrt(phi_I_l2_u0**2 + phi_R_l2_u0**2) + + # compare + run_fails = False + + tol = 1e-12 + htol = 1e-4 + + print(" ---- ---- ---- ---- ---- ---- ---- ") + print(" ---- ---- ---- ---- ---- ---- ---- ") + print(" rho L2 norms: ") + rho_l2_ref = rho_l2_ex.evalf() + print(f'rho_l2_ref = {rho_l2_ref}') + print(f'rho_l2 = {rho_l2}') + print(f'rho_l2_u0 = {rho_l2_u0}') + print(f'rho_l2_v0 = {rho_l2_v0}') + print(f'rho_h_l2 = {rho_h_l2}') + print(f'vrho_h_l2 = {vrho_h_l2}') + print(f'l2_urho = {l2_urho}') + print(f'l2_vrho = {l2_vrho}') + # assert abs(rho_l2 - rho_l2_ref) < tol # TODO + assert abs(rho_l2_u0 - rho_l2_ref) < tol + assert abs(rho_l2_v0 - rho_l2_ref) < tol + assert abs(rho_h_l2 - rho_l2_ref) < htol + assert abs(vrho_h_l2 - rho_l2_ref) < htol + assert abs(l2_urho - rho_l2_ref) < htol + assert abs(l2_vrho - rho_l2_ref) < htol + print(" ---- ---- ---- ---- ---- ---- ---- ") + print(" phi L2 norms: ") + phi_l2_ref = phi_l2_ex.evalf() + print(f'phi_l2_ref = {phi_l2_ref}') + print(f'phi_l2 = {phi_l2}') + print(f'phi_l2_v0 = {phi_l2_v0}') + print(f'phi_h_l2 = {phi_h_l2}') + print(f'l2_vphi = {l2_vphi}') + print(f'phi_IR_l2_u0 = {phi_IR_l2_u0}') + # assert abs(phi_l2 - phi_l2_ref) < tol # TODO + if run_fails: assert abs(phi_l2_v0 - phi_l2_ref) < tol # FAILS! + assert abs(phi_h_l2 - phi_l2_ref) < htol + if run_fails: assert abs(l2_vphi - phi_l2_ref) < htol # FAILS! + assert abs(phi_IR_l2_u0 - phi_l2_ref) < tol + + print(" ---- ---- ---- ---- ---- ---- ---- ") + print(" ---- ---- ---- ---- ---- ---- ---- ") + print(" rho H1-seminorms: ") + rho_h1s_ref = rho_h1s_ex.evalf() + print(f'rho_h1s_ref = {rho_h1s_ref}') + print(f'rho_h1s_u0 = {rho_h1s_u0}') + print(f'rho_h1s_v0 = {rho_h1s_v0}') + print(f'rho_h_h1s = {rho_h_h1s}') + print(f'vrho_h_h1s = {vrho_h_h1s}') + print(f'h1s_urho = {h1s_urho}') + print(f'h1s_vrho = {h1s_vrho}') + if run_fails: assert abs(rho_h1s_u0 - rho_h1s_ref) < tol # FAILS! + if run_fails: assert abs(rho_h1s_v0 - rho_h1s_ref) < tol # FAILS! + assert abs(rho_h_h1s - rho_h1s_ref) < htol + assert abs(vrho_h_h1s - rho_h1s_ref) < htol + assert abs(h1s_urho - rho_h1s_ref) < htol + assert abs(h1s_vrho - rho_h1s_ref) < htol + print(" ---- ---- ---- ---- ---- ---- ---- ") + print(" phi H1-seminorms: ") + phi_h1s_ref = phi_h1s_ex.evalf() + print(f'phi_h1s_ref = {phi_h1s_ref}') + print(f'phi_h1s_v0 = {phi_h1s_v0}') + print(f'phi_h_h1s = {phi_h_h1s}') + print(f'h1s_vphi = {h1s_vphi}') + if run_fails: assert abs(phi_h1s_v0 - phi_h1s_ref) < tol # FAILS! + assert abs(phi_h_h1s - phi_h1s_ref) < htol + if run_fails: assert abs(h1s_vphi - phi_h1s_ref) < htol # FAILS! + print(" ---- ---- ---- ---- ---- ---- ---- ") + print(" ---- ---- ---- ---- ---- ---- ---- ") + + + #============================================================================== if __name__ == '__main__': + + test_Norm_complex(None) + test_sympde_norm(None) + exit() test_field_and_constant(None) test_multiple_fields(None) test_math_imports(None)