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

Fix error in building PML when split across ranks

浏览代码 
The PML can only be split in two dimensions. I.e. an x PML layer can
only be split in the y and z dimension. The x dimension of the grid in
each rank should be greater than the thickness of the PML.
这个提交包含在：
nmannall
2024-06-25 09:51:47 +01:00
父节点 bc860d7d7c
当前提交 f376d472ab
共有 2 个文件被更改，包括 156 次插入 和 39 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

3
gprMax/grid/mpi_grid.py
查看文件

@@ -58,6 +58,9 @@ class MPIGrid(FDTDGrid):
super().__init__()
self.comm = comm
self.x_comm = comm.Sub([False, True, True])
self.y_comm = comm.Sub([True, False, True])
self.z_comm = comm.Sub([True, True, False])
self.mpi_tasks = np.array(self.comm.dims)

192
gprMax/pml.py
查看文件

@@ -20,6 +20,7 @@ import logging
from importlib import import_module
import numpy as np
from mpi4py import MPI
import gprMax.config as config
@@ -45,7 +46,15 @@ class CFSParameter:
}
scalingdirections = ["forward", "reverse"]
def __init__(self, ID=None, scaling="polynomial", scalingprofile=None, scalingdirection="forward", min=0, max=0):
def __init__(
self,
ID=None,
scaling="polynomial",
scalingprofile=None,
scalingdirection="forward",
min=0,
max=0,
):
"""
Args:
ID: string identifier for CFS parameter, can be: 'alpha', 'kappa' or
@@ -95,7 +104,9 @@ class CFS:
# Calculation of the maximum value of sigma from http://dx.doi.org/10.1109/8.546249
m = CFSParameter.scalingprofiles[self.sigma.scalingprofile]
self.sigma.max = (0.8 * (m + 1)) / (config.sim_config.em_consts["z0"] * d * np.sqrt(er * mr))
self.sigma.max = (0.8 * (m + 1)) / (
config.sim_config.em_consts["z0"] * d * np.sqrt(er * mr)
)
def scaling_polynomial(self, order, Evalues, Hvalues):
"""Applies the polynomial to be used for the scaling profile for
@@ -115,7 +126,9 @@ class CFS:
magnetic PML update.
"""
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]
Hvalues = tmp[1::2]
@@ -237,7 +250,9 @@ class PML:
kappamin = sum(cfs.kappa.min for cfs in self.CFS)
if kappamin < 1:
logger.exception(f"Sum of kappamin value(s) for PML is {kappamin} " "and must be greater than one.")
logger.exception(
f"Sum of kappamin value(s) for PML is {kappamin} " "and must be greater than one."
)
raise ValueError
def initialise_field_arrays(self):
@@ -245,42 +260,54 @@ class PML:
if self.direction[0] == "x":
self.EPhi1 = np.zeros(
(len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx + 1, self.ny, self.nz + 1),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.EPhi2 = np.zeros(
(len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx + 1, self.ny + 1, self.nz),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.HPhi1 = np.zeros(
(len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx, self.ny + 1, self.nz),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.HPhi2 = np.zeros(
(len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx, self.ny, self.nz + 1),
dtype=config.sim_config.dtypes["float_or_double"],
)
elif self.direction[0] == "y":
self.EPhi1 = np.zeros(
(len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx, self.ny + 1, self.nz + 1),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.EPhi2 = np.zeros(
(len(self.CFS), self.nx + 1, self.ny + 1, self.nz), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx + 1, self.ny + 1, self.nz),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.HPhi1 = np.zeros(
(len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx + 1, self.ny, self.nz),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.HPhi2 = np.zeros(
(len(self.CFS), self.nx, self.ny, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx, self.ny, self.nz + 1),
dtype=config.sim_config.dtypes["float_or_double"],
)
elif self.direction[0] == "z":
self.EPhi1 = np.zeros(
(len(self.CFS), self.nx, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx, self.ny + 1, self.nz + 1),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.EPhi2 = np.zeros(
(len(self.CFS), self.nx + 1, self.ny, self.nz + 1), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx + 1, self.ny, self.nz + 1),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.HPhi1 = np.zeros(
(len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx + 1, self.ny, self.nz),
dtype=config.sim_config.dtypes["float_or_double"],
)
self.HPhi2 = np.zeros(
(len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=config.sim_config.dtypes["float_or_double"]
(len(self.CFS), self.nx, self.ny + 1, self.nz),
dtype=config.sim_config.dtypes["float_or_double"],
)
def calculate_update_coeffs(self, er, mr):
@@ -291,14 +318,30 @@ class PML:
mr: float of average permeability of underlying material
"""
self.ERA = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.ERB = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.ERE = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.ERF = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.HRA = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.HRB = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.HRE = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.HRF = np.zeros((len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"])
self.ERA = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
self.ERB = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
self.ERE = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
self.ERF = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
self.HRA = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
self.HRB = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
self.HRE = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
self.HRF = np.zeros(
(len(self.CFS), self.thickness), dtype=config.sim_config.dtypes["float_or_double"]
)
for x, cfs in enumerate(self.CFS):
if not cfs.sigma.max:
@@ -310,20 +353,26 @@ class PML:
# Define different parameters depending on PML formulation
if self.G.pmls["formulation"] == "HORIPML":
# HORIPML electric update coefficients
tmp = (2 * config.sim_config.em_consts["e0"] * Ekappa) + self.G.dt * (Ealpha * Ekappa + Esigma)
tmp = (2 * config.sim_config.em_consts["e0"] * Ekappa) + self.G.dt * (
Ealpha * Ekappa + Esigma
)
self.ERA[x, :] = (2 * config.sim_config.em_consts["e0"] + self.G.dt * Ealpha) / tmp
self.ERB[x, :] = (2 * config.sim_config.em_consts["e0"] * Ekappa) / tmp
self.ERE[x, :] = (
(2 * config.sim_config.em_consts["e0"] * Ekappa) - self.G.dt * (Ealpha * Ekappa + Esigma)
(2 * config.sim_config.em_consts["e0"] * Ekappa)
- self.G.dt * (Ealpha * Ekappa + Esigma)
) / tmp
self.ERF[x, :] = (2 * Esigma * self.G.dt) / (Ekappa * tmp)
# HORIPML magnetic update coefficients
tmp = (2 * config.sim_config.em_consts["e0"] * Hkappa) + self.G.dt * (Halpha * Hkappa + Hsigma)
tmp = (2 * config.sim_config.em_consts["e0"] * Hkappa) + self.G.dt * (
Halpha * Hkappa + Hsigma
)
self.HRA[x, :] = (2 * config.sim_config.em_consts["e0"] + self.G.dt * Halpha) / tmp
self.HRB[x, :] = (2 * config.sim_config.em_consts["e0"] * Hkappa) / tmp
self.HRE[x, :] = (
(2 * config.sim_config.em_consts["e0"] * Hkappa) - self.G.dt * (Halpha * Hkappa + Hsigma)
(2 * config.sim_config.em_consts["e0"] * Hkappa)
- self.G.dt * (Halpha * Hkappa + Hsigma)
) / tmp
self.HRF[x, :] = (2 * Hsigma * self.G.dt) / (Hkappa * tmp)
@@ -332,14 +381,18 @@ class PML:
tmp = 2 * config.sim_config.em_consts["e0"] + self.G.dt * Ealpha
self.ERA[x, :] = Ekappa + (self.G.dt * Esigma) / tmp
self.ERB[x, :] = (2 * config.sim_config.em_consts["e0"]) / tmp
self.ERE[x, :] = ((2 * config.sim_config.em_consts["e0"]) - self.G.dt * Ealpha) / tmp
self.ERE[x, :] = (
(2 * config.sim_config.em_consts["e0"]) - self.G.dt * Ealpha
) / tmp
self.ERF[x, :] = (2 * Esigma * self.G.dt) / tmp
# MRIPML magnetic update coefficients
tmp = 2 * config.sim_config.em_consts["e0"] + self.G.dt * Halpha
self.HRA[x, :] = Hkappa + (self.G.dt * Hsigma) / tmp
self.HRB[x, :] = (2 * config.sim_config.em_consts["e0"]) / tmp
self.HRE[x, :] = ((2 * config.sim_config.em_consts["e0"]) - self.G.dt * Halpha) / tmp
self.HRE[x, :] = (
(2 * config.sim_config.em_consts["e0"]) - self.G.dt * Halpha
) / tmp
self.HRF[x, :] = (2 * Hsigma * self.G.dt) / tmp
def update_electric(self):
@@ -348,7 +401,9 @@ class PML:
"""
pmlmodule = "gprMax.cython.pml_updates_electric_" + self.G.pmls["formulation"]
func = getattr(import_module(pmlmodule), "order" + str(len(self.CFS)) + "_" + self.direction)
func = getattr(
import_module(pmlmodule), "order" + str(len(self.CFS)) + "_" + self.direction
)
func(
self.xs,
self.xf,
@@ -380,7 +435,9 @@ class PML:
"""
pmlmodule = "gprMax.cython.pml_updates_magnetic_" + self.G.pmls["formulation"]
func = getattr(import_module(pmlmodule), "order" + str(len(self.CFS)) + "_" + self.direction)
func = getattr(
import_module(pmlmodule), "order" + str(len(self.CFS)) + "_" + self.direction
)
func(
self.xs,
self.xf,
@@ -439,7 +496,11 @@ class CUDAPML(PML):
self.bpg = (
int(
np.ceil(
((self.EPhi1_dev.shape[1] + 1) * (self.EPhi1_dev.shape[2] + 1) * (self.EPhi1_dev.shape[3] + 1))
(
(self.EPhi1_dev.shape[1] + 1)
* (self.EPhi1_dev.shape[2] + 1)
* (self.EPhi1_dev.shape[3] + 1)
)
/ self.G.tpb[0]
)
),
@@ -673,17 +734,53 @@ def build_pml(G, pml_ID, thickness):
pml_type = OpenCLPML
if pml_ID == "x0":
pml = pml_type(G, ID=pml_ID, direction="xminus", xs=0, xf=thickness, ys=0, yf=G.ny, zs=0, zf=G.nz)
pml = pml_type(
G, ID=pml_ID, direction="xminus", xs=0, xf=thickness, ys=0, yf=G.ny, zs=0, zf=G.nz
)
elif pml_ID == "xmax":
pml = pml_type(G, ID=pml_ID, direction="xplus", xs=G.nx - thickness, xf=G.nx, ys=0, yf=G.ny, zs=0, zf=G.nz)
pml = pml_type(
G,
ID=pml_ID,
direction="xplus",
xs=G.nx - thickness,
xf=G.nx,
ys=0,
yf=G.ny,
zs=0,
zf=G.nz,
)
elif pml_ID == "y0":
pml = pml_type(G, ID=pml_ID, direction="yminus", xs=0, xf=G.nx, ys=0, yf=thickness, zs=0, zf=G.nz)
pml = pml_type(
G, ID=pml_ID, direction="yminus", xs=0, xf=G.nx, ys=0, yf=thickness, zs=0, zf=G.nz
)
elif pml_ID == "ymax":
pml = pml_type(G, ID=pml_ID, direction="yplus", xs=0, xf=G.nx, ys=G.ny - thickness, yf=G.ny, zs=0, zf=G.nz)
pml = pml_type(
G,
ID=pml_ID,
direction="yplus",
xs=0,
xf=G.nx,
ys=G.ny - thickness,
yf=G.ny,
zs=0,
zf=G.nz,
)
elif pml_ID == "z0":
pml = pml_type(G, ID=pml_ID, direction="zminus", xs=0, xf=G.nx, ys=0, yf=G.ny, zs=0, zf=thickness)
pml = pml_type(
G, ID=pml_ID, direction="zminus", xs=0, xf=G.nx, ys=0, yf=G.ny, zs=0, zf=thickness
)
elif pml_ID == "zmax":
pml = pml_type(G, ID=pml_ID, direction="zplus", xs=0, xf=G.nx, ys=0, yf=G.ny, zs=G.nz - thickness, zf=G.nz)
pml = pml_type(
G,
ID=pml_ID,
direction="zplus",
xs=0,
xf=G.nx,
ys=0,
yf=G.ny,
zs=G.nz - thickness,
zf=G.nz,
)
if pml_ID[0] == "x":
averageer, averagemr = pml_average_er_mr(
@@ -698,5 +795,22 @@ def build_pml(G, pml_ID, thickness):
G.nx, G.ny, config.get_model_config().ompthreads, G.solid[:, :, pml.zs], ers, mrs
)
# TODO: Use Grid type (currently will cause a circular dependency)
if config.sim_config.mpi:
if pml_ID[0] == "x":
comm = G.x_comm
elif pml_ID[0] == "y":
comm = G.y_comm
elif pml_ID[0] == "z":
comm = G.z_comm
comm.barrier()
print(f"[Rank {G.rank}] pml_ID: {pml_ID}, er: {averageer}, mr: {averagemr}")
averageer = comm.allreduce(averageer, MPI.SUM) / comm.size
averagemr = comm.allreduce(averagemr, MPI.SUM) / comm.size
print(f"[Rank {G.rank}] NEW pml_ID: {pml_ID}, er: {averageer}, mr {averagemr}")
else:
print(f"pml_ID: {pml_ID}, er: {averageer}, mr {averagemr}")
pml.calculate_update_coeffs(averageer, averagemr)
G.pmls["slabs"].append(pml)