Fix bugs in sod shocktube.

 - Before both solids and mirror bc are set which is incorrect.
 - Solid support is removed as mirror bc is sufficient for this problem.

pysph/examples/gas_dynamics/riemann_solver.py

@@ -130,7 +130,7 @@ def _star_pu(rho_l, u_l, p_l, c_l, rho_r, u_r, p_r, c_r, p_guess):
     fr = _flux_fsolve(p_guess, rho_r, c_r, p_r)
     f = lambda p: fl(p) + fr(p) + u_r - u_l
     from scipy.optimize import fsolve
-    p_star = fsolve(f, 0.0)
+    p_star = fsolve(f, 0.0)[0]
     u_star = (
         0.5 * (u_l + u_r + _flux_fsolve(p_star, rho_r, c_r, p_r)(p_star) -
                _flux_fsolve(p_star, rho_l, c_l, p_l)(p_star))

pysph/examples/gas_dynamics/sod_shocktube.py

@@ -4,7 +4,6 @@
 from pysph.sph.scheme import ADKEScheme, GasDScheme, GSPHScheme, SchemeChooser
 from pysph.sph.wc.crksph import CRKSPHScheme
 from pysph.base.nnps import DomainManager
-import numpy
 
 # Numerical constants
 dim = 1
@@ -50,41 +49,52 @@ def consume_user_options(self):
         self.ml = self.dxl * self.rhol
         self.h0 = self.hdx * self.dxr
         self.hdx = self.hdx
+        self.dt = dt
+        self.tf = tf
 
     def create_particles(self):
-        return self.generate_particles(xmin=self.xmin, xmax=self.xmax,
-                                       dxl=self.dxl, dxr=self.dxr,
-                                       m=self.ml, pl=self.pl, pr=self.pr,
-                                       h0=self.h0, bx=0.03, gamma1=gamma1,
-                                       ul=self.ul, ur=self.ur)
+        # Boundary particles are not needed as we are mirroring the particles
+        # using the domain manager. Hence, bx is set to 0.0.
+        f, b = self.generate_particles(
+            xmin=self.xmin, xmax=self.xmax, dxl=self.dxl, dxr=self.dxr,
+            m=self.ml, pl=self.pl, pr=self.pr, h0=self.h0, bx=0.00,
+            gamma1=gamma1, ul=self.ul, ur=self.ur
+        )
+        self.scheme.setup_properties([f, b])
+        return [f]
 
     def create_domain(self):
         return DomainManager(
-            xmin=self.xmin, xmax=self.xmax, mirror_in_x=True
+            xmin=self.xmin, xmax=self.xmax, mirror_in_x=True,
+            n_layers=2
         )
 
+    def configure_scheme(self):
+        scheme = self.scheme
+        if self.options.scheme in ['gsph', 'mpm']:
+            scheme.configure(kernel_factor=self.hdx)
+        scheme.configure_solver(tf=self.tf, dt=self.dt)
+
     def create_scheme(self):
-        self.dt = dt
-        self.tf = tf
         adke = ADKEScheme(
-            fluids=['fluid'], solids=['boundary'], dim=dim, gamma=gamma,
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
             alpha=1, beta=1.0, k=0.3, eps=0.5, g1=0.2, g2=0.4)
 
         mpm = GasDScheme(
-            fluids=['fluid'], solids=['boundary'], dim=dim, gamma=gamma,
-            kernel_factor=1.2, alpha1=1.0, alpha2=0.1,
-            beta=2.0, update_alpha1=True, update_alpha2=True
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            kernel_factor=None, alpha1=1.0, alpha2=0.1,
+            beta=2.0, update_alpha1=True, update_alpha2=True,
         )
         gsph = GSPHScheme(
-            fluids=['fluid'], solids=['boundary'], dim=dim, gamma=gamma,
-            kernel_factor=1.0,
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            kernel_factor=None,
             g1=0.2, g2=0.4, rsolver=2, interpolation=1, monotonicity=1,
-            interface_zero=True, hybrid=False, blend_alpha=2.0,
+            interface_zero=True, hybrid=True, blend_alpha=2.0,
             niter=20, tol=1e-6
         )
         crk = CRKSPHScheme(
-            fluids=['fluid'], dim=dim, rho0=0, c0=0, nu=0, h0=0, p0=0,
-            gamma=gamma, cl=3
+            fluids=['fluid'], dim=dim, rho0=0, c0=0,
+            nu=0, h0=0, p0=0, gamma=gamma, cl=3
         )
         s = SchemeChooser(
             default='adke', adke=adke, mpm=mpm, gsph=gsph, crk=crk

pysph/sph/scheme.py

@@ -1322,10 +1322,11 @@ def get_equations(self):
             )
         equations.append(Group(equations=group, update_nnps=True))
 
-        group = []
-        for solid in self.solids:
-            group.append(WallBoundary(solid, sources=self.fluids))
-        equations.append(Group(equations=group))
+        if self.solids:
+            group = []
+            for solid in self.solids:
+                group.append(WallBoundary(solid, sources=self.fluids))
+            equations.append(Group(equations=group))
 
         all_pa = self.fluids + self.solids
         group = []
@@ -1337,10 +1338,11 @@ def get_equations(self):
             )
         equations.append(Group(equations=group, update_nnps=False))
 
-        group = []
-        for solid in self.solids:
-            group.append(WallBoundary(solid, sources=self.fluids))
-        equations.append(Group(equations=group))
+        if self.solids:
+            group = []
+            for solid in self.solids:
+                group.append(WallBoundary(solid, sources=self.fluids))
+            equations.append(Group(equations=group))
 
         group = []
         for fluid in self.fluids:
@@ -1368,10 +1370,11 @@ def get_equations(self):
                                      gamma=self.gamma))
         equations.append(Group(equations=group))
 
-        group = []
-        for solid in self.solids:
-            group.append(WallBoundary(solid, sources=self.fluids))
-        equations.append(Group(equations=group))
+        if self.solids:
+            group = []
+            for solid in self.solids:
+                group.append(WallBoundary(solid, sources=self.fluids))
+            equations.append(Group(equations=group))
 
         g2 = []
         for fluid in self.fluids:
@@ -1477,10 +1480,11 @@ def get_equations(self):
         from pysph.sph.gas_dynamics.boundary_equations import WallBoundary
 
         equations = []
-        g1 = []
-        for solid in self.solids:
-            g1.append(WallBoundary(solid, sources=self.fluids))
-        equations.append(Group(equations=g1))
+        if self.solids:
+            g1 = []
+            for solid in self.solids:
+                g1.append(WallBoundary(solid, sources=self.fluids))
+            equations.append(Group(equations=g1))
 
         g2 = []
         for fluid in self.fluids:
@@ -1490,22 +1494,24 @@ def get_equations(self):
                     eps=self.eps
                 )
             )
-        equations.append(Group(g2, update_nnps=True, iterate=False))
+        equations.append(Group(g2, update_nnps=False, iterate=False))
 
-        g3 = []
-        for solid in self.solids:
-            g3.append(WallBoundary(solid, sources=self.fluids))
-        equations.append(Group(equations=g3))
+        if self.solids:
+            g3 = []
+            for solid in self.solids:
+                g3.append(WallBoundary(solid, sources=self.fluids))
+            equations.append(Group(equations=g3))
 
         g4 = []
         for fluid in self.fluids:
             g4.append(SummationDensity(fluid, self.fluids+self.solids))
-        equations.append(Group(g4))
+        equations.append(Group(g4,  update_nnps=True))
 
-        g5 = []
-        for solid in self.solids:
-            g5.append(WallBoundary(solid, sources=self.fluids))
-        equations.append(Group(equations=g5))
+        if self.solids:
+            g5 = []
+            for solid in self.solids:
+                g5.append(WallBoundary(solid, sources=self.fluids))
+            equations.append(Group(equations=g5))
 
         g6 = []
         for elem in self.fluids+self.solids: