你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 04:56:51 +08:00
hdf5 xdmf support
这个提交包含在:
@@ -19,13 +19,14 @@
|
|||||||
import sys
|
import sys
|
||||||
import numpy as np
|
import numpy as np
|
||||||
from struct import pack
|
from struct import pack
|
||||||
|
from gprMax.xdmf import write_output_file
|
||||||
|
|
||||||
from gprMax.utilities import round_value
|
from gprMax.utilities import round_value
|
||||||
|
|
||||||
|
|
||||||
class GeometryView:
|
class GeometryView:
|
||||||
"""Views of the geometry of the model."""
|
"""Views of the geometry of the model."""
|
||||||
|
|
||||||
if sys.byteorder == 'little':
|
if sys.byteorder == 'little':
|
||||||
byteorder = 'LittleEndian'
|
byteorder = 'LittleEndian'
|
||||||
else:
|
else:
|
||||||
@@ -39,7 +40,7 @@ class GeometryView:
|
|||||||
filename (str): Filename to save to.
|
filename (str): Filename to save to.
|
||||||
type (str): Either 'n' for a per cell geometry view, or 'f' for a per cell edge geometry view.
|
type (str): Either 'n' for a per cell geometry view, or 'f' for a per cell edge geometry view.
|
||||||
"""
|
"""
|
||||||
|
|
||||||
self.xs = xs
|
self.xs = xs
|
||||||
self.ys = ys
|
self.ys = ys
|
||||||
self.zs = zs
|
self.zs = zs
|
||||||
@@ -55,26 +56,30 @@ class GeometryView:
|
|||||||
self.filename = filename
|
self.filename = filename
|
||||||
self.type = type
|
self.type = type
|
||||||
|
|
||||||
|
def write_xdmf(self, modelrun, numbermodelruns, G):
|
||||||
|
filename = self.filename[:-4]
|
||||||
|
write_output_file(filename, G)
|
||||||
|
|
||||||
def write_vtk(self, modelrun, numbermodelruns, G):
|
def write_vtk(self, modelrun, numbermodelruns, G):
|
||||||
"""Writes the geometry information to a VTK file. Either ImageData (.vti) for a per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view.
|
"""Writes the geometry information to a VTK file. Either ImageData (.vti) for a per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view.
|
||||||
|
|
||||||
N.B. No Python 3 support for VTK at time of writing (03/2015)
|
N.B. No Python 3 support for VTK at time of writing (03/2015)
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
modelrun (int): Current model run number.
|
modelrun (int): Current model run number.
|
||||||
numbermodelruns (int): Total number of model runs.
|
numbermodelruns (int): Total number of model runs.
|
||||||
G (class): Grid class instance - holds essential parameters describing the model.
|
G (class): Grid class instance - holds essential parameters describing the model.
|
||||||
"""
|
"""
|
||||||
|
|
||||||
# Construct filename from user-supplied name and model run number
|
# Construct filename from user-supplied name and model run number
|
||||||
if numbermodelruns == 1:
|
if numbermodelruns == 1:
|
||||||
self.filename = G.inputdirectory + self.filename
|
self.filename = G.inputdirectory + self.filename
|
||||||
else:
|
else:
|
||||||
self.filename = G.inputdirectory + self.filename + str(modelrun)
|
self.filename = G.inputdirectory + self.filename + str(modelrun)
|
||||||
|
|
||||||
if self.type == 'n':
|
if self.type == 'n':
|
||||||
self.filename += '.vti'
|
self.filename += '.vti'
|
||||||
|
|
||||||
# Calculate number of cells according to requested sampling for geometry view
|
# Calculate number of cells according to requested sampling for geometry view
|
||||||
self.vtk_xscells = round_value(self.xs / self.dx)
|
self.vtk_xscells = round_value(self.xs / self.dx)
|
||||||
self.vtk_xfcells = round_value(self.xf / self.dx)
|
self.vtk_xfcells = round_value(self.xf / self.dx)
|
||||||
@@ -85,16 +90,16 @@ class GeometryView:
|
|||||||
self.vtk_nxcells = self.vtk_xfcells - self.vtk_xscells
|
self.vtk_nxcells = self.vtk_xfcells - self.vtk_xscells
|
||||||
self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells
|
self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells
|
||||||
self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells
|
self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells
|
||||||
|
|
||||||
# Create an array and add numeric IDs for PML, sources and receivers
|
# Create an array and add numeric IDs for PML, sources and receivers
|
||||||
self.srcs_rxs_pml = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8)
|
self.srcs_rxs_pml = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8)
|
||||||
for pml in G.pmls:
|
for pml in G.pmls:
|
||||||
self.srcs_rxs_pml[pml.xs:pml.xf, pml.ys:pml.yf, pml.zs:pml.zf] = 1
|
self.srcs_rxs_pml[pml.xs:pml.xf, pml.ys:pml.yf, pml.zs:pml.zf] = 1
|
||||||
for index, srcrx in enumerate(G.rxs + G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines):
|
for index, srcrx in enumerate(G.rxs + G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines):
|
||||||
self.srcs_rxs_pml[srcrx.xcoord, srcrx.ycoord, srcrx.zcoord] = index + 2
|
self.srcs_rxs_pml[srcrx.xcoord, srcrx.ycoord, srcrx.zcoord] = index + 2
|
||||||
|
|
||||||
vtk_srcs_rxs_pml_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize)
|
vtk_srcs_rxs_pml_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize)
|
||||||
|
|
||||||
with open(self.filename, 'wb') as f:
|
with open(self.filename, 'wb') as f:
|
||||||
f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
|
f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
|
||||||
f.write('<VTKFile type="ImageData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
|
f.write('<VTKFile type="ImageData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
|
||||||
@@ -105,7 +110,7 @@ class GeometryView:
|
|||||||
f.write('<DataArray type="Int8" Name="Sources_Receivers_PML" format="appended" offset="{}" />\n'.format(vtk_srcs_rxs_pml_offset).encode('utf-8'))
|
f.write('<DataArray type="Int8" Name="Sources_Receivers_PML" format="appended" offset="{}" />\n'.format(vtk_srcs_rxs_pml_offset).encode('utf-8'))
|
||||||
f.write('</CellData>\n'.encode('utf-8'))
|
f.write('</CellData>\n'.encode('utf-8'))
|
||||||
f.write('</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
f.write('</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
||||||
|
|
||||||
# Write material IDs
|
# Write material IDs
|
||||||
datasize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
datasize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||||
# Write number of bytes of appended data as UInt32
|
# Write number of bytes of appended data as UInt32
|
||||||
@@ -114,7 +119,7 @@ class GeometryView:
|
|||||||
for j in range(self.ys, self.yf, self.dy):
|
for j in range(self.ys, self.yf, self.dy):
|
||||||
for i in range(self.xs, self.xf, self.dx):
|
for i in range(self.xs, self.xf, self.dx):
|
||||||
f.write(pack('I', G.solid[i, j, k]))
|
f.write(pack('I', G.solid[i, j, k]))
|
||||||
|
|
||||||
# Write source/receiver IDs
|
# Write source/receiver IDs
|
||||||
datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
|
||||||
f.write(pack('I', datasize))
|
f.write(pack('I', datasize))
|
||||||
@@ -124,34 +129,34 @@ class GeometryView:
|
|||||||
f.write(pack('b', self.srcs_rxs_pml[i, j, k]))
|
f.write(pack('b', self.srcs_rxs_pml[i, j, k]))
|
||||||
|
|
||||||
f.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))
|
f.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))
|
||||||
|
|
||||||
self.write_gprmax_info(f, G)
|
self.write_gprmax_info(f, G)
|
||||||
|
|
||||||
elif self.type == 'f':
|
elif self.type == 'f':
|
||||||
self.filename += '.vtp'
|
self.filename += '.vtp'
|
||||||
|
|
||||||
vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
|
vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
|
||||||
vtk_numpoint_components = 3
|
vtk_numpoint_components = 3
|
||||||
vtk_numlines = 2 * self.nx * self.ny + 2 * self.ny * self.nz + 2 * self.nx * self.nz + 3 * self.nx * self.ny * self.nz + self.nx + self.ny + self.nz
|
vtk_numlines = 2 * self.nx * self.ny + 2 * self.ny * self.nz + 2 * self.nx * self.nz + 3 * self.nx * self.ny * self.nz + self.nx + self.ny + self.nz
|
||||||
vtk_numline_components = 2;
|
vtk_numline_components = 2;
|
||||||
|
|
||||||
vtk_connectivity_offset = round_value((vtk_numpoints * vtk_numpoint_components * np.dtype(np.float32).itemsize) + np.dtype(np.uint32).itemsize)
|
vtk_connectivity_offset = round_value((vtk_numpoints * vtk_numpoint_components * np.dtype(np.float32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||||
vtk_offsets_offset = round_value(vtk_connectivity_offset + (vtk_numlines * vtk_numline_components * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)
|
vtk_offsets_offset = round_value(vtk_connectivity_offset + (vtk_numlines * vtk_numline_components * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||||
vtk_materials_offset = round_value(vtk_offsets_offset + (vtk_numlines * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)
|
vtk_materials_offset = round_value(vtk_offsets_offset + (vtk_numlines * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||||
|
|
||||||
with open(self.filename, 'wb') as f:
|
with open(self.filename, 'wb') as f:
|
||||||
f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
|
f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
|
||||||
f.write('<VTKFile type="PolyData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
|
f.write('<VTKFile type="PolyData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
|
||||||
f.write('<PolyData>\n<Piece NumberOfPoints="{}" NumberOfVerts="0" NumberOfLines="{}" NumberOfStrips="0" NumberOfPolys="0">\n'.format(vtk_numpoints, vtk_numlines).encode('utf-8'))
|
f.write('<PolyData>\n<Piece NumberOfPoints="{}" NumberOfVerts="0" NumberOfLines="{}" NumberOfStrips="0" NumberOfPolys="0">\n'.format(vtk_numpoints, vtk_numlines).encode('utf-8'))
|
||||||
|
|
||||||
f.write('<Points>\n<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="0" />\n</Points>\n'.encode('utf-8'))
|
f.write('<Points>\n<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="0" />\n</Points>\n'.encode('utf-8'))
|
||||||
f.write('<Lines>\n<DataArray type="UInt32" Name="connectivity" format="appended" offset="{}" />\n'.format(vtk_connectivity_offset).encode('utf-8'))
|
f.write('<Lines>\n<DataArray type="UInt32" Name="connectivity" format="appended" offset="{}" />\n'.format(vtk_connectivity_offset).encode('utf-8'))
|
||||||
f.write('<DataArray type="UInt32" Name="offsets" format="appended" offset="{}" />\n</Lines>\n'.format(vtk_offsets_offset).encode('utf-8'))
|
f.write('<DataArray type="UInt32" Name="offsets" format="appended" offset="{}" />\n</Lines>\n'.format(vtk_offsets_offset).encode('utf-8'))
|
||||||
|
|
||||||
f.write('<CellData Scalars="Material">\n'.encode('utf-8'))
|
f.write('<CellData Scalars="Material">\n'.encode('utf-8'))
|
||||||
f.write('<DataArray type="UInt32" Name="Material" format="appended" offset="{}" />\n'.format(vtk_materials_offset).encode('utf-8'))
|
f.write('<DataArray type="UInt32" Name="Material" format="appended" offset="{}" />\n'.format(vtk_materials_offset).encode('utf-8'))
|
||||||
f.write('</CellData>\n'.encode('utf-8'))
|
f.write('</CellData>\n'.encode('utf-8'))
|
||||||
|
|
||||||
f.write('</Piece>\n</PolyData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
f.write('</Piece>\n</PolyData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
|
||||||
|
|
||||||
# Write points
|
# Write points
|
||||||
@@ -209,7 +214,7 @@ class GeometryView:
|
|||||||
for j in range(self.ys, self.yf + 1):
|
for j in range(self.ys, self.yf + 1):
|
||||||
for k in range(self.zs, self.zf + 1):
|
for k in range(self.zs, self.zf + 1):
|
||||||
f.write(pack('I', G.ID[0, i, j, k]))
|
f.write(pack('I', G.ID[0, i, j, k]))
|
||||||
|
|
||||||
for i in range(self.xs, self.xf + 1):
|
for i in range(self.xs, self.xf + 1):
|
||||||
for j in range(self.ys, self.yf):
|
for j in range(self.ys, self.yf):
|
||||||
for k in range(self.zs, self.zf + 1):
|
for k in range(self.zs, self.zf + 1):
|
||||||
@@ -227,7 +232,7 @@ class GeometryView:
|
|||||||
|
|
||||||
def write_gprmax_info(self, f, G, materialsonly=False):
|
def write_gprmax_info(self, f, G, materialsonly=False):
|
||||||
"""Writes gprMax specific information relating material, source, and receiver names to numeric identifiers.
|
"""Writes gprMax specific information relating material, source, and receiver names to numeric identifiers.
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
f (filehandle): VTK file.
|
f (filehandle): VTK file.
|
||||||
G (class): Grid class instance - holds essential parameters describing the model.
|
G (class): Grid class instance - holds essential parameters describing the model.
|
||||||
|
@@ -42,7 +42,7 @@ from gprMax.yee_cell_build import build_electric_components, build_magnetic_comp
|
|||||||
|
|
||||||
def main():
|
def main():
|
||||||
"""This is the main function for gprMax."""
|
"""This is the main function for gprMax."""
|
||||||
|
|
||||||
# Print gprMax logo, version, and licencing/copyright information
|
# Print gprMax logo, version, and licencing/copyright information
|
||||||
logo(gprMax.__version__ + ' (Bowmore)')
|
logo(gprMax.__version__ + ' (Bowmore)')
|
||||||
|
|
||||||
@@ -59,7 +59,7 @@ def main():
|
|||||||
numbermodelruns = args.n
|
numbermodelruns = args.n
|
||||||
inputdirectory = os.path.dirname(os.path.abspath(args.inputfile)) + os.sep
|
inputdirectory = os.path.dirname(os.path.abspath(args.inputfile)) + os.sep
|
||||||
inputfile = inputdirectory + os.path.basename(args.inputfile)
|
inputfile = inputdirectory + os.path.basename(args.inputfile)
|
||||||
|
|
||||||
# Create a separate namespace that users can access in any Python code blocks in the input file
|
# Create a separate namespace that users can access in any Python code blocks in the input file
|
||||||
usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'inputdirectory': inputdirectory}
|
usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'inputdirectory': inputdirectory}
|
||||||
|
|
||||||
@@ -92,7 +92,7 @@ def main():
|
|||||||
|
|
||||||
def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None):
|
def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None):
|
||||||
"""Run standard simulation - models are run one after another and each model is parallelised with OpenMP
|
"""Run standard simulation - models are run one after another and each model is parallelised with OpenMP
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
args (dict): Namespace with command line arguments
|
args (dict): Namespace with command line arguments
|
||||||
numbermodelruns (int): Total number of model runs.
|
numbermodelruns (int): Total number of model runs.
|
||||||
@@ -100,7 +100,7 @@ def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
|||||||
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
||||||
optparams (dict): Optional argument. For Taguchi optimisation it provides the parameters to optimise and their values.
|
optparams (dict): Optional argument. For Taguchi optimisation it provides the parameters to optimise and their values.
|
||||||
"""
|
"""
|
||||||
|
|
||||||
tsimstart = perf_counter()
|
tsimstart = perf_counter()
|
||||||
for modelrun in range(1, numbermodelruns + 1):
|
for modelrun in range(1, numbermodelruns + 1):
|
||||||
if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace
|
if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace
|
||||||
@@ -117,20 +117,20 @@ def run_std_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
|||||||
|
|
||||||
def run_benchmark_sim(args, inputfile, usernamespace):
|
def run_benchmark_sim(args, inputfile, usernamespace):
|
||||||
"""Run standard simulation in benchmarking mode - models are run one after another and each model is parallelised with OpenMP
|
"""Run standard simulation in benchmarking mode - models are run one after another and each model is parallelised with OpenMP
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
args (dict): Namespace with command line arguments
|
args (dict): Namespace with command line arguments
|
||||||
inputfile (str): Name of the input file to open.
|
inputfile (str): Name of the input file to open.
|
||||||
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
||||||
"""
|
"""
|
||||||
|
|
||||||
# Number of threads to test - start from max physical CPU cores and divide in half until 1
|
# Number of threads to test - start from max physical CPU cores and divide in half until 1
|
||||||
thread = psutil.cpu_count(logical=False)
|
thread = psutil.cpu_count(logical=False)
|
||||||
threads = [thread]
|
threads = [thread]
|
||||||
while not thread%2:
|
while not thread%2:
|
||||||
thread /= 2
|
thread /= 2
|
||||||
threads.append(int(thread))
|
threads.append(int(thread))
|
||||||
|
|
||||||
benchtimes = np.zeros(len(threads))
|
benchtimes = np.zeros(len(threads))
|
||||||
|
|
||||||
numbermodelruns = len(threads)
|
numbermodelruns = len(threads)
|
||||||
@@ -150,7 +150,7 @@ def run_benchmark_sim(args, inputfile, usernamespace):
|
|||||||
|
|
||||||
def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None):
|
def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None):
|
||||||
"""Run mixed mode MPI/OpenMP simulation - MPI task farm for models with each model parallelised with OpenMP
|
"""Run mixed mode MPI/OpenMP simulation - MPI task farm for models with each model parallelised with OpenMP
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
args (dict): Namespace with command line arguments
|
args (dict): Namespace with command line arguments
|
||||||
numbermodelruns (int): Total number of model runs.
|
numbermodelruns (int): Total number of model runs.
|
||||||
@@ -180,7 +180,7 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
|||||||
data = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, status=status)
|
data = comm.recv(source=MPI.ANY_SOURCE, tag=MPI.ANY_TAG, status=status)
|
||||||
source = status.Get_source()
|
source = status.Get_source()
|
||||||
tag = status.Get_tag()
|
tag = status.Get_tag()
|
||||||
|
|
||||||
if tag == tags.READY.value: # Worker is ready, so send it a task
|
if tag == tags.READY.value: # Worker is ready, so send it a task
|
||||||
if modelrun < numbermodelruns + 1:
|
if modelrun < numbermodelruns + 1:
|
||||||
comm.send(modelrun, dest=source, tag=tags.START.value)
|
comm.send(modelrun, dest=source, tag=tags.START.value)
|
||||||
@@ -188,10 +188,10 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
|||||||
modelrun += 1
|
modelrun += 1
|
||||||
else:
|
else:
|
||||||
comm.send(None, dest=source, tag=tags.EXIT.value)
|
comm.send(None, dest=source, tag=tags.EXIT.value)
|
||||||
|
|
||||||
elif tag == tags.DONE.value:
|
elif tag == tags.DONE.value:
|
||||||
print('Worker {}: completed.'.format(source))
|
print('Worker {}: completed.'.format(source))
|
||||||
|
|
||||||
elif tag == tags.EXIT.value:
|
elif tag == tags.EXIT.value:
|
||||||
print('Worker {}: exited.'.format(source))
|
print('Worker {}: exited.'.format(source))
|
||||||
closedworkers += 1
|
closedworkers += 1
|
||||||
@@ -202,7 +202,7 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
|||||||
comm.send(None, dest=0, tag=tags.READY.value)
|
comm.send(None, dest=0, tag=tags.READY.value)
|
||||||
modelrun = comm.recv(source=0, tag=MPI.ANY_TAG, status=status) # Receive a model number to run from the master
|
modelrun = comm.recv(source=0, tag=MPI.ANY_TAG, status=status) # Receive a model number to run from the master
|
||||||
tag = status.Get_tag()
|
tag = status.Get_tag()
|
||||||
|
|
||||||
# Run a model
|
# Run a model
|
||||||
if tag == tags.START.value:
|
if tag == tags.START.value:
|
||||||
if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace
|
if optparams: # If Taguchi optimistaion, add specific value for each parameter to optimise for each experiment to user accessible namespace
|
||||||
@@ -212,10 +212,10 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
|||||||
modelusernamespace.update({'optparams': tmp})
|
modelusernamespace.update({'optparams': tmp})
|
||||||
else:
|
else:
|
||||||
modelusernamespace = usernamespace
|
modelusernamespace = usernamespace
|
||||||
|
|
||||||
run_model(args, modelrun, numbermodelruns, inputfile, modelusernamespace)
|
run_model(args, modelrun, numbermodelruns, inputfile, modelusernamespace)
|
||||||
comm.send(None, dest=0, tag=tags.DONE.value)
|
comm.send(None, dest=0, tag=tags.DONE.value)
|
||||||
|
|
||||||
elif tag == tags.EXIT.value:
|
elif tag == tags.EXIT.value:
|
||||||
break
|
break
|
||||||
|
|
||||||
@@ -224,34 +224,34 @@ def run_mpi_sim(args, numbermodelruns, inputfile, usernamespace, optparams=None)
|
|||||||
|
|
||||||
def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
||||||
"""Runs a model - processes the input file; builds the Yee cells; calculates update coefficients; runs main FDTD loop.
|
"""Runs a model - processes the input file; builds the Yee cells; calculates update coefficients; runs main FDTD loop.
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
args (dict): Namespace with command line arguments
|
args (dict): Namespace with command line arguments
|
||||||
modelrun (int): Current model run number.
|
modelrun (int): Current model run number.
|
||||||
numbermodelruns (int): Total number of model runs.
|
numbermodelruns (int): Total number of model runs.
|
||||||
inputfile (str): Name of the input file to open.
|
inputfile (str): Name of the input file to open.
|
||||||
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
usernamespace (dict): Namespace that can be accessed by user in any Python code blocks in input file.
|
||||||
|
|
||||||
Returns:
|
Returns:
|
||||||
tsolve (int): Length of time (seconds) of main FDTD calculations
|
tsolve (int): Length of time (seconds) of main FDTD calculations
|
||||||
"""
|
"""
|
||||||
|
|
||||||
# Monitor memory usage
|
# Monitor memory usage
|
||||||
p = psutil.Process()
|
p = psutil.Process()
|
||||||
|
|
||||||
print('\n{}\n\nModel input file: {}\n'.format(68*'*', inputfile))
|
print('\n{}\n\nModel input file: {}\n'.format(68*'*', inputfile))
|
||||||
|
|
||||||
# Add the current model run to namespace that can be accessed by user in any Python code blocks in input file
|
# Add the current model run to namespace that can be accessed by user in any Python code blocks in input file
|
||||||
usernamespace['current_model_run'] = modelrun
|
usernamespace['current_model_run'] = modelrun
|
||||||
print('Constants/variables available for Python scripting: {}\n'.format(usernamespace))
|
print('Constants/variables available for Python scripting: {}\n'.format(usernamespace))
|
||||||
|
|
||||||
# Process any user input Python commands
|
# Process any user input Python commands
|
||||||
processedlines = process_python_include_code(inputfile, usernamespace)
|
processedlines = process_python_include_code(inputfile, usernamespace)
|
||||||
|
|
||||||
# Write a file containing the input commands after Python blocks have been processed
|
# Write a file containing the input commands after Python blocks have been processed
|
||||||
if args.write_processed:
|
if args.write_processed:
|
||||||
write_processed_file(inputfile, modelrun, numbermodelruns, processedlines)
|
write_processed_file(inputfile, modelrun, numbermodelruns, processedlines)
|
||||||
|
|
||||||
# Check validity of command names & that essential commands are present
|
# Check validity of command names & that essential commands are present
|
||||||
singlecmds, multicmds, geometry = check_cmd_names(processedlines)
|
singlecmds, multicmds, geometry = check_cmd_names(processedlines)
|
||||||
|
|
||||||
@@ -284,7 +284,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
|
|
||||||
# Build the PML and calculate initial coefficients
|
# Build the PML and calculate initial coefficients
|
||||||
build_pmls(G)
|
build_pmls(G)
|
||||||
|
|
||||||
# Build the model, i.e. set the material properties (ID) for every edge of every Yee cell
|
# Build the model, i.e. set the material properties (ID) for every edge of every Yee cell
|
||||||
tbuildstart = perf_counter()
|
tbuildstart = perf_counter()
|
||||||
build_electric_components(G.solid, G.rigidE, G.ID, G)
|
build_electric_components(G.solid, G.rigidE, G.ID, G)
|
||||||
@@ -295,7 +295,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
# Process any voltage sources (that have resistance) to create a new material at the source location
|
# Process any voltage sources (that have resistance) to create a new material at the source location
|
||||||
for voltagesource in G.voltagesources:
|
for voltagesource in G.voltagesources:
|
||||||
voltagesource.create_material(G)
|
voltagesource.create_material(G)
|
||||||
|
|
||||||
# Initialise arrays of update coefficients to pass to update functions
|
# Initialise arrays of update coefficients to pass to update functions
|
||||||
G.initialise_std_updatecoeff_arrays()
|
G.initialise_std_updatecoeff_arrays()
|
||||||
|
|
||||||
@@ -309,22 +309,22 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
print('ID\tName\t\tProperties')
|
print('ID\tName\t\tProperties')
|
||||||
print('{}'.format('-'*50))
|
print('{}'.format('-'*50))
|
||||||
for material in G.materials:
|
for material in G.materials:
|
||||||
|
|
||||||
# Calculate update coefficients for material
|
# Calculate update coefficients for material
|
||||||
material.calculate_update_coeffsE(G)
|
material.calculate_update_coeffsE(G)
|
||||||
material.calculate_update_coeffsH(G)
|
material.calculate_update_coeffsH(G)
|
||||||
|
|
||||||
# Store all update coefficients together
|
# Store all update coefficients together
|
||||||
G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
|
G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
|
||||||
G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm
|
G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm
|
||||||
|
|
||||||
# Store coefficients for any dispersive materials
|
# Store coefficients for any dispersive materials
|
||||||
if Material.maxpoles != 0:
|
if Material.maxpoles != 0:
|
||||||
z = 0
|
z = 0
|
||||||
for pole in range(Material.maxpoles):
|
for pole in range(Material.maxpoles):
|
||||||
G.updatecoeffsdispersive[material.numID, z:z+3] = e0 * material.eqt2[pole], material.eqt[pole], material.zt[pole]
|
G.updatecoeffsdispersive[material.numID, z:z+3] = e0 * material.eqt2[pole], material.eqt[pole], material.zt[pole]
|
||||||
z += 3
|
z += 3
|
||||||
|
|
||||||
if G.messages:
|
if G.messages:
|
||||||
if material.deltaer and material.tau:
|
if material.deltaer and material.tau:
|
||||||
tmp = 'delta_epsr={}, tau={} secs; '.format(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer), ', '.join('{:g}'.format(tau) for tau in material.tau))
|
tmp = 'delta_epsr={}, tau={} secs; '.format(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer), ', '.join('{:g}'.format(tau) for tau in material.tau))
|
||||||
@@ -340,7 +340,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
resolution = dispersion_check(G)
|
resolution = dispersion_check(G)
|
||||||
if resolution != 0 and max((G.dx, G.dy, G.dz)) > resolution:
|
if resolution != 0 and max((G.dx, G.dy, G.dz)) > resolution:
|
||||||
print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution should be less than {:g}m'.format(resolution))
|
print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution should be less than {:g}m'.format(resolution))
|
||||||
|
|
||||||
# Write files for any geometry views
|
# Write files for any geometry views
|
||||||
if not G.geometryviews and args.geometry_only:
|
if not G.geometryviews and args.geometry_only:
|
||||||
raise GeneralError('No geometry views found.')
|
raise GeneralError('No geometry views found.')
|
||||||
@@ -348,16 +348,17 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
tgeostart = perf_counter()
|
tgeostart = perf_counter()
|
||||||
for geometryview in G.geometryviews:
|
for geometryview in G.geometryviews:
|
||||||
geometryview.write_vtk(modelrun, numbermodelruns, G)
|
geometryview.write_vtk(modelrun, numbermodelruns, G)
|
||||||
|
geometryview.write_xdmf(modelrun, numbermodelruns, G)
|
||||||
tgeoend = perf_counter()
|
tgeoend = perf_counter()
|
||||||
print('\nGeometry file(s) written in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tgeoend - tgeostart))))
|
print('\nGeometry file(s) written in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tgeoend - tgeostart))))
|
||||||
|
|
||||||
# Run simulation if not doing only geometry
|
# Run simulation if not doing only geometry
|
||||||
if not args.geometry_only:
|
if not args.geometry_only:
|
||||||
|
|
||||||
# Prepare any snapshot files
|
# Prepare any snapshot files
|
||||||
for snapshot in G.snapshots:
|
for snapshot in G.snapshots:
|
||||||
snapshot.prepare_vtk_imagedata(modelrun, numbermodelruns, G)
|
snapshot.prepare_vtk_imagedata(modelrun, numbermodelruns, G)
|
||||||
|
|
||||||
# Adjust position of sources and receivers if required
|
# Adjust position of sources and receivers if required
|
||||||
if G.srcstepx > 0 or G.srcstepy > 0 or G.srcstepz > 0:
|
if G.srcstepx > 0 or G.srcstepy > 0 or G.srcstepz > 0:
|
||||||
for source in itertools.chain(G.hertziandipoles, G.magneticdipoles, G.voltagesources, G.transmissionlines):
|
for source in itertools.chain(G.hertziandipoles, G.magneticdipoles, G.voltagesources, G.transmissionlines):
|
||||||
@@ -396,10 +397,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
for timestep in range(G.iterations):
|
for timestep in range(G.iterations):
|
||||||
if timestep == 0:
|
if timestep == 0:
|
||||||
tstepstart = perf_counter()
|
tstepstart = perf_counter()
|
||||||
|
|
||||||
# Write field outputs to file
|
# Write field outputs to file
|
||||||
write_hdf5(f, timestep, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
|
write_hdf5(f, timestep, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
|
||||||
|
|
||||||
# Write any snapshots to file
|
# Write any snapshots to file
|
||||||
for snapshot in G.snapshots:
|
for snapshot in G.snapshots:
|
||||||
if snapshot.time == timestep + 1:
|
if snapshot.time == timestep + 1:
|
||||||
@@ -432,7 +433,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
|
|
||||||
# Increment absolute time value
|
# Increment absolute time value
|
||||||
abstime += 0.5 * G.dt
|
abstime += 0.5 * G.dt
|
||||||
|
|
||||||
# Update magnetic field components
|
# Update magnetic field components
|
||||||
update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
|
update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
|
||||||
|
|
||||||
@@ -447,7 +448,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
|
|
||||||
# Increment absolute time value
|
# Increment absolute time value
|
||||||
abstime += 0.5 * G.dt
|
abstime += 0.5 * G.dt
|
||||||
|
|
||||||
# Calculate time for two iterations, used to estimate overall runtime
|
# Calculate time for two iterations, used to estimate overall runtime
|
||||||
if timestep == 1:
|
if timestep == 1:
|
||||||
tstepend = perf_counter()
|
tstepend = perf_counter()
|
||||||
@@ -457,7 +458,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
|
|||||||
sys.stdout.flush()
|
sys.stdout.flush()
|
||||||
elif timestep > 1:
|
elif timestep > 1:
|
||||||
update_progress((timestep + 1) / G.iterations)
|
update_progress((timestep + 1) / G.iterations)
|
||||||
|
|
||||||
# Close output file
|
# Close output file
|
||||||
f.close()
|
f.close()
|
||||||
|
|
||||||
|
@@ -22,10 +22,33 @@ import matplotlib.pyplot as plt
|
|||||||
from gprMax.constants import c, floattype, complextype
|
from gprMax.constants import c, floattype, complextype
|
||||||
from gprMax.materials import Material
|
from gprMax.materials import Material
|
||||||
|
|
||||||
|
class Grid():
|
||||||
|
|
||||||
class FDTDGrid:
|
def __init__(self, grid):
|
||||||
|
self.nx = grid.shape[0]
|
||||||
|
self.ny = grid.shape[1]
|
||||||
|
self.nz = grid.shape[2]
|
||||||
|
self.grid = grid
|
||||||
|
|
||||||
|
def n_edges(self):
|
||||||
|
l = self.nx
|
||||||
|
m = self.ny
|
||||||
|
n = self.nz
|
||||||
|
e = (l * m * (n - 1)) + (m * n * (l - 1)) + (l * n * (m - 1))
|
||||||
|
return e
|
||||||
|
|
||||||
|
def n_nodes(self):
|
||||||
|
return self.nx * self.ny * self.nz
|
||||||
|
|
||||||
|
def n_cells(self):
|
||||||
|
return (self.nx - 1) * (self.ny - 1) * (self.nz - 1)
|
||||||
|
|
||||||
|
def get(self, i, j, k):
|
||||||
|
return self.grid[i, j, k]
|
||||||
|
|
||||||
|
class FDTDGrid(Grid):
|
||||||
"""Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters."""
|
"""Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters."""
|
||||||
|
|
||||||
def __init__(self):
|
def __init__(self):
|
||||||
self.inputdirectory = ''
|
self.inputdirectory = ''
|
||||||
self.title = ''
|
self.title = ''
|
||||||
@@ -61,7 +84,7 @@ class FDTDGrid:
|
|||||||
self.rxstepz = 0
|
self.rxstepz = 0
|
||||||
self.rxs = []
|
self.rxs = []
|
||||||
self.snapshots = []
|
self.snapshots = []
|
||||||
|
|
||||||
def initialise_std_arrays(self):
|
def initialise_std_arrays(self):
|
||||||
"""Initialise an array for volumetric material IDs (solid); boolean arrays for specifying whether materials can have dielectric smoothing (rigid);
|
"""Initialise an array for volumetric material IDs (solid); boolean arrays for specifying whether materials can have dielectric smoothing (rigid);
|
||||||
an array for cell edge IDs (ID); and arrays for the electric and magnetic field components. Solid and ID arrays are initialised to free_space (one); rigid arrays
|
an array for cell edge IDs (ID); and arrays for the electric and magnetic field components. Solid and ID arrays are initialised to free_space (one); rigid arrays
|
||||||
@@ -78,7 +101,7 @@ class FDTDGrid:
|
|||||||
self.Hx = np.zeros((self.nx + 1, self.ny, self.nz), dtype=floattype)
|
self.Hx = np.zeros((self.nx + 1, self.ny, self.nz), dtype=floattype)
|
||||||
self.Hy = np.zeros((self.nx, self.ny + 1, self.nz), dtype=floattype)
|
self.Hy = np.zeros((self.nx, self.ny + 1, self.nz), dtype=floattype)
|
||||||
self.Hz = np.zeros((self.nx, self.ny, self.nz + 1), dtype=floattype)
|
self.Hz = np.zeros((self.nx, self.ny, self.nz + 1), dtype=floattype)
|
||||||
|
|
||||||
def initialise_std_updatecoeff_arrays(self):
|
def initialise_std_updatecoeff_arrays(self):
|
||||||
"""Initialise arrays for storing update coefficients."""
|
"""Initialise arrays for storing update coefficients."""
|
||||||
self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=floattype)
|
self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=floattype)
|
||||||
@@ -94,25 +117,25 @@ class FDTDGrid:
|
|||||||
|
|
||||||
def dispersion_check(G):
|
def dispersion_check(G):
|
||||||
"""Check for potential numerical dispersion. Is the smallest wavelength present in the simulation discretised by at least a factor of 10
|
"""Check for potential numerical dispersion. Is the smallest wavelength present in the simulation discretised by at least a factor of 10
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
G (class): Grid class instance - holds essential parameters describing the model.
|
G (class): Grid class instance - holds essential parameters describing the model.
|
||||||
|
|
||||||
Returns:
|
Returns:
|
||||||
resolution (float): Potential numerical dispersion
|
resolution (float): Potential numerical dispersion
|
||||||
"""
|
"""
|
||||||
|
|
||||||
# Minimum number of spatial steps to resolve smallest wavelength
|
# Minimum number of spatial steps to resolve smallest wavelength
|
||||||
resolvedsteps = 10
|
resolvedsteps = 10
|
||||||
|
|
||||||
# Find maximum frequency
|
# Find maximum frequency
|
||||||
maxfreqs = []
|
maxfreqs = []
|
||||||
for waveform in G.waveforms:
|
for waveform in G.waveforms:
|
||||||
|
|
||||||
# User-defined waveform
|
# User-defined waveform
|
||||||
if waveform.uservalues is not None:
|
if waveform.uservalues is not None:
|
||||||
waveformvalues = waveform.uservalues
|
waveformvalues = waveform.uservalues
|
||||||
|
|
||||||
# Built-in waveform
|
# Built-in waveform
|
||||||
else:
|
else:
|
||||||
time = np.linspace(0, 1, G.iterations)
|
time = np.linspace(0, 1, G.iterations)
|
||||||
@@ -130,24 +153,24 @@ def dispersion_check(G):
|
|||||||
|
|
||||||
# Shift powers so that frequency with maximum power is at zero decibels
|
# Shift powers so that frequency with maximum power is at zero decibels
|
||||||
power -= np.amax(power)
|
power -= np.amax(power)
|
||||||
|
|
||||||
# Set maximum frequency to -60dB from maximum power
|
# Set maximum frequency to -60dB from maximum power
|
||||||
freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
|
freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
|
||||||
maxfreqs.append(freqs[freq])
|
maxfreqs.append(freqs[freq])
|
||||||
|
|
||||||
if maxfreqs:
|
if maxfreqs:
|
||||||
maxfreq = max(maxfreqs)
|
maxfreq = max(maxfreqs)
|
||||||
|
|
||||||
# Find minimum wavelength
|
# Find minimum wavelength
|
||||||
ers = [material.er for material in G.materials]
|
ers = [material.er for material in G.materials]
|
||||||
maxer = max(ers)
|
maxer = max(ers)
|
||||||
|
|
||||||
# Minimum velocity
|
# Minimum velocity
|
||||||
minvelocity = c / np.sqrt(maxer)
|
minvelocity = c / np.sqrt(maxer)
|
||||||
|
|
||||||
# Minimum wavelength
|
# Minimum wavelength
|
||||||
minwavelength = minvelocity / maxfreq
|
minwavelength = minvelocity / maxfreq
|
||||||
|
|
||||||
# Resolution of minimum wavelength
|
# Resolution of minimum wavelength
|
||||||
resolution = minwavelength / resolvedsteps
|
resolution = minwavelength / resolvedsteps
|
||||||
|
|
||||||
@@ -159,10 +182,10 @@ def dispersion_check(G):
|
|||||||
|
|
||||||
def get_other_directions(direction):
|
def get_other_directions(direction):
|
||||||
"""Return the two other directions from x, y, z given a single direction
|
"""Return the two other directions from x, y, z given a single direction
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
direction (str): Component x, y or z
|
direction (str): Component x, y or z
|
||||||
|
|
||||||
Returns:
|
Returns:
|
||||||
(tuple): Two directions from x, y, z
|
(tuple): Two directions from x, y, z
|
||||||
"""
|
"""
|
||||||
@@ -174,7 +197,7 @@ def get_other_directions(direction):
|
|||||||
|
|
||||||
def Ix(x, y, z, Hy, Hz, G):
|
def Ix(x, y, z, Hy, Hz, G):
|
||||||
"""Calculates the x-component of current at a grid position.
|
"""Calculates the x-component of current at a grid position.
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
x, y, z (float): Coordinates of position in grid.
|
x, y, z (float): Coordinates of position in grid.
|
||||||
Hy, Hz (memory view): numpy array of magnetic field values.
|
Hy, Hz (memory view): numpy array of magnetic field values.
|
||||||
@@ -184,14 +207,14 @@ def Ix(x, y, z, Hy, Hz, G):
|
|||||||
if y == 0 or z == 0:
|
if y == 0 or z == 0:
|
||||||
Ix = 0
|
Ix = 0
|
||||||
return Ix
|
return Ix
|
||||||
|
|
||||||
else:
|
else:
|
||||||
Ix = G.dy * (Hy[x, y, z - 1] - Hy[x, y, z]) + G.dz * (Hz[x, y, z] - Hz[x, y - 1, z])
|
Ix = G.dy * (Hy[x, y, z - 1] - Hy[x, y, z]) + G.dz * (Hz[x, y, z] - Hz[x, y - 1, z])
|
||||||
return Ix
|
return Ix
|
||||||
|
|
||||||
def Iy(x, y, z, Hx, Hz, G):
|
def Iy(x, y, z, Hx, Hz, G):
|
||||||
"""Calculates the y-component of current at a grid position.
|
"""Calculates the y-component of current at a grid position.
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
x, y, z (float): Coordinates of position in grid.
|
x, y, z (float): Coordinates of position in grid.
|
||||||
Hx, Hz (memory view): numpy array of magnetic field values.
|
Hx, Hz (memory view): numpy array of magnetic field values.
|
||||||
@@ -201,14 +224,14 @@ def Iy(x, y, z, Hx, Hz, G):
|
|||||||
if x == 0 or z == 0:
|
if x == 0 or z == 0:
|
||||||
Iy = 0
|
Iy = 0
|
||||||
return Iy
|
return Iy
|
||||||
|
|
||||||
else:
|
else:
|
||||||
Iy = G.dx * (Hx[x, y, z] - Hx[x, y, z - 1]) + G.dz * (Hz[x - 1, y, z] - Hz[x, y, z])
|
Iy = G.dx * (Hx[x, y, z] - Hx[x, y, z - 1]) + G.dz * (Hz[x - 1, y, z] - Hz[x, y, z])
|
||||||
return Iy
|
return Iy
|
||||||
|
|
||||||
def Iz(x, y, z, Hx, Hy, G):
|
def Iz(x, y, z, Hx, Hy, G):
|
||||||
"""Calculates the z-component of current at a grid position.
|
"""Calculates the z-component of current at a grid position.
|
||||||
|
|
||||||
Args:
|
Args:
|
||||||
x, y, z (float): Coordinates of position in grid.
|
x, y, z (float): Coordinates of position in grid.
|
||||||
Hx, Hy (memory view): numpy array of magnetic field values.
|
Hx, Hy (memory view): numpy array of magnetic field values.
|
||||||
@@ -218,7 +241,7 @@ def Iz(x, y, z, Hx, Hy, G):
|
|||||||
if x == 0 or y == 0:
|
if x == 0 or y == 0:
|
||||||
Iz = 0
|
Iz = 0
|
||||||
return Iz
|
return Iz
|
||||||
|
|
||||||
else:
|
else:
|
||||||
Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z])
|
Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z])
|
||||||
return Iz
|
return Iz
|
||||||
|
361
gprMax/xdmf.py
普通文件
361
gprMax/xdmf.py
普通文件
@@ -0,0 +1,361 @@
|
|||||||
|
import numpy as np
|
||||||
|
import h5py
|
||||||
|
from lxml import etree
|
||||||
|
from gprMax.grid import Grid
|
||||||
|
import copy
|
||||||
|
|
||||||
|
class Edges:
|
||||||
|
|
||||||
|
def __init__(self, grid):
|
||||||
|
|
||||||
|
"""
|
||||||
|
Class to define some connectivity of for an n x l x m
|
||||||
|
grid
|
||||||
|
"""
|
||||||
|
self.total_edges = grid.n_edges()
|
||||||
|
self.edges = np.zeros((self.total_edges, 2), np.float32)
|
||||||
|
self.edge_count = 0
|
||||||
|
self.grid = grid
|
||||||
|
"""
|
||||||
|
Adds the the edge specified by in_node and the i,j,k position of the outnode
|
||||||
|
"""
|
||||||
|
def add_edge(self, in_label, i, j, k):
|
||||||
|
|
||||||
|
out_label = self.grid.get(i, j, k)
|
||||||
|
self.edges[self.edge_count] = np.array([in_label, out_label])
|
||||||
|
self.edge_count += 1
|
||||||
|
|
||||||
|
class Coordinates:
|
||||||
|
|
||||||
|
def __init__(self, grid):
|
||||||
|
self.total_coordinates = grid.nx * grid.ny * grid.nz
|
||||||
|
self.coordinate_count = 0
|
||||||
|
self.coordinates = np.zeros((self.total_coordinates, 3), np.float32)
|
||||||
|
|
||||||
|
def add_coordinate(self, x, y, z):
|
||||||
|
self.coordinates[self.coordinate_count] = np.array([x, y, z])
|
||||||
|
self.coordinate_count += 1
|
||||||
|
|
||||||
|
def hexCellPicker(grid, i, j, k):
|
||||||
|
"""
|
||||||
|
This is the ordering of nodes in the hexahedron cell.
|
||||||
|
|
||||||
|
7 --------- 6
|
||||||
|
/ /|
|
||||||
|
4 --------- 5 2
|
||||||
|
| 3 | /
|
||||||
|
| / |/
|
||||||
|
0 --------- 1
|
||||||
|
|
||||||
|
0 1 2 3 4 5 6 7
|
||||||
|
"""
|
||||||
|
|
||||||
|
cell = [
|
||||||
|
grid[i][j][k],
|
||||||
|
# 1
|
||||||
|
grid[i + 1][j][k],
|
||||||
|
# 2
|
||||||
|
grid[i + 1][j + 1][k],
|
||||||
|
# 3
|
||||||
|
grid[i][j + 1][k],
|
||||||
|
# 4
|
||||||
|
grid[i][j][k + 1],
|
||||||
|
# 5
|
||||||
|
grid[i + 1][j][k + 1],
|
||||||
|
# 6
|
||||||
|
grid[i + 1][j + 1][k + 1],
|
||||||
|
# 7
|
||||||
|
grid[i][j + 1][k + 1]
|
||||||
|
]
|
||||||
|
|
||||||
|
return cell
|
||||||
|
|
||||||
|
class Solids:
|
||||||
|
|
||||||
|
def __init__(self, fdtd_grid):
|
||||||
|
self.count = 0
|
||||||
|
self.fdtd_grid = fdtd_grid
|
||||||
|
self.total_solids = fdtd_grid.n_cells()
|
||||||
|
self.solids = np.zeros((self.total_solids), np.float32)
|
||||||
|
|
||||||
|
def add_solid(self, i, j, k):
|
||||||
|
|
||||||
|
self.solids[self.count] = self.fdtd_grid.solid[i][j][k]
|
||||||
|
self.count += 1
|
||||||
|
|
||||||
|
class SolidLabels():
|
||||||
|
|
||||||
|
def __init__(self, label_grid):
|
||||||
|
self.count = 0
|
||||||
|
self.label_grid = label_grid
|
||||||
|
self.total_solids = label_grid.n_cells()
|
||||||
|
self.solid_labels = np.zeros((self.total_solids, 8), np.float32)
|
||||||
|
|
||||||
|
def add(self, i, j, k):
|
||||||
|
|
||||||
|
solid_labels = hexCellPicker(self.label_grid.grid, i, j, k)
|
||||||
|
self.solid_labels[self.count] = solid_labels
|
||||||
|
self.count += 1
|
||||||
|
|
||||||
|
|
||||||
|
class Materials:
|
||||||
|
|
||||||
|
def __init__(self, fdtd_grid):
|
||||||
|
self.fdtd_grid = fdtd_grid
|
||||||
|
self.n_edges = fdtd_grid.n_edges()
|
||||||
|
self.materials = np.zeros((self.n_edges), np.float32)
|
||||||
|
self.material_count = 0
|
||||||
|
|
||||||
|
# direction x->0 y->1 z->2
|
||||||
|
def add_material(self, i, j, k, direction):
|
||||||
|
|
||||||
|
material = self.fdtd_grid.ID[direction, i, j, k]
|
||||||
|
self.materials[self.material_count] = material
|
||||||
|
|
||||||
|
self.material_count += 1
|
||||||
|
|
||||||
|
def process_grid(fdtd_grid):
|
||||||
|
|
||||||
|
# Dimensions of the problem domain.
|
||||||
|
nx = fdtd_grid.nx
|
||||||
|
ny = fdtd_grid.ny
|
||||||
|
nz = fdtd_grid.nz
|
||||||
|
|
||||||
|
# label each node in the space
|
||||||
|
labels = np.arange(nx * ny * nz).reshape(nx, ny, nz)
|
||||||
|
|
||||||
|
label_grid = Grid(labels)
|
||||||
|
|
||||||
|
# Edges define the connectivity of the grid.
|
||||||
|
edges = Edges(label_grid)
|
||||||
|
|
||||||
|
# Material for each edge
|
||||||
|
edge_materials = Materials(fdtd_grid)
|
||||||
|
|
||||||
|
# Define coordinates for each node
|
||||||
|
coordinates = Coordinates(fdtd_grid)
|
||||||
|
|
||||||
|
# Material for each solid
|
||||||
|
solids = Solids(fdtd_grid)
|
||||||
|
|
||||||
|
# Connectivity for hexhahedron grid
|
||||||
|
solid_labels = SolidLabels(label_grid)
|
||||||
|
|
||||||
|
i_max = nx - 1
|
||||||
|
j_max = ny - 1
|
||||||
|
k_max = nz - 1
|
||||||
|
|
||||||
|
for i, ix in enumerate(labels):
|
||||||
|
for j, jx in enumerate(ix):
|
||||||
|
for k, kx in enumerate(jx):
|
||||||
|
|
||||||
|
label = labels[i][j][k]
|
||||||
|
|
||||||
|
# Each vertex can have varying numbers of edges
|
||||||
|
|
||||||
|
# Type 1 vertex
|
||||||
|
if i < i_max and j < j_max and k < k_max:
|
||||||
|
edges.add_edge(label, i + 1, j, k)
|
||||||
|
edges.add_edge(label, i, j + 1, k)
|
||||||
|
edges.add_edge(label, i, j, k + 1)
|
||||||
|
|
||||||
|
edge_materials.add_material(i, j, k, 0)
|
||||||
|
edge_materials.add_material(i, j, k, 1)
|
||||||
|
edge_materials.add_material(i, j, k, 2)
|
||||||
|
|
||||||
|
# Only this node can support a cell
|
||||||
|
solids.add_solid(i, j, k)
|
||||||
|
solid_labels.add(i, j, k)
|
||||||
|
|
||||||
|
# Type 2 vertex
|
||||||
|
elif i < i_max and j == j_max and k == k_max:
|
||||||
|
edges.add_edge(label, i + 1, j, k)
|
||||||
|
edge_materials.add_material(i, j, k, 0)
|
||||||
|
|
||||||
|
# Type 7 vertex
|
||||||
|
elif i < i_max and j == j_max and k < k_max:
|
||||||
|
edges.add_edge(label, i + 1, j, k)
|
||||||
|
edges.add_edge(label, i, j, k + 1)
|
||||||
|
edge_materials.add_material(i, j, k, 0)
|
||||||
|
edge_materials.add_material(i, j, k, 2)
|
||||||
|
|
||||||
|
# Type 6 vertex
|
||||||
|
elif i == i_max and j == j_max and k < k_max:
|
||||||
|
edges.add_edge(label, i, j, k + 1)
|
||||||
|
edge_materials.add_material(i, j, k, 2)
|
||||||
|
|
||||||
|
# Type 5 vertex
|
||||||
|
elif i == i_max and j < j_max and k < k_max:
|
||||||
|
edges.add_edge(label, i, j, k + 1)
|
||||||
|
edges.add_edge(label, i, j + 1, k)
|
||||||
|
edge_materials.add_material(i, j, k, 2)
|
||||||
|
edge_materials.add_material(i, j, k, 1)
|
||||||
|
|
||||||
|
# Type 4 vertex
|
||||||
|
elif i == i_max and j < j_max and k == k_max:
|
||||||
|
edges.add_edge(label, i, j + 1, k)
|
||||||
|
edge_materials.add_material(i, j, k, 1)
|
||||||
|
|
||||||
|
# Type 8 vertex
|
||||||
|
elif i < i_max and j < j_max and k == k_max:
|
||||||
|
edges.add_edge(label, i, j + 1, k)
|
||||||
|
edges.add_edge(label, i + 1, j, k)
|
||||||
|
edge_materials.add_material(i, j, k, 1)
|
||||||
|
edge_materials.add_material(i, j, k, 0)
|
||||||
|
|
||||||
|
# Type 3 vertex
|
||||||
|
# Has no new connectivity
|
||||||
|
elif i == i_max and j == j_max and k == k_max:
|
||||||
|
pass
|
||||||
|
else:
|
||||||
|
print('oh no')
|
||||||
|
|
||||||
|
# Add the coordinates
|
||||||
|
coordinates.add_coordinate(i, j, k)
|
||||||
|
|
||||||
|
#x = np.arange(fdtd_grid.nx)
|
||||||
|
#y = np.arange(fdtd_grid.ny)
|
||||||
|
#z = np.arange(fdtd_grid.nz)
|
||||||
|
|
||||||
|
return {
|
||||||
|
'coordinates': coordinates,
|
||||||
|
'solids': solids,
|
||||||
|
'solid_labels': solid_labels,
|
||||||
|
'edges': edges,
|
||||||
|
'edge_materials': edge_materials,
|
||||||
|
#'x': x,
|
||||||
|
#'y': y,
|
||||||
|
#'z': z
|
||||||
|
}
|
||||||
|
|
||||||
|
def write_output_file(filename, grid):
|
||||||
|
|
||||||
|
data = process_grid(grid)
|
||||||
|
data['filename'] = filename
|
||||||
|
data['xml_doc'] = create_xdmf_markup(data)
|
||||||
|
|
||||||
|
write_H5file(data)
|
||||||
|
write_xml_doc(data)
|
||||||
|
|
||||||
|
def write_xml_doc(options):
|
||||||
|
#write xml to file
|
||||||
|
with open(options['filename'] + '.xdmf', 'wb') as xdmf_f:
|
||||||
|
xdmf_f.write(options['xml_doc'])
|
||||||
|
|
||||||
|
def write_H5file(options):
|
||||||
|
|
||||||
|
f = h5py.File(options['filename'] + '.h5', "w")
|
||||||
|
coords = f.create_group("mesh")
|
||||||
|
coords.create_dataset('coordinates', data=options['coordinates'].coordinates)
|
||||||
|
coords.create_dataset('connectivity', data=options['edges'].edges)
|
||||||
|
coords.create_dataset('solid_connectivity', data=options['solid_labels'].solid_labels)
|
||||||
|
#coords.create_dataset('x', data=options['x'])
|
||||||
|
#coords.create_dataset('y', data=options['y'])
|
||||||
|
#coords.create_dataset('z', data=options['z'])
|
||||||
|
data = f.create_group("data")
|
||||||
|
data.create_dataset('materials', data=options['edge_materials'].materials)
|
||||||
|
data.create_dataset('solids', data=options['solids'].solids)
|
||||||
|
|
||||||
|
def create_xdmf_markup(options):
|
||||||
|
|
||||||
|
# Write the XDMF markup for edge style grid
|
||||||
|
xdmf_el = etree.Element("Xdmf", Version="2.0")
|
||||||
|
|
||||||
|
domain_el = etree.Element("Domain")
|
||||||
|
xdmf_el.append(domain_el)
|
||||||
|
|
||||||
|
grid_el = etree.Element("Grid", Name="Edges", GridType="Uniform")
|
||||||
|
domain_el.append(grid_el)
|
||||||
|
|
||||||
|
# Create the grid node
|
||||||
|
topology_el = etree.Element("Topology", TopologyType="Polyline", NumberOfElements=str(options['edges'].total_edges))
|
||||||
|
grid_el.append(topology_el)
|
||||||
|
|
||||||
|
topology_dimensions = "{} 2".format(options['edges'].total_edges)
|
||||||
|
top_data_el = etree.Element("DataItem", Dimensions=topology_dimensions, NumberType="Float", Precision="8", Format="HDF")
|
||||||
|
top_data_el.text = "{}:/mesh/connectivity".format(options['filename'] + '.h5')
|
||||||
|
topology_el.append(top_data_el)
|
||||||
|
|
||||||
|
# Create the Geometry node
|
||||||
|
geometry_el = etree.Element("Geometry", GeometryType="XYZ")
|
||||||
|
grid_el.append(geometry_el)
|
||||||
|
|
||||||
|
# Create the origin coordinates
|
||||||
|
coordinates_dimensions = "{} 3".format(options['coordinates'].total_coordinates)
|
||||||
|
origin_el = etree.Element("DataItem", Dimensions=coordinates_dimensions, NumberType="Float", Precision="8", Format="HDF")
|
||||||
|
origin_el.text = "{}:/mesh/coordinates".format(options['filename'] + '.h5')
|
||||||
|
geometry_el.append(origin_el)
|
||||||
|
|
||||||
|
# Create the materials attribute
|
||||||
|
attr_el = etree.Element("Attribute", Center="Cell", Name="Edge_Materials")
|
||||||
|
grid_el.append(attr_el)
|
||||||
|
|
||||||
|
materials_dimensions = "{} 1".format(options['edge_materials'].material_count)
|
||||||
|
materials_el = etree.Element("DataItem", Dimensions=materials_dimensions, NumberType="Float", Precision="8", Format="HDF")
|
||||||
|
materials_el.text = "{}:/data/materials".format(options['filename'] + '.h5')
|
||||||
|
attr_el.append(materials_el)
|
||||||
|
|
||||||
|
"""
|
||||||
|
# VOXEL style markup
|
||||||
|
v_grid_el = etree.Element("Grid", Name="Voxel", GridType="Uniform")
|
||||||
|
domain_el.append(v_grid_el)
|
||||||
|
|
||||||
|
noe = "{} {} {}".format(options['x'].size, options['y'].size, options['y'].size)
|
||||||
|
v_topology_el = etree.Element("Topology", TopologyType="3DRectMesh", NumberOfElements=noe)
|
||||||
|
v_grid_el.append(v_topology_el)
|
||||||
|
|
||||||
|
v_geometry = etree.Element("Geometry", GeometryType="VXVYVZ")
|
||||||
|
v_grid_el.append(v_geometry)
|
||||||
|
|
||||||
|
d1 = etree.Element("DataItem", Dimensions=str(options['x'].size), NumberType="Float", Precision="4", Format="HDF")
|
||||||
|
d1.text = "{}:/mesh/x".format(options['filename'] + '.h5')
|
||||||
|
v_geometry.append(d1)
|
||||||
|
|
||||||
|
d2 = etree.Element("DataItem", Dimensions=str(options['y'].size), NumberType="Float", Precision="4", Format="HDF")
|
||||||
|
d2.text = "{}:/mesh/y".format(options['filename'] + '.h5')
|
||||||
|
v_geometry.append(d2)
|
||||||
|
|
||||||
|
d3 = etree.Element("DataItem", Dimensions=str(options['z'].size), NumberType="Float", Precision="4", Format="HDF")
|
||||||
|
d3.text = "{}:/mesh/z".format(options['filename'] + '.h5')
|
||||||
|
v_geometry.append(d3)
|
||||||
|
|
||||||
|
v_attr = etree.Element("Attribute", Name="material-blocks", Center="Cell")
|
||||||
|
v_grid_el.append(v_attr)
|
||||||
|
|
||||||
|
d4 = etree.Element("DataItem", Format="HDF", NumberType="Float", Precision="4", Dimensions=str(options['solids'].solids.size))
|
||||||
|
d4.text = "{}:/data/solids".format(options['filename'] + '.h5')
|
||||||
|
v_attr.append(d4)
|
||||||
|
|
||||||
|
"""
|
||||||
|
|
||||||
|
v_grid_el = etree.Element("Grid", Name="Voxel", GridType="Uniform")
|
||||||
|
domain_el.append(v_grid_el)
|
||||||
|
|
||||||
|
n_solids = str(options['solids'].solids.size)
|
||||||
|
v_topology_el = etree.Element("Topology", TopologyType="Hexahedron", NumberOfElements=str(options['solids'].solids.size))
|
||||||
|
v_grid_el.append(v_topology_el)
|
||||||
|
|
||||||
|
solid_label_d = "{} {}".format(n_solids, 8)
|
||||||
|
solid_labels_el = etree.Element("DataItem", Dimensions=solid_label_d, Format="HDF")
|
||||||
|
solid_labels_el.text = "{}:/mesh/solid_connectivity".format(options['filename'] + '.h5')
|
||||||
|
v_topology_el.append(solid_labels_el)
|
||||||
|
|
||||||
|
# Same geometry as edges
|
||||||
|
v_grid_el.append(copy.deepcopy(geometry_el))
|
||||||
|
|
||||||
|
v_attr = etree.Element("Attribute", Name="Voxel_Materials", Center="Cell")
|
||||||
|
v_grid_el.append(v_attr)
|
||||||
|
|
||||||
|
d4 = etree.Element("DataItem", Format="HDF", NumberType="Float", Precision="4", Dimensions=str(options['solids'].solids.size))
|
||||||
|
d4.text = "{}:/data/solids".format(options['filename'] + '.h5')
|
||||||
|
v_attr.append(d4)
|
||||||
|
|
||||||
|
# Define a doctype - useful for parsers
|
||||||
|
doc_type = '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>'
|
||||||
|
|
||||||
|
# Serialize elements
|
||||||
|
xml_doc = etree.tostring(xdmf_el, xml_declaration=True,
|
||||||
|
encoding="utf-8", doctype=doc_type, pretty_print=True)
|
||||||
|
|
||||||
|
return xml_doc
|
||||||
|
|
在新工单中引用
屏蔽一个用户