- Generated by
1.9.8
|
Basilisk CFD
Adaptive Cartesian mesh PDE framework
|
Go to the source code of this file.
| 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 | ||
| ) |
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 \(x_{i+1} - x_i\) 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 \(df_b\) (resp. \(df_t\)) is non-zero, then Neumann boundary conditions are applied i.e.
\[ \partial_n P(0) = df_b \]
where \(n\) is the direction "normal" to the boundary i.e. \(+x\) for the bottom/left boundary and \(-x\) for the top/right boundary.
If \(df_b\) (resp. \(df_t\)) is zero, then Robin/Navier boundary conditions are applied i.e.
\[ P(0) = f_b + \lambda_b \partial_n P(0) \]
To impose a Neumann zero boundary condition one can set \(\lambda_b\) to HUGE and \(f_b\) 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.
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.
We go through each new interval.
We first check if the starting point (x) of the integration interval is beyond the current interval [xpos[k]:xpos[k+1]].
If this is the case, we need to consider the next interval and compute the corresponding polynomial coefficients.
The top layer is a special case.
If limiting is used together with Neumann zero fluxes, the value must be constant in the last layer.
Otherwise we use Neumann or Robin boundary conditions to compute the polynomial.
Since the functions above are written for the bottom layer, we need to apply symmetry conditions i.e. transform \(x\) into \(1 - x\). This gives the following polynomial coefficients.
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.
We use boundary conditions for the bottom layer.
Or central remapping, with optional limiting, using the left and right values for all the other layers.
The new left value is just the old right value.
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.
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.
The new integration interval starts at xe and concerns jnew.
Definition at line 340 of file remapc.h.
References a, assert, b, dx, HUGE, i, k, lambda_b, remap_central(), remap_neumann(), remap_robin(), right_value(), v, and x.
Referenced by vertical_remapping().
|
static |
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 \([x_k:x_{k+1}]\).
If limiting is used, we first limit the left and right values.
The polynomial fit is given by the linear system
\[ \begin{aligned} P(0) &= s_l, \\ \int_0^1P(x)dx &= f_k \\ P(1) &= s_r \end{aligned} \]
which has for solution
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 \(C_2\neq0\)) is at \(x=-\frac{C_1}{2C_2}\) which can be developed as
\[ x=-\frac{6f_k-2s_r-4s_l}{6(s_r+s_l-3f_k)} \]
If \(x\in[0,0.5]\), we change \(s_r\) such that \(x=0\). The new value is then \(s_r=3f_k-2s_l\).
If \(x\in[0.5,1]\), we change \(s_l\) such that \(x=1\). The new value is then \(s_l=3f_k-2s_r\).
We update the polynomial coefficients using these new bounds.
Definition at line 221 of file remapc.h.
References dx, f, k, minmodremap(), sigma, v, and x.
Referenced by remap_c().
This is inspired by the PPR Fortran90 library 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 within the multilayer framework.
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.
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
\[ P(x) = \sum_i C_i x^i \]
It must verify the following conditions
\[ \begin{aligned} P'(0) &= df_b, \\ \int_0^1P(x)dx &= f \\ P(1) &= s_r \end{aligned} \]
which leads to the code below.
Definition at line 45 of file remapc.h.
Referenced by remap_c().
This is a bit more complicated, with the following conditions
\[ \begin{aligned} P(0) &= f_b + \lambda_b P'(0), \\ \int_0^1P(x)dx &= f \\ P(1) &= s_r \end{aligned} \]
which leads to the solution of the linear system coded below.
Definition at line 68 of file remapc.h.
References a, f, i, j, lambda_b, s, and x.
Referenced by remap_c().
|
static |
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
\[ \int_{x_i}^{x_{i+1}} P(x') dx' = f_i \]
with \(k - 1 \leq i \leq k + 2\). The integration is normalized using
\[ x' = \frac{x}{x_{k+1} - x_k} \]
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 df_b is non-zero we apply a Neumann boundary condition, which changes one equation of the default linear system to
Otherwise we apply a Robin condition, as
This is the equivalent for the top layer. The systems are different at this location because \(x'=1\) for the top layer and \(x'=0\) for the bottom layer.
We invert the linear system to get the coefficients \(C_i\) of the polynomial and return the value on the "right" side of the interval by computing \(P(x'=1)\).
Definition at line 103 of file remapc.h.
References assert, c, dx, f, i, j, k, lambda_b, n, smatrix_inverse(), sq(), v, and x.
Referenced by remap_c().
1.9.8