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}
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.
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 :
Then we put a do
..while loop for the time iteration :
The time derivative discretization is done in black box; we simply wrote :
for the bilinear form, and :
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 :
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) :
Then compute the rhs according to the equation :
Eventually, compute the L2 and H1 error :
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 :
We have a very good convergence slope in P1. We should also get a good convergence in P2 and further.
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).