Aligned repo with sandbox

doc/centered

@@ -169,7 +169,7 @@ enforce the divergence-free condition on the predicted face velocity. To do so,
 we solve the Poisson equation:
 
 $$
-  \left(\dfrac{1}{\rho^{n+1/2}}\nabla p^{n+1/2}\right)
+  \nabla\cdot\left(\dfrac{1}{\rho^{n+1/2}}\nabla p^{n+1/2}\right)
   = \dfrac{\nabla\cdot\mathbf{u}_{p,f}^{n+1/2}}{\Delta t/2}
 $$
 
@@ -316,7 +316,7 @@ The remaining part of the momentum equation (with just the pressure gradient),
 and the continuity equation are combined in a Poisson equation:
 
 $$
-  \left(\dfrac{1}{\rho^{n+1/2}}\nabla p^{n+1}\right) =
+  \nabla\cdot\left(\dfrac{1}{\rho^{n+1/2}}\nabla p^{n+1}\right) =
   \dfrac{\nabla\cdot\mathbf{u}_f^{**}}{\Delta t}
 $$
 
@@ -353,7 +353,7 @@ The collocated $\mathbf{g}^{n+1}$ is updated from $\mathbf{g}_f^{n+1}$ using a
 linear interpolation:
 
 $$
-  \mathbf{g}^{n+1} = \dfrac{\mathbf{g}_f^{n+1}[1] + \mathbf{g}_f^{n+1}}{2}
+  \mathbf{g}^{n+1} = \dfrac{\mathbf{g}_f^{n+1}[1] + \mathbf{g}_f^{n+1}[]}{2}
 $$
 
 Finally, the collocated velocity at the end of the time step is obtained:

run/coalescence.c

@@ -0,0 +1,100 @@
+/**
+# Droplet Coalescence with Phase Change
+
+We modify the [coalescence.c](/sandbox/popinet/coalescence.c) file using the clsvof approach and including phase change with a simple isothermal evaporation model.
+
+![Evolution of the velocity field](coalescence/movie.mp4)
+
+![Evolution of the species mass fractions](coalescence/species.mp4)
+*/
+
+#include "grid/multigrid.h"
+#include "navier-stokes/centered-evaporation.h"
+#include "navier-stokes/centered-doubled.h"
+#include "two-phase-clsvof.h"
+#include "integral.h"
+#include "evaporation.h"
+#include "species-gradient.h"
+#include "view.h"
+
+double Dmix1, Dmix2, YIntVal;
+
+u.n[top] = neumann (0.);
+u.t[top] = neumann (0.);
+p[top] = dirichlet (0.);
+uext.n[top] = neumann (0.);
+uext.t[top] = neumann (0.);
+pext[top] = dirichlet (0.);
+
+u.n[bottom] = neumann (0.);
+u.t[bottom] = neumann (0.);
+p[bottom] = dirichlet (0.);
+uext.n[bottom] = neumann (0.);
+uext.t[bottom] = neumann (0.);
+pext[bottom] = dirichlet (0.);
+
+u.n[left] = neumann (0.);
+u.t[left] = neumann (0.);
+p[left] = dirichlet (0.);
+uext.n[left] = neumann (0.);
+uext.t[left] = neumann (0.);
+pext[left] = dirichlet (0.);
+
+u.n[right] = neumann (0.);
+u.t[right] = neumann (0.);
+p[right] = dirichlet (0.);
+uext.n[right] = neumann (0.);
+uext.t[right] = neumann (0.);
+pext[right] = dirichlet (0.);
+
+int main()
+{
+  size (4. [*]);
+  origin (-L0/2. [*], -L0/2. [*]);
+
+  rho1 = 1. [*], rho2 = 1. [*];
+  mu1 = 1e-2 [*], mu2 = 1e-2 [*];
+  Dmix1 = 0., Dmix2 = 1.e-2;
+  YIntVal = 0.6;
+
+  const scalar sigmav[] = 1.;
+  d.sigmaf = sigmav;
+
+  run();
+}
+
+event init (t = 0)
+{
+  fraction (f, max (- (sq(x + 1.) + sq(y) - sq(0.5)),
+		    - (sq(x - 1.) + sq(y) - sq(0.5))));
+  foreach()
+    d[] = max (- (sq(x + 1.) + sq(y) - sq(0.5)),
+        - (sq(x - 1.) + sq(y) - sq(0.5)));
+
+  foreach() {
+    u.x[] = - sign(x)*f[];
+    uext.x[] = u.x[];
+  }
+
+  foreach() {
+    YL[] = f[];
+    YG[] = 0.;
+    Y[] = YL[] + YG[];
+  }
+}
+
+event movie (t += 0.04; t <= 6.)
+{
+  clear();
+  squares ("u.x", spread = -1, linear = true);
+  draw_vof ("f");
+  box();
+  save ("movie.mp4");
+
+  clear();
+  squares ("Y", min = 0., max = 1., linear = true);
+  draw_vof ("f");
+  box();
+  save ("species.mp4");
+}
+

run/dimensionalteq.c

@@ -46,18 +46,12 @@ Problably it is not necessary to set the units for all these
 quantities, but it seems easier to me to just set the units of
 every input value. */
 
-//double T0 = 273. SI_Temp;
-//double TL = 373. SI_Temp;
-//double TR = 173. SI_Temp;
-//double rho = 1.  SI_Density;
-//double cp = 1.   SI_SpecificHeat;
-//double k = 1.e-2 SI_ThermalCond;
-double T0 = 273.  SI_Temp;
-double TL = 373.  SI_Temp;
-double TR = 173.  SI_Temp;
-double rho = 1.   SI_Density;
-double cp = 1.    SI_SpecificHeat;
-double k = 1.e-2  SI_ThermalCond;
+double T0 = 273. SI_Temp;
+double TL = 373. SI_Temp;
+double TR = 173. SI_Temp;
+double rho = 1.  SI_Density;
+double cp = 1.   SI_SpecificHeat;
+double k = 1.e-2 SI_ThermalCond;
 
 /**
 We create the temperature field and the thermal diffusivity,

run/embedfalldrop.c

@@ -88,7 +88,7 @@ event movie (t += 0.001; t <= tEnd) {
   clear();
   box();
   view (tx = -0.5, ty = -0.5);
-  draw_vof ("cs", "fs", filled = -1, fc = {1.,1.,1.});
+  draw_vof ("cs", "fs", filled = -1, fc = {0.3,0.3,0.3});
   draw_vof ("f", lw = 1.5);
   squares ("omega", linear = false, 
            min = -1000., max = 1000., 

run/epsteinplesset.c

@@ -12,10 +12,7 @@ Therefore, the Stefan flow is not considered in the theretical
 solution, obtained assuming quasi-static conditions. In this context,
 the chemical species concentration profile is assumed to be equal to
 the steady-state concentration profile, at any simulation time instants.
-This approximation is justified when the diffusion is much faster than
-the interface regression velocity. Using this approximation, the
-theretical solution describing the evolution of the droplet radius in
-time was obtained by [Epstein and Plesset in 1950](#epstein1950stability):
+This approximation is justified when the diffusion is much faster than the interface regression velocity. Using this approximation, the theretical solution describing the evolution of the droplet radius in time was obtained by [Epstein and Plesset in 1950](#epstein1950stability):
 
 $$
   \dfrac{dR}{dt} = -MW\dfrac{\mathcal{D} (\hat{c} - c_{bulk})}{\rho_g}

run/pureisothermal.c

@@ -289,7 +289,7 @@ event movie (t += 0.5e-5, t <= 1.5e-4) {
 }
 
 /**
-## References
+## Results
 
 The numerical results are compared with the results obtained
 by [Pathak et al., 2018](#pathak2018steady) using a radially

run/radialdiffusion.c

@@ -5,7 +5,7 @@ a 1D grid and the ([spherisym.h](/sandbox/ecipriano/src/spherisym.h)) metrics.
 */
 
 #include "grid/multigrid1D.h"
-#include "../src/spherisym.h"
+#include "spherisym.h"
 #include "run.h"
 #include "timestep.h"
 #include "diffusion.h"
@@ -150,4 +150,4 @@ p "outfile-12" index (STATS_blocks-2) u 1:3 w l t "Analytic", \
   "outfile-11" index (STATS_blocks-2) w l t "maxlevel 11", \
   "outfile-12" index (STATS_blocks-2) w l t "maxlevel 12"
 ~~~
-*/
+*/

run/suckingproblem.c

@@ -183,7 +183,7 @@ int main (void) {
   coefficient. */
 
   L0 = 10.e-3;
-  f.sigma = 0.059;
+  f.sigma = 0.0059;
 
   /**
   We run the simulation at four different levels of refinement. */

src/common-evaporation.h

@@ -14,7 +14,6 @@ for evaporation models. */
 #ifndef T_ERR
 # define T_ERR 0.
 #endif
-//#include "mass_refine_prolongation.h"
 
 #ifndef I_TOL
 # define I_TOL 1.e-8

src/navier-stokes/velocity-jump.h

@@ -23,35 +23,45 @@ The scheme implemented here is close to that used in Gerris ([Popinet,
 We will use the generic time loop, a CFL-limited timestep, the
 Bell-Collela-Glaz advection scheme and the implicit viscosity
 solver. If embedded boundaries are used, a different scheme is used
-for viscosity. This scheme is extended considering the phase change
-between two phases described by a VOF or CLSVOF approach. Using the
-jump condition is convenient because it does not require an explicit
-source term, localized at the gas-liquid interface, in the projection
-step, which causes oscillations in the velocity field.
-
-The method that we use here to enforce the velocity jump condition is
-the same proposed by [Tanguy et al. 2007](#tanguy2007level) for droplet
+for viscosity.
+
+We extend this scheme by including the velocity jump due to the phase
+change:
+$$
+  \left(\mathbf{u}_g - \mathbf{u}_l\right)\cdot\mathbf{n}_\Gamma
+  = \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)
+$$
+which is introduced into the velocity field directly, according to the
+ghost velocity approach proposed by [Nguyen et al. 2001](#nguyen2001boundary),
+that was modified by [Tanguy et al. 2007](#tanguy2007level) for droplet
 evaporation problems, and by [Tanguy et al. 2014](#tanguy2014benchmarks)
-for boiling simulations. It is based on the combination of the ghost fluid
-velocity method poposed by [Nguyen et al. 2001](#nguyen2001boundary),
-which defines ghost velocities as:
+for boiling simulations.
 
+The idea is to set the velocity jump in the "ghost regions" by exploiting the
+jump condition:
 $$
   \mathbf{u}^{ghost}_l = \mathbf{u}_g -
-  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}
+  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}_\Gamma
 $$
 
 $$
   \mathbf{u}^{ghost}_g = \mathbf{u}_l +
-  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}
+  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}_\Gamma
 $$
-
-The ghost velocities impose the continuity of the velocity fields
-across the interface. Two different momentum equations for the gas
-and liquid phase velocities are solved, a single projection step is
+Two different momentum equations for the gas
+and liquid phase velocities are solved, but a single projection step is
 used to update the velocities at the new time step. Additional
-velocity extensions are used to obtain a divergence-free velocity for
+projection steps are used to obtain a divergence-free velocity for
 the advection of the volume fraction field.
+
+This is a modified version of the method discussed by [Long et al. 2024](#long2024edge)
+for phase change using the EBIT approach, and implemented in basilisk by [Tian Long](/sandbox/tianlong/src/semushin/double-evaporation.h).
+Unlike the previous version, this approach has mainly been tested on evaporation
+problems rather than on boiling simulations, and it was designed for VOF simulations.
+
+The main advantage with respect to the [navier-stokes/centered-doubled.h](/sandbox/ecipriano/src/navier-stokes/centered-doubled.h)
+approach is that we avoid the oscillations that arise due to the introduction of
+the localized expansion term as a volumetric source in the projection step.
 */
 
 #define VELOCITY_JUMP
@@ -580,7 +590,7 @@ gas phase ghost velocity:
 
 $$
   \mathbf{u}^{ghost}_g = \mathbf{u}_l -
-  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}
+  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}_\Gamma
 $$
 
 while the initial liquid ghost velocity is set to zero. For boiling
@@ -590,7 +600,7 @@ the jump condition to set the value of the liquid phase ghost
 
 $$
   \mathbf{u}^{ghost}_l = \mathbf{u}_g +
-  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}
+  \dot{m}\left(\dfrac{1}{\rho_g} - \dfrac{1}{\rho_l}\right)\mathbf{n}_\Gamma
 $$
 */
 
@@ -600,13 +610,9 @@ void update_ghost_velocities (void) {
 #ifdef BOILING_SETUP
       u2g.x[] = u2g.x[];
       u1g.x[] = u2.x[] + mEvapTot1[]*n.x[];
-
-      //u2g.x[] = u1.x[] - mEvapTot2[]*n.x[]; // [Test]
 #else
       u1g.x[] = u1g.x[];
       u2g.x[] = u1.x[] - mEvapTot2[]*n.x[];
-
-      //u1g.x[] = u2.x[] + mEvapTot1[]*n.x[]; // [Test]
 #endif
       u1.x[] = (ls1[] < 0.) ? u1.x[] : u1g.x[];
       u2.x[] = (ls1[] > 0.) ? u2.x[] : u2g.x[];
@@ -618,16 +624,10 @@ void update_ghost_velocities (void) {
     double mEvapTot1f = 0.5*(mEvapTot1[] + mEvapTot1[-1]);
     uf2g.x[] = uf2g.x[];
     uf1g.x[] = uf2.x[] + mEvapTot1f*nf.x[]*fm.x[];
-
-    //double mEvapTot2f = 0.5*(mEvapTot2[] + mEvapTot2[-1]); // [Test]
-    //uf2g.x[] = uf1.x[] - mEvapTot2f*nf.x[]*fm.x[];             // [Test]
 #else
     double mEvapTot2f = 0.5*(mEvapTot2[] + mEvapTot2[-1]);
     uf1g.x[] = uf1g.x[];
     uf2g.x[] = uf1.x[] - mEvapTot2f*nf.x[]*fm.x[];
-
-    //double mEvapTot1f = 0.5*(mEvapTot1[] + mEvapTot1[-1]); // [Test]
-    //uf1g.x[] = uf2.x[] + mEvapTot1f*nf.x[]*fm.x[];             // [Test]
 #endif
     double lsf = 0.5*(ls1[] + ls1[-1]);
     uf1.x[] = (lsf < 0.) ? uf1.x[] : uf1g.x[];
@@ -1017,13 +1017,9 @@ event update_ghost (i++, last) {
 #ifdef BOILING_SETUP
       u2g.x[] = u2g.x[];
       u1g.x[] = u2.x[] + mEvapTot1[]*n.x[];
-
-      //u2g.x[] = u1.x[] - mEvapTot2[]*n.x[]; // [Test]
 #else
       u1g.x[] = u1g.x[];
       u2g.x[] = u1.x[] - mEvapTot2[]*n.x[];
-
-      //u1g.x[] = u2.x[] + mEvapTot1[]*n.x[]; // [Test]
 #endif
     }
   }
@@ -1033,16 +1029,10 @@ event update_ghost (i++, last) {
     double mEvapTot1f = 0.5*(mEvapTot1[] + mEvapTot1[-1]);
     uf2g.x[] = uf2g.x[];
     uf1g.x[] = uf2.x[] + mEvapTot1f*nf.x[]*fm.x[];
-
-    //double mEvapTot2f = 0.5*(mEvapTot2[] + mEvapTot2[-1]); // [Test]
-    //uf2g.x[] = uf1.x[] - mEvapTot2f*nf.x[]*fm.x[];             // [Test]
 #else
     double mEvapTot2f = 0.5*(mEvapTot2[] + mEvapTot2[-1]);
     uf1g.x[] = uf1g.x[];
     uf2g.x[] = uf1.x[] - mEvapTot2f*nf.x[]*fm.x[];
-
-    //double mEvapTot1f = 0.5*(mEvapTot1[] + mEvapTot1[-1]); // [Test]
-    //uf1g.x[] = uf2.x[] + mEvapTot1f*nf.x[]*fm.x[];             // [Test]
 #endif
   }
 }
@@ -1201,6 +1191,13 @@ event adapt (i++,last) {
   year={2001},
   publisher={Elsevier}
 }
+
+@article{long2024edge,
+  title={An Edge-based Interface Tracking (EBIT) Method for Multiphase Flows with Phase Change},
+  author={Long, Tian and Pan, Jieyun and Zaleski, St{\'e}phane},
+  journal={arXiv preprint arXiv:2402.13677},
+  year={2024}
+}
 ~~~
 */