From e831f0a713362caf03f2baa57fc86f45b23a3a1b Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 6 Oct 2016 17:08:52 +0100 Subject: [PATCH] Added more informative messages (via analysis) on numerical dispersion. --- gprMax/gprMax.py | 12 +++++++----- gprMax/grid.py | 44 ++++++++++++++++++++++++++++++-------------- 2 files changed, 37 insertions(+), 19 deletions(-) diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index 6b4540fa..03858de3 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -38,7 +38,7 @@ from .constants import c, e0, m0, z0 from .exceptions import GeneralError from .fields_outputs import write_hdf5_outputfile from .fields_update import update_electric, update_magnetic, update_electric_dispersive_multipole_A, update_electric_dispersive_multipole_B, update_electric_dispersive_1pole_A, update_electric_dispersive_1pole_B -from .grid import FDTDGrid, dispersion_check +from .grid import FDTDGrid, dispersion_analysis from .input_cmds_geometry import process_geometrycmds from .input_cmds_file import process_python_include_code, write_processed_file, check_cmd_names from .input_cmds_multiuse import process_multicmds @@ -46,7 +46,7 @@ from .input_cmds_singleuse import process_singlecmds from .materials import Material, process_materials from .pml import build_pmls from .receivers import store_outputs -from .utilities import logo, human_size, get_machine_cpu_os, get_terminal_width +from .utilities import logo, human_size, get_machine_cpu_os, get_terminal_width, round_value from .yee_cell_build import build_electric_components, build_magnetic_components @@ -418,9 +418,11 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace): print(materialstable.table) # Check to see if numerical dispersion might be a problem - resolution = dispersion_check(G) - if resolution and max((G.dx, G.dy, G.dz)) > resolution: - print(Fore.RED + '\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution <{:g}m'.format(resolution) + Style.RESET_ALL) + deltavp, N, material, maxfreq = dispersion_analysis(G) + if deltavp and np.abs(deltavp) > G.maxnumericaldisp: + print(Fore.RED + "\nWARNING: Potentially significant numerical dispersion. Largest physical phase-velocity error is {:.2f}% in material '{}' with wavelength sampled by {} cells (maximum significant frequency {:g}Hz)".format(deltavp, material.ID, round_value(N), maxfreq) + Style.RESET_ALL) + elif deltavp: + print("\nNumerical dispersion analysis: largest physical phase-velocity error is {:.2f}% in material '{}' with wavelength sampled by {} cells (maximum significant frequency {:g}Hz)".format(deltavp, material.ID, round_value(N), maxfreq)) # If geometry information to be reused between model runs else: diff --git a/gprMax/grid.py b/gprMax/grid.py index 9dafd0b0..525b60d6 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -75,6 +75,12 @@ class FDTDGrid(Grid): self.title = '' self.messages = True self.tqdmdisable = False + + # Threshold (dB) down from maximum power (0dB) of main frequency used to calculate highest frequency for disperion analysis + self.highestfreqthres = 60 + # Maximum allowable percentage physical phase-velocity phase error + self.maxnumericaldisp = 2 + self.nx = 0 self.ny = 0 self.nz = 0 @@ -136,19 +142,16 @@ class FDTDGrid(Grid): self.updatecoeffsdispersive = np.zeros((len(self.materials), 3 * Material.maxpoles), dtype=complextype) -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 +def dispersion_analysis(G): + """Analysis of numerical dispersion (Taflove et al, 2005, p112) - worse case of maximum frequency and minimum wavelength Args: G (class): Grid class instance - holds essential parameters describing the model. Returns: - resolution (float): Potential numerical dispersion + deltavp (float): Percentage difference between true and numerical phase velocity """ - # Minimum number of spatial steps to resolve smallest wavelength - resolvedsteps = 10 - # Find maximum frequency maxfreqs = [] for waveform in G.waveforms: @@ -185,7 +188,7 @@ def dispersion_check(G): power -= np.amax(power) # Set maximum frequency to -60dB from maximum power, ignoring DC value - freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1 + freq = np.where((np.amax(power[1::]) - power[1::]) > G.highestfreqthres)[0][0] + 1 maxfreqs.append(freqs[freq]) else: @@ -194,23 +197,36 @@ def dispersion_check(G): if maxfreqs: maxfreq = max(maxfreqs) - # Find minimum wavelength - ers = [material.er for material in G.materials] + # Find minimum wavelength (material with maximum permittivity) + ers = [x.er for x in G.materials if x.ID != 'pec'] maxer = max(ers) + material = next(x for x in G.materials if x.er == maxer and x.ID != 'pec') # Minimum velocity minvelocity = c / np.sqrt(maxer) # Minimum wavelength minwavelength = minvelocity / maxfreq - - # Resolution of minimum wavelength - resolution = minwavelength / resolvedsteps + + # Maximum spatial step + delta = max(G.dx, G.dy, G.dz) + + # Courant stability factor + S = (c * G.dt) / delta + + # Grid sampling density + N = minwavelength / delta + + # Numerical phase velocity + vp = np.pi / (N * np.arcsin((1 / S) * np.sin((np.pi * S) / N))) + + # Physical phase velocity error (percentage) + deltavp = (((vp * c) - c) / c) * 100 else: - resolution = False + deltavp, N, material, maxfreq = False - return resolution + return deltavp, N, material, maxfreq def get_other_directions(direction):