From b662fcd4479d0650009e7bfb3c82df89e2611ed4 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Tue, 28 May 2024 14:01:10 +0100 Subject: [PATCH 1/8] Wip PML Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake/pml.py | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py index 57ca9867..7b1665bc 100644 --- a/ngsPETSc/utils/firedrake/pml.py +++ b/ngsPETSc/utils/firedrake/pml.py @@ -2,16 +2,16 @@ This file contains all the functions related to creating a Perfectly Matched Layer (PML) in Firedrake. ''' - +from collections.abc import Iterable from firedrake import Constant, Function, solve, DirichletBC, TrialFunction, TestFunction -from firedrake import inner, grad, dx, conditional, real +from firedrake import inner, grad, dx, conditional, real, FunctionSpace from firedrake.__future__ import interpolate -from collections.abc import Iterable class PML: """ This class creates a Perfectly Matched Layer (PML) in Firedrake. :arg mesh: The Firedrake mesh. - :arg V: The function space of the Helmholtz problem. + :arg order: The function space order to solve the Poisson problem, when + constructing the PML. :arg pmls: A list of touple of the from: (pml_region_name, pml_inner_boundary_name, pml_outer_boundary_name). :arg k: The wavenumber of the Helmholtz problem, it can be a Constant, a Function, @@ -22,17 +22,17 @@ class PML: coefficient for one of the PML regions. :arg solver_parameters: The solver parameters for the PML problem. """ - def __init__(self, mesh, V, pml_regions, k, alpha=Constant(1j), solver_parameters=None): + def __init__(self, mesh, order, pml_regions, k, alpha=Constant(1j), solver_parameters=None): """ This function initializes the PML object. """ self.mesh = mesh - self.V = V + self.order = order self.pml_regions = pml_regions self.k = k self.alpha = alpha if solver_parameters is None: - solver_parameters = {"ksp_type":"preonly", "pc_type":"lu", + 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.") @@ -42,15 +42,13 @@ def __init__(self, mesh, V, pml_regions, k, alpha=Constant(1j), solver_parameter range(len(mesh.netgen_mesh.GetRegionNames(dim=2))))) if not isinstance(k, (Constant, Function, list)): raise ValueError("The wavenumber must be a Constant, a Function, or a list.") - else: - if not isinstance(k, Iterable): - k = [k]*len(pml_regions) + if not isinstance(k, Iterable): + k = [k]*len(pml_regions) if not isinstance(alpha, (Constant, Function, list)): - raise ValueError("The attenuation coefficient must be a Constant, a Function, or a list.") - else: - if not isinstance(alpha, Iterable): - alpha = [alpha]*len(pml_regions) - + 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.pmls = [] for region in self.pml_regions: @@ -61,6 +59,7 @@ def __init__(self, mesh, V, pml_regions, k, alpha=Constant(1j), solver_parameter 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 + V = FunctionSpace(self.mesh, "CG", self.order) u = TrialFunction(V) v = TestFunction(V) F = inner(grad(u), grad(v))*dx(labels2[region[0]]) @@ -73,6 +72,4 @@ def __init__(self, mesh, V, pml_regions, k, alpha=Constant(1j), solver_parameter solve(F == L, dalet, bcs=bcs, solver_parameters=solver_parameters) sigma = interpolate(1+(alpha/k)*dalet, V) self.pmls = self.pmls + [conditional(real(sigma) < 1e-12, 1, sigma)] - - - + \ No newline at end of file From a6a569527dd91410b1c055d118b89e17ac00271a Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Tue, 28 May 2024 14:09:51 +0100 Subject: [PATCH 2/8] PML weighted measure Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake/pml.py | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py index 7b1665bc..f34b7703 100644 --- a/ngsPETSc/utils/firedrake/pml.py +++ b/ngsPETSc/utils/firedrake/pml.py @@ -5,7 +5,20 @@ from collections.abc import Iterable from firedrake import Constant, Function, solve, DirichletBC, TrialFunction, TestFunction from firedrake import inner, grad, dx, conditional, real, FunctionSpace +from firedrake import ufl from firedrake.__future__ import interpolate + +class WeightedMeasure(ufl.Measure): + """ + This class creates a weighted measure in Firedrake. + """ + def __init__(self, *args, weight=None, **kwargs): + self.weight = weight + ufl.Measure.__init__(self, *args, **kwargs) + + def __rmul__(self, integrand): + return ufl.Measure.__rmul__(self, self.weight * integrand) + class PML: """ This class creates a Perfectly Matched Layer (PML) in Firedrake. @@ -59,17 +72,22 @@ def __init__(self, mesh, order, pml_regions, k, alpha=Constant(1j), solver_param 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 - V = FunctionSpace(self.mesh, "CG", self.order) - u = TrialFunction(V) - v = TestFunction(V) + self.V = FunctionSpace(self.mesh, "CG", self.order) + u = TrialFunction(self.V) + v = TestFunction(self.V) F = inner(grad(u), grad(v))*dx(labels2[region[0]]) for i in range(1, len(self.mesh.netgen_mesh.GetRegionNames(dim=2))+1): if i != labels2[region[0]]: F += inner(grad(u), grad(v))*dx(i) L = inner(Constant(1),v)*dx(2) - bcs = [DirichletBC(V, 1, labels1[region[1]]), DirichletBC(V, 0, labels1[region[2]])] - dalet = Function(V) + bcs = [DirichletBC(self.V, 1, labels1[region[1]]), + DirichletBC(self.V, 0, labels1[region[2]])] + dalet = Function(self.V) solve(F == L, dalet, bcs=bcs, solver_parameters=solver_parameters) - sigma = interpolate(1+(alpha/k)*dalet, V) + sigma = interpolate(1+(alpha/k)*dalet, self.V) self.pmls = self.pmls + [conditional(real(sigma) < 1e-12, 1, sigma)] - \ No newline at end of file + #Assembling the weighted measure + weight = Function(self.V) + for pml in self.pmls: + weight = weight+pml + self.dx = WeightedMeasure("dx", weight=weight) From 61550e6c9ed64081c71c3773e56a674db1852105 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Tue, 28 May 2024 14:41:31 +0100 Subject: [PATCH 3/8] No need to consturct V every pml Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake/pml.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py index f34b7703..6ab1b20d 100644 --- a/ngsPETSc/utils/firedrake/pml.py +++ b/ngsPETSc/utils/firedrake/pml.py @@ -64,6 +64,9 @@ def __init__(self, mesh, order, pml_regions, k, alpha=Constant(1j), solver_param alpha = [alpha]*len(pml_regions) #Construct the PML for each PML region self.pmls = [] + self.V = FunctionSpace(self.mesh, "CG", self.order) + u = TrialFunction(self.V) + v = 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.") @@ -72,9 +75,6 @@ def __init__(self, mesh, order, pml_regions, k, alpha=Constant(1j), solver_param 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 - self.V = FunctionSpace(self.mesh, "CG", self.order) - u = TrialFunction(self.V) - v = TestFunction(self.V) F = inner(grad(u), grad(v))*dx(labels2[region[0]]) for i in range(1, len(self.mesh.netgen_mesh.GetRegionNames(dim=2))+1): if i != labels2[region[0]]: From 924abdc3dcc64126f4a5acc99c84d51dd4917773 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Tue, 28 May 2024 14:42:47 +0100 Subject: [PATCH 4/8] Fix import --- ngsPETSc/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ngsPETSc/__init__.py b/ngsPETSc/__init__.py index 8ecdcd57..97edab7e 100644 --- a/ngsPETSc/__init__.py +++ b/ngsPETSc/__init__.py @@ -16,7 +16,8 @@ if firedrake: from ngsPETSc.utils.firedrake.meshes import * from ngsPETSc.utils.firedrake.hierarchies import * - __all__ = __all__ + ["FiredrakeMesh", "NetgenHierarchy"] + from ngsPETSc.utils.firedrake.pml import * + __all__ = __all__ + ["FiredrakeMesh", "NetgenHierarchy", "PML"] #FEniCSx try: From a0d5fa70e146706a1cf9b155497d34390bfe7c2b Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Tue, 28 May 2024 14:51:22 +0100 Subject: [PATCH 5/8] Import firedrake from the side Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake/pml.py | 43 +++++++++++++++------------------ 1 file changed, 20 insertions(+), 23 deletions(-) diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py index 6ab1b20d..ded7aba3 100644 --- a/ngsPETSc/utils/firedrake/pml.py +++ b/ngsPETSc/utils/firedrake/pml.py @@ -3,21 +3,18 @@ Perfectly Matched Layer (PML) in Firedrake. ''' from collections.abc import Iterable -from firedrake import Constant, Function, solve, DirichletBC, TrialFunction, TestFunction -from firedrake import inner, grad, dx, conditional, real, FunctionSpace -from firedrake import ufl -from firedrake.__future__ import interpolate +import firedrake as fd -class WeightedMeasure(ufl.Measure): +class WeightedMeasure(fd.ufl.Measure): """ This class creates a weighted measure in Firedrake. """ def __init__(self, *args, weight=None, **kwargs): self.weight = weight - ufl.Measure.__init__(self, *args, **kwargs) + fd.ufl.Measure.__init__(self, *args, **kwargs) def __rmul__(self, integrand): - return ufl.Measure.__rmul__(self, self.weight * integrand) + return fd.ufl.Measure.__rmul__(self, self.weight * integrand) class PML: """ @@ -35,7 +32,7 @@ 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, k, alpha=Constant(1j), solver_parameters=None): + def __init__(self, mesh, order, pml_regions, k, alpha=fd.Constant(1j), solver_parameters=None): """ This function initializes the PML object. """ @@ -53,20 +50,20 @@ def __init__(self, mesh, order, pml_regions, k, alpha=Constant(1j), solver_param 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(k, (Constant, Function, list)): + if not isinstance(k, (fd.Constant, fd.Function, list)): raise ValueError("The wavenumber must be a Constant, a Function, or a list.") if not isinstance(k, Iterable): k = [k]*len(pml_regions) - if not isinstance(alpha, (Constant, Function, list)): + 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.pmls = [] - self.V = FunctionSpace(self.mesh, "CG", self.order) - u = TrialFunction(self.V) - v = TestFunction(self.V) + 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.") @@ -75,19 +72,19 @@ def __init__(self, mesh, order, pml_regions, k, alpha=Constant(1j), solver_param 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 = inner(grad(u), grad(v))*dx(labels2[region[0]]) + 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 += inner(grad(u), grad(v))*dx(i) - L = inner(Constant(1),v)*dx(2) - bcs = [DirichletBC(self.V, 1, labels1[region[1]]), - DirichletBC(self.V, 0, labels1[region[2]])] - dalet = Function(self.V) - solve(F == L, dalet, bcs=bcs, solver_parameters=solver_parameters) - sigma = interpolate(1+(alpha/k)*dalet, self.V) - self.pmls = self.pmls + [conditional(real(sigma) < 1e-12, 1, sigma)] + F += fd.inner(fd.grad(u), fd.grad(v))*fd.dx(i) + L = fd.inner(fd.Constant(1),v)*fd.dx(labels2[region[0]]) + bcs = [fd.DirichletBC(self.V, 1, labels1[region[1]]), + fd.DirichletBC(self.V, 0, labels1[region[2]])] + dalet = fd.Function(self.V) + fd.solve(F == L, dalet, bcs=bcs, solver_parameters=solver_parameters) + sigma = fd.__future__.interpolate(1+(alpha/k)*dalet, self.V) + self.pmls = self.pmls + [fd.conditional(fd.real(sigma) < 1e-12, 1, sigma)] #Assembling the weighted measure - weight = Function(self.V) + weight = fd.Function(self.V) for pml in self.pmls: weight = weight+pml self.dx = WeightedMeasure("dx", weight=weight) From 6f3e1ce277c8364a3176bd38cdb69288a8c3171a Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 29 May 2024 01:01:25 +0100 Subject: [PATCH 6/8] Implementing adiabatic layer by Oskooi et al. (2008) Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake/pml.py | 49 +++++++++++++++++++-------------- 1 file changed, 28 insertions(+), 21 deletions(-) diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py index ded7aba3..2efb8a00 100644 --- a/ngsPETSc/utils/firedrake/pml.py +++ b/ngsPETSc/utils/firedrake/pml.py @@ -3,7 +3,13 @@ Perfectly Matched Layer (PML) in Firedrake. ''' from collections.abc import Iterable -import firedrake as fd + +try: + import firedrake as fd +except ImportError: + fd = None + +import numpy as np class WeightedMeasure(fd.ufl.Measure): """ @@ -24,22 +30,18 @@ class PML: constructing the PML. :arg pmls: A list of touple of the from: (pml_region_name, pml_inner_boundary_name, pml_outer_boundary_name). - :arg k: The wavenumber of the Helmholtz problem, it can be a Constant, a Function, - or a list, if it is a list we assume each element of the list the wavenumber - for one of the PML regions. :arg alpha: The attenuation coefficient of the PML, it can be a Constant, a Function, or a list, if it is a list we assume each element of the list the attenuation coefficient for one of the PML regions. :arg solver_parameters: The solver parameters for the PML problem. """ - def __init__(self, mesh, order, pml_regions, k, alpha=fd.Constant(1j), solver_parameters=None): + 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.k = k self.alpha = alpha if solver_parameters is None: solver_parameters = {"ksp_type":"preonly", "pc_type":"lu", @@ -50,17 +52,13 @@ def __init__(self, mesh, order, pml_regions, k, alpha=fd.Constant(1j), solver_pa 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(k, (fd.Constant, fd.Function, list)): - raise ValueError("The wavenumber must be a Constant, a Function, or a list.") - if not isinstance(k, Iterable): - k = [k]*len(pml_regions) 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.pmls = [] + self.dalets = [] self.V = fd.FunctionSpace(self.mesh, "CG", self.order) u = fd.TrialFunction(self.V) v = fd.TestFunction(self.V) @@ -75,16 +73,25 @@ def __init__(self, mesh, order, pml_regions, k, alpha=fd.Constant(1j), solver_pa 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(fd.grad(u), fd.grad(v))*fd.dx(i) + 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, 1, labels1[region[1]]), - fd.DirichletBC(self.V, 0, labels1[region[2]])] + 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) - sigma = fd.__future__.interpolate(1+(alpha/k)*dalet, self.V) - self.pmls = self.pmls + [fd.conditional(fd.real(sigma) < 1e-12, 1, sigma)] - #Assembling the weighted measure - weight = fd.Function(self.V) - for pml in self.pmls: - weight = weight+pml - self.dx = WeightedMeasure("dx", weight=weight) + 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 From 239247a36a65124537e19851dbe527dd1677df57 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 29 May 2024 01:12:36 +0100 Subject: [PATCH 7/8] 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 2efb8a00..c71b8804 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 From 184975b32b442121790896d8e62901362c68f9e1 Mon Sep 17 00:00:00 2001 From: Umberto Zerbinati Date: Wed, 29 May 2024 15:28:44 +0100 Subject: [PATCH 8/8] Using absorption length Signed-off-by: Umberto Zerbinati --- ngsPETSc/utils/firedrake/pml.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ngsPETSc/utils/firedrake/pml.py b/ngsPETSc/utils/firedrake/pml.py index c71b8804..f348280b 100644 --- a/ngsPETSc/utils/firedrake/pml.py +++ b/ngsPETSc/utils/firedrake/pml.py @@ -95,7 +95,8 @@ def adiabatic_layer(self, k, deg_absorb=2): 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.alpha*k*sigma0*(fd.real(d_absorb*dalet)**deg_absorb) + -(sigma0*fd.real(d_absorb*dalet))**2, self.V)) fd.File("absk.pvd").write(absk) return absk