/**
# The Grid Adaptation Algorithm Based on a Wavelet-Estimated Discretization Error
Basilisk can employ anisotropic meshes using tree-based grids. If you
are not familiar with this grid structure, [this
link](http://basilisk.fr/sandbox/Antoonvh/The_Tree-Grid_Structure_in_Basilisk)
might be helpfull. The tree structure facilitates a convienient basis
for refinement and coarsening, such that the mesh resolution can vary
over the course of the simulation based on the evolution of the
solution itself (i.e. adaptive). This requires a decision
algorithm. Ideally such algorithm is generally applicable for a
multitude of solvers. Basilisk's wavelet-based strategy is designed to
be such an adaptation algorithm. This page aims to discribe how it can
be employed and what goes on behind the scenes. The algorithm is
designed and implemented by [Stéphane
Popinet](/sandbox/popinet/README). It is consisely described in
[Popinet (2015)](https://hal.archives-ouvertes.fr/hal-01163101v2)
This page was written to guide the writing for another publication, by
[Van Hooft et
al.(2018)](https://doi.org/10.1007/s10546-018-0335-9). You may find a
another version of the content that is presented on this page via the
link.
## Syntax
The function `adapt_wavelet` is defined in
[grid/tree-common.h](http://basilisk.fr/src/grid/tree-common.h#151). In
short, the adaptive wavelet algorithm is based arround the estimation
of numerical errors in the representation of spatially-discretized
fields. The `adapt_wavelet()` function requires the user to define a
list of fields that it will analyze for the
refinement/coarsening. This list can consist of any combination of the
existing scalar fields. In general, it makes sense to use the fields
that appear in the equations that are being solved. Additionally the
maximum tolerated estimated error for each field needs to be
defined. Finally, a maximum level (i.e. resolution) that the algorithm
is allowed to employ should be provided. Optionally, two more inputs
can be defined. i.e. 1: the minimum level of refinement and 2: A list
of scalars that need to be updated. The default values of these
optional arguments are 1 and `all`, respectively. Adaptation can be
done according to the following example event, where two dimensional
scalar fields `A, B, C` and vector field `u` are present:
~~~literatec
...
event adapt(i++) {
adapt_wavelet((scalar*){A, u},(double[]){0.1,0.2,0.3},8,3,{A, u, B});}
....
~~~
In this example, adaption is done every timestep and is based on the
three fields associated with `A` and `u` (i.e. `u.x` & `u.y`), The
maximum tolerance of the estimated error is 0.1, 0.2 and 0.3 for `A`,
`u.x` and `u.y`, respectively. The grid is allowed to vary between 3
and 8 levels of refinement. Scalar field C is not updated after
adaption, meaning that this field will not have sensible values
initialized at locations where the grid is either refined or
coarsened. This could be for example a diagnostic field that does not
appear in any computation or output routine before its values are
computed again.
## Wavelet based error Estimation
The mathematical rigor of the used algorithm is rooted in the field of
Multi-resolution analysis. Consider a one dimensional signal $f$ that
is discretized using an even number ($n$) elements, we can write
$f_n$. The $i$-th entry of the signal $f_n$ is denoted as $f^i_n$. We
now define a down sampling operator ($D$) that coarsens the original
signal to a lower level solution such that,
$$f_{n/2}=D\left(f_n\right).$$
Alternatively, we define an upsampling operator ($U$) such that,
$$g_{n}=U(f_{n/2}).$$
Note that in general, $g^i_n \neq f^i_n$ and that the difference
defined as,
$$ \chi^i_n = |f^i_n-g^i_n|,$$
can be interpreted as an error associated with the subsequent
application of the downsampling and upsampling operators to the signal
$f_n$. For the downsamping operation we choose local volume-avereaging
of the fine resolution solution to obtain a coarse resolution
solution. In (Basilisk) practice this means that for a N-dimensional
grid, $2^N$ leaf cells that share a common parent cell are averaged
and the resulting value is assigned to the corresponding parent
cell. Notice that this operation is exact for a finite volume
formulation. For the upscaling operator however, a second order
accurate interpolatiton seems suitable as Basilisk (typically) employs
second-order-accurate formulations for its solvers. This entails using
the just, bi or tri linear interpolation techniques for 1,2 and 3
dimensional fields, respectively, using the coarse level
solution. Given the implementation of this formulation one can
evaluate $\chi(i)$ and destinguish tree cases regaring the grid cell'
resolution by defining a threshold on the allowed error ($\zeta$):
$$\text{The\ }i\text{-th Grid Cell is\ }
\begin{cases}
&\text{Too Coarse.} & \chi^i_n > \zeta,\\
&\mathrm{Too\ Fine.} & \chi^i_n < \frac{2\zeta}{3},\\
&\mathrm{Just\ Fine.} & \mathrm{Otherwise},
\end{cases}$$
$\zeta$ represents the so-called refinement criterion that has the
same units as $f$. Notice that since the down and up sampling operator
are defined using only local data, the signal $f$ is not required to
be defined on an globally equidistant grid. However to ensure a local
$3^N$ coarse grid stencil for the linear interpolation, the resolution
between cells is only allowed to vary one level of refinement. The
figure below demonstrates how the discribed algorithm accesses the
resolution of a signal ($f_n$) according to a refinement criterium
$\zeta$.
![A visual representation of the various steps done by the algorithm
to assess the resolution of a given
grid.](http://www.basilisk.fr/sandbox/Antoonvh/errorestimation.jpg){width=100%}
### An example application
This section aims to exemplify how the adaption algorithm assesses a
discretized signal and adapts the grid according to a refinement
criterion ζ. For this purpose, we apply the algorithm to a subset of
the data from the simulation of forced isotropic turbulence in Li et
al. (2008). The simulation is run using a fixed equidistant grid with
$1024^3$ nodes; in terms of the Kolmogorov length scale (η), the grid
spacing ($\Delta_i$) is $\Delta_i$=2.2η. For the analysis we assume
the data to be resolved well enough, and the results are kindly made
available via the [Johns Hopkins turbulence
databases](http://turbulence.pha.jhu.edu/). We analyze a 2D slice of
the data (i.e. $1024^2$ cells) and for simplicity, we only consider
the velocity component perpendicular to the sliced plane ($u_⊥$). The
data are presented in Fig. 5a (see below); using the algorithm
described in the previous section, we can evaluate the χ field
corresponding to the original $u_⊥$ field. A section of the resulting
field, indicated by the black box in Fig. 5a, is shown in Fig. 5b,
where we can clearly see that the estimated discretization error is
not distributed uniformly by the equidistant-grid approach that was
used in the simulation. Rather, it appears that there are anisotropic
structures present, visualized by relatively high χ values (in
yellow). These structures appear to correspond to vortex filaments
that characterize the dissipative structures of high-Reynolds-number
turbulence (Frisch 1995). This result motivates the application of the
grid refinement algorithm to the data sample shown. Note that we
cannot add new information by refinement and at this point we do not
make any claims regarding what χ values are reasonable for a
turbulence-resolving simulation (this will depend on the
numerical). As such, we only allow the algorithm to coarsen the field
with a maximum error threshold ζ. The number of grid cells resulting
from the application of the adaptation algorithm for a range of ζ
values is shown in Fig. 5c; as expected, the number of grid cells
decreases with an increasing ζ value. Note that the plot also shows
that even for the high ζ values, the grid still contains cells at the
maximum resolution.
The main concept of employing the described grid-adaption algorithm is
visualized in Fig.5d. Here histograms of the number of grid cells
within 512 equally-spaced χ bins are presented for the original data
and the data obtained from applying the grid adaptation technique with
three different refinement criteria. It appears that for the original
dataset, the histogram is monotonically decreasing with increasing
χ. This shows that many grid cells exist where the numerical solution
is relatively smooth compared to cells in the tail of the
histogram. Hence, if the grid is chosen such that the discretization
errors in the latter region do not affect the relevant statistics of
the flow evolution, then the grid must be over-refined elsewhere. The
histograms of the adapted grids show that the algorithm is able to
lower the number of grid cells with low χ values, such that fewer grid
cells are employed. Note that the grid coarsening does not introduce
new grid cells with χ>2ζ/3, as this part of the histogram remains
unaltered.
When grid cells with a small but finite χ value are coarsened, some of
the data are lost and in general cannot be exactly reconstructed by
interpolation techniques. In order to assess how the data from the
adapted grids compare with the original data, Fig. 5e presents the
corresponding power spectra. It appears that none of the adapted grid
data are able to exactly reproduce the original power spectrum; more
specifically, with increasing ζ values, the wavenumbers ($k$) that
show a significant deviation in $E(k)$ from the original appear to
decrease. We point out that in order to evaluate the spectrum we have
linearly interpolated the data from the non-uniform grids to an
equidistant grid with 1024×1024 data points. The choice of the
interpolation technique is arbitrary and will pollute the diagnosed
spectrum in a non-trivial manner. As such, we directly compare all
10242 $u_⊥$(x,y) samples in Fig. 5f, where we see that the deviation
of the data from the 1 : 1 line is a function of ζ.
![Fig. 5 Example of the adaption algorithm applied to a 2D slice of a
3D turbulent field. a Shows the data slice of the velocity component
in the plane-perpendicular direction ($u_⊥$, obtained from Li et
al. (2008). b Presents the χ field, evaluated using the method
described in the previous section Only the centre part of the slice,
indicated by the black box in a, is shown to reveal the small-scale
details in this simulation. c shows the grid-cell-number dependence
on the chosen refinement criterion (ζ), note the logarithmic vertical
axis. A histogram of the χ field with 512 bins for the original data,
and the data corresponding to three ζ values are presented in
d. Using the same colour coding as in d, power spectra and a direct
comparison of the $u_⊥$(y,z) field are shown in e, f,
respectively](http://www.basilisk.fr/sandbox/Antoonvh/exampleadaptation.png)
###The down and upsampling in Basilisk
In Basilisk various different downsampling operations are used. In the
code this is referred to as the $restriction$ operation. It is defined
as an attribute of a given scalar field. By default basilisk uses the
following definition when initializing a new scalar field $s$ (see
[here](http://basilisk.fr/src/grid/multigrid-common.h#199)):
~~~literatec
s.restriction = restriction_average;
~~~
The upsampling operation is referred to as $prolongation$. By default,
for a scalar field, this is defined as (see
[here](http://basilisk.fr/src/grid/multigrid-common.h#198)):
~~~literatec
s.prolongation = refine_bilinear;
~~~
One may change these attributes if the argumentation for the chosen
restriction and prolongation formulation does not apply for your
solver formulation. This is for example the case for a volume of fluid
(fraction) field $c$ (see [here](http://basilisk.fr/src/vof.h)).
~~~literatec
c.prolongation = fraction_refine;
~~~
Many attributes for prolongation and restriction are defined in the
[common multigrid header
file](http://basilisk.fr/src/grid/multigrid-common.h), but you are
free to employ your own definitions,
e.g. [here](http://basilisk.fr/sandbox/Antoonvh/mandelbrot3.c).
###So where are the wavelets?
Similar to the beter known Fourier transformation, a wavelet based
transformation entails the decomposition of a signal into a set of
orthorgonal functions. The main difference is that the used orthogonal
functions ($\psi$) are localized in both physical space ($x$) as in
the frequency domain ($f$). A wavelet function ($\psi(x)$) can be
rescaled with a factor $a$ (dyadic dilation), and shifted in space by
$b$ (dyadic positioning) to obtain a set of wavelet functions,
$$\psi_{ab}(x)=\psi \left(\frac{x-b}{a}\right),$$
that may be orthogonal for a prober choice of $\psi(x),a$ and
$b$. Once such choice results in the
[Haar-wavelets](https://en.wikipedia.org/wiki/Haar_wavelet):
$$\psi_H(x)=\begin{cases}
1 \quad & 0 \leq t < \frac{1}{2},\\
-1 & \frac{1}{2} \leq t < 1,\\
0 &\mathrm{otherwise.}
\end{cases}
$$
Supplemented with a pair of integers n,k that determine the dyadic
scaling and positioning, respectively, one obains an orthogonal set of
wavelets,
$$\psi_{nk}=2^{n/2}\psi_H(2^nx-k)$$
The associated wavelet transform can be carried out
[discretely](https://en.wikipedia.org/wiki/Discrete_wavelet_transform)
using the so-called Haar-transform matrix $H$ for a (1D), discretized
signal $f_n$, $$W_n=Hf_n,$$
where $W_n$ is a vector of the coefficients for the various
wavelet-components. The idea is that not all components contribute
equally to the original signal, so one may choose to do without the
components that are smaller than some threshold $\zeta$,
$$T_i =
\begin{cases}
W_i, & |W_i| \geq \zeta,\\
0, & |W_i| < \zeta.
\end{cases}$$
Here $T_i$ is the results of the so-called wavelet thresholding
processes of $W_n$. We can obtain an approximation of the orginal
signal f_n, F_n, $$F_n = H^{-1}T_n.$$ Let us look the result of such a
transformation:
![Haar transform of a signal containing 256 points using $\zeta$=0.1](http://www.basilisk.fr/sandbox/Antoonvh/haar1.jpg)
We see that at the locations of low variablility (i.e. low derivative) the signal is approxmimated coarsly, and vise versa.
...SOME TEXT THAT MAKES A CONNECTION...
Let us look what the wavelet-based algorithm of basilisk produces for
the same signal, using the following code:
~~~literatec
#include "grid/bitree.h"
#include "utils.h"
scalar f[];
int main(){
FILE * fp1 = fopen("original","w");
FILE * fp2 = fopen("adapted","w");
f.prolongation=refine_injection;
X0=0;
L0=1;
init_grid(256);
foreach(){
f[]= sin(((x)*M_PI*5)-11.5)* exp(-(sq((x-0.5)/0.15)));
fprintf(fp1,"%g\t%g\n",x,f[]);
}
while(adapt_wavelet({f},(double[]){0.1},8).nc){
foreach()
f[]= sin(((x)*M_PI*5)-11.5)* exp(-(sq((x-0.5)/0.15)));
boundary({f});
}
double xp=L0/512;
while(xp