From 26578f6b12d16351a281f27e55f53514815360ea Mon Sep 17 00:00:00 2001 From: Frederik Schnack Date: Wed, 2 Aug 2023 14:43:26 +0200 Subject: [PATCH] add new tests and last example --- .../examples/h1_source_pbms_conga_2d.py | 4 +- .../examples/hcurl_source_pbms_conga_2d.py | 5 +- .../examples_nc/h1_source_pbms_nc.py | 278 ++++++++++++++++++ .../examples_nc/hcurl_eigen_pbms_dg.py | 2 +- .../examples_nc/hcurl_eigen_pbms_nc.py | 4 +- .../examples_nc/hcurl_eigen_testcase.py | 8 +- .../examples_nc/hcurl_source_pbms_nc.py | 23 +- .../examples_nc/hcurl_source_testcase.py | 1 - .../tests/test_feec_maxwell_multipatch_2d.py | 162 +++++++++- .../tests/test_feec_poisson_multipatch_2d.py | 30 +- 10 files changed, 493 insertions(+), 24 deletions(-) create mode 100644 psydac/feec/multipatch/examples_nc/h1_source_pbms_nc.py diff --git a/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py b/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py index 92b2451ba..bbe83f525 100644 --- a/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py +++ b/psydac/feec/multipatch/examples/h1_source_pbms_conga_2d.py @@ -23,7 +23,7 @@ from psydac.feec.multipatch.operators import HodgeOperator from psydac.feec.multipatch.plotting_utilities import plot_field from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain -from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution +from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution_OBSOLETE from psydac.feec.multipatch.utilities import time_count from psydac.feec.multipatch.non_matching_operators import construct_V0_conforming_projection, construct_V1_conforming_projection @@ -165,7 +165,7 @@ def lift_u_bc(u_bc): # (not all the returned functions are useful here) N_diag = 200 method = 'conga' - f_scal, f_vect, u_bc, ph_ref, uh_ref, p_ex, u_ex, phi, grad_phi = get_source_and_solution( + f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE( source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name, refsol_params=[N_diag, method, source_proj], ) diff --git a/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py b/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py index 55dd506e6..94f4992ba 100644 --- a/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py +++ b/psydac/feec/multipatch/examples/hcurl_source_pbms_conga_2d.py @@ -24,7 +24,7 @@ from psydac.feec.multipatch.operators import HodgeOperator from psydac.feec.multipatch.plotting_utilities import plot_field from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain -from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution +from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution_OBSOLETE from psydac.feec.multipatch.utilities import time_count from psydac.linalg.utilities import array_to_psydac from psydac.fem.basic import FemField @@ -227,9 +227,8 @@ def lift_u_bc(u_bc): print('getting the source and ref solution...') N_diag = 200 method = 'conga' - f_scal, f_vect, u_bc, ph_ref, uh_ref, p_ex, u_ex, phi, grad_phi = get_source_and_solution( + f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE( source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name, - refsol_params=[N_diag, method, source_proj], ) # compute approximate source f_h diff --git a/psydac/feec/multipatch/examples_nc/h1_source_pbms_nc.py b/psydac/feec/multipatch/examples_nc/h1_source_pbms_nc.py new file mode 100644 index 000000000..9e81eb69d --- /dev/null +++ b/psydac/feec/multipatch/examples_nc/h1_source_pbms_nc.py @@ -0,0 +1,278 @@ +# coding: utf-8 + +from mpi4py import MPI + +import os +import numpy as np +from collections import OrderedDict + +from sympy import lambdify +from scipy.sparse.linalg import spsolve + +from sympde.expr.expr import LinearForm +from sympde.expr.expr import integral, Norm +from sympde.topology import Derham +from sympde.topology import element_of + + +from psydac.api.settings import PSYDAC_BACKENDS +from psydac.feec.multipatch.api import discretize +from psydac.feec.pull_push import pull_2d_h1 + +from psydac.feec.multipatch.fem_linear_operators import IdLinearOperator +from psydac.feec.multipatch.operators import HodgeOperator +from psydac.feec.multipatch.plotting_utilities import plot_field +from psydac.feec.multipatch.multipatch_domain_utilities import build_multipatch_domain +from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution_OBSOLETE +from psydac.feec.multipatch.utilities import time_count +from psydac.feec.multipatch.non_matching_operators import construct_V0_conforming_projection, construct_V1_conforming_projection +from psydac.api.postprocessing import OutputManager, PostProcessManager + +from psydac.linalg.utilities import array_to_psydac +from psydac.fem.basic import FemField + +def solve_h1_source_pbm_nc( + nc=4, deg=4, domain_name='pretzel_f', backend_language=None, source_proj='P_L2', source_type='manu_poisson', + eta=-10., mu=1., gamma_h=10., + plot_source=False, plot_dir=None, hide_plots=True +): + """ + solver for the problem: find u in H^1, such that + + A u = f on \Omega + u = u_bc on \partial \Omega + + where the operator + + A u := eta * u - mu * div grad u + + is discretized as Ah: V0h -> V0h in a broken-FEEC approach involving a discrete sequence on a 2D multipatch domain \Omega, + + V0h --grad-> V1h -—curl-> V2h + + Examples: + + - Helmholtz equation with + eta = -omega**2 + mu = 1 + + - Poisson equation with Laplace operator L = A, + eta = 0 + mu = 1 + + :param nc: nb of cells per dimension, in each patch + :param deg: coordinate degree in each patch + :param gamma_h: jump penalization parameter + :param source_proj: approximation operator for the source, possible values are 'P_geom' or 'P_L2' + :param source_type: must be implemented in get_source_and_solution() + """ + + ncells = nc + degree = [deg,deg] + + # if backend_language is None: + # if domain_name in ['pretzel', 'pretzel_f'] and nc > 8: + # backend_language='numba' + # else: + # backend_language='python' + # print('[note: using '+backend_language+ ' backends in discretize functions]') + + print('---------------------------------------------------------------------------------------------------------') + print('Starting solve_h1_source_pbm function with: ') + print(' ncells = {}'.format(ncells)) + print(' degree = {}'.format(degree)) + print(' domain_name = {}'.format(domain_name)) + print(' source_proj = {}'.format(source_proj)) + print(' backend_language = {}'.format(backend_language)) + print('---------------------------------------------------------------------------------------------------------') + + print('building the multipatch domain...') + domain = build_multipatch_domain(domain_name=domain_name) + mappings = OrderedDict([(P.logical_domain, P.mapping) for P in domain.interior]) + mappings_list = list(mappings.values()) + ncells_h = {patch.name: [ncells[i], ncells[i]] for (i,patch) in enumerate(domain.interior)} + domain_h = discretize(domain, ncells=ncells_h) + + print('building the symbolic and discrete deRham sequences...') + derham = Derham(domain, ["H1", "Hcurl", "L2"]) + derham_h = discretize(derham, domain_h, degree=degree) + + # multi-patch (broken) spaces + V0h = derham_h.V0 + V1h = derham_h.V1 + V2h = derham_h.V2 + print('dim(V0h) = {}'.format(V0h.nbasis)) + print('dim(V1h) = {}'.format(V1h.nbasis)) + print('dim(V2h) = {}'.format(V2h.nbasis)) + + print('broken differential operators...') + # broken (patch-wise) differential operators + bD0, bD1 = derham_h.broken_derivatives_as_operators + bD0_m = bD0.to_sparse_matrix() + # bD1_m = bD1.to_sparse_matrix() + + print('building the discrete operators:') + print('commuting projection operators...') + nquads = [4*(d + 1) for d in degree] + P0, P1, P2 = derham_h.projectors(nquads=nquads) + + I0 = IdLinearOperator(V0h) + I0_m = I0.to_sparse_matrix() + + print('Hodge operators...') + # multi-patch (broken) linear operators / matrices + H0 = HodgeOperator(V0h, domain_h, backend_language=backend_language) + H1 = HodgeOperator(V1h, domain_h, backend_language=backend_language) + + dH0_m = H0.get_dual_Hodge_sparse_matrix() # = mass matrix of V0 + H0_m = H0.to_sparse_matrix() # = inverse mass matrix of V0 + dH1_m = H1.get_dual_Hodge_sparse_matrix() # = mass matrix of V1 + # H1_m = H1.to_sparse_matrix() # = inverse mass matrix of V1 + + print('conforming projection operators...') + # conforming Projections (should take into account the boundary conditions of the continuous deRham sequence) + cP0_m = construct_V0_conforming_projection(V0h, domain_h, hom_bc=True) + # cP1_m = construct_V1_conforming_projection(V1h, domain_h, hom_bc=True) + + if not os.path.exists(plot_dir): + os.makedirs(plot_dir) + + def lift_u_bc(u_bc): + if u_bc is not None: + print('lifting the boundary condition in V0h... [warning: Not Tested Yet!]') + # note: for simplicity we apply the full P1 on u_bc, but we only need to set the boundary dofs + u_bc = lambdify(domain.coordinates, u_bc) + u_bc_log = [pull_2d_h1(u_bc, m.get_callable_mapping()) for m in mappings_list] + # it's a bit weird to apply P1 on the list of (pulled back) logical fields -- why not just apply it on u_bc ? + uh_bc = P0(u_bc_log) + ubc_c = uh_bc.coeffs.toarray() + # removing internal dofs (otherwise ubc_c may already be a very good approximation of uh_c ...) + ubc_c = ubc_c - cP0_m.dot(ubc_c) + else: + ubc_c = None + return ubc_c + + # Conga (projection-based) stiffness matrices: + # div grad: + pre_DG_m = - bD0_m.transpose() @ dH1_m @ bD0_m + + # jump penalization: + jump_penal_m = I0_m - cP0_m + JP0_m = jump_penal_m.transpose() * dH0_m * jump_penal_m + + pre_A_m = cP0_m.transpose() @ ( eta * dH0_m - mu * pre_DG_m ) # useful for the boundary condition (if present) + A_m = pre_A_m @ cP0_m + gamma_h * JP0_m + + print('getting the source and ref solution...') + # (not all the returned functions are useful here) + N_diag = 200 + method = 'conga' + f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE( + source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name, + refsol_params=[N_diag, method, source_proj], + ) + + # compute approximate source f_h + b_c = f_c = None + if source_proj == 'P_geom': + print('projecting the source with commuting projection P0...') + f = lambdify(domain.coordinates, f_scal) + f_log = [pull_2d_h1(f, m.get_callable_mapping()) for m in mappings_list] + f_h = P0(f_log) + f_c = f_h.coeffs.toarray() + b_c = dH0_m.dot(f_c) + + elif source_proj == 'P_L2': + print('projecting the source with L2 projection...') + v = element_of(V0h.symbolic_space, name='v') + expr = f_scal * v + l = LinearForm(v, integral(domain, expr)) + lh = discretize(l, domain_h, V0h) + b = lh.assemble() + b_c = b.toarray() + if plot_source: + f_c = H0_m.dot(b_c) + else: + raise ValueError(source_proj) + + if plot_source: + plot_field(numpy_coeffs=f_c, Vh=V0h, space_kind='h1', domain=domain, title='f_h with P = '+source_proj, filename=plot_dir+'/fh_'+source_proj+'.png', hide_plot=hide_plots) + + ubc_c = lift_u_bc(u_bc) + + if ubc_c is not None: + # modified source for the homogeneous pbm + print('modifying the source with lifted bc solution...') + b_c = b_c - pre_A_m.dot(ubc_c) + + # direct solve with scipy spsolve + print('solving source problem with scipy.spsolve...') + uh_c = spsolve(A_m, b_c) + + # project the homogeneous solution on the conforming problem space + print('projecting the homogeneous solution on the conforming problem space...') + uh_c = cP0_m.dot(uh_c) + + if ubc_c is not None: + # adding the lifted boundary condition + print('adding the lifted boundary condition...') + uh_c += ubc_c + + print('getting and plotting the FEM solution from numpy coefs array...') + title = r'solution $\phi_h$ (amplitude)' + params_str = 'eta={}_mu={}_gamma_h={}'.format(eta, mu, gamma_h) + plot_field(numpy_coeffs=uh_c, Vh=V0h, space_kind='h1', domain=domain, title=title, filename=plot_dir+params_str+'_phi_h.png', hide_plot=hide_plots) + + + if u_ex: + u = element_of(V0h.symbolic_space, name='u') + l2norm = Norm(u - u_ex, domain, kind='l2') + l2norm_h = discretize(l2norm, domain_h, V0h) + uh_c = array_to_psydac(uh_c, V0h.vector_space) + l2_error = l2norm_h.assemble(u=FemField(V0h, coeffs=uh_c)) + return l2_error + +if __name__ == '__main__': + + t_stamp_full = time_count() + + quick_run = True + # quick_run = False + + omega = np.sqrt(170) # source + roundoff = 1e4 + eta = int(-omega**2 * roundoff)/roundoff + # print(eta) + # source_type = 'elliptic_J' + source_type = 'manu_poisson' + + # if quick_run: + # domain_name = 'curved_L_shape' + # nc = 4 + # deg = 2 + # else: + # nc = 8 + # deg = 4 + + domain_name = 'pretzel_f' + # domain_name = 'curved_L_shape' + nc = np.array([8, 8, 16, 16, 8, 4, 4, 4, 4, 4, 2, 2, 4, 16, 16, 8, 2, 2, 2]) + deg = 2 + + # nc = 2 + # deg = 2 + + run_dir = '{}_{}_nc={}_deg={}/'.format(domain_name, source_type, nc, deg) + solve_h1_source_pbm_nc( + nc=nc, deg=deg, + eta=eta, + mu=1, #1, + domain_name=domain_name, + source_type=source_type, + backend_language='pyccel-gcc', + plot_source=True, + plot_dir='./plots/h1_tests_source_february/'+run_dir, + hide_plots=True, + ) + + time_count(t_stamp_full, msg='full program') diff --git a/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_dg.py b/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_dg.py index c2e986a5c..682ff3226 100644 --- a/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_dg.py +++ b/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_dg.py @@ -30,7 +30,7 @@ from psydac.feec.multipatch.non_matching_multipatch_domain_utilities import create_square_domain from psydac.api.postprocessing import OutputManager, PostProcessManager -def hcurl_solve_eigen_pbm_multipatch_dg(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0, +def hcurl_solve_eigen_pbm_dg(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0, generalized_pbm=False, sigma=None, ref_sigmas=[], nb_eigs_solve=8, nb_eigs_plot=5, skip_eigs_threshold=1e-7, plot_dir=None, hide_plots=True, m_load_dir="",): diff --git a/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_nc.py b/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_nc.py index 7efa56b15..a99d601a2 100644 --- a/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_nc.py +++ b/psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_nc.py @@ -33,9 +33,9 @@ from psydac.api.postprocessing import OutputManager, PostProcessManager -def hcurl_solve_eigen_pbm_multipatch_nc(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0, +def hcurl_solve_eigen_pbm_nc(ncells=[[2,2], [2,2]], degree=[3,3], domain=[[0, np.pi],[0, np.pi]], domain_name='refined_square', backend_language='pyccel-gcc', mu=1, nu=0, gamma_h=0, generalized_pbm=False, sigma=None, ref_sigmas=[], nb_eigs_solve=8, nb_eigs_plot=5, skip_eigs_threshold=1e-7, - plot_dir=None, hide_plots=True, m_load_dir="",): + plot_dir=None, hide_plots=True, m_load_dir=None,): diags = {} diff --git a/psydac/feec/multipatch/examples_nc/hcurl_eigen_testcase.py b/psydac/feec/multipatch/examples_nc/hcurl_eigen_testcase.py index e2ea765c9..853d33816 100644 --- a/psydac/feec/multipatch/examples_nc/hcurl_eigen_testcase.py +++ b/psydac/feec/multipatch/examples_nc/hcurl_eigen_testcase.py @@ -2,8 +2,8 @@ import numpy as np -from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_nc import hcurl_solve_eigen_pbm_multipatch_nc -from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_dg import hcurl_solve_eigen_pbm_multipatch_dg +from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_nc import hcurl_solve_eigen_pbm_nc +from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_dg import hcurl_solve_eigen_pbm_dg from psydac.feec.multipatch.utilities import time_count, get_run_dir, get_plot_dir, get_mat_dir, get_sol_dir, diag_fn from psydac.feec.multipatch.utils_conga_2d import write_diags_to_file @@ -235,7 +235,7 @@ # - we look for nb_eigs_solve eigenvalues close to sigma (skip zero eigenvalues if skip_zero_eigs==True) # - we plot nb_eigs_plot eigenvectors if method == 'feec': - diags, eigenvalues = hcurl_solve_eigen_pbm_multipatch_nc( + diags, eigenvalues = hcurl_solve_eigen_pbm_nc( ncells=ncells, degree=degree, gamma_h=gamma_h, generalized_pbm=generalized_pbm, @@ -253,7 +253,7 @@ m_load_dir=m_load_dir, ) elif method == 'dg': - diags, eigenvalues = hcurl_solve_eigen_pbm_multipatch_dg( + diags, eigenvalues = hcurl_solve_eigen_pbm_dg( ncells=ncells, degree=degree, gamma_h=gamma_h, generalized_pbm=generalized_pbm, diff --git a/psydac/feec/multipatch/examples_nc/hcurl_source_pbms_nc.py b/psydac/feec/multipatch/examples_nc/hcurl_source_pbms_nc.py index bbb35ce0f..7f6b314cf 100644 --- a/psydac/feec/multipatch/examples_nc/hcurl_source_pbms_nc.py +++ b/psydac/feec/multipatch/examples_nc/hcurl_source_pbms_nc.py @@ -29,6 +29,7 @@ from psydac.feec.multipatch.utilities import time_count #, export_sol, import_sol from psydac.linalg.utilities import array_to_psydac from psydac.fem.basic import FemField +from psydac.feec.multipatch.examples.ppc_test_cases import get_source_and_solution_OBSOLETE from psydac.feec.multipatch.non_matching_operators import construct_V0_conforming_projection, construct_V1_conforming_projection from psydac.api.postprocessing import OutputManager, PostProcessManager @@ -39,8 +40,8 @@ def solve_hcurl_source_pbm_nc( project_sol=False, plot_source=False, plot_dir=None, hide_plots=True, skip_plot_titles=False, cb_min_sol=None, cb_max_sol=None, - m_load_dir="", sol_filename="", sol_ref_filename="", - ref_nc=None, ref_deg=None, + m_load_dir=None, sol_filename="", sol_ref_filename="", + ref_nc=None, ref_deg=None, test=False ): """ solver for the problem: find u in H(curl), such that @@ -246,9 +247,14 @@ def lift_u_bc(u_bc): t_stamp = time_count(t_stamp) print() print(' -- getting source --') - f_vect, u_bc, u_ex, curl_u_ex, div_u_ex = get_source_and_solution_hcurl( + if source_type == 'manu_maxwell': + f_scal, f_vect, u_bc, p_ex, u_ex, phi, grad_phi = get_source_and_solution_OBSOLETE( source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name, - ) + ) + else: + f_vect, u_bc, u_ex, curl_u_ex, div_u_ex = get_source_and_solution_hcurl( + source_type=source_type, eta=eta, mu=mu, domain=domain, domain_name=domain_name, + ) # compute approximate source f_h t_stamp = time_count(t_stamp) tilde_f_c = f_c = None @@ -352,5 +358,14 @@ def lift_u_bc(u_bc): time_count(t_stamp) + if test: + u = element_of(V1h.symbolic_space, name='u') + l2norm = Norm(Matrix([u[0] - u_ex[0],u[1] - u_ex[1]]), domain, kind='l2') + l2norm_h = discretize(l2norm, domain_h, V1h) + uh_c = array_to_psydac(uh_c, V1h.vector_space) + l2_error = l2norm_h.assemble(u=FemField(V1h, coeffs=uh_c)) + print(l2_error) + return l2_error + return diags \ No newline at end of file diff --git a/psydac/feec/multipatch/examples_nc/hcurl_source_testcase.py b/psydac/feec/multipatch/examples_nc/hcurl_source_testcase.py index fe283a1b7..23ef31fb9 100644 --- a/psydac/feec/multipatch/examples_nc/hcurl_source_testcase.py +++ b/psydac/feec/multipatch/examples_nc/hcurl_source_testcase.py @@ -16,7 +16,6 @@ test_case = 'maxwell_hom_eta=170' # used in paper # test_case = 'maxwell_inhom' # used in paper -compute_ref_sol = False # (not needed for inhomogeneous test-case, as exact solution is known) # # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- diff --git a/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py b/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py index 51746c218..4d88e519f 100644 --- a/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py +++ b/psydac/feec/multipatch/tests/test_feec_maxwell_multipatch_2d.py @@ -3,6 +3,11 @@ import numpy as np from psydac.feec.multipatch.examples.hcurl_source_pbms_conga_2d import solve_hcurl_source_pbm +from psydac.feec.multipatch.examples_nc.hcurl_source_pbms_nc import solve_hcurl_source_pbm_nc + +from psydac.feec.multipatch.examples.hcurl_eigen_pbms_conga_2d import hcurl_solve_eigen_pbm +from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_nc import hcurl_solve_eigen_pbm_nc +from psydac.feec.multipatch.examples_nc.hcurl_eigen_pbms_dg import hcurl_solve_eigen_pbm_dg def test_time_harmonic_maxwell_pretzel_f(): nc,deg = 10,2 @@ -24,6 +29,161 @@ def test_time_harmonic_maxwell_pretzel_f(): assert abs(l2_error - 0.06246693595198972)<1e-10 +def test_time_harmonic_maxwell_pretzel_f_nc(): + deg = 2 + nc = np.array([20, 20, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 10, 10]) + + source_type = 'manu_maxwell' + domain_name = 'pretzel_f' + source_proj = 'tilde_Pi' + + omega = np.sqrt(170) # source + roundoff = 1e4 + eta = int(-omega**2 * roundoff)/roundoff + + l2_error = solve_hcurl_source_pbm_nc( + nc=nc, deg=deg, + eta=eta, + nu=0, + mu=1, + domain_name=domain_name, + source_type=source_type, + source_proj=source_proj, + plot_dir='./plots/th_maxell_nc', + backend_language='pyccel-gcc', + test=True) + + assert abs(l2_error - 0.04753982587323614)<1e-10 + +def test_maxwell_eigen_curved_L_shape(): + domain_name = 'curved_L_shape' + + nc = 10 + deg = 2 + + ref_sigmas = [ + 0.181857115231E+01, + 0.349057623279E+01, + 0.100656015004E+02, + 0.101118862307E+02, + 0.124355372484E+02, + ] + sigma = 7 + nb_eigs_solve = 7 + nb_eigs_plot = 7 + skip_eigs_threshold = 1e-7 + + eigenvalues, eigenvectors = hcurl_solve_eigen_pbm( + nc=nc, deg=deg, + gamma_h=0, + nu=0, + mu=1, + sigma=sigma, + skip_eigs_threshold=skip_eigs_threshold, + nb_eigs=nb_eigs_solve, + nb_eigs_plot=nb_eigs_plot, + domain_name=domain_name, + backend_language='pyccel-gcc', + plot_dir='./plots/eigen_maxell', + ) + + error = 0 + n_errs = min(len(ref_sigmas), len(eigenvalues)) + for k in range(n_errs): + error += (eigenvalues[k]-ref_sigmas[k])**2 + error = np.sqrt(error) + + assert abs(error - 0.023413963252245817)<1e-10 + +def test_maxwell_eigen_curved_L_shape_nc(): + domain_name = 'curved_L_shape' + domain=[[1, 3],[0, np.pi/4]] + + ncells = np.array([[None, 10], + [10, 20]]) + + degree = [2, 2] + + ref_sigmas = [ + 0.181857115231E+01, + 0.349057623279E+01, + 0.100656015004E+02, + 0.101118862307E+02, + 0.124355372484E+02, + ] + sigma = 7 + nb_eigs_solve = 7 + nb_eigs_plot = 7 + skip_eigs_threshold = 1e-7 + + diags, eigenvalues = hcurl_solve_eigen_pbm_nc( + ncells=ncells, degree=degree, + gamma_h=0, + generalized_pbm=True, + nu=0, + mu=1, + sigma=sigma, + ref_sigmas=ref_sigmas, + skip_eigs_threshold=skip_eigs_threshold, + nb_eigs_solve=nb_eigs_solve, + nb_eigs_plot=nb_eigs_plot, + domain_name=domain_name, domain=domain, + backend_language='pyccel-gcc', + plot_dir='./plots/eigen_maxell_nc', + ) + + error = 0 + n_errs = min(len(ref_sigmas), len(eigenvalues)) + for k in range(n_errs): + error += (eigenvalues[k]-ref_sigmas[k])**2 + error = np.sqrt(error) + + assert abs(error - 0.004289103786542442)<1e-10 + +def test_maxwell_eigen_curved_L_shape_dg(): + domain_name = 'curved_L_shape' + domain=[[1, 3],[0, np.pi/4]] + + ncells = np.array([[None, 10], + [10, 20]]) + + degree = [2, 2] + + ref_sigmas = [ + 0.181857115231E+01, + 0.349057623279E+01, + 0.100656015004E+02, + 0.101118862307E+02, + 0.124355372484E+02, + ] + sigma = 7 + nb_eigs_solve = 7 + nb_eigs_plot = 7 + skip_eigs_threshold = 1e-7 + + diags, eigenvalues = hcurl_solve_eigen_pbm_dg( + ncells=ncells, degree=degree, + gamma_h=0, + generalized_pbm=True, + nu=0, + mu=1, + sigma=sigma, + ref_sigmas=ref_sigmas, + skip_eigs_threshold=skip_eigs_threshold, + nb_eigs_solve=nb_eigs_solve, + nb_eigs_plot=nb_eigs_plot, + domain_name=domain_name, domain=domain, + backend_language='pyccel-gcc', + plot_dir='./plots/eigen_maxell_dg', + ) + + error = 0 + n_errs = min(len(ref_sigmas), len(eigenvalues)) + for k in range(n_errs): + error += (eigenvalues[k]-ref_sigmas[k])**2 + error = np.sqrt(error) + + assert abs(error - 0.004208158031148591)<1e-10 #============================================================================== # CLEAN UP SYMPY NAMESPACE @@ -35,4 +195,4 @@ def teardown_module(): def teardown_function(): from sympy.core import cache - cache.clear_cache() + cache.clear_cache() \ No newline at end of file diff --git a/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py b/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py index 7fd6eb040..c18cfe0a2 100644 --- a/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py +++ b/psydac/feec/multipatch/tests/test_feec_poisson_multipatch_2d.py @@ -1,4 +1,7 @@ +import numpy as np + from psydac.feec.multipatch.examples.h1_source_pbms_conga_2d import solve_h1_source_pbm +from psydac.feec.multipatch.examples_nc.h1_source_pbms_nc import solve_h1_source_pbm_nc def test_poisson_pretzel_f(): @@ -16,9 +19,27 @@ def test_poisson_pretzel_f(): backend_language='pyccel-gcc', plot_source=False, plot_dir='./plots/h1_tests_source_february/'+run_dir) - print(l2_error) - assert abs(l2_error-8.054935880021907e-05)<1e-10 + assert abs(l2_error-0.11860734907095004)<1e-10 + +def test_poisson_pretzel_f_nc(): + + source_type = 'manu_poisson_2' + domain_name = 'pretzel_f' + nc = np.array([20, 20, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10, 20, 20, 20, 10, 10]) + deg = 2 + run_dir = '{}_{}_nc={}_deg={}/'.format(domain_name, source_type, nc, deg) + l2_error = solve_h1_source_pbm_nc( + nc=nc, deg=deg, + eta=0, + mu=1, + domain_name=domain_name, + source_type=source_type, + backend_language='pyccel-gcc', + plot_source=False, + plot_dir='./plots/h1_tests_source_february/'+run_dir) + + assert abs(l2_error-0.04324704991715671)<1e-10 #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== @@ -29,7 +50,4 @@ def teardown_module(): def teardown_function(): from sympy.core import cache - cache.clear_cache() - -if __name__ == '__main__': - test_poisson_pretzel_f() + cache.clear_cache() \ No newline at end of file