你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 23:14:03 +08:00
255 行
11 KiB
Python
255 行
11 KiB
Python
# Copyright (C) 2015-2022: The University of Edinburgh, United Kingdom
|
|
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
|
#
|
|
# This file is part of gprMax.
|
|
#
|
|
# gprMax is free software: you can redistribute it and/or modify
|
|
# it under the terms of the GNU General Public License as published by
|
|
# the Free Software Foundation, either version 3 of the License, or
|
|
# (at your option) any later version.
|
|
#
|
|
# gprMax is distributed in the hope that it will be useful,
|
|
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
# GNU General Public License for more details.
|
|
#
|
|
# You should have received a copy of the GNU General Public License
|
|
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
import logging
|
|
import sys
|
|
from pathlib import Path
|
|
|
|
import gprMax
|
|
import h5py
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
|
|
from testing.analytical_solutions import hertzian_dipole_fs
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
if sys.platform == 'linux':
|
|
plt.switch_backend('agg')
|
|
|
|
|
|
"""Compare field outputs
|
|
|
|
Usage:
|
|
cd gprMax
|
|
python -m tests.test_models
|
|
"""
|
|
|
|
# Specify directoty with set of models to test
|
|
modelset = 'models_basic'
|
|
# modelset += 'models_advanced'
|
|
# modelset += 'models_pmls'
|
|
|
|
basepath = Path(__file__).parents[0] / modelset
|
|
|
|
|
|
# List of available basic test models
|
|
testmodels = ['hertzian_dipole_fs_analytical', '2D_ExHyHz', '2D_EyHxHz', '2D_EzHxHy',
|
|
'cylinder_Ascan_2D', 'hertzian_dipole_fs', 'hertzian_dipole_hs',
|
|
'hertzian_dipole_dispersive', 'magnetic_dipole_fs']
|
|
|
|
# List of available advanced test models
|
|
# testmodels = ['antenna_GSSI_1500_fs', 'antenna_MALA_1200_fs']
|
|
|
|
# List of available PML models
|
|
# testmodels = ['pml_x0', 'pml_y0', 'pml_z0', 'pml_xmax', 'pml_ymax', 'pml_zmax', 'pml_3D_pec_plate']
|
|
|
|
# Select a specific model if desired
|
|
# testmodels = testmodels[:-1]
|
|
# testmodels = [testmodels[6]]
|
|
testresults = dict.fromkeys(testmodels)
|
|
path = '/rxs/rx1/'
|
|
|
|
# Minimum value of difference to plot (dB)
|
|
plotmin = -160
|
|
|
|
for i, model in enumerate(testmodels):
|
|
|
|
testresults[model] = {}
|
|
|
|
# Run model
|
|
file = basepath / model / model
|
|
gprMax.run(inputfile=file.with_suffix('.in'), gpu=None)
|
|
|
|
# Special case for analytical comparison
|
|
if model == 'hertzian_dipole_fs_analytical':
|
|
# Get output for model file
|
|
filetest = h5py.File(file.with_suffix('.h5'), 'r')
|
|
testresults[model]['Test version'] = filetest.attrs['gprMax']
|
|
|
|
# Get available field output component names
|
|
outputstest = list(filetest[path].keys())
|
|
|
|
# Arrays for storing time
|
|
float_or_double = filetest[path + outputstest[0]].dtype
|
|
timetest = np.linspace(0, (filetest.attrs['Iterations'] - 1) * filetest.attrs['dt'],
|
|
num=filetest.attrs['Iterations']) / 1e-9
|
|
timeref = timetest
|
|
|
|
# Arrays for storing field data
|
|
datatest = np.zeros((filetest.attrs['Iterations'], len(outputstest)),
|
|
dtype=float_or_double)
|
|
for ID, name in enumerate(outputstest):
|
|
datatest[:, ID] = filetest[path + str(name)][:]
|
|
if np.any(np.isnan(datatest[:, ID])):
|
|
logger.exception('Test data contains NaNs')
|
|
raise ValueError
|
|
|
|
# Tx/Rx position to feed to analytical solution
|
|
rxpos = filetest[path].attrs['Position']
|
|
txpos = filetest['/srcs/src1/'].attrs['Position']
|
|
rxposrelative = ((rxpos[0] - txpos[0]),
|
|
(rxpos[1] - txpos[1]),
|
|
(rxpos[2] - txpos[2]))
|
|
|
|
# Analytical solution of a dipole in free space
|
|
dataref = hertzian_dipole_fs(filetest.attrs['Iterations'],
|
|
filetest.attrs['dt'],
|
|
filetest.attrs['dx_dy_dz'], rxposrelative)
|
|
|
|
filetest.close()
|
|
|
|
else:
|
|
# Get output for model and reference files
|
|
fileref = file.stem + '_ref'
|
|
fileref = file.parent / Path(fileref)
|
|
fileref = h5py.File(fileref.with_suffix('.h5'), 'r')
|
|
filetest = h5py.File(file.with_suffix('.h5'), 'r')
|
|
testresults[model]['Ref version'] = fileref.attrs['gprMax']
|
|
testresults[model]['Test version'] = filetest.attrs['gprMax']
|
|
|
|
# Get available field output component names
|
|
outputsref = list(fileref[path].keys())
|
|
outputstest = list(filetest[path].keys())
|
|
if outputsref != outputstest:
|
|
logger.exception('Field output components do not match reference solution')
|
|
raise ValueError
|
|
|
|
# Check that type of float used to store fields matches
|
|
if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype:
|
|
logger.warning(f'Type of floating point number in test model ' +
|
|
f'({filetest[path + outputstest[0]].dtype}) does not ' +
|
|
f'match type in reference solution ({fileref[path + outputsref[0]].dtype})\n')
|
|
float_or_doubleref = fileref[path + outputsref[0]].dtype
|
|
float_or_doubletest = filetest[path + outputstest[0]].dtype
|
|
|
|
# Arrays for storing time
|
|
timeref = np.zeros((fileref.attrs['Iterations']), dtype=float_or_doubleref)
|
|
timeref = np.linspace(0, (fileref.attrs['Iterations'] - 1) * fileref.attrs['dt'],
|
|
num=fileref.attrs['Iterations']) / 1e-9
|
|
timetest = np.zeros((filetest.attrs['Iterations']), dtype=float_or_doubletest)
|
|
timetest = np.linspace(0, (filetest.attrs['Iterations'] - 1) * filetest.attrs['dt'],
|
|
num=filetest.attrs['Iterations']) / 1e-9
|
|
|
|
# Arrays for storing field data
|
|
dataref = np.zeros((fileref.attrs['Iterations'], len(outputsref)),
|
|
dtype=float_or_doubleref)
|
|
datatest = np.zeros((filetest.attrs['Iterations'], len(outputstest)),
|
|
dtype=float_or_doubletest)
|
|
for ID, name in enumerate(outputsref):
|
|
dataref[:, ID] = fileref[path + str(name)][:]
|
|
datatest[:, ID] = filetest[path + str(name)][:]
|
|
if np.any(np.isnan(datatest[:, ID])):
|
|
logger.exception('Test data contains NaNs')
|
|
raise ValueError
|
|
|
|
fileref.close()
|
|
filetest.close()
|
|
|
|
# Diffs
|
|
datadiffs = np.zeros(datatest.shape, dtype=np.float64)
|
|
for i in range(len(outputstest)):
|
|
max = np.amax(np.abs(dataref[:, i]))
|
|
datadiffs[:, i] = np.divide(np.abs(dataref[:, i] - datatest[:, i]), max,
|
|
out=np.zeros_like(dataref[:, i]),
|
|
where=max != 0) # Replace any division by zero with zero
|
|
|
|
# Calculate power (ignore warning from taking a log of any zero values)
|
|
with np.errstate(divide='ignore'):
|
|
datadiffs[:, i] = 20 * np.log10(datadiffs[:, i])
|
|
# Replace any NaNs or Infs from zero division
|
|
datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0
|
|
|
|
# Store max difference
|
|
maxdiff = np.amax(np.amax(datadiffs))
|
|
testresults[model]['Max diff'] = maxdiff
|
|
|
|
# Plot datasets
|
|
fig1, ((ex1, hx1), (ey1, hy1), (ez1, hz1)) = plt.subplots(nrows=3, ncols=2,
|
|
sharex=False, sharey='col',
|
|
subplot_kw=dict(xlabel='Time [ns]'),
|
|
num=model + '.in',
|
|
figsize=(20, 10),
|
|
facecolor='w',
|
|
edgecolor='w')
|
|
ex1.plot(timetest, datatest[:, 0], 'r', lw=2, label=model)
|
|
ex1.plot(timeref, dataref[:, 0], 'g', lw=2, ls='--', label=model + '(Ref)')
|
|
ey1.plot(timetest, datatest[:, 1], 'r', lw=2, label=model)
|
|
ey1.plot(timeref, dataref[:, 1], 'g', lw=2, ls='--', label=model + '(Ref)')
|
|
ez1.plot(timetest, datatest[:, 2], 'r', lw=2, label=model)
|
|
ez1.plot(timeref, dataref[:, 2], 'g', lw=2, ls='--', label=model + '(Ref)')
|
|
hx1.plot(timetest, datatest[:, 3], 'r', lw=2, label=model)
|
|
hx1.plot(timeref, dataref[:, 3], 'g', lw=2, ls='--', label=model + '(Ref)')
|
|
hy1.plot(timetest, datatest[:, 4], 'r', lw=2, label=model)
|
|
hy1.plot(timeref, dataref[:, 4], 'g', lw=2, ls='--', label=model + '(Ref)')
|
|
hz1.plot(timetest, datatest[:, 5], 'r', lw=2, label=model)
|
|
hz1.plot(timeref, dataref[:, 5], 'g', lw=2, ls='--', label=model + '(Ref)')
|
|
ylabels = ['$E_x$, field strength [V/m]', '$H_x$, field strength [A/m]',
|
|
'$E_y$, field strength [V/m]', '$H_y$, field strength [A/m]',
|
|
'$E_z$, field strength [V/m]', '$H_z$, field strength [A/m]']
|
|
for i, ax in enumerate(fig1.axes):
|
|
ax.set_ylabel(ylabels[i])
|
|
ax.set_xlim(0, np.amax(timetest))
|
|
ax.grid()
|
|
ax.legend()
|
|
|
|
# Plot diffs
|
|
fig2, ((ex2, hx2), (ey2, hy2), (ez2, hz2)) = plt.subplots(nrows=3, ncols=2,
|
|
sharex=False, sharey='col',
|
|
subplot_kw=dict(xlabel='Time [ns]'),
|
|
num='Diffs: ' + model + '.in',
|
|
figsize=(20, 10),
|
|
facecolor='w',
|
|
edgecolor='w')
|
|
ex2.plot(timeref, datadiffs[:, 0], 'r', lw=2, label='Ex')
|
|
ey2.plot(timeref, datadiffs[:, 1], 'r', lw=2, label='Ey')
|
|
ez2.plot(timeref, datadiffs[:, 2], 'r', lw=2, label='Ez')
|
|
hx2.plot(timeref, datadiffs[:, 3], 'r', lw=2, label='Hx')
|
|
hy2.plot(timeref, datadiffs[:, 4], 'r', lw=2, label='Hy')
|
|
hz2.plot(timeref, datadiffs[:, 5], 'r', lw=2, label='Hz')
|
|
ylabels = ['$E_x$, difference [dB]', '$H_x$, difference [dB]',
|
|
'$E_y$, difference [dB]', '$H_y$, difference [dB]',
|
|
'$E_z$, difference [dB]', '$H_z$, difference [dB]']
|
|
for i, ax in enumerate(fig2.axes):
|
|
ax.set_ylabel(ylabels[i])
|
|
ax.set_xlim(0, np.amax(timetest))
|
|
ax.set_ylim([plotmin, np.amax(np.amax(datadiffs))])
|
|
ax.grid()
|
|
|
|
# Save a PDF/PNG of the figure
|
|
filediffs = file.stem + '_diffs'
|
|
filediffs = file.parent / Path(filediffs)
|
|
# fig1.savefig(file.with_suffix('.pdf'), dpi=None, format='pdf',
|
|
# bbox_inches='tight', pad_inches=0.1)
|
|
# fig2.savefig(savediffs.with_suffix('.pdf'), dpi=None, format='pdf',
|
|
# bbox_inches='tight', pad_inches=0.1)
|
|
fig1.savefig(file.with_suffix('.png'), dpi=150, format='png',
|
|
bbox_inches='tight', pad_inches=0.1)
|
|
fig2.savefig(filediffs.with_suffix('.png'), dpi=150, format='png',
|
|
bbox_inches='tight', pad_inches=0.1)
|
|
|
|
# Summary of results
|
|
for name, data in sorted(testresults.items()):
|
|
if 'analytical' in name:
|
|
logger.info(f"Test '{name}.in' using v.{data['Test version']} compared " +
|
|
f"to analytical solution. Max difference {data['Max diff']:.2f}dB.")
|
|
else:
|
|
logger.info(f"Test '{name}.in' using v.{data['Test version']} compared to " +
|
|
f"reference solution using v.{data['Ref version']}. Max difference " +
|
|
f"{data['Max diff']:.2f}dB.")
|