Merge pull request #225 from ananyo-work/crksph

CRKSPH scheme, MPM and mirrorbc fix

pysph/base/nnps_base.pyx

@@ -613,15 +613,17 @@ cdef class CPUDomainManager(DomainManagerBase):
             if mirror_in_x:
                 # x_low
                 copy = pa.extract_particles( x_low )
-                self._add_array_to_array(copy.get_carray('x'), xt_low)
-                self._mul_to_array(copy.get_carray('u'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('x'), xt_low)
+                    self._mul_to_array(copy.get_carray('u'), -1)
+                    added.append_parray(copy)
 
                 # x_high
                 copy = pa.extract_particles( x_high )
-                self._add_array_to_array(copy.get_carray('x'), xt_high)
-                self._mul_to_array(copy.get_carray('u'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('x'), xt_high)
+                    self._mul_to_array(copy.get_carray('u'), -1)
+                    added.append_parray(copy)
 
             if mirror_in_y:
                 # Now do the corners from the previous.
@@ -638,27 +640,31 @@ cdef class CPUDomainManager(DomainManagerBase):
                         high_translate.append(2*(ymax - yi))
 
                 copy = added.extract_particles(low)
-                self._add_array_to_array(copy.get_carray('y'), low_translate)
-                self._mul_to_array(copy.get_carray('v'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('y'), low_translate)
+                    self._mul_to_array(copy.get_carray('v'), -1)
+                    added.append_parray(copy)
 
                 copy = added.extract_particles(high)
-                self._add_array_to_array(copy.get_carray('y'), high_translate)
-                self._mul_to_array(copy.get_carray('v'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('y'), high_translate)
+                    self._mul_to_array(copy.get_carray('v'), -1)
+                    added.append_parray(copy)
 
                 # Add the actual y_high and y_low now.
                 # y_high
                 copy = pa.extract_particles( y_high )
-                self._add_array_to_array(copy.get_carray('y'), yt_high)
-                self._mul_to_array(copy.get_carray('v'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('y'), yt_high)
+                    self._mul_to_array(copy.get_carray('v'), -1)
+                    added.append_parray(copy)
 
                 # y_low
                 copy = pa.extract_particles( y_low )
-                self._add_array_to_array(copy.get_carray('y'), yt_low)
-                self._mul_to_array(copy.get_carray('v'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('y'), yt_low)
+                    self._mul_to_array(copy.get_carray('v'), -1)
+                    added.append_parray(copy)
 
             if mirror_in_z:
                 # Now do the corners from the previous.
@@ -675,27 +681,31 @@ cdef class CPUDomainManager(DomainManagerBase):
                         high_translate.append(2*(zmax - zi))
 
                 copy = added.extract_particles(low)
-                self._add_array_to_array(copy.get_carray('z'), low_translate)
-                self._mul_to_array(copy.get_carray('w'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('z'), low_translate)
+                    self._mul_to_array(copy.get_carray('w'), -1)
+                    added.append_parray(copy)
 
                 copy = added.extract_particles(high)
-                self._add_array_to_array(copy.get_carray('z'), high_translate)
-                self._mul_to_array(copy.get_carray('w'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('z'), high_translate)
+                    self._mul_to_array(copy.get_carray('w'), -1)
+                    added.append_parray(copy)
 
                 # Add the actual z_high and z_low now.
                 # z_high
                 copy = pa.extract_particles( z_high )
-                self._add_array_to_array(copy.get_carray('z'), zt_high)
-                self._mul_to_array(copy.get_carray('w'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('z'), zt_high)
+                    self._mul_to_array(copy.get_carray('w'), -1)
+                    added.append_parray(copy)
 
                 # z_low
                 copy = pa.extract_particles( z_low )
-                self._add_array_to_array(copy.get_carray('z'), zt_low)
-                self._mul_to_array(copy.get_carray('w'), -1)
-                added.append_parray(copy)
+                if copy.get_number_of_particles() > 0:
+                    self._add_array_to_array(copy.get_carray('z'), zt_low)
+                    self._mul_to_array(copy.get_carray('w'), -1)
+                    added.append_parray(copy)
 
 
             added.tag[:] = Ghost

pysph/examples/gas_dynamics/accuracy_test_2d.py

@@ -12,6 +12,7 @@
 from pysph.solver.application import Application
 
 from pysph.sph.scheme import GasDScheme, ADKEScheme, GSPHScheme, SchemeChooser
+from pysph.sph.wc.crksph import CRKSPHScheme
 
 # PySPH tools
 from pysph.tools import uniform_distribution as ud
@@ -46,13 +47,27 @@ def initialize(self):
         self.xmax = xmax
         self.ymin = ymin
         self.ymax = ymax
-        self.ny = 100
+        self.ny = 128
         self.nx = self.ny
         self.dx = (self.xmax - self.xmin) / (self.nx)
         self.hdx = 2.
         self.p = 1.
         self.u = 1
         self.v = -1
+        self.c_0 = 1.18
+        self.cfl = 0.1
+
+    def add_user_options(self, group):
+        group.add_argument(
+            "--nparticles", action="store", type=int, dest="nprt", default=256,
+            help="Number of particles in domain"
+        )
+
+    def consume_user_options(self):
+        self.nx = self.options.nprt
+        self.ny = self.nx
+        self.dx = (self.xmax - self.xmin) / (self.nx)
+        self.dt = self.cfl * self.dx / self.c_0
 
     def create_domain(self):
         return DomainManager(
@@ -100,29 +115,33 @@ def create_particles(self):
         return [fluid, ]
 
     def create_scheme(self):
-        self.dt = dt
         self.tf = tf
-
         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)
+            alpha=0, beta=0, k=1.5, eps=0., g1=0., g2=0.,
+            has_ghosts=True)
 
         mpm = GasDScheme(
             fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
-            kernel_factor=kernel_factor, alpha1=alpha1, alpha2=alpha2,
-            beta=beta
+            kernel_factor=kernel_factor, alpha1=0, alpha2=0,
+            beta=beta, has_ghosts=True
+        )
+
+        crksph = CRKSPHScheme(
+            fluids=['fluid'], dim=dim, rho0=0, c0=0, nu=0, h0=0, p0=0,
+            gamma=gamma, cl=2, has_ghosts=True
         )
 
         gsph = GSPHScheme(
             fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
             kernel_factor=1.,
-            g1=0., g2=0., rsolver=2, interpolation=0, monotonicity=0,
+            g1=0., g2=0., rsolver=7, interpolation=1, monotonicity=1,
             interface_zero=True, hybrid=False, blend_alpha=5.0,
-            niter=40, tol=1e-6
+            niter=40, tol=1e-6, has_ghosts=True
         )
 
         s = SchemeChooser(
-            default='gsph', adke=adke, mpm=mpm, gsph=gsph
+            default='gsph', adke=adke, mpm=mpm, gsph=gsph, crksph=crksph
         )
         return s
 
@@ -138,6 +157,9 @@ def configure_scheme(self):
         elif self.options.scheme == 'gsph':
             s.configure_solver(dt=self.dt, tf=self.tf,
                                adaptive_timestep=False, pfreq=50)
+        elif self.options.scheme == 'crksph':
+            s.configure_solver(dt=self.dt, tf=self.tf,
+                               adaptive_timestep=False, pfreq=50)
 
     def post_process(self):
         from pysph.solver.utils import load

pysph/examples/gas_dynamics/acoustic_wave.py

@@ -1,4 +1,4 @@
-r"""Diffusion of an acoustic wave in 1-d (30 secs)
+r"""Diffusion of an acoustic wave in 1-d (5 minutes)
 
 Propagation of acoustic wave
 particles have properties according
@@ -24,6 +24,7 @@
 from pysph.solver.application import Application
 from pysph.sph.scheme import \
     GSPHScheme, ADKEScheme, GasDScheme, SchemeChooser
+from pysph.sph.wc.crksph import CRKSPHScheme
 
 
 class AcousticWave(Application):
@@ -35,20 +36,31 @@ def initialize(self):
         self.p_0 = 1.
         self.c_0 = 1.
         self.delta_rho = 1e-6
-        self.n_particles = 400
+        self.n_particles = 8
         self.domain_length = self.xmax - self.xmin
-        self.dx = self.domain_length / (self.n_particles)
         self.k = -2 * numpy.pi / self.domain_length
-        self.hdx = 1.
+        self.cfl = 0.1
+        self.hdx = 1.0
         self.dt = 1e-3
-        self.tf = 10
+        self.tf = 5
         self.dim = 1
 
     def create_domain(self):
         return DomainManager(
             xmin=0, xmax=1, periodic_in_x=True
         )
 
+    def add_user_options(self, group):
+        group.add_argument(
+            "--nparticles", action="store", type=int, dest="nprt", default=256,
+            help="Number of particles in domain"
+        )
+
+    def consume_user_options(self):
+        self.n_particles = self.options.nprt
+        self.dx = self.domain_length / (self.n_particles)
+        self.dt = self.cfl * self.dx / self.c_0
+
     def create_particles(self):
         x = numpy.arange(
             self.xmin + self.dx*0.5, self.xmax, self.dx
@@ -68,7 +80,8 @@ def create_particles(self):
         m = numpy.ones_like(x) * self.dx * rho
         e = p / ((self.gamma - 1) * rho)
         fluid = gpa(
-            name='fluid', x=x, p=p, rho=rho, u=u, h=h, m=m, e=e, cs=cs
+            name='fluid', x=x, p=p, rho=rho, u=u, h=h, m=m, e=e, cs=cs,
+            h0=h.copy()
         )
         self.scheme.setup_properties([fluid])
 
@@ -80,17 +93,28 @@ def create_scheme(self):
             gamma=self.gamma, kernel_factor=1.0,
             g1=0., g2=0., rsolver=7, interpolation=1, monotonicity=1,
             interface_zero=True, hybrid=False, blend_alpha=5.0,
-            niter=40, tol=1e-6
+            niter=40, tol=1e-6, has_ghosts=True
         )
 
         mpm = GasDScheme(
             fluids=['fluid'], solids=[], dim=self.dim, gamma=self.gamma,
             kernel_factor=1.2, alpha1=0, alpha2=0,
-            beta=2.0, update_alpha1=False, update_alpha2=False
+            beta=2.0, update_alpha1=False, update_alpha2=False,
+            has_ghosts=True
         )
 
+        crksph = CRKSPHScheme(
+            fluids=['fluid'], dim=self.dim, rho0=0, c0=0, nu=0, h0=0, p0=0,
+            gamma=self.gamma, cl=2, has_ghosts=True
+        )
+
+        adke = ADKEScheme(
+            fluids=['fluid'], solids=[], dim=self.dim, gamma=self.gamma,
+            alpha=0, beta=0.0, k=1.5, eps=0.0, g1=0.0, g2=0.0,
+            has_ghosts=True)
+
         s = SchemeChooser(
-            default='gsph', gsph=gsph, mpm=mpm
+            default='gsph', gsph=gsph, mpm=mpm, crksph=crksph, adke=adke
         )
 
         return s
@@ -106,7 +130,15 @@ def configure_scheme(self):
         if self.options.scheme == 'mpm':
             s.configure(kernel_factor=1.2)
             s.configure_solver(dt=self.dt, tf=self.tf,
-                               adaptive_timestep=True, pfreq=50)
+                               adaptive_timestep=False, pfreq=50)
+
+        if self.options.scheme == 'crksph':
+            s.configure_solver(dt=self.dt, tf=self.tf,
+                               adaptive_timestep=False, pfreq=50)
+
+        if self.options.scheme == 'adke':
+            s.configure_solver(dt=self.dt, tf=self.tf,
+                               adaptive_timestep=False, pfreq=50)
 
     def post_process(self):
         from pysph.solver.utils import load
@@ -122,8 +154,21 @@ def post_process(self):
         l_inf = numpy.max(
             numpy.abs(u_c - u)
         )
+        l_1 = (numpy.sum(
+            numpy.abs(u_c - u)
+        ) / self.n_particles)
+        print("L_inf norm of velocity for the problem: %s" % (l_inf))
+        print("L_1 norm of velocity for the problem: %s" % (l_1))
+
+        rho = self.rho_0 + self.delta_rho *\
+            numpy.sin(self.k * x_c)
 
-        print("L_inf norm for the problem: %s" % (l_inf))
+        rho_c = pa.rho
+        l1 = numpy.sum(
+            numpy.abs(rho - rho_c)
+        )
+        l1 = l1 / self.n_particles
+        print("l_1 norm of density for the problem: %s" % (l1))
 
 
 if __name__ == "__main__":