Adding a new SPH solver that is more physically accurate. See patch #29681
authorAlex Fraser <alex@phatcore.com>
Fri, 14 Dec 2012 04:57:26 +0000 (04:57 +0000)
committerAlex Fraser <alex@phatcore.com>
Fri, 14 Dec 2012 04:57:26 +0000 (04:57 +0000)
    http://projects.blender.org/tracker/index.php?func=detail&aid=29681&group_id=9&atid=127

The solver was mostly implemented by John Mansour at VPAC, with help from me and with funding from the AutoCRC. The SPH formulation is due to Gingold and Monaghan, and the smoothing kernel is due to Wendland.

This solver does not replace the old one; it is available as an option. Note that the new solver uses different units than the old one. The patch page has a couple of attachments that can be used to test the new solver, particularly sphclassical_dam_s0.01_grav.blend (ignore the earlier tests). The simulation in that file compares well with a physical experimental dam break; details in a paper by Changhong Hu and Makoto Sueyoshi, also referred to on that page.

release/scripts/startup/bl_ui/properties_particle.py
source/blender/blenkernel/BKE_particle.h
source/blender/blenkernel/intern/particle_system.c
source/blender/makesdna/DNA_particle_types.h
source/blender/makesrna/intern/rna_particle.c

index 75154550bcba84f5279e39e8cefebf194a3d1a8b..2c2ced5db0c51cb6221a5eedc829b83975949328 100644 (file)
@@ -505,6 +505,10 @@ class PARTICLE_PT_physics(ParticleButtonsPanel, Panel):
             if part.physics_type == 'FLUID':
                 fluid = part.fluid
 
+                split = layout.split()
+                sub = split.row()
+                sub.prop(fluid, "solver", expand=True)
+
                 split = layout.split()
 
                 col = split.column()
@@ -516,13 +520,14 @@ class PARTICLE_PT_physics(ParticleButtonsPanel, Panel):
                 col = split.column()
                 col.label(text="Advanced:")
 
-                sub = col.row()
-                sub.prop(fluid, "repulsion", slider=fluid.factor_repulsion)
-                sub.prop(fluid, "factor_repulsion", text="")
+                if fluid.solver == 'DDR':
+                    sub = col.row()
+                    sub.prop(fluid, "repulsion", slider=fluid.factor_repulsion)
+                    sub.prop(fluid, "factor_repulsion", text="")
 
-                sub = col.row()
-                sub.prop(fluid, "stiff_viscosity", slider=fluid.factor_stiff_viscosity)
-                sub.prop(fluid, "factor_stiff_viscosity", text="")
+                    sub = col.row()
+                    sub.prop(fluid, "stiff_viscosity", slider=fluid.factor_stiff_viscosity)
+                    sub.prop(fluid, "factor_stiff_viscosity", text="")
 
                 sub = col.row()
                 sub.prop(fluid, "fluid_radius", slider=fluid.factor_radius)
@@ -532,27 +537,37 @@ class PARTICLE_PT_physics(ParticleButtonsPanel, Panel):
                 sub.prop(fluid, "rest_density", slider=fluid.factor_density)
                 sub.prop(fluid, "factor_density", text="")
 
-                split = layout.split()
-
-                col = split.column()
-                col.label(text="Springs:")
-                col.prop(fluid, "spring_force", text="Force")
-                col.prop(fluid, "use_viscoelastic_springs")
-                sub = col.column(align=True)
-                sub.active = fluid.use_viscoelastic_springs
-                sub.prop(fluid, "yield_ratio", slider=True)
-                sub.prop(fluid, "plasticity", slider=True)
-
-                col = split.column()
-                col.label(text="Advanced:")
-                sub = col.row()
-                sub.prop(fluid, "rest_length", slider=fluid.factor_rest_length)
-                sub.prop(fluid, "factor_rest_length", text="")
-                col.label(text="")
-                sub = col.column()
-                sub.active = fluid.use_viscoelastic_springs
-                sub.prop(fluid, "use_initial_rest_length")
-                sub.prop(fluid, "spring_frames", text="Frames")
+                if fluid.solver == 'CLASSICAL':
+                    # With the classical solver, it is possible to calculate the
+                    # spacing between particles when the fluid is at rest. This
+                    # makes it easier to set stable initial conditions.
+                    particle_volume = part.mass / fluid.rest_density
+                    spacing = pow(particle_volume, 1/3)
+                    sub = col.row()
+                    sub.label(text="Spacing: %g" % spacing)
+
+                elif fluid.solver == 'DDR':
+                    split = layout.split()
+
+                    col = split.column()
+                    col.label(text="Springs:")
+                    col.prop(fluid, "spring_force", text="Force")
+                    col.prop(fluid, "use_viscoelastic_springs")
+                    sub = col.column(align=True)
+                    sub.active = fluid.use_viscoelastic_springs
+                    sub.prop(fluid, "yield_ratio", slider=True)
+                    sub.prop(fluid, "plasticity", slider=True)
+
+                    col = split.column()
+                    col.label(text="Advanced:")
+                    sub = col.row()
+                    sub.prop(fluid, "rest_length", slider=fluid.factor_rest_length)
+                    sub.prop(fluid, "factor_rest_length", text="")
+                    col.label(text="")
+                    sub = col.column()
+                    sub.active = fluid.use_viscoelastic_springs
+                    sub.prop(fluid, "use_initial_rest_length")
+                    sub.prop(fluid, "spring_frames", text="Frames")
 
         elif part.physics_type == 'KEYED':
             split = layout.split()
index bdaffef6818e75cd1332fc89f8a48db2912c5baa..f15ad296e4aa5d91c7b0a81295409ce37414046a 100644 (file)
@@ -20,7 +20,9 @@
  *
  * The Original Code is: all of this file.
  *
- * Contributor(s): none yet.
+ * Adaptive time step
+ * Classical SPH
+ * Copyright 2011-2012 AutoCRC
  *
  * ***** END GPL LICENSE BLOCK *****
  */
@@ -58,6 +60,7 @@ struct RNG;
 struct SurfaceModifierData;
 struct BVHTreeRay;
 struct BVHTreeRayHit; 
+struct EdgeHash;
 
 #define PARTICLE_P              ParticleData * pa; int p
 #define LOOP_PARTICLES  for (p = 0, pa = psys->particles; p < psys->totpart; p++, pa++)
@@ -85,6 +88,24 @@ typedef struct ParticleSimulationData {
        float courant_num;
 } ParticleSimulationData;
 
+typedef struct SPHData {
+       ParticleSystem *psys[10];
+       ParticleData *pa;
+       float mass;
+       struct EdgeHash *eh;
+       float *gravity;
+       float hfac;
+       /* Average distance to neighbours (other particles in the support domain),
+          for calculating the Courant number (adaptive time step). */
+       int pass;
+       float element_size;
+       float flow[3];
+
+       /* Integrator callbacks. This allows different SPH implementations. */
+       void (*force_cb) (void *sphdata_v, ParticleKey *state, float *force, float *impulse);
+       void (*density_cb) (void *rangedata_v, int index, float squared_dist);
+} SPHData;
+
 typedef struct ParticleTexture {
        float ivel;                           /* used in reset */
        float time, life, exist, size;        /* used in init */
@@ -283,6 +304,10 @@ float psys_get_child_size(struct ParticleSystem *psys, struct ChildParticle *cpa
 void psys_get_particle_on_path(struct ParticleSimulationData *sim, int pa_num, struct ParticleKey *state, int vel);
 int psys_get_particle_state(struct ParticleSimulationData *sim, int p, struct ParticleKey *state, int always);
 
+void psys_sph_init(struct ParticleSimulationData *sim, struct SPHData *sphdata);
+void psys_sph_finalise(struct SPHData *sphdata);
+void psys_sph_density(struct BVHTree *tree, struct SPHData* data, float co[3], float vars[2]);
+
 /* for anim.c */
 void psys_get_dupli_texture(struct ParticleSystem *psys, struct ParticleSettings *part,
                             struct ParticleSystemModifierData *psmd, struct ParticleData *pa, struct ChildParticle *cpa,
index d5fabf600790920a683bc2aa47911e509dd120b4..7c95265b27b80c9d898505a94f57fb82f4f24f51 100644 (file)
@@ -23,7 +23,8 @@
  * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn.
  *
  * Adaptive time step
- * Copyright 2011 AutoCRC
+ * Classical SPH
+ * Copyright 2011-2012 AutoCRC
  *
  * ***** END GPL LICENSE BLOCK *****
  */
@@ -1889,7 +1890,6 @@ void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime,
                bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f;
        }
 
-
        if (part->type == PART_HAIR) {
                pa->lifetime = 100.0f;
        }
@@ -2390,33 +2390,34 @@ typedef struct SPHRangeData {
        SPHNeighbor neighbors[SPH_NEIGHBORS];
        int tot_neighbors;
 
-       float density, near_density;
-       float h;
+       float* data;
 
        ParticleSystem *npsys;
        ParticleData *pa;
 
+       float h;
        float massfac;
        int use_size;
 } SPHRangeData;
 
-typedef struct SPHData {
-       ParticleSystem *psys[10];
-       ParticleData *pa;
-       float mass;
-       EdgeHash *eh;
-       float *gravity;
-       /* Average distance to neighbors (other particles in the support domain),
-        * for calculating the Courant number (adaptive time step). */
-       int pass;
-       float element_size;
-       float flow[3];
-
-       /* Integrator callbacks. This allows different SPH implementations. */
-       void (*force_cb) (void *sphdata_v, ParticleKey *state, float *force, float *impulse);
-       void (*density_cb) (void *rangedata_v, int index, float squared_dist);
-} SPHData;
+static void sph_evaluate_func(BVHTree *tree, ParticleSystem **psys, float co[3], SPHRangeData *pfr, float interaction_radius, BVHTree_RangeQuery callback) {
+       int i;
+
+       pfr->tot_neighbors = 0;
+
+       for (i=0; i < 10 && psys[i]; i++) {
+               pfr->npsys    = psys[i];
+               pfr->massfac  = psys[i]->part->mass;
+               pfr->use_size = psys[i]->part->flag & PART_SIZEMASS;
 
+               if (tree) {
+                       BLI_bvhtree_range_query(tree, co, interaction_radius, callback, pfr);
+                       break;
+               } else {
+                       BLI_bvhtree_range_query(psys[i]->bvhtree, co, interaction_radius, callback, pfr);
+               }
+       }
+}
 static void sph_density_accum_cb(void *userdata, int index, float squared_dist)
 {
        SPHRangeData *pfr = (SPHRangeData *)userdata;
@@ -2446,8 +2447,8 @@ static void sph_density_accum_cb(void *userdata, int index, float squared_dist)
        if (pfr->use_size)
                q *= npa->size;
 
-       pfr->density += q*q;
-       pfr->near_density += q*q*q;
+       pfr->data[0] += q*q;
+       pfr->data[1] += q*q*q;
 }
 
 /*
@@ -2500,8 +2501,10 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa
 
        float inv_mass = 1.0f/mass;
        float spring_constant = fluid->spring_k;
-       
-       float h = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.f*pa->size : 1.f); /* 4.0 seems to be a pretty good value */
+
+       /* 4.0 seems to be a pretty good value */
+       float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
+       float h = interaction_radius * sphdata->hfac;
        float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f); /* 4.77 is an experimentally determined density factor */
        float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
 
@@ -2512,24 +2515,23 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa
        float vec[3];
        float vel[3];
        float co[3];
+       float data[2];
+       float density, near_density;
 
        int i, spring_index, index = pa - psys[0]->particles;
 
-       pfr.tot_neighbors = 0;
-       pfr.density = pfr.near_density = 0.f;
+       data[0] = data[1] = 0;
+       pfr.data = data;
        pfr.h = h;
        pfr.pa = pa;
 
-       for (i=0; i<10 && psys[i]; i++) {
-               pfr.npsys = psys[i];
-               pfr.massfac = psys[i]->part->mass*inv_mass;
-               pfr.use_size = psys[i]->part->flag & PART_SIZEMASS;
+       sph_evaluate_func( NULL, psys, state->co, &pfr, interaction_radius, sph_density_accum_cb);
 
-               BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sphdata->density_cb, &pfr);
-       }
+       density = data[0];
+       near_density = data[1];
 
-       pressure =  stiffness * (pfr.density - rest_density);
-       near_pressure = stiffness_near_fac * pfr.near_density;
+       pressure =  stiffness * (density - rest_density);
+       near_pressure = stiffness_near_fac * near_density;
 
        pfn = pfr.neighbors;
        for (i=0; i<pfr.tot_neighbors; i++, pfn++) {
@@ -2593,14 +2595,202 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa
        
        /* Artificial buoyancy force in negative gravity direction  */
        if (fluid->buoyancy > 0.f && gravity)
-               madd_v3_v3fl(force, gravity, fluid->buoyancy * (pfr.density-rest_density));
+               madd_v3_v3fl(force, gravity, fluid->buoyancy * (density-rest_density));
+
+       if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF)
+               sph_particle_courant(sphdata, &pfr);
+       sphdata->pass++;
+}
+
+/* powf is really slow for raising to integer powers. */
+MINLINE float pow2(float x) {
+       return x * x;
+}
+MINLINE float pow3(float x) {
+       return pow2(x) * x;
+}
+MINLINE float pow4(float x) {
+       return pow2(pow2(x));
+}
+MINLINE float pow7(float x) {
+       return pow2(pow3(x)) * x;
+}
+
+static void sphclassical_density_accum_cb(void *userdata, int index, float squared_dist)
+{
+       SPHRangeData *pfr = (SPHRangeData *)userdata;
+       ParticleData *npa = pfr->npsys->particles + index;
+       float q;
+       float qfac = 21.0f / (256.f * M_PI);
+       float rij, rij_h;
+       float vec[3];
+
+       /* Exclude particles that are more than 2h away. Can't use squared_dist here
+        * because it is not accurate enough. Use current state, i.e. the output of
+        * basic_integrate() - z0r */
+       sub_v3_v3v3(vec, npa->state.co, pfr->pa->state.co);
+       rij = len_v3(vec);
+       rij_h = rij / pfr->h;
+       if (rij_h > 2.0f)
+               return;
+
+       /* Smoothing factor. Utilise the Wendland kernel. gnuplot:
+        *     q1(x) = (2.0 - x)**4 * ( 1.0 + 2.0 * x)
+        *     plot [0:2] q1(x) */
+       q  = qfac / pow3(pfr->h) * pow4(2.0f - rij_h) * ( 1.0f + 2.0f * rij_h);
+       q *= pfr->massfac;
+
+       if (pfr->use_size)
+               q *= pfr->pa->size;
+
+       pfr->data[0] += q;
+       pfr->data[1] += q / npa->sphdensity;
+}
+
+static void sphclassical_neighbour_accum_cb(void *userdata, int index, float squared_dist)
+{
+       SPHRangeData *pfr = (SPHRangeData *)userdata;
+       ParticleData *npa = pfr->npsys->particles + index;
+       float rij, rij_h;
+       float vec[3];
+
+       if (pfr->tot_neighbors >= SPH_NEIGHBORS)
+               return;
+
+       /* Exclude particles that are more than 2h away. Can't use squared_dist here
+        * because it is not accurate enough. Use current state, i.e. the output of
+        * basic_integrate() - z0r */
+       sub_v3_v3v3(vec, npa->state.co, pfr->pa->state.co);
+       rij = len_v3(vec);
+       rij_h = rij / pfr->h;
+       if (rij_h > 2.0f)
+               return;
+
+       pfr->neighbors[pfr->tot_neighbors].index = index;
+       pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
+       pfr->tot_neighbors++;
+}
+static void sphclassical_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
+{
+       SPHData *sphdata = (SPHData *)sphdata_v;
+       ParticleSystem **psys = sphdata->psys;
+       ParticleData *pa = sphdata->pa;
+       SPHFluidSettings *fluid = psys[0]->part->fluid;
+       SPHRangeData pfr;
+       SPHNeighbor *pfn;
+       float *gravity = sphdata->gravity;
+
+       float dq, u, rij, dv[3];
+       float pressure, npressure;
+
+       float visc = fluid->viscosity_omega;
+
+       float interaction_radius;
+       float h, hinv;
+       /* 4.77 is an experimentally determined density factor */
+       float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.0f);
+
+       float stiffness = fluid->stiffness_k;
+
+       ParticleData *npa;
+       float vec[3];
+       float co[3];
+       float pressureTerm;
+
+       int i;
+
+       float qfac2 = 42.0f / (256.0f * M_PI);
+       float rij_h;
+
+       /* 4.0 here is to be consistent with previous formulation/interface */
+       interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
+       h = interaction_radius * sphdata->hfac;
+       hinv = 1.0f / h;
+
+       pfr.h = h;
+       pfr.pa = pa;
+
+       sph_evaluate_func(NULL, psys, state->co, &pfr, interaction_radius, sphclassical_neighbour_accum_cb);
+       pressure =  stiffness * (pow7(pa->sphdensity / rest_density) - 1.0f);
+
+       /* multiply by mass so that we return a force, not accel */
+       qfac2 *= sphdata->mass / pow3(pfr.h);
+
+       pfn = pfr.neighbors;
+       for (i = 0; i < pfr.tot_neighbors; i++, pfn++) {
+               npa = pfn->psys->particles + pfn->index;
+               if (npa == pa) {
+                       /* we do not contribute to ourselves */
+                       continue;
+               }
+
+               /* Find vector to neighbour. Exclude particles that are more than 2h
+                * away. Can't use current state here because it may have changed on
+                * another thread - so do own mini integration. Unlike basic_integrate,
+                * SPH integration depends on neighbouring particles. - z0r */
+               madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
+               sub_v3_v3v3(vec, co, state->co);
+               rij = normalize_v3(vec);
+               rij_h = rij / pfr.h;
+               if (rij_h > 2.0f)
+                       continue;
+
+               npressure = stiffness * (pow7(npa->sphdensity / rest_density) - 1.0f);
+
+               /* First derivative of smoothing factor. Utilise the Wendland kernel.
+                * gnuplot:
+                *     q2(x) = 2.0 * (2.0 - x)**4 - 4.0 * (2.0 - x)**3 * (1.0 + 2.0 * x)
+                *     plot [0:2] q2(x)
+                * Particles > 2h away are excluded above. */
+               dq = qfac2 * (2.0f * pow4(2.0f - rij_h) - 4.0f * pow3(2.0f - rij_h) * (1.0f + 2.0f * rij_h)  );
+
+               if (pfn->psys->part->flag & PART_SIZEMASS)
+                       dq *= npa->size;
+
+               pressureTerm = pressure / pow2(pa->sphdensity) + npressure / pow2(npa->sphdensity);
+
+               /* Note that 'minus' is removed, because vec = vecBA, not vecAB.
+                * This applies to the viscosity calculation below, too. */
+               madd_v3_v3fl(force, vec, pressureTerm * dq);
+
+               /* Viscosity */
+               if (visc > 0.0f) {
+                       sub_v3_v3v3(dv, npa->prev_state.vel, pa->prev_state.vel);
+                       u = dot_v3v3(vec, dv);
+                       /* Apply parameters */
+                       u *= -dq * hinv * visc / (0.5f * npa->sphdensity + 0.5f * pa->sphdensity);
+                       madd_v3_v3fl(force, vec, u);
+               }
+       }
+
+       /* Artificial buoyancy force in negative gravity direction  */
+       if (fluid->buoyancy > 0.f && gravity)
+               madd_v3_v3fl(force, gravity, fluid->buoyancy * (pa->sphdensity - rest_density));
 
        if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF)
                sph_particle_courant(sphdata, &pfr);
        sphdata->pass++;
 }
 
-static void sph_solver_init(ParticleSimulationData *sim, SPHData *sphdata)
+static void sphclassical_calc_dens(ParticleData *pa, float UNUSED(dfra), SPHData *sphdata){
+       ParticleSystem **psys = sphdata->psys;
+       SPHFluidSettings *fluid = psys[0]->part->fluid;
+       /* 4.0 seems to be a pretty good value */
+       float interaction_radius  = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
+       SPHRangeData pfr;
+       float data[2];
+
+       data[0] = 0;
+       data[1] = 0;
+       pfr.data = data;
+       pfr.h = interaction_radius * sphdata->hfac;
+       pfr.pa = pa;
+
+       sph_evaluate_func( NULL, psys, pa->state.co, &pfr, interaction_radius, sphclassical_density_accum_cb);
+       pa->sphdensity = MIN2(MAX2(data[0], fluid->rest_density * 0.9f), fluid->rest_density * 1.1f);
+}
+
+void psys_sph_init(ParticleSimulationData *sim, SPHData *sphdata)
 {
        ParticleTarget *pt;
        int i;
@@ -2621,17 +2811,44 @@ static void sph_solver_init(ParticleSimulationData *sim, SPHData *sphdata)
        sphdata->pa = NULL;
        sphdata->mass = 1.0f;
 
-       sphdata->force_cb = sph_force_cb;
-       sphdata->density_cb = sph_density_accum_cb;
+       if (sim->psys->part->fluid->solver == SPH_SOLVER_DDR) {
+               sphdata->force_cb = sph_force_cb;
+               sphdata->density_cb = sph_density_accum_cb;
+               sphdata->hfac = 1.0f;
+       } else {
+               /* SPH_SOLVER_CLASSICAL */
+               sphdata->force_cb = sphclassical_force_cb;
+               sphdata->density_cb = sphclassical_density_accum_cb;
+               sphdata->hfac = 0.5f;
+       }
+
 }
 
-static void sph_solver_finalise(SPHData *sphdata)
+void psys_sph_finalise(SPHData *sphdata)
 {
        if (sphdata->eh) {
                BLI_edgehash_free(sphdata->eh, NULL);
                sphdata->eh = NULL;
        }
 }
+/* Sample the density field at a point in space. */
+void psys_sph_density(BVHTree *tree, SPHData *sphdata, float co[3], float vars[2]) {
+       ParticleSystem **psys = sphdata->psys;
+       SPHFluidSettings *fluid = psys[0]->part->fluid;
+       /* 4.0 seems to be a pretty good value */
+       float interaction_radius  = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
+       SPHRangeData pfr;
+       float density[2];
+
+       density[0] = density[1] = 0.0f;
+       pfr.data = density;
+       pfr.h = interaction_radius*sphdata->hfac;
+
+       sph_evaluate_func(tree, psys, co, &pfr, interaction_radius, sphdata->density_cb);
+
+       vars[0] = pfr.data[0];
+       vars[1] = pfr.data[1];
+}
 
 static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, SPHData *sphdata)
 {
@@ -3815,15 +4032,14 @@ static float update_timestep(ParticleSystem *psys, ParticleSimulationData *sim,
        else
                dt_target = psys->dt_frac * (psys->part->courant_target / sim->courant_num);
 
-       /* Make sure the time step is reasonable.
-        * For some reason, the CLAMP macro doesn't work here. The time step becomes
-        * too large. - z0r */
+       /* Make sure the time step is reasonable. For some reason, the CLAMP macro
+        * doesn't work here. The time step becomes too large. - z0r */
        if (dt_target < MIN_TIMESTEP)
                dt_target = MIN_TIMESTEP;
        else if (dt_target > get_base_time_step(psys->part))
                dt_target = get_base_time_step(psys->part);
 
-       /* Decrease time step instantly, but expand slowly. */
+       /* Decrease time step instantly, but increase slowly. */
        if (dt_target > psys->dt_frac)
                psys->dt_frac = interpf(dt_target, psys->dt_frac, TIMESTEP_EXPANSION_FACTOR);
        else
@@ -3991,31 +4207,71 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra)
                case PART_PHYS_FLUID:
                {
                        SPHData sphdata;
-                       sph_solver_init(sim, &sphdata);
+                       ParticleSettings *part = sim->psys->part;
+                       psys_sph_init(sim, &sphdata);
 
-                       #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5)
-                       LOOP_DYNAMIC_PARTICLES {
-                               /* do global forces & effectors */
-                               basic_integrate(sim, p, pa->state.time, cfra);
+                       if (part->fluid->flag & SPH_SOLVER_DDR) {
+                               /* Apply SPH forces using double-density relaxation algorithm
+                                * (Clavat et. al.) */
+                               #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5)
+                               LOOP_DYNAMIC_PARTICLES {
+                                       /* do global forces & effectors */
+                                       basic_integrate(sim, p, pa->state.time, cfra);
 
-                               /* actual fluids calculations */
-                               sph_integrate(sim, pa, pa->state.time, &sphdata);
+                                       /* actual fluids calculations */
+                                       sph_integrate(sim, pa, pa->state.time, &sphdata);
 
-                               if (sim->colliders)
-                                       collision_check(sim, p, pa->state.time, cfra);
+                                       if(sim->colliders)
+                                               collision_check(sim, p, pa->state.time, cfra);
+
+                                       /* SPH particles are not physical particles, just interpolation
+                                        * particles,  thus rotation has not a direct sense for them */
+                                       basic_rotate(part, pa, pa->state.time, timestep);
+
+                                       #pragma omp critical
+                                       if (part->time_flag & PART_TIME_AUTOSF)
+                                               update_courant_num(sim, pa, dtime, &sphdata);
+                               }
+
+                               sph_springs_modify(psys, timestep);
+
+                       } else {
+                               /* SPH_SOLVER_CLASSICAL */
+                               /* Apply SPH forces using classical algorithm (due to Gingold
+                                * and Monaghan). Note that, unlike double-density relaxation,
+                                * this algorthim is separated into distinct loops. */
+
+                               #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5)
+                               LOOP_DYNAMIC_PARTICLES {
+                                       basic_integrate(sim, p, pa->state.time, cfra);
+                               }
+
+                               /* calculate summation density */
+                               #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5)
+                               LOOP_DYNAMIC_PARTICLES {
+                                       sphclassical_calc_dens(pa, pa->state.time, &sphdata);
+                               }
+
+                               /* do global forces & effectors */
+                               #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5)
+                               LOOP_DYNAMIC_PARTICLES {
+                                       /* actual fluids calculations */
+                                       sph_integrate(sim, pa, pa->state.time, &sphdata);
+
+                                       if(sim->colliders)
+                                               collision_check(sim, p, pa->state.time, cfra);
                                
-                               /* SPH particles are not physical particles, just interpolation
-                                * particles,  thus rotation has not a direct sense for them */
-                               basic_rotate(part, pa, pa->state.time, timestep);  
+                                       /* SPH particles are not physical particles, just interpolation
+                                        * particles,  thus rotation has not a direct sense for them */
+                                       basic_rotate(part, pa, pa->state.time, timestep);
 
-                               #pragma omp critical
-                               if (part->time_flag & PART_TIME_AUTOSF)
-                                       update_courant_num(sim, pa, dtime, &sphdata);
+                                       #pragma omp critical
+                                       if (part->time_flag & PART_TIME_AUTOSF)
+                                               update_courant_num(sim, pa, dtime, &sphdata);
+                               }
                        }
 
-                       sph_springs_modify(psys, timestep);
-
-                       sph_solver_finalise(&sphdata);
+                       psys_sph_finalise(&sphdata);
                        break;
                }
        }
@@ -4325,18 +4581,16 @@ static void system_step(ParticleSimulationData *sim, float cfra)
                if ((int)cfra == startframe && part->sta < startframe)
                        totframesback = (startframe - (int)part->sta);
 
-               /* Initialise time step. This may change if automatic subframes are
-                * enabled. */
                if (!(part->time_flag & PART_TIME_AUTOSF)) {
                        /* Constant time step */
                        psys->dt_frac = get_base_time_step(part);
                }
                else if ((int)cfra == startframe) {
-                       /* Variable time step; use a very conservative value at the start.
-                        * If it doesn't need to be so small, it will quickly grow. */
+                       /* Variable time step; initialise to subframes */
                        psys->dt_frac = get_base_time_step(part);
                }
                else if (psys->dt_frac < MIN_TIMESTEP) {
+                       /* Variable time step; subsequent frames */
                        psys->dt_frac = MIN_TIMESTEP;
                }
 
index 5952aa8afb03478c94b2ced00ab149e27e94dd4c..1cf16fb52b3a2e39188d8ff360011e0217c4bd95 100644 (file)
@@ -116,6 +116,9 @@ typedef struct ParticleData {
 
        float size;                             /* size and multiplier so that we can update size when ever */
 
+       float sphdensity;               /* density of sph particle */
+       int pad;
+
        int hair_index;
        short flag;
        short alive;                    /* the life state of a particle */
@@ -130,6 +133,8 @@ typedef struct SPHFluidSettings {
        float stiffness_k, stiffness_knear, rest_density;
        float buoyancy;
        int flag, spring_frames;
+       short solver;
+       short pad[3];
 } SPHFluidSettings;
 
 /* fluid->flag */
@@ -141,6 +146,10 @@ typedef struct SPHFluidSettings {
 #define SPH_FAC_VISCOSITY                      32
 #define SPH_FAC_REST_LENGTH                    64
 
+/* fluid->solver (numerical ID field, not bitfield) */
+#define SPH_SOLVER_DDR                                 0
+#define SPH_SOLVER_CLASSICAL                   1
+
 typedef struct ParticleSettings {
        ID id;
        struct AnimData *adt;
index 96b6197cf69c0148500af675c912fe5b0b0bab45..93f940b1aa30885360e688412737316d7d418e1b 100644 (file)
@@ -1112,12 +1112,25 @@ static void rna_def_fluid_settings(BlenderRNA *brna)
 {
        StructRNA *srna;
        PropertyRNA *prop;
-       
+
+       static EnumPropertyItem sph_solver_items[] = {
+               {SPH_SOLVER_DDR, "DDR", 0, "Double-Density", "An artistic solver with strong surface tension effects (original)"},
+               {SPH_SOLVER_CLASSICAL, "CLASSICAL", 0, "Classical", "A more physically-accurate solver"},
+               {0, NULL, 0, NULL, NULL}
+       };
+
        srna = RNA_def_struct(brna, "SPHFluidSettings", NULL);
        RNA_def_struct_path_func(srna, "rna_SPHFluidSettings_path");
        RNA_def_struct_ui_text(srna, "SPH Fluid Settings", "Settings for particle fluids physics");
-       
+
        /* Fluid settings */
+       prop = RNA_def_property(srna, "solver", PROP_ENUM, PROP_NONE);
+       RNA_def_property_enum_sdna(prop, NULL, "solver");
+       RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+       RNA_def_property_enum_items(prop, sph_solver_items);
+       RNA_def_property_ui_text(prop, "SPH Solver", "The code used to calculate internal forces on particles");
+       RNA_def_property_update(prop, 0, "rna_Particle_reset");
+
        prop = RNA_def_property(srna, "spring_force", PROP_FLOAT, PROP_NONE);
        RNA_def_property_float_sdna(prop, NULL, "spring_k");
        RNA_def_property_range(prop, 0.0f, 100.0f);
@@ -1186,7 +1199,7 @@ static void rna_def_fluid_settings(BlenderRNA *brna)
        /* Double density relaxation */
        prop = RNA_def_property(srna, "stiffness", PROP_FLOAT, PROP_NONE);
        RNA_def_property_float_sdna(prop, NULL, "stiffness_k");
-       RNA_def_property_range(prop, 0.0f, 100.0f);
+       RNA_def_property_range(prop, 0.0f, 100000.0f);
        RNA_def_property_ui_range(prop, 0.0f, 10.0f, 1, 3);
        RNA_def_property_ui_text(prop, "Stiffness", "How incompressible the fluid is");
        RNA_def_property_update(prop, 0, "rna_Particle_reset");
@@ -1201,7 +1214,7 @@ static void rna_def_fluid_settings(BlenderRNA *brna)
 
        prop = RNA_def_property(srna, "rest_density", PROP_FLOAT, PROP_NONE);
        RNA_def_property_float_sdna(prop, NULL, "rest_density");
-       RNA_def_property_range(prop, 0.0f, 100.0f);
+       RNA_def_property_range(prop, 0.0f, 10000.0f);
        RNA_def_property_ui_range(prop, 0.0f, 2.0f, 1, 3);
        RNA_def_property_ui_text(prop, "Rest Density", "Fluid rest density");
        RNA_def_property_update(prop, 0, "rna_Particle_reset");
@@ -2281,7 +2294,7 @@ static void rna_def_particle_settings(BlenderRNA *brna)
 
        /* physical properties */
        prop = RNA_def_property(srna, "mass", PROP_FLOAT, PROP_NONE);
-       RNA_def_property_range(prop, 0.001f, 100000.0f);
+       RNA_def_property_range(prop, 0.00000001f, 100000.0f);
        RNA_def_property_ui_range(prop, 0.01, 100, 1, 3);
        RNA_def_property_ui_text(prop, "Mass", "Mass of the particles");
        RNA_def_property_update(prop, 0, "rna_Particle_reset");