remap.h
📄 View in API Reference (Doxygen) · View on basilisk.fr
Test cases (13): bar, beach, conical, dispersion, dispersion_gravitocap, galilean_invariance, gaussian, kh, large, lock, overflow, stokes, wind-driven
Examples (4): breaking, held-suarez, lee, shoal
Vertical remapping
This implements a simple vertical remapping to "$\sigma$-coordinates" (equally-distributed by default).
We can optionally use the PPR Library of Engwirda and Kelley to perform the remapping. The default settings are using the Parabolic Piecewise Method without limiting.
By default we use the C PPR implementation without limiting.
#if USE_PPR
# include "ppr/ppr.h"
// int edge_meth = p1e_method, cell_meth = plm_method, cell_lim = null_limit;
int edge_meth = p3e_method, cell_meth = ppm_method;
// int edge_meth = p5e_method, cell_meth = pqm_method, cell_lim = null_limit;
#else // !USE_PPR
# include "remapc.h"
const bool null_limit = false, mono_limit = true;
#endif
int cell_lim = null_limit;The distribution of layers can be controlled using the *beta* array which defines the ratio of the thickness of each layer to the total depth $H$ (i.e. the relative thickness). By default all layers have the same relative thickness.
double * beta = NULL;
event defaults (i = 0)
{
beta = malloc (nl*sizeof(double));
for (int l = 0; l < nl; l++)
beta[l] = 1./nl;
}The default uniform layer distribution can be replaced with a geometric progression for the layer thicknesses. This needs to be called for example in the init() event. The rmin parameter specifies the minimum layer thickness relative to the uniform layer thickness (proportional to 1/nl). If the top parameter is set to true the minimum layer thickness is at the top (layer nl - 1), otherwise it is at the bottom (layer 0).
void geometric_beta (double rmin, bool top)
{
if (rmin <= 0. || rmin >= 1. || nl < 2)
return;
double r = 1. + 2.*(1./rmin - 1.)/(nl - 1.);
double hmin = (r - 1.)/(pow(r, nl) - 1.);
for (int l = 0; l < nl; l++)
beta[l] = hmin*pow(r, top ? nl - 1 - l : l);
}The *vertical_remapping()* function takes a (block) field of layer thicknesses and the corresponding list of tracer fields and performs the remapping (defined by *beta*).
trace
void vertical_remapping (scalar h, scalar * tracers)
{
const int npos = nl + 1, nvar = list_len(tracers);
#if USE_PPR
int ndof = 1;
#endif
foreach() {
#if HALF
double H0 = 0., H1 = 0., H;
foreach_layer() {
if (point.l < nl/2)
H0 += h[];
else
H1 += h[];
}
H = H0 + H1;
#else
double H = 0.;
foreach_layer()
H += h[];
#endif
if (H > dry) {
double zpos[npos], znew[npos];
#if USE_PPR
double fdat[nvar*nl], fnew[nvar*nl];
zpos[0] = znew[0] = 0.;
foreach_layer() {
zpos[point.l+1] = zpos[point.l] + max(h[],dry);
int i = nvar*point.l;
for (scalar s in tracers) {
dimensional (fnew[i] = s[]);
fdat[i++] = s[];
}
#if HALF
if (point.l < nl/2)
h[] = 2.*H0*beta[point.l];
else
h[] = 2.*H1*beta[point.l];
#else
h[] = H*beta[point.l];
#endif
znew[point.l+1] = znew[point.l] + h[];
}
my_remap (&npos, &npos, &nvar, &ndof, zpos, znew, fdat, fnew,
&edge_meth, &cell_meth, &cell_lim);
foreach_layer() {
int i = nvar*point.l;
for (scalar s in tracers)
s[] = fnew[i++];
}
#else // !USE_PPR
zpos[0] = znew[0] = 0.;
foreach_layer() {
zpos[point.l+1] = zpos[point.l] + max(h[],dry);
#if HALF
if (point.l < nl/2)
h[] = 2.*H0*beta[point.l];
else
h[] = 2.*H1*beta[point.l];
#else
h[] = H*beta[point.l];
#endif
znew[point.l+1] = znew[point.l] + h[];
}
znew[npos -1] = zpos[npos - 1];
double fdat[nl][nvar], fnew[nl][nvar];
int v = 0;
for (scalar s in tracers) {
foreach_layer() {
dimensional (fnew[point.l][v] = s[]);
fdat[point.l][v] = s[];
}
v++;
}For the moment we simply use Neumann zero top and bottom boundary conditions for all fields. Note that we also ignore the corresponding dimensions since this could be inconsistent when remapping quantities with different dimensions.
remap_c (npos, npos, zpos, znew, nvar, fdat, fnew,
0[*], HUGE[*], 0[*],
0[*], HUGE[*], 0[*],
cell_lim);
v = 0;
for (scalar s in tracers) {
foreach_layer()
s[] = fnew[point.l][v];
v++;
}
#endif // !USE_PPR
}
}
}The remapping is applied at every timestep.
event remap (i++) {
if (nl > 1)
vertical_remapping (h, tracers);
}The *beta* array is freed at the end of the run.
event cleanup (i = end)
{
free (beta), beta = NULL;
}