Fix/pressure

spatialpy/core/Domain.py

@@ -53,7 +53,7 @@ class Domain():
     """
 
 
-    def __init__(self, numpoints, xlim, ylim, zlim, rho0=1.0, c0=10, P0=10, gravity=None):
+    def __init__(self, numpoints, xlim, ylim, zlim, rho0=1.0, c0=10, P0=None, gravity=None):
         self.vertices = numpy.zeros((numpoints, 3), dtype=float)
         self.triangles = None
         self.tetrahedrons = None
@@ -73,7 +73,10 @@ def __init__(self, numpoints, xlim, ylim, zlim, rho0=1.0, c0=10, P0=10, gravity=
 
         self.rho0 = rho0
         self.c0 = c0
-        self.P0 = P0
+        if P0 != None:
+            self.P0 = P0
+        else:
+            self.P0 = rho0 * c0**2 # assume gamma of 1, full equation is rho0*c0**2/gamma
         self.gravity = gravity
 
         self.xlim = xlim

spatialpy/core/Model.py

@@ -464,10 +464,11 @@ def set_timesteps(self, output_interval, num_steps, timestep_size=None):
         # array of step numbers corresponding to the simulation times in the timespan
         output_steps = numpy.arange(0, self.num_timesteps + self.timestep_size, self.output_freq)
         self.output_steps = numpy.unique(numpy.round(output_steps).astype(int))
-        sim_steps = numpy.arange(0, self.num_timesteps + self.timestep_size, self.timestep_size)
         self.tspan = numpy.zeros((self.output_steps.size), dtype=float)
         for i, step in enumerate(self.output_steps):
-            self.tspan[i] = sim_steps[step]
+            #TODO: review conflict resolution here
+            self.tspan[i] = step*self.timestep_size  # changed in this branch
+            #self.tspan[i] = sim_steps[step]         # changed in other branch
 
     def timespan(self, time_span, timestep_size=None):
         """

spatialpy/solvers/Solver.py

@@ -24,7 +24,6 @@
 import tempfile
 import threading
 import time
-import tempfile
 import getpass
 import re
 
@@ -97,7 +96,7 @@ def __del__(self):
         try:
             if self.build_dir is not None:
                 try:
-                    shutil.rmtree(self.build_dir)
+                    shutil.rmtree(self.build_dir, ignore_errors=True)
                 except OSError as e:
                     print("Could not delete '{0}'".format(
                         self.solver_base_dir))

spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/model.c

@@ -169,6 +169,7 @@ namespace Spatialpy{
             //fflush(stdout);
         }
         //printf("pairwiseForce(id=%i) num_chem_rxns=%i\n",me->id,system->num_chem_rxns);
+        //printf("pairwiseForce(id=%i) F=[%e %e %e] Fbp=[%e %e %e]\n",me->id,me->F[0],me->F[1],me->F[2],me->Fbp[0],me->Fbp[1],me->Fbp[2]);
         //fflush(stdout);
 
         // after processing all neighbors
@@ -177,9 +178,6 @@ namespace Spatialpy{
         double vol = (me->mass / me->rho);
         double cur_time = system->current_step * system->dt;
         for(rxn=0; rxn < system->num_chem_rxns; rxn++){
-            //TODO: t (2nd arg) set to zero, fix
-            //TODO, vol (3rd arg) set to zero, fix
-            //TODO: data (4th arg) set to NULL, fix
             double flux = (*system->chem_rxn_rhs_functions[rxn])(me->C, cur_time, vol , me->data_fn, me->type);
             for(s=0; s< system->num_chem_species; s++){
                 int k = system->num_chem_rxns * rxn + s;
@@ -220,7 +218,7 @@ namespace Spatialpy{
           Wij = alpha*( (1+3*R) * pow( 1-R,3));
 
           //Compute numerator of Shepard filter
-          num += pt_j->rho * Wij;
+          num += pt_j->old_rho * Wij;
 
           // Compute denominator of Shepard filter
           den += Wij;

spatialpy/solvers/c_base/ssa_sdpd-c-simulation-engine/src/simulate.cpp

@@ -102,6 +102,9 @@ namespace Spatialpy{
             me->Q[i] = 0.0;
         }
 
+        // save the current per-particle rho for the Shepard filter
+        me->old_rho = me->rho;
+
 
     }
 
@@ -138,17 +141,17 @@ namespace Spatialpy{
             }
 
             // Update density using continuity equation and Shepard filter
-            //if (step % 20 == 0) {
-            //    filterDensity(me, system);
-            //}
+            if (step % 20 == 0) {
+                filterDensity(me, system);
+            }
             me->rho = me->rho + 0.5 * system->dt * me->Frho;
 
           // Solid (wall) particles should change density
         }else if (me->solidTag == 1 && system->static_domain == 0) {
             // Filter density field (for fixed solid particles)
-            //if (step % 20 == 0) {
-            //    filterDensity(me, system);
-            //}
+            if (step % 20 == 0) {
+                filterDensity(me, system);
+            }
         }