Updates to PML functionality

gprMax/__init__.py

@@ -28,8 +28,9 @@ from .cmds_multiuse import (PMLCFS, AddDebyeDispersion, AddDrudeDispersion,
                             Material, Rx, RxArray, Snapshot, SoilPeplinski,
                             TransmissionLine, VoltageSource, Waveform)
 from .cmds_singleuse import (Discretisation, Domain, ExcitationFile,
-                             OMPThreads, PMLCells, RxSteps, SrcSteps,
-                             TimeStepStabilityFactor, TimeWindow, Title)
+                             OMPThreads, PMLCells, PMLFormulation, 
+                             RxSteps, SrcSteps, TimeStepStabilityFactor, 
+                             TimeWindow, Title)
 from .gprMax import run as run
 from .hash_cmds_file import user_libs_fn_to_scene_obj
 from .scene import Scene

gprMax/cmds_multiuse.py

@@ -1582,17 +1582,10 @@ class PMLCFS(UserObjectMulti):
         sigmamax: float required for maximum value for the CFS sigma parameter.
     """
 
-    count = 0
-
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
         self.order = 16
         self.hash = '#pml_cfs'
-        PMLCFS.count += 1
-        if PMLCFS.count == 2:
-            logger.exception(self.params_str() + ' can only be used up to two '
-                             'times, for up to a 2nd order PML.')
-            raise ValueError
 
     def create(self, grid, uip):
         try:
@@ -1630,10 +1623,10 @@ class PMLCFS(UserObjectMulti):
             logger.exception(self.params_str() + ' minimum and maximum scaling '
                              'values must be greater than zero.')
             raise ValueError
-        if float(kappamin) < 1:
-            logger.exception(self.params_str() + ' minimum scaling value for '
-                             'kappa must be greater than or equal to one.')
-            raise ValueError
+        # if float(kappamin) < 1:
+        #     logger.exception(self.params_str() + ' minimum scaling value for '
+        #                      'kappa must be greater than or equal to one.')
+        #     raise ValueError
 
         cfsalpha = CFSParameter()
         cfsalpha.ID = 'alpha'
@@ -1652,9 +1645,7 @@ class PMLCFS(UserObjectMulti):
         cfssigma.scalingprofile = sigmascalingprofile
         cfssigma.scalingdirection = sigmascalingdirection
         cfssigma.min = float(sigmamin)
-        if sigmamax == 'None':
-            cfssigma.max = None
-        else:
+        if sigmamax is not None:
             cfssigma.max = float(sigmamax)
         cfs = CFS()
         cfs.alpha = cfsalpha
@@ -1669,9 +1660,14 @@ class PMLCFS(UserObjectMulti):
                     f'{cfskappa.max:g}), sigma (scaling: {cfssigma.scalingprofile}, '
                     f'scaling direction: {cfssigma.scalingdirection}, min: '
                     f'{cfssigma.min:g}, max: {cfssigma.max:g}) created.')
-
+        
         grid.cfs.append(cfs)
 
+        if len(grid.cfs) > 2:
+            logger.exception(self.params_str() + ' can only be used up to two '
+                             'times, for up to a 2nd order PML.')
+            raise ValueError
+            
 
 class Subgrid(UserObjectMulti):
     """"""

gprMax/cmds_singleuse.py

@@ -24,6 +24,7 @@ import gprMax.config as config
 import numpy as np
 from scipy import interpolate
 
+from .pml import PML
 from .utilities.host_info import set_omp_threads
 from .waveforms import Waveform
 
@@ -254,6 +255,30 @@ class TimeStepStabilityFactor(UserObjectSingle):
         logger.info(f'Time step (modified): {G.dt:g} secs')
 
 
+class PMLFormulation(UserObjectSingle):
+    """Allows you to specify the formulation (type) of the PML to be used.
+
+    Attributes:
+        pml: string specifying formulation of PML.
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.order = 8
+
+    def create(self, G, uip):
+        try:
+            pml = self.kwargs['pml']
+        except KeyError:
+            logger.exception(self.__str__() + ' requires exactly one parameter to specify the formulation of PML to use')
+            raise
+        if pml not in PML.formulations:
+            logger.exception(self.__str__() + f" requires the value to be one of {' '.join(PML.formulations)}")
+            raise ValueError
+
+        G.pmlformulation = pml
+
+
 class PMLCells(UserObjectSingle):
     """Allows you to control the number of cells (thickness) of PML that are used
     on the six sides of the model domain. Specify either single thickness or
@@ -277,7 +302,7 @@ class PMLCells(UserObjectSingle):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
-        self.order = 8
+        self.order = 9
 
     def create(self, G, uip):
         try:
@@ -317,7 +342,7 @@ class SrcSteps(UserObjectSingle):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
-        self.order = 9
+        self.order = 10
 
     def create(self, G, uip):
         try:
@@ -339,7 +364,7 @@ class RxSteps(UserObjectSingle):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
-        self.order = 10
+        self.order = 11
 
     def create(self, G, uip):
         try:
@@ -366,7 +391,7 @@ class ExcitationFile(UserObjectSingle):
 
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
-        self.order = 11
+        self.order = 12
 
     def create(self, G, uip):
         try:
@@ -446,7 +471,7 @@ class OutputDir(UserObjectSingle):
     """
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
-        self.order = 12
+        self.order = 13
 
     def create(self, grid, uip):
         config.get_model_config().set_output_file_path(self.kwargs['dir'])
@@ -461,7 +486,7 @@ class NumberOfModelRuns(UserObjectSingle):
     """
     def __init__(self, **kwargs):
         super().__init__(**kwargs)
-        self.order = 13
+        self.order = 14
 
     def create(self, grid, uip):
         try:

gprMax/grid.py

@@ -21,8 +21,8 @@ from collections import OrderedDict
 
 import gprMax.config as config
 import numpy as np
-from .pml import CFS, PML
-from .utilities.utilities import fft_power, human_size, round_value
+from .pml import PML
+from .utilities.utilities import fft_power, round_value
 
 np.seterr(invalid='raise')
 
@@ -53,7 +53,8 @@ class FDTDGrid:
         # same from model to model otherwise the numerical precision from adding
         # the PML corrections will be different.
         self.pmlthickness = OrderedDict((key, 10) for key in PML.boundaryIDs)
-        self.cfs = [CFS()]
+        # Default PML CFS parameters will be used if none are provided by user.
+        self.cfs = []
         self.pmls = []
         self.pmlformulation = 'HORIPML'
 

gprMax/hash_cmds_file.py

@@ -201,8 +201,7 @@ def check_cmd_names(processedlines, checkessential=True):
     # Commands that there should only be one instance of in a model
     singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', 
                                 '#title', '#cpu_threads', 
-                                '#time_step_stability_factor', 
-                                '#pml_formulation', '#pml_cells', 
+                                '#time_step_stability_factor', '#pml_cells', 
                                 '#excitation_file', '#src_steps', '#rx_steps', 
                                 '#output_dir'], None)
 

gprMax/pml.py

@@ -215,7 +215,10 @@ class PML:
             self.d = self.G.dz
             self.thickness = self.nz
 
-        self.CFS = self.G.cfs
+        if not self.G.cfs:
+            self.CFS = [CFS()]
+        else:
+            self.CFS = self.G.cfs
 
         self.initialise_field_arrays()
 
@@ -319,7 +322,7 @@ class PML:
                 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.sim_config.em_consts['e0']) 
+                self.HRE[x, :] = (((2 * config.sim_config.em_consts['e0']) 
                                   - self.G.dt * Halpha) / tmp)
                 self.HRF[x, :] = (2 * Hsigma * self.G.dt) / tmp