diff --git a/gprMax/grid/mpi_grid.py b/gprMax/grid/mpi_grid.py index 26662552..3d69f6ed 100644 --- a/gprMax/grid/mpi_grid.py +++ b/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) diff --git a/gprMax/pml.py b/gprMax/pml.py index 83d95ebc..781e4584 100644 --- a/gprMax/pml.py +++ b/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)