README

Organisation of the sandbox

Organisation of the sandbox

In this sandbox I present some work I’ve done using Basilisk during my PhD thesis (on mesh adaptation) and my post doc (on bubble dynamics).

Here is some work I’ve during my PhD on metric-based mesh adaptation.

If a message should be taken home, it would be: Every cursus learning computational methods should have an introduction on mesh and mesh adaptation, to enforce the idea that having a ‘good’ mesh is complicated, and that a carefully validated mesh adaptation method is a very very very interesting and powerfull tool. It is indeed too easy to manually create an inefficient mesh by hand (or, said differently: the element repartition of an optimal mesh (a mesh which minimizes the total error) is not obvious ! ).

A beautiful example which shows what is mesh adaptation is presented here

Basic theoretical principles

Let u be the exact solution of a PDE, u_h be the numerical solution of a discretized PDE and \Pi_h u the interpolation of the exact solution on the mesh. We try to reduce the total error ||u-u_h|| which contains local and non-local parts.

The exact solution interpolated on the mesh \Pi_h u doesn’t solves the discrete PDE. Thus, the numerical solution is modeled as u_h = \Pi_h u + u' with u' a correction. From that, the total error becomes ||u-\Pi_h u|| \leq ||u-\Pi_h u|| + ||u'|| where

||u-\Pi_h u|| is the interpolation error. This error is local (by definition). It is never null (except for linear equations): as soon as an exact solution solution is injected in a discretized mesh, the interpolation error can be defined.
||u'|| is a non necessarily local component of the total error.

We propose to reduce the total error by reducing the interpolation error. We observed the influence of ||u'|| and from that proposed a method which can restrict the minimal element size of a mesh. It is shown to reduce the total error when ||u'|| dominates it.

The interpolation error is derived from the metric-based theory and the the continuous mesh framework. It has been adapted to the grids of Basilisk.

Comparison with the previously existing AMR method in Basilisk

Basilisk already contains a wavelet-based AMR method. By construction, it only allows L^\infty-norm interpolation error estimation. However, L^2-norms errors are generally recommended when dealing with numerical solutions (sometimes other L^p-norms, with p<\infty).

Thus, we implemented a metric-based L^p-norm interpolation error. It is adapted from the anisotropic theory to the case of the isotropic elements of Basilisk quad/octree grids.

After that, we quantified the limits of interpolation error-based AMR criteria in the case of the 2nd order elliptic solver of Basilisk. It follows that:

interpolation error is not always the correct criterion to perform AMR. The total error of the numerical solution can propagate in the whole domain and this is linked with the cell size variation of the mesh.
by imposing a (non-trivial) minimal cell size, we slightly increase the interpolation error, but we can dramatically decrease the total numerical error, see this example.

Basic usage and example

The user simply need to add the new amr module (and all the required files available here)

       #include "./../../AMR_tools/amr.h"  // AMR

and to replace

      int maxl = 42;
      double eps=0.01;
      adapt_wavelet({psi}, (double[]){eps}, maxlevel=maxl);

      // restriction on the minimum grid size based on compression ratio (optional) 
      double etaopt = estimate_eta_opt(2, {psi});
      maxlevel = 0.5*log(Nobj/etaopt)/log(2.);   // maxlevel: global variable. Optional.

      // epsilon criteria for cell refinement/coarsening (mandatory) 
      AMReps = 0.01;                             // AMReps: global variable

      // AMR 
      adapt_metric( {psi} );                     // user interface similar to adapt_wavelet()

Note that the value for global variable AMReps can be eventually obtained in order to obtain iteratively a given number of elements using the function (this should apply to adapwavelts too)

AMReps = update_epsilon_control(Nobj);

The optimal converge rate for a given field can be simply obtained using the theoretical prefactor constants

struct PreFactorData cd = compute_prefactors (2, {psi}, NULL );

The minimum theoretical (interpolation) error only depends on the structure of the function and is \epsilon = c_d \; \Delta x^{-2}

double errtheory = cd.copt*pow(grid->tn,-(2./dimension)) )

A simple example to show that the user-interface is easy to use, and where we solve the laplace equation in axi-symmetric coordinates around a circular embedded boundary.

AMR algorithm files

The links to my AMR algorithm .h files.

AMR: performs AMR + user-interface

metric-based: tools to perform metric-based mesh adaptation and in particular compute metric-based interpolation errors.

linalg: tools for operation on tensors/matrices.

The last update was made in January 2024. It was necessary due to the struct Adapt being removed from src/ during the previous year, as the funtion calls using it were deemed to be done in an obsolete way. As a consequence, the current implementation of my code triggers a warning at compilation, but it has no impact on the computations.

AMR Examples

Boundary-Layer-Like Function

For this test case, the interpolation error alone provides nearly optimal meshes. Introducing the mesh-size constraint only slightly penalizes the interpolation error.

Exponential Function

For this function, the interpolation error alone does not lead to optimal meshes. The term ||u'|| dominates the adaptation process and becomes disconnected from the interpolation error when the mesh compression ratio is too large. The mesh-size constraint limits this compression ratio and enables the construction of nearly optimal meshes.

3D Exponential Function

Extension of the previous example to three dimensions.

Article and PhD thesis

AMR article: A metric-based adaptive mesh refinement criterion under constrain for solving elliptic problems on quad/octree grids. It is on HAL.

my PhD manuscript

Registered presentations

PMMH-\partial’Alembert interlab visit 2024: A 2 minutes popularizing science video presenting a numerical experiment of a cavitation bubble collapse. Both mesh adaptation and Arbitrary Lagrangian Eulerian methods are used. No equations are shown.

\partial’Alembertiennes 2023: A 2 minutes video about mesh adaptation and total numerical error minimization.

BGUM 2023

my PhD defense

Outputs

Usefull output functions.

vtk outputs: .vtk output for AMR grids.

Integral computation: compute the integral other the domain of either a scalar field or the difference between a scalar field and an analytical function using Gauss quadratures. In particular, I use these functions to compute the total error of numerical solutions when the reference analytical solution is known. NB: the current version of Basilisk triggers a warning at the compilation step due to the function calls I used being now deemed obsolete. However, it has no impact on the computation.

Miscellaneous

### 4th order Runge-Kutta scheme

4th order Runge-Kutta for a double using Basilisk events.

Rayleigh-Plesset equation time integration using the RK4 scheme.

Keller-Miksis equation time integration using the RK4 scheme.

Other

mathematica: C-functions as written by the C-export of mathematica

awk is wonderful <3: awk is so wonderful and so powerful that I want to share examples on how we can use it.