publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-08 07:24:19 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

Updated to account for new PML boundary ID naming.

浏览代码 
PML thickness now defined as N + 1 and therefore  '+1's removed from being passed to other functions.
Calculating which spatial step to use in sigmamax function simplified.
这个提交包含在：
Craig Warren
2017-01-30 18:18:45 +00:00
父节点 22fe57e3da
当前提交 6cc1f3c1d3
共有 1 个文件被更改，包括 42 次插入 和 38 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

80
gprMax/pml.py
查看文件

@@ -68,23 +68,16 @@ class CFS(object):
self.kappa = CFSParameter(ID='kappa', scalingprofile='constant', min=1, max=1) self.kappa = CFSParameter(ID='kappa', scalingprofile='constant', min=1, max=1)
self.sigma = CFSParameter(ID='sigma', scalingprofile='quartic', min=0, max=None) self.sigma = CFSParameter(ID='sigma', scalingprofile='quartic', min=0, max=None)
def calculate_sigmamax(self, direction, er, mr, G): def calculate_sigmamax(self, d, er, mr, G):
"""Calculates an optimum value for sigma max based on underlying material properties. """Calculates an optimum value for sigma max based on underlying material properties.
Args: Args:
direction (str): Direction of PML slab d (float): dx, dy, or dz in direction of PML.
er (float): Average permittivity of underlying material. er (float): Average permittivity of underlying material.
mr (float): Average permeability of underlying material. mr (float): Average permeability of underlying material.
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
if direction[0] == 'x':
d = G.dx
elif direction[0] == 'y':
d = G.dy
elif direction[0] == 'z':
d = G.dz
# Calculation of the maximum value of sigma from http://dx.doi.org/10.1109/8.546249 # Calculation of the maximum value of sigma from http://dx.doi.org/10.1109/8.546249
m = CFSParameter.scalingprofiles[self.sigma.scalingprofile] m = CFSParameter.scalingprofiles[self.sigma.scalingprofile]
self.sigma.max = (0.8 * (m + 1)) / (z0 * d * np.sqrt(er * mr)) self.sigma.max = (0.8 * (m + 1)) / (z0 * d * np.sqrt(er * mr))
@@ -105,6 +98,7 @@ class CFS(object):
tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2 * len(Evalues)) / (len(Evalues) - 1)) ** order tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2 * len(Evalues)) / (len(Evalues) - 1)) ** order
Evalues = tmp[0:-1:2] Evalues = tmp[0:-1:2]
Hvalues = tmp[1::2] Hvalues = tmp[1::2]
return Evalues, Hvalues return Evalues, Hvalues
def calculate_values(self, thickness, parameter): def calculate_values(self, thickness, parameter):
@@ -119,8 +113,8 @@ class CFS(object):
Hvalues (float): numpy array holding profile value for magnetic PML update. Hvalues (float): numpy array holding profile value for magnetic PML update.
""" """
Evalues = np.zeros(thickness + 1, dtype=floattype) Evalues = np.zeros(thickness, dtype=floattype)
Hvalues = np.zeros(thickness + 1, dtype=floattype) Hvalues = np.zeros(thickness, dtype=floattype)
if parameter.scalingprofile == 'constant': if parameter.scalingprofile == 'constant':
Evalues += parameter.max Evalues += parameter.max
@@ -147,16 +141,24 @@ class CFS(object):
class PML(object): class PML(object):
"""PML - the implementation comes from the derivation in: http://dx.doi.org/10.1109/TAP.2011.2180344""" """PML - the implementation comes from the derivation in: http://dx.doi.org/10.1109/TAP.2011.2180344"""
slabs = ['xminus', 'yminus', 'zminus', 'xplus', 'yplus', 'zplus'] # IDs for default PML slabs at boundaries of domain
boundaryIDs = ['x0', 'y0', 'z0', 'xmax', 'ymax', 'zmax']
def __init__(self, G, direction=None, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0): # Indicates direction of increasing absorption
# xminus, yminus, zminus - absorption increases in negative direction of x-axis, y-axis, or z-axis
# xplus, yplus, zplus - absorption increases in positive direction of x-axis, y-axis, or z-axis
directions = ['xminus', 'yminus', 'zminus', 'xplus', 'yplus', 'zplus']
def __init__(self, G, ID=None, direction=None, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0):
""" """
Args: Args:
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
direction (str): Identifier for PML slab. ID (str): Identifier for PML slab.
direction (str): Direction of increasing absorption.
xs, xf, ys, yf, zs, zf (float): Extent of the PML slab. xs, xf, ys, yf, zs, zf (float): Extent of the PML slab.
""" """
self.ID = ID
self.direction = direction self.direction = direction
self.xs = xs self.xs = xs
self.xf = xf self.xf = xf
@@ -168,15 +170,17 @@ class PML(object):
self.ny = yf - ys self.ny = yf - ys
self.nz = zf - zs self.nz = zf - zs
# Spatial discretisation and thickness
# (one extra cell of thickness required for interpolation of electric and magnetic scaling values)
if self.direction[0] == 'x': if self.direction[0] == 'x':
self.d = G.dx self.d = G.dx
self.thickness = self.nx self.thickness = self.nx + 1
elif self.direction[0] == 'y': elif self.direction[0] == 'y':
self.d = G.dy self.d = G.dy
self.thickness = self.ny self.thickness = self.ny + 1
elif self.direction[0] == 'z': elif self.direction[0] == 'z':
self.d = G.dz self.d = G.dz
self.thickness = self.nz self.thickness = self.nz + 1
self.CFS = G.cfs self.CFS = G.cfs
if not self.CFS: if not self.CFS:
@@ -212,18 +216,18 @@ class PML(object):
G (class): Grid class instance - holds essential parameters describing the model. G (class): Grid class instance - holds essential parameters describing the model.
""" """
self.ERA = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.ERA = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
self.ERB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.ERB = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
self.ERE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.ERE = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
self.ERF = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.ERF = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
self.HRA = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.HRA = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
self.HRB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.HRB = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
self.HRE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.HRE = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
self.HRF = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype) self.HRF = np.zeros((len(self.CFS), self.thickness), dtype=floattype)
for x, cfs in enumerate(self.CFS): for x, cfs in enumerate(self.CFS):
if not cfs.sigma.max: if not cfs.sigma.max:
cfs.calculate_sigmamax(self.direction, er, mr, G) cfs.calculate_sigmamax(self.d, er, mr, G)
Ealpha, Halpha = cfs.calculate_values(self.thickness, cfs.alpha) Ealpha, Halpha = cfs.calculate_values(self.thickness, cfs.alpha)
Ekappa, Hkappa = cfs.calculate_values(self.thickness, cfs.kappa) Ekappa, Hkappa = cfs.calculate_values(self.thickness, cfs.kappa)
Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma) Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma)
@@ -278,10 +282,10 @@ def build_pmls(G, pbar):
summr = 0 # Sum of relative permeabilities in PML slab summr = 0 # Sum of relative permeabilities in PML slab
if key[0] == 'x': if key[0] == 'x':
if key == 'xminus': if key == 'x0':
pml = PML(G, direction=key, xf=value, yf=G.ny, zf=G.nz) pml = PML(G, ID=key, direction='xminus', xf=value, yf=G.ny, zf=G.nz)
elif key == 'xplus': elif key == 'xmax':
pml = PML(G, direction=key, xs=G.nx - value, xf=G.nx, yf=G.ny, zf=G.nz) pml = PML(G, ID=key, direction='xplus', xs=G.nx - value, xf=G.nx, yf=G.ny, zf=G.nz)
G.pmls.append(pml) G.pmls.append(pml)
for j in range(G.ny): for j in range(G.ny):
for k in range(G.nz): for k in range(G.nz):
@@ -293,10 +297,10 @@ def build_pmls(G, pbar):
averagemr = summr / (G.ny * G.nz) averagemr = summr / (G.ny * G.nz)
elif key[0] == 'y': elif key[0] == 'y':
if key == 'yminus': if key == 'y0':
pml = PML(G, direction=key, yf=value, xf=G.nx, zf=G.nz) pml = PML(G, ID=key, direction='yminus', yf=value, xf=G.nx, zf=G.nz)
elif key == 'yplus': elif key == 'ymax':
pml = PML(G, direction=key, ys=G.ny - value, xf=G.nx, yf=G.ny, zf=G.nz) pml = PML(G, ID=key, direction='yplus', ys=G.ny - value, xf=G.nx, yf=G.ny, zf=G.nz)
G.pmls.append(pml) G.pmls.append(pml)
for i in range(G.nx): for i in range(G.nx):
for k in range(G.nz): for k in range(G.nz):
@@ -308,10 +312,10 @@ def build_pmls(G, pbar):
averagemr = summr / (G.nx * G.nz) averagemr = summr / (G.nx * G.nz)
elif key[0] == 'z': elif key[0] == 'z':
if key == 'zminus': if key == 'z0':
pml = PML(G, direction=key, zf=value, xf=G.nx, yf=G.ny) pml = PML(G, ID=key, direction='zminus', zf=value, xf=G.nx, yf=G.ny)
elif key == 'zplus': elif key == 'zmax':
pml = PML(G, direction=key, zs=G.nz - value, xf=G.nx, yf=G.ny, zf=G.nz) pml = PML(G, ID=key, direction='zplus', zs=G.nz - value, xf=G.nx, yf=G.ny, zf=G.nz)
G.pmls.append(pml) G.pmls.append(pml)
for i in range(G.nx): for i in range(G.nx):
for j in range(G.ny): for j in range(G.ny):