simple_discretization.h

Advection solvers
- Timestep
- Advection

Advection solvers

These functions are only redefinition of other functions (fluxes and advection) in order to ignore the presence of an embedded boundary.

#include "alex_functions.h"

Timestep

We need to define a timestep on the mesh without taking into accound the embedded boundary. We only take into account cells that will be updated, i.e. such that |dist[]| < maxd. We also don’t want the previous variable to be initialized at 0 for simple test cases, we want the timestep to be constant throughout the simulation.

double timestep_LS (vector u, double dtmax,
  scalar dist, double maxd)
{
  dtmax /= CFL;
  foreach(reduction(min:dtmax))
    foreach_dimension(){
      if (u.x[] != 0. && fabs(dist[]) < maxd) {
        double dt = Delta/fabs(u.x[]);
        if (dt < dtmax) dtmax = dt;
      }
    }

  return dtmax *= CFL;
}

Advection

We use a simple advection scheme with co-located variables.

#include "weno2.h"

void FE(scalar f, scalar fi, vector ndist, double dt, 
  double maxd){

Simple Forward Euler scheme.

  foreach(){
    if(fabs(f[]) < maxd)
      foreach_dimension(){
        double graplus  = WENOdiff_x(point, fi,-1);
        double graminus = WENOdiff_x(point, fi,1);
        f[] -= dt * ( max(ndist.x[],0.)*graplus
          + min(ndist.x[],0.)*graminus);
    }
  }
  boundary({f});
  restriction({f});
}

Runge Kutta 2 scheme

void RK2(scalar f, vector ndist, double dt, double maxd){
  scalar f1[],f2[];

  foreach(){
    f1[] = f[];
  }
  boundary({f1});
  restriction({f1});

  FE(f1, f, ndist, dt, maxd);
  
  foreach(){
    f2[] = f1[];
  }
  boundary({f2});
  restriction({f2});

  FE(f2, f1, ndist, dt, maxd);
  foreach(){
    f[] = (f2[] + f[])*0.5;
  }
  boundary({f});
  restriction({f});
}

Runge Kutta 3 advection scheme

void RK3(scalar f, vector ndist, double dt, double maxd){
  static const double coeff[2][2] = {
    {0.75, 0.25},
    {1./3.,2./3.}
  };
  scalar f1[],f2[];

  foreach(){
    f1[] = f[];
  }
  boundary({f1});
  restriction({f1});

  FE(f1, f, ndist, dt, maxd);
  foreach(){
    f2[] = f1[];
  }
  boundary({f2});
  restriction({f2});

  FE(f2, f1, ndist, dt, maxd);
  foreach(){
    f1[] = coeff[0][0]*f[]+coeff[0][1]*f2[];
    f2[] = f1[];
  }
  boundary({f1,f2});
  restriction({f1,f2});
  FE(f2, f1, ndist, dt, maxd);
  foreach(){
    f[] = coeff[1][0]*f[] + coeff[1][1]*f2[];
  }
  boundary({f});
  restriction({f});

}

Runge Kutta 3 advection scheme

#include "weno5.h"
void RK3_WENO5(scalar f, vector ndist, double dt, double maxd){
  static const double coeff[2][2] = {
    {0.75, 0.25},
    {1./3.,2./3.}
  };
  scalar f1[],f2[];

  foreach(){
    f1[] = f[];
  }
  boundary({f1});
  restriction({f1});

  FE_WENO5(f1, f, ndist, dt, maxd);
  foreach(){
    f2[] = f1[];
  }
  boundary({f2});
  restriction({f2});

  FE_WENO5(f2, f1, ndist, dt, maxd);
  foreach(){
    f1[] = coeff[0][0]*f[]+coeff[0][1]*f2[];
    f2[] = f1[];
  }
  boundary({f1,f2});
  restriction({f1,f2});
  FE_WENO5(f2, f1, ndist, dt, maxd);
  foreach(){
    f[] = coeff[1][0]*f[] + coeff[1][1]*f2[];
  }
  boundary({f});
  restriction({f});

}

void RK3_WENO3(scalar f, vector ndist, double dt, double maxd){
  static const double coeff[2][2] = {
    {0.75, 0.25},
    {1./3.,2./3.}
  };
  scalar f1[],f2[];

  foreach(){
    f1[] = f[];
  }
  boundary({f1});
  restriction({f1});

  FE_WENO3(f1, f, ndist, dt, maxd);
  foreach(){
    f2[] = f1[];
  }
  boundary({f2});
  restriction({f2});

  FE_WENO3(f2, f1, ndist, dt, maxd);
  foreach(){
    f1[] = coeff[0][0]*f[]+coeff[0][1]*f2[];
    f2[] = f1[];
  }
  boundary({f1,f2});
  restriction({f1,f2});
  FE_WENO3(f2, f1, ndist, dt, maxd);
  foreach(){
    f[] = coeff[1][0]*f[] + coeff[1][1]*f2[];
  }
  boundary({f});
  restriction({f});

}