LS_advection.h

Advection of a level-set function

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){
  void advection_LS(scalar dist, scalar cs, face vector fs, scalar TS, scalar TL,
                    vector vpcf, int itredist, double s_clean, double NB_width, scalar curve){

  // 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 = dist, dt = 0.2 * L0/(1 << grid->maxdepth), 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[
  

//

// We remove the overshoots that we might create with reinitialization. //

//   foreach(){
//     dist[] = clamp(dist[], -1.05*NB_width, 1.05*NB_width);
//   }
//   bounda]!=1.); // cell shouldn't be full
//       if (cs[] >= 1.)
//     continue; 
//      coord o = {0.,0.};
//       TL[] = embed_extrapolate_ls (point, TL, cs, o, false);
//     }
//   }

#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
}