sandbox/alimare/LS_advection.h
Advection of a level-set function
We will use a RK3 time integration scheme for collocated variables. With its own timestep (no restriction due to the presence of an embedded boundary).
#include "simple_discretization.h"
For this function we need the following files.
#include "LS_reinit.h"
#include "../ghigo/src/myquadratic.h" // best extrapolation so far
If we solve the Navier-Stokes equations, we use ghigo’s method to update the variables in emerged cells.
#if NS_emerged
#include "emerged_NS.h"
#endif
struct LSadv{
scalar dist;
scalar cs;
face vector fs;
scalar TS;
scalar TL;
vector vpcf;
int itredist;
double s_clean;
double NB_width;
scalar curve;
};
void advection_LS(struct LSadv p){
scalar dist = p.dist;
scalar cs = p.cs;
face vector fs = p.fs;
scalar TS = p.TS;
scalar TL = p.TL;
vector vpcf = p.vpcf;
int itredist = p.itredist;
double s_clean = p.s_clean;
double NB_width = p.NB_width;
scalar curve = p.curve;
Previous state of cs is saved into csm1
foreach(){
csm1[] = cs[];
}
boundary({csm1});
restriction({csm1});
face vector fsm1[];
foreach_face()
fsm1.x[] = fs.x[];
boundary((scalar *){fsm1});
restriction((scalar *){fsm1});
RK3_WENO5(dist,vpcf,dt, NB_width);
boundary ({dist});
restriction({dist});
LS_reinit(dist, it_max = itredist);
We remove the overshoots that we might create with reinitialization.
foreach(){
dist[] = clamp(dist[], -1.05*NB_width, 1.05*NB_width);
}
boundary({dist});
restriction({dist});
We update the volume and face fractions accordingly
LS2fractions(dist,cs,fs,s_clean);
#if CURVE_LS // mandatory in Gibbs Thomsmon cases
curvature_LS(dist,curve);
#else
curvature(cs,curve);
#endif
boundary({curve});
restriction({curve});
Sometimes, a new cell becomes uncovered/covered because the interface has moved. We must initialize the tracers field. We use Ghigo’s methodology in update_tracer() that can be found in this page. The basic idea is very similar to what can be found in the Dirichlet boundary condition calculation in embed.h
.
#if 1
foreach() {
if (cs[] > 0. && csm1[] <= 0.) { // Emerged cells
assert(cs[]!=1.); // cell shouldn't be full
coord o = {0.,0.};
TL[] = embed_extrapolate_ls (point, TL, cs, o, false);
}
}
We update the boundary condition for a.
boundary({TL});
restriction({TL});
#if NS_emerged
init_emerged_NS(cs,csm1);
#endif
invertcs(cs,fs);
invertcs(csm1,fsm1);
foreach() {
if (cs[] > 0. && csm1[] <= 0.) { // Emerged cells
assert(cs[]!=1.); // cell shouldn't be full
coord o = {0.,0.};
TS[] = embed_extrapolate_ls(point, TS,cs, o, false);
}
}
boundary ({TS});
invertcs(cs,fs);
invertcs(csm1,fsm1);
#endif
}