+
+ # <(I-P)v,(I-P)w> in rhs
+ generalized_pbm = True
+
+ # curl-curl operator
+ nu = 0
+ mu = 1
+
+ # reference eigenvalues for validation
+ if domain_name == 'refined_square':
+ assert domain == [[0, np.pi], [0, np.pi]]
+ ref_sigmas = [
+ 1, 1,
+ 2,
+ 4, 4,
+ 5, 5,
+ 8,
+ 9, 9,
+ ]
+ sigma = 5
+ nb_eigs_solve = 10
+ nb_eigs_plot = 10
+ skip_eigs_threshold = 1e-7
+
+ elif domain_name == 'curved_L_shape':
+ # ref eigenvalues from Monique Dauge benchmark page
+ assert domain == [[1, 3], [0, np.pi / 4]]
+ ref_sigmas = [
+ 0.181857115231E+01,
+ 0.349057623279E+01,
+ 0.100656015004E+02,
+ 0.101118862307E+02,
+ 0.124355372484E+02,
+ ]
+ sigma = 7
+ nb_eigs_solve = 5
+ nb_eigs_plot = 5
+ skip_eigs_threshold = 1e-7
+
+ eigenvalues = hcurl_solve_eigen_pbm_dg(
+ ncells=ncells, degree=degree,
+ nu=nu,
+ mu=mu,
+ sigma=sigma,
+ skip_eigs_threshold=skip_eigs_threshold,
+ nb_eigs_solve=nb_eigs_solve,
+ nb_eigs_plot=nb_eigs_plot,
+ domain_name=domain_name, domain=domain,
+ )
+
+ if ref_sigmas is not None:
+ n_errs = min(len(ref_sigmas), len(eigenvalues))
+ for k in range(n_errs):
+ print('error_{}: '.format(k), abs(eigenvalues[k] - ref_sigmas[k]))
diff --git a/examples/feec/hcurl_source_pbms_conga_2d.py b/examples/feec/hcurl_source_pbms_conga_2d.py
new file mode 100644
index 000000000..cfc31e331
--- /dev/null
+++ b/examples/feec/hcurl_source_pbms_conga_2d.py
@@ -0,0 +1,355 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+"""
+ solver for the problem: find u in H(curl), such that
+
+ A u = f on \\Omega
+ n x u = n x u_bc on \\partial \\Omega
+
+ where the operator
+
+ A u := eta * u + mu * curl curl u - nu * grad div u
+
+ is discretized as Ah: V1h -> V1h in a broken-FEEC approach involving a discrete sequence on a 2D multipatch domain \\Omega,
+
+ V0h --grad-> V1h -—curl-> V2h
+"""
+
+import os
+import numpy as np
+
+from sympde.topology import Derham
+
+from psydac.api.discretization import discretize
+from psydac.api.postprocessing import OutputManager, PostProcessManager
+
+from psydac.feec.multipatch_domain_utilities import build_multipatch_domain
+
+from psydac.fem.basic import FemField
+from psydac.fem.projectors import get_dual_dofs
+
+from psydac.linalg.basic import IdentityOperator
+from psydac.linalg.solvers import inverse
+
+#==============================================================================
+# Solver for H(curl) source problems
+#==============================================================================
+def solve_hcurl_source_pbm(
+ nc=4, deg=4, domain_name='pretzel_f', backend_language=None, source_type='manu_maxwell_inhom',
+ eta=-10., mu=1., nu=0., gamma_h=10.,
+ project_sol=True, plot_dir=None):
+ """
+ solver for the problem: find u in H(curl), such that
+
+ A u = f on \\Omega
+ n x u = n x u_bc on \\partial \\Omega
+
+ where the operator
+
+ A u := eta * u + mu * curl curl u - nu * grad div u
+
+ is discretized as Ah: V1h -> V1h in a broken-FEEC approach involving a discrete sequence on a 2D multipatch domain \\Omega,
+
+ V0h --grad-> V1h -—curl-> V2h
+
+ Examples:
+
+ - time-harmonic maxwell equation with
+ eta = -omega**2
+ mu = 1
+ nu = 0
+
+ - Hodge-Laplacian operator L = A with
+ eta = 0
+ mu = 1
+ nu = 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_type: must be implemented in get_source_and_solution()
+ """
+ degree = [deg, deg]
+
+
+ print('---------------------------------------------------------------------------------------------------------')
+ print('Starting solve_hcurl_source_pbm function with: ')
+ print(' ncells = {}'.format(nc))
+ print(' degree = {}'.format(degree))
+ print(' domain_name = {}'.format(domain_name))
+ print(' backend_language = {}'.format(backend_language))
+ print('---------------------------------------------------------------------------------------------------------')
+
+ print()
+ print(' -- building discrete spaces and operators --')
+
+ print(' .. multi-patch domain...')
+ domain = build_multipatch_domain(domain_name=domain_name)
+
+ if isinstance(nc, int):
+ ncells = [nc, nc]
+ else:
+ ncells = {patch.name: [nc[i], nc[i]]
+ for (i, patch) in enumerate(domain.interior)}
+
+ print(' .. derham sequence...')
+ derham = Derham(domain, ["H1", "Hcurl", "L2"])
+
+ print(' .. discrete domain...')
+ domain_h = discretize(domain, ncells=ncells)
+
+ print(' .. discrete derham sequence...')
+ derham_h = discretize(derham, domain_h, degree=degree)
+
+ print(' .. commuting projection operators...')
+ nquads = [10 * (d + 1) for d in degree]
+ P0, P1, P2 = derham_h.projectors(nquads=nquads)
+
+ print(' .. multi-patch spaces...')
+ V0h, V1h, V2h = derham_h.spaces
+ mappings = derham_h.callable_mapping
+
+ print('dim(V0h) = {}'.format(V0h.nbasis))
+ print('dim(V1h) = {}'.format(V1h.nbasis))
+ print('dim(V2h) = {}'.format(V2h.nbasis))
+
+
+ print(' .. Id operator and matrix...')
+ I1 = IdentityOperator(V1h.coeff_space)
+
+ print(' .. Hodge operators...')
+ # multi-patch (broken) linear operators / matrices
+ # other option: define as Hodge Operators:
+ H0, H1, H2 = derham_h.hodge_operators(kind='linop', backend_language=backend_language)
+ dH0, dH1, dH2 = derham_h.hodge_operators(kind='linop', dual=True, backend_language=backend_language)
+
+ print(' .. conforming Projection operators...')
+ # conforming Projections (should take into account the boundary conditions
+ # of the continuous deRham sequence)
+ cP0, cP1, cP2 = derham_h.conforming_projectors(kind='linop', hom_bc = True)
+
+ print(' .. broken differential operators...')
+ # broken (patch-wise) differential operators
+ bD0, bD1 = derham_h.derivatives(kind='linop')
+
+ # Conga (projection-based) stiffness matrices
+ # curl curl:
+ print(' .. curl-curl stiffness matrix...')
+ pre_CC = bD1.T @ H2 @ bD1
+
+ # grad div:
+ print(' .. grad-div stiffness matrix...')
+ pre_GD = - H1 @ bD0 @ cP0 @ dH0 @ cP0.T @ bD0.T @ H1
+
+ # jump stabilization:
+ print(' .. jump stabilization matrix...')
+ JS = (I1 - cP1).T @ H1 @ (I1 - cP1)
+
+ print(' .. full operator matrix...')
+ print('eta = {}'.format(eta))
+ print('mu = {}'.format(mu))
+ print('nu = {}'.format(nu))
+ print('STABILIZATION: gamma_h = {}'.format(gamma_h))
+ # useful for the boundary condition (if present)
+ pre_A = eta * cP1.T @ H1
+ if mu != 0:
+ pre_A += mu * cP1.T @ pre_CC
+ if nu != 0:
+ pre_A -= nu * cP1.T @ pre_GD
+
+ A = pre_A @ cP1 + gamma_h * JS
+
+ print()
+ print(' -- getting source --')
+
+ f_vect, u_bc, 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
+ # f_h = L2 projection of f_vect, with filtering if tilde_Pi
+ tilde_f = get_dual_dofs(Vh=V1h, f=f_vect, domain_h=domain_h, backend_language=backend_language)
+
+ print(' .. filtering the discrete source with P1.T ...')
+ tilde_f = cP1.T @ tilde_f
+
+ def lift_u_bc(u_bc):
+ if u_bc is not None:
+ ubc = P1(u_bc).coeffs
+ ubc -= cP1.dot(ubc)
+
+ else:
+ ubc = None
+
+ return ubc
+
+ ubc = lift_u_bc(u_bc)
+
+ if ubc is not None:
+ # modified source for the homogeneous pbm
+ print(' .. modifying the source with lifted bc solution...')
+ tilde_f -= pre_A.dot(ubc)
+
+ # direct solve with scipy spsolve
+ print('solving source problem with conjugate gradient...')
+ solver = inverse(A, solver='cg', tol=1e-8)
+ u = solver.solve(tilde_f)
+
+ # project the homogeneous solution on the conforming problem space
+ if project_sol:
+ print(' .. projecting the homogeneous solution on the conforming problem space...')
+ u = cP1.dot(u)
+
+ if ubc is not None:
+ # adding the lifted boundary condition
+ print(' .. adding the lifted boundary condition...')
+ u += ubc
+
+ uh = FemField(V1h, coeffs=u)
+ f = dH1.dot(tilde_f)
+ jh = FemField(V1h, coeffs=f)
+
+ print(' -- plots and diagnostics --')
+ if plot_dir:
+ if not os.path.exists(plot_dir):
+ os.makedirs(plot_dir)
+
+ OM = OutputManager(plot_dir + '/spaces.yml', plot_dir + '/fields.h5')
+ OM.add_spaces(V1h=V1h)
+ OM.set_static()
+ OM.export_fields(vh=uh)
+ OM.export_fields(jh=jh)
+ OM.export_space_info()
+ OM.close()
+
+ PM = PostProcessManager(
+ domain=domain,
+ space_file=plot_dir +
+ '/spaces.yml',
+ fields_file=plot_dir +
+ '/fields.h5')
+ PM.export_to_vtk(
+ plot_dir + "/sol",
+ grid=None,
+ npts_per_cell=[6] * 2,
+ snapshots='all',
+ fields='vh')
+ PM.export_to_vtk(
+ plot_dir + "/source",
+ grid=None,
+ npts_per_cell=[6] * 2,
+ snapshots='all',
+ fields='jh')
+
+ PM.close()
+
+
+ if u_ex:
+ u_ex_p = P1(u_ex).coeffs
+
+ err = u_ex_p - u
+ l2_error = np.sqrt( H1.dot_inner(err, err) / H1.dot_inner(u_ex_p, u_ex_p))
+ print("L2 error: ", l2_error)
+
+ return l2_error
+
+#==============================================================================
+# Test sources and exact solutions
+#==============================================================================
+def get_source_and_solution_hcurl(
+ source_type=None, eta=0, mu=0, nu=0,
+ domain=None, domain_name=None):
+ """
+ provide source, and exact solutions when available, for:
+
+ Find u in H(curl) such that
+
+ A u = f on \\Omega
+ n x u = n x u_bc on \\partial \\Omega
+
+ with
+
+ A u := eta * u + mu * curl curl u - nu * grad div u
+
+ see solve_hcurl_source_pbm()
+ """
+ from sympy import pi, cos, sin, Tuple, exp
+
+ # exact solutions (if available)
+ u_ex = None
+
+ # bc solution: describe the bc on boundary. Inside domain, values should
+ # not matter. Homogeneous bc will be used if None
+ u_bc = None
+
+ # source terms
+ f_vect = None
+
+ # auxiliary term (for more diagnostics)
+ grad_phi = None
+ phi = None
+
+ x, y = domain.coordinates
+
+ if source_type == 'manu_maxwell_inhom':
+ # used for Maxwell equation with manufactured solution
+ f_vect = Tuple(eta * sin(pi * y) - pi**2 * sin(pi * y) * cos(pi * x) + pi**2 * sin(pi * y),
+ eta * sin(pi * x) * cos(pi * y) + pi**2 * sin(pi * x) * cos(pi * y))
+ if nu == 0:
+ u_ex = Tuple(sin(pi * y), sin(pi * x) * cos(pi * y))
+ curl_u_ex = pi * (cos(pi * x) * cos(pi * y) - cos(pi * y))
+ div_u_ex = -pi * sin(pi * x) * sin(pi * y)
+ else:
+ raise NotImplementedError
+ u_bc = u_ex
+
+ elif source_type == 'elliptic_J':
+ # no manufactured solution for Maxwell pbm
+ x0 = 1.5
+ y0 = 1.5
+ s = (x - x0) - (y - y0)
+ t = (x - x0) + (y - y0)
+ a = (1 / 1.9)**2
+ b = (1 / 1.2)**2
+ sigma2 = 0.0121
+ tau = a * s**2 + b * t**2 - 1
+ phi = exp(-tau**2 / (2 * sigma2))
+ dx_tau = 2 * (a * s + b * t)
+ dy_tau = 2 * (-a * s + b * t)
+
+ f_x = dy_tau * phi
+ f_y = - dx_tau * phi
+ f_vect = Tuple(f_x, f_y)
+
+ else:
+ raise ValueError(source_type)
+
+ from sympy import lambdify
+ u_bc_x = lambdify(domain.coordinates, u_bc[0])
+ u_bc_y = lambdify(domain.coordinates, u_bc[1])
+
+ u_ex_x = lambdify(domain.coordinates, u_ex[0])
+ u_ex_y = lambdify(domain.coordinates, u_ex[1])
+
+ return f_vect, [u_bc_x, u_bc_y], [u_ex_x, u_ex_y]
+
+if __name__ == '__main__':
+ nc = 5
+ deg = 3
+
+ source_type = 'manu_maxwell_inhom'
+ domain_name = 'pretzel_f'
+
+ omega = np.pi
+ eta = -omega**2 # source
+
+ err = solve_hcurl_source_pbm(
+ nc=nc, deg=deg,
+ eta=eta,
+ nu=0,
+ mu=1,
+ domain_name=domain_name,
+ source_type=source_type,
+ backend_language='pyccel-gcc')
+
\ No newline at end of file
diff --git a/examples/feec/tests/__init__.py b/examples/feec/tests/__init__.py
new file mode 100644
index 000000000..419109b64
--- /dev/null
+++ b/examples/feec/tests/__init__.py
@@ -0,0 +1,5 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
diff --git a/examples/feec/timedomain_maxwell.py b/examples/feec/timedomain_maxwell.py
new file mode 100644
index 000000000..658dc8bfe
--- /dev/null
+++ b/examples/feec/timedomain_maxwell.py
@@ -0,0 +1,569 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+"""
+ solver for the TD Maxwell problem: find E(t) in H(curl), B in L2, such that
+
+ dt E - curl B = -J on \\Omega
+ dt B + curl E = 0 on \\Omega
+ n x E = n x E_bc on \\partial \\Omega
+
+ with Ampere discretized weakly and Faraday discretized strongly, in a broken-FEEC approach on a 2D multipatch domain \\Omega,
+
+ V0h --grad-> V1h -—curl-> V2h
+ (Eh) (Bh)
+"""
+import os
+import numpy as np
+
+from sympde.calculus import grad, dot, curl, cross
+from sympde.topology import NormalVector
+from sympde.topology import elements_of
+from sympde.topology import Derham
+from sympde.expr.expr import integral
+from sympde.expr.expr import BilinearForm
+
+from psydac.linalg.basic import IdentityOperator
+
+from psydac.api.settings import PSYDAC_BACKENDS
+from psydac.api.discretization import discretize
+from psydac.api.postprocessing import OutputManager, PostProcessManager
+
+from psydac.feec.multipatch_domain_utilities import build_multipatch_domain, build_cartesian_multipatch_domain
+
+from psydac.fem.basic import FemField
+from psydac.fem.projectors import get_dual_dofs
+
+#==============================================================================
+# Solver for the TD Maxwell problem
+#==============================================================================
+def solve_td_maxwell_pbm(*,
+ nc=4,
+ deg=4,
+ final_time=20,
+ cfl_max=0.8,
+ dt_max=None,
+ domain_name='pretzel_f',
+ backend='pyccel-gcc',
+ source_type='zero',
+ E0_type='pulse_2',
+ plot_dir=None,
+ domain_lims=None,
+ p_moments=-1,
+ ):
+ """
+ solver for the TD Maxwell problem: find E(t) in H(curl), B in L2, such that
+
+ dt E - curl B = -J on \\Omega
+ dt B + curl E = 0 on \\Omega
+ n x E = n x E_bc on \\partial \\Omega
+
+ with Ampere discretized weakly and Faraday discretized strongly, in a broken-FEEC approach on a 2D multipatch domain \\Omega,
+
+ V0h --grad-> V1h -—curl-> V2h
+ (Eh) (Bh)
+
+ Parameters
+ ----------
+ nc : int
+ Number of cells (same along each direction) in every patch.
+
+ deg : int
+ Polynomial degree (same along each direction) in every patch, for the
+ spline space V0 in H1.
+
+ final_time : float
+ Final simulation time. Given that the speed of light is set to c=1,
+ this can be easily chosen based on the wave transit time in the domain.
+
+ cfl_max : float
+ Maximum Courant parameter in the simulation domain, used to determine
+ the time step size.
+
+ dt_max : float
+ Maximum time step size, which has to be met together with cfl_max. This
+ additional constraint is useful to resolve a time-dependent source.
+
+ domain_name : str
+ Name of the multipatch geometry used in the simulation, to be chosen
+ among those available in the function `build_multipatch_domain`.
+
+ backend : str
+ Name of the backend used for acceleration of the computational kernels,
+ to be chosen among the available keys of the PSYDAC_BACKENDS dict.
+
+ source_type : str {'zero' | 'pulse' | 'cf_pulse' }
+ Name that identifies the space-time profile of the current source, to be
+ chosen among those available in the function get_source_and_solution().
+ Available options:
+ - 'zero' : no current source
+ - 'pulse' : div-free current source, time-harmonic
+ - 'cf_pulse': curl-free current source, time-harmonic
+
+ E0_type : str {'zero', 'pulse'}
+ Initial conditions for the electric field. Choose 'zero' for E0=0
+ and 'pulse' for a non-zero field localized in a small region.
+
+ plot_dir : str
+ Path to the directory where the figures will be saved.
+
+ domain_lims : list
+ If the domain_name is 'refined_square' or 'square_L_shape', this
+ parameter must be set to the list of the two intervals defining the
+ rectangular domain, i.e. `[[x_min, x_max], [y_min, y_max]]`.
+
+ p_moments : int
+ Degree of the polynomial moments used in the conforming projection.
+ """
+ degree = [deg, deg]
+
+
+ print('---------------------------------------------------------------------------------------------------------')
+ print('Starting solve_td_maxwell_pbm function with: ')
+ print(' ncells = {}'.format(nc))
+ print(' degree = {}'.format(degree))
+ print(' domain_name = {}'.format(domain_name))
+ print(' backend = {}'.format(backend))
+ print('---------------------------------------------------------------------------------------------------------')
+
+
+ print()
+ print(' -- building discrete spaces and operators --')
+
+ print(' .. multi-patch domain...')
+ if domain_name == 'refined_square' or domain_name == 'square_L_shape':
+ int_x, int_y = domain_lims
+ domain = build_cartesian_multipatch_domain(nc, int_x, int_y, mapping='identity')
+ else:
+ domain = build_multipatch_domain(domain_name=domain_name)
+
+ if isinstance(nc, int):
+ ncells = [nc, nc]
+ elif nc.ndim == 1:
+ ncells = {patch.name: [nc[i], nc[i]]
+ for (i, patch) in enumerate(domain.interior)}
+ elif nc.ndim == 2:
+ ncells = {patch.name: [nc[int(patch.name[2])][int(patch.name[4])],
+ nc[int(patch.name[2])][int(patch.name[4])]] for patch in domain.interior}
+
+
+ print(' .. derham sequence...')
+ derham = Derham(domain, ["H1", "Hcurl", "L2"])
+
+ print(' .. discrete domain...')
+ domain_h = discretize(domain, ncells=ncells)
+
+ print(' .. discrete derham sequence...')
+ derham_h = discretize(derham, domain_h, degree=degree)
+
+ print(' .. commuting projection operators...')
+ nquads = [4 * (d + 1) for d in degree]
+ P0, P1, P2 = derham_h.projectors(nquads=nquads)
+
+ print(' .. multi-patch spaces...')
+ V0h, V1h, V2h = derham_h.spaces
+
+ print(' .. Id operator and matrix...')
+ I1 = IdentityOperator(V1h.coeff_space)
+
+ print(' .. Hodge operators...')
+ H0, H1, H2 = derham_h.hodge_operators(kind='linop')
+ dH0, dH1, dH2 = derham_h.hodge_operators(kind='linop', dual=True)
+
+ print(' .. conforming Projection operators...')
+ cP0, cP1, cP2 = derham_h.conforming_projectors(kind='linop', p_moments = p_moments, hom_bc = False)
+
+ print(' .. broken differential operators...')
+ bD0, bD1 = derham_h.derivatives(kind='linop')
+
+ print(' .. matrix of the primal curl (in primal bases)...')
+ C = bD1 @ cP1
+
+ print(' .. matrix of the dual curl (also in primal bases)...')
+ dC = dH1 @ C.T @ H2
+
+ ### Silvermueller ABC
+ u, v = elements_of(derham.V1, names='u, v')
+ nn = NormalVector('nn')
+ boundary = domain.boundary
+ expr_b = cross(nn, u) * cross(nn, v)
+
+ a = BilinearForm((u, v), integral(boundary, expr_b))
+ ah = discretize(a, domain_h, [V1h, V1h], backend=PSYDAC_BACKENDS[backend],)
+ A_eps = ah.assemble()
+
+ # Compute stable time step size based on max CFL and max dt
+ dt = compute_stable_dt(C=C, dC=dC, cfl_max=cfl_max, dt_max=dt_max)
+ print(" Reduce time step to match the simulation final time:")
+ Nt = int(np.ceil(final_time / dt))
+ dt = final_time / Nt
+ print(f" . Time step size : dt = {dt}")
+ print(' total nb of time steps: Nt = {}, final time: T = {:5.4f}'.format(Nt, final_time))
+
+ H1A = H1 + dt * A_eps
+
+ # alternative inverse
+ # from psydac.linalg.solvers import inverse
+ # H1A_inv = inverse(H1A, solver='cg', tol=1e-8)
+ ###
+ M = H1A
+ from scipy.linalg import inv
+ from scipy.sparse import csr_matrix
+ from psydac.linalg.sparse import SparseMatrixLinearOperator
+ M_inv = inv(M.toarray())
+ M_inv = csr_matrix(M_inv)
+ H1A_inv = SparseMatrixLinearOperator(M.codomain, M.domain, M_inv)
+ ####
+
+ # Absorbing dC
+ dC = H1A_inv @ C.T @ H2
+ dCH1 = H1A_inv @ H1
+
+
+ # ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
+ # source
+ print()
+ print(' -- getting source --')
+
+ if source_type == 'zero':
+
+ f0 = None
+ f0_harmonic = None
+
+ elif source_type == 'pulse':
+
+ f0 = get_div_free_pulse(x_0=np.pi/2, y_0=np.pi/2, domain=domain)
+
+ elif source_type == 'cf_pulse':
+
+ f0 = get_curl_free_pulse(x_0=np.pi/2, y_0=np.pi/2, domain=domain)
+
+ else:
+
+ raise ValueError(source_type)
+
+
+ if f0 is not None:
+ print(' .. projecting the source f0 with L2 projection...')
+ tilde_f0_h = get_dual_dofs(Vh=V1h, f=f0, domain_h=domain_h, backend_language=backend)
+
+ print(' .. filtering the source...')
+ tilde_f0_h = cP1.T @ tilde_f0_h
+
+ f0_h = dH1.dot(tilde_f0_h)
+
+ else:
+
+ f0_h = V1h.coeff_space.zeros()
+
+ # ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
+ # initial solution
+
+ print(' -- initial solution --')
+
+ # initial B sol
+ B_h = V2h.coeff_space.zeros()
+ E_h = V1h.coeff_space.zeros()
+
+ # initial E sol
+ if E0_type == 'zero':
+ E_h = V1h.coeff_space.zeros()
+
+ elif E0_type == 'pulse':
+
+ E0 = get_div_free_pulse(x_0=np.pi/2, y_0=np.pi/2, domain=domain)
+
+ print(' .. projecting E0 with L2 projection...')
+ tilde_E0_h = get_dual_dofs(Vh=V1h, f=E0, domain_h=domain_h, backend_language=backend)
+ E_h = dH1.dot(tilde_E0_h)
+
+ elif E0_type == 'pulse_2':
+
+ E0, B0 = get_Gaussian_beam(y_0=np.pi/2, x_0=np.pi/2, domain=domain)
+
+ print(' .. projecting E0 with L2 projection...')
+ tilde_E0_h = get_dual_dofs(Vh=V1h, f=E0, domain_h=domain_h, backend_language=backend)
+ E_h = dH1.dot(tilde_E0_h)
+
+ tilde_B0_h = get_dual_dofs(Vh=V2h, f=B0, domain_h=domain_h, backend_language=backend)
+ B_h = dH2.dot(tilde_B0_h)
+
+ elif E0_type == 'Gaussian':
+
+ E0, B0 = get_Gaussian_beam(y_0=np.pi/2, x_0=np.pi/2, domain=domain)
+
+ print(' .. projecting E0 with L2 projection...')
+ tilde_E0_h = get_dual_dofs(Vh=V1h, f=E0, domain_h=domain_h, backend_language=backend)
+ E_h = dH1.dot(tilde_E0_h)
+
+ tilde_B0_h = get_dual_dofs(Vh=V2h, f=B0, domain_h=domain_h, backend_language=backend)
+ B_h = dH2.dot(tilde_B0_h)
+
+ else:
+ raise ValueError(E0_type)
+
+ # ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- -----
+ # time loop
+
+
+ if plot_dir is not None and not os.path.exists(plot_dir):
+ os.makedirs(plot_dir)
+
+ if plot_dir:
+ OM1 = OutputManager(plot_dir + '/spaces1.yml', plot_dir + '/fields1.h5')
+ OM1.add_spaces(V1h=V1h)
+ OM1.export_space_info()
+
+ OM2 = OutputManager(plot_dir + '/spaces2.yml', plot_dir + '/fields2.h5')
+ OM2.add_spaces(V2h=V2h)
+ OM2.export_space_info()
+
+ Eh = FemField(V1h, coeffs=cP1 @ E_h)
+ OM1.add_snapshot(t=0, ts=0)
+ OM1.export_fields(Eh=Eh)
+
+ Bh = FemField(V2h, coeffs=B_h)
+ OM2.add_snapshot(t=0, ts=0)
+ OM2.export_fields(Bh=Bh)
+
+
+ f_h = f0_h.copy()
+ Btemp_h = B_h.copy()
+ Etemp_h = E_h.copy()
+
+ print(" -- time loop --")
+ for nt in range(Nt):
+ print(' .. nt+1 = {}/{}'.format(nt+1, Nt))
+
+ # 1/2 faraday: Bn -> Bn+1/2
+ # B_h -= (dt/2) * C @ E_h
+ # E_h = A_eps @ E_h + dt * dC @ B_h
+ # B_h -= (dt/2) * C @ E_h
+
+ C.dot(E_h, out=Btemp_h)
+ B_h -= (dt/2) * Btemp_h
+
+ dCH1.dot(E_h, out=E_h)
+ dC.dot(B_h, out=Etemp_h)
+ E_h += dt * (Etemp_h - f_h)
+
+ C.dot(E_h, out=Btemp_h)
+ B_h -= (dt/2) * Btemp_h
+
+ if plot_dir:
+ Eh = FemField(V1h, coeffs=cP1 @ E_h)
+ OM1.add_snapshot(t=nt*dt, ts=nt)
+ OM1.export_fields(Eh = Eh)
+
+ Bh = FemField(V2h, coeffs=B_h)
+ OM2.add_snapshot(t=nt*dt, ts=nt)
+ OM2.export_fields(Bh=Bh)
+
+ if plot_dir:
+ OM1.close()
+
+ print("Post process fields")
+ PM = PostProcessManager(
+ domain=domain,
+ space_file=plot_dir + '/spaces1.yml',
+ fields_file=plot_dir + '/fields1.h5')
+ PM.export_to_vtk(
+ plot_dir + "/Eh",
+ grid=None,
+ npts_per_cell=4,
+ snapshots='all',
+ fields='Eh')
+ PM.close()
+
+ PM = PostProcessManager(
+ domain=domain,
+ space_file=plot_dir + '/spaces2.yml',
+ fields_file=plot_dir + '/fields2.h5')
+ PM.export_to_vtk(
+ plot_dir + "/Bh",
+ grid=None,
+ npts_per_cell=4,
+ snapshots='all',
+ fields='Bh')
+ PM.close()
+
+# ==============================================================================
+# Compute stable time step size
+# ==============================================================================
+def compute_stable_dt(*, C, dC, cfl_max, dt_max=None):
+ """
+ Compute a stable time step size based on the maximum CFL parameter in the
+ domain. To this end we estimate the operator norm of
+
+ `dC @ C: V1h -> V1h`,
+
+ find the largest stable time step compatible with Strang splitting, and
+ rescale it by the provided `cfl_max`. Setting `cfl_max = 1` would run the
+ scheme exactly at its stability limit, which is not safe because of the
+ unavoidable round-off errors. Hence we require `0 < cfl_max < 1`.
+
+ Optionally the user can provide a maximum time step size in order to
+ properly resolve some time scales of interest (e.g. a time-dependent
+ current source).
+
+ Parameters
+ ----------
+ C : LinearOperator
+ Matrix of the Curl operator.
+
+ dC : LinearOperator
+ Matrix of the dual Curl operator.
+
+ cfl_max : float
+ Maximum Courant parameter in the domain, intended as a stability
+ parameter (=1 at the stability limit). Must be `0 < cfl_max < 1`.
+
+ dt_max : float, optional
+ If not None, restrict the computed dt by this value in order to
+ properly resolve time scales of interest. Must be > 0.
+
+ Returns
+ -------
+ dt : float
+ Largest stable dt which satisfies the provided constraints.
+
+ """
+
+ print(" .. compute_stable_dt by estimating the operator norm of ")
+ print(" .. dC_m @ C_m: V1h -> V1h ")
+ print(" .. with dim(V1h) = {} ...".format(C.domain.dimension))
+
+ if not (0 < cfl_max < 1):
+ print(' ****** ****** ****** ****** ****** ****** ')
+ print(' WARNING !!! cfl = {} '.format(cfl))
+ print(' ****** ****** ****** ****** ****** ****** ')
+
+ V = C.domain
+ from psydac.linalg.utilities import array_to_psydac
+ vv = array_to_psydac(np.random.rand(V.dimension), V)
+
+ norm_vv = np.sqrt(vv.inner(vv))
+
+ max_ncfl = 500
+ ncfl = 0
+ spectral_rho = 1
+ conv = False
+ CC = dC @ C
+
+ while not (conv or ncfl > max_ncfl):
+
+ vv *= (1. / norm_vv)
+ ncfl += 1
+ CC.dot(vv, out=vv)
+
+ norm_vv = np.sqrt(vv.inner(vv))
+ old_spectral_rho = spectral_rho
+ spectral_rho = norm_vv # approximation
+ conv = abs((spectral_rho - old_spectral_rho) / spectral_rho) < 0.001
+ print(" ... spectral radius iteration: spectral_rho( dC @ C ) ~= {}".format(spectral_rho))
+
+ norm_op = np.sqrt(spectral_rho)
+ c_dt_max = 2. / norm_op
+
+ light_c = 1
+ dt = cfl_max * c_dt_max / light_c
+
+ if dt_max is not None:
+ dt = min(dt, dt_max)
+
+ print(" Time step dt computed for Maxwell solver:")
+ print(f" Based on cfl_max = {cfl_max} and dt_max = {dt_max}, we set dt = {dt}")
+ print(f" -- note that c*Dt = {light_c*dt} and c_dt_max = {c_dt_max}, thus c * dt / c_dt_max = {light_c*dt/c_dt_max}")
+ print(f" -- and spectral_radius((c*dt)**2* dC @ C ) = {(light_c * dt * norm_op)**2} (should be < 4).")
+
+ return dt
+
+# ==============================================================================
+# Test Sources
+# ==============================================================================
+def get_div_free_pulse(x_0, y_0, domain=None):
+
+ from sympy import pi, cos, sin, Tuple, exp
+
+ x, y = domain.coordinates
+ ds2_0 = (0.02)**2
+ sigma_0 = (x - x_0)**2 + (y - y_0)**2
+ phi_0 = exp(-sigma_0**2 / (2 * ds2_0))
+ dx_sig_0 = 2 * (x - x_0)
+ dy_sig_0 = 2 * (y - y_0)
+ dx_phi_0 = - dx_sig_0 * sigma_0 / ds2_0 * phi_0
+ dy_phi_0 = - dy_sig_0 * sigma_0 / ds2_0 * phi_0
+ f_x = dy_phi_0
+ f_y = - dx_phi_0
+ f_vect = Tuple(f_x, f_y)
+
+ return f_vect
+
+
+def get_curl_free_pulse(x_0, y_0, domain=None, pp=False):
+
+ from sympy import pi, cos, sin, Tuple, exp
+
+ # return -grad phi_0
+ x, y = domain.coordinates
+ if pp:
+ # psi=phi
+ ds2_0 = (0.02)**2
+ else:
+ ds2_0 = (0.1)**2
+ sigma_0 = (x - x_0)**2 + (y - y_0)**2
+ phi_0 = exp(-sigma_0**2 / (2 * ds2_0))
+ dx_sig_0 = 2 * (x - x_0)
+ dy_sig_0 = 2 * (y - y_0)
+ dx_phi_0 = - dx_sig_0 * sigma_0 / ds2_0 * phi_0
+ dy_phi_0 = - dy_sig_0 * sigma_0 / ds2_0 * phi_0
+ f_x = -dx_phi_0
+ f_y = -dy_phi_0
+ f_vect = Tuple(f_x, f_y)
+
+ return f_vect
+
+def get_Gaussian_beam(x_0, y_0, domain=None):
+
+ from sympy import pi, cos, sin, Tuple, exp
+
+ # return E = cos(k*x) exp( - x^2 + y^2 / 2 sigma^2) v
+ x, y = domain.coordinates
+
+ x = x - x_0
+ y = y - y_0
+
+ sigma = 0.1
+
+ xy = x**2 + y**2
+ ef = 1 / (sigma**2) * exp(- xy / (2 * sigma**2))
+
+ # E = curl exp
+ E = Tuple(y * ef, -x * ef)
+
+ # B = curl E
+ B = (xy / (sigma**2) - 2) * ef
+
+ return E, B
+
+if __name__ == '__main__':
+ domain_name = 'refined_square'
+ domain_lims = [[0, np.pi], [0, np.pi]]
+
+ nc = 20
+ ncells = np.array([[nc, nc, nc],
+ [nc, 2*nc, nc],
+ [nc, nc, nc]])
+
+ deg = 3
+ p_moments = deg+1
+
+ final_time = 2
+
+ plot_dir = './td_maxwell_pulse/'
+
+ solve_td_maxwell_pbm(nc=ncells, deg=deg, p_moments=p_moments, final_time=final_time,
+ domain_name=domain_name, domain_lims=domain_lims,
+ source_type='zero', E0_type='pulse', plot_dir=plot_dir)
diff --git a/examples/notebooks/feec_curlcurl_eigenvalue.ipynb b/examples/notebooks/feec_curlcurl_eigenvalue.ipynb
new file mode 100644
index 000000000..d77f7ef9e
--- /dev/null
+++ b/examples/notebooks/feec_curlcurl_eigenvalue.ipynb
@@ -0,0 +1,342 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "6be9f483",
+ "metadata": {},
+ "source": [
+ "# Example of solving the curl-curl eigenvalue problem\n",
+ "\n",
+ "In this notebook we show how to compute approximate eigenvalues for the curl-curl eigenvalue problem using the FEEC (Finite Element Exterior Calculus) API of PSYDAC. The problem reads as follows:\n",
+ "\n",
+ "\n",
+ "Let $\\Omega=[0, \\pi] \\times [0, \\pi]$. We want to find $\\lambda >0$ and $\\boldsymbol{E} \\in H_0(\\mathrm{curl}, \\Omega)$ such that\n",
+ "\n",
+ "$$ \\mathbf{curl} \\ \\mathrm{curl} \\ \\boldsymbol{E} = \\lambda \\boldsymbol{E}. $$\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Step 1 : Building the domain\n",
+ "\n",
+ "This is similar to what is done in other examples. We leave the option to use a single patch or multiple patches.\n",
+ "The FEEC interface is compatible with both. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c9d51f8e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "\n",
+ "single_patch = False\n",
+ "\n",
+ "# degree of the splines in each direction\n",
+ "degree = (3, 3)\n",
+ "\n",
+ "if single_patch:\n",
+ " #number of cells in each direction\n",
+ " ncells = [10, 10]\n",
+ "\n",
+ " from sympde.topology import Square, IdentityMapping\n",
+ " logical_domain = Square('Omega', bounds1=(0, np.pi), bounds2=(0, np.pi))\n",
+ " mapping = IdentityMapping('M1', dim=2)\n",
+ " domain = mapping(logical_domain)\n",
+ "\n",
+ "else:\n",
+ " # We can give the shape of the multipatch domain \n",
+ " # through the shape of the array containing the number of cells in each direction\n",
+ " ncells = np.array([[5, 5], \n",
+ " [5, 5]])\n",
+ " # The multipatch domain then has the following structure:\n",
+ " # A | B\n",
+ " # -----\n",
+ " # C | D\n",
+ "\n",
+ " from psydac.feec.multipatch_domain_utilities import build_cartesian_multipatch_domain\n",
+ " domain = build_cartesian_multipatch_domain(ncells, (0, np.pi), (0, np.pi), mapping='identity')\n",
+ "\n",
+ " # We now need to convert the ncells array into a dictionary, \n",
+ " # where the keys are the patches\n",
+ " ncells = {patch.name: [ncells[int(patch.name[2])][int(patch.name[4])], \n",
+ " ncells[int(patch.name[2])][int(patch.name[4])]] for patch in domain.interior}\n",
+ "\n",
+ "\n",
+ "from sympde.utilities.utils import plot_domain\n",
+ "plot_domain(domain, isolines=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6bcee6a7",
+ "metadata": {},
+ "source": [
+ "## Step 2: Create and discretize the de Rham complex"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "dc4a3e7a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.api.discretization import discretize\n",
+ "from sympde.topology import Derham\n",
+ "\n",
+ "domain_h = discretize(domain, ncells=ncells) \n",
+ "\n",
+ "derham = Derham(domain, [\"H1\", \"Hcurl\", \"L2\"])\n",
+ "derham_h = discretize(derham, domain_h, degree=degree)\n",
+ "\n",
+ "V0h, V1h, V2h = derham_h.spaces"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "540b8e43",
+ "metadata": {},
+ "source": [
+ "## Step 3: Use the FEEC interface to get the operators\n",
+ "\n",
+ "The (broken) FEEC operators are given by: \n",
+ "The conforming projections $P^\\ell: V^\\ell_\\mathrm{pw} \\rightarrow V^\\ell_h$, mapping from broken to continuous spaces. The broken derivatives $D^\\ell: V^\\ell_\\mathrm{pw} \\rightarrow V^{\\ell + 1}_\\mathrm{pw}$ and the Hodge-operators $H^\\ell$ which can be seen as mass matrices. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "18134cef",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# conforming projectors with homogeneous boundary conditions\n",
+ "# (In the single patch case, they only involve the zero BCs)\n",
+ "cP0, cP1, cP2 = derham_h.conforming_projectors(kind='linop', hom_bc = True)\n",
+ "\n",
+ "# broken (patch-wise) derivative operators\n",
+ "bD0, bD1 = derham_h.derivatives(kind='linop')\n",
+ "\n",
+ "# broken Hodge operators (mass matrices) and their inverses\n",
+ "H0, H1, H2 = derham_h.hodge_operators(kind='linop')\n",
+ "dH0, dH1, dH2 = derham_h.hodge_operators(kind='linop', dual=True)\n",
+ "\n",
+ "from psydac.linalg.basic import IdentityOperator\n",
+ "I1 = IdentityOperator(V1h.coeff_space)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "21c450cf",
+ "metadata": {},
+ "source": [
+ "## Step 4: Set up the system matrices\n",
+ "\n",
+ "The $\\mathrm{curl}$ operators are discretized by the conforming curl operator $D^1 P^1$ and its transpose. The right-hand-side is given by the conforming mass matrix $(P^1)^T H^1 P^1$ where we also have to add a penalization term corresponding to the discontinuous jumps in the broken spaces: $(I - P^1)^T H^1 (I - P^1)$. In total, the discrete equation reads: \n",
+ "\n",
+ "$$ \\left( D^1 P^1 \\right)^T H^2 D^1 P^1 \\boldsymbol{E}_h = \\lambda \\left((P^1)^T H^1 P^1 + (I - P^1)^T H^1 (I - P^1)\\right) \\boldsymbol{E}_h$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4d3c55c1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "CC = cP1.T @ bD1.T @ H2 @ bD1 @ cP1\n",
+ "RHS = cP1.T @ H1 @ cP1 + (I1 - cP1).T @ H1 @ (I1 - cP1)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "17f466e6",
+ "metadata": {},
+ "source": [
+ "## Step 5: Scipy eigenvalue solver"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f02eab35",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_eigenvalues(nb_eigs, sigma, A_m, M_m):\n",
+ " \"\"\"\n",
+ " Compute the eigenvalues of the matrix A close to sigma and right-hand-side M\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " nb_eigs : int\n",
+ " Number of eigenvalues to compute\n",
+ " sigma : float\n",
+ " Value close to which the eigenvalues are computed\n",
+ " A_m : sparse matrix\n",
+ " Matrix A\n",
+ " M_m : sparse matrix\n",
+ " Matrix M\n",
+ " \"\"\"\n",
+ "\n",
+ " from scipy.sparse.linalg import spilu, lgmres\n",
+ " from scipy.sparse.linalg import LinearOperator, eigsh, minres\n",
+ " from scipy.linalg import norm\n",
+ "\n",
+ " print('----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ----- ')\n",
+ " print(\n",
+ " 'computing {0} eigenvalues (and eigenvectors) close to sigma={1} with scipy.sparse.eigsh...'.format(\n",
+ " nb_eigs,\n",
+ " sigma))\n",
+ " mode = 'normal'\n",
+ " which = 'LM'\n",
+ " # from eigsh docstring:\n",
+ " # ncv = number of Lanczos vectors generated ncv must be greater than k and smaller than n;\n",
+ " # it is recommended that ncv > 2*k. Default: min(n, max(2*k + 1, 20))\n",
+ " ncv = 4 * nb_eigs\n",
+ " try_lgmres = True\n",
+ " max_shape_splu = 24000 # OK for nc=20, deg=6 on pretzel_f\n",
+ " if A_m.shape[0] < max_shape_splu:\n",
+ " print('(via sparse LU decomposition)')\n",
+ " OPinv = None\n",
+ " tol_eigsh = 0\n",
+ " else:\n",
+ "\n",
+ " OP_m = A_m - sigma * M_m\n",
+ " tol_eigsh = 1e-7\n",
+ " if try_lgmres:\n",
+ " print(\n",
+ " '(via SPILU-preconditioned LGMRES iterative solver for A_m - sigma*M1_m)')\n",
+ " OP_spilu = spilu(OP_m, fill_factor=15, drop_tol=5e-5)\n",
+ " preconditioner = LinearOperator(\n",
+ " OP_m.shape, lambda x: OP_spilu.solve(x))\n",
+ " tol = tol_eigsh\n",
+ " OPinv = LinearOperator(\n",
+ " matvec=lambda v: lgmres(OP_m, v, x0=None, tol=tol, atol=tol, M=preconditioner,\n",
+ " callback=lambda x: print(\n",
+ " 'cg -- residual = ', norm(OP_m.dot(x) - v))\n",
+ " )[0],\n",
+ " shape=M_m.shape,\n",
+ " dtype=M_m.dtype\n",
+ " )\n",
+ "\n",
+ " else:\n",
+ " # from https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigsh.html:\n",
+ " # the user can supply the matrix or operator OPinv, which gives x = OPinv @ b = [A - sigma * M]^-1 @ b.\n",
+ " # > here, minres: MINimum RESidual iteration to solve Ax=b\n",
+ " # suggested in https://github.com/scipy/scipy/issues/4170\n",
+ " print('(with minres iterative solver for A_m - sigma*M1_m)')\n",
+ " OPinv = LinearOperator(\n",
+ " matvec=lambda v: minres(\n",
+ " OP_m,\n",
+ " v,\n",
+ " tol=1e-10)[0],\n",
+ " shape=M_m.shape,\n",
+ " dtype=M_m.dtype)\n",
+ "\n",
+ " eigenvalues, eigenvectors = eigsh(\n",
+ " A_m, k=nb_eigs, M=M_m, sigma=sigma, mode=mode, which=which, ncv=ncv, tol=tol_eigsh, OPinv=OPinv)\n",
+ "\n",
+ " print(\"done: eigenvalues found: \" + repr(eigenvalues))\n",
+ " return eigenvalues, eigenvectors"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c21fcd3e",
+ "metadata": {},
+ "source": [
+ "## Step 6: Solve the system\n",
+ "\n",
+ "We solve for the first `nb_eigs` many eigenvalue closest to the reference value `sigma`.\n",
+ "In the case `nb_eigs=10` and `sigma=5`, the exact eigenvalues are the first 10 sum of squares, i.e.\n",
+ "`ref_sigmas = [1, 1, 2, 4, 4, 5, 5, 8, 9, 9]`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5fd5ab7b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ref_sigmas = [\n",
+ " 1, 1,\n",
+ " 2,\n",
+ " 4, 4,\n",
+ " 5, 5,\n",
+ " 8,\n",
+ " 9, 9,\n",
+ " ]\n",
+ "sigma = 5\n",
+ "nb_eigs = 10\n",
+ "\n",
+ "eigenvalues, eigenvectors_transp = get_eigenvalues(nb_eigs, sigma, CC.tosparse(), RHS.tosparse())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "f2b1242e",
+ "metadata": {},
+ "source": [
+ "## Step 7: Plot the eigenvalues"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5f452f48",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "fig,ax = plt.subplots(1,1)\n",
+ "ax.set_title(f\"The first {nb_eigs} eigenvalues\")\n",
+ "im = ax.scatter(np.arange(eigenvalues.size), eigenvalues)\n",
+ "plt.xlabel('index');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a978966d",
+ "metadata": {},
+ "source": [
+ "### Testing the notebook"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7f4f0c33",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import ipytest\n",
+ "ipytest.autoconfig(raise_on_error=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2e38bb10",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%ipytest\n",
+ "\n",
+ "l2_error = np.sqrt(np.sum((eigenvalues - ref_sigmas)**2)) \n",
+ "\n",
+ "def test_l2error():\n",
+ " assert l2_error < 5e-04"
+ ]
+ }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/notebooks/feec_time_harmonic_Maxwell.ipynb b/examples/notebooks/feec_time_harmonic_Maxwell.ipynb
new file mode 100644
index 000000000..cd4913bd9
--- /dev/null
+++ b/examples/notebooks/feec_time_harmonic_Maxwell.ipynb
@@ -0,0 +1,326 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "31838f69",
+ "metadata": {},
+ "source": [
+ "# Solving the time-harmonic Maxwells equations on a multipatch domain\n",
+ "\n",
+ "For $\\omega \\in \\mathbb{R}$ and $\\boldsymbol{J} \\in L^2(\\Omega)$, we solve\n",
+ "$$- \\omega^2 \\boldsymbol{E} + \\mathbf{curl} \\mathrm{curl} \\boldsymbol{E} = \\boldsymbol{J}, $$\n",
+ "where $\\boldsymbol{E} \\in H(\\mathrm{curl}, \\Omega).$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d7cf3108",
+ "metadata": {},
+ "source": [
+ "## Step 1: Discretize the domain. \n",
+ "\n",
+ "For the domain $\\Omega$, we choose the pretzel domain:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c18b57f5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.feec.multipatch_domain_utilities import build_multipatch_domain\n",
+ "\n",
+ "domain = build_multipatch_domain(domain_name='pretzel_f')\n",
+ "\n",
+ "from sympde.utilities.utils import plot_domain\n",
+ "plot_domain(domain, isolines=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ee044477",
+ "metadata": {},
+ "source": [
+ "## Step 2: Discretize the de Rham sequence"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "63def2e4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.api.discretization import discretize\n",
+ "from sympde.topology import Derham\n",
+ "\n",
+ "degree = [2, 2]\n",
+ "# For simplicity, we use the same number of cells for all patches\n",
+ "ncells = [4, 4]\n",
+ "\n",
+ "derham = Derham(domain, [\"H1\", \"Hcurl\", \"L2\"])\n",
+ "\n",
+ "domain_h = discretize(domain, ncells=ncells)\n",
+ "derham_h = discretize(derham, domain_h, degree=degree)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d451401e",
+ "metadata": {},
+ "source": [
+ "## Step 3: Obtain the FEEC operators"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d783864a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Finite Element spaces\n",
+ "V0h, V1h, V2h = derham_h.spaces\n",
+ "\n",
+ "# Geometric global projectors\n",
+ "nquads = [(d + 1) for d in degree]\n",
+ "P0, P1, P2 = derham_h.projectors(nquads=nquads)\n",
+ "\n",
+ "# Identity operator\n",
+ "from psydac.linalg.basic import IdentityOperator\n",
+ "I1 = IdentityOperator(V1h.coeff_space)\n",
+ "\n",
+ "# Hodge and dual Hodge operators\n",
+ "H0, H1, H2 = derham_h.hodge_operators(kind='linop')\n",
+ "dH0, dH1, dH2 = derham_h.hodge_operators(kind='linop', dual=True)\n",
+ "\n",
+ "# Conforming projectors with homogeneous boundary conditions\n",
+ "cP0, cP1, cP2 = derham_h.conforming_projectors(kind='linop', hom_bc=True)\n",
+ "\n",
+ "# Broken derivative operators\n",
+ "bD0, bD1 = derham_h.derivatives(kind='linop')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "aa3148da",
+ "metadata": {},
+ "source": [
+ "## Step 4: Assemble the system matrices"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f03ccb61",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "omega = np.pi\n",
+ "\n",
+ "CC = cP1.T @ bD1.T @ H2 @ bD1 \n",
+ "M = -omega**2 * cP1.T @ H1 \n",
+ "pre_A = M + CC\n",
+ "\n",
+ "JS = 10 * (I1 - cP1).T @ H1 @ (I1 - cP1)\n",
+ "\n",
+ "A = pre_A @ cP1 + JS"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7e61a0da",
+ "metadata": {},
+ "source": [
+ "## Step 5: Source and exact solution"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2aee193c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_source_and_solution_hcurl(eta=0, domain=None):\n",
+ " from sympy import pi, cos, sin, Tuple, exp\n",
+ "\n",
+ " x, y = domain.coordinates\n",
+ "\n",
+ " # used for Maxwell equation with manufactured solution\n",
+ " f_vect = Tuple(eta * sin(pi * y) - pi**2 * sin(pi * y) * cos(pi * x) + pi**2 * sin(pi * y),\n",
+ " eta * sin(pi * x) * cos(pi * y) + pi**2 * sin(pi * x) * cos(pi * y))\n",
+ "\n",
+ " u_ex = Tuple(sin(pi * y), sin(pi * x) * cos(pi * y))\n",
+ "\n",
+ "\n",
+ " from sympy import lambdify\n",
+ " u_ex_x = lambdify(domain.coordinates, u_ex[0])\n",
+ " u_ex_y = lambdify(domain.coordinates, u_ex[1])\n",
+ "\n",
+ " return f_vect, [u_ex_x, u_ex_y]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6959cf37",
+ "metadata": {},
+ "source": [
+ "## Step 6: Get and project source"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6f41db72",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "f_vect, u_ex = get_source_and_solution_hcurl(eta=-omega**2, domain=domain)\n",
+ "\n",
+ "from psydac.fem.projectors import get_dual_dofs\n",
+ "tilde_f = get_dual_dofs(Vh=V1h, f=f_vect, domain_h=domain_h)\n",
+ "f = dH1.dot(tilde_f)\n",
+ "\n",
+ "# filtering the discrete source \n",
+ "tilde_f = cP1.T @ tilde_f"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "63fc7530",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.fem.plotting_utilities import plot_field_2d as plot_field\n",
+ "plot_field(stencil_coeffs=f, Vh=V1h, space_kind='hcurl', domain=domain, title='tilde f_h: Geometric projection of f', hide_plot=False)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a4a2d3f6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Manipulate the source to account for inhomogeneous boundary conditions in u\n",
+ "ubc = P1(u_ex).coeffs\n",
+ "ubc -= cP1 @ ubc\n",
+ "\n",
+ "tilde_f -= pre_A @ ubc"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "25860d86",
+ "metadata": {},
+ "source": [
+ "## Step 7: Solve the linear system using CG"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3e83202a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.linalg.solvers import inverse\n",
+ "\n",
+ "solver = inverse(A, solver='cg', tol=1e-5)\n",
+ "u = solver.solve(tilde_f)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "33576e21",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Project back to the full space and add the boundary conditions\n",
+ "u = cP1.dot(u)\n",
+ "u += ubc"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "90124e42",
+ "metadata": {},
+ "source": [
+ "## Step 8: Look at the exact solution and compute the error"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2b31445c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "uex = P1(u_ex).coeffs\n",
+ "err = uex - u\n",
+ "l2_error = np.sqrt(H1.dot_inner(err, err) / H1.dot_inner(uex, uex))\n",
+ "print('L2 relative error = ', l2_error)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ab5d7502",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_field(stencil_coeffs=u, Vh=V1h, space_kind='hcurl', domain=domain, title='Numerical solution u', hide_plot=False)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "44ab3744",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_field(stencil_coeffs=uex, Vh=V1h, space_kind='hcurl', domain=domain, title='Exact solution u_ex', hide_plot=False)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c89c62c3",
+ "metadata": {},
+ "source": [
+ "### Test the notebook"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "81370f97",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import ipytest\n",
+ "ipytest.autoconfig(raise_on_error=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e47d2bb5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%ipytest\n",
+ "\n",
+ "def test_l2error():\n",
+ " assert l2_error < 8e-03"
+ ]
+ }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/notebooks/feec_vector_potential_torus.ipynb b/examples/notebooks/feec_vector_potential_torus.ipynb
new file mode 100644
index 000000000..9cf5e26a2
--- /dev/null
+++ b/examples/notebooks/feec_vector_potential_torus.ipynb
@@ -0,0 +1,534 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "0c5b846b",
+ "metadata": {},
+ "source": [
+ "# Computing the vector potential of a solenoidal field on a Torus\n",
+ "A torus is a 3D domain with zero cavities and one tunnel. \n",
+ "\n",
+ "As such, on a Torus, there exists only the trivial normal harmonic field, but a one-dimensional space of tangential harmonic fields.\n",
+ "\n",
+ "Hence, any solenoidal field $B\\in H_0(div0)$ on a torus can be decomposed (Helmholtz decomposition) as\n",
+ "\n",
+ "$B = curl(A_0) + \\lambda B_H$\n",
+ "\n",
+ "where $B_H$ is a tangential harmonic field, and $A_0\\in H_0(curl)$ is the vector potential of the \"harmonic-free\" solenoidal field $B - \\lambda B_H$.\n",
+ "\n",
+ "In this notebook, we compute a vector potential of the solenoidal magnetic field\n",
+ "\n",
+ "$B(x, y, z) = \\left(\\begin{matrix} \n",
+ "(1/(x^2 + y^2)) \\cdot ( (-C_0 y) + x(2z-1)(\\sqrt{x^2 + y^2}-r)(\\sqrt{x^2 + y^2}-R) ) \\\\ \n",
+ "(1/(x^2 + y^2)) \\cdot ( ( C_0 x) + y(2z-1)(\\sqrt{x^2 + y^2}-r)(\\sqrt{x^2 + y^2}-R) ) \\\\ \n",
+ "(1/(x^2 + y^2)) \\cdot ( \\sqrt{x^2 + y^2}((r+R)-2\\sqrt{x^2 + y^2})z(z-1) ) \n",
+ "\\end{matrix}\\right),\\quad\\quad C_0=-10,\\ r=0.5,\\ R=1.$\n",
+ "\n",
+ "where $r$ is the inner and $R$ is the outer radius. We use this magnetic field as initial guess for a relaxation method for the computation of an MHD equilibrium.\n",
+ "\n",
+ "We compute $A_0$ and a vector potential of the tangential harmonic field $B_H$ separately."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "6c5d3a3e",
+ "metadata": {},
+ "source": [
+ "## Step 1 : Computation of the \"harmonic-free\" vector potential $A_0$\n",
+ "S.th. $curl(A_0) = B - \\lambda B_H$.\n",
+ "\n",
+ "A (weak-divergence-free) vector potential $A_0\\in H_0(curl)$ of the \"harmonic-free\" part of $B$ can be computed as the solution of the Hodge-Laplace problem\n",
+ "\n",
+ "$\\left( \\widetilde{curl}_h\\ curl\\ - grad\\ \\widetilde{div}_h \\right) A_0 = \\widetilde{curl}_h\\ B.$\n",
+ "\n",
+ "### Define and discretize the domain and the de Rham sequence"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9a7d4d77",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from mpi4py import MPI\n",
+ "import numpy as np\n",
+ "\n",
+ "from sympde.topology import Cube, Mapping, Derham\n",
+ "\n",
+ "from psydac.api.discretization import discretize\n",
+ "\n",
+ "comm = MPI.COMM_WORLD # MPI communicator\n",
+ "ncells = [16, 16, 16] # number of cells in each direction\n",
+ "degree = [3, 3, 3] # B-spline degree in each direction\n",
+ "periodic = [False, True, False] # periodicity of the domain\n",
+ "\n",
+ "r = 0.5 # inner radius of the torus\n",
+ "R = 1. # outer radius of the torus\n",
+ "\n",
+ "logical_domain = Cube('C', bounds1=(0.5,1), bounds2=(0,2*np.pi), bounds3=(0,1)) # 3D cubic domain - to be mapped to a toroidal domain (with square cross section)\n",
+ "\n",
+ "class SquareTorus(Mapping):\n",
+ "\n",
+ " _expressions = {'x': 'x1 * cos(x2)',\n",
+ " 'y': 'x1 * sin(x2)',\n",
+ " 'z': 'x3'}\n",
+ " \n",
+ " _ldim = 3\n",
+ " _pdim = 3\n",
+ "\n",
+ "mapping = SquareTorus('ST')\n",
+ "\n",
+ "domain = mapping(logical_domain) # the mapped domain\n",
+ "derham = Derham(domain)\n",
+ "domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm)\n",
+ "derham_h = discretize(derham, domain_h, degree=degree)\n",
+ "\n",
+ "V0, V1, V2, V3 = derham.spaces\n",
+ "V0h, V1h, V2h, V3h = derham_h.spaces\n",
+ "V0cs, V1cs, V2cs, V3cs = [Vh.coeff_space for Vh in derham_h.spaces]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d3cd1ca2",
+ "metadata": {},
+ "source": [
+ "### Assembly of the system matrix\n",
+ "\n",
+ "corresponding to the bilinear form, and additionally taking care of homogeneous boundary conditions,\n",
+ "\n",
+ "$V^1_h\\times V^1_h \\ni (u, v) \\mapsto \\int curl(u)\\cdot curl(v) + \\widetilde{div}_h(u) \\widetilde{div}_h(v)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ef34b831",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.calculus import inner\n",
+ "from sympde.expr import integral, BilinearForm\n",
+ "from sympde.topology import elements_of\n",
+ "\n",
+ "from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL\n",
+ "from psydac.linalg.basic import IdentityOperator\n",
+ "from psydac.linalg.solvers import inverse\n",
+ "\n",
+ "backend = PSYDAC_BACKEND_GPYCCEL\n",
+ "\n",
+ "G, C, D = derham_h.derivatives(kind='linop') # gradient, curl and divergence as psydac.linalg.basic.LinearOperator objects\n",
+ "\n",
+ "# symbolic functions of the symbolic function spaces - to be used as trial and test functions\n",
+ "u0, v0 = elements_of(V0, names='u0, v0')\n",
+ "u1, v1 = elements_of(V1, names='u1, v1')\n",
+ "u2, v2 = elements_of(V2, names='u2, v2')\n",
+ "u3, v3 = elements_of(V3, names='u3, v3')\n",
+ "\n",
+ "# Bilinear Forms corresponding to mass matrices of all four function spaces\n",
+ "m0 = BilinearForm((u0, v0), integral(domain, u0*v0))\n",
+ "m1 = BilinearForm((u1, v1), integral(domain, inner(u1, v1)))\n",
+ "m2 = BilinearForm((u2, v2), integral(domain, inner(u2, v2)))\n",
+ "m3 = BilinearForm((u3, v3), integral(domain, u3*v3))\n",
+ "\n",
+ "# Discretization of the mass-matrix bilinear forms\n",
+ "m0h = discretize(m0, domain_h, (V0h, V0h), backend=backend)\n",
+ "m1h = discretize(m1, domain_h, (V1h, V1h), backend=backend)\n",
+ "m2h = discretize(m2, domain_h, (V2h, V2h), backend=backend)\n",
+ "m3h = discretize(m3, domain_h, (V3h, V3h), backend=backend)\n",
+ "\n",
+ "# Assembly of the mass matrices\n",
+ "M0 = m0h.assemble()\n",
+ "M1 = m1h.assemble()\n",
+ "M2 = m2h.assemble()\n",
+ "M3 = m3h.assemble()\n",
+ "\n",
+ "# Dirichlet projectors in order to apply the projection method\n",
+ "DP0, DP1, _, _ = derham_h.dirichlet_projectors(kind='linop')\n",
+ "I0 = IdentityOperator(V0cs)\n",
+ "\n",
+ "# Modified mass matrix for the proper computation of the weak divergence\n",
+ "M0_0 = DP0 @ M0 @ DP0 + (I0 - DP0)\n",
+ "\n",
+ "# Conjugate Gradient inverse of M0_0, preconditioned using an LST (Loli, Sangalli, Tani) preconditioner\n",
+ "M0_0_pc, = derham_h.LST_preconditioners(M0=M0, hom_bc=True)\n",
+ "M0_0_inv = inverse(M0_0, 'CG', pc=M0_0_pc, maxiter=1000, tol=1e-15)\n",
+ "\n",
+ "# System matrix\n",
+ "S = C.T @ M2 @ C + M1 @ G @ M0_0_inv @ DP0 @ G.T @ M1"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1dec8028",
+ "metadata": {},
+ "source": [
+ "### Assembly of the right-hand-side\n",
+ "\n",
+ "corresponding to the linear form\n",
+ "\n",
+ "$V^1_h \\ni v \\mapsto \\int v\\cdot\\widetilde{curl}_h B$.\n",
+ "\n",
+ "We also check whether\n",
+ "\n",
+ "$div(B) = 0\\text{ in }\\Omega\\quad\\quad\\text{ and }\\quad\\quad n\\cdot B=0\\text{ on }\\partial\\Omega$\n",
+ "\n",
+ "hold."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7cf36d2f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.topology import Union, NormalVector\n",
+ "\n",
+ "# Define callable function corresponding to B\n",
+ "def get_B(r, R, C0=-10):\n",
+ " rad = lambda x, y : np.sqrt(x**2 + y**2)\n",
+ "\n",
+ " B1 = lambda x, y, z : (1/rad(x, y)**2) * ( (-C0 * y) + x * (2*z-1) * (rad(x, y)-r) * (rad(x, y)-R) )\n",
+ " B2 = lambda x, y, z : (1/rad(x, y)**2) * ( ( C0 * x) + y * (2*z-1) * (rad(x, y)-r) * (rad(x, y)-R) )\n",
+ " B3 = lambda x, y, z : (1/rad(x, y)**2) * ( rad(x, y) * ((r+R)-2*rad(x, y)) * z * (z-1) )\n",
+ "\n",
+ " B = (B1, B2, B3)\n",
+ " return B\n",
+ "\n",
+ "# Project B into V2h\n",
+ "_, _, P2, _ = derham_h.projectors()\n",
+ "r, R = domain.logical_domain.bounds1\n",
+ "\n",
+ "B = P2(get_B(r, R))\n",
+ "b = B.coeffs\n",
+ "\n",
+ "# Check whether B is indead divergence free\n",
+ "div_b = D @ b\n",
+ "div_b_norm = np.sqrt( M3.dot_inner(div_b, div_b) )\n",
+ "print(f'|| div(B) ||_L2(domain) = {div_b_norm:.3g}')\n",
+ "\n",
+ "# Check whether B is indead tangential to the boundary\n",
+ "def get_boundaries(*args):\n",
+ "\n",
+ " if not args:\n",
+ " return ()\n",
+ " else:\n",
+ " assert all(1 <= a <= 6 for a in args)\n",
+ " assert len(set(args)) == len(args)\n",
+ "\n",
+ " boundaries = {1: {'axis': 0, 'ext': -1},\n",
+ " 2: {'axis': 0, 'ext': 1},\n",
+ " 3: {'axis': 1, 'ext': -1},\n",
+ " 4: {'axis': 1, 'ext': 1},\n",
+ " 5: {'axis': 2, 'ext': -1},\n",
+ " 6: {'axis': 2, 'ext': 1}}\n",
+ "\n",
+ " return tuple(boundaries[i] for i in args)\n",
+ "\n",
+ "dir_zero_boundary = get_boundaries(1, 2, 5, 6) # we exclude the \"periodic boundary\" in y-direction as it is not technically a boundary anymore\n",
+ "boundary = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary])\n",
+ "\n",
+ "# Assembly of \"tangential trace mass matrix\"\n",
+ "nn = NormalVector('nn')\n",
+ "m2_bd = BilinearForm((u2, v2), integral(boundary, inner(nn, u2) * inner(nn, v2)))\n",
+ "m2_bd_h = discretize(m2_bd, domain_h, (V2h, V2h), backend=backend, sum_factorization=False)\n",
+ "M2_bd = m2_bd_h.assemble()\n",
+ "\n",
+ "dirichlet_boundary_norm = np.sqrt( M2_bd.dot_inner(b, b) )\n",
+ "print(f'|| n \\\\cdot B ||_L2(boundary) = {dirichlet_boundary_norm:.3g}')\n",
+ "\n",
+ "# Assembly of the right-hand-side\n",
+ "rhs = C.T @ M2 @ b"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "890ee988",
+ "metadata": {},
+ "source": [
+ "### Employ the projection method and solve for $A_0$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7b3b6552",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "I1 = IdentityOperator(V1cs)\n",
+ "\n",
+ "# Modify system matrix and right-hand-side in order to take care of boundary conditions\n",
+ "S_0 = DP1 @ S @ DP1 + (I1 - DP1)\n",
+ "rhs_0 = DP1 @ rhs\n",
+ "\n",
+ "# Solve for A_0\n",
+ "maxiter = 1000\n",
+ "tol = 1e-10\n",
+ "\n",
+ "S_0_inv = inverse(S_0, 'cg', maxiter=maxiter, tol=tol)\n",
+ "\n",
+ "a_0 = S_0_inv @ rhs_0"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "48eb3ff5",
+ "metadata": {},
+ "source": [
+ "## Step 2 : Computation of the vector potential $A_H$ of $B_H$\n",
+ "S.th. $curl(A_H) = B_H$ and $\\|curl(A_H)\\| = 1$.\n",
+ "\n",
+ "The theoretical background for the computation of vector potentials of tangential harmonic fields has been recently published\n",
+ "by Martin Campos Pinto and Julian Owezarek (https://arxiv.org/abs/2508.16822).\n",
+ "\n",
+ "The potential is decomposed into $A_H = A_H^b + A_H^0$. We call $A_H^b\\in H(curl)$ the lifting potential and $A_H^0\\in H_0(curl)$ the correction potential.\n",
+ "\n",
+ "### Construct a suitable choice of $A^b_H$ for the torus"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "83d2264e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.fem.basic import FemField\n",
+ "from psydac.linalg.utilities import array_to_psydac\n",
+ "\n",
+ "def get_lifting_field():\n",
+ " dim = V1cs.dimension\n",
+ " array = np.zeros(dim)\n",
+ "\n",
+ " nx, ny, nz = ncells\n",
+ " px, py, pz = degree\n",
+ "\n",
+ " i3 = int(np.floor((nz+pz-1)/2))\n",
+ " Nz = nz + pz\n",
+ "\n",
+ " start = (nx+px-1)*ny*(nz+pz) + (nx+px)*ny*(nz+pz)\n",
+ "\n",
+ " for i2 in range(ny):\n",
+ " index = i2*(Nz-1) + i3\n",
+ " array[start + index] = 1\n",
+ " \n",
+ " coeffs = array_to_psydac(array, V1cs)\n",
+ " return coeffs\n",
+ "\n",
+ "a_H_b = get_lifting_field()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "45008c1a",
+ "metadata": {},
+ "source": [
+ "### Visualization of the lifting field\n",
+ "\n",
+ "\n",
+ " ![](\"./paraview_images/vector_potential_A_H_b.png\")
\n",
+ " ![](\"./paraview_images/vector_potential_A_H_b_lines.png\")
\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "60cf5046",
+ "metadata": {},
+ "source": [
+ "### Compute the correction potential $A^0_H$\n",
+ "\n",
+ "It's yet again a solution to the previously introduced Hodge-Laplace problem, but now with a different right-hand-side\n",
+ "\n",
+ "$\\left( \\widetilde{curl}_h\\ curl\\ - grad\\ \\widetilde{div}_h \\right) A^0_H = - \\widetilde{curl}_h\\ curl A^b_H.$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7a862b0d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "rhs_aH0 = - C.T @ M2 @ C @ a_H_b\n",
+ "rhs_aH0_0 = DP1 @ rhs_aH0\n",
+ "\n",
+ "a_H_0 = S_0_inv @ rhs_aH0_0"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "d5fc3b98",
+ "metadata": {},
+ "source": [
+ "### Assemble the (normalized) harmonic potential $A_H$ and compute $\\lambda$\n",
+ "\n",
+ "s.th. $B_H = curl(\\lambda\\ A_H)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "eefee648",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a_H = a_H_0 + a_H_b\n",
+ "a_H /= np.sqrt( M2.dot_inner(C @ a_H, C @ a_H) )\n",
+ "\n",
+ "b_H = C @ a_H\n",
+ "lam = M2.dot_inner(b, b_H)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3030f8ef",
+ "metadata": {},
+ "source": [
+ "## Step 3 : Assemble the entire vector potential and verify the solution\n",
+ "\n",
+ "$A = A_0 + \\lambda A_H$\n",
+ "\n",
+ "$\\epsilon := \\| curl(A) - B \\|_{L^2}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d27f99d5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "a = a_0 + lam * a_H\n",
+ "\n",
+ "diff = C @ a - b\n",
+ "err_l2 = np.sqrt( M2.dot_inner(diff, diff) )\n",
+ "print(f'|| curl(A) - B ||_L2 = {err_l2:.3g}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "701a6557",
+ "metadata": {},
+ "source": [
+ "## Step 4 : Save results"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8e573154",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.api.postprocessing import OutputManager\n",
+ "from psydac.api.postprocessing import PostProcessManager\n",
+ "import os\n",
+ "\n",
+ "# B = curl(A_0) + curl(lam * A_H) =: B_0 + B_H = curl(A_0 + lam * A_H) =: curl(A)\n",
+ "B_H = FemField(V2h, C @ (lam * a_H))\n",
+ "B_0 = FemField(V2h, b - C @ (lam * a_H))\n",
+ "curlA = FemField(V2h, C @ a)\n",
+ "\n",
+ "# A = A_0 + lam * A_H, here as above we write A_H and mean lam * A_H\n",
+ "A = FemField(V1h, a)\n",
+ "A_H = FemField(V1h, lam * a_H)\n",
+ "A_0 = FemField(V1h, a_0)\n",
+ "\n",
+ "# lifting field A_H_b\n",
+ "A_H_b = FemField(V1h, a_H_b)\n",
+ "\n",
+ "# Save the results using OutputManager\n",
+ "# Export the results to VTK using PostProcessManager\n",
+ "\n",
+ "os.makedirs('results_vector_potential', exist_ok=True)\n",
+ "\n",
+ "Om = OutputManager(\n",
+ " 'results_vector_potential/space_info.yml',\n",
+ " 'results_vector_potential/field_info.h5',\n",
+ " comm=comm,\n",
+ " save_mpi_rank=True, \n",
+ " mode = 'w' \n",
+ ")\n",
+ "\n",
+ "Om.add_spaces(V1=V1h)\n",
+ "Om.add_spaces(V2=V2h)\n",
+ "Om.export_space_info()\n",
+ "\n",
+ "Om.set_static() # the fields do not depend on time\n",
+ "\n",
+ "save_fields = {\n",
+ " 'B' : B,\n",
+ " 'B_H' : B_H,\n",
+ " 'B_0' : B_0,\n",
+ " 'curlA' : curlA,\n",
+ " 'A' : A,\n",
+ " 'A_H' : A_H,\n",
+ " 'A_0' : A_0,\n",
+ " 'A_H_b' : A_H_b\n",
+ "}\n",
+ "\n",
+ "Om.export_fields(**save_fields) # saves the coefficients of the fields\n",
+ "\n",
+ "Om.close()\n",
+ "\n",
+ "# This generates a YAML file containing the space information, and a HDF5 file that stores coefficients. \n",
+ "# At this point the solution has not been evaluated on any grid.\n",
+ "\n",
+ "Pm = PostProcessManager(\n",
+ " domain=domain,\n",
+ " space_file=f'results_vector_potential/space_info.yml',\n",
+ " fields_file=f'results_vector_potential/field_info.h5',\n",
+ " comm=comm\n",
+ ")\n",
+ "\n",
+ "Pm.export_to_vtk(\n",
+ " f'results_vector_potential/visu',\n",
+ " grid=None,\n",
+ " npts_per_cell=[6,6,6],\n",
+ " fields=save_fields.keys() # evaluate all the fields that were saved\n",
+ ")\n",
+ "\n",
+ "Pm.close()\n",
+ "# If this file is run in sequential, then this step will generate a single .vtu file that can be opened with ParaView.\n",
+ "# If it is run in parallel, it generates a .vtu file per process (containing the local grid) and a .pvtu file. To visualize it in ParaView, open the .pvtu file."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "83437c67",
+ "metadata": {},
+ "source": [
+ "## Visualization in ParaView\n",
+ "\n",
+ "### Compare $B$ and $curl(A)$\n",
+ "\n",
+ "\n",
+ " ![](\"./paraview_images/vector_potential_B.png\")
\n",
+ " ![](\"./paraview_images/vector_potential_curlA.png\")
\n",
+ "
\n",
+ "\n",
+ "### Inverstigate decomposition of $B$ into $B_0$, the harmonic-free part, and $B_H$, the harmonic part\n",
+ "\n",
+ "\n",
+ " ![](\"./paraview_images/vector_potential_B_0.png\")
\n",
+ " ![](\"./paraview_images/vector_potential_B_H.png\")
\n",
+ "
\n",
+ "\n",
+ "### Invetsigate decomposition of $A$ into $A_0$, the potential of $B_0$, and $A_H$, the potential of $B_H$\n",
+ "\n",
+ "\n",
+ " ![](\"./paraview_images/vector_potential_A.png\")
\n",
+ " ![](\"./paraview_images/vector_potential_A_0.png\")
\n",
+ " ![](\"./paraview_images/vector_potential_A_H.png\")
\n",
+ "
"
+ ]
+ }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/notebooks/fem_L2_projection.ipynb b/examples/notebooks/fem_L2_projection.ipynb
new file mode 100644
index 000000000..96b0d585b
--- /dev/null
+++ b/examples/notebooks/fem_L2_projection.ipynb
@@ -0,0 +1,498 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "2443e734",
+ "metadata": {},
+ "source": [
+ "# Finite element projections of scalar and vector-valued functions\n",
+ "\n",
+ "*Warning: this notebook is a draft, meant to discuss which interfaces should be exposed to users*\n",
+ "\n",
+ "In this notebook we show how to: \n",
+ "- define FEM spaces over a given domain, \n",
+ "- plot some FEM functions with a simple 2D plotting function,\n",
+ "- apply H1 and Hcurl conforming projections to discrete multipatch fields,\n",
+ "- assemble mass matrices and compute the L2 projections of some given functions."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "455febb2",
+ "metadata": {},
+ "source": [
+ "## Step 1 : define a domain with sympde"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3f2bc988",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from sympde.topology import Square, PolarMapping\n",
+ "\n",
+ "# Define each patch as a mapped subdomain\n",
+ "rmin, rmax = 0.3, 1.\n",
+ "\n",
+ "domain_log_1 = Square('A_1', bounds1=(0., 1.), bounds2=(0., np.pi/2))\n",
+ "F_1 = PolarMapping('F_1', dim=2, c1=0., c2=0., rmin=rmin, rmax=rmax)\n",
+ "Omega_1 = F_1(domain_log_1)\n",
+ "\n",
+ "domain_log_2 = Square('A_2', bounds1=(0., 1.), bounds2=(np.pi, 2 * np.pi))\n",
+ "F_2 = PolarMapping('F_2', dim=2, c1=rmin+rmax, c2=0., rmin=rmax, rmax=rmin)\n",
+ "Omega_2 = F_2(domain_log_2)\n",
+ "\n",
+ "domain_log_3 = Square('A_3', bounds1=(0., 1.), bounds2=(np.pi/2, np.pi))\n",
+ "F_3 = PolarMapping('F_3', dim=2, c1=2*(rmin+rmax), c2=0., rmin=rmin, rmax=rmax)\n",
+ "Omega_3 = F_3(domain_log_3)\n",
+ "\n",
+ "# Join the patches with the Domain.join function from sympde\n",
+ "from sympde import Domain\n",
+ "connectivity = [((Omega_1,1,-1), (Omega_2,1,-1), 1), ((Omega_2,1,1), (Omega_3,1,1), 1)]\n",
+ "patches = [Omega_1, Omega_2, Omega_3]\n",
+ "Omega = Domain.join(patches, connectivity, 'Omega')\n",
+ "\n",
+ "# Simple visualization of the domain\n",
+ "from sympde.utilities.utils import plot_domain\n",
+ "npatches = len(Omega.subdomains)\n",
+ "print(f'A domain named {Omega}, with {npatches} patches:')\n",
+ "plot_domain(Omega, isolines=True)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b1f65c13",
+ "metadata": {},
+ "source": [
+ "## Step 2 : define scalar and vector FEM spaces, plot some basis functions\n",
+ "\n",
+ "We first define a scalar-valued FEM space V"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "39c21151",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace\n",
+ "\n",
+ "# the kind argument specifies the pushforward defining the FEM space on the physical domain\n",
+ "# kind=None results in a space defined by a simple change of variable\n",
+ "V = ScalarFunctionSpace('V', Omega, kind=None)\n",
+ "\n",
+ "ncells = [10, 10]\n",
+ "degree = [2, 2]\n",
+ "\n",
+ "from psydac.api.discretization import discretize\n",
+ "Omega_h = discretize(Omega, ncells=ncells)\n",
+ "Vh = discretize(V, Omega_h, degree=degree)\n",
+ "print(f\"we have defined a FEM space of type {type(Vh)}, with dimension {Vh.nbasis}\")\n",
+ "\n",
+ "# plot some basis function from the FEM space, using a simple wrapper for matplotlib\n",
+ "from psydac.fem.basic import FemField\n",
+ "from psydac.linalg.utilities import array_to_psydac\n",
+ "\n",
+ "# here we set the coefficients by hand using numpy arrays\n",
+ "bf_c = np.zeros(Vh.nbasis)\n",
+ "i0 = 108\n",
+ "bf_c[i0] = 1\n",
+ "i1 = 211\n",
+ "bf_c[i1] = 1\n",
+ "i2 = 215\n",
+ "bf_c[i2] = 1\n",
+ "bf = FemField(Vh, coeffs=array_to_psydac(bf_c, Vh.coeff_space))\n",
+ "\n",
+ "from psydac.fem.plotting_utilities import plot_field_2d\n",
+ "plot_field_2d(fem_field=bf, Vh=Vh, domain=Omega, cmap='viridis', title=fr'three scalar-valued basis functions: $\\phi_{ {i0} }$, $\\phi_{ {i1} }$, $\\phi_{ {i2} }$', filename='')\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "79eb29ab",
+ "metadata": {},
+ "source": [
+ "Note: the multipatch space Vh is made of discontinuous functions. A conforming projection may be used to map them to continuous functions"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "76d8fa2a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.feec.conforming_projectors import construct_h1_conforming_projection\n",
+ "cP = construct_h1_conforming_projection(Vh, hom_bc=False)\n",
+ "Pbf = FemField(Vh, coeffs=array_to_psydac(cP@bf_c, Vh.coeff_space))\n",
+ "plot_field_2d(fem_field=Pbf, Vh=Vh, domain=Omega, cmap='viridis', title='H1-conforming projection of the above basis functions', filename='')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "1de55bc2",
+ "metadata": {},
+ "source": [
+ "Then we define a vector-valued FEM space W"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d4b778e1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# again, kind=None will result in a FEM space defined by a simple change of variable on the physical domain\n",
+ "W = VectorFunctionSpace('W', Omega, kind=None) \n",
+ "Wh = discretize(W, Omega_h, degree=degree)\n",
+ "print(f\"we have defined a FEM space of type {type(Wh)}, with dimension {Wh.nbasis}\")\n",
+ "\n",
+ "# plot some basis function from the space \n",
+ "# [remark] here I'm using numpy arrays, but one could illustrate the psydac indexing\n",
+ "dim_Wh = Wh.nbasis\n",
+ "bf2_c = np.zeros(dim_Wh)\n",
+ "i0 = 108; bf2_c[i0] = 1\n",
+ "i1 = 314; bf2_c[i1] = 1\n",
+ "i2 = 479; bf2_c[i2] = 1\n",
+ "bf2 = FemField(Wh, coeffs=array_to_psydac(bf2_c, Wh.coeff_space))\n",
+ "\n",
+ "plot_field_2d(fem_field=bf2, Vh=Wh, domain=Omega, \n",
+ " plot_type='components', cmap='viridis', title=fr'three vector-valued basis functions: $\\Phi_{ {i0} }$, $\\Phi_{ {i1} }$, $\\Phi_{ {i2} }$', filename='')\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e2793980",
+ "metadata": {},
+ "source": [
+ "A stream plot of these basis functions is also possible, to better vizualise the vector fields"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5249428d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plot_field_2d(fem_field=bf2, Vh=Wh, domain=Omega, \n",
+ " plot_type='vector_field', vf_skip=1, title=fr'three vector-valued basis functions: $\\Phi_{ {i0} }$, $\\Phi_{ {i1} }$, $\\Phi_{ {i2} }$', filename='')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "21ca10a3",
+ "metadata": {},
+ "source": [
+ "Note: again, the multipatch space Wh is made of discontinuous functions. A conforming projection could be used, but continuity is harder to impose for vector-valued functions: this will be discussed in a further example."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "36a9632b",
+ "metadata": {},
+ "source": [
+ "## Step 3 : compute the mass matrices for Vh and Wh\n",
+ "\n",
+ "We assemble the mass matrices and visualize them. Again we start with the scalar-valued FEM space"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "091d4560",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "from matplotlib import cm, colors\n",
+ "\n",
+ "from sympde.topology import elements_of\n",
+ "from sympde.expr.expr import BilinearForm\n",
+ "from sympde.expr.expr import integral \n",
+ "\n",
+ "from psydac.api.settings import PSYDAC_BACKENDS\n",
+ "\n",
+ "backend_language = 'python'\n",
+ "\n",
+ "u, v = elements_of(V, names='u, v')\n",
+ "a = BilinearForm((u, v), integral(Omega, u * v))\n",
+ "ah = discretize(a, Omega_h, [Vh, Vh], backend=PSYDAC_BACKENDS[backend_language])\n",
+ "\n",
+ "# Mass matrix in PSYDAC 'stencil' (linear operator) format \n",
+ "M_V = ah.assemble() \n",
+ "\n",
+ "# visualize the mass matrix by plotting its numpy conversion\n",
+ "mat = M_V.toarray()\n",
+ "fig,ax = plt.subplots(1,1)\n",
+ "ax.set_title(f\"Mass matrix of Vh, with {npatches} blocks (1 per patch)\")\n",
+ "im = ax.matshow(mat, norm=colors.LogNorm(), cmap='hot_r')\n",
+ "cb = fig.colorbar(im, ax=ax)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3a3db420",
+ "metadata": {},
+ "source": [
+ "We next assemble the mass matrix of the vector-valued FEM space"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "aa017e79",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.calculus import dot\n",
+ "\n",
+ "uu, vv = elements_of(W, names='uu, vv')\n",
+ "a = BilinearForm((uu, vv), integral(Omega, dot(uu,vv)))\n",
+ "ah_W = discretize(a, Omega_h, [Wh, Wh], backend=PSYDAC_BACKENDS[backend_language])\n",
+ "\n",
+ "# Mass matrix in PSYDAC 'stencil' (linear operator) format \n",
+ "M_W = ah_W.assemble() \n",
+ "\n",
+ "# visualize the mass matrix by plotting its numpy conversion\n",
+ "mat_W = M_W.toarray()\n",
+ "fig,ax = plt.subplots(1,1)\n",
+ "ax.set_title(f\"Mass matrix of Wh, with {2*npatches} blocks (2 per patch)\")\n",
+ "im = ax.matshow(mat_W, norm=colors.LogNorm(), cmap='hot_r')\n",
+ "cb = fig.colorbar(im, ax=ax)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "dae32f1c",
+ "metadata": {},
+ "source": [
+ "## Step 4 : compute the moments of some given target function, and their L2 projections\n",
+ "\n",
+ "we first define a scalar function f and visualize it on the domain"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "de936090",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympy import cos, exp\n",
+ "x,y = Omega.coordinates\n",
+ "ref_f = exp(-(y-0.7*cos(3*x))**2 / (2 * 0.15**2))\n",
+ "\n",
+ "from sympy import lambdify\n",
+ "f_call = lambdify(Omega.coordinates, ref_f) # make the expression a callable function\n",
+ "\n",
+ "from psydac.fem.plotting_utilities import get_plotting_grid, my_small_plot\n",
+ "etas, xx, yy = get_plotting_grid(Omega.mappings, N=100)\n",
+ "ref_f_vals = [f_call(xx_k, yy_k) for xx_k, yy_k in zip(xx, yy)]\n",
+ "\n",
+ "my_small_plot(\n",
+ " title=r'a reference function $f$ (scalar valued)',\n",
+ " vals=[ref_f_vals],\n",
+ " xx=xx,\n",
+ " yy=yy,\n",
+ " save_fig='',\n",
+ " cmap='viridis',\n",
+ ")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "75014e6a",
+ "metadata": {},
+ "source": [
+ "we next compute the moments of f in Vh and invert the mass matrix to get the L2 projection"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "db7f74ca",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.linalg.solvers import inverse\n",
+ "from psydac.fem.projectors import get_dual_dofs\n",
+ "\n",
+ "# the dual dofs are the moments, i.e. the products _{L2(Omega)}\n",
+ "tilde_f = get_dual_dofs(Vh=Vh, f=ref_f, domain_h=Omega_h, backend_language=backend_language)\n",
+ "\n",
+ "inv_M = inverse(M_V, solver='cg')\n",
+ "f_coeffs = inv_M @ tilde_f\n",
+ "fh = FemField(Vh, f_coeffs)\n",
+ "\n",
+ "plot_field_2d(fem_field=fh, domain=Omega, title=r'$f_h$: L2 projection of $f$ in Vh', cmap='viridis', filename='')\n",
+ "\n",
+ "print('note the lower accuracy in the regions where the grid is not aligned with the anisotropic smoothness')\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0de90ef8",
+ "metadata": {},
+ "source": [
+ "compute the norm and the error"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a23c225a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.expr.expr import Norm \n",
+ "error = ref_f - u\n",
+ "err_l2norm = Norm(error, Omega, kind='l2')\n",
+ "nquads = [p + 1 for p in degree]\n",
+ "err_l2norm_h = discretize(err_l2norm, Omega_h, Vh, nquads=nquads, backend=PSYDAC_BACKENDS[backend_language])\n",
+ "err_l2 = err_l2norm_h.assemble(u=fh)\n",
+ "print(f'> L2 projection error : {err_l2}')\n",
+ "\n",
+ "zero_h = FemField(Vh, coeffs=array_to_psydac(np.zeros(Vh.nbasis), Vh.coeff_space))\n",
+ "f_l2 = err_l2norm_h.assemble(u=zero_h)\n",
+ "print(f'> L2 norm of f : {f_l2}')\n",
+ "\n",
+ "rel_err_f = err_l2 / f_l2\n",
+ "print(f'> L2 relative error : {rel_err_f}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8258971a",
+ "metadata": {},
+ "source": [
+ "we then repeat the above steps for the vector-valued space Wh"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "772df7b3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympy import pi, cos, Tuple\n",
+ "\n",
+ "ref_g_x = ref_f * cos(2*pi*y)\n",
+ "ref_g_y = ref_f * 1\n",
+ "ref_g = Tuple(ref_g_x, ref_g_y)\n",
+ "\n",
+ "# make the expression a callable function\n",
+ "g_vect = [lambdify(Omega.coordinates, ref_g[d]) for d in [0,1]]\n",
+ "\n",
+ "ref_vals_g = [[g_vect[d](xx_k, yy_k) for xx_k, yy_k in zip(xx, yy)] for d in [0,1]]\n",
+ "\n",
+ "from psydac.fem.plotting_utilities import my_small_streamplot\n",
+ "my_small_streamplot(\n",
+ " title=r'a reference function $g$ (vector valued)',\n",
+ " vals_x=ref_vals_g[0],\n",
+ " vals_y=ref_vals_g[1],\n",
+ " xx=xx,\n",
+ " yy=yy,\n",
+ " skip=8,\n",
+ " amp_factor=1,\n",
+ " hide_plot=False,\n",
+ " dpi=200,\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1f0f103d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# again, the dual dofs are the moments relative to the basis of Wh, i.e. the products _{L2(Omega)}\n",
+ "tilde_g = get_dual_dofs(Vh=Wh, f=ref_g, domain_h=Omega_h, backend_language=backend_language)\n",
+ "\n",
+ "inv_M_W = inverse(M_W, solver='cg')\n",
+ "g_coeffs = inv_M_W @ tilde_g\n",
+ "\n",
+ "gh = FemField(Wh, g_coeffs)\n",
+ "\n",
+ "plot_field_2d(fem_field=gh, domain=Omega, plot_type='vector_field', title='g_h: L2 projection of g', filename='')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3f9d6e2b",
+ "metadata": {},
+ "source": [
+ "compute the norm and the error"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d6e6b7b4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympy import Matrix \n",
+ "\n",
+ "error = Matrix([ref_g[0]-uu[0], ref_g[1]-uu[1]])\n",
+ "err_l2norm = Norm(error, Omega, kind='l2')\n",
+ "nquads = [p + 1 for p in degree]\n",
+ "err_l2norm_h = discretize(err_l2norm, Omega_h, Wh, nquads=nquads, backend=PSYDAC_BACKENDS[backend_language])\n",
+ "err_l2 = err_l2norm_h.assemble(uu=gh)\n",
+ "print(f'> L2 projection error : {err_l2}')\n",
+ "\n",
+ "zero_h = FemField(Wh, coeffs=array_to_psydac(np.zeros(Wh.nbasis), Wh.coeff_space))\n",
+ "g_l2 = err_l2norm_h.assemble(uu=zero_h)\n",
+ "print(f'> L2 norm of g : {g_l2}')\n",
+ "\n",
+ "rel_err_g = err_l2 / g_l2\n",
+ "print(f'> L2 relative error : {rel_err_g}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a5b15eef",
+ "metadata": {},
+ "source": [
+ "### Testing the notebook"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2a7634d4",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import ipytest\n",
+ "ipytest.autoconfig(raise_on_error=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "caa039a8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%ipytest\n",
+ "\n",
+ "def test_l2error_f():\n",
+ " assert rel_err_f < 5e-02\n",
+ "\n",
+ "def test_l2error_g():\n",
+ " assert rel_err_g < 5e-02"
+ ]
+ }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/notebooks/paraview_images/regularized_curlcurl_3D_fx.png b/examples/notebooks/paraview_images/regularized_curlcurl_3D_fx.png
new file mode 100644
index 000000000..00d948219
Binary files /dev/null and b/examples/notebooks/paraview_images/regularized_curlcurl_3D_fx.png differ
diff --git a/examples/notebooks/paraview_images/regularized_curlcurl_3D_fy.png b/examples/notebooks/paraview_images/regularized_curlcurl_3D_fy.png
new file mode 100644
index 000000000..c77febe72
Binary files /dev/null and b/examples/notebooks/paraview_images/regularized_curlcurl_3D_fy.png differ
diff --git a/examples/notebooks/paraview_images/regularized_curlcurl_3D_fz.png b/examples/notebooks/paraview_images/regularized_curlcurl_3D_fz.png
new file mode 100644
index 000000000..0967959c4
Binary files /dev/null and b/examples/notebooks/paraview_images/regularized_curlcurl_3D_fz.png differ
diff --git a/examples/notebooks/paraview_images/regularized_curlcurl_3D_ux.png b/examples/notebooks/paraview_images/regularized_curlcurl_3D_ux.png
new file mode 100644
index 000000000..0b3410b33
Binary files /dev/null and b/examples/notebooks/paraview_images/regularized_curlcurl_3D_ux.png differ
diff --git a/examples/notebooks/paraview_images/regularized_curlcurl_3D_uy.png b/examples/notebooks/paraview_images/regularized_curlcurl_3D_uy.png
new file mode 100644
index 000000000..169100636
Binary files /dev/null and b/examples/notebooks/paraview_images/regularized_curlcurl_3D_uy.png differ
diff --git a/examples/notebooks/paraview_images/regularized_curlcurl_3D_uz.png b/examples/notebooks/paraview_images/regularized_curlcurl_3D_uz.png
new file mode 100644
index 000000000..e7b7e01e0
Binary files /dev/null and b/examples/notebooks/paraview_images/regularized_curlcurl_3D_uz.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_A.png b/examples/notebooks/paraview_images/vector_potential_A.png
new file mode 100644
index 000000000..f1079991b
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_A.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_A_0.png b/examples/notebooks/paraview_images/vector_potential_A_0.png
new file mode 100644
index 000000000..729b8d3a5
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_A_0.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_A_H.png b/examples/notebooks/paraview_images/vector_potential_A_H.png
new file mode 100644
index 000000000..c3157b11b
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_A_H.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_A_H_b.png b/examples/notebooks/paraview_images/vector_potential_A_H_b.png
new file mode 100644
index 000000000..0002a300d
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_A_H_b.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_A_H_b_lines.png b/examples/notebooks/paraview_images/vector_potential_A_H_b_lines.png
new file mode 100644
index 000000000..1b191c032
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_A_H_b_lines.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_B.png b/examples/notebooks/paraview_images/vector_potential_B.png
new file mode 100644
index 000000000..a3fbdc580
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_B.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_B_0.png b/examples/notebooks/paraview_images/vector_potential_B_0.png
new file mode 100644
index 000000000..0d0c5eec0
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_B_0.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_B_H.png b/examples/notebooks/paraview_images/vector_potential_B_H.png
new file mode 100644
index 000000000..e2e13f442
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_B_H.png differ
diff --git a/examples/notebooks/paraview_images/vector_potential_curlA.png b/examples/notebooks/paraview_images/vector_potential_curlA.png
new file mode 100644
index 000000000..def527bc3
Binary files /dev/null and b/examples/notebooks/paraview_images/vector_potential_curlA.png differ
diff --git a/examples/notebooks/petsc_eigenvalues_regularized_curlcurl.ipynb b/examples/notebooks/petsc_eigenvalues_regularized_curlcurl.ipynb
new file mode 100644
index 000000000..119a01872
--- /dev/null
+++ b/examples/notebooks/petsc_eigenvalues_regularized_curlcurl.ipynb
@@ -0,0 +1,282 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "2443e734",
+ "metadata": {},
+ "source": [
+ "# Example of computing approximative eigenvalues with PETSc\n",
+ "\n",
+ "In this notebook we show how to compute approximative eigenvalues using PETSc. We use Conjugate Gradient to solve a 3D problem in parallel. Besides, we use a Jacobi preconditioner to illustrate how to set up a preconditioned Krylov solver with PETSc. The eigenvalues are estimated using the Krylov method and the number of estimated eigenvalues is the number of iterations.\n",
+ "\n",
+ "Similar to the 3D VTK ParaView test, we solve the regularized curl curl problem in 3D with a callable function for the source. In practice, the source could be the interpolation of some experimental data. Besides, we impose periodic boundary conditions everywhere in the boundary.\n",
+ "\n",
+ "Let $\\Omega=[0, 10] \\times [-2\\pi, 2\\pi]\\times [-3, 1]$. We want to find $u\\in H(\\mathrm{curl}, \\Omega)$ such that\n",
+ "\n",
+ "$$ \\mathrm{curl} \\mathrm{curl} u + \\mu u = f, $$\n",
+ "\n",
+ "where $\\mu = \\pi$, and $f$ is a callable function in terms of the domain coordinates, the expression of which is a priori not known. \n",
+ "\n",
+ "Only for illustration here we use the callable function\n",
+ "$$f(x,y,z) = (10e^{-0.5(x-4)^2 - 1.5z^2 - 0.1y^2}, 0, \\cos(y)\\sin(\\pi x / 5)),$$\n",
+ "\n",
+ "and project it from $H(\\mathrm{curl}, \\Omega)$ to the corresponding discretized space. This process is analogous with any callable vectorial function.\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Step 1 : define the domain and spaces and discretize them\n",
+ "\n",
+ "This is similar to what is done in other examples"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "39c21151",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from mpi4py import MPI\n",
+ "from sympde.topology import Cube\n",
+ "from sympde.topology import VectorFunctionSpace\n",
+ "from psydac.api.discretization import discretize\n",
+ "\n",
+ "comm = MPI.COMM_WORLD # MPI communicator\n",
+ "\n",
+ "Omega = Cube('Omega', bounds1=(0, 10), bounds2=(-2*np.pi, 2*np.pi), bounds3=(-3, 1)) # 3D cubic domain\n",
+ "\n",
+ "W = VectorFunctionSpace('W', Omega, kind='hcurl') # Hcurl space\n",
+ "\n",
+ "ncells = [25, 30, 11] # number of cells in each direction\n",
+ "degree = [3, 2, 2] # B-spline degree in each direction\n",
+ "\n",
+ "periodic = [True, True, True] # set periodic BC in every direction\n",
+ "\n",
+ "Omega_h = discretize(Omega, ncells=ncells, periodic=periodic, comm=comm)\n",
+ "Wh = discretize(W, Omega_h, degree=degree)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "36a9632b",
+ "metadata": {},
+ "source": [
+ "## Step 2 : compute the matrix\n",
+ "\n",
+ "We assemble the mass and curl curl matrices."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "091d4560",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.topology import elements_of, element_of\n",
+ "from sympde.expr.expr import BilinearForm, LinearForm, integral\n",
+ "from sympde.calculus import dot, curl\n",
+ "from psydac.api.settings import PSYDAC_BACKENDS\n",
+ "from psydac.fem.basic import FemField\n",
+ "\n",
+ "backend_language = 'pyccel-gcc'\n",
+ "backend = PSYDAC_BACKENDS[backend_language]\n",
+ "\n",
+ "u, v = elements_of(W, names='u, v') #trial and test functions\n",
+ "\n",
+ "m = BilinearForm((u, v), integral(Omega, dot(u, v)))\n",
+ "mh = discretize(m, Omega_h, [Wh, Wh], backend=backend)\n",
+ "M = mh.assemble() # Mass matrix\n",
+ "\n",
+ "k = BilinearForm((u, v), integral(Omega, dot(curl(u), curl(v))))\n",
+ "kh = discretize(k, Omega_h, [Wh, Wh], backend=backend)\n",
+ "K = kh.assemble() # Curl curl matrix\n",
+ "\n",
+ "mu = 1e-1\n",
+ "\n",
+ "A = K + mu*M # system matrix, positive-definite"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8927f88c",
+ "metadata": {},
+ "source": [
+ "## Step 3 : compute the right-hand size\n",
+ "\n",
+ "We project a callable function to the vector space and use it to assemble the right-hand side of the system."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a13044fa",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.feec.global_geometric_projectors import GlobalGeometricProjectorHcurl\n",
+ "f_callable = [lambda x, y, z: 10*np.exp(-0.5*(x-4)**2 - 1.5*z**2 - 1e-1*y**2),\n",
+ " lambda x, y, z: 0,\n",
+ " lambda x, y, z: np.cos(y)*np.sin(x*np.pi/5)]\n",
+ "\n",
+ "proj_Wh = GlobalGeometricProjectorHcurl(Wh, nquads=[p+1 for p in degree]) # get Hcurl projector, nquads is the number of points for the Gauss quadrature.\n",
+ "fh = proj_Wh(f_callable) # project callable to field in Wh\n",
+ "\n",
+ "f = element_of(W, name='f') # free variable for the FemField\n",
+ "\n",
+ "l = LinearForm(v, integral(Omega, dot(f,v)))\n",
+ "lh = discretize(l, Omega_h, Wh, backend=backend)\n",
+ "rhs = lh.assemble(f=fh) # assmble RHS, f is a free parameter"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8258971a",
+ "metadata": {},
+ "source": [
+ "## Step 4 : solve system with PETSc\n",
+ "First we convert matrix to a PETSc.Mat object and the right-hand side to a PETSc.Vec object. Then we solve the system with conjugate gradient and Jacobi preconditioner. We also compute the approximative eigenvalues during the Krylov method."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1f0f103d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from petsc4py import PETSc\n",
+ "from psydac.linalg.utilities import petsc_to_psydac\n",
+ "\n",
+ "Ap = A.topetsc()\n",
+ "rhsp = rhs.topetsc()\n",
+ "\n",
+ "# initial solution\n",
+ "u = rhsp.copy()\n",
+ "u.zeroEntries() # set all entries to zero\n",
+ "\n",
+ "iterative_solver = PETSc.KSP().create(comm=comm) # create Krylov solver object\n",
+ "iterative_solver.setType(PETSc.KSP.Type.CG)\n",
+ "iterative_solver.setOperators(A=Ap, P=Ap) # set operator and preconditioner\n",
+ "preconditioner = iterative_solver.getPC()\n",
+ "preconditioner.setType(PETSc.PC.Type.JACOBI) # set Jacobi preconditioner, i.e. diagonal of matrix Ap\n",
+ "preconditioner.setUp()\n",
+ "\n",
+ "iterative_solver.setNormType(PETSc.KSP.NormType.UNPRECONDITIONED) #use unpreconditioned norm for convergence\n",
+ "\n",
+ "iterative_solver.setComputeEigenvalues(True) # compute eigenvalues during Krylov solver\n",
+ "def callback(ksp, iter, rnorm):\n",
+ " if ksp.getComm().Get_rank() == 0:\n",
+ " print(f'Iter {iter} | res = {\"{:.2e}\".format(rnorm)}')\n",
+ "\n",
+ "iterative_solver.setMonitor(callback) # set a callback function to print iteration and residual.\n",
+ "iterative_solver.setTolerances(atol=1e-5, rtol=1e-2, max_it=100)\n",
+ "iterative_solver.setInitialGuessNonzero(True) #if True it uses the solution vector as initial guess, otherwise it takes zero.\n",
+ "\n",
+ "iterative_solver.solve(rhsp, u)\n",
+ "\n",
+ "eigenvalues_approx = iterative_solver.computeEigenvalues() # the array size is the number of CG iterations\n",
+ "\n",
+ "# convert the solution back to PSYDAC\n",
+ "u_psydac = petsc_to_psydac(u, Wh.coeff_space)\n",
+ "uh = FemField(Wh, coeffs=u_psydac)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ab581ba7",
+ "metadata": {},
+ "source": [
+ "## Step 5 : Plot approximative eigenvalues\n",
+ "The eigenvalues computed with the Krylov method are only approximative. The number of computed \"eigenvalues\" is the number of iterations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1a04a463",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "\n",
+ "if comm is None or comm.Get_rank() == 0:\n",
+ " fig,ax = plt.subplots(1,1)\n",
+ " ax.set_title(f\"Eigenvalues estimated with CG {iterative_solver.getIterationNumber()} iterations\")\n",
+ " im = ax.scatter(np.arange(eigenvalues_approx.size), eigenvalues_approx)\n",
+ " plt.xlabel('index')\n",
+ " plt.show()\n",
+ " \n",
+ "if comm is not None:\n",
+ " comm.Barrier()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e0c73790",
+ "metadata": {},
+ "source": [
+ "The estimated eigenvalues are shown above.\n",
+ "We can see that the eigenavalues corresponding to the kernel of the curlcurl operator are shifted thanks to the regularization.\n",
+ "\n",
+ "\n",
+ "## Step 6 : Save and export solution\n",
+ "This step is similar to what is explained in other tests."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6c77d25a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Save the results using OutputManager\n",
+ "from psydac.api.postprocessing import OutputManager\n",
+ "import os\n",
+ "\n",
+ "os.makedirs('results_eigenvalues_regularized_curlcurl', exist_ok=True)\n",
+ "\n",
+ "Om = OutputManager(\n",
+ " f'results_eigenvalues_regularized_curlcurl/space_info_{Omega.name}.yml',\n",
+ " f'results_eigenvalues_regularized_curlcurl/field_info_{Omega.name}.h5',\n",
+ " comm=comm,\n",
+ " save_mpi_rank=True, \n",
+ " mode = 'w' \n",
+ ")\n",
+ "\n",
+ "Om.add_spaces(V=Wh)\n",
+ "Om.export_space_info()\n",
+ "\n",
+ "Om.set_static() # the fields do not depend on time\n",
+ "\n",
+ "Om.export_fields(u=uh) # saves the coefficients of the solution\n",
+ "\n",
+ "Om.close()\n",
+ "\n",
+ "# Export the results to VTK using PostProcessManager\n",
+ "from psydac.api.postprocessing import PostProcessManager\n",
+ "\n",
+ "Pm = PostProcessManager(\n",
+ " domain=Omega,\n",
+ " space_file=f'results_eigenvalues_regularized_curlcurl/space_info_{Omega.name}.yml',\n",
+ " fields_file=f'results_eigenvalues_regularized_curlcurl/field_info_{Omega.name}.h5',\n",
+ " comm=comm\n",
+ ")\n",
+ "\n",
+ "Pm.export_to_vtk(\n",
+ " f'results_eigenvalues_regularized_curlcurl/visu_{Omega.name}',\n",
+ " grid=None,\n",
+ " npts_per_cell=[4,2,4],\n",
+ " fields='u'\n",
+ ")\n",
+ "\n",
+ "Pm.close()"
+ ]
+ }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/notebooks/poisson_2d_square.ipynb b/examples/notebooks/poisson_2d_square.ipynb
new file mode 100644
index 000000000..8a2e6f3e0
--- /dev/null
+++ b/examples/notebooks/poisson_2d_square.ipynb
@@ -0,0 +1,333 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Solving Poisson's equation on a square domain with a mapped grid\n",
+ "\n",
+ "In this notebook, we will show how to solve Poisson's equation on a single-patch square domain.\n",
+ "Thanks to the simple geometry, will use the method of manufactured solution to build a test case where the exact solution is known a priori.\n",
+ "As we are running this in a notebook, everything will be run in serial and hence we are limiting ourselves to a fairly coarse discretization to avoid taking too much time.\n",
+ "However, PSYDAC allows for hybrid MPI + OpenMP parallelization with barely any changes to the code.\n",
+ "The lines that are impacted by that will be preceeded by their MPI equivalent."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 1: Problem definition with SymPDE\n",
+ "\n",
+ "PSYDAC uses the powerful topological tools of SymPDE to build a large variety of domains.\n",
+ "Here we use SymPDE to map a square domain $\\widehat{\\Omega}$ to another square $\\Omega := F(\\widehat{\\Omega})$ with non-orthogonal coordinates.\n",
+ "On such a simple domain, we define an exact solution $u$ which satisfies homogeneous Dirichlet boundary conditions, and we compute the right-hand side $f$ of the equation.\n",
+ "The strong form of Poisson's equation reads simply\n",
+ "\n",
+ "$$\n",
+ "\\left\\{\n",
+ "\\begin{aligned}\n",
+ " -\\nabla^2 u &= f && \\text{in $\\Omega$}, \\\\\n",
+ " u &= 0 && \\text{on $\\partial\\Omega$}.\n",
+ "\\end{aligned}\n",
+ "\\right.\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympy import sin, pi\n",
+ "\n",
+ "from sympde.topology import Square, CollelaMapping2D\n",
+ "from sympde.calculus import laplace\n",
+ "from sympde.utilities.utils import plot_domain\n",
+ "\n",
+ "# Problem definition\n",
+ "OmegaP = Square('Omega') # parametric domain\n",
+ "F = CollelaMapping2D('F', eps=0.1, k1=1, k2=1) # mapping\n",
+ "Omega = F(OmegaP) # physical domain\n",
+ "x, y = Omega.coordinates # physical coordinates\n",
+ "u_ex = sin(pi * x) * sin(pi * y) # manufactured solution\n",
+ "f = -laplace(u_ex) # right-hand side\n",
+ "\n",
+ "# Simple visualization of the domain\n",
+ "plot_domain(Omega, isolines=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 2: Abstract weak formulation with SymPDE\n",
+ "\n",
+ "The weak formulation of Poisson's equation reads:\n",
+ "\n",
+ "$$\n",
+ "\\text{Given $f \\in L^2(\\Omega)$},\n",
+ "\\quad\n",
+ "\\text{find $u \\in H^1_0(\\Omega)$ s.t.}\n",
+ "\\quad\n",
+ "\\underbrace{\\int_{\\Omega} \\nabla u \\cdot \\nabla v\\,\\mathrm{d}\\Omega}_{a(u, v)} =\n",
+ "\\underbrace{\\int_{\\Omega} f\\,v \\,\\mathrm{d}\\Omega}_{l(v)}\n",
+ "\\quad\n",
+ "\\forall v \\in H^1_0(\\Omega).\n",
+ "$$\n",
+ "\n",
+ "Here $H^1_0(\\Omega)$ satisfies homogeneous Dirichlet boundary conditions.\n",
+ "Currently, SymPDE does not allow attaching this information to the function space.\n",
+ "(This is consistent with the fact that the FEM fields in PSYDAC are in general non-zero at the boundary.)\n",
+ "Instead, an additional constraint is attached to the weak formulation explicitly, in the form of an `EssentialBC` object.\n",
+ "\n",
+ "Since we are using the method of manufactured solutions, we can compute the $L^2$ and $H^1$ norms of the error.\n",
+ "The norms are defined symbolically in SymPDE using `Norm` objects.\n",
+ "The error is $u-u_{ex}$, where the trial function $u$ is a SymPDE `ScalarFunction` object, while the exact solution $u_{ex}$ is a SymPy expression of the physical coordinates.\n",
+ "Since $u$ represents a free parameter, it will have to be replaced by a PSYDAC FEM field when assemblying the norm."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.topology import ScalarFunctionSpace\n",
+ "from sympde.topology import elements_of\n",
+ "from sympde.topology import NormalVector\n",
+ "from sympde.calculus import grad, inner\n",
+ "from sympde.expr import LinearForm, BilinearForm, Norm\n",
+ "from sympde.expr import integral\n",
+ "from sympde.expr import find, EssentialBC\n",
+ "\n",
+ "# Space, trial function u, and test function v\n",
+ "V = ScalarFunctionSpace('V', Omega)\n",
+ "u, v = elements_of(V, names='u, v')\n",
+ "\n",
+ "# Bilinear form a(u, v) and linear form l(v)\n",
+ "a = BilinearForm((u, v), integral(Omega, inner(grad(u), grad(v))))\n",
+ "l = LinearForm(v, integral(Omega, f * v))\n",
+ "\n",
+ "# Homogeneous Dirichlet boundary conditions: u=0 on ∂Ω\n",
+ "bc = EssentialBC(u, 0, Omega.boundary)\n",
+ "\n",
+ "# Weak formulation of Poisson's eqn. with boundary conditions as constraint\n",
+ "equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc)\n",
+ "\n",
+ "# L^2 and H^1 norms of the error\n",
+ "l2norm = Norm(u - u_ex, Omega, kind='l2')\n",
+ "h1norm = Norm(u - u_ex, Omega, kind='h1')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 3: Galerkin method with PSYDAC\n",
+ "\n",
+ "Starting from the weak formulation of Poisson's equation, a Galerkin method of solution can be formulated as follows:\n",
+ "\n",
+ "- Let $u_h \\approx u$ in a finite dimensional subspace $V_h \\subset H^1_0$\n",
+ "- Choose a basis for $V_h = span(\\phi_1, \\phi_2, \\dots \\phi_N)$\n",
+ "- Search for $u_h(x,y) := \\sum_{j=1}^N w_j\\, \\phi_j(x,y)$\n",
+ "- Choose test functions $v = \\phi_i$\n",
+ "- A linear PDE yields the linear system $\\mathbb{A}~\\mathsf{w} = \\mathsf{b}$\n",
+ "\n",
+ "$$\n",
+ "\\text{Given $f \\in L^2(\\Omega)$},\n",
+ "\\quad\n",
+ "\\text{find $\\mathsf{w}\\in\\mathbb{R}^N$ s.t.}\n",
+ "\\quad\n",
+ "\\sum_{j=1}^N \\underbrace{\\left(\\int_{\\Omega_h} \\nabla\\phi_i \\cdot \\nabla\\phi_j\\,\\mathrm{d}\\Omega \\right)}_{A_{ij}} w_j\n",
+ "= \\underbrace{\\int_{\\Omega_h} f\\,\\phi_i\\,\\mathrm{d}\\Omega}_{b_i},\n",
+ "\\quad\n",
+ "\\text{for $i = 1, \\dots, N$}.\n",
+ "$$\n",
+ "\n",
+ "PSYDAC provides the function `discretize` which maps SymPDE abstract objects to their discrete PSYDAC counterpart.\n",
+ "For a Poisson problem over a single-patch domain, this means:\n",
+ "\n",
+ "| SymPDE | PSYDAC |\n",
+ "| --------------------- | ---------------------- |\n",
+ "| `Domain` | `Geometry` |\n",
+ "| `ScalarFunctionSpace` | `TensorFemSpace` |\n",
+ "| `Equation` | `DiscreteEquation` |\n",
+ "| `Functional` (including `Norm` and `SemiNorm`) | `DiscreteFunctional` |\n",
+ "| `LinearForm` | `DiscreteLinearForm` |\n",
+ "| `BilinearForm` | `DiscreteBilinearForm` |\n",
+ "\n",
+ "In addition, SymPDE's `ScalarFunction` objects correspond to PSYDAC's `FemField` objects.\n",
+ "\n",
+ "When discretizing an object of type `Functional`, `LinearForm`, or `BilinearForm`, PSYDAC automatically generates Python code for the corresponding assembly functions.\n",
+ "This Python code can be accelerated to Fortran or C, and even include OpenMP instructions, depending on the \"backend\" chosen by the user.\n",
+ "The available backends can be accessed through the `PSYDAC_BACKENDS` dictionary defined in the module `psydac.api.settings`, where the dictionary key is one the following strings: `['python', 'pyccel-gcc', 'pyccel-intel', 'pyccel-pgi', 'pyccel-nvidia']`.\n",
+ "Each of the dictionary values is itself a dictionary containing various options which can be modified by the user.\n",
+ "\n",
+ "`PSYDAC_BACKENDS['python']` defines a pure Python backend; this means that the generated Python assembly files will be used as they are.\n",
+ "The other options are based on accelerating the Python assembly files using the Pyccel transpiler and one of the compilers: GCC, Intel, PGI, or NVIDIA.\n",
+ "For all but the Python backend, OpenMP parallelization can be activated by setting `backend['omp'] = True`, and compilation flags can be accessed and changed using `backend['flags']`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.api.discretization import discretize\n",
+ "from psydac.api.settings import PSYDAC_BACKENDS\n",
+ "\n",
+ "# Modify to use MPI and/or OpenMP\n",
+ "use_mpi = False\n",
+ "use_openmp = False\n",
+ "\n",
+ "# MPI communicator, if any\n",
+ "if use_mpi:\n",
+ " from mpi4py import MPI\n",
+ " comm = MPI.COMM_WORLD\n",
+ "else:\n",
+ " comm = None\n",
+ "\n",
+ "# Choose backend\n",
+ "backend = PSYDAC_BACKENDS['pyccel-gcc']\n",
+ "if use_openmp:\n",
+ " import os\n",
+ " os.environ['OMP_NUM_THREADS'] = \"4\"\n",
+ " backend['omp'] = True \n",
+ "\n",
+ "# Discretization parameters\n",
+ "ncells = [9, 9]\n",
+ "degree = [3, 3]\n",
+ "periodic = [False, False]\n",
+ "nquads = [p + 1 for p in degree]\n",
+ "\n",
+ "# Discretize domain and function space\n",
+ "Omega_h = discretize(Omega, ncells=ncells, periodic=periodic)\n",
+ "V_h = discretize(V, Omega_h, degree=degree)\n",
+ "\n",
+ "# Discretize weak formulation and error norms.\n",
+ "# Python code is generated for the assembly kernels.\n",
+ "# The number of quadrature points should be provided.\n",
+ "# If a Pyccel backend is chosen, the Python code is translated to Fortran.\n",
+ "equation_h = discretize(equation, Omega_h, [V_h, V_h], nquads=nquads, backend=backend)\n",
+ "l2norm_h = discretize(l2norm, Omega_h, V_h, nquads=nquads, backend=backend)\n",
+ "h1norm_h = discretize(h1norm, Omega_h, V_h, nquads=nquads, backend=backend)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 4: Solve the equation and compute the error norms\n",
+ "\n",
+ "When calling the `assemble` method of a `DiscreteBilinearForm`, `DiscreteLinearForm`, or `DiscreteFunctional`, PSYDAC assembles a matrix, vector, or scalar quantity, respectively.\n",
+ "Objects of type `DiscreteEquation` do not have an `assemble` method, but they have the methods `set_solver` and `solve` instead:\n",
+ "- The `set_solver` method specifies the iterative algorithm (and its parameters) for solving the linear system $\\mathbb{A} \\mathrm{w} = \\mathrm{b}$ to obtain the spline coefficients $\\mathrm{w}$ of the Galerkin solution $u_h$.\n",
+ "- The `solve` methods performs a few steps.\n",
+ " First, it assembles both the left-hand-side matrix $\\mathbb{A}$ and right-hand-side vector $\\mathrm{b}$.\n",
+ " Then, it solves the linear system $\\mathbb{A} \\mathrm{w} = \\mathrm{b}$ with the chosen iterative algorithm.\n",
+ " Finally, it creates a callable `FemField` with spline coefficients $\\mathrm{w}$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import time\n",
+ "from psydac.fem.plotting_utilities import plot_field_2d\n",
+ "\n",
+ "# Set the solver parameters (CG = Conjugate Gradient)\n",
+ "equation_h.set_solver('CG', info=True, tol=1e-14)\n",
+ "\n",
+ "# Compute the numerical solution u_h:\n",
+ "# 1. assemble (distributed) sparse matrix A\n",
+ "# 2. assemble (distributed) dense vector b\n",
+ "# 3. solve linear system A w = b\n",
+ "# 4. create callable field u_h(x,y)\n",
+ "t0_s = time.time()\n",
+ "u_h, info = equation_h.solve()\n",
+ "t1_s = time.time()\n",
+ "\n",
+ "# Compute the L^2 and H^1 error norms\n",
+ "t0_d = time.time()\n",
+ "l2_error = l2norm_h.assemble(u=u_h)\n",
+ "h1_error = h1norm_h.assemble(u=u_h)\n",
+ "t1_d = time.time()\n",
+ "\n",
+ "# Print convergence information, error norms, and timings\n",
+ "print('> CG info :: ',info )\n",
+ "print('> L2 error :: {:.2e}'.format(l2_error))\n",
+ "print('> H1 error :: {:.2e}'.format(h1_error))\n",
+ "print('> Solution time :: {:.2e} s'.format(t1_s - t0_s))\n",
+ "print('> Evaluat. time :: {:.2e} s'.format(t1_d - t0_d))\n",
+ "\n",
+ "# Plot the results\n",
+ "plot_field_2d(u_h,\n",
+ " domain = Omega,\n",
+ " title = '$u_h(x,y)$: computed FEM solution',\n",
+ " plot_type = 'components',\n",
+ " N_vis = 50,\n",
+ " hide_plot = False,\n",
+ " filename = None,\n",
+ " cmap = 'jet',\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Testing the notebook"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import ipytest\n",
+ "ipytest.autoconfig(raise_on_error=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "%%ipytest\n",
+ "\n",
+ "def test_l2error():\n",
+ " assert l2_error < 8.68e-03\n",
+ "\n",
+ "def test_h1error():\n",
+ " assert h1_error < 1.35e-01"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/examples/notebooks/regularized_curlcurl_3D_VTK.ipynb b/examples/notebooks/regularized_curlcurl_3D_VTK.ipynb
new file mode 100644
index 000000000..e8ef51c29
--- /dev/null
+++ b/examples/notebooks/regularized_curlcurl_3D_VTK.ipynb
@@ -0,0 +1,276 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "2443e734",
+ "metadata": {},
+ "source": [
+ "# 3D example with VTK and ParaView \n",
+ "\n",
+ "In this notebook we show how to solve a 3D problem in parallel, save the solution in HDF5 format and export it to VTK. \n",
+ "\n",
+ "We solve the regularized curl-curl problem in 3D with a callable function for the source. In practice, the source could be the interpolation of some experimental data. Besides, we impose periodic boundary conditions everywhere in the boundary.\n",
+ "\n",
+ "Let $\\Omega=[0, 10] \\times [-2\\pi, 2\\pi]\\times [-3, 1]$. We want to find $u\\in H(\\mathrm{curl}, \\Omega)$ such that\n",
+ "\n",
+ "$$ \\mathrm{curl} \\mathrm{curl} u + \\mu u = f, $$\n",
+ "\n",
+ "where $\\mu = \\pi$, and $f$ is a callable function (of the domain coordinates) for which we may not have a symbolic expression. \n",
+ "\n",
+ "Since $\\mu>0$ the resulting operator is positive definite and we can use conjugate gradient to solve the system. In the case of the time-harmonic Maxwell problem we have that $\\mu<0$, which yields an indefinite operator and conjugate gradient is no longer guaranteed to converge.\n",
+ "\n",
+ "Only for illustration here we use the callable function\n",
+ "$$f(x,y,z) = (10e^{-0.5(x-4)^2 - 1.5z^2 - 0.1y^2}, 0, \\cos(y)\\sin(\\pi x / 5)),$$\n",
+ "\n",
+ "and project it from $H(\\mathrm{curl}, \\Omega)$ to the corresponding discretized space. This process is analogous with any callable vectorial function.\n",
+ "\n",
+ "\n",
+ "\n",
+ "## Step 1 : define the domain and spaces and discretize them\n",
+ "\n",
+ "This is similar to what is done in other examples"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "39c21151",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "from mpi4py import MPI\n",
+ "from sympde.topology import Cube\n",
+ "from sympde.topology import VectorFunctionSpace\n",
+ "from psydac.api.discretization import discretize\n",
+ "\n",
+ "comm = MPI.COMM_WORLD # MPI communicator\n",
+ "\n",
+ "Omega = Cube('Omega', bounds1=(0, 10), bounds2=(-2*np.pi, 2*np.pi), bounds3=(-3, 1)) # 3D cubic domain\n",
+ "\n",
+ "W = VectorFunctionSpace('W', Omega, kind='hcurl') # Hcurl space\n",
+ "\n",
+ "ncells = [25, 30, 11] # number of cells in each direction\n",
+ "degree = [3, 2, 2] # B-spline degree in each direction\n",
+ "\n",
+ "periodic = [True, True, True] # set periodic BC in every direction\n",
+ "\n",
+ "Omega_h = discretize(Omega, ncells=ncells, periodic=periodic, comm=comm)\n",
+ "Wh = discretize(W, Omega_h, degree=degree)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "36a9632b",
+ "metadata": {},
+ "source": [
+ "## Step 2 : compute the matrix\n",
+ "\n",
+ "We assemble the mass and curl curl matrices."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "091d4560",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympde.topology import elements_of, element_of\n",
+ "from sympde.expr.expr import BilinearForm, LinearForm, integral\n",
+ "from sympde.calculus import dot, curl\n",
+ "from psydac.api.settings import PSYDAC_BACKENDS\n",
+ "from psydac.fem.basic import FemField\n",
+ "\n",
+ "backend_language = 'pyccel-gcc'\n",
+ "backend = PSYDAC_BACKENDS[backend_language]\n",
+ "\n",
+ "u, v = elements_of(W, names='u, v') #trial and test functions\n",
+ "\n",
+ "m = BilinearForm((u, v), integral(Omega, dot(u, v)))\n",
+ "mh = discretize(m, Omega_h, [Wh, Wh], backend=backend)\n",
+ "M = mh.assemble() # Mass matrix\n",
+ "\n",
+ "k = BilinearForm((u, v), integral(Omega, dot(curl(u), curl(v))))\n",
+ "kh = discretize(k, Omega_h, [Wh, Wh], backend=backend)\n",
+ "K = kh.assemble() # Curl curl matrix\n",
+ "\n",
+ "mu = np.pi\n",
+ "\n",
+ "A = K + mu*M # system matrix, positive-definite"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8927f88c",
+ "metadata": {},
+ "source": [
+ "## Step 3 : compute the right-hand size\n",
+ "\n",
+ "We project a callable function to the vector space and use it to assemble the right-hand side of the system."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a13044fa",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.feec.global_geometric_projectors import GlobalGeometricProjectorHcurl\n",
+ "f_callable = [lambda x, y, z: 10*np.exp(-0.5*(x-4)**2 - 1.5*z**2 - 1e-1*y**2),\n",
+ " lambda x, y, z: 0,\n",
+ " lambda x, y, z: np.cos(y)*np.sin(x*np.pi/5)]\n",
+ "\n",
+ "proj_Wh = GlobalGeometricProjectorHcurl(Wh, nquads=[p+1 for p in degree]) # get Hcurl projector, nquads is the number of points for the Gauss quadrature.\n",
+ "fh = proj_Wh(f_callable) # project callable to field in Wh\n",
+ "\n",
+ "f = element_of(W, name='f') # free variable for the FemField\n",
+ "\n",
+ "l = LinearForm(v, integral(Omega, dot(f,v)))\n",
+ "lh = discretize(l, Omega_h, Wh, backend=backend)\n",
+ "rhs = lh.assemble(f=fh) # assmble RHS, f is a free parameter"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8258971a",
+ "metadata": {},
+ "source": [
+ "## Step 4 : solve the system\n",
+ "Solve the system using Conjugate Gradient. Notice that since $\\mu>0$, the matrix is positive definite."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1f0f103d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from psydac.linalg.solvers import inverse\n",
+ "\n",
+ "inv_A = inverse(A, solver='cg', tol=1e-2, maxiter=100) # invert using Conjugate Gradient\n",
+ "\n",
+ "u = inv_A @ rhs # array of coefficients\n",
+ "\n",
+ "uh = FemField(Wh, coeffs=u) # finite element field, callable"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "e0c73790",
+ "metadata": {},
+ "source": [
+ "## Step 5 : Save results to HDF5 format\n",
+ "We save the coefficients of the solution and the field used in the right-hand size. This step generates a YAML file containing the space information, and a HDF5 file that stores the coefficients. The fields are saved in parallel. This step does not include an evaluation of the fields on a grid."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6c77d25a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Save the results using OutputManager\n",
+ "from psydac.api.postprocessing import OutputManager\n",
+ "import os\n",
+ "\n",
+ "os.makedirs('results_regularized_curlcurl', exist_ok=True)\n",
+ "\n",
+ "Om = OutputManager(\n",
+ " f'results_regularized_curlcurl/space_info_{Omega.name}.yml',\n",
+ " f'results_regularized_curlcurl/field_info_{Omega.name}.h5',\n",
+ " comm=comm,\n",
+ " save_mpi_rank=True, \n",
+ " mode = 'w' \n",
+ ")\n",
+ "\n",
+ "Om.add_spaces(V=Wh)\n",
+ "Om.export_space_info()\n",
+ "\n",
+ "Om.set_static() # the fields do not depend on time\n",
+ "\n",
+ "save_fields = {\n",
+ " 'u': uh,\n",
+ " 'f': fh\n",
+ "}\n",
+ "\n",
+ "Om.export_fields(**save_fields) # saves the coefficients of the fields\n",
+ "\n",
+ "Om.close()\n",
+ "\n",
+ "# This generates a YAML file containing the space information, and a HDF5 file that stores the solution coefficients. \n",
+ "# At this point the solution has not been evaluated on any grid."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a619a838",
+ "metadata": {},
+ "source": [
+ "## Step 6 : Evaluate fields on a grid and export to VTK\n",
+ "Evaluate the saved fields on a grid and export the result to VTK format."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "86a1c531",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Export the results to VTK using PostProcessManager\n",
+ "from psydac.api.postprocessing import PostProcessManager\n",
+ "\n",
+ "Pm = PostProcessManager(\n",
+ " domain=Omega,\n",
+ " space_file=f'results_regularized_curlcurl/space_info_{Omega.name}.yml',\n",
+ " fields_file=f'results_regularized_curlcurl/field_info_{Omega.name}.h5',\n",
+ " comm=comm\n",
+ ")\n",
+ "\n",
+ "Pm.export_to_vtk(\n",
+ " f'results_regularized_curlcurl/visu_{Omega.name}',\n",
+ " grid=None,\n",
+ " npts_per_cell=[4,2,4],\n",
+ " fields=save_fields.keys() # evaluate all the fields that were saved\n",
+ ")\n",
+ "\n",
+ "Pm.close()\n",
+ "# If this file is run in sequential, then this step will generate a single .vtu file that can be opened with ParaView.\n",
+ "# If it is run in parallel, it generates a .vtu file per process (containing the local grid) and a .pvtu file. To visualize it in ParaView, open the .pvtu file."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "af41d09b",
+ "metadata": {},
+ "source": [
+ "## Visualization of the solution in ParaView\n",
+ "\n",
+ " ![](\"./paraview_images/regularized_curlcurl_3D_ux.png\")
\n",
+ " ![](\"./paraview_images/regularized_curlcurl_3D_uy.png\")
\n",
+ " ![](\"./paraview_images/regularized_curlcurl_3D_uz.png\")
\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "0bb786a7",
+ "metadata": {},
+ "source": [
+ "## Visualization of the source in ParaView\n",
+ "\n",
+ " ![](\"./paraview_images/regularized_curlcurl_3D_fx.png\")
\n",
+ " ![](\"./paraview_images/regularized_curlcurl_3D_fy.png\")
\n",
+ " ![](\"./paraview_images/regularized_curlcurl_3D_fz.png\")
\n",
+ "
"
+ ]
+ }
+ ],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/performance/compare_3d_matrix_assembly_speed.py b/examples/performance/compare_3d_matrix_assembly_speed.py
similarity index 100%
rename from performance/compare_3d_matrix_assembly_speed.py
rename to examples/performance/compare_3d_matrix_assembly_speed.py
diff --git a/performance/matrix_assembly_speed_log.md b/examples/performance/matrix_assembly_speed_log.md
similarity index 100%
rename from performance/matrix_assembly_speed_log.md
rename to examples/performance/matrix_assembly_speed_log.md
diff --git a/feectools/api/ast/tests/test_poisson.py b/feectools/api/ast/tests/test_poisson.py
new file mode 100644
index 000000000..00778483d
--- /dev/null
+++ b/feectools/api/ast/tests/test_poisson.py
@@ -0,0 +1,91 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import pytest
+import sys
+
+# Skip entire module if sympde is not available
+sympde_available = True
+try:
+ import sympde
+except ImportError:
+ sympde_available = False
+
+if not sympde_available:
+ pytest.skip("sympde not installed", allow_module_level=True)
+
+import os
+from pathlib import Path
+
+from sympy import symbols
+from sympy import sin, pi
+
+from sympde.calculus import grad, dot
+from sympde.topology import ScalarFunctionSpace
+from sympde.topology import elements_of, LogicalExpr
+from sympde.topology import Square
+from sympde.topology import Mapping, IdentityMapping, PolarMapping
+from sympde.expr import integral
+from sympde.expr import LinearForm
+from sympde.expr import BilinearForm
+from sympde.expr import Norm, SemiNorm
+from sympde.expr.evaluation import TerminalExpr
+
+from psydac.api.ast.fem import AST
+from psydac.api.ast.parser import parse
+from psydac.api.discretization import discretize
+from psydac.api.printing.pycode import pycode
+from psydac.api.settings import PSYDAC_BACKENDS
+
+#==============================================================================
+# Get the mesh directory
+import psydac.cad.mesh as mesh_mod
+
+mesh_dir = Path(mesh_mod.__file__).parent
+filename = os.path.join(mesh_dir, 'identity_2d.h5')
+
+# Choose backend
+backend = PSYDAC_BACKENDS['python']
+
+#==============================================================================
+def test_codegen():
+ domain = Square()
+ M = Mapping('M',2)
+ domain = M(domain)
+ V = ScalarFunctionSpace('V', domain)
+ u,v = elements_of(V, names='u,v')
+
+ x,y = symbols('x, y')
+
+ b = LogicalExpr(BilinearForm((u,v), integral(domain, dot(grad(u), grad(v)))), domain)
+ l = LogicalExpr(LinearForm(v, integral(domain, v*2*pi**2*sin(pi*x)*sin(pi*y))), domain)
+
+ # Create computational domain from topological domain
+ domain_h = discretize(domain, filename=filename)
+
+ # Discrete spaces
+ Vh = discretize(V, domain_h)
+
+ error = u - sin(pi*x)*sin(pi*y)
+ l2norm = LogicalExpr( Norm(error, domain, kind='l2'), domain)
+ h1norm = LogicalExpr(SemiNorm(error, domain, kind='h1'), domain)
+
+ ast_b = AST(b, TerminalExpr(b, domain)[0], [Vh, Vh], nquads=(2, 2), backend=backend)
+ ast_b = parse(ast_b.expr, settings={'dim':2, 'nderiv':1, 'mapping':M, 'target':domain.logical_domain}, backend=backend)
+ print(pycode(ast_b))
+
+ print('==============================================================================================================')
+ ast_l = AST(l, TerminalExpr(l, domain)[0], Vh, nquads=(2, 2), backend=backend)
+ ast_l = parse(ast_l.expr, settings={'dim':2, 'nderiv':1, 'mapping':M, 'target':domain.logical_domain}, backend=backend)
+ print(pycode(ast_l))
+
+ print('==============================================================================================================')
+ ast_l2norm = AST(l2norm, TerminalExpr(l2norm, domain)[0], Vh, nquads=(2, 2), backend=backend)
+ ast_l2norm = parse(ast_l2norm.expr, settings={'dim':2, 'nderiv':1, 'mapping':M, 'target':domain.logical_domain}, backend=backend)
+ print(pycode(ast_l2norm))
+
+#==============================================================================
+if __name__ == '__main__':
+ test_codegen()
diff --git a/feectools/api/tests/test_sum_factorization_assembly_3d.py b/feectools/api/tests/test_sum_factorization_assembly_3d.py
new file mode 100644
index 000000000..e0b8385f9
--- /dev/null
+++ b/feectools/api/tests/test_sum_factorization_assembly_3d.py
@@ -0,0 +1,600 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import os
+from pathlib import Path
+
+import pytest
+import time
+
+# Skip entire module if sympde is not available
+try:
+ import sympde
+except ImportError:
+ pytest.skip("sympde not installed", allow_module_level=True)
+
+import numpy as np
+from sympy import sin, sqrt, pi, Abs, cos, tan
+
+from sympde.calculus import dot, cross, grad, curl, div, laplace
+from sympde.expr import BilinearForm, integral
+from sympde.topology import element_of, elements_of, Cube, Mapping, ScalarFunctionSpace, VectorFunctionSpace, Domain, Derham
+
+from psydac.api.discretization import discretize
+from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL
+from psydac.linalg.block import BlockVectorSpace
+from psydac.fem.basic import FemField
+
+# Get the mesh directory
+import psydac.cad.mesh as mesh_mod
+mesh_dir = Path(mesh_mod.__file__).parent
+
+# With PR #448, matrices corresponding to bilinear forms on 3D domains are being assembled using a so called sum factorization algorithm.
+# Unless explicitely using the old algorithm, this happens automatically, and hence all old tests passing should indicate that the implementation of the sum factorization algorithm has been successful.
+# Nonetheless, there are various difficulties in the implementation, and possibly not all of them are accounted for in the existing tests.
+
+# This file is designed to test such "difficult" edge cases - for mapped (Bspline & analytical) and parametric domains:
+
+# Such "difficult" edge cases are:
+
+# 1. bilinear forms on different spaces
+# 2. (FemField / analytical / ...) weight functions
+# 3. high derivatives (>=2)
+# 4. (multiple) free FemFields
+# 5. complicated expressions
+# 6. uncommen (numpy) functions that need be imported correctly in the assembly file (here in the weight function)
+
+# These tests also return old and new discretization and assembly times.
+
+# Most of the time, being close to the "old matrix" (generated using the old assembly algorithm) will be the requirement to pass a test,
+# as the old implementation has not caused problems in a long time and is considered to function properly.
+
+# Update: Instead of testing all mapping options all of the time, we now rather randomly test one of the three options!
+#@pytest.mark.parametrize('mapping', ('None', 'Analytical', 'Bspline'))
+def test_assembly(): # mapping):
+
+ rng = np.random.default_rng() # (seed=42)
+ mapping_options = ['None', 'Analytical', 'Bspline']
+ mapping = mapping_options[int(np.floor(rng.random()*3))]
+
+ ncells = [7, 5, 6]
+ degree = [2, 4, 3]
+ periodic = [int(np.floor(rng.random()*2))==True for _ in range(3)]
+ print(f'Random periodicity: {periodic}')
+
+ trial_multiplicity = [1, 3, 2]
+ test_multiplicity = [2, 2, 3]
+
+ backend = PSYDAC_BACKEND_GPYCCEL
+
+ if mapping == 'None':
+
+ domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1))
+ derham = Derham(domain)
+
+ domain_h = discretize(domain, ncells=ncells, periodic=periodic)
+ derham_h = discretize(derham, domain_h, degree=degree, multiplicity=trial_multiplicity)
+ derham_test_h = discretize(derham, domain_h, degree=degree, multiplicity=test_multiplicity)
+
+ elif mapping == 'Bspline':
+
+ filename = os.path.join(mesh_dir, 'identity_3d.h5')
+
+ domain = Domain.from_file(filename=filename)
+ derham = Derham(domain)
+
+ domain_h = discretize(domain, filename=filename)
+ derham_h = discretize(derham, domain_h, degree=domain.mapping.get_callable_mapping().space.degree, multiplicity=trial_multiplicity)
+ derham_test_h = discretize(derham, domain_h, degree=domain.mapping.get_callable_mapping().space.degree, multiplicity=test_multiplicity)
+
+ elif mapping == 'Analytical':
+
+ class HalfSquareTorusMapping3D(Mapping):
+ _expressions = {'x': 'x1 * cos(x2)',
+ 'y': 'x1 * sin(x2)',
+ 'z': 'x3'}
+
+ _ldim = 3
+ _pdim = 3
+
+ M = HalfSquareTorusMapping3D('M')
+ logical_domain = Cube('C', bounds1=(0.3,1), bounds2=(0,np.pi), bounds3=(0,1))
+
+ domain = M(logical_domain)
+ derham = Derham(domain)
+
+ domain_h = discretize(domain, ncells=ncells, periodic=periodic)
+ derham_h = discretize(derham, domain_h, degree=degree, multiplicity=trial_multiplicity)
+ derham_test_h = discretize(derham, domain_h, degree=degree, multiplicity=test_multiplicity)
+
+ x, y, z = domain.coordinates
+ weight = 1 + sqrt(Abs(x*y**2 + z)) + abs(sin(x-2*y + pi + np.pi))
+
+ V0 = derham.V0
+ V0h = derham_h.V0
+ V1 = derham.V1
+ V1h = derham_h.V1
+ V2 = derham.V2
+ V2h = derham_h.V2
+ V3 = derham.V3
+ V3h = derham_h.V3
+
+ V0h_test = derham_test_h.V0
+ V1h_test = derham_test_h.V1
+ V2h_test = derham_test_h.V2
+ V3h_test = derham_test_h.V3
+
+ Vs = ScalarFunctionSpace('Vsh', domain)
+ Vsh = discretize(Vs, domain_h, degree=degree)
+ Vvc = VectorFunctionSpace('Vvc', domain, kind='hcurl')
+ Vvch = discretize(Vvc, domain_h, degree=degree)
+ Vvd = VectorFunctionSpace('Vvd', domain, kind='hdiv')
+ Vvdh = discretize(Vvd, domain_h, degree=degree)
+
+ Vsh_test = discretize(Vs, domain_h, degree=degree, multiplicity=test_multiplicity)
+ Vvch_test = discretize(Vs, domain_h, degree=degree, multiplicity=test_multiplicity)
+ Vvdh_test = discretize(Vs, domain_h, degree=degree, multiplicity=test_multiplicity)
+
+ u1, u2, F01, F02, F03 = elements_of(V0, names='u1, u2, F01, F02, F03')
+ v1, v2, F11, F12, F13 = elements_of(V1, names='v1, v2, F11, F12, F13')
+ w1, w2, F21, F22, F23 = elements_of(V2, names='w1, w2, F21, F22, F23')
+ f1, f2, F31, F32, F33 = elements_of(V3, names='f1, f2, F31, F32, F33')
+
+ fs1, fs2, Fs1, Fs2, Fs3 = elements_of(Vs, names='fs1, fs2, Fs1, Fs2, Fs3')
+ fvc1, fvc2, Fvc1, Fvc2, Fvc3 = elements_of(Vvc, names='fvc1, fvc2, Fvc1, Fvc2, Fvc3')
+ fvd1, fvd2, Fvd1, Fvd2, Fvd3 = elements_of(Vvd, names='fvd1, fvd2, Fvd1, Fvd2, Fvd3')
+
+ trial_spaces = {'V0': {'Vh':V0h, 'funcs':[u1, u2]},
+ 'V1': {'Vh':V1h, 'funcs':[v1, v2]},
+ 'V2': {'Vh':V2h, 'funcs':[w1, w2]},
+ 'V3': {'Vh':V3h, 'funcs':[f1, f2]},
+ 'Vs': {'Vh':Vsh, 'funcs':[fs1, fs2]},
+ 'Vvc':{'Vh':Vvch, 'funcs':[fvc1, fvc2]},
+ 'Vvd':{'Vh':Vvdh, 'funcs':[fvd1, fvd2]}}
+
+ test_spaces = {'V0': {'Vh':V0h_test, 'funcs':[u1, u2]},
+ 'V1': {'Vh':V1h_test, 'funcs':[v1, v2]},
+ 'V2': {'Vh':V2h_test, 'funcs':[w1, w2]},
+ 'V3': {'Vh':V3h_test, 'funcs':[f1, f2]},
+ 'Vs': {'Vh':Vsh_test, 'funcs':[fs1, fs2]},
+ 'Vvc':{'Vh':Vvch_test, 'funcs':[fvc1, fvc2]},
+ 'Vvd':{'Vh':Vvdh_test, 'funcs':[fvd1, fvd2]}}
+
+ F01_coeffs = V0h.coeff_space.zeros()
+ rng.random(size=F01_coeffs._data.shape, dtype='float64', out=F01_coeffs._data)
+ F11_coeffs = V1h.coeff_space.zeros()
+ for block in F11_coeffs.blocks:
+ rng.random(size=block._data.shape, dtype='float64', out=block._data)
+ F12_coeffs = V1h.coeff_space.zeros()
+ for block in F12_coeffs.blocks:
+ rng.random(size=block._data.shape, dtype='float64', out=block._data)
+ F21_coeffs = V2h.coeff_space.zeros()
+ for block in F21_coeffs.blocks:
+ rng.random(size=block._data.shape, dtype='float64', out=block._data)
+ Fvc1_coeffs = Vvch.coeff_space.zeros()
+ for block in Fvc1_coeffs.blocks:
+ rng.random(size=block._data.shape, dtype='float64', out=block._data)
+ Fs1_coeffs = Vsh.coeff_space.zeros()
+ rng.random(size=Fs1_coeffs._data.shape, dtype='float64', out=Fs1_coeffs._data)
+ Fs2_coeffs = Vsh.coeff_space.zeros()
+ rng.random(size=Fs2_coeffs._data.shape, dtype='float64', out=Fs2_coeffs._data)
+ Fvd1_coeffs = Vvdh.coeff_space.zeros()
+ for block in Fvd1_coeffs.blocks:
+ rng.random(size=block._data.shape, dtype='float64', out=block._data)
+
+ F01_field = FemField(V0h, F01_coeffs)
+ F11_field = FemField(V1h, F11_coeffs)
+ F12_field = FemField(V1h, F12_coeffs)
+ F21_field = FemField(V2h, F21_coeffs)
+ Fvc1_field = FemField(Vvch, Fvc1_coeffs)
+ Fs1_field = FemField(Vsh, Fs1_coeffs)
+ Fs2_field = FemField(Vsh, Fs2_coeffs)
+ Fvd1_field = FemField(Vvdh, Fvd1_coeffs)
+
+ bilinear_forms = { # one and two free FemFields without derivatives (with derivatives in seperate test)
+ # complicated expressions
+ 'Q' :{'trial' :'V1', 'test':'V1',
+ 'expr' :dot(cross(F11, v1), cross(F11, v2)),
+ 'fields':[F11_field, ]},
+ 'equilibrium' :{'trial' :'V1', 'test':'V1',
+ 'expr' :dot(cross(F11, v1), cross(v2, F12)),
+ 'fields':[F11_field, F12_field]},
+ 'Elena' :{'trial' :'V1', 'test':'V1',
+ 'expr' :dot(F01*v1, v2),
+ 'fields':[F01_field, ]},
+ # weight function, free FemField, different spaces
+ 'dot(grad(u),v)':{'trial' :'V0', 'test':'V1',
+ 'expr' :dot(grad(u1), v2)*F01*weight,
+ 'fields':[F01_field, ]},
+ 'dot(curl(v),w)':{'trial' :'V1', 'test':'V2',
+ 'expr' :dot(curl(v1), F21)*div(w2)*weight,
+ 'fields':[F21_field, ]},
+ # among other difficulties: multiple scalar FemFields
+ 'ScalarFields' :{'trial' :'V0', 'test':'V1',
+ 'expr' :dot(grad(u1), curl(Fvc1))*dot(grad(Fs1), curl(v2))*Fs2*div(Fvd1),
+ 'fields':[Fvc1_field, Fs1_field, Fs2_field, Fvd1_field]},
+ # high derivatives, not FEEC
+ 'bilaplace' :{'trial' :'Vs', 'test':'Vs',
+ 'expr' :laplace(fs1)*laplace(fs2)}
+ }
+
+ # test all BFs
+ # bilinear_form_strings_to_test = list(bilinear_forms.keys())
+
+ # or only a subset
+ bilinear_form_strings_to_test = list(bilinear_forms.keys())[:-1] # exclude expensive bilaplace test
+
+ bilinear_forms_to_test = {}
+ for name in bilinear_form_strings_to_test:
+ bilinear_forms_to_test[name] = bilinear_forms[name]
+
+ int_0 = lambda expr: integral(domain, expr)
+ print()
+
+ for bf_name, bf_data in bilinear_forms_to_test.items():
+
+ trial_space = trial_spaces[bf_data['trial']]
+ Vh = trial_space['Vh']
+ u = trial_space['funcs'][0]
+ test_space = test_spaces[bf_data['test']]
+ Wh = test_space ['Vh']
+ v = test_space['funcs'][1]
+ expr = bf_data['expr']
+ if 'fields' in bf_data.keys():
+ fields = bf_data['fields']
+
+ a = BilinearForm((u, v), int_0(expr))
+
+ t0 = time.time()
+ ah_old = discretize(a, domain_h, (Vh, Wh), backend=backend, sum_factorization=False)
+ t1 = time.time()
+ discretization_time_old = t1 - t0
+
+ t0 = time.time()
+ ah = discretize(a, domain_h, (Vh, Wh), backend=backend)
+ t1 = time.time()
+ discretization_time = t1 - t0
+
+ if bf_name == 'Q':
+ t0_old = time.time()
+ A_old = ah_old.assemble(F11=fields[0])
+ t1_old = time.time()
+
+ t0 = time.time()
+ A = ah.assemble(F11=fields[0])
+ t1 = time.time()
+ elif bf_name == 'equilibrium':
+ t0_old = time.time()
+ A_old = ah_old.assemble(F11=fields[0], F12=fields[1])
+ t1_old = time.time()
+
+ t0 = time.time()
+ A = ah.assemble(F11=fields[0], F12=fields[1])
+ t1 = time.time()
+ elif bf_name in ('Elena', 'dot(grad(u),v)'):
+ t0_old = time.time()
+ A_old = ah_old.assemble(F01=fields[0])
+ t1_old = time.time()
+
+ t0 = time.time()
+ A = ah.assemble(F01=fields[0])
+ t1 = time.time()
+ elif bf_name == 'dot(curl(v),w)':
+ t0_old = time.time()
+ A_old = ah_old.assemble(F21=fields[0])
+ t1_old = time.time()
+
+ t0 = time.time()
+ A = ah.assemble(F21=fields[0])
+ t1 = time.time()
+ elif bf_name == 'ScalarFields':
+ t0_old = time.time()
+ A_old = ah_old.assemble(Fvc1=fields[0], Fs1=fields[1], Fs2=fields[2], Fvd1=fields[3])
+ t1_old = time.time()
+
+ t0 = time.time()
+ A = ah.assemble(Fvc1=fields[0], Fs1=fields[1], Fs2=fields[2], Fvd1=fields[3])
+ t1 = time.time()
+ else:
+ t0_old = time.time()
+ A_old = ah_old.assemble()
+ t1_old = time.time()
+
+ t0 = time.time()
+ A = ah.assemble()
+ t1 = time.time()
+
+ assembly_time_old = t1_old - t0_old
+ assembly_time = t1 - t0
+
+ # Testing whether two linear operators are identical by comparing their arrays is quite expensive.
+ # Thus we instead test whether three random domain vectors applied to both
+ # the old and the new matrix produce the same codomain vector.
+
+ domain_vector1 = Vh.coeff_space.zeros()
+ domain_vector2 = Vh.coeff_space.zeros()
+ domain_vector3 = Vh.coeff_space.zeros()
+ domain_vectors = [domain_vector1, domain_vector2, domain_vector3]
+
+ if isinstance(Vh.coeff_space, BlockVectorSpace):
+ for domain_vector in domain_vectors:
+ for block in domain_vector.blocks:
+ rng.random(size=block._data.shape, dtype='float64', out=block._data)
+ else:
+ for domain_vector in domain_vectors:
+ rng.random(size=domain_vector._data.shape, dtype='float64', out=domain_vector._data)
+
+ err = []
+ rel_err = []
+
+ for domain_vector in domain_vectors:
+ codomain_vector = A @ domain_vector
+ codomain_vector_old = A_old @ domain_vector
+
+ norm_old = np.sqrt(codomain_vector_old.inner(codomain_vector_old))
+
+ diff = codomain_vector - codomain_vector_old
+
+ err.append(np.sqrt(diff.inner(diff)))
+ rel_err.append(err[-1] / norm_old)
+
+ print(f' >>> Mapping: {mapping}')
+ print(f' >>> BilinearForm: {bf_name}')
+ print(f' >>> Discretization in: Old {discretization_time_old:.3g} \t\t || New {discretization_time:.3g} \t\t || Old/New {discretization_time_old/discretization_time:.3g}')
+ print(f' >>> Assembly in: Old {assembly_time_old:.3g} \t \t || New {assembly_time:.3g} \t\t || Old/New {assembly_time_old/assembly_time:.3g}')
+ print(f' >>> Error: {max(err):.3g}')
+ print(f' >>> Rel. Error: {max(rel_err):.3g}')
+ print()
+
+ assert max(rel_err) < 1e-12 # arbitrary rel. error bound (How to test better?)
+
+
+# fixed by PR #507
+#@pytest.mark.xfail
+def test_allocate_matrix_bug():
+ """
+ This test is related to Issue #504.
+
+ The bilinear form
+ (V0 x V3) ni (u, f) mapsto int_{Omega} u*f
+ should be the transpose of the bilinear form
+ (V3 x V0) ni (f, u) mapsto int_{Omega} u*f
+ but is not.
+ """
+
+ ncells = [15, 16, 17]
+ degree = [4, 3, 2]
+ periodic = [False, True, False]
+
+ backend = PSYDAC_BACKEND_GPYCCEL
+
+ domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1))
+ derham = Derham(domain)
+
+ domain_h = discretize(domain, ncells=ncells, periodic=periodic)
+ derham_h = discretize(derham, domain_h, degree=degree)
+
+ P0, _, _, P3 = derham_h.projectors()
+
+ V0 = derham.V0
+ V0h = derham_h.V0
+ V3 = derham.V3
+ V3h = derham_h.V3
+
+ u = element_of(V0, name='u')
+ f = element_of(V3, name='f')
+
+ fun = lambda x, y, z : 1
+ u_coeffs = P0(fun).coeffs
+ f_coeffs = P3(fun).coeffs
+
+ a0 = BilinearForm((u, f), integral(domain, u*f))
+ a1 = BilinearForm((f, u), integral(domain, u*f))
+
+ a0h = discretize(a0, domain_h, (V0h, V3h), backend=backend, sum_factorization=False)
+ a1h = discretize(a1, domain_h, (V3h, V0h), backend=backend, sum_factorization=False)
+
+ A0 = a0h.assemble()
+ A1 = a1h.assemble()
+ A1T = A1.T
+
+ # Clearly, it should hold A1T = A0, and further ||A0|| = ||A1||.
+ A0arr = A0.toarray()
+ A1arr = A1.toarray()
+ A1Tarr = A1T.toarray()
+
+ diff1 = np.linalg.norm(A0arr - A1Tarr)
+ diff2 = np.linalg.norm(A0arr) - np.linalg.norm(A1arr)
+
+ print(f' || A0 - A1.T || = {diff1:.3g}')
+ print(f' ||A0|| - ||A1|| = {diff2:.3g}')
+
+ # Further, the following integral should evaluate to 1.
+ # This however is only the case for the second integral,
+ # independent on whether one uses the new or old assembly algorithm.
+
+ print(f' 1 =? {A0.dot_inner(u_coeffs, f_coeffs)}')
+ print(f' 1 =? {A1.dot_inner(f_coeffs, u_coeffs)}')
+
+ assert diff1 <= 1e-12 # arbitrary error bound
+ assert diff2 <= 1e-12 # arbitrary error bound
+
+#@pytest.mark.xfail
+def test_free_FemField_derivatives():
+ """
+ These particular bilinear forms, when using a constant 1-vector coefficient vector for the free FemFields,
+ causes problems in a different test file of mine.
+ In particular, the assembled matrices used to have really small norms (~e-13).
+ That is probably due to the constant 1-vector corresponding to a constant function,
+ which means that all appearing derivatives of free FemFields are 0.
+
+ When using meaningful free FemFields, these dubious observations disappeared.
+
+ """
+
+ ncells = [5, 2, 4]
+ degree = [2, 1, 3]
+ periodic = [False, False, False]
+
+ backend = PSYDAC_BACKEND_GPYCCEL
+
+ domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1))
+ derham = Derham(domain)
+
+ domain_h = discretize(domain, ncells=ncells, periodic=periodic)
+ derham_h = discretize(derham, domain_h, degree=degree)
+
+ P0, P1, P2, P3 = derham_h.projectors()
+
+ V0 = derham.V0
+ V0h = derham_h.V0
+ V1 = derham.V1
+ V1h = derham_h.V1
+ V2 = derham.V2
+ V2h = derham_h.V2
+
+ u = element_of (V0, name= 'u')
+ v, F1 = elements_of(V1, names='v, F1')
+ w1, w2, F2 = elements_of(V2, names='w1, w2, F2')
+
+ u_func = lambda x, y, z: x + y + z
+ u_coeffs = P0(u_func).coeffs
+
+ v_1 = lambda x, y, z: 1
+ v_2 = lambda x, y, z: 1
+ v_3 = lambda x, y, z: 1
+ v_func = (v_1, v_2, v_3)
+ v_coeffs = P1(v_func).coeffs
+
+ w_1 = lambda x, y, z: x
+ w_2 = lambda x, y, z: y
+ w_3 = lambda x, y, z: z
+ w_func = (w_1, w_2, w_3)
+ w_coeffs = P2(w_func).coeffs
+
+ dubious_observations = False
+
+ if dubious_observations:
+ F1_coeffs = V1h.coeff_space.zeros()
+ F2_coeffs = V2h.coeff_space.zeros()
+ for block in F1_coeffs.blocks:
+ block._data = np.ones(block._data.shape, dtype='float64')
+ for block in F2_coeffs.blocks:
+ block._data = np.ones(block._data.shape, dtype='float64')
+ F1_FF = FemField(V1h, F1_coeffs)
+ F2_FF = FemField(V2h, F2_coeffs)
+ else:
+ F1_1 = lambda x, y, z: z
+ F1_2 = lambda x, y, z: x
+ F1_3 = lambda x, y, z: y
+ F1_func = (F1_1, F1_2, F1_3)
+ F1_FF = P1(F1_func)
+
+ F2_FF = FemField(V2h, w_coeffs.copy())
+
+ # with the above choices (dubious_observation = False):
+ # a0 reduces to 3*int_{Omega}x+y+z with Omega being the unit square. Expected value: 4.5
+ a0 = BilinearForm((u, v), integral(domain, dot(grad(u), curl(F1)) * dot(v, F1)))
+ # a1 reduces to 9* ----------------------------------------- " -------------------- 13.5
+ a1 = BilinearForm((u, w2), integral(domain, dot(grad(u), F2)*div(w2)*div(F2)))
+ # a2 reduces to 3* ----------------------------------------- " --------------------- 4.5
+ a2 = BilinearForm((w1, w2), integral(domain, dot(curl(F1), w1)*div(w2)))
+
+ a0h_old = discretize(a0, domain_h, (V0h, V1h), backend=backend, sum_factorization=False)
+ a1h_old = discretize(a1, domain_h, (V0h, V2h), backend=backend, sum_factorization=False)
+ a2h_old = discretize(a2, domain_h, (V2h, V2h), backend=backend, sum_factorization=False)
+
+ a0h = discretize(a0, domain_h, (V0h, V1h), backend=backend)
+ a1h = discretize(a1, domain_h, (V0h, V2h), backend=backend)
+ a2h = discretize(a2, domain_h, (V2h, V2h), backend=backend)
+
+ bfs = [(a0h_old, a0h), (a1h_old, a1h), (a2h_old, a2h)]
+ print()
+
+ for i, (ah_old, ah) in enumerate(bfs):
+
+ if i in (0, 2):
+ A_old = ah_old.assemble(F1=F1_FF)
+ A = ah.assemble(F1=F1_FF)
+ else:
+ A_old = ah_old.assemble(F2=F2_FF)
+ A = ah.assemble(F2=F2_FF)
+
+ if i == 0:
+ value_old = A_old.dot_inner(u_coeffs, v_coeffs)
+ value = A.dot_inner(u_coeffs, v_coeffs)
+ elif i == 1:
+ value_old = A_old.dot_inner(u_coeffs, w_coeffs)
+ value = A.dot_inner(u_coeffs, w_coeffs)
+ else:
+ value_old = A_old.dot_inner(w_coeffs, w_coeffs)
+ value = A.dot_inner(w_coeffs, w_coeffs)
+
+ A_old_arr = A_old.toarray()
+ A_arr = A.toarray()
+ A_old_norm = np.linalg.norm(A_old_arr)
+ A_norm = np.linalg.norm(A_arr)
+
+ err = np.linalg.norm(A_old_arr - A_arr)
+ rel_err = err / A_old_norm
+
+ print(f' i = {i}')
+ print(f' >>> Error: {err:.3g}')
+ print(f' >>> Rel. Error: {rel_err:.3g}')
+ print(f' >>> Norms: ||A_old|| = {A_old_norm:.3g} \t\t ||A|| = {A_norm:.3g}')
+ print()
+
+ # arbitrary tolerance
+ tol = 1e-12
+ if not dubious_observations:
+ assert abs(value-value_old) < tol
+ assert rel_err < tol
+
+def test_assembly_free_FemFields():
+
+ backend = PSYDAC_BACKEND_GPYCCEL
+
+ domain = Cube(bounds1=(2626,3179), bounds2=(-138, 138), bounds3=(-760.3, 69))
+ derham = Derham(domain)
+ V0 = derham.V0
+ V1 = derham.V1
+
+ u = element_of(V1, name='u')
+ v = element_of(V1, name='v')
+ p = element_of(V0, name='p')
+
+ p_call = lambda xi,yi,zi: np.cos(2*np.pi*(xi-2626)/553) + np.tan(2*np.pi*(zi+760.3)/(5*829.3))
+ x,y,z = domain.coordinates
+ p_sym = cos(2*np.pi*(x-2626)/553) + tan(2*np.pi*(z+760.3)/(5*829.3))
+
+ a_sym = BilinearForm((u, v), integral(domain, p_sym*dot(u, v)))
+ a_fem = BilinearForm((u, v), integral(domain, p*dot(u, v)))
+
+ ncells = (11, 1, 17)
+ degree = (3, 1, 4)
+ domain_h = discretize(domain, ncells=ncells, periodic=(False, True, False))
+ derham_h = discretize(derham, domain_h, degree=degree)
+ V1_h = derham_h.V1
+ P0 = derham_h.projectors()[0]
+
+ p_fem = P0(p_call)
+
+ a_sym_h = discretize(a_sym, domain_h, (V1_h, V1_h), backend=backend)
+ a_fem_h = discretize(a_fem, domain_h, (V1_h, V1_h), backend=backend)
+
+ A_sym = a_sym_h.assemble()
+ A_fem = a_fem_h.assemble(p=p_fem)
+
+ A_sym_sp = A_sym.tosparse()
+ A_fem_sp = A_fem.tosparse()
+ diff = A_sym_sp - A_fem_sp
+ norm_sym = np.linalg.norm(A_sym_sp.data)
+ norm_fem = np.linalg.norm(A_fem_sp.data)
+ error = np.linalg.norm(diff.data)
+ rel_err = error / norm_sym
+
+ #print(norm_sym, norm_fem, error, rel_err)
+
+ assert rel_err < 1e-4
diff --git a/feectools/cad/mesh/__init__.py b/feectools/cad/mesh/__init__.py
new file mode 100644
index 000000000..419109b64
--- /dev/null
+++ b/feectools/cad/mesh/__init__.py
@@ -0,0 +1,5 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
diff --git a/feectools/cad/mesh/generate_pipe.py b/feectools/cad/mesh/generate_pipe.py
new file mode 100644
index 000000000..5c8ba3e3b
--- /dev/null
+++ b/feectools/cad/mesh/generate_pipe.py
@@ -0,0 +1,25 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import numpy as np
+from igakit.cad import circle, ruled, bilinear, join
+
+from psydac.cad.geometry import export_nurbs_to_hdf5, refine_nurbs
+
+# create pipe geometry
+C0 = circle(center=(-1,0),angle=(-np.pi/3,0))
+C1 = circle(radius=2,center=(-1,0),angle=(-np.pi/3,0))
+annulus = ruled(C0,C1).transpose()
+square = bilinear(np.array([[[0,0],[0,3]],[[1,0],[1,3]]]) )
+pipe = join(annulus, square, axis=1)
+
+# refine the nurbs object
+ncells = [2**5,2**5]
+degree = [2,2]
+multiplicity = [2,2]
+
+new_pipe = refine_nurbs(pipe, ncells=ncells, degree=degree, multiplicity=multiplicity)
+filename = "pipe.h5"
+export_nurbs_to_hdf5(filename, new_pipe)
diff --git a/feectools/cad/mesh/multipatch/__init__.py b/feectools/cad/mesh/multipatch/__init__.py
new file mode 100644
index 000000000..419109b64
--- /dev/null
+++ b/feectools/cad/mesh/multipatch/__init__.py
@@ -0,0 +1,5 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
diff --git a/feectools/cmd/argparse_helpers.py b/feectools/cmd/argparse_helpers.py
new file mode 100644
index 000000000..9b0cf5626
--- /dev/null
+++ b/feectools/cmd/argparse_helpers.py
@@ -0,0 +1,71 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import sys
+import argparse
+
+from termcolor import colored
+
+from psydac import __version__ as psydac_version
+from psydac import __path__ as psydac_path
+
+__all__ = (
+ 'add_help_flag',
+ 'add_version_flag',
+ 'exit_with_error_message',
+)
+
+#------------------------------------------------------------------------------
+def add_help_flag(parser: argparse.ArgumentParser) -> None:
+ """
+ Add `-h/--help` flag to argument parser.
+
+ Add `-h/--help` flag to argument parser.
+
+ Parameters
+ ----------
+ parser : argparse.ArgumentParser
+ The parser to be modified.
+ """
+ message = 'Show this help message and exit.'
+ parser.add_argument('-h', '--help', action='help', help=message)
+
+#------------------------------------------------------------------------------
+def add_version_flag(parser: argparse.ArgumentParser) -> None:
+ """
+ Add `-V/--version` flag to argument parser.
+
+ Add `-V/--version` flag to argument parser.
+
+ Parameters
+ ----------
+ parser : argparse.ArgumentParser
+ The parser to be modified.
+ """
+ version = psydac_version
+ libpath = psydac_path[0]
+ python = f'python {sys.version_info.major}.{sys.version_info.minor}'
+ message = f'psydac {version} from {libpath} ({python})'
+
+ parser.add_argument('-V', '--version', action='version',
+ help='Show version and exit.', version=message)
+
+#------------------------------------------------------------------------------
+def exit_with_error_message(msg: str) -> None:
+ """
+ Print a colored error message and exit with status code 2.
+
+ Print a colored error message and exit with status code 2.
+
+ Parameters
+ ----------
+ msg : str
+ The error message to be printed.
+ """
+ err = colored('ERROR', color='magenta', attrs=['bold'])
+ sep = colored(': ', color='magenta')
+ msg = colored(msg, color='magenta')
+ print(f'{err}{sep}{msg}')
+ sys.exit(2)
diff --git a/feectools/cmd/main.py b/feectools/cmd/main.py
new file mode 100644
index 000000000..0b8f10892
--- /dev/null
+++ b/feectools/cmd/main.py
@@ -0,0 +1,58 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import argparse
+import sys
+
+from psydac.cmd.argparse_helpers import add_help_flag, add_version_flag
+from psydac.cmd.psydac_compile import setup_psydac_compile_parser, psydac_compile, PSYDAC_COMPILE_DESCR
+from psydac.cmd.psydac_mesh import setup_psydac_mesh_parser, psydac_mesh, PSYDAC_MESH_DESCR
+from psydac.cmd.psydac_test import setup_psydac_test_parser, psydac_test, PSYDAC_TEST_DESCR
+
+#==============================================================================
+def psydac_command() -> None:
+ """
+ Main entry point for the `psydac` command line interface.
+
+ Main entry point for the `psydac` command line interface.
+ Parses the command line arguments and calls the appropriate sub-command.
+ """
+ parser = argparse.ArgumentParser(
+ description = 'PSYDAC: Python Spline librarY for Differential equations with Automatic Code generation',
+ add_help = False,
+ )
+
+ group = parser.add_argument_group("Options")
+ add_help_flag(group)
+ add_version_flag(group)
+
+ sub_commands = {
+ 'compile': (setup_psydac_compile_parser, psydac_compile, PSYDAC_COMPILE_DESCR),
+ 'mesh' : (setup_psydac_mesh_parser , psydac_mesh , PSYDAC_MESH_DESCR ),
+ 'test' : (setup_psydac_test_parser , psydac_test , PSYDAC_TEST_DESCR ),
+ }
+
+ subparsers = parser.add_subparsers(
+ required=True, title='Subcommands', metavar='COMMAND')
+
+ for key, (parser_setup, exe_func, descr) in sub_commands.items():
+ sparser = subparsers.add_parser(
+ key,
+ help = descr.splitlines()[0],
+ description = f"PSYDAC's CLI: {descr}",
+ formatter_class = argparse.RawDescriptionHelpFormatter,
+ add_help = False,
+ )
+ parser_setup(sparser)
+ sparser.set_defaults(func = exe_func)
+
+ argv = sys.argv[1:]
+ if len(argv) == 0:
+ parser.print_help()
+ sys.exit(2)
+
+ kwargs = vars(parser.parse_args())
+ func = kwargs.pop('func')
+ func(**kwargs)
diff --git a/feectools/cmd/psydac_compile.py b/feectools/cmd/psydac_compile.py
new file mode 100644
index 000000000..1198250e9
--- /dev/null
+++ b/feectools/cmd/psydac_compile.py
@@ -0,0 +1,82 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+"""
+The purpose of this module is to pyccelize all PSYDAC kernels, in the case
+that these were modified after an editable installation of PSYDAC.
+"""
+from psydac.cmd.argparse_helpers import (
+ add_help_flag,
+ add_version_flag,
+ exit_with_error_message,
+)
+
+__all__ = (
+ 'setup_psydac_compile_parser',
+ 'psydac_compile',
+ 'PSYDAC_COMPILE_DESCR',
+)
+
+PSYDAC_COMPILE_DESCR = "Accelerate all computational kernels in PSYDAC using Pyccel (editable install only)."
+
+#==============================================================================
+def setup_psydac_compile_parser(parser):
+ """
+ Add the `psydac compile` arguments to the parser.
+
+ Add the `psydac compile` arguments to the parser for command line arguments.
+
+ Parameters
+ ----------
+ parser : argparse.ArgumentParser
+ The parser to be modified.
+ """
+ parser.add_argument('--language',
+ type = str,
+ choices = ('fortran', 'c'),
+ default = 'fortran',
+ help = 'Language used to pyccelize all the kernel files (default: fortran).'
+ )
+ add_help_flag(parser)
+ add_version_flag(parser)
+
+#==============================================================================
+def psydac_compile(*, language):
+ """
+ Accelerate all computational kernels in PSYDAC using Pyccel (editable install only).
+
+ Accelerate all computational kernels in PSYDAC using Pyccel (editable install only).
+
+ Parameters
+ ----------
+ language : str
+ Language used to pyccelize the kernels files.
+ """
+ from pathlib import Path
+ import shutil
+ import subprocess
+ import psydac
+
+ # Absolute path to the psydac directory
+ psydac_path = Path(psydac.__file__).parent
+
+ # Glob pattern of all kernel files
+ glob = '**/*_kernels.py'
+
+ # Command to be executed
+ cmd = (
+ shutil.which('pyccel'),
+ 'make',
+ '--language', language,
+ '--glob', (psydac_path / glob).as_posix(),
+ '--openmp',
+ )
+
+ print('Executing command:')
+ print(f' {" ".join(cmd)}\n')
+ result = subprocess.run(cmd, shell=False)
+
+ if result.returncode != 0:
+ exit_with_error_message('failed to compile PSYDAC kernels.')
diff --git a/feectools/cmd/psydac_mesh.py b/feectools/cmd/psydac_mesh.py
new file mode 100644
index 000000000..dcf80892b
--- /dev/null
+++ b/feectools/cmd/psydac_mesh.py
@@ -0,0 +1,147 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import sys
+
+from psydac.cmd.argparse_helpers import add_help_flag, add_version_flag, exit_with_error_message
+from psydac.mapping.discrete_gallery import available_mappings_2d, available_mappings_3d
+
+__all__ = (
+ 'setup_psydac_mesh_parser',
+ 'psydac_mesh',
+ 'PSYDAC_MESH_DESCR',
+)
+
+PSYDAC_MESH_DESCR = f"""Generate an HDF5 geometry file with a discrete mapping.
+
+Available 2D mappings: {', '.join(available_mappings_2d)}.
+Available 3D mappings: {', '.join(available_mappings_3d)}.
+"""
+#==============================================================================
+def setup_psydac_mesh_parser(parser):
+ """
+ Add the `psydac mesh` arguments to the parser.
+
+ Add the `psydac mesh` arguments to the parser for command line arguments.
+
+ Parameters
+ ----------
+ parser : argparse.ArgumentParser
+ The parser to be modified.
+ """
+ group = parser.add_argument_group('Mapping')
+ megrp = group.add_mutually_exclusive_group(required=True)
+ megrp.add_argument('--map-2d',
+ type = str,
+ metavar = 'MAP-NAME',
+ choices = available_mappings_2d,
+ dest = 'map_2d',
+ help = '2D analytical mapping to be interpolated by a spline.'
+ )
+ megrp.add_argument('--map-3d',
+ type = str,
+ metavar = 'MAP-NAME',
+ choices = available_mappings_3d,
+ dest = 'map_3d',
+ help = '3D analytical mapping to be interpolated by a spline.'
+ )
+
+ group = parser.add_argument_group('Discretization')
+ group.add_argument('-n',
+ required = True,
+ type = int,
+ nargs = '+',
+ dest = 'ncells',
+ metavar = ('N1','N2'),
+ help = 'Number of grid cells (elements) along each dimension.'
+ )
+ group.add_argument('-d',
+ required = True,
+ type = int,
+ nargs = '+',
+ dest = 'degree',
+ metavar = ('P1','P2'),
+ help = 'Spline degree along each dimension.'
+ )
+
+ group = parser.add_argument_group('Other options')
+ group.add_argument( '-o',
+ type = str,
+ default = 'out.h5',
+ dest = 'filename',
+ help = 'Name of output geometry file (default: out.h5).'
+ )
+ add_help_flag(group)
+ add_version_flag(group)
+
+#==============================================================================
+def psydac_mesh(*, map_2d, map_3d, ncells, degree, filename):
+ """
+ Generate an HDF5 geometry file with a discrete mapping.
+ """
+ if map_2d:
+ ndim = 2
+ error = len(degree) != 2 or len(ncells) != 2
+ elif map_3d:
+ ndim = 3
+ error = len(degree) != 3 or len(ncells) != 3
+
+ if error:
+ exit_with_error_message(
+ f'ncells and degree must have length {ndim} for a {ndim}D mapping'
+ )
+
+ map_name = map_2d or map_3d
+ export_analytical_mapping(map_name, ncells, degree, filename)
+
+#==============================================================================
+def export_analytical_mapping(mapping, ncells, degree, filename, **kwargs):
+
+ from psydac.cad.geometry import Geometry
+ from psydac.mapping.discrete_gallery import discrete_mapping
+
+ # create the discrete mapping from an analytical one
+ mapping = discrete_mapping(mapping, ncells=ncells, degree=degree)
+
+ # create a geometry from a discrete mapping
+ geometry = Geometry.from_discrete_mapping(mapping)
+
+ # export the geometry
+ geometry.export(filename)
+
+#==============================================================================
+# TODO [YG 13.01.2026] fix this using psydac.cad.gallery
+#def psydac_mesh(*, filename, geo_name, map_name, degree, ncells):
+#
+# if len(degree) != len(ncells):
+# raise ValueError('> ncells and degree must have same dimension')
+#
+# if geo_name:
+# export_caid_geometry(geo_name, ncells, degree, filename)
+#
+# elif map_name:
+# export_analytical_mapping(map_name, ncells, degree, filename)
+#
+#==============================================================================
+# TODO [YG 13.01.2026] fix this using psydac.cad.gallery
+#def export_caid_geometry(name, ncells, degree, filename):
+#
+# import os.path
+#
+# from caid.cad_geometry import cad_geometry
+# from caid.cad_geometry import line
+# from caid.cad_geometry import square
+# from caid.cad_geometry import circle
+# from caid.cad_geometry import cube
+#
+# constructor = eval(name)
+# geo = constructor(n=[i-1 for i in ncells], p=degree)
+# # ...
+#
+# extension = os.path.splitext(filename)[1]
+# if not extension == '.h5':
+# raise ValueError('> Only h5 extension is allowed for filename')
+#
+# geo.save(filename)
diff --git a/feectools/cmd/psydac_test.py b/feectools/cmd/psydac_test.py
new file mode 100644
index 000000000..8b6655274
--- /dev/null
+++ b/feectools/cmd/psydac_test.py
@@ -0,0 +1,178 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+"""
+The purpose of this module is to pyccelize all PSYDAC kernels, in the case
+that these were modified after an editable installation of PSYDAC.
+"""
+from psydac.cmd.argparse_helpers import (
+ add_help_flag,
+ add_version_flag,
+ exit_with_error_message,
+)
+
+__all__ = (
+ 'setup_psydac_test_parser',
+ 'psydac_test',
+ 'PSYDAC_TEST_DESCR',
+)
+
+PSYDAC_TEST_DESCR = "Run the PSYDAC test suite."
+
+#==============================================================================
+def setup_psydac_test_parser(parser):
+ """
+ Add the `psydac test` arguments to the parser.
+
+ Add the `psydac test` arguments to the parser for command line arguments.
+
+ Parameters
+ ----------
+ parser : argparse.ArgumentParser
+ The parser to be modified.
+ """
+
+ group = parser.add_argument_group('Test selection')
+ group.add_argument('--mod',
+ type = str,
+ help = 'Only run tests from the specified module (e.g. psydac.linalg).'
+ )
+ group.add_argument('--mpi',
+ action = 'store_true',
+ help = 'Only run parallel tests with 4 MPI processes (default: run serial tests).',
+ )
+ group.add_argument('--petsc',
+ action = 'store_true',
+ help = 'Only run tests using PETSc and petsc4py (default: run the other tests).',
+ )
+
+ group = parser.add_argument_group('Pytest options')
+ group.add_argument('-v', '--verbose',
+ action = 'store_true',
+ help = 'Increase verbosity of Pytest output.'
+ )
+ group.add_argument('-x', '--exitfirst',
+ action = 'store_true',
+ help = 'Exit instantly on first error or failed test.'
+ )
+
+ group = parser.add_argument_group('Other options')
+ add_help_flag(group)
+ add_version_flag(group)
+
+#==============================================================================
+def psydac_test(*, mod, mpi, petsc, verbose, exitfirst):
+ """
+ Run the PSYDAC test suite.
+ """
+ if mod is None:
+ mod = 'psydac'
+ else:
+ submods = mod.split('.')
+ if len(submods) == 0:
+ exit_with_error_message("module name cannot be empty")
+ elif submods[0] != 'psydac':
+ exit_with_error_message("module name must start with 'psydac'")
+ try:
+ modname = mod.split('::')[0]
+ import importlib
+ importlib.import_module(modname)
+ except ImportError:
+ exit_with_error_message(f"module '{modname}' not found")
+
+ # Import modules here to speed up parser
+ import os
+ import shutil
+ import subprocess
+ import time
+
+ # Clear Pytest cache from the current working directory
+ cache_dir = '.pytest_cache'
+ if os.path.isdir(cache_dir):
+ print(f'Removing existing Pytest cache directory: {cache_dir}\n', flush=True)
+ shutil.rmtree(cache_dir)
+
+ # If no pytest.toml file exists in the current working directory, copy it
+ # from the parent directory of this script (which is installed with PSYDAC)
+ if not os.path.isfile('pytest.toml'):
+ script_dir = os.path.dirname(os.path.abspath(__file__))
+ parent_dir = os.path.dirname(script_dir)
+ pytest_cfg = os.path.join(parent_dir, 'pytest.toml')
+ if not os.path.isfile(pytest_cfg):
+ exit_with_error_message(f'could not find pytest.toml file in {parent_dir}')
+ else:
+ print(f'Copying pytest.toml from: {parent_dir}\n', flush=True)
+ shutil.copy(pytest_cfg, os.getcwd())
+
+ # Build the list of flags for pytest
+ flags = []
+
+ # Set up MPI execution command, if needed
+ if mpi:
+ mpirun = shutil.which('mpirun')
+ mpi_exe = [mpirun, '-n', '4']
+ flags.append('--with-mpi')
+
+ # Determine if we are using OpenMPI, MPICH, or Intel MPI
+ result = subprocess.run([mpirun, '--version'], capture_output=True, text=True)
+ output = result.stdout.lower() + result.stderr.lower()
+ if 'open mpi' in output:
+ mpi_implementation = 'OpenMPI'
+ oversubscribe_flag = '--oversubscribe'
+ elif 'mpich' in output:
+ mpi_implementation = 'MPICH'
+ oversubscribe_flag = ''
+ elif 'intel mpi' in output:
+ mpi_implementation = 'Intel MPI'
+ oversubscribe_flag = '--oversubscribe' # to be verified
+ else:
+ exit_with_error_message("cannot determine MPI implementation from output of 'mpirun --version'")
+
+ print(f'MPI implementation detected: {mpi_implementation}')
+ if oversubscribe_flag:
+ mpi_exe.append(oversubscribe_flag)
+
+ else:
+ mpi_exe = []
+ flags.extend(['-n', 'auto', '--dist', 'loadgroup']) # for pytest-xdist
+
+ # If PETSc tests are requested, check that petsc4py is installed
+ if petsc:
+ try:
+ import petsc4py # noqa: F401
+ except ImportError:
+ exit_with_error_message("petsc4py is not installed, cannot run PETSc tests")
+
+ # Add appropriate markers for test selection
+ flags.extend(['--pyargs', mod])
+ mpi_mark = 'mpi' if mpi else 'not mpi'
+ petsc_mark = 'petsc' if petsc else 'not petsc'
+ flags.extend(['-m', f'({mpi_mark} and {petsc_mark})'])
+
+ # Default flags for pytest
+ flags.append('-ra') # show extra test summary info for skipped, failed, etc.
+
+ # Additional flags
+ if verbose:
+ flags.append('-v')
+ if exitfirst:
+ flags.append('-x')
+
+ # Command to be executed
+ cmd = [*mpi_exe, shutil.which('pytest'), *flags]
+
+ # Print command
+ cmd_to_print = [a.replace('(', '"(',).replace(')', ')"') for a in cmd]
+ print('Executing command:')
+ print(f' {" ".join(cmd_to_print)}', end='\n\n', flush=True)
+ time.sleep(0.1) # ensure the print is shown before subprocess output
+
+ # Execute the command
+ result = subprocess.run(cmd, shell=False, env=os.environ)
+
+ if result.returncode != 0:
+ msg = 'the PSYDAC test suite failed. '\
+ 'Please check the output above for details.'
+ exit_with_error_message(msg)
diff --git a/feectools/cmd/tests/__init__.py b/feectools/cmd/tests/__init__.py
new file mode 100644
index 000000000..419109b64
--- /dev/null
+++ b/feectools/cmd/tests/__init__.py
@@ -0,0 +1,5 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
diff --git a/feectools/cmd/tests/failing_test.py b/feectools/cmd/tests/failing_test.py
new file mode 100644
index 000000000..6c27828ab
--- /dev/null
+++ b/feectools/cmd/tests/failing_test.py
@@ -0,0 +1,26 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+"""
+This module contains only tests which are designed to fail, to check that the
+error handling in the `psydac test` command works correctly. This will be used
+in the CI to verify that the test suite correctly reports failures, by running
+`psydac test` on this file and verifying that the CI fails as expected.
+"""
+import pytest
+
+
+msg_tmp = "This {}test is designed to fail to check error handling in the test suite."
+
+def test_failure():
+ assert False, msg_tmp.format("")
+
+@pytest.mark.mpi
+def test_failure_mpi():
+ assert False, msg_tmp.format("MPI ")
+
+@pytest.mark.petsc
+def test_failure_petsc():
+ assert False, msg_tmp.format("PETSc ")
diff --git a/feectools/core/bsplines.py b/feectools/core/bsplines.py
index 900178871..369d2f895 100644
--- a/feectools/core/bsplines.py
+++ b/feectools/core/bsplines.py
@@ -1,7 +1,8 @@
-# coding: utf-8
-#
-# Copyright 2018 Yaman Güçlü
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
"""
Basic module that provides the means for evaluating the B-Splines basis
functions and their derivatives. In order to simplify automatic Fortran code
diff --git a/feectools/core/bsplines_kernels.py b/feectools/core/bsplines_kernels.py
index c8ba08e92..0ad57da3c 100644
--- a/feectools/core/bsplines_kernels.py
+++ b/feectools/core/bsplines_kernels.py
@@ -1,3 +1,9 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+
# This file holds the pyccelisable versions of the functions in bsplines.py
# This will be changed once pyccel can return arrays and can get out=None arguments
# like Numpy functions.
diff --git a/feectools/core/tests/test_bsplines.py b/feectools/core/tests/test_bsplines.py
index fe588c071..da6a213e3 100644
--- a/feectools/core/tests/test_bsplines.py
+++ b/feectools/core/tests/test_bsplines.py
@@ -1,5 +1,8 @@
-#coding: utf-8
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
import numpy as np
diff --git a/feectools/core/tests/test_bsplines_pyccel.py b/feectools/core/tests/test_bsplines_pyccel.py
index 768ca08e8..7ac56b3f7 100644
--- a/feectools/core/tests/test_bsplines_pyccel.py
+++ b/feectools/core/tests/test_bsplines_pyccel.py
@@ -1,3 +1,8 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import numpy as np
import pytest
diff --git a/feectools/ddm/TODO b/feectools/ddm/TODO
deleted file mode 100644
index c222a9629..000000000
--- a/feectools/ddm/TODO
+++ /dev/null
@@ -1,3 +0,0 @@
-TODO
-====
-* Update Poisson test
diff --git a/feectools/ddm/basic.py b/feectools/ddm/basic.py
index 5c217ad39..8d1d5c23a 100644
--- a/feectools/ddm/basic.py
+++ b/feectools/ddm/basic.py
@@ -1,7 +1,11 @@
-# coding: utf-8
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
from abc import ABC, abstractmethod
+
__all__ = ('CartDataExchanger',)
#===============================================================================
class CartDataExchanger(ABC):
diff --git a/feectools/ddm/blocking_data_exchanger.py b/feectools/ddm/blocking_data_exchanger.py
index 6a7e8adf7..8e3a46662 100644
--- a/feectools/ddm/blocking_data_exchanger.py
+++ b/feectools/ddm/blocking_data_exchanger.py
@@ -1,11 +1,15 @@
-# coding: utf-8
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import numpy as np
from feectools.ddm.mpi import mpi as MPI
from .cart import CartDecomposition, find_mpi_type
from .basic import CartDataExchanger
+
__all__ = ('BlockingCartDataExchanger',)
class BlockingCartDataExchanger(CartDataExchanger):
diff --git a/feectools/ddm/partition.py b/feectools/ddm/partition.py
index fb8519430..39b1981b5 100644
--- a/feectools/ddm/partition.py
+++ b/feectools/ddm/partition.py
@@ -1,8 +1,14 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import numpy as np
import numpy.ma as ma
from sympy.ntheory import factorint
+
__all__ = ('compute_dims', 'partition_procs_per_patch')
#==============================================================================
diff --git a/feectools/ddm/petsc.py b/feectools/ddm/petsc.py
index 88a3b8abe..538766413 100644
--- a/feectools/ddm/petsc.py
+++ b/feectools/ddm/petsc.py
@@ -1,7 +1,11 @@
-# coding: utf-8
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+from itertools import product
import numpy as np
-from itertools import product
from .cart import CartDecomposition
diff --git a/feectools/ddm/tests/test_cart_1d.py b/feectools/ddm/tests/test_cart_1d.py
index 6afc2c09b..1aa761c8e 100644
--- a/feectools/ddm/tests/test_cart_1d.py
+++ b/feectools/ddm/tests/test_cart_1d.py
@@ -1,5 +1,8 @@
-# Contents of test_cart_1d.py
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import numpy as np
from feectools.ddm.blocking_data_exchanger import BlockingCartDataExchanger
@@ -111,7 +114,7 @@ def run_cart_1d( data_exchanger_type, verbose=False ):
import pytest
@pytest.mark.parametrize( 'data_exchanger_type', [BlockingCartDataExchanger, NonBlockingCartDataExchanger] )
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_cart_1d( data_exchanger_type ):
namespace = run_cart_1d( data_exchanger_type )
diff --git a/feectools/ddm/tests/test_multicart_2d.py b/feectools/ddm/tests/test_multicart_2d.py
index a66dae0f3..3e6103694 100644
--- a/feectools/ddm/tests/test_multicart_2d.py
+++ b/feectools/ddm/tests/test_multicart_2d.py
@@ -1,4 +1,8 @@
-# File test_multicart_2d.py
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
#===============================================================================
# TEST MultiCartDecomposition in 2D
diff --git a/feectools/ddm/tests/test_partition.py b/feectools/ddm/tests/test_partition.py
index 1e20e5c82..808207d62 100644
--- a/feectools/ddm/tests/test_partition.py
+++ b/feectools/ddm/tests/test_partition.py
@@ -1,3 +1,8 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
from feectools.ddm.partition import compute_dims
@@ -119,6 +124,3 @@ def test_partition_3d_dims_mask( mpi_size, npts, mask ):
# test_partition_2d_dims_mask(10, [58, 64], [True, False])
test_partition_3d_dims_mask(10, [32, 64, 128], [True, False, True])
-
-
-
diff --git a/feectools/ddm/utilities.py b/feectools/ddm/utilities.py
index 35cfbc8b3..34b6dac36 100644
--- a/feectools/ddm/utilities.py
+++ b/feectools/ddm/utilities.py
@@ -1,10 +1,14 @@
-# coding: utf-8
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
from .cart import CartDecomposition, InterfaceCartDecomposition
from .blocking_data_exchanger import BlockingCartDataExchanger
from .nonblocking_data_exchanger import NonBlockingCartDataExchanger
from .interface_data_exchanger import InterfaceCartDataExchanger
+
__all__ = ('get_data_exchanger',)
def get_data_exchanger(cart, dtype, *, coeff_shape=(), assembly=False, axis=None, shape=None, blocking=True):
diff --git a/feectools/fem/grid.py b/feectools/fem/grid.py
index 58bfd1726..522839792 100644
--- a/feectools/fem/grid.py
+++ b/feectools/fem/grid.py
@@ -1,7 +1,8 @@
-# coding: utf-8
-#
-# Copyright 2018 Yaman Güçlü
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import numpy as np
from feectools.core.bsplines import elements_spans
diff --git a/feectools/fem/lst_preconditioner.py b/feectools/fem/lst_preconditioner.py
new file mode 100644
index 000000000..630ddc4ba
--- /dev/null
+++ b/feectools/fem/lst_preconditioner.py
@@ -0,0 +1,292 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+from functools import lru_cache
+from scipy.sparse import dia_matrix
+import numpy as np
+
+from sympde.topology import elements_of, Line, ScalarFunctionSpace
+from sympde.topology.datatype import SpaceType
+from sympde.expr import integral, BilinearForm
+
+from psydac.linalg.basic import IdentityOperator, LinearOperator
+from psydac.linalg.block import BlockVectorSpace, BlockLinearOperator
+from psydac.linalg.direct_solvers import BandedSolver
+from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix
+from psydac.linalg.stencil import StencilVectorSpace
+from psydac.fem.projectors import DirichletProjector
+
+
+@lru_cache
+def construct_LST_preconditioner(M, domain_h, fem_space, hom_bc=False, kind=None):
+ """
+ LST (Loli, Sangalli, Tani) preconditioners [1] are mass matrix preconditioners of the form
+ pc = D_inv_sqrt @ D_log_sqrt @ M_log_kron_solver @ D_log_sqrt @ D_inv_sqrt, where
+
+ D_inv_sqrt is the diagonal matrix of the square roots of the inverse diagonal entries of the mass matrix M,
+ D_log_sqrt is the diagonal matrix of the square roots of the diagonal entries of the mass matrix on the logical domain,
+ M_log_kron_solver is the Kronecker Solver of the mass matrix on the logical domain.
+
+ These preconditioners work very well even on complex domains as numerical experiments have shown.
+ Upon choosing hom_bc=True, a preconditioner for the modified mass matrix M_0 is returned.
+ M_0 is a mass matrix of the form
+ M_0 = DP @ M @ DP + (I - DP)
+ where DP and I are the corresponding DirichletProjector and IdentityOperator.
+ See examples/vector_potential_3d.
+
+ Parameters
+ ----------
+ M : psydac.linalg.stencil.StencilMatrix | psydac.linalg.block.BlockLinearOperator
+ Mass matrix corresponding to fem_space
+
+ domain_h : psydac.cad.geometry.Geometry
+ discretized physical domain used to discretize fem_space
+
+ fem_space : psydac.fem.basic.FemSpace
+ discretized Scalar- or VectorFunctionSpace. M is the corresponding mass matrix
+
+ hom_bc : bool
+ If True, return LST preconditioner for modified M_0 = DP @ M @ DP + (I - DP) mass matrix.
+ The argument M in that case remains the same (M, not M_0). DP and I are DirichletProjector and IdentityOperator.
+ Default: False.
+
+ kind : str | None
+ Optional. Must be passed if fem_space has no kind. Must match the kind of fem_space if fem_space has a kind.
+ Relevant as we must know whether M is a H1, Hcurl, Hdiv or L2 mass matrix.
+
+ Returns
+ -------
+ psydac.linalg.stencil.StencilMatrix | psydac.linalg.block.BlockLinearOperator
+ LST preconditioner for M (hom_bc=False) or M_0 (hom_b=True).
+
+ References
+ ----------
+ [1] Gabriele Loli, Giancarlo Sangalli, Mattia Tani. “Easy and efficient preconditioning of the isogeometric mass
+ matrix”. In: Computers & Mathematics with Applications 116 (2022), pp. 245–264
+ """
+
+ # to avoid circular import
+ from psydac.api.discretization import discretize
+
+ dim = fem_space.ldim
+ # In 1D one can solve the linear system directly (instead of using this preconditioner)
+ assert dim in (2, 3)
+
+ if hom_bc == True:
+ def toarray_1d(A):
+ """
+ Obtain a numpy array representation of a (1D) LinearOperator (which has not implemented toarray()).
+
+ We fill an empty numpy array row by row by repeatedly applying unit vectors
+ to the transpose of A. In order to obtain those unit vectors in Stencil format,
+ we make use of an auxiliary function that takes periodicity into account.
+ """
+
+ assert isinstance(A, LinearOperator)
+ W = A.codomain
+ assert isinstance(W, StencilVectorSpace)
+
+ def get_unit_vector_1d(v, periodic, n1, npts1, pads1):
+
+ v *= 0.0
+ v._data[pads1+n1] = 1.
+
+ if periodic:
+ if n1 < pads1:
+ v._data[-pads1+n1] = 1.
+ if n1 >= npts1-pads1:
+ v._data[n1-npts1+pads1] = 1.
+
+ return v
+
+ periods = W.periods
+ periodic = periods[0]
+
+ w = W.zeros()
+ At = A.T
+
+ A_arr = np.zeros(A.shape, dtype=A.dtype)
+
+ npts1, = W.npts
+ pads1, = W.pads
+ for n1 in range(npts1):
+ e_n1 = get_unit_vector_1d(w, periodic, n1, npts1, pads1)
+ A_n1 = At @ e_n1
+ A_arr[n1, :] = A_n1.toarray()
+
+ return A_arr
+
+ def M_0_1d_to_bandsolver(A):
+ """
+ Converts the M0_0_1d StencilMatrix to a BandedSolver.
+
+ Closely resembles BandedSolver.from_stencil_mat_1d,
+ the difference being that M0_0_1d neither has a
+ remove_spurious_entries() nor a toarray() function.
+
+ """
+
+ dmat = dia_matrix(toarray_1d(A), dtype=A.dtype)
+ la = abs(dmat.offsets.min())
+ ua = dmat.offsets.max()
+ cmat = dmat.tocsr()
+
+ A_bnd = np.zeros((1+ua+2*la, cmat.shape[1]), A.dtype)
+
+ for i,j in zip(*cmat.nonzero()):
+ A_bnd[la+ua+i-j, j] = cmat[i,j]
+
+ return BandedSolver(ua, la, A_bnd)
+
+ domain = domain_h.domain
+
+ ncells, = domain_h.ncells.values()
+ degree = fem_space.degree
+ periodic, = domain_h.periodic.values()
+
+ V_cs = fem_space.coeff_space
+
+ logical_domain = domain.logical_domain
+
+ # ----- Compute D_inv_sqrt
+
+ D_inv_sqrt = M.diagonal(inverse=True, sqrt=True)
+
+ # ----- Compute M_log_kron_solver
+
+ logical_domain_1d_x = Line('L', bounds=logical_domain.bounds1)
+ logical_domain_1d_y = Line('L', bounds=logical_domain.bounds2)
+ if dim == 3:
+ logical_domain_1d_z = Line('L', bounds=logical_domain.bounds3)
+
+ logical_domain_1d_list = [logical_domain_1d_x, logical_domain_1d_y]
+ if dim == 3:
+ logical_domain_1d_list += [logical_domain_1d_z]
+
+ # We gather the 1D mass matrices.
+ # Those will be used to obtain D_log_sqrt using the new
+ # diagonal function for KroneckerStencilMatrices.
+ M_1d_solvers = [[],[]]
+ Ms_1d = [[],[]]
+ if dim == 3:
+ M_1d_solvers += [[]]
+ Ms_1d += [[]]
+
+ # Mark 1D 'h1' spaces built using B-splines.
+ # 1D spaces for which (i, j) \notin keys are 'l2' spaces built using M-splines.
+ fem_space_kind = fem_space.symbolic_space.kind.name
+
+ if kind is not None:
+ if isinstance(kind, str):
+ kind = kind.lower()
+ assert(kind in ['h1', 'hcurl', 'hdiv', 'l2'])
+ elif isinstance(kind, SpaceType):
+ kind = kind.name
+ else:
+ raise TypeError(f'Expecting kind {kind} to be a str or of SpaceType')
+
+ # If fem_space has a kind, it must be compatible with kind
+ if fem_space_kind != 'undefined':
+ assert fem_space_kind == kind, f'fem_space and space_kind are not compatible.'
+ else:
+ kind = fem_space_kind
+
+ if kind == 'h1':
+ keys = ((0, 0), (1, 0), (2, 0))
+ elif kind == 'hcurl':
+ keys = ((0, 1), (0, 2), (1, 0), (1, 2), (2, 0), (2, 1))
+ elif kind == 'hdiv':
+ keys = ((0, 0), (1, 1), (2, 2))
+ elif kind == 'l2':
+ keys = ()
+ else:
+ raise ValueError(f'kind {kind} must be either h1, hcurl, hdiv or l2.')
+
+ for i, (ncells_1d, periodic_1d, logical_domain_1d) in enumerate(zip(ncells, periodic, logical_domain_1d_list)):
+
+ logical_domain_1d_h = discretize(logical_domain_1d, ncells=[ncells_1d, ], periodic=[periodic_1d, ])
+
+ degrees_1d = [degree_dir[i] for degree_dir in degree] if isinstance(fem_space.coeff_space, BlockVectorSpace) else [degree[i], ]
+
+ for j, d in enumerate(degrees_1d):
+
+ kind_1d = 'h1' if (i, j) in keys else 'l2'
+ basis = 'B' if (i, j) in keys else 'M'
+
+ if basis == 'M':
+ d += 1
+
+ V_1d = ScalarFunctionSpace('V', logical_domain_1d, kind =kind_1d)
+ Vh_1d = discretize(V_1d, logical_domain_1d_h, degree=[d,], basis=basis)
+
+ u, v = elements_of(V_1d, names='u, v')
+ a_1d = BilinearForm((u, v), integral(logical_domain_1d, u*v))
+ ah_1d = discretize(a_1d, logical_domain_1d_h, (Vh_1d, Vh_1d))
+ M_1d = ah_1d.assemble()
+ Ms_1d[j].append(M_1d)
+
+ if (hom_bc == True) and ((i, j) in keys):
+ DP = DirichletProjector(Vh_1d, space_kind='h1')
+ if DP.bcs != ():
+ I = IdentityOperator(Vh_1d.coeff_space)
+ M_0_1d = DP @ M_1d @ DP + (I - DP)
+
+ M_0_1d_solver = M_0_1d_to_bandsolver(M_0_1d)
+ M_1d_solvers[j].append(M_0_1d_solver)
+ else:
+ M_1d_solver = BandedSolver.from_stencil_mat_1d(M_1d)
+ M_1d_solvers[j].append(M_1d_solver)
+ else:
+ M_1d_solver = BandedSolver.from_stencil_mat_1d(M_1d)
+ M_1d_solvers[j].append(M_1d_solver)
+
+ if isinstance(V_cs, StencilVectorSpace):
+ M_log_kron_solver = KroneckerLinearSolver(V_cs, V_cs, M_1d_solvers[0])
+
+ else:
+ M_0_log_kron_solver = KroneckerLinearSolver(V_cs[0], V_cs[0], M_1d_solvers[0])
+ M_1_log_kron_solver = KroneckerLinearSolver(V_cs[1], V_cs[1], M_1d_solvers[1])
+ if dim == 3:
+ M_2_log_kron_solver = KroneckerLinearSolver(V_cs[2], V_cs[2], M_1d_solvers[2])
+
+ if dim == 2:
+ blocks = [[M_0_log_kron_solver, None],
+ [None, M_1_log_kron_solver]]
+ else:
+ blocks = [[M_0_log_kron_solver, None, None],
+ [None, M_1_log_kron_solver, None],
+ [None, None, M_2_log_kron_solver]]
+
+ M_log_kron_solver = BlockLinearOperator (V_cs, V_cs, blocks)
+
+ # ----- Compute D_log_sqrt
+
+ if isinstance(V_cs, StencilVectorSpace):
+ M_log = KroneckerStencilMatrix(V_cs, V_cs, *Ms_1d[0])
+
+ D_log_sqrt = M_log.diagonal (inverse=False, sqrt=True)
+
+ else:
+ M_0_log = KroneckerStencilMatrix(V_cs[0], V_cs[0], *Ms_1d[0])
+ M_1_log = KroneckerStencilMatrix(V_cs[1], V_cs[1], *Ms_1d[1])
+ if dim == 3:
+ M_2_log = KroneckerStencilMatrix(V_cs[2], V_cs[2], *Ms_1d[2])
+
+ if dim == 2:
+ blocks = [[M_0_log, None],
+ [None, M_1_log]]
+ else:
+ blocks = [[M_0_log, None, None],
+ [None, M_1_log, None],
+ [None, None, M_2_log]]
+
+ M_log = BlockLinearOperator(V_cs, V_cs, blocks=blocks)
+ D_log_sqrt = M_log.diagonal(inverse=False, sqrt=True)
+
+ # --------------------------------
+
+ M_pc = D_inv_sqrt @ D_log_sqrt @ M_log_kron_solver @ D_log_sqrt @ D_inv_sqrt
+
+ return M_pc
diff --git a/feectools/fem/tests/analytical_profiles_1d.py b/feectools/fem/tests/analytical_profiles_1d.py
index 7bdfe2991..fb16f77b0 100644
--- a/feectools/fem/tests/analytical_profiles_1d.py
+++ b/feectools/fem/tests/analytical_profiles_1d.py
@@ -1,13 +1,15 @@
-# coding: utf-8
-# Copyright 2018 Yaman Güçlü
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import math
import numpy as np
from feectools.fem.tests.analytical_profiles_base import AnalyticalProfile
from feectools.fem.tests.utilities import horner, falling_factorial
-__all__ = ('AnalyticalProfile_1d_cos', 'AnalyticalProfile1D_Poly')
+__all__ = ('AnalyticalProfile1D_Cos', 'AnalyticalProfile1D_Poly')
#===============================================================================
class AnalyticalProfile1D_Cos( AnalyticalProfile ):
diff --git a/feectools/fem/tests/analytical_profiles_base.py b/feectools/fem/tests/analytical_profiles_base.py
index 4d0ba4a0d..ae1ecf170 100644
--- a/feectools/fem/tests/analytical_profiles_base.py
+++ b/feectools/fem/tests/analytical_profiles_base.py
@@ -1,6 +1,8 @@
-# coding: utf-8
-# Copyright 2018 Yaman Güçlü
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
from abc import ABCMeta, abstractmethod
__all__ = ('AnalyticalProfile',)
diff --git a/feectools/fem/tests/splines_error_bounds.py b/feectools/fem/tests/splines_error_bounds.py
index 46b641a3b..49f14b0da 100644
--- a/feectools/fem/tests/splines_error_bounds.py
+++ b/feectools/fem/tests/splines_error_bounds.py
@@ -1,6 +1,9 @@
-# coding: utf-8
-# Copyright 2018 Yaman Güçlü
-#
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+
# This file is the Python translation of a Selalib Fortran module:
# 'selalib/src/splines/tests/m_splines_error_bounds.F90'
diff --git a/feectools/fem/tests/test_dirichlet_projectors.py b/feectools/fem/tests/test_dirichlet_projectors.py
new file mode 100644
index 000000000..02def9574
--- /dev/null
+++ b/feectools/fem/tests/test_dirichlet_projectors.py
@@ -0,0 +1,517 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import numpy as np
+import pytest
+
+# Skip entire module if sympde is not available
+try:
+ import sympde
+except ImportError:
+ pytest.skip("sympde not installed", allow_module_level=True)
+
+from mpi4py import MPI
+from sympy import sin, pi, sqrt, Tuple
+
+from sympde.calculus import inner, cross
+from sympde.expr import integral, LinearForm, BilinearForm
+from sympde.topology import elements_of, Derham, Mapping, Line, Square, Cube, Union, NormalVector, ScalarFunctionSpace, VectorFunctionSpace
+from sympde.topology.datatype import H1Space, HcurlSpace
+
+from psydac.api.discretization import discretize
+from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL
+from psydac.fem.projectors import DirichletProjector
+from psydac.linalg.basic import LinearOperator, IdentityOperator
+from psydac.linalg.block import BlockVectorSpace
+from psydac.linalg.solvers import inverse
+from psydac.linalg.tests.utilities import check_linop_equality_using_rng, SinMapping1D, Annulus, SquareTorus
+
+#===============================================================================
+@pytest.mark.parametrize('dim', [1, 2])
+
+def test_function_space_boundary_projector(dim):
+
+ tol = 1e-15
+
+ ncells_3d = [8, 8, 8]
+ degree_3d = [2, 2, 2]
+ periodic_3d = [False, True, False]
+
+ comm = None
+ backend = PSYDAC_BACKEND_GPYCCEL
+
+ logical_domain_1d = Line ('L', bounds= (0, 1))
+ logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi))
+ logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1))
+ logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d]
+
+ mapping_1d = SinMapping1D('LM')
+ mapping_2d = Annulus ('A' )
+ mapping_3d = SquareTorus ('ST')
+ mappings = [mapping_1d, mapping_2d, mapping_3d]
+
+ rng = np.random.default_rng(42)
+
+ print()
+ print(f' ----- Test projectors in dimension {dim} -----')
+ print()
+
+ domain = mappings[dim-1](logical_domains[dim-1])
+ from sympde.utilities.utils import plot_domain
+ #plot_domain(domain, draw=True, isolines=True)
+
+ # Obtain "true" boundary, i.e., remove periodic y-direction boundary
+ if dim == 1:
+ boundary = domain.boundary
+ elif dim == 2:
+ boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1))
+ else:
+ boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1),
+ domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1))
+
+ ncells = [ncells_3d[0], ] if dim == 1 else ncells_3d [0:dim]
+ degree = [degree_3d[0], ] if dim == 1 else degree_3d [0:dim]
+ periodic = [periodic_3d[0], ] if dim == 1 else periodic_3d[0:dim]
+
+ domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm)
+
+ nn = NormalVector('nn')
+
+ for i in range(dim):
+ print(f' - Test DP{i}')
+
+ # The function defined here satisfy the corresponding homogeneous Dirichlet BCs
+ if dim == 1:
+ x = domain.coordinates
+ V = ScalarFunctionSpace('V', domain, kind='H1') # testing various kind arguments
+ f = sin(2*pi*x)
+ if dim == 2:
+ x, y = domain.coordinates
+ if i == 0:
+ V = ScalarFunctionSpace('V', domain, kind=H1Space) # testing various kind arguments
+ f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1)
+ else:
+ V = VectorFunctionSpace('V', domain, kind='hCuRl') # testing various kind arguments
+ f1 = x
+ f2 = y
+ f = Tuple(f1, f2)
+ if dim == 3:
+ x, y, z = domain.coordinates
+ if i == 0:
+ V = ScalarFunctionSpace('V', domain, kind='h1') # testing various kind arguments
+ f = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1) * z * (z-1)
+ elif i == 1:
+ V = VectorFunctionSpace('V', domain, kind=HcurlSpace) # testing various kind arguments
+ f1 = z * (z - 1) * x
+ f2 = z * (z - 1) * y
+ f3 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1)
+ f = Tuple(f1, f2, f3)
+ else:
+ V = VectorFunctionSpace('V', domain, kind='Hdiv') # testing various kind arguments
+ f1 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1)
+ f2 = (sqrt(x**2 + y**2)-0.5) * (sqrt(x**2 + y**2)-1)
+ f3 = z * (z-1) * sin(x*y)
+ f = Tuple(f1, f2, f3)
+
+ u, v = elements_of(V, names='u, v')
+ if i == 0:
+ boundary_expr = u*v
+ if (i == 1) and (dim == 2):
+ boundary_expr = cross(nn, u) * cross(nn, v)
+ if (i == 1) and (dim == 3):
+ boundary_expr = inner(cross(nn, u), cross(nn, v))
+ if i == 2:
+ boundary_expr = inner(nn, u) * inner(nn, v)
+
+ Vh = discretize(V, domain_h, degree=degree)
+ expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v
+
+ a = BilinearForm((u, v), integral(domain, expr))
+ ab = BilinearForm((u, v), integral(boundary, boundary_expr))
+
+ ah = discretize(a, domain_h, (Vh, Vh), backend=backend)
+ abh = discretize(ab, domain_h, (Vh, Vh), backend=backend, sum_factorization=False)
+
+ I = IdentityOperator(Vh.coeff_space)
+ DP = DirichletProjector(Vh)
+
+ M = ah.assemble()
+ M_0 = DP @ M @ DP + (I - DP)
+ Mb = abh.assemble()
+
+ # We project f into the conforming discrete space using a penalization method. It's coefficients are stored in fc
+ lexpr = inner(v, f) if isinstance(Vh.coeff_space, BlockVectorSpace) else v*f
+ l = LinearForm(v, integral(domain, lexpr))
+ lh = discretize(l, domain_h, Vh, backend=backend)
+ rhs = lh.assemble()
+ A = M + 1e30*Mb
+ A_inv = inverse(A, 'cg', maxiter=1000, tol=1e-10)
+ fc = A_inv @ rhs
+
+ # 1.
+ # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet
+ # boundary conditions should not change under application of the corresponding projector
+ fc2 = DP @ fc
+ diff = fc - fc2
+ err_sqr = diff.inner(diff)
+ print(f' || f - P @ f ||^2 = {err_sqr}')
+ assert err_sqr < tol**2
+
+ # 2.1
+ # After applying a projector to a random vector, we want to verify that the
+ # corresponding boundary integral vanishes
+ rdm_coeffs = Vh.coeff_space.zeros()
+ print(' Random boundary integrals:')
+ for _ in range(3):
+ if isinstance(rdm_coeffs.space, BlockVectorSpace):
+ for block in rdm_coeffs.blocks:
+ rng.random(size=block._data.shape, dtype="float64", out=block._data)
+ else:
+ rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data)
+ rdm_coeffs2 = DP @ rdm_coeffs
+ scaled_boundary_int_rdm_sqr = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2
+ scaled_boundary_int_proj_rdm_sqr = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2
+ print(f' rdm: {scaled_boundary_int_rdm_sqr} proj. rdm: {scaled_boundary_int_proj_rdm_sqr}')
+ assert scaled_boundary_int_proj_rdm_sqr < tol**2
+
+ # 2.2
+ # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray())
+ DP_arr = DP.toarray()
+ rdm_coeffs_arr = rdm_coeffs.toarray()
+ diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray()
+ err_sqr = diff_arr.dot(diff_arr)
+ assert err_sqr < tol**2
+
+ # 3.
+ # We want to verify that applying a projector twice does not change the vector twice
+ fc3 = DP @ fc2
+ diff = fc2 - fc3
+ err_sqr = diff.inner(diff)
+ print(f' || P @ f - P @ P @ f ||^2 = {err_sqr}')
+ assert err_sqr < tol**2
+
+ # 4.
+ # Finally, the modified mass matrix should still compute inner products correctly
+ l2_norm_sqr = M.dot_inner (fc, fc)
+ l2_norm2_sqr = M_0.dot_inner(fc, fc)
+ err_sqr = abs(l2_norm_sqr - l2_norm2_sqr)
+ print(f' || P @ f ||^2 = {l2_norm_sqr} should be equal to')
+ print(f' || P @ f ||^2 (alt) = {l2_norm2_sqr}')
+ # M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same only up to order 1e-15.
+ # Hence, we can't expect err_sqr to be less than tol**2, but only less than tol.
+ assert err_sqr < tol
+
+ print()
+
+#===============================================================================
+@pytest.mark.parametrize('dim', [1, 3])
+@pytest.mark.mpi
+
+def test_discrete_derham_boundary_projector(dim):
+
+ tol = 1e-15
+
+ ncells = [8, 8, 8]
+ degree = [2, 2, 2]
+ periodic = [False, True, False]
+
+ comm = MPI.COMM_WORLD
+ backend = PSYDAC_BACKEND_GPYCCEL
+
+ logical_domain_1d = Line ('L', bounds= (0, 1))
+ logical_domain_2d = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi))
+ logical_domain_3d = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1))
+ logical_domains = [logical_domain_1d, logical_domain_2d, logical_domain_3d]
+
+ mapping_1d = SinMapping1D('LM')
+ mapping_2d = Annulus ('A' )
+ mapping_3d = SquareTorus ('ST')
+ mappings = [mapping_1d, mapping_2d, mapping_3d]
+
+ rng = np.random.default_rng(42)
+
+ # The following are functions (1D, 2D & 3D) satisfying homogeneous Dirichlet BCs
+
+ f11 = lambda x : np.sin(2*np.pi*x)
+
+ r2 = lambda x, y : np.sqrt(x**2 + y**2)
+ f21 = lambda x, y : (r2(x, y) - 0.5) * (r2(x, y) - 1)
+ f22_1 = lambda x, y : x
+ f22_2 = lambda x, y : y
+ f22 = (f22_1, f22_2)
+
+ f31 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1) * z * (z - 1)
+ f32_1 = lambda x, y, z : z * (z - 1) * x
+ f32_2 = lambda x, y, z : z * (z - 1) * y
+ f32_3 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1)
+ f32 = (f32_1, f32_2, f32_3)
+ f33_1 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1)
+ f33_2 = lambda x, y, z : (r2(x, y) - 0.5) * (r2(x, y) - 1)
+ f33_3 = lambda x, y, z : z * (z - 1) * np.sin(x*y)
+ f33 = (f33_1, f33_2, f33_3)
+
+ funs = [[f11], [f21, f22], [f31, f32, f33]]
+
+ print()
+ print(f' ----- Test projectors in dimension {dim} -----')
+ print()
+
+ domain = mappings[dim-1](logical_domains[dim-1])
+ from sympde.utilities.utils import plot_domain
+ #plot_domain(domain, draw=True, isolines=True)
+
+ # Obtain "true" boundary, i.e., remove periodic y-direction boundary
+ if dim == 1:
+ boundary = domain.boundary
+ elif dim == 2:
+ boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1))
+ else:
+ boundary = Union(domain.get_boundary(axis=0, ext=-1), domain.get_boundary(axis=0, ext=1),
+ domain.get_boundary(axis=2, ext=-1), domain.get_boundary(axis=2, ext=1))
+
+ derham = Derham(domain) if dim in (1, 3) else Derham(domain, sequence=['h1', 'hcurl', 'l2'])
+
+ ncells_dim = [ncells[0], ] if dim == 1 else ncells[0:dim]
+ degree_dim = [degree[0], ] if dim == 1 else degree[0:dim]
+ periodic_dim = [periodic[0], ] if dim == 1 else periodic[0:dim]
+
+ domain_h = discretize(domain, ncells=ncells_dim, periodic=periodic_dim, comm=comm)
+ derham_h = discretize(derham, domain_h, degree=degree_dim)
+
+ d_projectors = derham_h.dirichlet_projectors(kind='linop')
+
+ if dim == 2:
+ conf_projectors = derham_h.conforming_projectors(kind='linop', hom_bc=True)
+
+ nn = NormalVector('nn')
+
+ for i in range(dim):
+ print(f' - Test DP{i}')
+
+ u, v = elements_of(derham.spaces[i], names='u, v')
+
+ if i == 0:
+ boundary_expr = u*v
+ if (i == 1) and (dim == 2):
+ boundary_expr = cross(nn, u) * cross(nn, v)
+ if (i == 1) and (dim == 3):
+ boundary_expr = inner(cross(nn, u), cross(nn, v))
+ if i == 2:
+ boundary_expr = inner(nn, u) * inner(nn, v)
+
+ expr = inner(u, v) if isinstance(derham_h.spaces[i].coeff_space, BlockVectorSpace) else u*v
+
+ a = BilinearForm((u, v), integral(domain, expr))
+ ab = BilinearForm((u, v), integral(boundary, boundary_expr))
+
+ ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend)
+ abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False)
+
+ I = IdentityOperator(derham_h.spaces[i].coeff_space)
+ DP = d_projectors[i]
+
+ if dim == 2:
+ CP = conf_projectors[i]
+ check_linop_equality_using_rng(DP, CP)
+
+ M = ah.assemble()
+ M_0 = DP @ M @ DP + (I - DP)
+ Mb = abh.assemble()
+
+ f = funs[dim-1][i]
+ fc = derham_h.projectors()[i](f).coeffs
+
+ # 1.
+ # In 1D, 2D, 3D, the coefficients of functions satisfying homogeneous Dirichlet
+ # boundary conditions should not change under application of the corresponding projector
+ fc2 = DP @ fc
+ diff = fc - fc2
+ err_sqr = diff.inner(diff)
+ print(f' || f - P @ f ||^2 = {err_sqr}')
+ assert err_sqr < tol**2
+
+ # 2.1
+ # After applying a projector to a random vector, we want to verify that the
+ # corresponding boundary integral vanishes
+ rdm_coeffs = derham_h.spaces[i].coeff_space.zeros()
+ print(' Random boundary integrals:')
+ for _ in range(3):
+ if isinstance(rdm_coeffs.space, BlockVectorSpace):
+ for block in rdm_coeffs.blocks:
+ rng.random(size=block._data.shape, dtype="float64", out=block._data)
+ else:
+ rng.random(size=rdm_coeffs._data.shape, dtype="float64", out=rdm_coeffs._data)
+ rdm_coeffs2 = DP @ rdm_coeffs
+ scaled_boundary_int_rdm_sqr = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2
+ scaled_boundary_int_proj_rdm_sqr = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2
+ print(f' rdm: {scaled_boundary_int_rdm_sqr} proj. rdm: {scaled_boundary_int_proj_rdm_sqr}')
+ assert scaled_boundary_int_proj_rdm_sqr < tol**2
+
+ # 2.2
+ # Test tosparse(): (DP @ rdm_coeffs).toarray() should be equal to DP.tosparse().dot(rdm_coeffs.toarray())
+ DP_spr = DP.tosparse()
+ rdm_coeffs_arr = rdm_coeffs.toarray()
+ diff_arr = DP_spr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray()
+ err_sqr = diff_arr.dot(diff_arr)
+ assert err_sqr < tol**2
+
+
+ # 3.
+ # We want to verify that applying a projector twice does not change the vector twice
+ fc3 = DP @ fc2
+ diff = fc2 - fc3
+ err_sqr = diff.inner(diff)
+ print(f' || P @ f - P @ P @ f ||^2 = {err_sqr}')
+ assert err_sqr < tol**2
+
+ # 4.
+ # Finally, the modified mass matrix should still compute inner products correctly
+ l2_norm_sqr = M.dot_inner (fc, fc)
+ l2_norm2_sqr = M_0.dot_inner(fc, fc)
+ err_sqr = abs(l2_norm_sqr - l2_norm2_sqr)
+ print(f' || P @ f ||^2 = {l2_norm_sqr} should be equal to')
+ print(f' || P @ f ||^2 (alt) = {l2_norm2_sqr}')
+ # M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same only up to order 1e-15.
+ # Hence, we can't expect err_sqr to be less than tol**2, but only less than tol.
+ assert err_sqr < tol
+
+ print()
+
+#===============================================================================
+def test_discrete_derham_boundary_projector_multipatch():
+
+ tol = 1e-15
+
+ ncells = [8, 8]
+ degree = [2, 2]
+
+ comm = None
+ backend = PSYDAC_BACKEND_GPYCCEL
+
+ from psydac.feec.multipatch_domain_utilities import build_multipatch_domain
+ domain = build_multipatch_domain(domain_name='annulus_3')
+
+ rng = np.random.default_rng(42)
+
+ # The following are functions satisfying homogeneous Dirichlet BCs
+ r = lambda x, y : np.sqrt(x**2 + y**2)
+ f1 = lambda x, y : (r(x, y) - 0.5) * (r(x, y) - 1)
+ f2_1 = lambda x, y : x
+ f2_2 = lambda x, y : y
+ f2 = (f2_1, f2_2)
+ funs = [f1, f2]
+ print()
+
+ boundary = domain.boundary
+
+ derham = Derham(domain, sequence=['h1', 'hcurl', 'l2'])
+
+ ncells_h = {}
+ for D in domain.interior:
+ ncells_h[D.name] = ncells
+
+ domain_h = discretize(domain, ncells=ncells_h, comm=comm)
+ derham_h = discretize(derham, domain_h, degree=degree)
+
+ projectors = derham_h.projectors(nquads=[(d + 1) for d in degree])
+
+ d_projectors = derham_h.dirichlet_projectors(kind='linop')
+
+ nn = NormalVector('nn')
+
+ for i in range(2):
+ print(f' - Test DP{i}')
+
+ u, v = elements_of(derham.spaces[i], names='u, v')
+
+ if i == 0:
+ boundary_expr = u*v
+ expr = u*v
+ if (i == 1):
+ boundary_expr = cross(nn, u) * cross(nn, v)
+ expr = inner(u,v)
+
+ a = BilinearForm((u, v), integral(domain, expr))
+ ab = BilinearForm((u, v), integral(boundary, boundary_expr))
+
+ ah = discretize(a, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend)
+ abh = discretize(ab, domain_h, (derham_h.spaces[i], derham_h.spaces[i]), backend=backend, sum_factorization=False)
+
+ I = IdentityOperator(derham_h.spaces[i].coeff_space)
+ DP = d_projectors[i]
+
+ M = ah.assemble()
+ M_0 = DP @ M @ DP + (I - DP)
+ Mb = abh.assemble()
+
+ f = funs[i]
+ fc = projectors[i](f).coeffs
+
+ # 1.
+ # The coefficients of functions satisfying homogeneous Dirichlet
+ # boundary conditions should not change under application of the corresponding projector
+ fc2 = DP @ fc
+ diff = fc - fc2
+ err_sqr = diff.inner(diff)
+ print(f' || f - P @ f ||^2 = {err_sqr}')
+ assert err_sqr < tol**2
+
+ # 2.1
+ # After applying a projector to a random vector, we want to verify that the
+ # corresponding boundary integral vanishes
+ rdm_coeffs = derham_h.spaces[i].coeff_space.zeros()
+ print(' Random boundary integrals:')
+ for _ in range(3):
+ for patch in rdm_coeffs.blocks:
+
+ if isinstance(patch.space, BlockVectorSpace):
+ for block in patch.blocks:
+ rng.random(size=block._data.shape, dtype="float64", out=block._data)
+ else:
+ rng.random(size=patch._data.shape, dtype="float64", out=patch._data)
+
+ rdm_coeffs2 = DP @ rdm_coeffs
+ scaled_boundary_int_rdm_sqr = Mb.dot_inner(rdm_coeffs, rdm_coeffs) / rdm_coeffs.space.dimension**2
+ scaled_boundary_int_proj_rdm_sqr = Mb.dot_inner(rdm_coeffs2, rdm_coeffs2) / rdm_coeffs.space.dimension**2
+ print(f' rdm: {scaled_boundary_int_rdm_sqr} proj. rdm: {scaled_boundary_int_proj_rdm_sqr}')
+ assert scaled_boundary_int_proj_rdm_sqr < tol**2
+
+ # 2.2
+ # Test toarray(): (DP @ rdm_coeffs).toarray() should be equal to DP.toarray().dot(rdm_coeffs.toarray())
+ DP_arr = DP.toarray()
+ rdm_coeffs_arr = rdm_coeffs.toarray()
+ diff_arr = DP_arr.dot(rdm_coeffs_arr) - rdm_coeffs2.toarray()
+ err_sqr = diff_arr.dot(diff_arr)
+ assert err_sqr < tol**2
+
+ # 3.
+ # We want to verify that applying a projector twice does not change the vector twice
+ fc3 = DP @ fc2
+ diff = fc2 - fc3
+ err_sqr = diff.inner(diff)
+ print(f' || P @ f - P @ P @ f ||^2 = {err_sqr}')
+ assert err_sqr < tol**2
+
+ # 4.
+ # Finally, the modified mass matrix should still compute inner products correctly
+ l2_norm_sqr = M.dot_inner (fc, fc)
+ l2_norm2_sqr = M_0.dot_inner(fc, fc)
+ err_sqr = abs(l2_norm_sqr - l2_norm2_sqr)
+ print(f' || P @ f ||^2 = {l2_norm_sqr} should be equal to')
+ print(f' || P @ f ||^2 (alt) = {l2_norm2_sqr}')
+ # M.dot_inner(fc, fc) and M_0.dot_inner(fc, fc) are the same only up to order 1e-15.
+ # Hence, we can't expect err_sqr to be less than tol**2, but only less than tol.
+ assert err_sqr < tol
+
+ print()
+
+
+# ===============================================================================
+# SCRIPT FUNCTIONALITY
+#===============================================================================
+
+if __name__ == "__main__":
+ import sys
+ pytest.main( sys.argv )
diff --git a/feectools/fem/tests/test_spline_histopolation.py b/feectools/fem/tests/test_spline_histopolation.py
index 85872551c..a6e798029 100644
--- a/feectools/fem/tests/test_spline_histopolation.py
+++ b/feectools/fem/tests/test_spline_histopolation.py
@@ -1,3 +1,8 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
import numpy as np
import matplotlib.pyplot as plt
diff --git a/feectools/fem/tests/test_spline_interpolation.py b/feectools/fem/tests/test_spline_interpolation.py
index 406097a0e..1f8a9d0f4 100644
--- a/feectools/fem/tests/test_spline_interpolation.py
+++ b/feectools/fem/tests/test_spline_interpolation.py
@@ -1,11 +1,14 @@
-# coding: utf-8
-# Copyright 2018 Yaman Güçlü
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import time
from feectools.ddm.mpi import mpi as MPI
import numpy as np
import pytest
-import time
from feectools.core.bsplines import make_knots
from feectools.fem.basic import FemField
@@ -17,10 +20,8 @@
from feectools.fem.tests.splines_error_bounds import spline_1d_error_bound
from feectools.fem.tests.analytical_profiles_1d import (AnalyticalProfile1D_Cos, AnalyticalProfile1D_Poly)
#===============================================================================
-@pytest.mark.serial
@pytest.mark.parametrize( "ncells", [1,5,10,23] )
@pytest.mark.parametrize( "degree", range(1,11) )
-
def test_SplineInterpolation1D_exact( ncells, degree ):
domain = [-1.0, 1.0]
@@ -53,10 +54,9 @@ def args_SplineInterpolation1D_cosine():
for degree in range(1,pmax+1):
yield (ncells, degree, periodic)
-@pytest.mark.serial
+
@pytest.mark.parametrize( "ncells,degree,periodic",
args_SplineInterpolation1D_cosine() )
-
def test_SplineInterpolation1D_cosine( ncells, degree, periodic ):
f = AnalyticalProfile1D_Cos()
@@ -78,12 +78,11 @@ def test_SplineInterpolation1D_cosine( ncells, degree, periodic ):
assert max_norm_err < err_bound
#===============================================================================
-@pytest.mark.parallel
+@pytest.mark.mpi
@pytest.mark.parametrize( "nc1", [7,10,23] )
@pytest.mark.parametrize( "nc2", [7,10,23] )
@pytest.mark.parametrize( "deg1", range(1,5) )
@pytest.mark.parametrize( "deg2", range(1,5) )
-
def test_SplineInterpolation2D_parallel_exact( nc1, nc2, deg1, deg2 ):
# Communicator, size, rank
diff --git a/feectools/fem/tests/test_splines.py b/feectools/fem/tests/test_splines.py
index 815ba8233..c56adf5ff 100644
--- a/feectools/fem/tests/test_splines.py
+++ b/feectools/fem/tests/test_splines.py
@@ -1,4 +1,9 @@
-# -*- coding: UTF-8 -*-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+from numpy import linspace
from feectools.fem.basic import FemField
from feectools.fem.splines import SplineSpace
@@ -6,8 +11,6 @@
from feectools.fem.vector import VectorFemSpace
from feectools.ddm.cart import DomainDecomposition
-from numpy import linspace
-
def test_1d_1():
print ('>>> test_1d_1')
diff --git a/feectools/fem/tests/test_tensor.py b/feectools/fem/tests/test_tensor.py
new file mode 100644
index 000000000..60fd8e8d6
--- /dev/null
+++ b/feectools/fem/tests/test_tensor.py
@@ -0,0 +1,232 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+import os
+import contextlib
+from pathlib import Path
+
+import pytest
+
+# Skip entire module if sympde is not available
+try:
+ import sympde
+except ImportError:
+ pytest.skip("sympde not installed", allow_module_level=True)
+
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from PIL import Image
+from mpi4py import MPI
+
+from sympde.topology.domain import Square
+from sympde.topology.space import ScalarFunctionSpace
+from sympde.topology.analytical_mapping import TargetMapping
+
+from psydac.mapping.discrete_gallery import discrete_mapping
+from psydac.api.discretization import discretize
+
+
+#==============================================================================
+# Machinery for comparing a PNG image with a reference on the root MPI process
+#==============================================================================
+
+def similar_images(file1, file2, tolerance=0.01):
+ """
+ Compare two PNG images and check if they are similar within tolerance.
+
+ Parameters
+ ----------
+ file1 : str
+ Path to first PNG file.
+ file2 : str
+ Path to second PNG file.
+ tolerance: float
+ Maximum allowed average difference between pixel values (0-1).
+
+ Returns
+ -------
+ bool :
+ True if images are similar enough, False otherwise.
+ """
+ # Load images and convert to numpy arrays
+ img1 = np.array(Image.open(file1)).astype(float)
+ img2 = np.array(Image.open(file2)).astype(float)
+
+ # Check dimensions match
+ if img1.shape != img2.shape:
+ return False
+
+ # Normalize pixel values to 0-1
+ img1 = img1 / 255.0
+ img2 = img2 / 255.0
+
+ # Calculate mean absolute difference
+ diff = np.mean(np.abs(img1 - img2))
+
+ return diff <= tolerance
+
+
+@contextlib.contextmanager
+def consistent_png_rendering():
+ """
+ Context manager for consistent Matplotlib rendering across platforms.
+ """
+ # Store original settings
+ orig_backend = mpl.get_backend()
+ orig_settings = {
+ 'text.usetex': mpl.rcParams['text.usetex'],
+ 'font.family': mpl.rcParams['font.family'],
+ }
+
+ try:
+ # Use Agg backend (pure python, no GUI)
+ mpl.use('Agg')
+ # Configure settings for consistent rendering
+ mpl.rcParams.update({
+ 'text.usetex': False,
+ 'font.family': 'DejaVu Sans',
+ })
+ yield # Control returns to the with block
+ finally:
+ # Restore original settings
+ mpl.use(orig_backend)
+ mpl.rcParams.update(orig_settings)
+
+
+def compare_figure_to_reference(fig, filename, *, dpi, tol, folder, comm, root):
+ """
+ Compare a matplotlib figure to a reference PNG file on the root MPI process.
+
+ Parameters
+ ----------
+ fig : matplotlib.figure.Figure
+ The figure to compare with the reference image.
+ filename : str
+ Name of the PNG file to save and compare.
+ dpi : int
+ Dots per inch resolution for saving the figure.
+ tol : float
+ Tolerance for image comparison (between 0 and 1).
+ folder : str
+ Name of the folder containing reference images.
+ comm : mpi4py.MPI.Comm
+ MPI communicator object.
+ root : int
+ Rank of the MPI process that should perform the comparison.
+
+ Returns
+ -------
+ bool
+ True if the images are similar enough within the tolerance,
+ False otherwise. The result is broadcast to all MPI processes.
+
+ Notes
+ -----
+ The function saves the figure to a temporary file, compares it with
+ the reference image, and then removes the temporary file. Only the
+ root process performs the actual comparison, but the result is
+ broadcast to all processes.
+ """
+ if comm.rank == root:
+ test_dir = Path(__file__).parent.absolute()
+ file1 = test_dir / filename
+ file2 = test_dir / folder / filename
+ with consistent_png_rendering():
+ fig.savefig(file1, dpi=dpi)
+ close_enough = similar_images(file1, file2, tol)
+ # Clean up the temporary file
+ os.remove(file1)
+ else:
+ close_enough = None
+
+ # Broadcast the boolean result from the root process to all others
+ close_enough = comm.bcast(close_enough, root=root)
+
+ # All MPI processes return the same result
+ return close_enough
+
+#==============================================================================
+# Unit tests
+#==============================================================================
+@pytest.mark.mpi
+@pytest.mark.parametrize('root', ['first', 'last'])
+@pytest.mark.parametrize('kind', ['spline', 'analytical'])
+def test_plot_2d_decomposition(kind, root):
+
+ # MPI communicator
+ mpi_comm = MPI.COMM_WORLD
+ mpi_size = mpi_comm.size
+ mpi_rank = mpi_comm.rank
+
+ # MPI rank which should make the plot
+ if root == 'first':
+ mpi_root = 0
+ elif root == 'last':
+ mpi_root = mpi_size - 1
+ else:
+ raise ValueError(f'root argument has wrong value {root}')
+
+ # Parameters of tensor-product 2D spline space
+ ncells = (6, 9)
+ degree = (2, 2)
+
+ if kind == 'spline':
+ # 2D spline mapping and tensor FEM space (distributed)
+ F, Vh = discrete_mapping('target', ncells=ncells, degree=degree,
+ comm=mpi_comm, return_space=True)
+ elif kind == 'analytical':
+ Omega = Square('Omega', bounds1=(0, 1), bounds2=(0, 2 * np.pi))
+ params = dict(c1=0, c2=0, k=0.3, D=0.2)
+ M = TargetMapping('M', dim=2, **params)
+ domain = M(Omega)
+ V = ScalarFunctionSpace('V', domain)
+
+ # 2D Geometry object
+ domain_h = discretize(domain, ncells=ncells, periodic=(False, True),
+ comm=mpi_comm)
+
+ # 2D spline tensor FEM space (distributed)
+ Vh = discretize(V, domain_h, degree=degree)
+
+ # 2D callable mapping (analytical)
+ F = M.get_callable_mapping()
+ else:
+ raise ValueError(f'kind argument has wrong value {kind}')
+
+ # Name of temporary image file to be compared with reference one
+ filename = f'decomp_{kind}_{mpi_size}_procs.png'
+
+ # Relative tolerance for image comparison
+ RTOL = 0.02
+
+ # Plot 2D decomposition
+ # [1] Run without passing (fig, ax)
+ fig = Vh.plot_2d_decomposition(F, refine=5, mpi_root=mpi_root)
+ assert compare_figure_to_reference(fig, filename, folder='data', dpi=100,
+ tol=RTOL, comm=mpi_comm, root=mpi_root)
+
+ # [2] Run with given (fig, ax), compatible
+ fig2, ax2 = plt.subplots(1, 1) if mpi_rank == mpi_root else (None, None)
+ Vh.plot_2d_decomposition(F, refine=5, fig=fig2, ax=ax2, mpi_root=mpi_root)
+ assert compare_figure_to_reference(fig2, filename, folder='data', dpi=100,
+ tol=RTOL, comm=mpi_comm, root=mpi_root)
+
+ # [3] Run with given (fig, ax), incompatible
+ if mpi_rank == mpi_root:
+ fig3, ax3 = plt.subplots(1, 1)
+ with pytest.raises(AssertionError) as excinfo:
+ Vh.plot_2d_decomposition(F, refine=5, fig=fig2, ax=ax3, mpi_root=mpi_root)
+ assert "Argument `ax` must be in `fig.axes`" in str(excinfo.value)
+ plt.close(fig3)
+ else:
+ Vh.plot_2d_decomposition(F, refine=5, fig=None, ax=None, mpi_root=mpi_root)
+
+#==============================================================================
+if __name__ == '__main__':
+
+ test_plot_2d_decomposition('spline', 'first')
+ test_plot_2d_decomposition('analytical', 'last')
+ plt.show()
diff --git a/feectools/fem/tests/test_vector_spaces.py b/feectools/fem/tests/test_vector_spaces.py
index adac30752..ff1d48507 100644
--- a/feectools/fem/tests/test_vector_spaces.py
+++ b/feectools/fem/tests/test_vector_spaces.py
@@ -1,5 +1,8 @@
-# -*- coding: UTF-8 -*-
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
from numpy import linspace
diff --git a/feectools/fem/tests/utilities.py b/feectools/fem/tests/utilities.py
index 838212942..84be02bd7 100644
--- a/feectools/fem/tests/utilities.py
+++ b/feectools/fem/tests/utilities.py
@@ -1,6 +1,8 @@
-# coding: utf-8
-# Copyright 2018 Yaman Güçlü
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import numpy as np
#===============================================================================
diff --git a/feectools/linalg/direct_solvers.py b/feectools/linalg/direct_solvers.py
index c8f06b809..f9d515a9b 100644
--- a/feectools/linalg/direct_solvers.py
+++ b/feectools/linalg/direct_solvers.py
@@ -1,15 +1,32 @@
-# coding: utf-8
-# Copyright 2018 Jalal Lakhlili, Yaman Güçlü
-
-from abc import abstractmethod
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import numpy as np
from scipy.linalg.lapack import dgbtrf, dgbtrs, sgbtrf, sgbtrs, cgbtrf, cgbtrs, zgbtrf, zgbtrs
-from scipy.sparse import spmatrix
+from scipy.sparse import spmatrix, dia_matrix
from scipy.sparse.linalg import splu
from feectools.linalg.basic import LinearSolver
-__all__ = ('BandedSolver', 'SparseSolver')
+__all__ = ('to_bnd', 'BandedSolver', 'SparseSolver')
+
+#===============================================================================
+def to_bnd(A):
+ """Converts a 1D StencilMatrix to a band matrix"""
+
+ dmat = dia_matrix(A.toarray(), dtype=A.dtype)
+ la = abs(dmat.offsets.min())
+ ua = dmat.offsets.max()
+ cmat = dmat.tocsr()
+
+ A_bnd = np.zeros((1+ua+2*la, cmat.shape[1]), A.dtype)
+
+ for i,j in zip(*cmat.nonzero()):
+ A_bnd[la+ua+i-j, j] = cmat[i,j]
+
+ return A_bnd, la, ua
#===============================================================================
class BandedSolver(LinearSolver):
@@ -58,6 +75,14 @@ def __init__(self, u, l, bmat, transposed=False):
self._space = np.ndarray
self._dtype = bmat.dtype
+ @staticmethod
+ def from_stencil_mat_1d(A):
+ """Converts a 1D StencilMatrix to a BandedSolver."""
+
+ A.remove_spurious_entries()
+ A_bnd, la, ua = to_bnd(A)
+ return BandedSolver(ua, la, A_bnd)
+
@property
def finfo(self):
return self._finfo
diff --git a/feectools/linalg/stencil.py b/feectools/linalg/stencil.py
index 7c9f83e9b..d32504760 100644
--- a/feectools/linalg/stencil.py
+++ b/feectools/linalg/stencil.py
@@ -1,13 +1,13 @@
-# coding: utf-8
-#
-# Copyright 2018 Yaman Güçlü
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import os
import warnings
+from types import MappingProxyType
import numpy as np
-
-from types import MappingProxyType
from scipy.sparse import coo_matrix, diags as sp_diags
from feectools.ddm.mpi import mpi as MPI
@@ -1429,11 +1429,12 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None):
Returns
-------
StencilDiagonalMatrix
- The matrix which contains the main diagonal of self (or its inverse).
+ The matrix which contains the main diagonal of self (or its inverse (square root)).
"""
- # Check `inverse` argument
+ # Check `inverse` and `sqrt` argument
assert isinstance(inverse, bool)
+ assert isinstance(sqrt, bool)
# Determine domain and codomain of the StencilDiagonalMatrix
V, W = self.domain, self.codomain
@@ -1446,7 +1447,6 @@ def diagonal(self, *, inverse = False, sqrt = False, out = None):
assert out.domain is V
assert out.codomain is W
-
# Extract diagonal data from self and identify output array
diagonal_indices = self._get_diagonal_indices()
diag = self._data[diagonal_indices]
@@ -2085,22 +2085,15 @@ def transpose(self, *, conjugate=False, out=None):
assert isinstance(out, StencilDiagonalMatrix)
assert out.domain is self.codomain
assert out.codomain is self.domain
-
- if not (conjugate and self.dtype is complex):
-
- if out is None:
- data = self._data.copy()
+ if conjugate and self.dtype is complex:
+ np.conjugate(self._data, out=out._data, casting='no')
else:
np.copyto(out._data, self._data, casting='no')
-
else:
-
- if out is None:
+ if conjugate and self.dtype is complex:
data = np.conjugate(self._data, casting='no')
else:
- np.conjugate(self._data, out=out._data, casting='no')
-
- if out is None:
+ data = self._data.copy()
out = StencilDiagonalMatrix(self.codomain, self.domain, data)
return out
@@ -2124,16 +2117,21 @@ def copy(self, *, out=None):
return out
- def diagonal(self, *, inverse = False, out = None):
+ def diagonal(self, *, inverse = False, sqrt = False, out = None):
"""
Get the coefficients on the main diagonal as a StencilDiagonalMatrix object.
- In the default case (inverse=False, out=None) self is returned.
+ In the default case (inverse=False, sqrt=False, out=None) self is returned.
Parameters
----------
inverse : bool
If True, get the inverse of the diagonal. (Default: False).
+ Can be combined with sqrt to get the inverse square root.
+
+ sqrt : bool
+ If True, get the square root of the diagonal. (Default: False).
+ Can be combined with inverse to get the inverse square root.
out : StencilDiagonalMatrix
If provided, write the diagonal entries into this matrix. (Default: None).
@@ -2141,11 +2139,12 @@ def diagonal(self, *, inverse = False, out = None):
Returns
-------
StencilDiagonalMatrix
- Either self, or another StencilDiagonalMatrix with the diagonal inverse.
+ Either self, or another StencilDiagonalMatrix with the diagonal (or its inverse (square root)).
"""
- # Check `inverse` argument
+ # Check `inverse` and `sqrt` argument
assert isinstance(inverse, bool)
+ assert isinstance(sqrt, bool)
# Determine domain and codomain of the `out` matrix
V, W = self.domain, self.codomain
@@ -2161,10 +2160,16 @@ def diagonal(self, *, inverse = False, out = None):
assert out.codomain is W
data = out._data
+ diag = self._data
+
# Calculate entries, or set `out=self` in default case
if inverse:
data = np.divide(1, diag, out=data)
- elif out:
+ if sqrt:
+ data = np.sqrt(data, out=data)
+ elif sqrt:
+ data = np.sqrt(diag, out=data)
+ elif out not in (None, self):
np.copyto(data, diag)
else:
out = self
diff --git a/feectools/linalg/tests/test_fft.py b/feectools/linalg/tests/test_fft.py
index f33475f18..935d0818c 100644
--- a/feectools/linalg/tests/test_fft.py
+++ b/feectools/linalg/tests/test_fft.py
@@ -1,3 +1,8 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
import scipy.fft as scifft
import numpy as np
@@ -98,7 +103,7 @@ def test_kron_fft_ser(seed, params, ffttype):
@pytest.mark.parametrize( 'seed', [0, 2] )
@pytest.mark.parametrize( 'params', [([16], [2], [False]), ([16,18], [2,3], [False,True]), ([16,18,37], [2,3,7], [False,True,False])] )
@pytest.mark.parametrize( 'ffttype', ['fft', 'ifft', 'dct', 'idct', 'dst', 'idst'] )
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_kron_fft_par(seed, params, ffttype):
method_test(seed, MPI.COMM_WORLD, params, *decode_fft_type(ffttype), verbose=False)
diff --git a/feectools/linalg/tests/test_stencil_interface_matrix.py b/feectools/linalg/tests/test_stencil_interface_matrix.py
index b3427d029..bfafd71b9 100644
--- a/feectools/linalg/tests/test_stencil_interface_matrix.py
+++ b/feectools/linalg/tests/test_stencil_interface_matrix.py
@@ -1,8 +1,10 @@
-# -*- coding: UTF-8 -*-
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
import numpy as np
-from random import random
from feectools.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix, StencilInterfaceMatrix
from feectools.api.settings import *
@@ -223,7 +225,7 @@ def test_stencil_interface_matrix_3d_serial_init(dtype, n1, n2, n3, p1, p2, p3,
(12,12,1,1, 3023467041788.0),
(12,12,2,2, 19555497680544.0),
(12,12,3,3, 62573623909332.0)])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_interface_matrix_2d_parallel_dot(n1, n2, p1, p2, expected):
from feectools.ddm.mpi import mpi as MPI
diff --git a/feectools/linalg/tests/test_stencil_vector.py b/feectools/linalg/tests/test_stencil_vector.py
index 40f9b25e7..15901fe3c 100644
--- a/feectools/linalg/tests/test_stencil_vector.py
+++ b/feectools/linalg/tests/test_stencil_vector.py
@@ -1,5 +1,8 @@
-# coding: utf-8
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
import numpy as np
@@ -526,7 +529,7 @@ def test_stencil_vector_2d_serial_update_ghost_region_interior(dtype, n1, n2, p1
@pytest.mark.parametrize('n1', [12, 22])
@pytest.mark.parametrize('p1', [1, 3])
@pytest.mark.parametrize('s1', [1, 2])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_1d_parallel_init(dtype, n1, p1, s1, P1=True):
comm = MPI.COMM_WORLD
@@ -561,7 +564,7 @@ def test_stencil_vector_1d_parallel_init(dtype, n1, p1, s1, P1=True):
@pytest.mark.parametrize('p2', [3])
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('s2', [2])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_2d_parallel_init(dtype, n1, n2, p1, p2, s1, s2, P1=True, P2=False):
comm = MPI.COMM_WORLD
@@ -598,7 +601,7 @@ def test_stencil_vector_2d_parallel_init(dtype, n1, n2, p1, p2, s1, s2, P1=True,
@pytest.mark.parametrize('s2', [2])
@pytest.mark.parametrize('P1', [True, False])
@pytest.mark.parametrize('P2', [True])
-@pytest.mark.parallel
+@pytest.mark.mpi
@pytest.mark.petsc
def test_stencil_vector_2d_parallel_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P2):
@@ -641,7 +644,7 @@ def test_stencil_vector_2d_parallel_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P
@pytest.mark.parametrize('p1', [1, 3])
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('P1', [True, False])
-@pytest.mark.parallel
+@pytest.mark.mpi
@pytest.mark.petsc
def test_stencil_vector_1d_parallel_topetsc(dtype, n1, p1, s1, P1):
@@ -692,7 +695,7 @@ def test_stencil_vector_1d_parallel_topetsc(dtype, n1, p1, s1, P1):
@pytest.mark.parametrize('P2', [True])
@pytest.mark.parametrize('P3', [False])
-@pytest.mark.parallel
+@pytest.mark.mpi
@pytest.mark.petsc
def test_stencil_vector_3d_parallel_topetsc(dtype, n1, n2, n3, p1, p2, p3, s1, s2, s3, P1, P2, P3):
@@ -742,7 +745,7 @@ def test_stencil_vector_3d_parallel_topetsc(dtype, n1, n2, n3, p1, p2, p3, s1, s
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('s2', [3])
@pytest.mark.parametrize('s3', [1])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_3d_parallel_init(dtype, n1, n2, n3, p1, p2, p3, s1, s2, s3, P1=True, P2=False, P3=True):
comm = MPI.COMM_WORLD
@@ -778,7 +781,7 @@ def test_stencil_vector_3d_parallel_init(dtype, n1, n2, n3, p1, p2, p3, s1, s2,
@pytest.mark.parametrize('p2', [2])
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('s2', [2])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_2d_parallel_toarray(dtype, n1, n2, p1, p2, s1, s2, P1=True, P2=False):
# Create domain decomposition
comm = MPI.COMM_WORLD
@@ -845,7 +848,7 @@ def test_stencil_vector_2d_parallel_toarray(dtype, n1, n2, p1, p2, s1, s2, P1=Tr
@pytest.mark.parametrize('s2', [1])
@pytest.mark.parametrize('P1', [True, False])
@pytest.mark.parametrize('P2', [True])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_2d_parallel_array_to_psydac(dtype, n1, n2, p1, p2, s1, s2, P1, P2):
npts = [n1, n2]
@@ -902,7 +905,7 @@ def test_stencil_vector_2d_parallel_array_to_psydac(dtype, n1, n2, p1, p2, s1, s
@pytest.mark.parametrize('p2', [2])
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('s2', [1])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_2d_parallel_dot(dtype, n1, n2, p1, p2, s1, s2, P1=True, P2=False):
comm = MPI.COMM_WORLD
@@ -969,7 +972,7 @@ def test_stencil_vector_2d_parallel_dot(dtype, n1, n2, p1, p2, s1, s2, P1=True,
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('s2', [1, 2])
@pytest.mark.parametrize('s3', [1])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_3d_parallel_dot(dtype, n1, n2, n3, p1, p2, p3, s1, s2, s3, P1=True, P2=False, P3=True):
comm = MPI.COMM_WORLD
diff --git a/feectools/linalg/tests/test_stencil_vector_space.py b/feectools/linalg/tests/test_stencil_vector_space.py
index 3bb9122b9..59ab5e75d 100644
--- a/feectools/linalg/tests/test_stencil_vector_space.py
+++ b/feectools/linalg/tests/test_stencil_vector_space.py
@@ -1,3 +1,8 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
import pytest
import numpy as np
@@ -276,7 +281,7 @@ def test_stencil_vector_space_2D_serial_set_interface(dtype, n1, n2, p1, p2, s1,
@pytest.mark.parametrize('p1', [1, 2])
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('P1', [True, False])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_space_1d_parallel_init(dtype, n1, p1, s1, P1):
@@ -321,7 +326,7 @@ def test_stencil_vector_space_1d_parallel_init(dtype, n1, p1, s1, P1):
@pytest.mark.parametrize('s2', [2])
@pytest.mark.parametrize('P1', [True, False])
@pytest.mark.parametrize('P2', [True])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_space_2d_parallel_init(dtype, n1, n2, p1, p2, s1, s2, P1, P2):
@@ -366,7 +371,7 @@ def test_stencil_vector_space_2d_parallel_init(dtype, n1, n2, p1, p2, s1, s2, P1
@pytest.mark.parametrize('s1', [1, 2])
@pytest.mark.parametrize('s2', [1, 2])
@pytest.mark.parametrize('s3', [1])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_space_3d_parallel_init(dtype, n1, n2, n3, p1, p2, p3, s1, s2, s3, P1=True, P2=False, P3=True):
@@ -404,7 +409,7 @@ def test_stencil_vector_space_3d_parallel_init(dtype, n1, n2, n3, p1, p2, p3, s1
@pytest.mark.parametrize('dtype', [float, complex])
@pytest.mark.parametrize('n1', [15])
@pytest.mark.parametrize('n2', [20])
-@pytest.mark.parallel
+@pytest.mark.mpi
def test_stencil_vector_space_2D_parallel_parent(dtype, n1, n2, P1=True, P2=False):
diff --git a/feectools/linalg/tests/utilities.py b/feectools/linalg/tests/utilities.py
new file mode 100644
index 000000000..755e22197
--- /dev/null
+++ b/feectools/linalg/tests/utilities.py
@@ -0,0 +1,89 @@
+import numpy as np
+
+from sympde.topology import Mapping
+
+from psydac.linalg.basic import LinearOperator
+from psydac.linalg.block import BlockVectorSpace
+
+__all__ = (
+ 'SquareTorus',
+ 'Annulus',
+ 'SinMapping1D',
+ 'check_linop_equality_using_rng'
+)
+
+class SquareTorus(Mapping):
+ _expressions = {'x': 'x1 * cos(x2)',
+ 'y': 'x1 * sin(x2)',
+ 'z': 'x3'}
+ _ldim = 3
+ _pdim = 3
+
+
+class Annulus(Mapping):
+ _expressions = {'x': 'x1 * cos(x2)',
+ 'y': 'x1 * sin(x2)'}
+ _ldim = 2
+ _pdim = 2
+
+
+class SinMapping1D(Mapping):
+ _expressions = {'x': 'sin((pi/2)*x1)'}
+ _ldim = 1
+ _pdim = 1
+
+
+def check_linop_equality_using_rng(A, B, tol=1e-15):
+ """
+ A simple tool to check with almost certainty that two linear operators are
+ identical, by applying them repeatedly to random vectors with entries drawn
+ uniformly in [0, 1).
+
+ Parameters
+ ----------
+ A : LinearOperator
+ First linear operator to compare.
+ B : LinearOperator
+ Second linear operator to compare.
+ tol : float, optional
+ Tolerance for the comparison on the relative 2-norm of the error.
+ The default is 1e-15.
+
+ Raises
+ ------
+ AssertionError
+ If the two operators are incompatible, or found to be different. This
+ behavior is intended for using this function in unit tests with Pytest.
+
+ """
+
+ assert isinstance(A, LinearOperator)
+ assert isinstance(B, LinearOperator)
+ assert A.domain is B.domain
+ assert A.codomain is B.codomain
+
+ rng = np.random.default_rng(42)
+
+ x = A.domain.zeros()
+ y1 = A.codomain.zeros()
+ y2 = A.codomain.zeros()
+
+ n = 10
+
+ for _ in range(n):
+
+ x *= 0.
+
+ if isinstance(A.domain, BlockVectorSpace):
+ for block in x.blocks:
+ rng.random(size=block._data.shape, dtype="float64", out=block._data)
+ else:
+ rng.random(size=x._data.shape, dtype="float64", out=x._data)
+
+ A.dot(x, out=y1)
+ B.dot(x, out=y2)
+
+ diff = y1 - y2
+ scaled_err_sqr = diff.inner(diff) / diff.space.dimension**2
+
+ assert scaled_err_sqr < tol**2
diff --git a/feectools/pytest.toml b/feectools/pytest.toml
new file mode 100644
index 000000000..7b525ecca
--- /dev/null
+++ b/feectools/pytest.toml
@@ -0,0 +1,52 @@
+# this file shows to pytest what to collect
+# here we avoid collecting test classes (which are not used in PSYDAC)
+# this is to avoid getting the warning on TestFunction
+# TODO can we exlude TestFunction from python_classes pattern?
+[pytest]
+minversion = "9.0"
+addopts = ["--strict-markers"]
+markers = [
+ "mpi: parallel test to be run using 'mpiexec' or 'mpirun'",
+ "petsc: test requiring a working PETSc installation with petsc4py Python bindings",
+ "pyccel: test for checking Pyccel setup on machine",
+]
+python_files = ["test_*.py"]
+python_classes = []
+python_functions = ["test_*"]
+
+[coverage.run]
+branch = true
+omit = [
+ # Exclude pyccelised kernels
+ "*/__psydac__/*",
+
+ # Examples don't need to be covered
+ "*/examples/*",
+
+ # Unit tests shouldn't be included
+ "*/tests/*",
+]
+
+[coverage.report]
+# Regexes for lines to exclude from consideration
+exclude_also = [
+ # Don't complain about missing debug-only code:
+ "def __repr__",
+ "if self\\.debug",
+
+ # Don't complain if tests don't hit defensive assertion code:
+ "raise AssertionError",
+ "raise NotImplementedError",
+
+ # Don't complain if non-runnable code isn't run:
+ "if False:",
+ "if __name__ == .__main__.:",
+
+ # Don't complain about abstract methods, they aren't run:
+ "@(abc\\.)?abstractmethod",
+]
+# Ignore source code that can’t be found, emitting a warning instead of an exception.
+ignore_errors = true
+
+[coverage.html]
+directory = "coverage_html_report"
diff --git a/feectools/utilities/quadratures.py b/feectools/utilities/quadratures.py
index fa99a273f..535561105 100644
--- a/feectools/utilities/quadratures.py
+++ b/feectools/utilities/quadratures.py
@@ -1,17 +1,17 @@
-# -*- coding: UTF-8 -*-
-#! /usr/bin/python
-
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
"""
This module contains some routines to generate quadrature points in 1D
it has also a routine uniform, which generates uniform points
with weights equal to 1
"""
-import numpy as np
-
from math import cos, pi
-from numpy import zeros
+import numpy as np
__all__ = ('gauss_legendre', 'gauss_lobatto', 'quadrature')
@@ -55,8 +55,8 @@ def legendre(t, m):
dp = m*(p0 - t*p1)/(1.0 - t**2)
return p1, dp
- A = zeros(m)
- x = zeros(m)
+ A = np.zeros(m)
+ x = np.zeros(m)
nRoots = (m + 1) // 2 # Number of non-neg. roots
for i in range(nRoots):
t = cos(pi*(i + 0.75)/(m + 0.5)) # Approx. root
diff --git a/feectools/utilities/utils.py b/feectools/utilities/utils.py
index e91b45184..6389d2f34 100644
--- a/feectools/utilities/utils.py
+++ b/feectools/utilities/utils.py
@@ -1,9 +1,11 @@
-# coding: utf-8
-#
-# Copyright 2018 Yaman Güçlü
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+from numbers import Number
import numpy as np
-from numbers import Number
__all__ = (
'refine_array_1d',
diff --git a/feectools/version.py b/feectools/version.py
index a4e2017f0..3652431d0 100644
--- a/feectools/version.py
+++ b/feectools/version.py
@@ -1 +1,35 @@
-__version__ = "0.1"
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+def extract_version() -> str:
+ """
+ Returns either the version of the installed package (e.g. "0.17.1") or the
+ one found in pyproject.toml (e.g. "0.17.1-dev (at /project/location)").
+
+ Returns
+ -------
+ str
+ The package version.
+
+ See also
+ --------
+ https://stackoverflow.com/a/76206192
+
+ """
+ from contextlib import suppress
+ from pathlib import Path
+
+ with suppress(FileNotFoundError, StopIteration):
+ root_dir = Path(__file__).parent.parent
+ with open(root_dir / "pyproject.toml", encoding="utf-8") as pyproject_toml:
+ version_line = next(line for line in pyproject_toml if line.startswith("version"))
+ version = version_line.split("=")[1].strip("'\"\n ")
+ return f"{version}-dev (at {root_dir})"
+
+ import importlib.metadata
+ return importlib.metadata.version(__package__)
+
+
+__version__ = extract_version()
diff --git a/meson.build b/meson.build
new file mode 100644
index 000000000..c5c9c6760
--- /dev/null
+++ b/meson.build
@@ -0,0 +1,52 @@
+project(
+ 'psydac',
+ 'c', 'fortran',
+ meson_version: '>=1.1.0',
+)
+
+pyccel_main_language = get_option('pyccel_language')
+
+message('Translating kernels to ' + pyccel_main_language)
+
+run_command('pyccel', 'make', '-g', 'psydac/**/*_kernels.py', '--convert-only', '--language=' + pyccel_main_language,
+ check: true)
+
+find = find_program('find', required: true)
+
+py = import('python').find_installation(modules: ['numpy'], pure: false)
+
+# Math dependencies
+cc = meson.get_compiler('c')
+m_dep = cc.find_library('m', required : false)
+
+pyccel_meson_math_file = run_command(find, '__pyccel__/math', '-name', 'meson.build',
+ check: true).stdout().strip()
+pyccel_meson_cwrapper_file = run_command(find, '__pyccel__/cwrapper', '-name', 'meson.build',
+ check: true).stdout().strip()
+pyccel_meson_files = run_command(find, '__pyccel__/psydac', '-name', 'meson.build',
+ check: true).stdout().split()
+
+sed = find_program('sed', required: true)
+
+run_command(sed, '-i.swp', 's/install: true//g', pyccel_meson_math_file, pyccel_meson_cwrapper_file,
+ check: true)
+
+foreach f : pyccel_meson_files
+ f_parts = f.split('/')
+ rel_f = f_parts.get(1)
+ foreach i : range(2, f_parts.length()-1)
+ rel_f = rel_f / f_parts.get(i)
+ endforeach
+ run_command(sed, '-i.swp', 's/install_dir: .*/subdir: \'' + rel_f.replace('/','\/') + '\',/g', f,
+ check: true)
+endforeach
+
+
+igakit_proj = subproject('igakit')
+
+subdir('__pyccel__/math')
+subdir('__pyccel__/cwrapper')
+subdir('__pyccel__/psydac')
+
+# Install the pure-Python package itself
+install_subdir('psydac', install_dir: py.get_install_dir())
diff --git a/meson.options b/meson.options
new file mode 100644
index 000000000..77647aafb
--- /dev/null
+++ b/meson.options
@@ -0,0 +1,2 @@
+option('pyccel_language', type : 'string', value : 'fortran', description : 'The language that Pyccel translates to [default: Fortran].')
+
diff --git a/psydac_accelerate.py b/psydac_accelerate.py
deleted file mode 100644
index 6840a2773..000000000
--- a/psydac_accelerate.py
+++ /dev/null
@@ -1,58 +0,0 @@
-
-import argparse
-import os
-from subprocess import run as sub_run
-import shutil
-
-
-
-#The purpose of this file is to be launched after an editable installation of Psydac, to pyccelise all the Psydac kernels.
-#This file is useless during a classic installation because the kernels are already pyccelised in the construction folder.
-
-
-
-parser = argparse.ArgumentParser(
- formatter_class=argparse.ArgumentDefaultsHelpFormatter,
- description="Accelerate all computational kernels in Psydac using Pyccel (editable install only)"
-)
-
-# Add Argument --language at the pyccel command
-parser.add_argument('--language',
- type=str,
- default='fortran',
- choices=['fortran', 'c'],
- action='store',
- dest='language',
- help='Language used to pyccelize all the _kernels files'
- )
-
-# Add flag --openmp at the pyccel command
-parser.add_argument('--openmp',
- default=False,
- action='store_true',
- dest='openmp',
- help="Use OpenMP multithreading in generated code."
- )
-
-# Read input arguments
-args = parser.parse_args()
-
-# get the absolute path to the psydac directory
-psydac_path = os.path.dirname(os.path.abspath(__file__))+'/feectools'
-
-print("\nNOTE: This script should only be used if Psydac was installed in editable mode.\n")
-
-# Define all the parameters of the command in the parameters array
-parameters = ['--language', args.language]
-
-# check if the flag --openmp is passed and add it to the argument if it's the case
-if args.openmp:
- parameters.append('--openmp')
-
-# search in psydac/psydac folder all the files ending with the tag _kernels.py
-for path, subdirs, files in os.walk(psydac_path):
- for name in files:
- if name.endswith('_kernels.py'):
- print(' Pyccelize file: ' + os.path.join(path, name))
- sub_run([shutil.which('pyccel'), os.path.join(path, name), *parameters], shell=False)
-print()
diff --git a/pyproject.toml b/pyproject.toml
index dafc3f469..43ea27331 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -107,3 +107,15 @@ ignore_errors = true
[tool.coverage.html]
directory = "coverage_html_report"
+
+[tool.pytest.ini_options]
+markers = [
+ "mpi: mark test as requiring MPI",
+ "petsc: mark test as requiring PETSc",
+ "parallel: mark test as parallel",
+]
+testpaths = ["feectools"]
+norecursedirs = [".git", ".venv", "env", "__pycache__", "*.egg-info"]
+python_files = "test_*.py"
+python_classes = "Test*"
+python_functions = "test_*"
diff --git a/pytest.ini b/pytest.ini
index 91ef13cfa..c4a9a600f 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -1,16 +1,9 @@
-# this file shows to pytest what to collect
-# here we avoid collecting test classes (which are not used in PSYDAC)
-# this is to avoid getting the warning on TestFunction
-# TODO can we exlude TestFunction from python_classes pattern?
[pytest]
-minversion = 4.5
-addopts = --strict-markers
markers =
- serial: single-process test,
- parallel: test to be run using 'mpiexec',
- petsc: test requiring a working PETSc installation with petsc4py Python bindings
- pyccel: test for checking Pyccel setup on machine
-
+ mpi: mark test as requiring MPI
+ petsc: mark test as requiring PETSc
+ parallel: mark test as parallel
+testpaths = feectools
python_files = test_*.py
-python_classes =
+python_classes = Test*
python_functions = test_*
diff --git a/scripts/add_header.py b/scripts/add_header.py
new file mode 100644
index 000000000..18e49dd29
--- /dev/null
+++ b/scripts/add_header.py
@@ -0,0 +1,47 @@
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+
+def add_header_to_python_files(directory, string_to_add, *, dry_run=True):
+
+ print(f'Adding header to all Python files:')
+ print(f'{header}')
+
+ import os
+ for root, dirs, files in os.walk(directory):
+ if os.path.basename(root) in ['__pycache__', '__pyccel__', '__psydac__']:
+ print(f'- Skipping folder: {root}')
+ continue
+ for file in files:
+ if file.endswith('.py'):
+ file_path = os.path.join(root, file)
+ print(f'+ Processing file: {file_path}')
+ with open(file_path, 'r+') as f:
+ content = f.read()
+ f.seek(0, 0)
+ if not dry_run:
+ f.write(string_to_add.rstrip('\r\n') + '\n' + content)
+ print()
+
+#==============================================================================
+if __name__ == '__main__':
+
+ # Test the function
+# header = '# This is a test string'
+
+ # Add license header
+ header = """
+#---------------------------------------------------------------------------#
+# This file is part of PSYDAC which is released under MIT License. See the #
+# LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE #
+# for full license details. #
+#---------------------------------------------------------------------------#
+"""
+ header = header.strip()
+
+ dry_run = True
+ add_header_to_python_files('../psydac' , header, dry_run=dry_run)
+ add_header_to_python_files('../examples' , header, dry_run=dry_run)
+ add_header_to_python_files('../mesh' , header, dry_run=dry_run)
diff --git a/mpi_tester.py b/scripts/mpi_tester.py
similarity index 100%
rename from mpi_tester.py
rename to scripts/mpi_tester.py
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 9f5944ec7..000000000
--- a/setup.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import logging
-import os
-from shutil import which
-from subprocess import PIPE, STDOUT # nosec B404
-from subprocess import run as sub_run
-from setuptools import setup
-from setuptools.command.build_py import build_py
-
-
-class BuildPyCommand(build_py):
- """Custom build command to pyccelise _kernels files in the build directory."""
-
- # Rewrite the build_module function to copy each module in the build
- # repository and pyccelise the modules ending with _kernels
- def build_module(self, module, module_file, package):
- outfile, copied = super().build_module(module, module_file, package)
-
- # This part check if the module is pyccelisable and pyccelise it in
- # case
- # if module.endswith('_kernels'):
- # self.announce(f"\nPyccelising [{module}] ...", level=logging.INFO)
- # # TODO --openmp doesn't work with c-compilation
- # pyccel = sub_run([which('pyccel'), outfile, '--language', 'c'],
- # stdout=PIPE, stderr=STDOUT,
- # text=True, shell=False, check=True) # nosec B603
- # self.announce(pyccel.stdout, level=logging.INFO)
-
- return outfile, copied
-
- def run(self):
- super().run()
-
- # Remove __pyccel__ directories
- sub_run([which('pyccel-clean'), self.build_lib], shell=False, check=True) # nosec B603, B607
-
- # Remove useless .lock files
- for path, subdirs, files in os.walk(self.build_lib):
- for name in files:
- if name == '.lock_acquisition.lock':
- os.remove(os.path.join(path, name))
-
-
-setup(
- cmdclass={
- 'build_py': BuildPyCommand,
- },
-)
diff --git a/subprojects/igakit b/subprojects/igakit
new file mode 160000
index 000000000..8f8a64991
--- /dev/null
+++ b/subprojects/igakit
@@ -0,0 +1 @@
+Subproject commit 8f8a64991b930c62f0ab7d459a55cc17e0a12a22