add time harmonic and time domain maxwell examples

psydac/feec/multipatch/examples/ppc_test_cases.py

@@ -90,15 +90,15 @@ def get_Gaussian_beam(x_0, y_0, domain=None):
     x = x - x_0
     y = y - y_0
 
-    k = 2*pi 
-    sigma = 0.7
+    k = pi 
+    sigma = 0.5
 
     xy = x**2 + y**2
     ef = exp( - xy/(2*sigma**2) )
 
     E = cos(k * y) * ef
-    B = y/(sigma**2) * E - sin(k * y) * ef 
-
+    B = -y/(sigma**2) * E 
+    
     return Tuple(E, 0), B
 
 def get_Gaussian_beam2(x_0, y_0, domain=None):

psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_dg.py

@@ -148,7 +148,10 @@ def hcurl_solve_eigen_pbm_multipatch_dg(ncells=[[2,2], [2,2]], degree=[3,3], dom
 
     t_stamp = time_count(t_stamp)
     print('plotting the eigenmodes...')     
-
+
+    if not os.path.exists(plot_dir):
+        os.makedirs(plot_dir)
+
    # OM = OutputManager('spaces.yml', 'fields.h5')
    # OM.add_spaces(V1h=V1h)
 

psydac/feec/multipatch/examples_nc/hcurl_eigen_pbms_nc.py

@@ -32,25 +32,6 @@
 
 from psydac.api.postprocessing import OutputManager, PostProcessManager
 
-#from said
-from scipy.sparse.linalg import spsolve, inv
-
-from sympde.calculus      import grad, dot, curl, cross
-from sympde.calculus      import minus, plus
-from sympde.topology      import VectorFunctionSpace
-from sympde.topology      import elements_of
-from sympde.topology      import NormalVector
-from sympde.topology      import Square
-from sympde.topology      import IdentityMapping, PolarMapping
-from sympde.expr.expr     import LinearForm, BilinearForm
-from sympde.expr.expr     import integral
-from sympde.expr.expr     import Norm
-from sympde.expr.equation import find, EssentialBC
-
-from psydac.api.tests.build_domain   import build_pretzel
-from psydac.fem.basic                import FemField
-from psydac.api.settings             import PSYDAC_BACKEND_GPYCCEL
-from psydac.feec.pull_push           import pull_2d_hcurl
 
 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,
                           generalized_pbm=False, sigma=None, ref_sigmas=[], nb_eigs_solve=8, nb_eigs_plot=5, skip_eigs_threshold=1e-7,
@@ -253,112 +234,7 @@ def hcurl_solve_eigen_pbm_multipatch_nc(ncells=[[2,2], [2,2]], degree=[3,3], dom
 
         t_stamp = time_count(t_stamp)
 
-    ### Saids Code
-
-    V  = VectorFunctionSpace('V', domain, kind='hcurl')
-
-    u, v, F  = elements_of(V, names='u, v, F')
-    nn       = NormalVector('nn')
-
-    I        = domain.interfaces
-    boundary = domain.boundary
-
-    kappa   = 10
-    k       = 1
-
-    jump = lambda w:plus(w)-minus(w)
-    avr  = lambda w:0.5*plus(w) + 0.5*minus(w)
-
-    expr1_I  =  cross(nn, jump(v))*curl(avr(u))\
-               +k*cross(nn, jump(u))*curl(avr(v))\
-               +kappa*cross(nn, jump(u))*cross(nn, jump(v))
-
-    expr1   = curl(u)*curl(v) 
-    expr1_b = -cross(nn, v) * curl(u) -k*cross(nn, u)*curl(v)  + kappa*cross(nn, u)*cross(nn, v)
-    ## curl curl u = - omega**2 u 
-
-    expr2   = dot(u,v)
-    #expr2_I = kappa*cross(nn, jump(u))*cross(nn, jump(v))
-    #expr2_b = -k*cross(nn, u)*curl(v) + kappa * cross(nn, u) * cross(nn, v)
-
-    # Bilinear form a: V x V --> R
-    a      = BilinearForm((u,v),  integral(domain, expr1) + integral(I, expr1_I) + integral(boundary, expr1_b))
-
-    # Linear form l: V --> R
-    b     = BilinearForm((u,v), integral(domain, expr2))# + integral(I, expr2_I) + integral(boundary, expr2_b))
-
-    #+++++++++++++++++++++++++++++++
-    # 2. Discretization
-    #+++++++++++++++++++++++++++++++
-
-    domain_h = discretize(domain, ncells=ncells_h)
-    Vh       = discretize(V, domain_h, degree=degree)
-
-    ah = discretize(a, domain_h, [Vh, Vh])
-    Ah_m = ah.assemble().tosparse()
-
-    bh = discretize(b, domain_h, [Vh, Vh])
-    Bh_m = bh.assemble().tosparse()
-
-    all_eigenvalues_2, all_eigenvectors_transp_2 = get_eigenvalues(nb_eigs_solve, sigma, Ah_m, Bh_m)
-
-    #Eigenvalue processing
-    t_stamp = time_count(t_stamp)
-    print('sorting out eigenvalues...')
-    zero_eigenvalues2 = []
-    if skip_eigs_threshold is not None:
-        eigenvalues2 = []
-        eigenvectors2 = []
-        for val, vect in zip(all_eigenvalues_2, all_eigenvectors_transp_2.T):
-            if abs(val) < skip_eigs_threshold: 
-                zero_eigenvalues2.append(val)
-                # we skip the eigenvector
-            else:
-                eigenvalues2.append(val)
-                eigenvectors2.append(vect)
-    else:
-        eigenvalues2 = all_eigenvalues_2
-        eigenvectors2 = all_eigenvectors_transp_2.T
-    diags['DG'] = True 
-    for k, val in enumerate(eigenvalues2):
-        diags['eigenvalue2_{}'.format(k)] = val #eigenvalues[k]
-
-    for k, val in enumerate(zero_eigenvalues2):
-        diags['skipped eigenvalue2_{}'.format(k)] = val
-
-    t_stamp = time_count(t_stamp)
-    print('plotting the eigenmodes...')     
-
-   # OM = OutputManager('spaces.yml', 'fields.h5')
-   # OM.add_spaces(V1h=V1h)
-
-    nb_eigs = len(eigenvalues2)
-    for i in range(min(nb_eigs_plot, nb_eigs)):
-        OM = OutputManager(plot_dir+'/spaces2.yml', plot_dir+'/fields2.h5')
-        OM.add_spaces(V1h=Vh)
-        print('looking at emode i = {}... '.format(i))
-        lambda_i  = eigenvalues2[i]
-        emode_i = np.real(eigenvectors2[i])
-        norm_emode_i = np.dot(emode_i,Bh_m.dot(emode_i))
-        eh_c = emode_i/norm_emode_i 
-        stencil_coeffs = array_to_psydac(eh_c, Vh.vector_space)
-        vh = FemField(Vh, coeffs=stencil_coeffs)
-        OM.set_static()
-        #OM.add_snapshot(t=i , ts=0) 
-        OM.export_fields(vh = vh)
-
-        #print('norm of computed eigenmode: ', norm_emode_i)
-        # plot the broken eigenmode:
-        OM.export_space_info()
-        OM.close()
-
-        PM = PostProcessManager(domain=domain, space_file=plot_dir+'/spaces2.yml', fields_file=plot_dir+'/fields2.h5' )
-        PM.export_to_vtk(plot_dir+"/eigen2_{}".format(i),grid=None, npts_per_cell=[6]*2,snapshots='all', fields='vh' )
-        PM.close()
-
-        t_stamp = time_count(t_stamp)
-
-    return diags, eigenvalues, eigenvalues2
+    return diags, eigenvalues
 
 
 def get_eigenvalues(nb_eigs, sigma, A_m, M_m):

psydac/feec/multipatch/examples_nc/hcurl_eigen_testcase.py

@@ -1,9 +1,9 @@
 import os
 import numpy as np
 
-#from psydac.feec.multipatch.examples_nc.multipatch_non_conf_examples import hcurl_solve_eigen_pbm_multipatch_nc
-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_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.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
@@ -15,6 +15,8 @@
 # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 
 #
 # test-case and numerical parameters:
+# method = 'feec'
+method = 'dg'
 
 operator = 'curl-curl' 
 degree = [3,3] # shared across all patches
@@ -69,8 +71,8 @@
 domain=[[1, 3],[0, np.pi/4]] # interval in x- and y-direction
 
 
-ncells = np.array([[None, 10],
-                   [10,   5]])
+ncells = np.array([[None, 5],
+                   [5,   10]])
 
 
 
@@ -111,7 +113,7 @@
 else:
     raise ValueError(operator)
 
-case_dir = 'talk_eigenpbm_'+operator
+case_dir = 'eigenpbm_'+operator+'_'+method
 ref_case_dir = case_dir
 
 cb_min_sol = None
@@ -232,8 +234,26 @@
 # note:
 #   - we look for nb_eigs_solve eigenvalues close to sigma (skip zero eigenvalues if skip_zero_eigs==True)
 #   - we plot nb_eigs_plot eigenvectors
-
-diags, eigenvalues, eigenvalues2 = hcurl_solve_eigen_pbm_multipatch_nc(
+if method == 'feec':
+    diags, eigenvalues = hcurl_solve_eigen_pbm_multipatch_nc(
+        ncells=ncells, degree=degree,
+        gamma_h=gamma_h,
+        generalized_pbm=generalized_pbm,
+        nu=nu,
+        mu=mu,
+        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=backend_language,
+        plot_dir=plot_dir,
+        hide_plots=True,
+        m_load_dir=m_load_dir,
+    )
+elif method == 'dg': 
+    diags, eigenvalues = hcurl_solve_eigen_pbm_multipatch_dg(
     ncells=ncells, degree=degree,
     gamma_h=gamma_h,
     generalized_pbm=generalized_pbm,
@@ -249,14 +269,14 @@
     plot_dir=plot_dir,
     hide_plots=True,
     m_load_dir=m_load_dir,
-)
+    )
+
 
 if ref_sigmas is not None:
         errors = []
         n_errs = min(len(ref_sigmas), len(eigenvalues))
         for k in range(n_errs):
             diags['error_{}'.format(k)] = abs(eigenvalues[k]-ref_sigmas[k])
-            diags['error2_{}'.format(k)] = abs(eigenvalues2[k]-ref_sigmas[k])
 #
 # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----