remapc.h

Go to the documentation of this file.

/** @file remapc.h
 */
/**
# Piecewise Polynomial Reconstruction
 
This is inspired by the [PPR Fortran90
library](https://github.com/dengwirda/PPR) of Darren Engwirda and
Maxwell Kelley, which is included in [/src/ppr]().
 
The initial version was written by Antoine Aubert.
 
It is used for [layer remapping](/src/layered/remap.h) within the
[multilayer framework](/src/layered/README).
 
The basic idea is to construct a second-order polynomial (a parabola)
fitting the integral over each layer and layer-interface values
obtained with fourth-order polynomial fits. Monotonic limiting can
also be applied.
 
A significant difference with the PPR library is that explicit boundary
conditions (Neumann or Robin/Navier) are imposed on the left and right
(or top and bottom) boundaries. Also, being written in C99, this code
runs transparently on GPUs.
 
## Neumann boundary conditions
 
We start with the simplest case i.e. the reconstruction of a parabola
with Neumann conditions on the left boundary and a Dirichlet condition
on the right boundary.
 
The polynomial is given by
\f[
P(x) = \sum_i C_i x^i
\f]
It must verify the following conditions
\f[
\begin{aligned}
P'(0) &= df_b, \\
\int_0^1P(x)dx &= f \\
P(1) &= s_r
\end{aligned}
\f]
which leads to the code below. */
 
static void remap_neumann (double s_r,
                           double f,
                           double df_b,
                           double C[3])
{
  C[1] = df_b;
  C[2] = 3.*(s_r - df_b/2. - f)/2.;
  C[0] = s_r - df_b - C[2];
}
 
/**
## Robin/Navier boundary conditions
 
This is a bit more complicated, with the following conditions
\f[
\begin{aligned}
P(0) &= f_b + \lambda_b P'(0), \\
\int_0^1P(x)dx &= f \\
P(1) &= s_r
\end{aligned}
\f]
which leads to the solution of the linear system coded below. */
 
static void remap_robin (double s_r,
                         double f,
                         double f_b, double lambda_b,
                         double C[3])
{
  double a = - lambda_b, s = - a + (a - 1.)/3. + 1/2. [0];
  double M[3][3] = {
    {    - a,    1./6.,     a/3.    },
    {     1.,  - 2./3.,  - 1./3.[0] },
    { a - 1.,    1./2,  1/2. - a    }
  };
  double B[3] = { f, f_b, s_r };
  for (int i = 0; i < 3; i++) {
    C[i] = 0.;
    for (int j = 0; j < 3; j++)
      C[i] += M[i][j]*B[j];
    C[i] /= s;
  }
}
 
/**
## Fourth-order polynomial for interface values
 
The values between layers are approximated by fitting a fourth-order
polynomial to the two layers above and below the interface. This gives
the 4x4 linear system
\f[
\int_{x_i}^{x_{i+1}} P(x') dx' = f_i
\f]
with \f$k - 1 \leq i \leq k + 2\f$. The integration is normalized using
\f[
x' = \frac{x}{x_{k+1} - x_k}
\f]
*/
 
static void right_value (int nvar, double s_right[nvar],
                         int k, int n,
                         const double x[n], double f[n-1][nvar],
                         double f_b, double lambda_b, double df_b,
                         double f_t, double lambda_t, double df_t)
{
  double xk = x[k], dx = x[k+1] - xk;
  double M[4][4];
  for (int i = (k == 0); i < (k == n - 3 ? 3 : 4); i++) {
    double zeta1 = (x[i+k-1] - xk)/dx;
    double z1n = zeta1;
    double zeta0 = (x[i+k] - xk)/dx;
    double z0n = zeta0;
    for (int j = 0; j < 4; j++, z1n *= zeta1, z0n *= zeta0)
      M[i][j] = (z0n - z1n)/(j + 1);
  }
 
  /**
  The left (or bottom) and right (or top) layers must use the left and
  right boundary conditions to close the system. 
 
  For the bottom layer we have: */
  
  if (k == 0) {
 
    /**
    If `df_b` is non-zero we apply a Neumann boundary condition, which
    changes one equation of the default linear system to */
    
    if (df_b) {
      M[0][0] = 0.;
      M[0][1] = 1.;
      M[0][2] = 0.;
      M[0][3] = 0.;
    }
 
    /**
    Otherwise we apply a Robin condition, as */
    
    else {
      M[0][0] = 1. [0];
      M[0][1] = - lambda_b/dx;
      M[0][2] = 0.;
      M[0][3] = 0.;
    }
  }
 
  /**
  This is the equivalent for the top layer. The systems are different
  at this location because \f$x'=1\f$ for the top layer and \f$x'=0\f$ for the
  bottom layer. */
  
  if (k == n - 3) {
    double zetab = (x[k+2] - xk)/dx;
    if (df_t) {
      M[3][0] = 0.;
      M[3][1] = 1.;
      M[3][2] = 2.*zetab;
      M[3][3] = 3.*sq(zetab);
    }
    else {
      M[3][0] = 1.;
      M[3][1] = zetab - lambda_t/dx;
      M[3][2] = zetab*(zetab - 2.*lambda_t/dx);
      M[3][3] = sq(zetab)*(zetab - 3.*lambda_t/dx);
    }
  }
 
  /**
  ### Solution of the linear system
 
  We invert the linear system to get the coefficients \f$C_i\f$ of the
  polynomial and return the value on the "right" side of the interval
  by computing \f$P(x'=1)\f$. */
  
  assert (smatrix_inverse (4, M, 1.e-30) > 1.e-20);
 
  for (int v = 0; v < nvar; v++) {
    double B[4];
    for (int i = (k == 0); i < (k == n - 3 ? 3 : 4); i++)
      B[i] = (x[k+i] - x[k-1+i])/dx*f[k-1+i][v];
    if (k == 0) {
      if (df_b)
        B[0] = df_b*dx;
      else
        B[0] = f_b;
    }
    if (k == n - 3) {
      if (df_t)
        B[3] = df_t*dx;
      else
        B[3] = f_t;
    }
    double c[4] = {0.,0.,0.,0.};
    for (int i = 0; i < 4; i++)
      for (int j = 0; j < 4; j++)
        c[i] += M[i][j]*B[j];
    s_right[v] = c[3];
    for (int i = 0; i < 3; i++)
      s_right[v] += c[i];
  }
}
 
/**
## Polynomial reconstruction for "central" layers
 
We first define a simple monotonic limiter. */
 
static double minmodremap (double a, double b) {
  return a*b <= 0. ? 0. :
    fabs(a) < fabs(b) ? a : b;
}
 
/**
The reconstruction function take the left and right values, the layers
positions `x` and corresponding average values `f`, a limiting option
and returns the polynomial `C` on the interval \f$[x_k:x_{k+1}]\f$. */
 
static void remap_central (int nvar,
                           const double s_l[nvar], const double s_r[nvar],
                           const int n, const double x[n], double f[n-1][nvar],
                           const int k,
                           double C[3][nvar],
                           const bool limiter)
{
  double xk = x[k], dx = x[k+1] - xk;
  for (int v = 0; v < nvar; v++) {
    double sl = s_l[v], sr = s_r[v];
 
    /**
    If limiting is used, we first limit the left and right values. */
  
    if (limiter) {
      if ((f[k+1][v] - f[k][v])*(f[k][v] - f[k-1][v]) < 0.) {
        sl = f[k][v];
        sr = f[k][v];
      }
      else {
        double sigma_left = 2.*(f[k][v] - f[k-1][v])/dx;
        double sigma_center = 2.*(f[k+1][v] - f[k-1][v])/(x[k+2] - x[k-1] + dx);
        double sigma_right = 2.*(f[k+1][v] - f[k][v])/dx;
        double sigma = minmodremap(sigma_center, minmodremap(sigma_left, sigma_right));
        if ((sl - f[k-1][v])*(f[k][v] - sl) < 0.)
          sl = f[k][v] - 1./2.*dx*sigma;
        if ((sr - f[k][v])*(f[k+1][v] - sr) < 0.)
          sr = f[k][v] + 1./2.*dx*sigma;
      }
    }
 
    /**
    The polynomial fit is given by the linear system
    \f[
    \begin{aligned}
    P(0) &= s_l, \\
    \int_0^1P(x)dx &= f_k \\
    P(1) &= s_r
    \end{aligned}
    \f]
    which has for solution */
  
    C[0][v] = sl;
    C[1][v] = 6.*f[k][v] - 2.*sr - 4.*sl;
    C[2][v] = 3.*(sr + sl - 2.*f[k][v]);
 
    /**
    If limiting is used, we check that the extrema of the polynomial do
    not over- or undershoot. For a second-order polynomial, the extrema (if
    \f$C_2\neq0\f$) is at
    \f$x=-\frac{C_1}{2C_2}\f$ which can be developed as
    \f[
    x=-\frac{6f_k-2s_r-4s_l}{6(s_r+s_l-3f_k)}
    \f]  
    If \f$x\in[0,0.5]\f$, we change \f$s_r\f$ such that
    \f$x=0\f$. The new value is then \f$s_r=3f_k-2s_l\f$.
    
    If \f$x\in[0.5,1]\f$, we change \f$s_l\f$ such that
    \f$x=1\f$. The new value is then \f$s_l=3f_k-2s_r\f$.
    */
  
    if (limiter && C[1][v]*C[2][v] < 0. && C[1][v]/C[2][v] > - 2.) {
      if (C[1][v]/C[2][v] > - 1.)
        sr = 3.*f[k][v] - 2.*sl;
      else
        sl = 3.*f[k][v] - 2.*sr;
 
      /**
      We update the polynomial coefficients using these new bounds. */
    
      C[0][v] = sl;
      C[1][v] = 6.*f[k][v] - 2.*sr - 4.*sl;
      C[2][v] = 3.*(sr + sl - 2.*f[k][v]);
    }
  }
}
 
/**
## Remapping function
 
This is the remapping interface. It takes as arguments the number of
initial positions `npos`, the number of "remapped" positions `nnew`
and the corresponding arrays of initial and remapped positions `xpos`
and `xnew` respectively, as well as the initial averaged values on the
corresponding intervals `fdat`. The remapped averaged values will be
stored in `fnew`.
 
The `fdat` and `fnew` arrays are also indexed by `nvar` which allows
to remap several fields simultaneously, for performance. Note however
that the interface does not allow to apply different boundary
conditions for each of these fields.
 
Note that the `xpos` and `xnew` must be stored in strict increasing
order i.e. negative interval values \f$x_{i+1} - x_i\f$ are not allowed.
 
The "bottom" and "top" boundary conditions are given by the
$(f_b,\lambda_b,df_b)$ parameters (resp. $(f_t,\lambda_t,df_t)$). If
\f$df_b\f$ (resp. \f$df_t\f$) is non-zero, then Neumann boundary conditions
are applied i.e.
\f[
\partial_n P(0) = df_b
\f]
where \f$n\f$ is the direction "normal" to the boundary i.e. \f$+x\f$ for the
bottom/left boundary and \f$-x\f$ for the top/right boundary.
 
If \f$df_b\f$ (resp. \f$df_t\f$) is zero, then Robin/Navier boundary conditions
are applied i.e.
\f[
P(0) = f_b + \lambda_b \partial_n P(0)
\f]
 
To impose a Neumann zero boundary condition one can set \f$\lambda_b\f$ to
`HUGE` and \f$f_b\f$ to zero.
 
If limiter is `true` then monotonic limiters are applied. In this case
the Neumann zero condition leads to a constant profile in the top or
bottom layer, which guarantees boundedness of the remapping. Note that
other boundary conditions will not guarantee boundedness. */
 
void remap_c (int npos, int nnew,
              const double xpos[npos], const double xnew[nnew],
              const int nvar,
              double fdat[npos-1][nvar], double fnew[nnew-1][nvar],
              double f_b, double lambda_b, double df_b,
              double f_t, double lambda_t, double df_t,
              bool limiter)
{
 
  /**
  We can only remap on the same interval. It would be possible to
  remap on a smaller interval (a subset) but there is no need for now
  and this restriction makes things simpler. */
    
  assert (xnew[0] == xpos[0] && xnew[nnew-1] == xpos[npos - 1]);
 
  /**
  We go through each new interval. */
  
  double x = xnew[0], s_left[nvar], C[3][nvar];
  for (int v = 0; v < nvar; v++)
    fnew[0][v] = 0., s_left[v] = 0.;
  for (int inew = 0, k = -1; inew < nnew - 1 && x < xpos[npos - 1];) {
 
    /**
    ## Integration interval
 
    We first check if the starting point (x) of the integration
    interval is beyond the current interval [xpos[k]:xpos[k+1]]. */
    
    if (x >= xpos[k + 1]) {
 
      /**
      If this is the case, we need to consider the next interval and
      compute the corresponding polynomial coefficients. */
      
      k++;
 
      /**
      The top layer is a special case. */
      
      if (k == npos - 2) {
 
        /**
        If limiting is used together with Neumann zero fluxes, the
        value must be constant in the last layer. */
        
        if (limiter && lambda_t == HUGE)
          for (int v = 0; v < nvar; v++) {
            C[0][v] = fdat[k][v];
            C[1][v] = 0.;
            C[2][v] = 0.;
          }
 
        /**
        Otherwise we use Neumann or Robin
        boundary conditions to compute the polynomial. */
        
        else
          for (int v = 0; v < nvar; v++) {
            double B[3];
            if (df_t)
              remap_neumann (s_left[v], fdat[k][v], - df_t*(xpos[k+1] - xpos[k]), B);
            else
              remap_robin (s_left[v], fdat[k][v], f_t, - lambda_t/(xpos[k+1] - xpos[k]), B);
 
            /**
            Since the functions above are written for the bottom layer, we
            need to apply symmetry conditions i.e. transform \f$x\f$ into \f$1 -
            x\f$. This gives the following polynomial coefficients. */
            
            C[0][v] = B[0] + B[1] + B[2];
            C[1][v] = - B[1] - 2.*B[2];
            C[2][v] = B[2];
          }
      }
      else {
 
        /**
        For all other layers, we need to compute the right value. If
        limiting is used together with Neumann zero fluxes we impose a
        constant profile in the bottom or top layer. */
 
        double s_right[nvar];
        if (limiter && lambda_b == HUGE && k == 0)
          for (int v = 0; v < nvar; v++)
            s_right[v] = fdat[0][v];
        else if (limiter && lambda_t == HUGE && k == npos - 3)
          for (int v = 0; v < nvar; v++)
            s_right[v] = fdat[npos - 2][v];
        else
          right_value (nvar, s_right, k, npos, xpos, fdat, f_b, lambda_b, df_b, f_t, lambda_t, df_t);
        
        /**
        We use boundary conditions for the bottom layer. */
 
        if (k == 0)
          for (int v = 0; v < nvar; v++) {
            double B[3];
            if (df_b)
              remap_neumann (s_right[v], fdat[0][v], df_b*(xpos[1] - xpos[0]), B);
            else
              remap_robin (s_right[v], fdat[0][v], f_b, lambda_b/(xpos[1] - xpos[0]), B);
            for (int i = 0; i < 3; i++)
              C[i][v] = B[i];
          }
 
        /**
        Or central remapping, with optional limiting, using the left
        and right values for all the other layers. */
        
        else
          remap_central (nvar, s_left, s_right, npos, xpos, fdat, k, C, limiter);
        
        /**
        The new left value is just the old right value. */
 
        for (int v = 0; v < nvar; v++)
          s_left[v] = s_right[v];
      }
    }
 
    /**
    ## End of integration interval
 
    We determine the value of `xe`, the end of the integration
    interval, by comparing the end of the new interval, `xnew[inew+1]`
    with the end of the current interval `xpos[k+1]`. `jnew` records
    whether we need to change interval after the integration. */
    
    double xe;
    int jnew = inew;
    if (xnew[inew + 1] < xpos[k + 1]) {
      jnew = inew + 1;
      xe = xnew[jnew];
      for (int v = 0; v < nvar; v++)
        fnew[jnew][v] = 0.;
    }
    else
      xe = xpos[k + 1];
 
    /**
    ## Integration
 
    We compute the integral of the current polynomial between `x` and
    `xe`. Since the polynomial is defined on the normalized interval
    [0:1], we first normalize the values of `x` and `xe`, in `a` and
    `b` respectively. We also need to de-normalize the integral with
    the ratio `dx1`. Note that it is particularly important to perform
    the integration in normalized space to minimize round-off errors
    which can cause instabilities when using single-precision on
    GPUs. */
 
    double xk = xpos[k], dx = xpos[k+1] - xk, dx1 = dx/(xnew[inew + 1] - xnew[inew]);
    double a = (x - xk)/dx, b = (xe - xk)/dx;
    assert (a >= 0. && a <= b && b <= 1.); // xpos and xnew must be ordered properly
    for (int v = 0; v < nvar; v++) {
      double an = a, bn = b;
      for (int i = 0; i < 3; i++, an *= a, bn *= b)
        fnew[inew][v] += C[i][v]/(i + 1)*(bn - an)*dx1;
    }
 
    /**
    ## Update of integration bounds
 
    The new integration interval starts at `xe` and concerns `jnew`. */
 
    x = xe;
    inew = jnew;
  }
}