redistance.h

📄 View in API Reference (Doxygen) · View on basilisk.fr

Redistancing of a distance field

The original implementation was by Alexandre Limare used in particular in Limare et al., 2022.

This file implements redistancing of a distance field with subcell correction, see the work of Russo & Smereka, 2000 with corrections by Min & Gibou, 2007 and by Min, 2010.

Let $\phi$ be a function close to a signed function that has been perturbed by numerical diffusion (more precisely, a non-zero tangential velocity). By iterating on the eikonal equation

\left\{\begin{array}{ll} \phi_t + \text{sign}(\phi^{0}) \left(\left| \nabla \phi\right| - 1 \right) = 0\\ \phi(x,0) = \phi^0(x) \end{array} \right.

we can correct or redistance $\phi$ to make it a signed function.

We use a Godunov Hamiltonian approximation for $\left| \nabla \phi\right|$

\left| \nabla \phi \right|_{ij} = H_G(D_x^+\phi_{ij}, D_x^-\phi_{ij}, D_y^+\phi_{ij}, D_y^-\phi_{ij})

where $D^\pm\phi_{ij}$ denotes the one-sided ENO difference finite difference in the x- direction

D_x^+ = \dfrac{\phi_{i+1,j}-\phi_{i,j}}{\Delta} - \dfrac{\Delta}{2}\text{minmod}(D_{xx}\phi_{ij}, D_{xx}\phi_{i+1,j})

D_x^- = \dfrac{\phi_{i,j}-\phi_{i-1,j}}{\Delta} + \dfrac{\Delta}{2}\text{minmod}(D_{xx}\phi_{ij}, D_{xx}\phi_{i+1,j})

here $D_{xx}\phi_{ij} = (\phi_{i-1,j} - 2\phi{ij} + \phi_{i+1,j})/\Delta^2$.

The minmod function is zero when the two arguments have different signs, and takes the argument with smaller absolute value when the two have the same sign.

The Godunov Hamiltonian $H_G$ is given as

H_G(a,b,c,d) = \left\{ \begin{array}{ll} \sqrt{\text{max}((a^-)^2,(b^+)^2 + (c^-)^2,(d^+)^2)} \text { when } \text{sign}(\phi^0_{ij}) \geq 0\\ \sqrt{\text{max}((a^+)^2,(b^-)^2 + (c^+)^2,(d^-)^2)} \text { when } \text{sign}(\phi^0_{ij}) < 0 \end{array} \right.

with

x^+ = \text{max}(0, x)\\ x^- = \text{min}(0, x)\\

We use a minmod limiter.

static inline double minmod3 (double a, double b)
{
  if (a == b || a*b <= 0.)
    return 0.;
  return fabs(a) < fabs(b) ? a : b;
}

#define BGHOSTS 2

Time derivative

This function fills *dphi* with $\partial_t\phi$ .

static
void dphidt (scalar phi, scalar dphi, scalar phi0, const double cfl, const double phixxmin)
{
  foreach() {
    double dt = cfl*Delta;

We first calculate the inputs of the Hamiltonian

D^+\phi_{ij} = \dfrac{\phi_{i+1,j} - \phi_{ij}}{\Delta x} - \dfrac{\Delta} {2}\text{minmod}(D_{xx}\phi_{ij},D_{xx}\phi_{i+1,j})

D^-\phi_{ij} = \dfrac{\phi_{i,j} - \phi_{i-1,j}}{\Delta x} - \dfrac{\Delta} {2}\text{minmod}(D_{xx}\phi_{ij},D_{xx}\phi{i-1,j})

where

D_{xx}\phi_{ij} = \dfrac{\phi_{i-1,j} - 2\phi_{ij} + \phi_{i+1,j}}{\Delta x^2}

coord gra[2];
    for (int i = 0, j = 1; i < 2; j = 1 - 2*++i)
      foreach_dimension() {
	double s1 = (phi[2*j] + phi[] - 2.*phi[j])/Delta; 
	double s2 = (phi[1] + phi[-1] - 2.*phi[])/Delta;
	gra[i].x = j*((phi[j] - phi[])/Delta - minmod3(s1, s2)/2.);
      }

We check for interfacial cells.

bool interfacial = false;
    foreach_dimension()
      if (phi0[-1]*phi0[] < 0 || phi0[1]*phi0[] < 0)
	interfacial = true;
    
    dphi[] = - sign2(phi0[]);
    
    if (interfacial) {

Near the interface, *i.e.* for cells where

\phi^0_i\phi^0_{i+1} \leq 0 \text{ or } \phi^0_i\phi^0_{i-1} \leq 0

The scheme must stay truly upwind, meaning that the movement of the 0 level-set of the function must be as small as possible. Therefore the upwind numerical scheme is modified to

D_x^+ = \dfrac{0-\phi_{ij}}{\Delta x^+} - \dfrac{\Delta x^+}{2} \text{minmod}(D_ {xx}\phi_{ij},D_{xx}\phi_{i+1,j}) \text{ if } \phi_{ij}\phi_{i+1,j} < 0

D_x^- = \dfrac{\phi_{ij}-0}{\Delta x^-} + \dfrac{\Delta x^-}{2} \text{minmod}(D_ {xx}\phi_{ij},D_{xx}\phi_{i-1,j}) \text{ if } \phi_{ij}\phi_{i+1,j} < 0

which is the correction by Min & Gibou 2007.

double size = HUGE;
      foreach_dimension()
	for (int i = 0, j = 1; i < 2; j = 1 - 2*++i)
	  if (phi0[]*phi0[j] < 0.) {

We compute the subcell fix near the interface (see section 2.2 in Min, 2010).

\Delta x^+ = \left\{ \begin{array}{ll} \Delta x \cdot \left( \frac{1}{2} + \dfrac{\phi^0_{i,j}-\phi^0_{i+1,j} - \text{sgn}(\phi^0_{i,j}-\phi^0_{i+1,j})\sqrt{D}}{}\right) \text{ if } \left| \phi^0_{xx}\right| >\epsilon \\ \Delta x \cdot \dfrac{\phi^0_{ij}}{\phi^0_{i,j}-\phi^0_{i+1,j}} \text{ else.} \\ \end{array} \right.

with

\phi_{xx}^0 = \text{minmod}(\phi^0_{i-1,j}-2\phi^0_{ij}+\phi^0_{i+1,j}, \phi^0_{i,j}-2\phi^0_{i+1j}+\phi^0_{i+2,j}) \\ D = \left( \phi^0_{xx}/2 - \phi_{ij}^0 - \phi_{i+1,j} \right)^2 - 4\phi_{ij}^0\phi_{i+1,j}^0

For the $\Delta x^-$ calculation, replace all the $+$ subscript by $-$, this is dealt with properly with the j parameter below.

double dx = Delta;
	    double phixx = minmod3 (phi0[2*j] + phi0[] - 2.*phi0[j],
				    phi0[1] + phi0[-1] - 2.*phi0[]);
	    if (fabs(phixx) > phixxmin) {
	      double D = sq(phixx/2. - phi0[] - phi0[j]) - 4.*phi0[]*phi0[j];
	      dx *= 1/2. + (phi0[] - phi0[j] - sign2(phi0[] - phi0[j])*sqrt(D))/phixx;
	    }
	    else
	      dx *= phi0[]/(phi0[] - phi0[j]);
	    
	    if (dx != 0.) {
	      double sxx1 = phi[2 - 4*i] + phi[] - 2.*phi[1 - 2*i];
	      double sxx2 = phi[1] + phi[-1] - 2.*phi[];
	      gra[i].x = (2*i - 1)*(phi[]/dx + dx*minmod3(sxx1, sxx2)/(2.*sq(Delta)));
	    }
	    else 
	      gra[i].x = 0.;
	    size = min(size, dx);
	  }
      dphi[] *= min(dt, fabs(size)/2.)/dt;
    }

The Godunov Hamiltonian is

H_G(a,b,c,d) = \left\{ \begin{array}{ll} \sqrt{\text{max}((a^-)^2,(b^+)^2 + (c^-)^2,(d^+)^2)} \text { when } \text{sign}(\phi^0_{ij}) \geq 0\\ \sqrt{\text{max}((a^+)^2,(b^-)^2 + (c^+)^2,(d^-)^2)} \text { when } \text{sign}(\phi^0_{ij}) < 0 \end{array} \right.

double H_G = 0;
    if (phi0[] > 0)
      foreach_dimension() {
	double a = min(0., gra[0].x); 
	double b = max(0., gra[1].x);
	H_G += max(sq(a), sq(b));
      }
    else
      foreach_dimension() {
	double a = max(0., gra[0].x);
	double b = min(0., gra[1].x);
	H_G += max(sq(a), sq(b));
      }
    
    dphi[] *= sqrt(H_G) - 1.;
  }
}

The redistance() function

trace
int redistance (scalar phi,
		int imax = 1,       // The maximum number of iterations
		double cfl = 0.5,   // The CFL number
		int order = 3,      // The order of time integration
		double eps = 1e-6,  // The maximum error on $|\nabla\phi| - 1$
		double band = HUGE, // The width of the band in which to compute the error
		scalar resf = {-1}, // The residual $|\nabla\phi| - 1$
		const double phixxmin = 1./HUGE) // The threshold for second-order subcell correction
{

We create phi0[], a copy of the initial level-set function before the iterations.

scalar phi0[];
  foreach()
    phi0[] = phi[] ;

Time integration iteration loop.

for (int i = 1; i <= imax; i++) {
    double maxres = 0.;

We use either a RK2 scheme...

if (order == 2) {
      scalar tmp1[];
      dphidt (phi, tmp1, phi0, cfl/2., phixxmin);
      foreach()
	tmp1[] = phi[] + Delta*cfl/2.*tmp1[];
      scalar tmp2[];
      dphidt (tmp1, tmp2, phi0, cfl, phixxmin);
      foreach (reduction(max:maxres)) {
	double res = tmp2[];
	if (resf.i >= 0) resf[] = res;
	if (fabs (res) > maxres && fabs(phi0[]) < band*Delta) maxres = fabs (res);
	phi[] += Delta*cfl*tmp2[];
      }
    }

... or a RK3 compact scheme from Shu and Osher, 1988.

else {
      scalar tmp1[];
      dphidt (phi, tmp1, phi0, cfl, phixxmin);
      foreach()
	tmp1[] = phi[] + Delta*cfl*tmp1[];
      scalar tmp2[];
      dphidt (tmp1, tmp2, phi0, cfl, phixxmin);
      foreach()
	tmp1[] = (3.*phi[] + tmp1[] + Delta*cfl*tmp2[])/4.;
      dphidt (tmp1, tmp2, phi0, cfl, phixxmin);
      foreach (reduction(max:maxres)) {
	double res = 2./3.*((phi[] - tmp1[])/(Delta*cfl) - tmp2[]);
	if (resf.i >= 0) resf[] = res;
	if (fabs (res) > maxres && fabs(phi0[]) < band*Delta) maxres = fabs (res);
	phi[] = (phi[] + 2.*(tmp1[] + Delta*cfl*tmp2[]))/3.;
      }
    }
    
    if (maxres < eps)
      return i;
  }
  return imax;
}

References

~~~bib @Article{shu1988, author = {Chi-Wang Shu and Stanley Osher}, title = {Efficient implementation of essentially non-oscillatory shock-capturing schemes}, journal = {Journal of Computational Physics}, year = {1988}, volume = {77}, pages = {439-471}, issn = {0021-9991}, doi = {10.1016/0021-9991(88)90177-5}, }

@article{russo2000, title = {A remark on computing distance functions}, volume = {163}, number = {1}, journal = {Journal of Computational Physics}, author = {Russo, Giovanni and Smereka, Peter}, year = {2000}, pages = {51--67} }

@article{min2007, author = {Chohong Min and Frédéric Gibou}, title = {A second order accurate level set method on non-graded adaptive cartesian grids}, journal = {Journal of Computational Physics}, year = {2007}, volume = {225}, pages = {300-321}, issn = {0021-9991}, doi = {10.1016/j.jcp.2006.11.034}, }

@article{min2010, author = {Chohong Min}, title = {On reinitializing level set functions}, journal = {Journal of Computational Physics}, year = {2010}, volume = {229}, pages = {2764-2772}, issn = {0021-9991}, doi = {10.1016/j.jcp.2009.12.032}, }

@hal{limare2022, hal-03889680} ~~~

Basilisk CFD

Navier–Stokes

Saint-Venant / Shallow Water

Multilayer

Compressible

Advection & VOF

Interfacial Forces

Two-Phase

Elliptic Solvers

Electrohydrodynamics

Viscoelasticity

Grid System

Core

Coordinates & Metrics

Input / Output

Utilities

Other

redistance.h

Redistancing of a distance field

Time derivative

The redistance() function

References