Per-Particle artificial speed of sound

spatialpy/Domain.py

@@ -67,6 +67,7 @@ def __init__(self, numpoints, xlim, ylim, zlim, rho0=1.0, c0=10, P0=10, gravity=
         self.mass = numpy.zeros((numpoints), dtype=float)
         self.type = numpy.zeros((numpoints), dtype=int)
         self.nu = numpy.zeros((numpoints), dtype=float)
+        self.c = numpy.zeros((numpoints), dtype=float)
         self.rho = numpy.zeros((numpoints), dtype=float)
         self.fixed = numpy.zeros((numpoints), dtype=bool)
 
@@ -88,6 +89,8 @@ def __str__(self):
         domain_strs.extend(["", "Paritcles", ""])
         for i, vertex in enumerate(self.vertices):
             v_str = f"{pad}{i+1}: {vertex}\n{pad}   Volume:{self.vol[i]}, Mass: {self.mass[i]}, "
+            v_str += f"Type: {self.type[i]}, Viscosity: {self.nu[i]}, Density: {self.rho[i]}, "
+            v_str += f"Artificial Speed of Sound: {self.c[i]}, Fixed: {self.fixed[i]}"
             v_str += f"Type: {self.type[i]}, Viscosity: {self.nu[i]}, Density: {self.rho[i]}, Fixed: {self.fixed[i]}"
             domain_strs.append(v_str)
         if self.triangles is not None:
@@ -104,7 +107,7 @@ def __str__(self):
     def _ipython_display_(self, use_matplotlib=False):
         self.plot_types(width="auto", height="auto", use_matplotlib=use_matplotlib)
 
-    def add_point(self, x, vol, mass, type, nu, fixed, rho=None):
+    def add_point(self, x, vol, mass, type, nu, fixed, rho=None, c=0):
         """ Add a single point particle to the domain space.
 
             :param x: Spatial coordinate vertices of point to be added
@@ -125,6 +128,9 @@ def add_point(self, x, vol, mass, type, nu, fixed, rho=None):
             :param fixed: True if particle is spatially fixed, else False
             :type fixed: bool
 
+            :param c: Default artificial speed of sound of particle to be created
+            :type c: float
+
             :param rho: Default density of particle to be created
             :type rho: float
         """
@@ -139,6 +145,7 @@ def add_point(self, x, vol, mass, type, nu, fixed, rho=None):
         self.mass = numpy.append(self.mass, mass)
         self.type = numpy.append(self.type, type)
         self.nu = numpy.append(self.nu, nu)
+        self.c = numpy.append(self.c, c)
         self.rho = numpy.append(self.rho, rho)
         self.fixed = numpy.append(self.fixed, fixed)
 
@@ -569,16 +576,18 @@ def read_stochss_domain(cls, filename):
             for particle in domain['particles']:
                 # StochSS backward compatability check for rho
                 rho = None if "rho" not in particle.keys() else particle['rho']
+                # StochSS backward compatability check for c
+                c = 0 if "c" not in particle.keys() else particle['c']
                 obj.add_point(particle['point'], particle['volume'], particle['mass'],
-                              particle['type'], particle['nu'], particle['fixed'], rho=rho)
+                              particle['type'], particle['nu'], particle['fixed'], rho=rho, c=c)
 
             return obj
         except KeyError as e:
             raise DomainError("The file is not a StochSS Domain (.domn) or a StochSS Spatial Model (.smdl).")
 
 
     @classmethod
-    def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=1.0, rho=None, fixed=False, **kwargs):
+    def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=1.0, rho=None, c=0, fixed=False, **kwargs):
         """ Create a filled 3D domain
 
             :param xlim: highest and lowest coordinate in the x dimension
@@ -608,6 +617,9 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=
             :param nu: default viscosity of particles to be created. Defaults to 1.0
             :type nu: float
 
+            :param c: default artificial speed of sound of particles to be created. Defaults to 0.0.
+            :type c: float
+
             :param rho: default density of particles to be created.
             :type rho:
 
@@ -650,6 +662,7 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=
                     obj.type[ndx] = type_id
                     obj.mass[ndx] = mass
                     obj.nu[ndx] = nu
+                    obj.c[ndx] = c
                     obj.rho[ndx] = rho
                     obj.fixed[ndx] = fixed
                     ndx+=1
@@ -658,7 +671,7 @@ def create_3D_domain(cls, xlim, ylim, zlim, nx, ny, nz, type_id=1, mass=1.0, nu=
         return obj
 
     @classmethod
-    def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, rho=None, fixed=False, **kwargs):
+    def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, rho=None, c=0, fixed=False, **kwargs):
         """ Create a filled 2D domain
 
             :param xlim: highest and lowest coordinate in the x dimension
@@ -682,6 +695,9 @@ def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, rho=N
             :param nu: default viscosity of particles to be created. Defaults to 1.0
             :type nu: float
 
+            :param c: default artificial speed of sound of particles to be created. Defaults to 0.0.
+            :type c: float
+
             :param rho: default density of particles to be created.
             :type rho:
 
@@ -722,6 +738,7 @@ def create_2D_domain(cls, xlim, ylim, nx, ny, type_id=1, mass=1.0, nu=1.0, rho=N
                 obj.type[ndx] = type_id
                 obj.mass[ndx] = mass
                 obj.nu[ndx] = nu
+                obj.c[ndx] = c
                 obj.rho[ndx] = rho
                 obj.fixed[ndx] = fixed
                 ndx+=1

spatialpy/Model.py

@@ -237,7 +237,7 @@ def timespan(self, time_span, timestep_size=None):
 
 
 
-    def set_type(self, geometry_ivar, type_id, mass=None, nu=None, rho=None, fixed=False):
+    def set_type(self, geometry_ivar, type_id, mass=None, nu=None, rho=None, c=None, fixed=False):
         """ Add a type definition to the model.  By default, all regions are set to
         type 0. Returns the number of domain points that were tagged with this type_id
 
@@ -252,6 +252,8 @@ def set_type(self, geometry_ivar, type_id, mass=None, nu=None, rho=None, fixed=F
             :type rho: float
             :param nu: The viscosity of each particle in the type
             :type nu: float
+            :param c: The artificial speed of sound of each particle in the type
+            :type c: float
             :param fixed: Are the particles in this type immobile
             :type fixed: bool
 
@@ -275,6 +277,8 @@ def set_type(self, geometry_ivar, type_id, mass=None, nu=None, rho=None, fixed=F
                     self.domain.rho[v_ndx] = rho
                 if nu is not None:
                     self.domain.nu[v_ndx] = nu
+                if c is not None:
+                    self.domain.c[v_ndx] = c
                 self.domain.fixed[v_ndx] = fixed
                 count +=1
         if count == 0:

spatialpy/Solver.py

@@ -430,9 +430,9 @@ def __create_propensity_file(self, file_name=None):
             if self.model.domain.type[i] == 0:
                 raise SimulationError(
                     "Not all particles have been defined in a type. Mass and other properties must be defined")
-            init_particles += "init_create_particle(sys,id++,{0},{1},{2},{3},{4},{5},{6},{7},{8});".format(
+            init_particles += "init_create_particle(sys,id++,{0},{1},{2},{3},{4},{5},{6},{7},{8},{9});".format(
                 self.model.domain.coordinates()[i,0],self.model.domain.coordinates()[i,1],self.model.domain.coordinates()[i,2],
-                self.model.domain.type[i],self.model.domain.nu[i],self.model.domain.mass[i],
+                self.model.domain.type[i],self.model.domain.nu[i],self.model.domain.mass[i],self.model.domain.c[i],
                 self.model.domain.rho[i],int(self.model.domain.fixed[i]),num_chem_species )+"\n"
         propfilestr = propfilestr.replace("__INIT_PARTICLES__", init_particles)
 

spatialpy/ssa_sdpd-c-simulation-engine/include/particle.hpp

@@ -42,7 +42,7 @@ namespace Spatialpy{
 
     struct Particle{
         Particle(ParticleSystem *sys, unsigned int id=0, double x=0, double y=0, double z=0,
-                    int type=0, double nu=0.01, double mass=1,
+                    int type=0, double nu=0.01, double mass=1, double c=0,
                     double rho=1, int solidTag=0);
         ParticleSystem *sys;
         unsigned int id;
@@ -56,6 +56,7 @@ namespace Spatialpy{
         double mass;
         double rho;
         double nu;
+        double c;
         int solidTag;
         double bvf_phi;
         double normal[3];

spatialpy/ssa_sdpd-c-simulation-engine/propensity_file_template.cpp

@@ -70,8 +70,8 @@ namespace Spatialpy{
 
     //dsfmt_t dsfmt;
 
-    void init_create_particle(ParticleSystem *sys, unsigned int id, double x, double y, double z, int type, double nu, double mass, double rho, int solidTag, int num_chem_species){
-        sys->particles.emplace_back(Particle(sys, id, x, y, z, type, nu, mass, rho,
+    void init_create_particle(ParticleSystem *sys, unsigned int id, double x, double y, double z, int type, double nu, double mass, double c, double rho, int solidTag, int num_chem_species){
+        sys->particles.emplace_back(Particle(sys, id, x, y, z, type, nu, mass, c, rho,
             solidTag));
         Particle *this_particle = &sys->particles.back() ;
         if(num_chem_species > 0){

spatialpy/ssa_sdpd-c-simulation-engine/src/particle.cpp

@@ -66,9 +66,9 @@ namespace Spatialpy{
 
     Particle::Particle(ParticleSystem *sys, unsigned int id, 
                         double xl, double yl, double zl, int type, double nu, 
-                        double mass, double rho, int solidTag) : 
-                        sys(sys), id(id), type(type), 
-                        nu(nu), mass(mass), rho(rho), solidTag(solidTag)
+                        double mass, double c, double rho, int solidTag) : 
+                        sys(sys), id(id), type(type), nu(nu),
+                        c(c), mass(mass), rho(rho), solidTag(solidTag)
     {
     	x[0] = xl ;
         x[1] = yl ;