Lesson 4

Lesson 4 : A time-dependent problem

We introduce, in this lesson, new aspects of OFELI programmation :

We consider a time-dependent heat transfer problem that we solve by Backward Euler time stepping scheme.
We consider the case of Neumann Boundary conditions.
Data (problem parameters) are introduced by a data file using IPF format.

The Finite Element Code

The main program

Here is a description of the source code. The source file can be found in Example 4 in the examples directory.

We start, as usual, by including required headers and naming the appropriate namespace.
```
#include "OFELI.h"
#include "Therm.h"
using namespace OFELI;
```
The program will have as argument the name of the parameter data file :
```
int main(int argc, char *argv[])
{
```
We next declare pointers to classes Element and Side that will be used later.
```
   Element *el;
   Side    *sd;
```

We expand program arguments and declare an instance of class IPF for parameter file :

   if (argc <= 1) {
      cout << " Usage : ex4 <parameter_file>" << endl;
      exit(1)
   }
   IPF data(argv[1]);

Parameters max_time (maximum time value) and deltat (time step) are retrieved as IPF class members.
```
   double max_time = data.getMaxTime();
   double deltat = data.getTimeStep();
```
The mesh data file is also obtained by the same method.
```
   ms.Get(data.getMeshFile());
```
In the present example, we introduce boundary conditions through a user defined class. This may be optional for Dirichlet conditions but necessary for Neumann ones.
```
   User ud(ms);
```
Implementation of class User will be given later.
We declare matrix and vector data: first, the matrix a is declared as instance of class SkSMatrix. The vectors b, u and bc will contain respectively, alternatively the right-hand side and the current solution, the previous solution and prescribed Dirichlet boundary conditions.
```
   SkSMatrix<double> a(ms);
   Vect<double> b(ms.getNbDOF()), u(ms.getNbDOF()), bc(ms.getNbDOF());
```
Since, we are dealing with a transient problem, we need initial data. This is retrieved from class member setInitialData of class User :
```
   ud.setInitialData(u);
```
Before starting time stepping loop, we calculate the number of time steps and initialize time :
```
   int nb_step = int(max_time/deltat);
   double time = 0;
```

We start a loop over time steps :

   for (int step=1; step<=nb_step; step++) {

The first thing to do here is to update time value and initialize the right-hand side to zero since this one will be assembled.
```
      time += deltat;
      b = 0;
```
We next transmit the user data class instance ud the time value :
```
      ud.setTime(time);
```
In order to deal with a problem with time-dependent boundary condition we re-fill vector bc at this level.
```
      ud.setDBC(bc);
```
We write a loop over finite elements as in the previous lessons :
```
      for (ms.topElement(); (el=ms.getElement());) {
```
We use here class DC2DT3 with the constructor that involves time. Instance eq will then be used to build matrix and right-hand side.
```
         DC2DT3 eq(el,u,time);
```
The element matrix is constructed with capacity term (chosen here to be lumped) and diffusion term :
```
         eq.LCapacity(1./deltat);
         eq.Diffusion();
```
Note that capacity matrix is multiplied by the inverse of time step. This is necessary to implement the backward Euler scheme.
We assemble matrix and right-hand side (useless for the present example). Note that, since the matrix does not depend on time, it is assembled once and factorized once.
```
         if (step==1)
            a.Assembly(el,eq.A());
         b.Assembly(el,eq.b());
      }
```
The loop on elements is closed.
To deal with Neumann boundary conditions (involving boundary integrals), we have to loop over given sides. The loop looks like the one over elements :
```
      for (ms.topSide(); (sd=ms.getSide());) {
```
For each side (pointed by eq) we invoke a constructor that involves sides.
```
         DC2DT3 eq(sd,u,time);
```
We fill the side vector using the instance ud of class User. The function BoundaryRHS calculates the side integral.
```
         eq.BoundaryRHS(ud);
```
We assemble side vectors just like for elements and close the loop.
```
         b.Assembly(sd,eq.b());
      }
```
Once the linear system is assembled, we impose Dirichlet boundary conditions by a penalty techniques implemented in member function Prescribe :
```
      a.Prescribe(ms,b,bc,step-1);
```
As said before, factorization is carried out at the first time step only. Obviously, solution is called each time step.
```
      if (step == 1)
         a.Factor();
      a.Solve(b);
```
Now, vector b contains the solution. We copy it to u to store it as a previous solution.
```
      u = b;
```
We may want to output the solution each time step :
```
      cout << "\nSolution for time : " << time << endl << u;
   }
   return 0;
}
```
and then close the time stepping loop and the program.

A User defined class

We have now to implement class User that defines boundary conditions, initial conditions. The class is defined in file User.h in the Example 4 in the tutorial directory.

Of course, we start by including file OFELI.h and invoking the namespace :
```
#include "OFELI.h"
using namespace OFELI;
```
Class User derives from abstract class UserData.
```
class User : public UserData<double> {
```
This class has only public members and no attributes.
```
public :
```
We have a constructor that provides the mesh to the class : nothing to do, the parent class does the job for you.
```
   User(Mesh &mesh) : UserData<double>(mesh) { }
```
We define a membre class to gives a value to prescribe for boundary condition in function of node code, node coordinates, time value and degree of freedom : Here, we impose that a code 2 imposes the value 1.0. Any other code will impose the default value 0.0.
```
   double BoundaryCondition(const Point<double> &x, int code, double time=0., size_t dof=1)
   {
      double ret = 0.0;
      if (code == 2)
         ret = 1.0;
      return ret;
   }
```

The same scheme works for Neumann boundary condition :

   double SurfaceForce(const Point<double> &x, int code, double time, size_t dof)
   {
      if (code)
         return 1.0;
      else
         return 0.0;
   }

The class definition ends here
```
};
```
Let us finally note that since no implementation is given for initial condition, the default one is 0.0 for each degree of freedom.

A test

We use here exactly the same mesh file as in the previous lesson. Of course, we obtain the same solution. The convergence is obtained after 3 iterations.