Adding a new SPH solver that is more physically accurate. See patch #29681

    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.

2 files changed, 345 insertions, 66 deletions

diff --git a/source/blender/blenkernel/BKE_particle.h b/source/blender/blenkernel/BKE_particle.h
index bdaffef6818..f15ad296e4a 100644
--- a/source/blender/blenkernel/BKE_particle.h
+++ b/source/blender/blenkernel/BKE_particle.h
@@ -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,
diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c
index d5fabf60079..7c95265b27b 100644
--- a/source/blender/blenkernel/intern/particle_system.c
+++ b/source/blender/blenkernel/intern/particle_system.c
@@ -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;
 		}