Merge pull request #207 from ananyo-work/gas-dynamics-fixes

Gas dynamics fixes

pysph/examples/gas_dynamics/accuracy_test_2d.py

@@ -0,0 +1,166 @@
+"""Acoustic wave diffusion in 2-d (2 mins)
+
+Two Dimensional constant pressure accuracy test
+particles should simply advect in a periodic domain
+"""
+
+# NumPy and standard library imports
+import numpy
+
+from pysph.base.nnps import DomainManager
+from pysph.base.utils import get_particle_array as gpa
+from pysph.solver.application import Application
+
+from pysph.sph.scheme import GasDScheme, ADKEScheme, GSPHScheme, SchemeChooser
+
+# PySPH tools
+from pysph.tools import uniform_distribution as ud
+
+# Numerical constants
+dim = 2
+gamma = 1.4
+gamma1 = gamma - 1.0
+
+# solution parameters
+dt = 5e-3
+tf = 1.
+
+
+# domain size
+xmin = 0.
+xmax = 1.
+ymin = 0.
+ymax = 1.
+
+
+# scheme constants
+alpha1 = 1.0
+alpha2 = 0.1
+beta = 2.0
+kernel_factor = 1.5
+
+
+class AccuracyTest2D(Application):
+    def initialize(self):
+        self.xmin = xmin
+        self.xmax = xmax
+        self.ymin = ymin
+        self.ymax = ymax
+        self.ny = 100
+        self.nx = self.ny
+        self.dx = (self.xmax - self.xmin) / (self.nx)
+        self.hdx = 2.
+        self.p = 1.
+        self.u = 1
+        self.v = -1
+
+    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
+        data = ud.uniform_distribution_cubic2D(
+            self.dx, xmin, xmax, ymin, ymax
+            )
+
+        x = data[0].ravel()
+        y = data[1].ravel()
+        dx = data[2]
+
+        volume = dx * dx
+
+        rho = 1 + 0.2 * numpy.sin(
+            numpy.pi * (x + y)
+        )
+
+        p = numpy.ones_like(x) * self.p
+
+        # const h and mass
+        h = numpy.ones_like(x) * self.hdx * dx
+        m = numpy.ones_like(x) * volume * rho
+
+        # u = 1
+        u = numpy.ones_like(x) * self.u
+
+        # v = -1
+        v = numpy.ones_like(x) * self.v
+
+        # thermal energy from the ideal gas EOS
+        e = p/(gamma1*rho)
+
+        fluid = gpa(name='fluid', x=x, y=y, rho=rho, p=p, e=e, h=h, m=m,
+                    h0=h.copy(), u=u, v=v)
+        self.scheme.setup_properties([fluid])
+
+        print("2D Accuracy Test with %d particles"
+              % (fluid.get_number_of_particles()))
+
+        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)
+
+        mpm = GasDScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            kernel_factor=kernel_factor, alpha1=alpha1, alpha2=alpha2,
+            beta=beta
+        )
+
+        gsph = GSPHScheme(
+            fluids=['fluid'], solids=[], dim=dim, gamma=gamma,
+            kernel_factor=1.,
+            g1=0., g2=0., rsolver=2, interpolation=0, monotonicity=0,
+            interface_zero=True, hybrid=False, blend_alpha=5.0,
+            niter=40, tol=1e-6
+        )
+
+        s = SchemeChooser(
+            default='gsph', adke=adke, mpm=mpm, gsph=gsph
+        )
+        return s
+
+    def configure_scheme(self):
+        s = self.scheme
+        if self.options.scheme == 'mpm':
+            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)
+
+    def post_process(self):
+        from pysph.solver.utils import load
+        if len(self.output_files) < 1:
+            return
+        outfile = self.output_files[-1]
+        data = load(outfile)
+        pa = data['arrays']['fluid']
+        x_c = pa.x
+        y_c = pa.y
+        rho_c = pa.rho
+        rho_e = 1 + 0.2 * numpy.sin(
+            numpy.pi * (x_c + y_c)
+        )
+        num_particles = rho_c.size
+        l1_norm = numpy.sum(
+            numpy.abs(rho_c - rho_e)
+        ) / num_particles
+
+        print(l1_norm)
+
+
+if __name__ == '__main__':
+    app = AccuracyTest2D()
+    app.run()
+    app.post_process()

pysph/examples/gas_dynamics/acoustic_wave.py

@@ -0,0 +1,130 @@
+"""Diffusion of an acoustic wave in 1-d (30 secs)
+
+Propagation of acoustic wave
+particles have properties according
+to the following distribuion
+\rho = \rho_0 + \Delta\rho sin(kx)
+p = p_0 + c_0^2\Delta\rho sin(kx)
+u = c_0\rho_0^{-1}\Delta\rho sin(kx)
+
+with \Delta\rho = 1e-6 and k = 2\pi/\lambda
+where \lambda is the domain length.
+\rho_0 = \gamma = 1.4 and p_0 = 1.0
+"""
+
+
+# standard library and numpy imports
+import numpy
+
+# pysph imports
+from pysph.base.utils import get_particle_array as gpa
+from pysph.base.nnps import DomainManager
+from pysph.solver.application import Application
+from pysph.sph.scheme import \
+    GSPHScheme, ADKEScheme, GasDScheme, SchemeChooser
+
+
+class AcousticWave(Application):
+    def initialize(self):
+        self.xmin = 0.
+        self.xmax = 1.
+        self.gamma = 1.4
+        self.rho_0 = self.gamma
+        self.p_0 = 1.
+        self.c_0 = 1.
+        self.delta_rho = 1e-6
+        self.n_particles = 400
+        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.dt = 1e-3
+        self.tf = 10
+        self.dim = 1
+
+    def create_domain(self):
+        return DomainManager(
+            xmin=0, xmax=1, periodic_in_x=True
+        )
+
+    def create_particles(self):
+        x = numpy.arange(
+            self.xmin + self.dx*0.5, self.xmax, self.dx
+        )
+        rho = self.rho_0 + self.delta_rho *\
+            numpy.sin(self.k * x)
+
+        p = self.p_0 + self.c_0**2 *\
+            self.delta_rho * numpy.sin(self.k * x)
+
+        u = self.c_0 * self.delta_rho * numpy.sin(self.k * x) /\
+            self.rho_0
+        cs = numpy.sqrt(
+            self.gamma * p / rho
+        )
+        h = numpy.ones_like(x) * self.dx * self.hdx
+        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
+        )
+        self.scheme.setup_properties([fluid])
+
+        return [fluid, ]
+
+    def create_scheme(self):
+        gsph = GSPHScheme(
+            fluids=['fluid'], solids=[], dim=self.dim,
+            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
+        )
+
+        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
+        )
+
+        s = SchemeChooser(
+            default='gsph', gsph=gsph, mpm=mpm
+        )
+
+        return s
+
+    def configure_scheme(self):
+        s = self.scheme
+        if self.options.scheme == 'gsph':
+            s.configure_solver(
+                dt=self.dt, tf=self.tf,
+                adaptive_timestep=True, pfreq=50
+            )
+
+        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)
+
+    def post_process(self):
+        from pysph.solver.utils import load
+        if len(self.output_files) < 1:
+            return
+        outfile = self.output_files[-1]
+        data = load(outfile)
+        pa = data['arrays']['fluid']
+        x_c = pa.x
+        u = self.c_0 * self.delta_rho * numpy.sin(self.k * x_c) /\
+            self.rho_0
+        u_c = pa.u
+        l_inf = numpy.max(
+            numpy.abs(u_c - u)
+        )
+
+        print("L_inf norm for the problem: %s" % (l_inf))
+
+
+if __name__ == "__main__":
+    app = AcousticWave()
+    app.run()
+    app.post_process()

pysph/examples/gas_dynamics/cheng_shu_1d.py

@@ -0,0 +1,99 @@
+"""Cheng and Shu's 1d acoustic wave propagation in 1d (1 min)
+
+particles have properties according
+to the following distribuion
+\rho = \rho_0 + \Delta\rho sin(kx)
+p = 1.0
+u = 1 + 0.1sin(kx)
+
+with \Delta\rho = 1 and k = 2\pi/\lambda
+where \lambda is the domain length.
+\rho_0 = 2, \gamma = 1.4 and p_0 = 1.0
+"""
+
+
+# standard library and numpy imports
+import numpy
+
+# pysph imports
+from pysph.base.utils import get_particle_array as gpa
+from pysph.base.nnps import DomainManager
+from pysph.solver.application import Application
+from pysph.sph.scheme import \
+    GSPHScheme, ADKEScheme, GasDScheme, SchemeChooser
+
+
+class ChengShu(Application):
+    def initialize(self):
+        self.xmin = 0.
+        self.xmax = 1.
+        self.gamma = 1.4
+        self.p_0 = 1.
+        self.c_0 = 1.
+        self.delta_rho = 1
+        self.n_particles = 1000
+        self.domain_length = self.xmax - self.xmin
+        self.dx = self.domain_length / (self.n_particles - 1)
+        self.k = 2 * numpy.pi / self.domain_length
+        self.hdx = 2.
+        self.dt = 1e-4
+        self.tf = 1.0
+        self.dim = 1
+
+    def create_domain(self):
+        return DomainManager(
+            xmin=self.xmin, xmax=self.xmax, periodic_in_x=True
+        )
+
+    def create_particles(self):
+        x = numpy.linspace(
+            self.xmin, self.xmax, self.n_particles
+            )
+        rho = 2 + numpy.sin(2 * numpy.pi * x)*self.delta_rho
+
+        p = numpy.ones_like(x)
+
+        u = 1 + 0.1 * numpy.sin(2 * numpy.pi * x)
+
+        cs = numpy.sqrt(
+            self.gamma * p / rho
+        )
+        h = numpy.ones_like(x) * self.dx * self.hdx
+        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
+        )
+
+        self.scheme.setup_properties([fluid])
+
+        return [fluid, ]
+
+    def create_scheme(self):
+        gsph = GSPHScheme(
+            fluids=['fluid'], solids=[], dim=self.dim,
+            gamma=self.gamma, kernel_factor=1.,
+            g1=0., g2=0., rsolver=3, interpolation=1, monotonicity=1,
+            interface_zero=True, hybrid=False, blend_alpha=5.0,
+            niter=200, tol=1e-6
+        )
+
+        s = SchemeChooser(
+            default='gsph', gsph=gsph
+        )
+
+        return s
+
+    def configure_scheme(self):
+        s = self.scheme
+        if self.options.scheme == 'gsph':
+            s.configure_solver(
+                dt=self.dt, tf=self.tf,
+                adaptive_timestep=False, pfreq=1000
+            )
+
+
+if __name__ == "__main__":
+    app = ChengShu()
+    app.run()