Merge pull request #373 from nauaneed/minor-fixes

Minor fixes and enhancements

Author: prabhuramachandran

pysph/examples/gas_dynamics/shocktube.py

@@ -36,18 +36,15 @@
 # domain size
 xmin = 0.
 xmax = 1
-dx = 0.002
-ny = 50
 ymin = 0
-ymax = ny * dx
-x0 = 0.5  # initial discontuinity
+ymax = 0.1
+x0 = 0.5  # initial discontinuity
 
 # scheme constants
 alpha1 = 1.0
 alpha2 = 1.0
 beta = 2.0
 kernel_factor = 1.5
-h0 = kernel_factor * dx
 
 
 class ShockTube2D(Application):
@@ -56,10 +53,8 @@ def initialize(self):
         self.xmax = xmax
         self.ymin = ymin
         self.ymax = ymax
-        self.dx = dx
         self.hdx = 1.7
         self.x0 = x0
-        self.ny = ny
         self.pl = 1000
         self.pr = 0.01
         self.rhol = 1.0
@@ -72,17 +67,22 @@ def initialize(self):
     def add_user_options(self, group):
         add_bool_argument(group, 'smooth-ic', dest='smooth_ic', default=False,
                           help="Smooth the initial condition.")
+        group.add_argument(
+            "--dx", action="store", type=float, dest="dx",
+            default=0.002, help="Particle spacing"
+        )
 
     def consume_user_options(self):
         self.smooth_ic = self.options.smooth_ic
+        self.dx = self.options.dx
 
     def create_domain(self):
         return DomainManager(
             xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
             periodic_in_x=True, periodic_in_y=True)
 
     def create_particles(self):
-        global dx
+        dx = self.dx
         data = ud.uniform_distribution_cubic2D(dx, xmin, xmax, ymin, ymax)
 
         x = data[0]
@@ -110,7 +110,7 @@ def create_particles(self):
             p[right_indices] = self.pr
 
         # const h and mass
-        h = numpy.ones_like(x) * self.hdx * self.dx
+        h = numpy.ones_like(x) * self.hdx * dx
         m = numpy.ones_like(x) * volume * rho
 
         # ul = ur = 0

pysph/sph/gas_dynamics/magma2.py

@@ -708,7 +708,7 @@ def loop(self, d_idx, s_m, s_idx, XIJ, s_rho, d_cm, WI):
                 drowcol = dsi2 + row * dim + col
                 d_cm[drowcol] += mbbyrhob * XIJ[row] * XIJ[col] * WI
 
-    def post_loop(self, d_idx, d_dv, d_divv, d_cm):
+    def post_loop(self, d_idx, d_cm):
         invcm, cm, idmat = declare('matrix(9)', 3)
         augcm = declare('matrix(18)')
         dsi2, row, col, rowcol, dim, dimsq = declare('int', 6)
@@ -926,7 +926,7 @@ def post_loop(self, d_h, d_idx, d_cs, d_alpha, d_s, d_p, d_rho, dt,
 class WallBoundary(Equation):
     """:class:`WallBoundary
     <pysph.sph.gas_dynamics.boundary_equations.WallBoundary>` modified
-    for GADGET2.
+    for MAGMA2.
     """
 
     def __init__(self, dest, sources, dim):
@@ -1058,7 +1058,7 @@ def post_loop(self, d_idx, d_p, d_rho, d_e, d_m, d_cs, d_divv, d_h, d_u,
 class UpdateGhostProps(Equation):
     """
     :class:`MPMUpdateGhostProps
-    <pysph.sph.gas_dynamics.basic.MPMUpdateGhostProps>` modified for GADGET2.
+    <pysph.sph.gas_dynamics.basic.MPMUpdateGhostProps>` modified for MAGMA2.
     """
 
     def __init__(self, dest, dim, sources=None):

pysph/sph/gas_dynamics/psph.py

@@ -201,9 +201,9 @@ def get_equations(self):
                 dest=fluid, sources=all_pa, hfact=self.hfact,
                 density_iterations=True, dim=self.dim,
                 htol=self.density_iteration_tolerance, gamma=self.gamma))
-            equations.append(
-                Group(equations=g1, update_nnps=True, iterate=True,
-                      max_iterations=self.max_density_iterations))
+        equations.append(
+            Group(equations=g1, update_nnps=True, iterate=True,
+                  max_iterations=self.max_density_iterations))
 
         g2 = []
         for fluid in self.fluids:

pysph/sph/gas_dynamics/tsph.py

@@ -161,9 +161,9 @@ def get_equations(self):
                 dest=fluid, sources=all_pa, hfact=self.hfact,
                 density_iterations=True, dim=self.dim,
                 htol=self.density_iteration_tolerance))
-            equations.append(Group(
-                equations=g1, update_nnps=True, iterate=True,
-                max_iterations=self.max_density_iterations))
+        equations.append(Group(
+            equations=g1, update_nnps=True, iterate=True,
+            max_iterations=self.max_density_iterations))
 
         g2 = []
         for fluid in self.fluids:

pysph/sph/scheme.py

@@ -1171,8 +1171,8 @@ def __init__(self, fluids, solids, dim, gamma, kernel_factor, g1=0.0,
         monotonicity : int
             Type of monotonicity algorithm to use:
             0 : First order GSPH
-            1 : I02 algorithm
-            2 : IwIn algorithm
+            1 : I02 algorithm # https://doi.org/10.1006/jcph.2002.7053
+            2 : IwIn algorithm # https://doi.org/10.1111/j.1365-2966.2011.19588.x  # noqa: E501
         interface_zero : bool
             Set Interface position s^*_{ij} = 0 for the Riemann problem.
         hybrid, blend_alpha : bool, double
@@ -1204,22 +1204,42 @@ def __init__(self, fluids, solids, dim, gamma, kernel_factor, g1=0.0,
         self.niter = niter
         self.tol = tol
         self.has_ghosts = has_ghosts
+        self.rsolver_choices = {'non_diffusive': 0,
+                                'van_leer': 1,
+                                'exact': 2,
+                                'hllc': 3,
+                                'ducowicz': 4,
+                                'hlle': 5,
+                                'roe': 6,
+                                'llxf': 7,
+                                'hllc_ball': 8,
+                                'hll_ball': 9,
+                                'hllsy': 10}
+        self.interpolation_choices = {'delta': 0,
+                                      'linear': 1,
+                                      'cubic': 2}
+        self.monotonicity_choices = {'first_order': 0,
+                                     'i02': 1,
+                                     'iwin': 2}
 
     def add_user_options(self, group):
         group.add_argument(
-            "--rsolver", action="store", type=int, dest="rsolver",
-            default=None,
-            help="Riemann solver to use."
+            "--rsolver", action="store", type=str, dest="rsolver",
+            default=None, choices=set(self.rsolver_choices.keys()),
+            help=f"Riemann solver to use, one of :"
+                 f"{set(self.rsolver_choices.keys())}"
         )
         group.add_argument(
-            "--interpolation", action="store", type=int, dest="interpolation",
-            default=None,
-            help="Interpolation algorithm to use."
+            "--interpolation", action="store", type=str, dest="interpolation",
+            default=None, choices=set(self.interpolation_choices.keys()),
+            help=f"Interpolation algorithm to use, one of :"
+                 f"{set(self.interpolation_choices.keys())}"
         )
         group.add_argument(
-            "--monotonicity", action="store", type=int, dest="monotonicity",
-            default=None,
-            help="Monotonicity algorithm to use."
+            "--monotonicity", action="store", type=str, dest="monotonicity",
+            default=None, choices=set(self.monotonicity_choices.keys()),
+            help=f"Monotonicity algorithm to use, one of :"
+                 f"{set(self.monotonicity_choices.keys())}"
         )
         group.add_argument(
             "--g1", action="store", type=float, dest="g1",
@@ -1253,13 +1273,24 @@ def add_user_options(self, group):
         )
 
     def consume_user_options(self, options):
-        vars = ['gamma', 'g1', 'g2', 'rsolver', 'interpolation',
-                'monotonicity', 'interface_zero', 'hybrid',
-                'blend_alpha']
+        vars = ['gamma', 'g1', 'g2', 'interface_zero',
+                'hybrid', 'blend_alpha']
         data = dict((var, self._smart_getattr(options, var))
                     for var in vars)
+        map_vars = ['monotonicity', 'rsolver', 'interpolation']
+        for var in map_vars:
+            data[var] = self._smart_getattr_mapped(options, var)
         self.configure(**data)
 
+    def _smart_getattr_mapped(self, obj, var):
+        res = getattr(obj, var)
+        if res is None:
+            return getattr(self, var)
+        else:
+            choices = getattr(self, f'{var}_choices')
+            return choices[res]
+
+
     def configure_solver(self, kernel=None, integrator_cls=None,
                          extra_steppers=None, **kw):
         """Configure the solver to be generated.