Commit 7aa80c51 authored by Praetorius, Simon's avatar Praetorius, Simon

Initial commit

parents
Pipeline #1262 failed with stages
in 9 seconds
*.pdf filter=lfs diff=lfs merge=lfs -text
*.png filter=lfs diff=lfs merge=lfs -text
*.jpg filter=lfs diff=lfs merge=lfs -text
[submodule "MathJax"]
path = MathJax
url = https://github.com/mathjax/MathJax.git
Subproject commit 419b0a6eee7eefc0f85e47f7d4f8227ec28b8e57
<!DOCTYPE html>
<html>
<head>
<title>AMDiS - Adaptive Multi-Dimensional Simulations</title>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"/>
<style type="text/css">
@import url(https://fonts.googleapis.com/css?family=Yanone+Kaffeesatz);
@import url(https://fonts.googleapis.com/css?family=Raleway);
@import url(https://fonts.googleapis.com/css?family=Ubuntu);
@import url(https://fonts.googleapis.com/css?family=Droid+Serif:400,700,400italic);
@import url(https://fonts.googleapis.com/css?family=Ubuntu+Mono:400,700,400italic);
</style>
<!--<link rel="stylesheet" type="text/css" href="style_display.css" />-->
<link rel="stylesheet" type="text/css" href="style_print.css" />
</head>
<body>
<textarea id="source">
class: center, middle
# Session 3
## Adaptivity and systems of equation
---
# Agenda
### Monday
- Scalar linear second order PDEs
- Handling data on unstructured grids
- **Adaptivity and systems of equations**
- Introduction to Student Projects
---
# Motivation
### Mesh adaptivity
- Two-phase flow problem
- Steep transition between two phases
<img src="images/capillary_ellipse.png" width="40%" /><img src="images/capillary_ellipse2.png" width="50%" />
---
# Adaptivity and error estimators: (`poisson3.cc`)
<table width="100%"><tr><td width="60%">
Adative solution strategy:<ul>
<li>Refine elements to reduce global error</li>
<li>Local refinement instead of global refinement</li>
<li>Repeated solution with adapted mesh:</li>
</ul></td><td><img src="images/bisection.png" width="100%" /></td>
</tr></table>
--
.center[
<img src="images/adaption-loop.png" width="90%" />
]
---
# Adaptivity and error estimators: (`poisson3.cc`)
```
prob.assemble(&amp;adaptInfo);
prob.solve(&amp;adaptInfo);
prob.estimate(&amp;adaptInfo); // local error indicator
Flag adapted;
for (int i = 0; i < MAX\_ITER; ++i)
{
if (adaptInfo.getEstSum(0) < TOL)
break;
Flag markFlag = prob.markElements(&amp;adaptInfo); // mesh adaption}
adapted = prob.refineMesh(&amp;adaptInfo);
adapted|= prob.coarsenMesh(&amp;adaptInfo);
prob.buildAfterCoarsen(&amp;adaptInfo, markFlag); // assemble}
prob.solve(&amp;adaptInfo);
prob.estimate(&amp;adaptInfo);
}
```
---
class: center, middle
# Better: `AdaptStationary`
--
.left[
### Encapsulation into class:
```
AdaptStationary adaptStat("adapt", prob, adaptInfo);
adaptStat.adapt(); // start adaptive solution process
```
]
---
# Adaptivity and error estimators: (`poisson3.cc`)
Solution strategy: `assemble->solve->estimate->mark->refine/coarsen`
- Estimator type can be set in init-file
--
- e.g. Residual-Estimator: Element error-indicator `\(\eta_T\)`:
\\[
\eta\_T := C\_0 h\_T^2 \||f + \nabla\cdot\mathbb{A}\nabla u\_h\||\_{L\_2(T)} + C\_1 \left(\sum\_{e\subset T} h\_e \left[\left[\mathbb{A}\nabla u\_h\cdot\underline{\nu\_e}\right]\right]^2\right)^{1/2}
\\]
--
- e.g. Recovery Estimator:
\\[
\eta\_T := \||\nabla u\_h^\ast - \nabla u\_h\||\_{L\_2(T)}
\\]
--
- Various *marking strategies*: global refinement, maximum strategy, equidistribution strategy, guaranteed error reduction strategy
---
# Stationary equation with adaptivity
Source-code: see `poisson4.cc`
```
ProblemStat prob("poisson5");
prob.initialize(INIT_ALL);
Operator opL( prob.getFeSpace(), prob.getFeSpace() );
addSOT( opL, A ); // e.g. A = 1.0
Operator opF( prob.getFeSpace() );
addZOT( opF, f ); // e.g. f = X()*X() + 1
prob.addMatrixOperator( opL, 0, 0 );
prob.addVectorOperator( opF, 0 );
prob.addDirichletBC( nr, 0, 0, new G );
AdaptInfo adaptInfo("adapt");
AdaptStationary("adapt", prob, adaptInfo).adapt();
prob.writeFiles(&amp;adaptInfo, true);
```
---
# Stationary equation with adaptivity
Init-file: see `poisson4.dat.2d`
```matlab
poisson->estimator[0]: residual % or recovery
poisson->estimator[0]->C0: 0.1 % weight of element residual
poisson->estimator[0]->C1: 0.1 % weight of jump residual
poisson->marker[0]->strategy: 2
% 1... global refinement
% 2... maximum strategy
% 3... equidistribution strategy
% 4... guaranteed error reduction strategy
adapt[0]->tolerance: 1e-4
adapt->max iteration: 5
```
---
# Stationary equation with adaptivity
Visualization of the locally refined grid: TODO: change images!!
<img src="images/local_refinement_0.png" width="30%" /><img src="images/local_refinement_3.png" width="30%" /><img src="images/local_refinement_4.png" width="30%" />
<img src="images/local_refinement_5.png" width="30%" /><img src="images/local_refinement_6.png" width="30%" /><img src="images/local_refinement_10.png" width="30%" />
---
class: center, middle
# System of equations
---
# System of equations
Motivation:
- Navier-Stokes equations
- Phase-Field (Crystal) equation
- Willmore flow
<img src="images/120particle_f3.png" width="30%" />
<img src="images/rho_big.png" width="30%" />
<img src="images/cell.png" width="30%" />
---
# Example: Biharmonic equation
We consider the stream-function formulation for the Stokes equation (2D):
\\[
-\nu\Delta^2\psi = f,\quad\text{ in }\Omega,\qquad\big(\mathbf{v} = -\nabla\times\mathbf{e}\_z \psi\big)
\\]
with `\(\psi|_{\partial\Omega}=\Delta\psi|_{\partial\Omega} = 0\)` and `\(f=(\nabla\times\mathbf{f})_z\)`.
--
Reformulated as system of second order equation, we obtain:
\\[
-\nu\Delta\phi = f,\quad\text{ in }\Omega \\\\
-\Delta\psi + \phi = 0
\\]
--
or in weak form:
\\[
\langle\nu\nabla\phi,\nabla\theta\rangle -\langle\nu\partial\_n\phi,\theta\rangle\_{\partial\Omega} = \langle f, \theta\rangle,\quad\forall\theta\in V \\\\
\langle\nabla\psi,\nabla\theta'\rangle -\langle\partial\_n\psi,\theta'\rangle\_{\partial\Omega} + \langle\phi,\theta'\rangle = 0,\qquad\forall\theta'\in V'
\\]
with solution components `\((\psi\in V, \phi\in V')\)`.
---
# Example: Biharmonic equation
We consider the stream-function formulation for the Stokes equation (2D):
\\[
-\nu\Delta^2\psi = f,\quad\text{ in }\Omega,\qquad\big(\mathbf{v} = -\nabla\times\mathbf{e}\_z \psi\big)
\\]
with **`\(\psi|_{\partial\Omega}=\Delta\psi|_{\partial\Omega} = 0\)`** and `\(f=(\nabla\times\mathbf{f})_z\)`.
Reformulated as system of second order equation, we obtain:
\\[
-\nu\Delta\phi = f,\quad\text{ in }\Omega \\\\
-\Delta\psi + \phi = 0,\quad\text{ with }\psi|\_{\partial\Omega}=0\text{ and }\phi|\_{\partial\Omega}=0
\\]
or in weak form:
\\[
\langle\nu\nabla\phi,\nabla\theta\rangle = \langle f, \theta\rangle,\quad\forall\theta\in V\_0 \\\\
\langle\nabla\psi,\nabla\theta'\rangle + \langle\phi,\theta'\rangle = 0,\qquad\forall\theta'\in V'\_0
\\]
with solution components `\((\psi\in V_0, \phi\in V'_0)\)`.
---
# The biharmonic equation
Source-code: see `biharmonic1.cc`
\\[
\overbrace{\langle\nu\nabla\phi,\nabla\theta\rangle}^{\text{opL\_phi}} = \overbrace{\langle f, \theta\rangle}^{\text{opF}}, \\\\
\underbrace{\langle\nabla\psi,\nabla\theta'\rangle}\_{\text{opL\_psi}} + \underbrace{\langle\phi,\theta'\rangle}\_{\text{opM\_phi}} = 0
\\]
```
Operator opL_phi(prob.getFeSpace(0), prob.getFeSpace(1));
addSOT(opL_phi, nu);
Operator opF(prob.getFeSpace(0));
addZOT(opF, f);
Operator opL_psi(prob.getFeSpace(1), prob.getFeSpace(0));
addSOT(opL_psi, 1.0);
Operator opM_phi(prob.getFeSpace(1), prob.getFeSpace(1));
addZOT(opM_phi, 1.0);
// ===== add operators to problem =====
prob.addMatrixOperator(opL_phi, 0, 1); // -nu*laplace(phi)
prob.addMatrixOperator(opL_psi, 1, 0); // -laplace(psi)
prob.addMatrixOperator(opM_phi, 1, 1); // phi
prob.addVectorOperator(opF, 0); // f(x)
/// ...
```
---
# The biharmonic equation
Source-code: see `biharmonic1.cc`
```
/// ...
BoundaryType nr = 1;
prob.addDirichletBC(nr, 1, 0, new Constant(0.0));
prob.addDirichletBC(nr, 0, 1, new Constant(0.0));
AdaptInfo adaptInfo("adapt");
AdaptStationary("adapt", prob, adaptInfo).adapt();
prob.writeFiles(&amp;adaptInfo, true);
```
Init-file: see `biharmonic1.dat.2d`
```matlab
biharmonic->mesh: mesh
biharmonic->components: 2
biharmonic->feSpace[0]: P1
biharmonic->feSpace[1]: P1
biharmonic->dim: 2
biharmonic->name: [psi,phi]
```
---
# The biharmonic equation
Velocity field `\(\mathbf{v}\)` and stream function `\(\psi\)` for `\(\mathbf{f}=(y, -x)^\top\)`
.center[
<img src="images/biharmonic.png" width="80%" />
]
---
class: center, middle
# Exercise3
## Biharmonic equation
---
# Exercise3: System of equations
Implement the 4th order PDE
\\[
u - \Delta(u -\epsilon\Delta u) = f(x)\quad\text{ in }\Omega, \\\\
\partial\_n u|\_{\partial\Omega} = \partial\_n(u - \epsilon\Delta u)|\_{\partial\Omega} = 0
\\]
in a rectangular domain `\(\Omega\)`, with `\(\epsilon<1\)` and `\(f(x)\in[-1,1]\)` (random values).
1. Formulate Splitting of equation in system of equations.
2. Assemble and solve the system.
3. Vary the value `\(\epsilon\)` from 1 to `\(10^{-3}\)`. Introduce a parameter and read it from init-file. Change the mesh resolution accordingly.
---
# Advanced Exercise3: Error estimators
1. Modify exercise 1 to allow local mesh adaption by an error estimator.
2. Change the init-file correspondingly, i.e. add an entry for estimator and marker.
3. Choose a residual estimator/ recovery estimator and vary the estimator parameters
4. Vary the marking strategies and visualize the effect.
---
# Some hints
### Used functions/classes:
```
// functor that returns random values in mean +- amplitude/2
Random(mean, amplitude);
EXPRESSION: eval(&amp;functor)
```
### Parameters to modify:
```
poisson->estimator[0]: residual / recovery
poisson->marker[0]->strategy: 1/2/3/4
```
### References:
- ALBERT manual, Section 1.5.2 (marker strategies)
</textarea>
<script src="lib/remark.js" type="text/javascript"></script>
<script type="text/javascript" src="../MathJax/MathJax.js?config=TeX-AMS_HTML"></script>
<script type="text/javascript">
var slideshow = remark.create({
ratio: "4:3",
highlightLanguage: "cpp"
});
// Setup MathJax
MathJax.Hub.Config({
tex2jax: {
skipTags: ['script', 'noscript', 'style', 'textarea', 'pre']
}
});
MathJax.Hub.Queue(function() {
$(MathJax.Hub.getAllJax()).map(function(index, elem) {
return(elem.SourceElement());
}).parent().addClass('has-jax');
});
MathJax.Hub.Configured();
</script>
</body>
</html>
<!DOCTYPE html>
<html>
<head>
<title>AMDiS - Adaptive Multi-Dimensional Simulations</title>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"/>
<style type="text/css">
@import url(https://fonts.googleapis.com/css?family=Yanone+Kaffeesatz);
@import url(https://fonts.googleapis.com/css?family=Raleway);
@import url(https://fonts.googleapis.com/css?family=Ubuntu);
@import url(https://fonts.googleapis.com/css?family=Droid+Serif:400,700,400italic);
@import url(https://fonts.googleapis.com/css?family=Ubuntu+Mono:400,700,400italic);
</style>
<!--<link rel="stylesheet" type="text/css" href="style_display.css" />-->
<link rel="stylesheet" type="text/css" href="style_print.css" />
</head>
<body>
<textarea id="source">
class: center, middle
# Session 7
## Boundary conditions and composite FEM
---
# Boundary regions
Parts of the boundary can be identified by boundary-nr given in the macro-file
<table width="100%"><tr>
<td width="60%"><pre><code>element vertices:
0 1 4
1 2 4...
element boundaries:
0 0 1
0 0 2
0 0 3
0 0 4
vertex coordinates:
0.0 0.0
1.0 0.0...
</code></pre>
</td><td width="40%" style="text-align:center;">
<img src="images/ellipt_macro_boundary.png" width="80%" alt="A macro mresh" />
</td>
</tr></table>
Order of the values? i'th entry identifies the face opposite to the i'th vertex of the element.
Nr. **0 means** *inner face*, **negative numbers** identify *periodic faces*, **positive numbers** all other.
--
```
prob.addDirichletBC(boundaryNr, 0, 0, functor);
```
Boundary number is stored in macro-mesh and is propagated during refinement.
---
# Set Boundary regions manually
Using an indicator function **`b:WorldVector<double> -> int`** for the face-centers
```
const int dim = prob.getMesh()->getDim();
for (auto* m : prob.getMesh()->getMacroElements())
{
for (int i = 0; i < dim+1; ++i) {
WorldVector<double> center; center.set(0.0);
for (int j = 0; j < dim+1; ++j) {
if (j == i)
continue;
center += m->getCoord(j);
}
center /= dim;
int boundaryNr = b(center);
if (m->getBoundary(i) != 0)
m->setBoundary(i, boundaryNr);
}
}
```
---
# Set Boundary regions manually
Using an element-vector
```
map<int, vector<int> > elementBoundaries = /*...*/;
for (auto* m : prob.getMesh()->getMacroElements())
{
auto const& boundaryNr = elementBoundaries[m->getIndex()];
for (int i = 0; i < boundaryNr.size(); ++i)
m->setBoundary(i, boundaryNr[i]);
}
```
*Problem*: Macro-mesh only indentifies boundary faces, not boundary vertices directly.
---
# DirichletBC based on indicator function
Introduce AbstractFunction wrapper
```
template <class ReturnType, class... ArgumentTypes>
class LambdaFunction
: public AbstractFunction<ReturnType, ArgumentTypes...>
{
public:
template <class F>
LambdaFunction(F&& f) : fct(std::forward<F>(f)) {}
virtual ReturnType operator()(ArgumentTypes const&... args) const override
{
return fct(args...);
}
private:
std::function<ReturnType(ArgumentTypes...)> fct;
};
```
---
# DirichletBC based on indicator function
Introduce AbstractFunction wrapper
```
template <class ReturnType, class... ArgumentTypes>
class LambdaFunction;
```
Specify boundary belonging / or boundary value based on coordinates:
```
prob.addDirichletBC(new LambdaFunction<bool, WorldVector<double>>(
[](WorldVector<double> const& x)
{
return x[0] < 1.e-8 || x[0] > 1.0-1.e-8;
}),
0, 0, new LambdaFunction<double, WorldVector<double>>(
[](WorldVector<double> const& x)
{
return std::sin(x[1]);
})
);
```
--
Here, `prob` must be of type `ExtendedProblemStat`.
---
# DirichletBC based on indicator function
Introduce AbstractFunction wrapper
```
template <class ReturnType, class... ArgumentTypes>
class LambdaFunction;
```
Specify boundary value based on coordinates:
```
prob.addDirichletBC(boundaryNr,
0, 0, [](WorldVector<double> const& x)
{
return std::sin(x[1]);
}
);
```
Here, no wrapper-functor necessary and standard `ProblemStat` possible.
---
# Periodic boundary conditions
Mesh is part of an infinite domain, with local periodic solution.
.center[
<img src="images/periodic_boundaries.png" width="90%" alt="Periodic boundaries" />
]
---
# Periodic boundary conditions
Identify boundary face of periodic boundaries, e.g. left-right periodic:
.center[
<img src="images/periodic_macro.png" width="80%" alt="Periodic boundaries" />
]
---
# Periodic boundary conditions
Identify boundary face of periodic boundaries, e.g. left-right periodic:
<table width="100%"><tr>
<td width="60%"><pre><code>element boundaries:
-1 1 0
0 0 0
0 1 0
-1 0 0 ...
element neighbors:
3 -1 1
2 4 0
1 -1 3
0 6 2 ...
</code></pre>
</td><td width="40%" style="text-align:center;">