Fix T73988: Mantaflow fluid simulation - Particles for Spray, Foam and Bubbles are one frame ahead of Mesh

Fixes an issue with secondary particles being out of sync with the main simulation. Cleaned up the secondary particle code in general too (making sure that all solver attributes - timestep, framelength, etc. - are set correctly).

3 files changed, 62 insertions, 42 deletions

diff --git a/intern/mantaflow/intern/MANTA_main.cpp b/intern/mantaflow/intern/MANTA_main.cpp
index 5d35de7898f..c169d242583 100644
--- a/intern/mantaflow/intern/MANTA_main.cpp
+++ b/intern/mantaflow/intern/MANTA_main.cpp
@@ -837,6 +837,8 @@ std::string MANTA::getRealValue(const std::string &varName, FluidModifierData *m
     ss << mmd->domain->flame_smoke_color[2];
   else if (varName == "CURRENT_FRAME")
     ss << mmd->time;
+  else if (varName == "START_FRAME")
+    ss << mmd->domain->cache_frame_start;
   else if (varName == "END_FRAME")
     ss << mmd->domain->cache_frame_end;
   else if (varName == "CACHE_DATA_FORMAT")
diff --git a/intern/mantaflow/intern/strings/fluid_script.h b/intern/mantaflow/intern/strings/fluid_script.h
index c442dd56c09..0aad0546aea 100644
--- a/intern/mantaflow/intern/strings/fluid_script.h
+++ b/intern/mantaflow/intern/strings/fluid_script.h
@@ -123,6 +123,11 @@ cflCond_s$ID$      = $CFL$\n\
 timestepsMin_s$ID$ = $TIMESTEPS_MIN$\n\
 timestepsMax_s$ID$ = $TIMESTEPS_MAX$\n\
 \n\
+# Start and stop for simulation\n\
+current_frame_s$ID$ = $CURRENT_FRAME$\n\
+start_frame_s$ID$   = $START_FRAME$\n\
+end_frame_s$ID$     = $END_FRAME$\n\
+\n\
 # Fluid diffusion / viscosity\n\
 domainSize_s$ID$ = $FLUID_DOMAIN_SIZE$ # longest domain side in meters\n\
 viscosity_s$ID$ = $FLUID_VISCOSITY$ / (domainSize_s$ID$*domainSize_s$ID$) # kinematic viscosity in m^2/s\n\
@@ -208,6 +213,8 @@ def fluid_adapt_time_step_$ID$():\n\
     # time params are animatable\n\
     s$ID$.frameLength = frameLength_s$ID$\n\
     s$ID$.cfl         = cflCond_s$ID$\n\
+    s$ID$.timestepMin  = s$ID$.frameLength / max(1, timestepsMax_s$ID$)\n\
+    s$ID$.timestepMax  = s$ID$.frameLength / max(1, timestepsMin_s$ID$)\n\
     \n\
     # ensure that vel grid is full (remember: adaptive domain can reallocate solver)\n\
     copyRealToVec3(sourceX=x_vel_s$ID$, sourceY=y_vel_s$ID$, sourceZ=z_vel_s$ID$, target=vel_s$ID$)\n\
@@ -225,7 +232,7 @@ const std::string fluid_alloc =
 mantaMsg('Fluid alloc data')\n\
 flags_s$ID$       = s$ID$.create(FlagGrid)\n\
 vel_s$ID$         = s$ID$.create(MACGrid)\n\
-velC_s$ID$        = s$ID$.create(MACGrid)\n\
+velTmp_s$ID$      = s$ID$.create(MACGrid)\n\
 x_vel_s$ID$       = s$ID$.create(RealGrid)\n\
 y_vel_s$ID$       = s$ID$.create(RealGrid)\n\
 z_vel_s$ID$       = s$ID$.create(RealGrid)\n\
@@ -247,7 +254,7 @@ phiIn_s$ID$.setConst(9999)\n\
 phiOut_s$ID$.setConst(9999)\n\
 \n\
 # Keep track of important objects in dict to load them later on\n\
-fluid_data_dict_final_s$ID$  = dict(vel=vel_s$ID$)\n\
+fluid_data_dict_final_s$ID$  = dict(vel=vel_s$ID$, velTmp=velTmp_s$ID$)\n\
 fluid_data_dict_resume_s$ID$ = dict(phiObs=phiObs_s$ID$, phiIn=phiIn_s$ID$, phiOut=phiOut_s$ID$, flags=flags_s$ID$)\n";
 
 const std::string fluid_alloc_obstacle =
@@ -493,7 +500,8 @@ def bake_fluid_process_data_$ID$(framenr, format_data, format_particles, format_
     mantaMsg('Bake fluid data')\n\
     \n\
     s$ID$.frame = framenr\n\
-    # Must not set 'timeTotal' here. Remember, this function is called from manta.c while-loop\n\
+    s$ID$.frameLength = frameLength_s$ID$\n\
+    s$ID$.timeTotal = timeTotal_s$ID$\n\
     \n\
     start_time = time.time()\n\
     if using_smoke_s$ID$:\n\
@@ -514,9 +522,9 @@ def bake_noise_process_$ID$(framenr, format_data, format_noise, path_data, path_
     mantaMsg('Bake fluid noise')\n\
     \n\
     sn$ID$.frame = framenr\n\
-    sn$ID$.timeTotal = (framenr-1) * frameLength_s$ID$\n\
-    sn$ID$.timestep  = dt0_s$ID$\n\
-    mantaMsg('sn$ID$.timeTotal: ' + str(sn$ID$.timeTotal))\n\
+    sn$ID$.frameLength = frameLength_s$ID$\n\
+    sn$ID$.timeTotal = abs(framenr - start_frame_s$ID$) * frameLength_s$ID$\n\
+    sn$ID$.timestep = frameLength_s$ID$ # no adaptive timestep for noise\n\
     \n\
     smoke_step_noise_$ID$(framenr)\n\
     smoke_save_noise_$ID$(path_noise, framenr, format_noise, resumable)\n\
@@ -533,8 +541,9 @@ def bake_mesh_process_$ID$(framenr, format_data, format_mesh, format_particles,
     mantaMsg('Bake fluid mesh')\n\
     \n\
     sm$ID$.frame = framenr\n\
-    sm$ID$.timeTotal = (framenr-1) * frameLength_s$ID$\n\
-    sm$ID$.timestep  = dt0_s$ID$\n\
+    sm$ID$.frameLength = frameLength_s$ID$\n\
+    sm$ID$.timeTotal = abs(framenr - start_frame_s$ID$) * frameLength_s$ID$\n\
+    sm$ID$.timestep = frameLength_s$ID$ # no adaptive timestep for mesh\n\
     \n\
     #if using_smoke_s$ID$:\n\
         # TODO (sebbas): Future update could include smoke mesh (vortex sheets)\n\
@@ -556,8 +565,9 @@ def bake_particles_process_$ID$(framenr, format_data, format_particles, path_dat
     mantaMsg('Bake secondary particles')\n\
     \n\
     sp$ID$.frame = framenr\n\
-    sp$ID$.timeTotal = (framenr-1) * frameLength_s$ID$\n\
-    sp$ID$.timestep  = dt0_s$ID$\n\
+    sp$ID$.frameLength = frameLength_s$ID$\n\
+    sp$ID$.timeTotal = abs(framenr - start_frame_s$ID$) * frameLength_s$ID$\n\
+    sp$ID$.timestep = frameLength_s$ID$ # no adaptive timestep for particles\n\
     \n\
     #if using_smoke_s$ID$:\n\
         # TODO (sebbas): Future update could include smoke particles (e.g. fire sparks)\n\
@@ -710,28 +720,24 @@ file_format_noise     = '$CACHE_NOISE_FORMAT$'\n\
 file_format_particles = '$CACHE_PARTICLE_FORMAT$'\n\
 file_format_mesh      = '$CACHE_MESH_FORMAT$'\n\
 \n\
-# Start and stop for simulation\n\
-current_frame  = $CURRENT_FRAME$\n\
-end_frame      = $END_FRAME$\n\
-\n\
 # How many frame to load from cache\n\
 from_cache_cnt = 100\n\
 \n\
 loop_cnt = 0\n\
-while current_frame <= end_frame:\n\
+while current_frame_s$ID$ <= end_frame_s$ID$:\n\
     \n\
     # Load already simulated data from cache:\n\
     if loop_cnt < from_cache_cnt:\n\
-        load(current_frame, cache_resumable)\n\
+        load(current_frame_s$ID$, cache_resumable)\n\
     \n\
     # Otherwise simulate new data\n\
     else:\n\
-        while(s$ID$.frame <= current_frame):\n\
+        while(s$ID$.frame <= current_frame_s$ID$):\n\
             if using_adaptTime_s$ID$:\n\
                 fluid_adapt_time_step_$ID$()\n\
-            step(current_frame)\n\
+            step(current_frame_s$ID$)\n\
     \n\
-    current_frame += 1\n\
+    current_frame_s$ID$ += 1\n\
     loop_cnt += 1\n\
     \n\
     if gui:\n\
diff --git a/intern/mantaflow/intern/strings/liquid_script.h b/intern/mantaflow/intern/strings/liquid_script.h
index 23ec16b9d84..2376d49aace 100644
--- a/intern/mantaflow/intern/strings/liquid_script.h
+++ b/intern/mantaflow/intern/strings/liquid_script.h
@@ -133,9 +133,8 @@ pLifeSnd_pp$ID$      = ppSnd_sp$ID$.create(PdataReal)\n\
 vel_sp$ID$           = sp$ID$.create(MACGrid)\n\
 flags_sp$ID$         = sp$ID$.create(FlagGrid)\n\
 phi_sp$ID$           = sp$ID$.create(LevelsetGrid)\n\
-phiIn_sp$ID$         = sp$ID$.create(LevelsetGrid)\n\
 phiObs_sp$ID$        = sp$ID$.create(LevelsetGrid)\n\
-phiObsIn_sp$ID$      = sp$ID$.create(LevelsetGrid)\n\
+phiOut_sp$ID$        = sp$ID$.create(LevelsetGrid)\n\
 normal_sp$ID$        = sp$ID$.create(VecGrid)\n\
 neighborRatio_sp$ID$ = sp$ID$.create(RealGrid)\n\
 trappedAir_sp$ID$    = sp$ID$.create(RealGrid)\n\
@@ -144,9 +143,8 @@ kineticEnergy_sp$ID$ = sp$ID$.create(RealGrid)\n\
 \n\
 # Set some initial values\n\
 phi_sp$ID$.setConst(9999)\n\
-phiIn_sp$ID$.setConst(9999)\n\
 phiObs_sp$ID$.setConst(9999)\n\
-phiObsIn_sp$ID$.setConst(9999)\n\
+phiOut_sp$ID$.setConst(9999)\n\
 \n\
 # Keep track of important objects in dict to load them later on\n\
 liquid_particles_dict_final_s$ID$  = dict(ppSnd=ppSnd_sp$ID$, pVelSnd=pVelSnd_pp$ID$, pLifeSnd=pLifeSnd_pp$ID$)\n\
@@ -229,13 +227,15 @@ def liquid_step_$ID$():\n\
     mantaMsg('Pushing particles out of obstacles')\n\
     pushOutofObs(parts=pp_s$ID$, flags=flags_s$ID$, phiObs=phiObs_s$ID$)\n\
     \n\
+    # save original states for later (used during mesh / secondary particle creation)\n\
+    phiTmp_s$ID$.copyFrom(phi_s$ID$)\n\
+    velTmp_s$ID$.copyFrom(vel_s$ID$)\n\
+    \n\
     mantaMsg('Advecting phi')\n\
     advectSemiLagrange(flags=flags_s$ID$, vel=vel_s$ID$, grid=phi_s$ID$, order=1) # first order is usually enough\n\
     mantaMsg('Advecting velocity')\n\
     advectSemiLagrange(flags=flags_s$ID$, vel=vel_s$ID$, grid=vel_s$ID$, order=2)\n\
     \n\
-    phiTmp_s$ID$.copyFrom(phi_s$ID$) # save original phi for later use in mesh creation\n\
-    \n\
     # create level set of particles\n\
     gridParticleIndex(parts=pp_s$ID$, flags=flags_s$ID$, indexSys=pindex_s$ID$, index=gpi_s$ID$)\n\
     unionParticleLevelset(parts=pp_s$ID$, indexSys=pindex_s$ID$, flags=flags_s$ID$, index=gpi_s$ID$, phi=phiParts_s$ID$, radiusFactor=radiusFactor_s$ID$)\n\
@@ -308,7 +308,13 @@ const std::string liquid_step_mesh =
 def liquid_step_mesh_$ID$():\n\
     mantaMsg('Liquid step mesh')\n\
     \n\
-    interpolateGrid(target=phi_sm$ID$, source=phiTmp_s$ID$)\n\
+    # no upres: just use the loaded grids\n\
+    if upres_sm$ID$ <= 1:\n\
+        phi_sm$ID$.copyFrom(phiTmp_s$ID$)\n\
+    \n\
+    # with upres: recreate grids\n\
+    else:\n\
+        interpolateGrid(target=phi_sm$ID$, source=phiTmp_s$ID$)\n\
     \n\
     # create surface\n\
     pp_sm$ID$.readParticles(pp_s$ID$)\n\
@@ -342,30 +348,36 @@ def liquid_step_particles_$ID$():\n\
     \n\
     # no upres: just use the loaded grids\n\
     if upres_sp$ID$ <= 1:\n\
-        flags_sp$ID$.copyFrom(flags_s$ID$)\n\
-        vel_sp$ID$.copyFrom(vel_s$ID$)\n\
+        vel_sp$ID$.copyFrom(velTmp_s$ID$)\n\
         phiObs_sp$ID$.copyFrom(phiObs_s$ID$)\n\
-        phi_sp$ID$.copyFrom(phi_s$ID$)\n\
+        phi_sp$ID$.copyFrom(phiTmp_s$ID$)\n\
+        phiOut_sp$ID$.copyFrom(phiOut_s$ID$)\n\
     \n\
     # with upres: recreate grids\n\
     else:\n\
         # create highres grids by interpolation\n\
-        interpolateMACGrid(target=vel_sp$ID$, source=vel_s$ID$)\n\
-        interpolateGrid(target=phi_sp$ID$, source=phi_s$ID$)\n\
-        flags_sp$ID$.initDomain(boundaryWidth=boundaryWidth_s$ID$, phiWalls=phiObs_sp$ID$, outflow=boundConditions_s$ID$)\n\
-        flags_sp$ID$.updateFromLevelset(levelset=phi_sp$ID$)\n\
+        interpolateMACGrid(target=vel_sp$ID$, source=velTmp_s$ID$)\n\
+        interpolateGrid(target=phiObs_sp$ID$, source=phiObs_s$ID$)\n\
+        interpolateGrid(target=phi_sp$ID$, source=phiTmp_s$ID$)\n\
+        interpolateGrid(target=phiOut_sp$ID$, source=phiOut_s$ID$)\n\
+    \n\
+    flags_sp$ID$.initDomain(boundaryWidth=1 if using_fractions_s$ID$ else 0, phiWalls=phiObs_sp$ID$, outflow=boundConditions_s$ID$)\n\
+    setObstacleFlags(flags=flags_sp$ID$, phiObs=phiObs_sp$ID$, phiOut=None, phiIn=None) # phiIn not needed\n\
+    flags_sp$ID$.updateFromLevelset(levelset=phi_sp$ID$)\n\
     \n\
-    # actual secondary simulation\n\
-    #extrapolateLsSimple(phi=phi_sp$ID$, distance=radius+1, inside=True)\n\
+    # Actual secondary particle simulation\n\
     flipComputeSecondaryParticlePotentials(potTA=trappedAir_sp$ID$, potWC=waveCrest_sp$ID$, potKE=kineticEnergy_sp$ID$, neighborRatio=neighborRatio_sp$ID$, flags=flags_sp$ID$, v=vel_sp$ID$, normal=normal_sp$ID$, phi=phi_sp$ID$, radius=pot_radius_sp$ID$, tauMinTA=tauMin_ta_sp$ID$, tauMaxTA=tauMax_ta_sp$ID$, tauMinWC=tauMin_wc_sp$ID$, tauMaxWC=tauMax_wc_sp$ID$, tauMinKE=tauMin_k_sp$ID$, tauMaxKE=tauMax_k_sp$ID$, scaleFromManta=scaleFromManta_sp$ID$)\n\
-    flipSampleSecondaryParticles(mode='single', flags=flags_sp$ID$, v=vel_sp$ID$, pts_sec=ppSnd_sp$ID$, v_sec=pVelSnd_pp$ID$, l_sec=pLifeSnd_pp$ID$, lMin=lMin_sp$ID$, lMax=lMax_sp$ID$, potTA=trappedAir_sp$ID$, potWC=waveCrest_sp$ID$, potKE=kineticEnergy_sp$ID$, neighborRatio=neighborRatio_sp$ID$, c_s=c_s_sp$ID$, c_b=c_b_sp$ID$, k_ta=k_ta_sp$ID$, k_wc=k_wc_sp$ID$, dt=s$ID$.frameLength)\n\
-    flipUpdateSecondaryParticles(mode='linear', pts_sec=ppSnd_sp$ID$, v_sec=pVelSnd_pp$ID$, l_sec=pLifeSnd_pp$ID$, f_sec=pForceSnd_pp$ID$, flags=flags_sp$ID$, v=vel_sp$ID$, neighborRatio=neighborRatio_sp$ID$, radius=update_radius_sp$ID$, gravity=gravity_s$ID$, k_b=k_b_sp$ID$, k_d=k_d_sp$ID$, c_s=c_s_sp$ID$, c_b=c_b_sp$ID$, dt=s$ID$.frameLength)\n\
+    flipSampleSecondaryParticles(mode='single', flags=flags_sp$ID$, v=vel_sp$ID$, pts_sec=ppSnd_sp$ID$, v_sec=pVelSnd_pp$ID$, l_sec=pLifeSnd_pp$ID$, lMin=lMin_sp$ID$, lMax=lMax_sp$ID$, potTA=trappedAir_sp$ID$, potWC=waveCrest_sp$ID$, potKE=kineticEnergy_sp$ID$, neighborRatio=neighborRatio_sp$ID$, c_s=c_s_sp$ID$, c_b=c_b_sp$ID$, k_ta=k_ta_sp$ID$, k_wc=k_wc_sp$ID$, dt=sp$ID$.timestep)\n\
+    flipUpdateSecondaryParticles(mode='linear', pts_sec=ppSnd_sp$ID$, v_sec=pVelSnd_pp$ID$, l_sec=pLifeSnd_pp$ID$, f_sec=pForceSnd_pp$ID$, flags=flags_sp$ID$, v=vel_sp$ID$, neighborRatio=neighborRatio_sp$ID$, radius=update_radius_sp$ID$, gravity=gravity_s$ID$, k_b=k_b_sp$ID$, k_d=k_d_sp$ID$, c_s=c_s_sp$ID$, c_b=c_b_sp$ID$, dt=sp$ID$.timestep)\n\
     if $SNDPARTICLE_BOUNDARY_PUSHOUT$:\n\
-        pushOutofObs(parts = ppSnd_sp$ID$, flags = flags_sp$ID$, phiObs = phiObs_sp$ID$, shift = 1.0)\n\
-    flipDeleteParticlesInObstacle(pts=ppSnd_sp$ID$, flags=flags_sp$ID$)\n\
-    #debugGridInfo(flags = flags_sp$ID$, grid = trappedAir_sp$ID$, name = 'Trapped Air')\n\
-    #debugGridInfo(flags = flags_sp$ID$, grid = waveCrest_sp$ID$, name = 'Wave Crest')\n\
-    #debugGridInfo(flags = flags_sp$ID$, grid = kineticEnergy_sp$ID$, name = 'Kinetic Energy')\n";
+        pushOutofObs(parts=ppSnd_sp$ID$, flags=flags_sp$ID$, phiObs=phiObs_sp$ID$, shift=1.0)\n\
+    flipDeleteParticlesInObstacle(pts=ppSnd_sp$ID$, flags=flags_sp$ID$) # delete particles inside obstacle and outflow cells\n\
+    \n\
+    # Print debug information in the console\n\
+    if 0:\n\
+        debugGridInfo(flags=flags_sp$ID$, grid=trappedAir_sp$ID$, name='Trapped Air')\n\
+        debugGridInfo(flags=flags_sp$ID$, grid=waveCrest_sp$ID$, name='Wave Crest')\n\
+        debugGridInfo(flags=flags_sp$ID$, grid=kineticEnergy_sp$ID$, name='Kinetic Energy')\n";
 
 //////////////////////////////////////////////////////////////////////
 // IMPORT