Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Pages
Parabolic equation example
Author
Christophe Prud'homme chris.nosp@m.toph.nosp@m.e.pru.nosp@m.dhom.nosp@m.me@fe.nosp@m.elpp.nosp@m..org
Date
2013-06-19

Description

This section is about another example easy to learn and understand : the parabolic equation. The aim is to solve equations like :

\begin{equation} \left\{ \begin{aligned} \dfrac{\partial u}{\partial t} - nu*\Delta u = f & \text{on}\; \Omega \;, \ u & = 0 & \text{on}\;\partial\Omega \;,\ \end{aligned} \right. \end{equation}

where \(u\in\Omega\) is the unknown "trial" function and \(\Omega\) the domain.

However this LaplacianParabolic application can also be used in the stationnary form, that is to say, solve the equation :

\begin{equation} \left\{ \begin{aligned} - nu*\Delta u = f & \text{on}\; \Omega \;, \ u & = 0 & \text{on}\;\partial\Omega \;,\ \end{aligned} \right. \end{equation}

Time Discretization

Implementation

The overall code is based on the laplacian.cpp application to which we added time discretization and other useful options. The time derivative will be discretized following a backward method called {Backward Differentiation Formula}; it is a finite differences method, used a lot in time discretization. Here, we don't go deeply in the implementation of BDF; a how-to about it can be get in Feel++ pdf tutorial.

BDF set up

The instationnary solving of this equation is not far different from solving the laplacian equation : in fact the laplacian equation is solved at each time iteration (called time-step). First we instantiate BDF objects :

// create the BDF structure
bdf_ptrtype M_bdf = bdf( _space=Xh, _name="u" );
// start the time at time-initial
M_bdf->start();
// create the finite difference polynome of the unknown
M_bdf->initialize(u);

Then we put a do..while loop for the time iteration :

do{
//set the linear form
form1=...
//set the bilinear form
form2=...
//solve
backend(_rebuild=true)->solve(...);
}while( M-bdf->isFinished() == false );

The time derivative discretization is done in black box; we simply wrote :

if( !steady )
a += integrate( _range=elements( mesh ), _expr = cst(M_bdf->polyDerivCoefficient( 0 ))*idt(u)*id(v) );

for the bilinear form, and :

if( !steady )
{
auto bdf_poly = M_bdf->polyDeriv();
l += integrate( _range = elements(mesh), _expr = idv( bdf_poly )*id(v));
}

for the linear form. In steady mode, those terms will not be added.

At the begining, one should have put all terms inside the time loop. But here, in our case, there are terms which are not time depending ones; so it is useless to have them inside the time loop and performing useless computations. The application will compute once and for all before the begining of the time loop and will add them to the temporal terms :

auto D = backend()->newMatrix( _test=Xh, _trial=Xh );
auto a = form2( _test=Xh, _trial=Xh, _matrix=D );
a = integrate( _range=elements( mesh ), _expr=nu*gradt( u )*trans( grad( v ) ) );
a += integrate( _range=markedfaces( mesh,"Dirichlet" ),
_expr= nu * ( -( gradt( u )*vf::N() )*id( v )
-( grad( v )*vf::N() )*idt( u )
+penaldir*id( v )*idt( u )/hFace() ) );
Warning
At the end of each solve, we should re-initialised $^{n}$, $^{n-1}$...; so don't forget the shifting operation :
M_bdf->shiftRight(u);
and also, don't forget to increment the time :
M_bdf->next();

Error class

The Error class allows users to define a test solution that has to be found by the application, to compute the Right Hand Side term of the equation and to compute L2 and H1 errors.

First we set the exact solution (if it is given in the .cfg file) and the associated parameter(s) :

if ( !params.empty() )
cvg->setParams ( params );
cvg->setSolution(exact, params);

Then compute the rhs according to the equation :

// if rhs should be computed
if(cvg->computedrhs())
{
// if it is not defined
if( rhs.empty() )
{
// compute it from the exact function
FEELPP_ASSERT( exact.size() ).error( "no exact function defined, no rhs function defined. Aborting...");
// with the edp of the parabolic equation
if ( !steady )
{
typedef typename boost::function<ex (ex u, std::vector<symbol> vars, std::vector<symbol> p)> t_edp_type;
t_edp_type t_edp = transient_edp();
cvg->setRhs(&t_edp);
}
// or with the normal laplacian equation
else
{
typedef typename boost::function<ex (ex u, std::vector<symbol> vars)> edp_type;
edp_type edp = steady_edp();
cvg->setRhs(&edp);
}
LOG(INFO) << "computed rhs is : " << cvg->getRhs() << "\n";
}
// else it is defined, so ok
else
{
cvg->setRhs();
LOG(INFO) << "rhs is : " << cvg->getRhs() << "\n";
}
std::cout << "rhs is : " << cvg->getRhs() << "\n";
}

Eventually, compute the L2 and H1 error :

double L2error = normL2( _range=elements( mesh ),_expr=( idv( u )-idv(gproj) ) );
double H1seminorm = math::sqrt( integrate( elements(mesh), (gradv(u) - gradg)*trans(gradv(u) - gradg) ).evaluate()(0,0) );
double H1error = math::sqrt( L2error*L2error + H1seminorm*H1seminorm);
Warning
Be careful, do NOT use the L2 and H1 error computation of the Error class, because they are written in cylindrical coordinates !

Results

To verify and validate the implementation, we made convergence study in stationnary and in temporal modes. The exact solution in the input is :

\begin{equation} \sin(\Pi(x-1)) \sin(\Pi\dfrac{y-0.05}{0.1})e^{-t} \end{equation}

(in the border, it is indeed equal to 0).

For the time error, we implemented this formula :

\begin{equation} E_{r} = \left( \Delta t \: \sum\limits_{t=t_i}^{t_f} \| u - u^n \|_{L^2(\Omega^{tn})}^2 \right)^{\frac{1}{2}} \end{equation}

Here we go with the graphs :

Stationnary study

We have a very good convergence slope in P1. We should also get a good convergence in P2 and further.

Laplacian_parabolic_L2_error.png
Laplacian_parabolic_H1_error.png

Temporal study

Here is something very interesting : depending the problem (the exact solution, the rhs, boundary conditions etc...), the theorical slope is not reached for "big" time-steps, but only when they are decreasing (BDF1, BDF2, BDF3); then, if the slope begins to behave wrongly, try to increase the precision of the geometry (BDF2 and BDF3).

Laplacian_parabolic_BDF1.png
Laplacian_parabolic_BDF2.png
Laplacian_parabolic_BDF3.png