Restructured CFS and CFS parameters for PML. Corrected implementation of kappa scaling when not constant.

gprMax/pml.py

@@ -21,35 +21,46 @@ import numpy as np
 from gprMax.constants import e0, z0, floattype
 
 
-class CFS():
+class CFSParameter():
-    """PML CFS parameters."""
+    """CFS parameter for PML."""
     
-    # Allowable scaling types
+    # Allowable scaling profiles and directions
-    scalingtypes = {'constant': 0, 'linear': 1, 'inverselinear': -1, 'quadratic': 2, 'cubic': 3, 'quartic': 4}
+    scalingprofiles = {'constant': 0, 'linear': 1, 'quadratic': 2, 'cubic': 3, 'quartic': 4}
+    scalingdirections = ['forward', 'reverse']
     
-    def __init__(self, alphascaling='constant', alphamin=0, alphamax=0, kappascaling='constant', kappamin=1, kappamax=1, sigmascaling='quartic', sigmamin=0, sigmamax=None):
+    def __init__(self, ID=None, scaling='polynomial', scalingprofile=None, scalingdirection='forward', min=0, max=0):
         """
         Args:
-           alphascaling (str): Type of scaling used for alpha parameter. Can be: 'constant', 'linear', 'inverselinear', 'quadratic', 'cubic', 'quartic'.
+            ID (str): Identifier for CFS parameter, can be: 'alpha', 'kappa' or 'sigma'.
-           alphamin (float): Minimum value for alpha parameter.
+            scaling (str): Type of scaling, can be: 'polynomial'.
-           alphamax (float): Maximum value for alpha parameter.
+            scalingprofile (str): Type of scaling profile from scalingprofiles.
-           kappascaling (str): Type of scaling used for kappa parameter. Can be: 'constant', 'linear', 'inverselinear', 'quadratic', 'cubic', 'quartic'.
+            scalingdirection (str): Direction of scaling profile from scalingdirections.
-           kappamin (float): Minimum value for kappa parameter.
+            min (float): Minimum value for parameter.
-           kappamax (float): Maximum value for kappa parameter.
+            max (float): Maximum value for parameter.
-           sigmascaling (str): Type of scaling used for sigma parameter. Can be: 'constant', 'linear', 'inverselinear', 'quadratic', 'cubic', 'quartic'.
-           sigmamin (float): Minimum value for sigma parameter.
-           sigmamax (float): Maximum value for sigma parameter.
         """
         
-        self.alphascaling = alphascaling
+        self.ID = ID
-        self.alphamin = alphamin
+        self.scaling = scaling
-        self.alphamax = alphamax
+        self.scalingprofile = scalingprofile
-        self.kappascaling = kappascaling
+        self.scalingdirection = scalingdirection
-        self.kappamin = kappamin
+        self.min = min
-        self.kappamax = kappamax
+        self.max = max
-        self.sigmascaling = sigmascaling
+
-        self.sigmamin = sigmamin
+
-        self.sigmamax = sigmamax
+class CFS():
+    """CFS term for PML."""
+    
+    def __init__(self):
+        """
+        Args:
+            alpha (CFSParameter): alpha parameter for CFS.
+            kappa (CFSParameter): kappa parameter for CFS.
+            sigma (CFSParameter): sigma parameter for CFS.
+        """
+        
+        self.alpha = CFSParameter(ID='alpha', scalingprofile='constant')
+        self.kappa = CFSParameter(ID='kappa', scalingprofile='constant', min=1, max=1)
+        self.sigma = CFSParameter(ID='sigma', scalingprofile='quartic', min=0, max=None)
     
     def calculate_sigmamax(self, direction, er, mr, G):
         """Calculates an optimum value for sigma max based on underlying material properties.
@@ -61,63 +72,69 @@ class CFS():
             G (class): Grid class instance - holds essential parameters describing the model.
         """
         
-        # Get general direction from first letter of PML direction
+        if direction == 'xminus' or direction == 'xplus':
-        if direction[0] == 'x':
             d = G.dx
-        elif direction[0] == 'y':
+        elif direction == 'yminus' or direction == 'yplus':
             d = G.dy
-        elif direction[0] == 'z':
+        elif direction == 'zminus' or direction == 'zplus':
             d = G.dz
-        # Calculation of the maximum value of sigma from http://dx.doi.org/10.1109/8.546249
-        m = CFS.scalingtypes[self.sigmascaling]
-        self.sigmamax = (0.8 * (m + 1)) / (z0 * d * np.sqrt(er * mr))
         
-    def scaling_polynomial(self, min, max, order, Evalues, Hvalues):
+        # Calculation of the maximum value of sigma from http://dx.doi.org/10.1109/8.546249
-        """Applies the polynomial to be used for scaling for electric and magnetic PML updates based on scaling type and minimum and maximum values.
+        m = CFSParameter.scalingprofiles[self.sigma.scalingprofile]
+        self.sigma.max = (0.8 * (m + 1)) / (z0 * d * np.sqrt(er * mr))
+
+    def scaling_polynomial(self, order, Evalues, Hvalues):
+        """Applies the polynomial to be used for the scaling profile for electric and magnetic PML updates.
             
         Args:
-            min (float): Minimum value for scaling.
+            order (int): Order of polynomial for scaling profile.
-            max (float): Maximum value for scaling.
+            Evalues (float): numpy array holding scaling profile values for electric PML update.
-            order (int): Order of polynomial for scaling.
+            Hvalues (float): numpy array holding scaling profile values for magnetic PML update.
-            Evalues (float): numpy array holding scaling value for electric PML update.
-            Hvalues (float): numpy array holding scaling value for magnetic PML update.
             
         Returns:
-            Evalues (float): numpy array holding scaling value for electric PML update.
+            Evalues (float): numpy array holding scaling profile values for electric PML update.
-            Hvalues (float): numpy array holding scaling value for magnetic PML update.
+            Hvalues (float): numpy array holding scaling profile values for magnetic PML update.
         """
         
-        tmp = max * ((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]
         return Evalues, Hvalues
 
-    def calculate_values(self, min, max, scaling, Evalues, Hvalues):
+    def calculate_values(self, thickness, parameter):
-        """Calculates values for electric and magnetic PML updates based on scaling type and minimum and maximum values.
+        """Calculates values for electric and magnetic PML updates based on profile type and minimum and maximum values.
             
         Args:
-            min (float): Minimum value for scaling.
+            thickness (int): Thickness of PML in cells.
-            max (float): Maximum value for scaling.
+            parameter (CFSParameter): Instance of CFSParameter
-            scaling (int): Type of scaling, can be: 'constant', 'linear', 'inverselinear', 'quadratic', 'cubic', 'quartic'.
-            Evalues (float): numpy array holding scaling value for electric PML update.
-            Hvalues (float): numpy array holding scaling value for magnetic PML update.
             
         Returns:
-            Evalues (float): numpy array holding scaling value for electric PML update.
+            Evalues (float): numpy array holding profile value for electric PML update.
-            Hvalues (float): numpy array holding scaling value for magnetic PML update.
+            Hvalues (float): numpy array holding profile value for magnetic PML update.
         """
         
-        if scaling == 'constant':
+        Evalues= np.zeros(thickness + 1, dtype=floattype)
-            Evalues += max
+        Hvalues = np.zeros(thickness + 1, dtype=floattype)
-            Hvalues += max
-        else:
-            Evalues, Hvalues = self.scaling_polynomial(min, max, CFS.scalingtypes[scaling], Evalues, Hvalues)
         
-        if scaling == 'inverselinear':
+        if parameter.scalingprofile == 'constant':
+            Evalues += parameter.max
+            Hvalues += parameter.max
+        
+        elif parameter.scaling == 'polynomial':
+            Evalues, Hvalues = self.scaling_polynomial(CFSParameter.scalingprofiles[parameter.scalingprofile], Evalues, Hvalues)
+            if parameter.ID == 'alpha':
+                pass
+            elif parameter.ID == 'kappa':
+                Evalues = Evalues * (self.kappa.max - 1) + 1
+                Hvalues = Hvalues * (self.kappa.max - 1) + 1
+            elif parameter.ID == 'sigma':
+                Evalues *= self.sigma.max
+                Hvalues *= self.sigma.max
+            
+        if parameter.scalingdirection == 'reverse':
             Evalues = Evalues[::-1]
             Hvalues = Hvalues[::-1]
-#        print('Evalues: scaling {}, {}'.format(scaling, Evalues))
+        
-#        print('Hvalues: scaling {}, {}'.format(scaling, Hvalues))
         return Evalues, Hvalues
 
 
@@ -186,17 +203,11 @@ class PML():
         """
         
         for x, cfs in enumerate(self.CFS):
-            Ealpha = np.zeros(self.thickness + 1, dtype=floattype)
+            if not cfs.sigma.max:
-            Halpha = np.zeros(self.thickness + 1, dtype=floattype)
-            Ekappa = np.zeros(self.thickness + 1, dtype=floattype)
-            Hkappa = np.zeros(self.thickness + 1, dtype=floattype)
-            Esigma = np.zeros(self.thickness + 1, dtype=floattype)
-            Hsigma = np.zeros(self.thickness + 1, dtype=floattype)
-            if not cfs.sigmamax:
                 cfs.calculate_sigmamax(self.direction, er, mr, G)
-            Ealpha, Halpha = cfs.calculate_values(cfs.alphamin, cfs.alphamax, cfs.alphascaling, Ealpha, Halpha)
+            Ealpha, Halpha = cfs.calculate_values(self.thickness, cfs.alpha)
-            Ekappa, Hkappa = cfs.calculate_values(cfs.kappamin, cfs.kappamax, cfs.kappascaling, Ekappa, Hkappa)
+            Ekappa, Hkappa = cfs.calculate_values(self.thickness, cfs.kappa)
-            Esigma, Hsigma = cfs.calculate_values(cfs.sigmamin, cfs.sigmamax, cfs.sigmascaling, Esigma, Hsigma)
+            Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma)
             
 #            print('Ealpha {}'.format(Ealpha))
 #            print('Halpha {}'.format(Halpha))
@@ -273,7 +284,7 @@ def build_pml(G):
 
 
 def calculate_initial_pml_params(G):
-    """ This function calculates the initial parameters and coefficients for PML including setting scaling
+    """ This function calculates the initial parameters and coefficients for PML including setting profile
         (based on underlying material er and mr from solid array).
     """