Adaptivity

From Medusa: Coordinate Free Mehless Method implementation
Revision as of 14:07, 11 June 2019 by Jureslak (talk | contribs) (Helmholz equation)

Jump to: navigation, search

Go back to Examples.

Solutions to many physical problems governed by partial differential equations (PDE) often significantly vary in magnitude throughout the problem domain. Although in some special cases the areas with high error are known in advance, in general the error distribution is unknown beforehand. Adaptive techniques for solving PDEs are a standard way of dealing with this problem, where problematic regions are iteratively refined. A step further is automatic adaptivity, where problematic regions are chosen automatically using an error indicator and then refined, until certain error threshold is reached. Below, we show some examples of fully automatic adaptivity in Medusa.


Basic concept

The adaptive methodology in this paper behaves similarly to "remeshing" used commonly in FEM. Some initial (possibly variable) nodal spacing h^0 is chosen, as well as its lower and upper bounds h_L and h_U, respectively. 3 Domain \Omega is filled with nodes, conforming to h^0 and the solution u^0 is obtained. An error indicator is employed to determine which nodes should be (de)refined and the nodal density h^0 is altered appropriately. This adaptive cycle below is repeated until the convergence criterion is met. The procedure on j-th iteration is written in more detail below:

  1. Fill \Omega with nodes conforming to h^j.
  2. Solve the problem to obtain u^j.
  3. Compute the error indicator values \varepsilon_i^j for each node p_i.
  4. If the mean of \varepsilon_i^j is below some tolerance \varepsilon return u^j as the solution and stop.
  5. Adapt h^j to obtain h^{j+1}.

More details can be found in our paper: https://arxiv.org/abs/1811.10368

Node density adaptation

The existing nodal spacing function h^j is evaluated at nodes p_i to obtain values h_{i,j} = h^j(p_i) . These values h_{i,j} are modified by a density factor f_i as h_{i,j+1} = \min(\max(h_{i,j} / f_i, h_L(p_i)), h_U(p_i)),

where density factor f_i is computed as f_i = \begin{cases} 1 + \frac{\eta - \varepsilon_i}{\eta - m} (\frac{1}{\beta} - 1), & \varepsilon_i \leq \eta, \quad \text{i.e. decrease the density} \\ 1, & \eta < \varepsilon_i < \varepsilon, \quad \text{i.e. no change in density}\\ 1 + \frac{\varepsilon_i - \varepsilon}{M - \varepsilon} (\alpha - 1), & \varepsilon_i \geq \varepsilon, \quad \text{i.e. increase the density} \end{cases}
and \alpha represents the refine aggressiveness, \beta the derefine aggressiveness, \varepsilon the refinement threshold, \eta the derenfinement threshold, and m = \min_i \varepsilon_i is the minimal and M = \max_i \varepsilon_i is the maximal value of the error indicator. Note that setting \alpha=1 or \beta=1 disables refinement and derefinement, respectively.


Error indicators

The work on error indicators is ongoing. For now, we use an ad hoc error indicator \varepsilon_i = \operatorname{std}_{j \in I_i}(u_j),

which represents the standard deviation of function values over all stencil nodes of a given node p_i.

Numerical exmaples

Below are several numerical examples where adaptivity has been tested or used to obtain solutions.

The errors e_1, e_2, e_\infty and e_E refer to relative discrete L^1, L^2, L_\infty and energy norm errors, respectively. These are evaluated in the computation nodes or on a denser grid by reinterpolation.

L shaped domain

The L shaped domain problem is defined on \Omega = [-1, 1]^2 \setminus [0, 1] \times [-1, 0]. The Laplace problem \nabla^2u = 0 with the solution u = r^{\frac{2}{3}} \sin(\frac{2}{3}\theta) given in polar coordinates.

BF-FD method with Polyharmonic splines augmented with monomials up to and including 2nd order was used to approximate the differential operators. The stencils for each node were chosen by simply selecting the closest n=15 nodes. The resulting sparse system was solved using the Intel ® MKL Pardiso sparse solver. Both uniform and fully adaptive refinement was tested. The adaptive procedure was run with \alpha=3, \varepsilon = 10^-2, \beta=1 and \eta=0.

The errors under uniform (left) and adaptive (right) refinement are shown below.

L shape uniform error.pngL shape adaptive error.png

The error (left) and the nodal density (right) during the adaptive iteration are shown below.

L shape progress.png

Disk under stress

Hertzian contact

Fretting fatigue contact

Bousinesq problem

Helmholz equation

Osc2d uniform error.pngOsc2d adaptive error.png

Osc2d adaptive sol.pngOsc2d adaptive den.png

Osc3d adaptive error.png

Osc3d adaptive sol.pngOsc3d adaptive error.png


Go back to Examples.