From 239247a36a65124537e19851dbe527dd1677df57 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 29 May 2024 01:12:36 +0100 Subject: [PATCH] Fix CI Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake/pml.py | 122 +++++++++++++++++--------------- 1 file changed, 63 insertions(+), 59 deletions(-) diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py index 2efb8a0..c71b880 100644 --- a/ngsPETSc/utils/firedrake/pml.py +++ b/ngsPETSc/utils/firedrake/pml.py @@ -35,63 +35,67 @@ class PML: coefficient for one of the PML regions. :arg solver_parameters: The solver parameters for the PML problem. """ - def __init__(self, mesh, order, pml_regions, alpha=fd.Constant(1j), solver_parameters=None): - """ - This function initializes the PML object. - """ - self.mesh = mesh - self.order = order - self.pml_regions = pml_regions - self.alpha = alpha - if solver_parameters is None: - solver_parameters = {"ksp_type":"preonly", "pc_type":"lu", - "pc_factor_mat_solver_type":"mumps"} - if not hasattr(self.mesh, "netgen_mesh"): - raise ValueError("The mesh must be a Netgen mesh.") - labels1 = dict(map(lambda i,j : (i,j+1) , mesh.netgen_mesh.GetRegionNames(dim=1), - range(len(mesh.netgen_mesh.GetRegionNames(dim=1))))) - labels2 = dict(map(lambda i,j : (i,j+1) , mesh.netgen_mesh.GetRegionNames(dim=2), - range(len(mesh.netgen_mesh.GetRegionNames(dim=2))))) - if not isinstance(alpha, (fd.Constant, fd.Function, list)): - raise ValueError("The attenuation coefficient must be a \ - Constant, a Function, or a list.") - if not isinstance(alpha, Iterable): - alpha = [alpha]*len(pml_regions) - #Construct the PML for each PML region - self.dalets = [] - self.V = fd.FunctionSpace(self.mesh, "CG", self.order) - u = fd.TrialFunction(self.V) - v = fd.TestFunction(self.V) - for region in self.pml_regions: - if region[0] not in labels2: - raise ValueError("The PML region name must be a valid region name.") - if region[1] not in labels1: - raise ValueError("The PML inner boundary name must be a valid region name.") - if region[2] not in labels1: - raise ValueError("The PML outer boundary name must be a valid region name.") - #Construct the weight function for the PML - F = fd.inner(fd.grad(u), fd.grad(v))*fd.dx(labels2[region[0]]) - for i in range(1, len(self.mesh.netgen_mesh.GetRegionNames(dim=2))+1): - if i != labels2[region[0]]: - F += fd.inner(u, v)*fd.dx(i) - L = fd.inner(fd.Constant(1),v)*fd.dx(labels2[region[0]]) - bcs = [fd.DirichletBC(self.V, 0, labels1[region[1]]), - fd.DirichletBC(self.V, 1, labels1[region[2]])] - dalet = fd.Function(self.V) - fd.solve(F == L, dalet, bcs=bcs, solver_parameters=solver_parameters) - self.dalets = self.dalets + [dalet] + if fd: + if fd.firedrake_configuration.get_config()['options']['complex']: + def __init__(self, mesh, order, pml_regions, alpha=fd.Constant(1j), + solver_parameters=None): + """ + This function initializes the PML object. + """ + self.mesh = mesh + self.order = order + self.pml_regions = pml_regions + self.alpha = alpha + if solver_parameters is None: + solver_parameters = {"ksp_type":"preonly", "pc_type":"lu", + "pc_factor_mat_solver_type":"mumps"} + if not hasattr(self.mesh, "netgen_mesh"): + raise ValueError("The mesh must be a Netgen mesh.") + labels1 = dict(map(lambda i,j : (i,j+1) , mesh.netgen_mesh.GetRegionNames(dim=1), + range(len(mesh.netgen_mesh.GetRegionNames(dim=1))))) + labels2 = dict(map(lambda i,j : (i,j+1) , mesh.netgen_mesh.GetRegionNames(dim=2), + range(len(mesh.netgen_mesh.GetRegionNames(dim=2))))) + if not isinstance(alpha, (fd.Constant, fd.Function, list)): + raise ValueError("The attenuation coefficient must be a \ + Constant, a Function, or a list.") + if not isinstance(alpha, Iterable): + alpha = [alpha]*len(pml_regions) + #Construct the PML for each PML region + self.dalets = [] + self.V = fd.FunctionSpace(self.mesh, "CG", self.order) + u = fd.TrialFunction(self.V) + v = fd.TestFunction(self.V) + for region in self.pml_regions: + if region[0] not in labels2: + raise ValueError("The PML region name must be a valid region name.") + if region[1] not in labels1: + raise ValueError("The PML inner boundary name must be a valid region name.") + if region[2] not in labels1: + raise ValueError("The PML outer boundary name must be a valid region name.") + #Construct the weight function for the PML + F = fd.inner(fd.grad(u), fd.grad(v))*fd.dx(labels2[region[0]]) + for i in range(1, len(self.mesh.netgen_mesh.GetRegionNames(dim=2))+1): + if i != labels2[region[0]]: + F += fd.inner(u, v)*fd.dx(i) + L = fd.inner(fd.Constant(1),v)*fd.dx(labels2[region[0]]) + bcs = [fd.DirichletBC(self.V, 0, labels1[region[1]]), + fd.DirichletBC(self.V, 1, labels1[region[2]])] + dalet = fd.Function(self.V) + fd.solve(F == L, dalet, bcs=bcs, solver_parameters=solver_parameters) + self.dalets = self.dalets + [dalet] - def adiabatic_layer(self, k, deg_absorb=2): - """ - This function returns a complex absorption coefficient, simulating a PML - """ - RT = 1.0e-6 - wave_len = 2*np.pi / k # wavelength - d_absorb = 2 * wave_len # depth of absorber - sigma0 = -(deg_absorb + 1) * np.log(RT) / (2.0 * d_absorb) - absk = fd.assemble(fd.interpolate(fd.Constant(0), self.V)) - for dalet in self.dalets: - absk += fd.assemble(fd.interpolate(self.alpha*k*sigma0*(fd.real(dalet)**deg_absorb) - -(sigma0*fd.real(dalet))**2, self.V)) - fd.File("absk.pvd").write(absk) - return absk + def adiabatic_layer(self, k, deg_absorb=2): + """ + This function returns a complex absorption coefficient, simulating a PML + """ + RT = 1.0e-6 + wave_len = 2*np.pi / k # wavelength + d_absorb = 2 * wave_len # depth of absorber + sigma0 = -(deg_absorb + 1) * np.log(RT) / (2.0 * d_absorb) + absk = fd.assemble(fd.interpolate(fd.Constant(0), self.V)) + for dalet in self.dalets: + absk += fd.assemble(fd.interpolate( + self.alpha*k*sigma0*(fd.real(dalet)**deg_absorb)-(sigma0*fd.real(dalet))**2, + self.V)) + fd.File("absk.pvd").write(absk) + return absk