Learning Feel++
:leveloffset:+1 = Generic Partial Differential Equations
:leveloffset:+1
The Laplacian
Problem statement
We are interested in this section in the conforming finite element approximation of the following problem:
\$\partial \Omega_D\$, \$\partial \Omega_N\$ and \$\partial \Omega_R\$ can be empty sets. In the case \$\partial \Omega_D =\partial \Omega_R = \emptyset\$, then the solution is known up to a constant. 
In the implementation presented later, \$\partial \Omega_D =\partial \Omega_N = \partial \Omega_R = \emptyset\$, then we set Dirichlet boundary conditions all over the boundary. The problem then reads like a standard laplacian with inhomogeneous Dirichlet boundary conditions: 
Variational formulation
We assume that \$f, h, l \in L^2(\Omega)\$. The weak formulation of the problem then reads:
Conforming Approximation
We now turn to the finite element approximation using Lagrange finite element. We assume \$\Omega\$ to be a segment in 1D, a polygon in 2D or a polyhedron in 3D. We denote \$V_\delta \subset H^1(\Omega)\$ an approximation space such that \$V_{g,\delta} \equiv P^k_{c,\delta}\cap H^1_{g,\Gamma_D}(\Omega)\$.
The weak formulation reads:
from now on, we omit \$\delta\$ to lighten the notations. Be careful that it appears both the geometrical and approximation level. 
Feel++ Implementation
In Feel++, \$V_{g,\delta}\$ is not built but rather \$P^k_{c,\delta}\$.
The Dirichlet boundary conditions can be treated using different techniques and we use from now on the elimination technique. 
We start with the mesh
auto mesh = loadMesh(_mesh=new Mesh<Simplex<FEELPP_DIM,1>>);
the keyword auto enables type inference, for more details see Wikipedia C++11 page.

Next the discretization setting by first defining Vh=Pch<k>(mesh)
\$\equiv P^k_{c,h}\$, then elements of Vh
and expressions f
, n
and g
given by command line options or configuration file.
auto Vh = Pch<2>( mesh );
auto u = Vh>element("u");
auto mu = doption(_name="mu");
auto f = expr( soption(_name="functions.f"), "f" );
auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression
auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression
auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression
auto g = expr( soption(_name="functions.g"), "g" );
auto v = Vh>element( g, "g" );
at the following line

the variational formulation is implemented below, we define the
bilinear form a
and linear form l
and we set strongly the
Dirichlet boundary conditions with the keyword on
using
elimination. If we don’t find Dirichlet
, Neumann
or Robin
in the
list of physical markers in the mesh data structure then we impose
Dirichlet boundary conditions all over the boundary.
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=f*id(v));
l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v));
l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v));
toc("l");
tic();
auto a = form2( _trial=Vh, _test=Vh);
a = integrate(_range=elements(mesh),
_expr=mu*gradt(u)*trans(grad(v)) );
a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v));
a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g );
//! if no markers Robin Neumann or Dirichlet are present in the mesh then
//! impose Dirichlet boundary conditions over the entire boundary
if ( !mesh>hasAnyMarker({"Robin", "Neumann","Dirichlet"}) )
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g );
toc("a");
tic();
//! solve the linear system, find u s.t. a(u,v)=l(v) for all v
if ( !boption( "nosolve" ) )
a.solve(_rhs=l,_solution=u);
toc("a.solve");
cout << "u_hg_L2=" << normL2(_range=elements(mesh), _expr=idv(u)g) << std::endl;
tic();
auto e = exporter( _mesh=mesh );
e>addRegions();
e>add( "u", u );
e>add( "g", v );
e>save();
toc("Exporter");
return 0;
}
We have the following correspondance:

next we solve the algebraic problem
//! solve the linear system, find u s.t. a(u,v)=l(v) for all v
if ( !boption( "nosolve" ) )
a.solve(_rhs=l,_solution=u);
next we compute the \$L^2\$ norm of \$u_\deltag\$, it could serve as an \$L^2\$ error if \$g\$ was manufactured to be the exact solution of the Laplacian problem.
cout << "u_hg_L2=" << normL2(_range=elements(mesh), _expr=idv(u)g) << std::endl;
and finally we export the results, by default it is in the ensight gold format and the files can be read with Paraview and Ensight. We save both \$u\$ and \$g\$.
auto e = exporter( _mesh=mesh );
e>addRegions();
e>add( "u", u );
e>add( "g", v );
e>save();
Testcases
The Feel++ Implementation comes with testcases in 2D and 3D.
circle
circle
is a 2D testcase where \$\Omega\$ is a disk whose boundary
has been split such that \$\partial \Omega=\partial \Omega_D \cup
\partial \Omega_N \cup \partial \Omega_R\$.
Here are some results we can observe after use the following command
cd Testcases/quickstart/circle
mpirun np 4 /usr/local/bin/feelpp_qs_laplacian_2d configfile circle.cfg
This give us some data such as solution of our problem or the mesh used in the application.
Solution \$u_\delta\$ 
Mesh 
feelpp2d and feelpp3d
This testcase solves the Laplacian problem in \$\Omega\$ an quadrangle or hexadra containing the letters of Feel++
feelpp2d
After running the following command
cd Testcases/quickstart/feelpp2d
mpirun np 4 /usr/local/bin/feelpp_qs_laplacian_2d configfile feelpp2d.cfg
we obtain the result \$u_\delta\$ and also the mesh
/images/Laplacian/TestCases/Feelpp2d/meshfeelpp2d.png[] 

Solution \$u_\delta\$ 
Mesh 
feelpp3d
We can launch this application with the current line
cd Testcases/quickstart/feelpp3d
mpirun np 4 /usr/local/bin/feelpp_qs_laplacian_3d configfile feelpp3d.cfg
When it’s finish, we can extract some informations
Solution \$u_\delta\$ 
Mesh 
Levelset
Having the possibility to determine where two regions meeting can be really useful in some scientific domains. That’s why the levelset method is born.
Levelset introduction
Levelset function
By using a scalar function \phi, define on all regions as the null value is obtained when it’s placed on an interface of two domains.
We denote \Omega_1 and \Omega_2 two domains with \Gamma the interface betwen them. Then \phi can be define as
\phi(\boldsymbol{x}) = \left\{ \begin{array}{cccc} \text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in &\Omega_1 \\ 0, & \boldsymbol{x}& \in &\Gamma\\ \text{dist}(\boldsymbol{x}, \Gamma), & \boldsymbol{x}& \in & \Omega_2 \end{array} \right.
with \text{dist}(\boldsymbol{x}, \Gamma ) = \underset{\boldsymbol{y} \; \in \; \Gamma}{\min}( \boldsymbol{x}  \boldsymbol{y} ).
This function \phi had also the following property \nabla\phi=1.
Moreover, the unit normal vector \boldsymbol{n} outgoing from the interface and the curvature \mathcal{\kappa} can be obtained from the levelset function.
\boldsymbol{n}=\frac{\nabla\phi}{\nabla\phi} \\ \mathcal{\kappa}=\nabla \cdot \boldsymbol{n}= \nabla \cdot \frac{\nabla\phi}{\nabla\phi}
Now we have exposed the levelset function, we need to define how the levelset will evolve and will spread into all the space. To do this, we use the following advection equation :
\partial_t\phi+\boldsymbol{u}\cdot\nabla\phi=0
where \boldsymbol{u} is an incompressible velocity field.
Heaviside and Dirac functions
We define also the regularized Heaviside function H_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq  \varepsilon,\\ \dfrac{1}{2} \left(1+\dfrac{\phi}{\varepsilon}+\dfrac{\sin\left(\dfrac{\pi \phi}{\varepsilon}\right)}{\pi}\right), & \varepsilon \leq \phi \leq \varepsilon, \\ 1, & \phi \geq \varepsilon \end{array} \right.
and the regularized Dirac function \delta_\varepsilon(\phi) = \left\{ \begin{array}{cc} 0, & \phi \leq  \varepsilon,\\ \displaystyle\dfrac{1}{2 \varepsilon} \left(1+\cos\left(\dfrac{\pi \phi}{\varepsilon}\right)\right), & \varepsilon \leq \phi \leq \varepsilon, \\ 0, & \phi \geq \varepsilon. \end{array} \right.
The first one gives a different value to each side of the interface ( here 0 in and 1 out ). The second one allow us to define quantities, with value different from 0 at the interface. A typical value of \varepsilon in literature is 1.5h where h is the mesh step size.
It should be noted that these functions allow us to determine respectively the volume and the surface of the interface by V^+_{\varepsilon} = \int_{\Omega} H_\varepsilon \\ S^{\Gamma}_{\varepsilon} = \int_{\Omega} \delta_\varepsilon
Solid rotation of a slotted disk
We describe the benchmark proposed by Zalesak.
Computer codes, used for the acquisition of results, are from Vincent Doyeux.
Problem Description
In order to test our interface propagation method, i.e. the levelset method \phi, we will study the rotation of a slotted disk into a square domain. The geometry can be represented as
We denote \Omega, the square domain [0,1]\times[0,1]. The center of the slotted disk is placed at (0.5,0.75).
To model the rotation, we will apply an angular velocity, centered in (0.5,0.5), as the disk is back to its initial position after t_f=628.
During this test, we observe three different errors to measure the quality of our method. With these values, two kinds of convergence will be studied : the time convergence, with different time step on an imposed grid and the space one, where the space discretization and the time step are linked by a relation. Several stabilization methods are used such as CIP ( Continuous Interior Penalty ) or SUPG ( StreamlineUpwing/PetrovGalerkin ).
Boundary conditions
We set a Neumann boundary condition on the boundary of the domain.
Initial conditions
The velocity is imposed as \boldsymbol{u}=\left( \frac{\pi}{314} (50y),\frac{\pi}{314} (x50) \right)
Here is the velocity look in the square domain
Inputs
The following table displays the various fixed and variables parameters of this testcase.
Name 
Description 
Nominal Value 
Units 
r 
disk radius 
0.15 
m 
l 
slot base 
0.05 
m 
h 
slot height 
0.25 
m 
t_f 
slotted disk rotation period 
628 
s 
Outputs
We observe during this benchmark three different errors.
First at all, the mass error, define by
e_{\text{m}} = \frac{ \left m_{\phi_f}  m_{\phi_0}\right }{m_{\phi_0}} = \frac{ \left \displaystyle \int_{\Omega} \chi( \phi_f < 0 )  \displaystyle \int_{\Omega} \chi( \phi_0 < 0 ) \right }{ \displaystyle \int_{\Omega} \chi( \phi_0 < 0 )}
where \chi is the characteristic function. \chi( f( \phi ) ) = \left\{ \begin{array}{rcl} 1 & \text{ if } & f( \phi ) \neq 0 \\ 0 & \text{ if } & f( \phi ) = 0 \end{array} \right.
However, mass is gain and loose at different emplacements on the mesh, and at the same time, with the level set method.
Secondly, the sign change error e_{\text{sc}} = \sqrt{ \int_\Omega \left( (1H_0)  (1H_f) \right)^2 }
with H_0=H_\epsilon(\phi_0) and H_f=H_\epsilon(\phi_f), H_\epsilon the smoothed Heaviside function of thickness 2ε.
This error is better to define the interface displacement. In fact, we can determine where \phi_0\phi_f<0, in other words where the interface has moved.
Finally, we define the classical L^2 error at the interface, as e_{L^2} = \sqrt{ \frac{1}{\displaystyle \int_\Omega \chi( \delta(\phi_0) > 0 ) } \int_\Omega (\phi_0  \phi_f)^2 \chi( \delta(\phi_0) > 0 ) }.
Discretization
Time convergence
For this case, we set a fixed grid with mesh step size h=0.04, and so 72314 degree of freedom on a \mathbb{P}^1.
Then, after the disk made one round, we measure the errors obtained from two different discretizations ( BDF2 and Euler ) and compared them.
We repeat this with several time step dt\in \{2.14, 1, 0.5, 0.25, 0.20\}.
Only one stabilization method is used : SUPG
Space convergence
We define the following relation, between time step and mesh step size : dt=C\frac{h}{U_{max}}
where C<1 constant and U_{max} the maximum velocity of \Omega.
From the definition of our velocity, U_{max} is reached at the farthest point from the center of \Omega. In this case, we have U_{max}=0.007, and we set C=0.8.
We use the BDF2 method for time discretization. As in time convergence, we wait one round of the disk to measure the errors and we repeat this test for these values of h: 0.32, 0.16, 0.08, 0.04.
We compare the results from different stabilization methods : CIP, SUPG, GLS ( GalerkinLeastSquares ) and SGS ( SubGrid Scale ).
Implementation
Results
Time convergence
dt 
e_{L^2} 
e_{sc} 
e_m 
2.14 
0.0348851 
0.202025 
0.202025 
1.00 
0.0187567 
0.147635 
0.147635 
0.5 
0.0098661 
0.10847 
0.10847 
0.25 
0.008791 
0.0782569 
0.0782569 
0.20 
0.00803373 
0.0670677 
0.0670677 
dt 
e_{L^2} 
e_{sc} 
e_m 
2.14 
0.0118025 
0.0906617 
0.0492775 
1.00 
0.00436957 
0.0445275 
0.0163494 
0.5 
0.00173637 
0.0216359 
0.0100621 
0.25 
0.001003 
0.0125971 
0.00354644 
0.20 
0.000949343 
0.0117449 
0.00317368 
Space convergence
stab 
h 
e_{L^2} 
e_{sc} 
e_m 
CIP 
0.32 
0.0074 
0.072 
0.00029 
0.16 
0.0046 
0.055 
0.00202 

0.08 
0.0025 
0.033 
0.00049 

0.04 
0.0023 
0.020 
0.00110 

SUPG 
0.32 
0.012 
0.065 
0.01632 
0.16 
0.008 
0.049 
0.07052 

0.08 
0.004 
0.030 
0.00073 

0.04 
0.001 
0.018 
0.00831 

GLS 
0.32 
0.013 
0.066 
0.02499 
0.16 
0.008 
0.051 
0.05180 

0.08 
0.004 
0.031 
0.00805 

0.04 
0.001 
0.019 
0.00672 

SGS 
0.32 
0.012 
0.065 
0.01103 
0.16 
0.008 
0.050 
0.07570 

0.08 
0.004 
0.030 
0.00084 

0.04 
0.001 
0.018 
0.00850 
Conclusion
Let’s begin with time convergence results. Tables shows us that sign change error is better to define the quality of the chosen scheme than the mass error. In fact, the loss of mass somewhere can be nullify by a gain of mass elsewhere. Sign change error shows half an order gain from Euler scheme to BDF2 one, as L^2 errors show us a gain of one order. For the slotted disk shape, BDF2 uses the two previous iterations to obtain the current result, while Euler only need the previous iteration. This explain why we can see an asymmetrical tendency in the first one.
As for space convergence, the different stabilization methods we used give us the same convergence rate equals to 0.6, with close error values, for the sign change error. For the L^2 error case, it’s not as evident as the previous one. Aside the CIP stabilization method with a 0.6 convergence rate, the others show us a convergence rate of 0.9.
Bibliography

[Zalesak] Steven T. Zalesak, Fully multidimensional fluxcorrected transport algorithms for fluids, Journal of Computational Physics, 1979.

[Doyeux] Vincent Doyeux, Modelisation et simulation de systemes multifluides. Application aux ecoulements sanguins., Physique Numérique [physics.compph], Université de Grenoble, 2014
:leveloffset:1
Computational Fluid Dynamics Benchmarks
Turek & Hron CFD Benchmark
Introduction
We implement the benchmark proposed by Turek and Hron, on the behavior of drag and lift forces of a flow around an object composed by a pole and a bar, see Figure Initial Geometry..
The software and the numerical results were initially obtained from Vincent Chabannes.
This benchmark is linked to the Turek Hron CSM and Turek Hron FSI benchmarks. 
Problem Description
We consider a 2D model representative of a laminar incompressible flow around an obstacle. The flow domain, named \$\Omega_f\$, is contained into the rectangle \$ \lbrack 0,2.5 \rbrack \times \lbrack 0,0.41 \rbrack \$. It is characterised, in particular, by its dynamic viscosity \$\mu_f\$ and by its density \$\rho_f\$. In this case, the fluid material we used is glycerine.
In order to describe the flow, the incompressible NavierStokes model is chosen for this case, define by the conservation of momentum equation and the conservation of mass equation. At them, we add the material constitutive equation, that help us to define \$\boldsymbol{\sigma}_f\$
The goal of this benchmark is to study the behavior of lift forces \$F_L\$ and drag forces \$F_D\$, with three different fluid dynamics applied on the obstacle, i.e on \$\Gamma_{obst}\$, we made rigid by setting specific structure parameters. To simulate these cases, different mean inflow velocities, and thus different Reynolds numbers, will be used.
Boundary conditions
We set

on \$\Gamma_{in}\$, an inflow Dirichlet condition : \$ \boldsymbol{u}_f=(v_{in},0) \$

on \$\Gamma_{wall}\$ and \$\Gamma_{obst}\$, a homogeneous Dirichlet condition : \$ \boldsymbol{u}_f=\boldsymbol{0} \$

on \$\Gamma_{out}\$, a Neumann condition : \$ \boldsymbol{\sigma}_f\boldsymbol{ n }_f=\boldsymbol{0} \$
Initial conditions
We use a parabolic velocity profile, in order to describe the flow inlet by \$ \Gamma_{in} \$, which can be express by
v_{cst} = 1.5 \bar{U} \frac{4}{0.1681} y \left(0.41y\right)
where \$\bar{U}\$ is the mean inflow velocity.
However, we want to impose a progressive increase of this velocity profile. That’s why we define
v_{in} = \left\{ \begin{aligned} & v_{cst} \frac{1\cos\left( \frac{\pi}{2} t \right) }{2} \quad & \text{ if } t < 2 \\ & v_{cst} \quad & \text{ otherwise } \end{aligned} \right.
With t the time.
Moreover, in this case, there is no source term, so \$f_f\equiv 0\$.
Inputs
The following table displays the various fixed and variables parameters of this testcase.
Name  Description  Nominal Value  Units  

\$l\$ 
elastic structure length 
\$0.35\$ 
\$m\$ 

\$h\$ 
elastic structure height 
\$0.02\$ 
\$m\$ 

\$r\$ 
cylinder radius 
\$0.05\$ 
\$m\$ 

\$C\$ 
cylinder center coordinates 
\$(0.2,0.2)\$ 
\$m\$ 

\$\nu_f\$ 
kinematic viscosity 
\$1\times 10^{3}\$ 
\$m^2/s\$ 

\$\mu_f\$ 
dynamic viscosity 
\$1\$ 
\$kg/(m \times s)\$ 

\$\rho_f\$ 
density 
\$1000\$ 
\$kg/m^3\$ 

\$f_f\$ 
source term 
0 
\$kg/(m^3 \times s)\$ 

\$\bar{U}\$ 
characteristic inflow velocity 

\$m/s\$ 
Outputs
As defined above, the goal of this benchmark is to measure the drag and lift forces, \$F_D\$ and \$F_L\$, to control the fluid solver behavior. They can be obtain from
(F_D,F_L)=\int_{\Gamma_{obst}}\boldsymbol{\sigma}_f \boldsymbol{ n }_f
where \$\boldsymbol{n}_f\$ the outer unit normal vector from \$\partial \Omega_f\$.
Discretization
To realize these tests, we made the choice to used \$P_N\$\$P_{N1}\$ TaylorHood finite elements, described by Chabannes, to discretize space. With the time discretization, we use BDF, for Backward Differentation Formulation, schemes at different orders \$q\$.
Solvers
Here are the different solvers ( linear and nonlinear ) used during results acquisition.
type 
gmres 
relative tolerance 
1e13 
max iteration 
1000 
reuse preconditioner 
false 
relative tolerance 
1e8 
steps tolerance 
1e8 
max iteration 
CFD1/CFD2 : 100  CFD3 : 50 
max iteration with reuse 
CFD1/CFD2 : 100  CFD3 : 50 
reuse jacobian 
false 
reuse jacobian rebuild at first Newton step 
true 
relative tolerance 
1e5 
max iteration 
1000 
max iteration with reuse 
CFD1/CFD2 : 100  CFD3 : 1000 
reuse preconditioner 
false 
reuse preconditioner rebuild at first Newton step 
false 
type 
lu 
package 
mumps 
Running the model
The configuration files are in toolboxes/fluid/TurekHron
. The different cases are implemented in the corresponding .cfg
files e.g. cfd1.cfg
, cfd2.cfg
and cfd3.cfg
.
The command line in feelpptoolboxes docker reads
$ mpirun np 4 /usr/local/bin/feelpp_toolbox_fluid_2d configfile cfd1.cfg
The result files are then stored by default in
feel/applications/models/fluid/TurekHron/"case_name"/"velocity_space""pression_space""Geometric_order"/"processor_used"
For example, for CFD2 case executed on \$12\$ processors, with a \$P_2\$ velocity approximation space, a \$P_1\$ pressure approximation space and a geometric order of \$1\$, the path is
feel/toolboxes/fluid/TurekHron/cfd2/P2P1G1/np_12
Results
Here are results from the different cases studied in this benchmark.
CFD1
\$\mathbf{N_{geo}}\$  \$\mathbf{N_{elt}}\$  \$\mathbf{N_{dof}}\$  Drag  Lift 

Reference Turek and Hron 
14.29 
1.119 

1 
9874 
45533 (\$P_2/P_1\$) 
14.217 
1.116 
1 
38094 
173608 (\$P_2/P_1\$) 
14.253 
1.120 
1 
59586 
270867 (\$P_2/P_1\$) 
14.262 
1.119 
2 
7026 
78758 (\$P_3/P_2\$) 
14.263 
1.121 
2 
59650 
660518 (\$P_3/P_2\$) 
14.278 
1.119 
3 
7026 
146057 (\$P_4/P_3\$) 
14.270 
1.120 
3 
59650 
1228831 (\$P_4/P_3\$) 
14.280 
1.119 
All the files used for this case can be found in this rep [geo file, config file, json file]
CFD2
\$\mathbf{N_{geo}}\$  \$\mathbf{N_{elt}}\$  \$\mathbf{N_{dof}}\$  Drag  Lift 

Reference Turek and Hron 
136.7 
10.53 

1 
7020 
32510 (\$P_2/P_1\$) 
135.33 
10.364 
1 
38094 
173608 (\$P_2/P_1\$) 
136.39 
10.537 
1 
59586 
270867 (\$P_2/P_1\$) 
136.49 
10.531 
2 
7026 
78758 (\$P_3/P_2\$) 
136.67 
10.548 
2 
59650 
660518 (\$P_3/P_2\$) 
136.66 
10.532 
3 
7026 
146057 (\$P_4/P_3\$) 
136.65 
10.539 
3 
59650 
1228831 (\$P_4/P_3\$) 
136.66 
10.533 
All the files used for this case can be found in this rep [geo file, config file, json file]
CFD3
As CFD3 is timedependent ( from BDF use ), results will be expressed as
mean ± amplitude [frequency]
where

mean is the average of the min and max values at the last period of oscillations.
mean=\frac{1}{2}(max+min)

amplitude is the difference of the max and the min at the last oscillation.
amplitude=\frac{1}{2}(maxmin)

frequency can be obtain by Fourier analysis on periodic data and retrieve the lowest frequency or by the following formula, if we know the period time T.
frequency=\frac{1}{T}
\$\mathbf{\Delta t}\$  \$\mathbf{N_{geo}}\$  \$\mathbf{N_{elt}}\$  \$\mathbf{N_{dof}}\$  \$\mathbf{N_{bdf}}\$  Drag  Lift 

0.005 
Reference Turek and Hron 
439.45 ± 5.6183[4.3956] 
−11.893 ± 437.81[4.3956] 
0.01 
1 
8042 
37514 (\$P_2/P_1\$) 
2 
437.47 ± 5.3750[4.3457] 
9.786 ± 437.54[4.3457] 
2 
2334 
26706 (\$P_3/P_2\$) 
2 
439.27 ± 5.1620[4.3457] 
8.887 ± 429.06[4.3457] 

2 
7970 
89790 (\$P_2/P_2\$) 
2 
439.56 ± 5.2335[4.3457] 
11.719 ± 425.81[4.3457] 
0.005 
1 
3509 
39843\$(P_3/P_2)\$ 
2 
438.24 ± 5.5375[4.3945] 
11.024 ± 433.90[4.3945] 
1 
8042 
90582 (\$P_3/P_2\$) 
2 
439.25 ± 5.6130[4.3945] 
10.988 ± 437.70[4.3945] 

2 
2334 
26706 (\$P_3/P_2\$) 
2 
439.49 ± 5.5985[4.3945] 
10.534 ± 441.02[4.3945] 

2 
7970 
89790 (\$P_3/P_2\$) 
2 
439.71 ± 5.6410[4.3945] 
11.375 ± 438.37[4.3945] 

3 
3499 
73440 (\$P_4/P_3\$) 
3 
439.93 ± 5.8072[4.3945] 
14.511 ± 440.96[4.3945] 

4 
2314 
78168 (\$P_5/P_4\$) 
2 
439.66 ± 5.6412[4.3945] 
11.329 ± 438.93[4.3945] 
0.002 
2 
7942 
89482 (\$P_3/P_2)\$ 
2 
439.81 ± 5.7370[4.3945] 
13.730 ± 439.30[4.3945] 
3 
2340 
49389 (\$P_4/P_3\$) 
2 
440.03 ± 5.7321[4.3945] 
13.250 ± 439.64[4.3945] 

3 
2334 
49266 (\$P_4/P_3\$) 
3 
440.06 ± 5.7773[4.3945] 
14.092 ± 440.07[4.3945] 
All the files used for this case can be found in this rep [geo file, config file, json file].
Geometrical Order
Add a section on geometrical order. 
Conclusion
The reference results, Turek and Hron, have been obtained with a time step \$\Delta t=0.05\$. When we compare our results, with the same step and \$\mathrm{BDF}_2\$, we observe that they are in accordance with the reference results.
With a larger \$\Delta t\$, a discrepancy is observed, in particular for the drag force. It can also be seen at the same time step, with a higher order \$\mathrm{BDF}_n\$ ( e.g. \$\mathrm{BDF}_3\$ ). This suggests that the couple \$\Delta t=0.05\$ and \$\mathrm{BDF}_2\$ isn’t enough accurate.
Bibliography

[TurekHron] S. Turek and J. Hron, Proposal for numerical benchmarking of fluidstructure interaction between an elastic object and laminar incompressible flow, Lecture Notes in Computational Science and Engineering, 2006.

[Chabannes] Vincent Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles [math.AP], Universitée de Grenoble, 2013.
Unresolved directive in CFD/README.adoc  include::MultiFluid/README.adoc[]
2D Drops Benchmark
This benchmark has been proposed and realised by Hysing. It allows us to verify our level set code, our NavierStokes solver and how they couple together.
Computer codes, used for the acquisition of results, are from Vincent Doyeux, with the use of Chabannes's NavierStokes code.
Problem Description
We want to simulate the rising of a 2D bubble in a Newtonian fluid. The bubble, made of a specific fluid, is placed into a second one, with a higher density. Like this, the bubble, due to its lowest density and by the action of gravity, rises.
The equations used to define fluid bubble rising in an other are the NavierStokes for the fluid and the advection one for the level set method. As for the bubble rising, two forces are defined :

The gravity force : \$\boldsymbol{f}_g=\rho_\phi\boldsymbol{g}\$

The surface tension force : \$\boldsymbol{f}_{st}=\int_\Gamma\sigma\kappa\boldsymbol{ n } \$
We denote \$ \Omega\times\rbrack0,3\rbrack \$ the interest domain with \$ \Omega=(0,1)\times(0,2) \$. \$\Omega\$ can be decompose into \$\Omega_1\$, the domain outside the bubble and \$\Omega_2\$ the domain inside the bubble and \$\Gamma\$ the interface between these two.
Durig this benchmark, we will study two different cases : the first one with a ellipsoidal bubble and the second one with a squirted bubble.
Boundary conditions

On the lateral walls, we imposed slip conditions

On the horizontal walls, no slip conditions are imposed : \$\boldsymbol{u}=0 \$
Initial conditions
In order to let the bubble rise, its density must be inferior to the density of the exterior fluid, so \$\rho_1>\rho_2\$
Inputs
The following table displays the various fixed and variables parameters of this testcase.
Name  Description  Nominal Value  Units 

\$\boldsymbol{g}\$ 
gravity acceleration 
\$(0,0.98)\$ 
\$m/s^2\$ 
\$l\$ 
length domain 
\$1\$ 
\$m\$ 
\$h\$ 
height domain 
\$2\$ 
\$m\$ 
\$r\$ 
bubble radius 
\$0.25\$ 
\$m\$ 
\$B_c\$ 
bubble center 
\$(0.5,0.5)\$ 
\$m\$ 
Outputs
In the first place, the quantities we want to measure are \$X_c\$ the position of the center of the mass of the bubble, the velocity of the center of the mass \$U_c\$ and the circularity \$c\$, define as the ratio between the perimeter of a circle and the perimeter of the bubble. They can be expressed by
After that, we interest us to quantitative points for comparison as \$c_{min}\$, the minimum of the circularity and \$t_{c_{min}}\$, the time needed to obtain this minimum, \$u_{c_{max}}\$ and \$t_{u_{c_{max}}}\$ the maximum velocity and the time to attain it, or \$y_c(t=3)\$ the position of the bubble at the final time step. We add a second maximum velocity \$u_{max}\$ and \$u_{c_{max_2}}\$ and its time \$t_{u_{c_{max_2}}}\$ for the second test on the squirted bubble.
Discretization
This is the parameters associate to the two cases, which interest us here.
Case 
\$\rho_1\$ 
\$\rho_2\$ 
\$\mu_1\$ 
\$\mu_2\$ 
\$\sigma\$ 
Re 
\$E_0\$ 
ellipsoidal bubble (1) 
1000 
100 
10 
1 
24.5 
35 
10 
squirted bubble (2) 
1000 
1 
10 
0.1 
1.96 
35 
125 
Implementation
Results
Test 1
Test 2
We describe the different quantitative results for the two studied cases.
h 
\$c_{min}\$ 
\$t_{c_{min}}\$ 
\$u_{c_{max}}\$ 
\$t_{u_{c_{max}}}\$ 
\$y_c(t=3)\$ 
lower bound 
0.9011 
1.8750 
0.2417 
0.9213 
1.0799 
upper bound 
0.9013 
1.9041 
0.2421 
0.9313 
1.0817 
0.02 
0.8981 
1.925 
0.2400 
0.9280 
1.0787 
0.01 
0.8999 
1.9 
0.2410 
0.9252 
1.0812 
0.00875 
0.89998 
1.9 
0.2410 
0.9259 
1.0814 
0.0075 
0.9001 
1.9 
0.2412 
0.9251 
1.0812 
0.00625 
0.8981 
1.9 
0.2412 
0.9248 
1.0815 
h 
\$c_{min}\$ 
\$t_{c_{min}}\$ 
\$u_{c_{max_1}}\$ 
\$t_{u_{c_{max_1}}}\$ 
\$u_{c_{max_2}}\$ 
\$t_{u_{c_{max_2}}}\$ 
\$y_c(t=3)\$ 
lower bound 
0.4647 
2.4004 
0.2502 
0.7281 
0.2393 
1.9844 
1.1249 
upper bound 
0.5869 
3.0000 
0.2524 
0.7332 
0.2440 
2.0705 
1.1380 
0.02 
0.4744 
2.995 
0.2464 
0.7529 
0.2207 
1.8319 
1.0810 
0.01 
0.4642 
2.995 
0.2493 
0.7559 
0.2315 
1.8522 
1.1012 
0.00875 
0.4629 
2.995 
0.2494 
0.7565 
0.2324 
1.8622 
1.1047 
0.0075 
0.4646 
2.995 
0.2495 
0.7574 
0.2333 
1.8739 
1.1111 
0.00625 
0.4616 
2.995 
0.2496 
0.7574 
0.2341 
1.8828 
1.1186 
Conclusion
Bibliography

[Hysing] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, and L. Tobiska, Quantitative benchmark computations of twodimensional bubble dynamics, International Journal for Numerical Methods in Fluids, 2009.

[Chabannes] V. Chabannes, Vers la simulation numérique des écoulements sanguins, Équations aux dérivées partielles. PhD thesis, Université de Grenoble, 2013.

[Doyeux] V. Doyeux, Modélisation et simulation de systèmes multifluides, Application aux écoulements sanguins, PhD thesis, Université de Grenoble, 2014.
3D Drop benchmark
The previous section described the strategy we used to track the interface. We couple it now to the Navier Stokes equation solver described in \cite{chabannes11:_high}. In the current section, we present a 3D extension of the 2D benchmark introduced in \cite{Hysing2009} and realised using Feel++ in \cite{Doyeux2012}.
Benchmark problem
The benchmark objective is to simulate the rise of a 3D bubble in a Newtonian fluid. The equations solved are the incompressible Navier Stokes equations for the fluid and the advection for the level set:
where \$\rho\$ is the density of the fluid, \$\nu\$ its viscosity, and \$\mathbf{g} \approx (0, 0.98)^T\$ is the gravity acceleration.
The computational domain is \$\Omega \times \rbrack0, T\rbrack \$ where \$\Omega\$ is a cylinder which has a radius \$R\$ and a heigth \$H\$ so that \$R=0.5\$ and \$H=2\$ and \$T=3\$. We denote \$\Omega_1\$ the domain outside the bubble \$ \Omega_1= \{\mathbf{x}  \phi(\mathbf{x})>0 \} \$, \$\Omega_2\$ the domain inside the bubble \$ \Omega_2 = \{\mathbf{x}  \phi(\mathbf{x})<0 \} stem:[ and stem:[\Gamma\$ the interface \$ \Gamma = \{\mathbf{x}  \phi(\mathbf{x})=0 \} \$. On the lateral walls and on the bottom walls, noslip boundary conditions are imposed, i.e. \$\mathbf{u} = 0\$ and \$\mathbf{t} \cdot (\nabla \mathbf{u} + (\nabla \mathbf{u})^T) \cdot \mathbf{n}=0\$ where \$\mathbf{n}\$ is the unit normal to the interface and \$\mathbf{t}\$ the unit tangent. Neumann condition is imposed on the top wall i.e. \$\dfrac{\partial \mathbf{u}}{\partial \mathbf{n}}=\mathbf{0}\$. The initial bubble is circular with a radius \$r_0 = 0.25\$ and centered on the point \$(0.5, 0.5, 0.)\$. A surface tension force \$\mathbf{f}_{st}\$ is applied on \$\Gamma\$, it reads : \$\mathbf{f}_{st} = \int_{\Gamma} \sigma \kappa \mathbf{n} \simeq \int_{\Omega} \sigma \kappa \mathbf{n} \delta_{\varepsilon}(\phi)\$ where \$\sigma\$ stands for the surface tension between the twofluids and \$\kappa = \nabla \cdot (\frac{\nabla \mathbf{\phi}}{\nabla \phi})\$ is the curvature of the interface. Note that the normal vector \$\mathbf{n}\$ is defined here as \$\mathbf{n}=\frac{\nabla \phi}{\nabla \phi}\$.
We denote with indices \$1\$ and \$2\$ the quantities relative to the fluid in respectively in \$\Omega_1\$ and \$\Omega_2\$. The parameters of the benchmark are \$\rho_1\$, \$\rho_2\$, \$\nu_1\$, \$\nu_2\$ and \$\sigma\$ and we define two dimensionless numbers: first, the Reynolds number which is the ratio between inertial and viscous terms and is defined as : \$Re = \dfrac{\rho_1 \sqrt{\mathbf{g} (2r_0)^3}}{\nu_1}\$; second, the E\"otv\"os number which represents the ratio between the gravity force and the surface tension \$E_0 = \dfrac{4 \rho_1 \mathbf{g} r_0^2}{\sigma}\$. The table below reports the values of the parameters used for two different test cases proposed in~\cite{Hysing2009}.
Tests 
\$\rho_1\$ 
\$\rho_2\$ 
\$\nu_1\$ 
\$\nu_2\$ 
\$\sigma\$ 
Re 
\$E_0\$ 
Test 1 (ellipsoidal bubble) 
1000 
100 
10 
1 
24.5 
35 
10 
Test 2 (skirted bubble) 
1000 
1 
10 
0.1 
1.96 
35 
125 
The quantities measured in \cite{Hysing2009} are \$\mathbf{X_c}\$ the center of mass of the bubble, \$\mathbf{U_c}\$ its velocity and the circularity. For the 3D case we extend the circularity to the sphericity defined as the ratio between the surface of a sphere which has the same volume and the surface of the bubble which reads \$\Psi(t) = \dfrac{4\pi\left(\dfrac{3}{4\pi} \int_{\Omega_2} 1 \right)^{\frac{2}{3}}}{\int_{\Gamma} 1}\$.
Simulations parameters
The simulations have been performed on the supercomputer SUPERMUC using 160 or 320 processors. The number of processors was chosen depending on the memory needed for the simulations. The table Numerical parameters used for the test 1 simulations: Mesh size, Number of processors, Time step, Average time per iteration, Total time of the simulation. summarize for the test 1 the different simulation properties and the table Mesh caracteristics: mesh size given, number of Tetrahedra, number of points, number of order 1 degrees of freedom, number of order 2 degrees of freedom give the carachteristics of each mesh.
h 
Number of processors 
\$\Delta t\$ 
Time per iteration (s) 
Total Time (h) 
0.025 
360 
0.0125 
18.7 
1.25 
0.02 
360 
0.01 
36.1 
3.0 
0.0175 
180 
0.00875 
93.5 
8.9 
0.015 
180 
0.0075 
163.1 
18.4 
0.0125 
180 
0.00625 
339.7 
45.3 
h 
Tetrahedra 
Points 
Order 1 
Order 2 
0.025 
73010 
14846 
67770 
1522578 
0.02 
121919 
23291 
128969 
2928813 
0.0175 
154646 
30338 
187526 
4468382 
0.015 
217344 
41353 
292548 
6714918 
0.0125 
333527 
59597 
494484 
11416557 
The NavierStokes equations are linearized using the Newton’s method and we used a KSP method to solve the linear system. We use an Additive Schwarz Method for the preconditioning (GASM) and a LU method as a sub preconditionner. We run the simulations looking for solutions in finite element spaces spanned by Lagrange polynomials of order \$(2,1,1)\$ for respectively the velocity, the pressure and the level set.
Results Test 1: Ellipsoidal bubble
Accordind to the 2D results we expect that the drop would became ellipsoid. The figure~\ref{subfig:elli_sh} shows the shape of the drop at the final time step. The contour is quite similar to the one we obtained in the two dimensions simulations. The shapes are similar and seems to converge when the mesh size is decreasing. The drop reaches a stationary circularity and its topology does not change. The velocity increases until it attains a constant value. Figure~\ref{subfig:elli_uc} shows the results we obtained with different mesh sizes.
Bibliography
Computational Solid Mechanics
Computational solid mechanics ( or CSM ) is a part of mechanics that studies solid deformations or motions under applied forces or other parameters.
Second Newton’s law
The second Newton’s law allows us to define the fundamental equation of the solid mechanic, as follows
It’s define here into a Lagrangian frame.
Variables, symbols and units
Notation 
Quantity 
Unit 
\$\rho_s^*\$ 
strucure density 
\$kg/m^3\$ 
\$\boldsymbol{\eta}_s\$ 
displacement 
\$m\$ 
\$\boldsymbol{F}_s\$ 
deformation gradient 

\$\boldsymbol{\Sigma}_s\$ 
second PiolaKirchhoff tensor 
\$N/m^2\$ 
\$f_s^t\$ 
body force 
\$N/m^2\$ 
Lamé coefficients
The Lamé coefficients are deducing from the Young’s modulus \$E_s\$ and the Poisson’s ratio \$\nu_s\$ of the material we work on and can be express
Axisymmetric reduced model
We interest us here to a 1D reduced model, named generalized string.
The axisymmetric form, which will interest us here, is a tube of length \$L\$ and radius \$R_0\$. It is oriented following the axis \$z\$ and \$r\$ represent the radial axis. The reduced domain, named \$\Omega_s^*\$ is represented by the dotted line. So the domain, where radial displacement \$\eta_s\$ is calculated, is \$\Omega_s^*=\lbrack0,L\rbrack\$.
We introduce then \$\Omega_s^{'*}\$, where we also need to estimate a radial displacement as before. The unique variance is this displacement direction.
The mathematical problem associated to this reduce model can be describe as
where \$\eta_s\$, the radial displacement that satisfy this equation, \$k\$ is the Timoshenko’s correction factor, and \$\gamma_v\$ is a viscoelasticity parameter. The material is defined by its density \$\rho_s^*\$, its Young’s modulus \$E_s\$, its Poisson’s ratio \$\nu_s\$ and its shear modulus \$G_s\$
At the end, we take \$ \eta_s=0\text{ on }\partial\Omega_s^*\$ as a boundary condition, which will fixed the wall to its extremities.
inlude::FSI/README.adoc[]
Thermal Building Environment
ISO 10211:2007 Thermal bridges in building construction
Introduction
ISO 10211:2007 sets out the specifications for a threedimensional and a twodimensional geometrical model of a thermal bridge for the numerical calculation of:

heat flows, in order to assess the overall heat loss from a building or part of it;

minimum surface temperatures, in order to assess the risk of surface condensation.
These specifications include the geometrical boundaries and subdivisions of the model, the thermal boundary conditions, and the thermal values and relationships to be used.
ISO 10211:2007 is based upon the following assumptions:

all physical properties are independent of temperature;

there are no heat sources within the building element.
ISO 10211:2007 can also be used for the derivation of linear and point thermal transmittances and of surface temperature factors. More information here.
Implementation
Only the 2D specifications have been implemented.
Running the testcase
$ mpirun np 4 /usr/local/bin/feelpp_thermodyn_2d configfile thermo2dCase2.cfg
inlude::HeatFluid/README.adoc[]
Optimization problems
This section presents some mathematical optimization problems. The following examples source codes are located in "doc/manual/opt/".
:leveloffset:1