aslam.h

Aslam Extrapolations

Aslam Extrapolations

This module implements the PDE-based Aslam extrapolation methods (Aslam 2003). Using these functions, a discontinuous field can be extended from the liquid phase to the gas-phase and vice-versa using constant and extrapolations.

These methods are expensive, since the PDE must be solved at steady-state. However, most of the times they are used just for a few steps in order to extend the field for a few cells across the interface, without covering the whole domain.

#ifndef ASLAM_H
#define ASLAM_H

#include "mapregion.h"
#include "redistance.h"
#include "fractions.h"

void vof_to_ls (scalar f, scalar ls, int imax = 3) {
#if 0
  double weight = 0.1;
  foreach()
    if (f[] > 1e-6 && f[] < 1. - 1e-6) {
      coord n = interface_normal (point, f);
      normalize (&n);
      double alpha = plane_alpha (f[], n);
      ls[] = (1. - weight)*ls[] + weight*Delta*alpha;
      ls[] *= -1;
    }

  redistance (ls, imax = imax);
#endif
  double deltamin = L0/(1 << grid->maxdepth);
  foreach()
    ls[] = -(2.*f[] - 1.)*deltamin*0.75;
#if TREE
  restriction({ls});
#endif
  redistance (ls, imax);
}

Constant Extrapolation

The constant extrapolation of a field f, existing only in a portion of space defined by ls, can be extrapolated in the whole domain solving the following PDE:

\dfrac{\partial f}{\partial t} + H(\phi)\hat{\mathbf{n}}\cdot\nabla f = 0

where H(\phi) is the Heaviside function, used to define the region where the field should be extrapolated and where it should not be modified:

H(\phi) = \begin{cases} 1 & \text{if } \phi > 0,\\ 0 & \text{if } \phi \leq 0. \end{cases}

while \hat{\mathbf{n}} is the unit normal.

trace
void constant_extrapolation (
  scalar f,                   // field to extrapolate
  scalar ls,                  // level set field
  double cfl = 0.5,           // CFL number
  int nmax,                   // number of maximum time steps
  (const) scalar s = {-1},    // source term, default zero
  (const) scalar c = {-1},    // vof field, optional
  int nl = 0,                 // from which layer of cells (optional, min=0, max=2)
  int nointerface = 0,        // remove the interface from the heaviside (default false)
  int inverse = 0,            // the vof field if = 1 if gas (default false)
  double tol = 1e-10          // tolerance which defines the interface
)
{
  scalar H[];
  vector n[], gf[];

  if (s.i < 0)
    s = zeroc;

We compute the gradients of the level set for the calcation of the interface normals.

  gradients ({ls}, {n});

We compute the normals and the heaviside function H, which is non-null in the region where the field must be extrapolated. In case of vof field, the user can decide to extrapolate the field from a layer of cells which in not adjacent to the interface. This can be specified setting the variables nl, 0 by default, it can be set to 1 or 2.

  if (c.i > 0)
    mapregion (H, c, nl=nl, nointerface=nointerface, inverse=inverse,
        tol=tol);
  else
    foreach()
      if (inverse)
        H[] = (-ls[] <= 0.) ? 0. : 1.;
      else
        H[] = (ls[] <= 0.) ? 0. : 1.;

  foreach() {
    double maggf = 0.;
    foreach_dimension()
      maggf += sq (n.x[]);
    maggf = sqrt (maggf);
    foreach_dimension()
      n.x[] /= inverse ? -(maggf + 1.e-10) : (maggf + 1.e-10);
  }

We solve the PDE inside a loop over the maximum number of time steps that, in principle, should ensure to arrive at steady-state. In practice, less time-steps will be sufficients to extrapolate the fields over the interface.

  int ts = 0;

  while (ts < nmax) {

The gradients of the extrapolated functions are updated using an upwind scheme.

    foreach()
      foreach_dimension()
        gf.x[] = (n.x[] <= 0.) ? (f[1] - f[])/Delta : (f[] - f[-1])/Delta;

We solve a single step of the PDE.

    foreach() {
      double dt = cfl*Delta;
      double nscalargf = 0.;
      foreach_dimension()
        nscalargf += n.x[]*gf.x[];
      f[] -= dt*H[]*(nscalargf - s[]);
    }
    ts++;
  }
}

Linear Extrapolation

The linear extrapolation of the field f along the normal direction, can be achieved from the solution of the following PDE:

\dfrac{\partial f}{\partial t} + H(\phi)\left(\hat{\mathbf{n}}\cdot\nabla f - f_n \right) = 0

Which is similar to the equation resolved with the constant extrapolation procedure, except for the source term f_n, which is the directional derivative of f in the normal direction, defined as:

f_n = \hat{\mathbf{n}}\cdot \nabla f

Since this function is defined only in the region where f is defined, we apply the constant extrapolation in order to obtain a field f_n defined on the whole domain:

\dfrac{\partial f_n}{\partial t} + H(\phi)\hat{\mathbf{n}}\cdot\nabla f_n = 0

Therefore, the linear interpolation requires the solution of two different extrapolation PDEs, and the same logic applies if higher order extrapolations have to be solved.

void linear_extrapolation (
  scalar f,                   // field to extrapolate
  scalar ls,                  // level set field
  double cfl,                 // CFL number
  int nmax,                   // number of maximum time steps
  (const) scalar s = {-1},    // source term, default zero
  (const) scalar c = {-1},    // vof field, optional
  int nl = 0,                 // from which layer of cells (optional, min=0, max=2)
  int nointerface = 0,        // remove the interface from the heaviside (default false)
  int inverse = 0,            // the vof field if = 1 if gas (default false)
  double tol = 1e-10          // tolerance which defines the interface
)
{
  scalar H[], fn[];
  vector n[], gf[];

  if (s.i < 0)
    s = zeroc;

We compute the gradients of the level set for the calcation of the interface normals, and the gradient of f for the calculation of the directional derivative.

  gradients ({ls, f}, {n, gf});

  if (c.i > 0)
    mapregion (H, c, nl = (nl == 0.) ? 1. : min (nl+1, 2),
        nointerface=nointerface, inverse=inverse, tol=tol);
  else
    foreach()
      if (inverse)
        H[] = (ls[]+Delta <= 0.) ? 0. : 1.;
      else
        H[] = (-ls[]+Delta <= 0.) ? 0. : 1.;

  foreach() {
    double maggf = 0.;
    foreach_dimension()
      maggf += sq (n.x[]);
    maggf = sqrt (maggf);
    foreach_dimension()
      n.x[] /= inverse ? -(maggf + 1.e-10) : (maggf + 1.e-10);
  }

We compute the directional derviative fn.

  foreach()
    foreach_dimension()
      fn[] += n.x[]*gf.x[];

We solve the constant extrapolation for extending the directional derivative.

  int ts = 0;

  while (ts < nmax) {

The gradients of the extrapolated functions are updated using an upwind scheme.

    foreach()
      foreach_dimension()
        gf.x[] = (n.x[] <= 0.) ? (fn[1] - fn[])/Delta : (fn[] - fn[-1])/Delta;

We solve a single step of the PDE.

    foreach() {
      double dt = cfl*Delta;
      double nscalargf = 0.;
      foreach_dimension()
        nscalargf += n.x[]*gf.x[];
      fn[] -= dt*H[]*(nscalargf - s[]);
    }
    ts++;
  }

We solve a constant extrapolation with source equal to the directional derivative.

  if (c.i > 0)
    constant_extrapolation (f, ls, cfl, nmax, fn, c, nl, nointerface, tol=tol);
  else
    constant_extrapolation (f, ls, cfl, nmax, fn, tol=tol);
}

#endif

References

[aslam2004partial]

Tariq D Aslam. A partial differential equation approach to multidimensional extrapolation. Journal of Computational Physics, 193(1):349–355, 2004.