diff --git a/examples/vector_potential_3d.py b/examples/vector_potential_3d.py new file mode 100644 index 000000000..5be586037 --- /dev/null +++ b/examples/vector_potential_3d.py @@ -0,0 +1,287 @@ +#---------------------------------------------------------------------------# +# 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 +import numpy as np + +from sympde.calculus import inner +from sympde.expr import integral, BilinearForm +from sympde.topology import elements_of, Derham, Mapping, Cube + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.linalg.basic import IdentityOperator +from psydac.linalg.solvers import inverse + +def compute_vector_potential_3d(b, derham_h): + """ + Computes a weak-divergence-free vector potential A in a subspace of H_0(curl) of a divergence-free function B in a subspace of H_0(div). + This example highlights the usage of LST preconditioners and Dirichlet projectors. + + Parameters + ---------- + b : psydac.linalg.block.BlockVector + Coefficient vector of a divergence-free function belonging to the third space (H_0(div)) + of the discrete de Rham sequence derham_h. + + derham_h : psydac.api.feec.DiscreteDeRham + b belongs to the coefficient space of the third space V2_h of derham_h. + + Returns + ------- + a : psydac.linalg.block.BlockVector + Coefficient vector of the weak-divergence-free vector potential of the function corresponding + to the coefficient vector b. + + Notes + ----- + We denote the function spaces of a 3D de Rham sequence satisfying + homogeneous Dirichlet boundary conditions (DBCs) as V0h, V1h, V2h and V3h. + Then the problem + + Given B in V2h s.th. div(B)=0, find A in V1h s.th. + curl(A) = B, + weak-div(A) = 0, + + can be equivalently stated as the Hodge-Laplace problem + + Given B in V2h s.th. div(B)=0, find A in V1h s.th. + weak-curl(curl(A)) - grad(weak-div(A)) = weak-curl(B). + + The corresponding variational formulation reads + + Given B in V2h s.th. div(B)=0, find A in V1h s.th. + (curl(v), curl(A)) + (weak-div(v), weak-div(A)) = (curl(v), B) for all v in V1h. + + The weak-divergence is defined implicitly by the set of equations + + (grad(u), v) = - (u, weak-div(v)) for all u in V0h, v in V1h. + + Hence, in terms of matrices (G for grad, wD for weak-divergence) + + u^T G^T M1 v = - u^T M0 wD v <=> wD = - (M0)^-1 G^T M1. + + Thus, not taking into account boundary conditions, the full system of equations reads + + v^T C^T M2 C A + v^T M1 G (M0)^-1 M0 (M0)^-1 G^T M1 A = v^T C^T M2 B for all v in V1h + + <=> ( C^T M2 C + M1 G (M0)^-1 G^T M1 ) A = C^T M2 B. + + Note: (M0)^-1 here is the inverse "mass matrix of functions satisfying hom. DBCs", + not the inverse of the entire M0 "mass matrix of all functions". + + This example highlights the importance of choosing the correct inverse, + choosing the correct preconditioner for this correct inverse, + and applying the projection method using `DirichletProjector`s to solve this discrete problem. + """ + + assert b.space is derham_h.spaces[2].coeff_space + + # ----- Obtain standard objects: domain, domain_h, derham, function spaces, mass matrices, derivative operators, ... + + domain_h = derham_h.domain_h + domain = domain_h.domain + derham = Derham(domain) + + V0, V1, V2, V3 = derham.spaces + V0h, V1h, V2h, V3h = derham_h.spaces + V0cs, V1cs, V2cs, V3cs = [Vh.coeff_space for Vh in derham_h.spaces] + + u0, v0 = elements_of(V0, names='u0, v0') + u1, v1 = elements_of(V1, names='u1, v1') + u2, v2 = elements_of(V2, names='u2, v2') + u3, v3 = elements_of(V3, names='u3, v3') + + G, C, D = derham_h.derivatives(kind='linop') + + a0 = BilinearForm((u0, v0), integral(domain, u0*v0)) + a1 = BilinearForm((u1, v1), integral(domain, inner(u1, v1))) + a2 = BilinearForm((u2, v2), integral(domain, inner(u2, v2))) + a3 = BilinearForm((u3, v3), integral(domain, u3*v3)) + + t0 = time.time() + a0h = discretize(a0, domain_h, (V0h, V0h), backend=backend) + a1h = discretize(a1, domain_h, (V1h, V1h), backend=backend) + a2h = discretize(a2, domain_h, (V2h, V2h), backend=backend) + a3h = discretize(a3, domain_h, (V3h, V3h), backend=backend) + t1 = time.time() + print() + print(f'Mass matrix bilinear forms discretized in {t1-t0:.3g}') + + t0 = time.time() + M0 = a0h.assemble() + M1 = a1h.assemble() + M2 = a2h.assemble() + M3 = a3h.assemble() + t1 = time.time() + print(f'Mass matrices assembled in {t1-t0:.3g}') + + # ----- Sanity check: b should be divergence-free. The method will convergence even if it is not, + # but b = curl(a) won't hold. + divB = D @ b + l2norm_divB = np.sqrt(M3.dot_inner(divB, divB)) + assert l2norm_divB < 1e-10, 'The coefficient vector passed to this function should belong to a divergence-free function.' + + # ----- We need to obtain an inverse M0 object in order to assemble the system matrix + DP0, DP1, _, _ = derham_h.dirichlet_projectors(kind='linop') # <- Dirichlet boundary projectors for the projection method + I0 = IdentityOperator(V0cs) + + # Option1: Modified mass matrix of functions satisfying hom. DBCs + M0_0 = DP0 @ M0 @ DP0 + (I0 - DP0) + t0 = time.time() + M0_0_pc, = derham_h.LST_preconditioners(M0=M0, hom_bc=True) + t1 = time.time() + M0_0_inv = inverse(M0_0, 'pcg', pc=M0_0_pc, maxiter=1000, tol=1e-15) + + # Option2: Inverse of entire mass matrix + t2 = time.time() + M0_pc, = derham_h.LST_preconditioners(M0=M0, hom_bc=False) + t3 = time.time() + M0_inv = inverse(M0, 'pcg', pc=M0_pc, maxiter=1000, tol=1e-15) + + print(f'M0 and M0_0 preconditioner obtained in {((t3-t2)+(t1-t0)):.3g}') + print() + + # ----- There are now (at least) three ways to assemble the system matrix + + # Option1: The correct option. Using M0_0_inv and a projection operator immediately before M0_0_inv + S1 = C.T @ M2 @ C + M1 @ G @ M0_0_inv @ DP0 @ G.T @ M1 + # Note that not including the projector here, i.e., using the following stiffness matrix + #S1 = C.T @ M2 @ C + M1 @ G @ M0_0_inv @ G.T @ M1 + # results in a wrong solution (and is super slow) + + # Option2: The entirely wrong option. Using M0_inv and no additional projector + S2 = C.T @ M2 @ C + M1 @ G @ M0_inv @ G.T @ M1 + + # Option3: The better of the wrong options. Using M0_inv, but an additional projector + S3 = C.T @ M2 @ C + M1 @ G @ M0_inv @ DP0 @ G.T @ M1 + + # ----- We assemble the rhs vector + rhs = C.T @ M2 @ b + + # ----- We can now solve the three resulting systems + I1 = IdentityOperator(V1cs) + + maxiter = 5000 + tol = 1e-10 + + vector_potential_list = [] + timings_list = [] + info_list = [] + + for S in [S1, S2, S3]: + # Apply the projection method + S_bc = DP1 @ S @ DP1 + (I1 - DP1) + rhs_bc = DP1 @ rhs + + # Obtain a solver for the system matrix + S_bc_inv = inverse(S_bc, 'cg', maxiter=maxiter, tol=tol) + + # Solve the system + t0 = time.time() + a = S_bc_inv @ rhs_bc + t1 = time.time() + info = S_bc_inv.get_info() + + vector_potential_list.append(a) + timings_list.append(t1-t0) + info_list.append(info) + + # ----- Analyse the results + + # Option1 is the fastest, and accurate. + a1 = vector_potential_list[0] + info1 = info_list[0] + diff = b - C @ a1 + l2norm_diff = np.sqrt(M2.dot_inner(diff, diff)) + print(f' ----- Option1: C.T @ M2 @ C + M1 @ G @ M0_0_inv @ DP0 @ G.T @ M1 -----') + print() + print(f' || B - curl(A) || = {l2norm_diff:.3g}') + print(f' Time : {timings_list[0]:.3g}') + print(f' Niter : {info1["niter"]}') + print(f' Convergence : {info1["success"]}') + print(f' Res. Norm : {info1["res_norm"]:.3g}') + print() + + # Option2 is not much slower, but delivers an entirely wrong solution + a2 = vector_potential_list[1] + info2 = info_list[1] + diff = b - C @ a2 + l2norm_diff = np.sqrt(M2.dot_inner(diff, diff)) + print(f' ----- Option2: C.T @ M2 @ C + M1 @ G @ M0_inv @ G.T @ M1 -----') + print() + print(f' || B - curl(A) || = {l2norm_diff:.3g}') + print(f' Time : {timings_list[1]:.3g}') + print(f' Niter : {info2["niter"]}') + print(f' Convergence : {info2["success"]}') + print(f' Res. Norm : {info2["res_norm"]:.3g}') + print() + + # Option3 is takes forever to converge, but delivers an accurate solution + a3 = vector_potential_list[2] + info3 = info_list[2] + diff = b - C @ a3 + l2norm_diff = np.sqrt(M2.dot_inner(diff, diff)) + print(f' ----- Option3: C.T @ M2 @ C + M1 @ G @ M0_inv @ DP0 @ G.T @ M1 -----') + print() + print(f' || B - curl(A) || = {l2norm_diff:.3g}') + print(f' Time : {timings_list[2]:.3g}') + print(f' Niter : {info3["niter"]}') + print(f' Convergence : {info3["success"]}') + print(f' Res. Norm : {info3["res_norm"]:.3g}') + print() + print(f'Note that increasing the problem size, or decreasing the tolerance, renders Option3 unafforadable.') + + return a1 + +#============================================================================== +if __name__ == '__main__': + + ncells = [16, 16, 16] + degree = [3, 3, 3] + periodic = [False, False, False] + + comm = None + backend = PSYDAC_BACKEND_GPYCCEL + + logical_domain = Cube('C', bounds1=(0,1), bounds2=(0,1), bounds3=(0,1)) + + class CollelaMap3D(Mapping): + _expressions = {'x': 'x1 + (a/2)*sin(2.*pi*(x1-0.5))*sin(2.*pi*(x2-0.5))', + 'y': 'x2 + (a/2)*sin(2.*pi*(x1-0.5))*sin(2.*pi*(x2-0.5))', + 'z': 'x3'} + + _ldim = 3 + _pdim = 3 + + mapping = CollelaMap3D('C', a=0.1) + + domain = mapping(logical_domain) + derham = Derham(domain) + + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree) + + G, C, D = derham_h.derivatives(kind='linop') + P0, P1, P2, P3 = derham_h.projectors() + + from sympde.utilities.utils import plot_domain + #plot_domain(domain, draw=True, isolines=True) + + A1 = lambda x, y, z: np.sin(2*np.pi*y) * np.sin(2*np.pi*z) + A2 = lambda x, y, z: np.sin(2*np.pi*z) * np.sin(2*np.pi*x) + A3 = lambda x, y, z: np.sin(2*np.pi*x) * np.sin(2*np.pi*y) + A = (A1, A2, A3) + + a_ex = P1(A).coeffs + # a should already satisfy hom. DBCs (n \times a = 0 on the boundary) + # But we can make sure that this holds exactly by projecting it + _, DP1, _, _ = derham_h.dirichlet_projectors(kind='linop') + a_ex = DP1 @ a_ex + + # Now b belongs to H_0(div0) as required + b = C @ a_ex + + a = compute_vector_potential_3d(b, derham_h) diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 1611786b9..6527bbdde 100644 --- a/psydac/api/feec.py +++ b/psydac/api/feec.py @@ -30,6 +30,7 @@ from psydac.feec.pull_push import pull_3d_hdiv, pull_3d_l2, pull_3d_h1vec from psydac.fem.basic import FemSpace, FemLinearOperator +from psydac.fem.lst_preconditioner import construct_LST_preconditioner from psydac.fem.vector import VectorFemSpace from psydac.fem.projectors import DirichletProjector, MultipatchDirichletProjector @@ -332,6 +333,62 @@ def dirichlet_projectors(self, kind='femlinop'): return self._dirichlet_proj elif kind == 'linop': return tuple(femlinop.linop for femlinop in self._dirichlet_proj) + + #-------------------------------------------------------------------------- + def LST_preconditioners(self, *, M0=None, M1=None, M2=None, M3=None, hom_bc=False): + """ + 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, preconditioners for the modified mass matrices M{i}_0 are being returned. + The preconditioner for the last mass matrix of the sequence remains identical as there are no BCs to take care of. + M{i}_0 is a mass matrix of the form + M{i}_0 = DP @ M{i} @ DP + (I - DP) + where DP and I are the corresponding DirichletProjector and IdentityOperator. + See examples/vector_potential_3d. + + Parameters + ---------- + M0, M1, M2, M3 : psydac.linalg.stencil.StencilMatrix | psydac.linalg.block.BlockLinearOperator | None + H1, Hcurl/Hdiv, L2 (2D) or H1, Hcurl, Hdiv, L2 mass matrices or None. + Returns only preconditioners for passed mass matrices. + + hom_bc : bool + If True, return LST preconditioners for modified M{i}_0 = DP @ M{i} @ DP + (I - DP) mass matrices (i=0,1 (2D), i=0,1,2 (3D)). + The arguments M{i} in that case remain the same (M{i}, not M{i}_0). DP and I are DirichletProjector and IdentityOperator. + Default: False. + + Returns + ------- + tuple + tuple of psydac.linalg.stencil.StencilMatrix and/or psydac.linalg.block.BlockLinearOperator + LST preconditioner(s) for the passed mass matrices. + + 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 + + """ + + domain_h = self.domain_h + Ms = (M0, M1, M2, M3) + + M_pc_arr = [] + + for M, Vh in zip(Ms, self.spaces): + if M is not None: + + M_pc = construct_LST_preconditioner(M, domain_h, Vh, hom_bc=hom_bc) + M_pc_arr.append(M_pc) + + return tuple(M_pc_arr) #-------------------------------------------------------------------------- def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False): diff --git a/psydac/fem/lst_preconditioner.py b/psydac/fem/lst_preconditioner.py new file mode 100644 index 000000000..630ddc4ba --- /dev/null +++ b/psydac/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/psydac/fem/tests/test_dirichlet_projectors.py b/psydac/fem/tests/test_dirichlet_projectors.py index bea727244..0298ad321 100644 --- a/psydac/fem/tests/test_dirichlet_projectors.py +++ b/psydac/fem/tests/test_dirichlet_projectors.py @@ -135,7 +135,7 @@ def test_function_space_boundary_projector(dim): nn = NormalVector('nn') for i in range(dim): - print(f' - Test DBP{i}') + print(f' - Test DP{i}') # The function defined here satisfy the corresponding homogeneous Dirichlet BCs if dim == 1: @@ -343,7 +343,7 @@ def test_discrete_derham_boundary_projector(dim): nn = NormalVector('nn') for i in range(dim): - print(f' - Test DBP{i}') + print(f' - Test DP{i}') u, v = elements_of(derham.spaces[i], names='u, v') @@ -477,7 +477,7 @@ def test_discrete_derham_boundary_projector_multipatch(): nn = NormalVector('nn') for i in range(2): - print(f' - Test DBP{i}') + print(f' - Test DP{i}') u, v = elements_of(derham.spaces[i], names='u, v') diff --git a/psydac/linalg/direct_solvers.py b/psydac/linalg/direct_solvers.py index ef2783162..db553a5b2 100644 --- a/psydac/linalg/direct_solvers.py +++ b/psydac/linalg/direct_solvers.py @@ -5,12 +5,28 @@ #---------------------------------------------------------------------------# 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 psydac.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): @@ -59,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/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index a080067a2..a47f4ae2d 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -13,8 +13,8 @@ from psydac.utilities.utils import is_real from psydac.linalg.utilities import _sym_ortho -from psydac.linalg.basic import (Vector, LinearOperator, - InverseLinearOperator, IdentityOperator, ScaledLinearOperator) +from psydac.linalg.basic import (Vector, LinearOperator, InverseLinearOperator, + IdentityOperator, ScaledLinearOperator) __all__ = ( 'inverse', diff --git a/psydac/linalg/tests/test_kron_direct_solver.py b/psydac/linalg/tests/test_kron_direct_solver.py index d11618312..eb157a3b4 100644 --- a/psydac/linalg/tests/test_kron_direct_solver.py +++ b/psydac/linalg/tests/test_kron_direct_solver.py @@ -8,7 +8,7 @@ import pytest import numpy as np from mpi4py import MPI -from scipy.sparse import csc_matrix, dia_matrix, kron +from scipy.sparse import csc_matrix, kron from scipy.sparse.linalg import splu from sympde.calculus import dot @@ -58,28 +58,6 @@ def kron_solve_seq_ref(Y, A, transposed): X = C_op.solve(Y.flatten()) return X.reshape(Y.shape) -# ... - -# ... convert a 1D stencil matrix to band matrix -def to_bnd(A): - - 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 -# ... - -def matrix_to_bandsolver(A): - A.remove_spurious_entries() - A_bnd, la, ua = to_bnd(A) - return BandedSolver(ua, la, A_bnd) def matrix_to_sparse(A): A.remove_spurious_entries() @@ -248,9 +226,9 @@ def get_M1_block_kron_solver(V1, ncells, degree, periodic): B2_mat = [M0_matrices[0], M1_matrices[1], M0_matrices[2]] B3_mat = [M0_matrices[0], M0_matrices[1], M1_matrices[2]] - B1_solvers = [matrix_to_bandsolver(Ai) for Ai in B1_mat] - B2_solvers = [matrix_to_bandsolver(Ai) for Ai in B2_mat] - B3_solvers = [matrix_to_bandsolver(Ai) for Ai in B3_mat] + B1_solvers = [BandedSolver.from_stencil_mat_1d(Ai) for Ai in B1_mat] + B2_solvers = [BandedSolver.from_stencil_mat_1d(Ai) for Ai in B2_mat] + B3_solvers = [BandedSolver.from_stencil_mat_1d(Ai) for Ai in B3_mat] B1_kron_inv = KroneckerLinearSolver(V1_1, V1_1, B1_solvers) B2_kron_inv = KroneckerLinearSolver(V1_2, V1_2, B2_solvers) @@ -321,7 +299,7 @@ def get_inverse_mass_matrices(derham_h, domain_h): B_mat_V0 = [M0_matrices[0], M0_matrices[1]] - B_solvers_V0 = [matrix_to_bandsolver(Ai) for Ai in B_mat_V0] + B_solvers_V0 = [BandedSolver.from_stencil_mat_1d(Ai) for Ai in B_mat_V0] M0_kron_solver = KroneckerLinearSolver(V0h, V0h, B_solvers_V0) @@ -333,8 +311,8 @@ def get_inverse_mass_matrices(derham_h, domain_h): B1_mat_V1 = [M0_matrices[0], M1_matrices[1]] B2_mat_V1 = [M1_matrices[0], M0_matrices[1]] - B1_solvers_V1 = [matrix_to_bandsolver(Ai) for Ai in B1_mat_V1] - B2_solvers_V1 = [matrix_to_bandsolver(Ai) for Ai in B2_mat_V1] + B1_solvers_V1 = [BandedSolver.from_stencil_mat_1d(Ai) for Ai in B1_mat_V1] + B2_solvers_V1 = [BandedSolver.from_stencil_mat_1d(Ai) for Ai in B2_mat_V1] B1_kron_inv_V1 = KroneckerLinearSolver(V1_1, V1_1, B1_solvers_V1) B2_kron_inv_V1 = KroneckerLinearSolver(V1_2, V1_2, B2_solvers_V1) @@ -346,7 +324,7 @@ def get_inverse_mass_matrices(derham_h, domain_h): B_mat_V2 = [M1_matrices[0], M1_matrices[1]] - B_solvers_V2 = [matrix_to_bandsolver(Ai) for Ai in B_mat_V2] + B_solvers_V2 = [BandedSolver.from_stencil_mat_1d(Ai) for Ai in B_mat_V2] M2_kron_solver = KroneckerLinearSolver(V2h, V2h, B_solvers_V2) @@ -361,7 +339,7 @@ def get_inverse_mass_matrices(derham_h, domain_h): @pytest.mark.parametrize( 'p', [1, 3] ) @pytest.mark.parametrize( 'P', [True, False] ) @pytest.mark.parametrize( 'nrhs', [1, 3] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) @pytest.mark.parametrize( 'transposed', [True, False] ) def test_direct_solvers(dtype, seed, n, p, P, nrhs, direct_solver, transposed): @@ -416,7 +394,7 @@ def test_direct_solvers(dtype, seed, n, p, P, nrhs, direct_solver, transposed): @pytest.mark.parametrize( 'dtype', [float, complex] ) @pytest.mark.parametrize( 'seed', [0, 2] ) @pytest.mark.parametrize( 'params', [([8], [2], [False]), ([8,9], [2,3], [False,True])] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) def test_kron_solver_nompi(seed, params, direct_solver, dtype): compare_solve(seed, None, params[0], params[1], params[2], direct_solver, dtype=dtype, transposed=False, verbose=False) @@ -432,7 +410,7 @@ def test_kron_solver_nompi(seed, params, direct_solver, dtype): @pytest.mark.parametrize( 'n1', [8, 17] ) @pytest.mark.parametrize( 'p1', [1, 2, 3] ) @pytest.mark.parametrize( 'P1', [True, False] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) def test_kron_solver_1d_ser(dtype, seed, n1, p1, P1, direct_solver): compare_solve(seed, MPI.COMM_SELF, [n1], [p1], [P1], direct_solver, dtype=dtype, transposed=False, verbose=False) #=============================================================================== @@ -446,7 +424,7 @@ def test_kron_solver_1d_ser(dtype, seed, n1, p1, P1, direct_solver): @pytest.mark.parametrize( 'p2', [1, 2] ) @pytest.mark.parametrize( 'P1', [True, False] ) @pytest.mark.parametrize( 'P2', [True, False] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) def test_kron_solver_2d_ser(dtype, seed, n1, n2, p1, p2, P1, P2, direct_solver): compare_solve(seed, MPI.COMM_SELF, [n1,n2], [p1,p2], [P1,P2], direct_solver, dtype=dtype, transposed=False, verbose=False) #=============================================================================== @@ -459,7 +437,7 @@ def test_kron_solver_2d_ser(dtype, seed, n1, n2, p1, p2, P1, P2, direct_solver): @pytest.mark.parametrize( 'p2', [1, 2] ) @pytest.mark.parametrize( 'P1', [True, False] ) @pytest.mark.parametrize( 'P2', [True, False] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) def test_kron_solver_2d_transposed_ser(seed, n1, n2, p1, p2, P1, P2, direct_solver, dtype): compare_solve(seed, MPI.COMM_SELF, [n1,n2], [p1,p2], [P1,P2], direct_solver, dtype=dtype, transposed=True, verbose=False) #=============================================================================== @@ -497,7 +475,7 @@ def test_kron_solver_nd_ser(seed, dim, dtype): @pytest.mark.parametrize( 'n1', [8, 17] ) @pytest.mark.parametrize( 'p1', [1, 2, 3] ) @pytest.mark.parametrize( 'P1', [True, False] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) @pytest.mark.parallel def test_kron_solver_1d_par(seed, n1, p1, P1, direct_solver, dtype): # we take n1*p1 here to prevent MPI topology problems @@ -512,7 +490,7 @@ def test_kron_solver_1d_par(seed, n1, p1, P1, direct_solver, dtype): @pytest.mark.parametrize( 'p2', [1, 2] ) @pytest.mark.parametrize( 'P1', [True, False] ) @pytest.mark.parametrize( 'P2', [True, False] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) @pytest.mark.parallel def test_kron_solver_2d_par(seed, n1, n2, p1, p2, P1, P2, direct_solver, dtype): compare_solve(seed, MPI.COMM_WORLD, [n1,n2], [p1,p2], [P1,P2], direct_solver, dtype=dtype, transposed=False, verbose=False) @@ -526,7 +504,7 @@ def test_kron_solver_2d_par(seed, n1, n2, p1, p2, P1, P2, direct_solver, dtype): @pytest.mark.parametrize( 'p2', [1, 2] ) @pytest.mark.parametrize( 'P1', [True, False] ) @pytest.mark.parametrize( 'P2', [True, False] ) -@pytest.mark.parametrize( 'direct_solver', [matrix_to_bandsolver, matrix_to_sparse] ) +@pytest.mark.parametrize( 'direct_solver', [BandedSolver.from_stencil_mat_1d, matrix_to_sparse] ) @pytest.mark.parallel def test_kron_solver_2d_transposed_par(seed, n1, n2, p1, p2, P1, P2, direct_solver, dtype): compare_solve(seed, MPI.COMM_WORLD, [n1,n2], [p1,p2], [P1,P2], direct_solver, dtype=dtype, transposed=True, verbose=False) @@ -775,5 +753,5 @@ def get_A_fun_scalar(n=1, m=1, A0=1e04): if __name__ == '__main__': # showcase testcase - compare_solve(0, MPI.COMM_WORLD, [4,4,5], [1,2,3], [False,True,False], matrix_to_bandsolver, dtype=[float, complex], transposed=False, verbose=True) + compare_solve(0, MPI.COMM_WORLD, [4,4,5], [1,2,3], [False,True,False], BandedSolver.from_stencil_mat_1d, dtype=[float, complex], transposed=False, verbose=True) #compare_solve(0, MPI.COMM_WORLD, [2]*10, [1]*10, [False]*10, matrix_to_sparse, verbose=True) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 912a6aea8..c3c9f1e0b 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -3,12 +3,32 @@ # LICENSE file or go to https://github.com/pyccel/psydac/blob/devel/LICENSE # # for full license details. # #---------------------------------------------------------------------------# +import time import numpy as np import pytest +from mpi4py import MPI -from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.linalg.solvers import inverse -from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector +from sympde.topology import elements_of, Mapping, Derham, Square, Cube +from sympde.calculus import inner +from sympde.expr import integral, BilinearForm + +from psydac.api.discretization import discretize +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL +from psydac.ddm.cart import DomainDecomposition, CartDecomposition +from psydac.fem.lst_preconditioner import construct_LST_preconditioner +from psydac.linalg.basic import IdentityOperator +from psydac.linalg.block import BlockVectorSpace +from psydac.linalg.solvers import inverse +from psydac.linalg.stencil import StencilVectorSpace, StencilMatrix, StencilVector +from psydac.fem.tests.test_dirichlet_projectors import _test_LO_equality_using_rng, Annulus, SquareTorus + + +class SinMapping1D(Mapping): + + _expressions = {'x': 'sin((pi/2)*x1)'} + + _ldim = 1 + _pdim = 1 def define_data_hermitian(n, p, dtype=float): @@ -207,6 +227,189 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): assert errh_norm < tol assert solver == 'pcg' or errc_norm < tol +#=============================================================================== +def test_LST_preconditioner(comm=None): + + ncells_3d = [16, 7, 11] + degree_3d = [1, 4, 2] + periodic_3d = [False, True, False] + + prin = True if ((comm is None) or (comm.rank == 0)) else False + backend = PSYDAC_BACKEND_GPYCCEL + + dimensions = [2, 3] + + maxiter = 20000 + tol = 1e-13 + + if prin: + print() + # Test both in 2D and 3D + for dim in dimensions: + if prin: + print(f' ----- Start {dim}D test -----') + + ncells = ncells_3d [0:2] if dim == 2 else ncells_3d + degree = degree_3d [0:2] if dim == 2 else degree_3d + periodic = periodic_3d[0:2] if dim == 2 else periodic_3d + + if dim == 2: + logical_domain = Square('S', bounds1=(0.5, 1), bounds2=(0, 2*np.pi)) + mapping = Annulus('A') + sequence = ['h1', 'hcurl', 'l2'] + else: + logical_domain = Cube ('C', bounds1=(0.5, 1), bounds2=(0, 2*np.pi), bounds3=(0, 1)) + mapping = SquareTorus('ST') + + domain = mapping(logical_domain) + + derham = Derham(domain, sequence=sequence) if dim == 2 else Derham(domain) + + domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=comm) + derham_h = discretize(derham, domain_h, degree=degree) + + Vs = derham.spaces + Vhs = derham_h.spaces + + d_projectors = derham_h.dirichlet_projectors(kind='linop') + + mass_matrices = [] + mass_0_matrices = [] + + for i, (V, Vh) in enumerate(zip(Vs, Vhs)): + u, v = elements_of(V, names='u, v') + expr = inner(u, v) if isinstance(Vh.coeff_space, BlockVectorSpace) else u*v + a = BilinearForm((u, v), integral(domain, expr)) + ah = discretize(a, domain_h, (Vh, Vh), backend=backend) + M = ah.assemble() + mass_matrices.append(M) + if i < dim: + DP = d_projectors[i] + I = IdentityOperator(Vhs[i].coeff_space) + M_0 = DP @ M @ DP + (I - DP) + mass_0_matrices.append(M_0) + + if dim == 2: + M0, M1, M2 = mass_matrices + else: + M0, M1, M2, M3 = mass_matrices + + if dim == 2: + mass_matrix_preconditioners = derham_h.LST_preconditioners(M0=M0, M1=M1, M2=M2 ) + mass_0_matrix_preconditioners = derham_h.LST_preconditioners(M0=M0, M1=M1, hom_bc=True) + else: + mass_matrix_preconditioners = derham_h.LST_preconditioners(M0=M0, M1=M1, M2=M2, M3=M3, ) + mass_0_matrix_preconditioners = derham_h.LST_preconditioners(M0=M0, M1=M1, M2=M2, hom_bc=True) + + # Prepare testing whether obtaining only a subset of preconditioners works + M1_pc, = derham_h.LST_preconditioners(M1=M1 ) + M0_pc, M2_pc = derham_h.LST_preconditioners(M0=M0, M2=M2) + if dim == 3: + M3_pc, = derham_h.LST_preconditioners(M3=M3 ) + + test_pcs = [M0_pc, M1_pc, M2_pc] + if dim == 3: + test_pcs += [M3_pc] + + # Test whether obtaining only a subset of all possible preconditioners works + for pc, test_pc in zip(mass_matrix_preconditioners, test_pcs): + _test_LO_equality_using_rng(pc, test_pc)#, tol=1e-13) + + if prin: + print(f' Accessing a subset of all possible preconditioners works.') + + rng = np.random.default_rng(42) + + # For comparison and testing: Number of iterations required, not using and using a preconditioner + # More information via " -s" when running the test + # dim 2 dim 3 + # M0 M1 M2 M0_0 M1_0 M0 M1 M2 M3 M0_0 M1_0 M2_0 + true_cg_niter = [[90, 681, 62, 77, 600], [486, 7970, 5292, 147, 356, 5892, 4510]] + true_pcg_niter = [[ 6, 6, 2, 5, 5], [ 6, 7, 6, 2, 5, 5, 5]] + # M{i}_0 matrices preconditioned with a LST preconditioner designed for M{i} instead: + # M0_0 M1_0 M0_0 M1_0 M2_0 + true_pcg_niter2 = [[ 23, 24], [ 367, 2867, 220]] + + mass_matrices += mass_0_matrices + mass_matrix_preconditioners += mass_0_matrix_preconditioners + extended_fem_spaces = Vhs + Vhs[:-1] + + for i, (M, Mpc, Vh) in enumerate(zip(mass_matrices, mass_matrix_preconditioners, extended_fem_spaces)): + + cg = False # Set to True to compare iterations and time with not-preconditioned Conjugate Gradient solver + + # hom_bc = False for M0 M1 M2 (M3), then hom_bc = True for M0_0 M1_0 (M2_0) + hom_bc = True if i > dim else False + + # In order to obtain an LST for M{i}_0, we still have to pass M{i} to `construct_LST_preconditioner`.` + # M2 = M{i} if M = M{i}_0 and hence can be used to obtain the pc for M{i}_0 + M2 = mass_matrices[i-dim-1] if i > dim else M + Mpc2 = construct_LST_preconditioner(M2, domain_h, Vh, hom_bc=hom_bc) + _test_LO_equality_using_rng(Mpc, Mpc2, tol=1e-12) + if prin: + print(' The LST pc obtained using derham_h.LST_preconditioners is the same as the one obtained from construct_LST_preconditioner.') + + if cg: + M_inv_cg = inverse(M, 'cg', maxiter=maxiter, tol=tol) + M_inv_pcg = inverse(M, 'pcg', pc=Mpc, maxiter=maxiter, tol=tol) + + y = M.codomain.zeros() + if isinstance(M.codomain, BlockVectorSpace): + for block in y.blocks: + rng.random(size=block._data.shape, dtype="float64", out=block._data) + else: + rng.random(size=y._data.shape, dtype="float64", out=y._data) + + if (i > dim): + if prin: + print(f' Projecting rhs vector into space of functions satisfying hom. DBCs') + DP = d_projectors[i-(dim+1)] + y = DP @ y + + if cg: + t0 = time.time() + x_cg = M_inv_cg @ y + t1 = time.time() + + y_cg = M @ x_cg + diff_cg = y - y_cg + err_cg = np.sqrt(M.codomain.inner(diff_cg, diff_cg)) + time_cg = t1 - t0 + info_cg = M_inv_cg.get_info() + + t0 = time.time() + x_pcg = M_inv_pcg @ y + t1 = time.time() + + y_pcg = M @ x_pcg + diff_pcg = y - y_pcg + err_pcg = np.sqrt(M.codomain.inner(diff_pcg, diff_pcg)) + time_pcg = t1 - t0 + info_pcg = M_inv_pcg.get_info() + + if dim == 2: + mat_txt = f'M{i}' if i <= 2 else f'M{i-3}_0' + else: + mat_txt = f'M{i}' if i <= 3 else f'M{i-4}_0' + + if prin: + print(f' - {mat_txt} test -') + if cg: + print(f' CG : {info_cg} in {time_cg:.3g}s - err.: {err_cg:.3g}') + print(f' PCG: {info_pcg} in {time_pcg:.3g}s - err.: {err_pcg:.3g}') + + if dim == 2: + assert info_pcg['niter'] == true_pcg_niter[0][i] + else: + assert info_pcg['niter'] == true_pcg_niter[1][i] + + print() + +#=============================================================================== +@pytest.mark.parallel +def test_LST_preconditioner_parallel(): + comm = MPI.COMM_WORLD + test_LST_preconditioner(comm=comm) # =============================================================================== # SCRIPT FUNCTIONALITY