Merge pull request #334 from nauaneed/TSPH-PSPH-KISS

Hopkins' TSPH and PSPH with squashed commits.

docs/source/reference/equations.rst

@@ -59,6 +59,10 @@ Gas Dynamics
    :members:
    :undoc-members:
 
+.. automodule:: pysph.sph.gas_dynamics.boundary_equations
+   :members:
+   :undoc-members:
+
 Surface tension
 ----------------
 
@@ -73,6 +77,24 @@ Implicit Incompressible SPH
    :members:
    :undoc-members:
 
+Hopkins' ‘Traditional’ SPH (TSPH)
+--------------
+
+.. automodule:: pysph.sph.gas_dynamics.tsph
+   :members: TSPHScheme, SummationDensity, IdealGasEOS, VelocityGradDivC1,
+             BalsaraSwitch, WallBoundary, UpdateGhostProps, MomentumAndEnergy
+   :undoc-members:
+   :member-order: bysource
+
+Hopkins' ‘Modern’ SPH (PSPH)
+--------------
+
+.. automodule:: pysph.sph.gas_dynamics.psph
+   :members: PSPHScheme, PSPHSummationDensityAndPressure, GradientKinsfolkC1,
+             VelocityGradDivC1, LimiterAndAlphas, WallBoundary,
+             UpdateGhostProps, MomentumAndEnergy
+   :undoc-members:
+   :member-order: bysource
 
 Rigid body motion
 -----------------
@@ -92,12 +114,11 @@ Miscellaneous
    :members:
    :undoc-members:
 
-
 Group of equations
 -------------------
 
 .. autoclass:: pysph.sph.equation.Group
    :special-members:
 
 .. autoclass:: pysph.sph.equation.MultiStageEquations
-   :special-members:
+   :special-members:

pysph/examples/gas_dynamics/accuracy_test_2d.py

@@ -13,6 +13,8 @@
 
 from pysph.sph.scheme import GasDScheme, ADKEScheme, GSPHScheme, SchemeChooser
 from pysph.sph.wc.crksph import CRKSPHScheme
+from pysph.sph.gas_dynamics.psph import PSPHScheme
+from pysph.sph.gas_dynamics.tsph import TSPHScheme
 
 # PySPH tools
 from pysph.tools import uniform_distribution as ud
@@ -140,8 +142,19 @@ def create_scheme(self):
             niter=40, tol=1e-6, has_ghosts=True
         )
 
+        psph = PSPHScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            hfact=kernel_factor
+        )
+
+        tsph = TSPHScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            hfact=kernel_factor
+        )
+
         s = SchemeChooser(
-            default='gsph', adke=adke, mpm=mpm, gsph=gsph, crksph=crksph
+            default='gsph', adke=adke, mpm=mpm, gsph=gsph, crksph=crksph,
+            psph=psph, tsph=tsph
         )
         return s
 
@@ -151,13 +164,7 @@ def configure_scheme(self):
             s.configure(kernel_factor=kernel_factor)
             s.configure_solver(dt=self.dt, tf=self.tf,
                                adaptive_timestep=True, pfreq=50)
-        elif self.options.scheme == 'adke':
-            s.configure_solver(dt=self.dt, tf=self.tf,
-                               adaptive_timestep=False, pfreq=50)
-        elif self.options.scheme == 'gsph':
-            s.configure_solver(dt=self.dt, tf=self.tf,
-                               adaptive_timestep=False, pfreq=50)
-        elif self.options.scheme == 'crksph':
+        elif self.options.scheme in ['adke', 'gsph', 'crksph', 'tsph', 'psph']:
             s.configure_solver(dt=self.dt, tf=self.tf,
                                adaptive_timestep=False, pfreq=50)
 

pysph/examples/gas_dynamics/acoustic_wave.py

@@ -25,6 +25,8 @@
 from pysph.sph.scheme import \
     GSPHScheme, ADKEScheme, GasDScheme, SchemeChooser
 from pysph.sph.wc.crksph import CRKSPHScheme
+from pysph.sph.gas_dynamics.psph import PSPHScheme
+from pysph.sph.gas_dynamics.tsph import TSPHScheme
 
 
 class AcousticWave(Application):
@@ -113,8 +115,18 @@ def create_scheme(self):
             alpha=0, beta=0.0, k=1.5, eps=0.0, g1=0.0, g2=0.0,
             has_ghosts=True)
 
+        psph = PSPHScheme(
+            fluids=['fluid'], solids=[], dim=self.dim, gamma=self.gamma,
+            hfact=1.2
+        )
+
+        tsph = TSPHScheme(
+            fluids=['fluid'], solids=[], dim=self.dim, gamma=self.gamma,
+            hfact=1.2)
+
         s = SchemeChooser(
-            default='gsph', gsph=gsph, mpm=mpm, crksph=crksph, adke=adke
+            default='gsph', gsph=gsph, mpm=mpm, crksph=crksph, adke=adke,
+            psph=psph, tsph=tsph
         )
 
         return s
@@ -127,20 +139,26 @@ def configure_scheme(self):
                 adaptive_timestep=True, pfreq=50
             )
 
-        if self.options.scheme == 'mpm':
+        elif self.options.scheme == 'mpm':
             s.configure(kernel_factor=1.2)
             s.configure_solver(dt=self.dt, tf=self.tf,
                                adaptive_timestep=False, pfreq=50)
 
-        if self.options.scheme == 'crksph':
+        elif self.options.scheme == 'crksph':
+            s.configure_solver(dt=self.dt, tf=self.tf,
+                               adaptive_timestep=False, pfreq=50)
+
+        elif self.options.scheme == 'adke':
             s.configure_solver(dt=self.dt, tf=self.tf,
                                adaptive_timestep=False, pfreq=50)
 
-        if self.options.scheme == 'adke':
+        elif self.options.scheme in ['tsph', 'psph']:
+            s.configure(hfact=1.2)
             s.configure_solver(dt=self.dt, tf=self.tf,
                                adaptive_timestep=False, pfreq=50)
 
     def post_process(self):
+        import os
         from pysph.solver.utils import load
         if len(self.output_files) < 1:
             return
@@ -169,6 +187,8 @@ def post_process(self):
         )
         l1 = l1 / self.n_particles
         print("l_1 norm of density for the problem: %s" % (l1))
+        fname = os.path.join(self.output_dir, 'norms.npz')
+        numpy.savez(fname, linf_vel=l_inf, l1_vel=l_1, l1_rho=l1)
 
 
 if __name__ == "__main__":

pysph/examples/gas_dynamics/blastwave.py

@@ -1,10 +1,9 @@
 """Simulate a 1D blast wave problem (30 seconds).
 """
-import os
-import numpy
 from pysph.examples.gas_dynamics.shocktube_setup import ShockTubeSetup
 from pysph.sph.scheme import ADKEScheme, GasDScheme, GSPHScheme, SchemeChooser
-
+from pysph.sph.gas_dynamics.psph import PSPHScheme
+from pysph.sph.gas_dynamics.tsph import TSPHScheme
 
 # Numerical constants
 dim = 1
@@ -82,7 +81,17 @@ def create_scheme(self):
             niter=20, tol=1e-6
         )
 
-        s = SchemeChooser(default='adke', adke=adke, gsph=gsph)
+        psph = PSPHScheme(
+            fluids=['fluid'], solids=['boundary'], dim=dim, gamma=gamma,
+            hfact=1.2
+        )
+
+        tsph = TSPHScheme(
+            fluids=['fluid'], solids=['boundary'], dim=dim, gamma=gamma,
+            hfact=1.2)
+
+        s = SchemeChooser(default='adke', adke=adke, gsph=gsph,
+                          psph=psph, tsph=tsph)
         return s
 
 

pysph/examples/gas_dynamics/cheng_shu_1d.py

@@ -22,6 +22,8 @@
 from pysph.base.nnps import DomainManager
 from pysph.solver.application import Application
 from pysph.sph.scheme import GSPHScheme, SchemeChooser
+from pysph.sph.gas_dynamics.psph import PSPHScheme
+from pysph.sph.gas_dynamics.tsph import TSPHScheme
 
 
 class ChengShu(Application):
@@ -80,8 +82,17 @@ def create_scheme(self):
             niter=200, tol=1e-6
         )
 
+        psph = PSPHScheme(
+            fluids=['fluid'], solids=[], dim=self.dim, gamma=self.gamma,
+            hfact=1.2
+        )
+
+        tsph = TSPHScheme(
+            fluids=['fluid'], solids=[], dim=self.dim, gamma=self.gamma,
+            hfact=1.2)
+
         s = SchemeChooser(
-            default='gsph', gsph=gsph
+            default='gsph', gsph=gsph, psph=psph, tsph=tsph
         )
 
         return s
@@ -93,6 +104,10 @@ def configure_scheme(self):
                 dt=self.dt, tf=self.tf,
                 adaptive_timestep=False, pfreq=1000
             )
+        elif self.options.scheme in ['tsph', 'psph']:
+            s.configure(hfact=1.2)
+            s.configure_solver(dt=self.dt, tf=self.tf,
+                               adaptive_timestep=False, pfreq=1000)
 
 
 if __name__ == "__main__":

pysph/examples/gas_dynamics/hydrostatic_box.py

@@ -8,6 +8,8 @@
 from pysph.sph.scheme import (
     ADKEScheme, SchemeChooser, GasDScheme, GSPHScheme
 )
+from pysph.sph.gas_dynamics.psph import PSPHScheme
+from pysph.sph.gas_dynamics.tsph import TSPHScheme
 from pysph.base.utils import get_particle_array as gpa
 from pysph.base.nnps import DomainManager
 from pysph.solver.application import Application
@@ -84,9 +86,19 @@ def create_scheme(self):
         adke = ADKEScheme(
             fluids=['fluid'], solids=[], dim=2, gamma=self.gamma,
             alpha=0.1, beta=0.1, k=1.5, eps=0., g1=0.1, g2=0.1,
-            has_ghosts=True)
+            has_ghosts=True
+        )
+        psph = PSPHScheme(
+            fluids=['fluid'], solids=[], dim=2, gamma=self.gamma,
+            hfact=1.2
+        )
+        tsph = TSPHScheme(
+            fluids=['fluid'], solids=[], dim=2, gamma=self.gamma,
+            hfact=1.2
+        )
         s = SchemeChooser(
-            default='crksph', crksph=crk, adke=adke, mpm=mpm, gsph=gsph
+            default='crksph', crksph=crk, adke=adke, mpm=mpm, gsph=gsph,
+            psph=psph, tsph=tsph
         )
         return s
 
@@ -108,6 +120,10 @@ def configure_scheme(self):
         elif self.options.scheme == 'adke':
             s.configure_solver(dt=self.dt, tf=self.tf,
                                adaptive_timestep=False, pfreq=50)
+        elif self.options.scheme in ['tsph', 'psph']:
+            s.configure(hfact=1.2)
+            s.configure_solver(dt=self.dt, tf=self.tf,
+                               adaptive_timestep=False, pfreq=50)
 
 
 if __name__ == "__main__":

pysph/examples/gas_dynamics/kelvin_helmholtz_instability.py

@@ -8,6 +8,8 @@
 from pysph.base.utils import get_particle_array as gpa
 from pysph.solver.application import Application
 from pysph.sph.scheme import GasDScheme, SchemeChooser, ADKEScheme, GSPHScheme
+from pysph.sph.gas_dynamics.psph import PSPHScheme
+from pysph.sph.gas_dynamics.tsph import TSPHScheme
 from pysph.sph.wc.crksph import CRKSPHScheme
 from pysph.base.nnps import DomainManager
 from pysph.tools import uniform_distribution as ud
@@ -149,9 +151,17 @@ def create_scheme(self):
             interface_zero=True, hybrid=False, blend_alpha=2.0,
             niter=40, tol=1e-6, has_ghosts=True
         )
-
+        psph = PSPHScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            hfact=1.2
+        )
+        tsph = TSPHScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            hfact=1.2
+        )
         s = SchemeChooser(
-            default='crksph', crksph=crk, gsph=gsph, adke=adke, mpm=mpm
+            default='crksph', crksph=crk, gsph=gsph, adke=adke, mpm=mpm,
+            psph=psph, tsph=tsph
         )
 
         return s
@@ -172,6 +182,10 @@ def configure_scheme(self):
         elif self.options.scheme == 'gsph':
             s.configure_solver(dt=dt, tf=tf,
                                adaptive_timestep=False, pfreq=50)
+        elif self.options.scheme in ['tsph', 'psph']:
+            s.configure(hfact=1.2)
+            s.configure_solver(dt=dt, tf=tf,
+                               adaptive_timestep=False, pfreq=50)
 
 
 if __name__ == "__main__":

pysph/examples/gas_dynamics/noh.py

@@ -8,6 +8,8 @@
 from pysph.base.utils import get_particle_array as gpa
 from pysph.solver.application import Application
 from pysph.sph.scheme import GasDScheme, SchemeChooser, ADKEScheme, GSPHScheme
+from pysph.sph.gas_dynamics.psph import PSPHScheme
+from pysph.sph.gas_dynamics.tsph import TSPHScheme
 from pysph.sph.wc.crksph import CRKSPHScheme
 from pysph.base.nnps import DomainManager
 
@@ -104,10 +106,22 @@ def create_scheme(self):
         adke = ADKEScheme(
             fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
             alpha=1, beta=1, k=1.0, eps=0.8, g1=0.5, g2=0.5,
-            has_ghosts=True)
+            has_ghosts=True
+        )
+
+        psph = PSPHScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            hfact=1.2
+        )
+
+        tsph = TSPHScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            hfact=1.2
+        )
 
         s = SchemeChooser(
-            default='crksph', crksph=crksph, mpm=mpm, adke=adke, gsph=gsph
+            default='crksph', crksph=crksph, mpm=mpm, adke=adke, gsph=gsph,
+            psph=psph, tsph=tsph
         )
         s.configure_solver(dt=dt, tf=tf, adaptive_timestep=False)
         return s
@@ -131,6 +145,11 @@ def configure_scheme(self):
             s.configure_solver(
                 dt=dt, tf=tf, adaptive_timestep=False, pfreq=50
             )
+        elif self.options.scheme in ['tsph', 'psph']:
+            s.configure(hfact=1.2)
+            s.configure_solver(
+                dt=dt, tf=tf, adaptive_timestep=False, pfreq=50
+            )
 
     def post_process(self):
         try:
@@ -194,6 +213,12 @@ def post_process(self):
         pyplot.savefig(fname, dpi=300)
         pyplot.close('all')
 
+        fname = os.path.join(self.output_dir, '1dresults.npz')
+        numpy.savez(fname, r=r, rho=rho, p=p)
+
+        fname = os.path.join(self.output_dir, '1dexact.npz')
+        numpy.savez(fname, r=re, rho=rho_e, p=p_e)
+
 
 if __name__ == '__main__':
     app = NohImplosion()