Demo: Class TimeStepping: Solution of the time dependent heat equation

We present a program to solve the time dependent 2-D heat equation using the P1 finite element method for space discretization and one of the available time integration schemes for time discretization. To validate the code, we test the exact solution 

   u(x,y) = exp(-t)*exp(x+y)
This gives as a source term, in the case where all coefficients are equal to 1: 
   f(x,y) = -3*exp(-t)*exp(x+y)
  • We start, as usual, by including the main header file of the library, and the header for thermal problems. The program has, as argument the mesh file and the time step. We construct a Mesh instance with the flag that enables eliminating imposed boundary conditions. Note that the final time value is also set by using the OFELI global variable theFinalTime. 

    #include "OFELI.h"
#include "Therm.h"
using namespace OFELI;
                          
int main(int argc, char *argv[])
{
   theFinalTime = 1.;
   Mesh ms(argv[1],true);
   theTimeStep = atof(argv[2]);

  • We declare the vectors u, f and bc that are respectively the solution vector, the vector containing the source, and the vector containing imposed boundary values. We then introduce the initial condition by assigning the corresponding algebraic expression to the solution vector.
    We construct an instance of class DC2DT3 using the Mesh instance and choose the terms of the equation. Here the equation contains the (consistent) capacity, the diffusion and the source terms.
    We can then declare the time integration scheme by instantiating the class TimeStepping. Clearly, the backward (implicit) Euler scheme is used. To this instance, we transmit the equation reference and set initial solution vector. Note that the solution will be stored in this vector. We also choose the conjugate gradient method preconditioned by the diagonal ILU preconditioner. 

       Vect<double> u(ms), f(ms), bc(ms);
   u.set("exp(x+y)");
                          
   DC2DT3 eq(ms);
   eq.setTerms(CAPACITY|DIFFUSION|SOURCE);
                          
   TimeStepping ts(BACKWARD_EULER,theTimeStep,theFinalTime);
   ts.setPDE(eq);
   ts.setInitial(u);
   ts.setLinearSolver(CG_SOLVER,DILU_PREC);

  • We can now loop over time steps. For each time step, we give the source terms, which is time dependent and transfer it to the instance of TimeStepping. The same treatment is realized for the Dirichlet boundary condition. Finally we run the time step. Note that the vector u contains now the solution at the current time. 

       TimeLoop {
      f.setTime(theTime);
      f.set("-3*exp(-t)*exp(x+y)");
      ts.setRHS(f);                    
      bc.setTime(theTime);
      bc.setNodeBC(1,"exp(-t)*exp(x+y)");
      ts.setBC(bc);
      ts.runOneTimeStep();
   }

  • We end by printing the TimeStepping instance information. To compute error, we consider the vector U containing the exact solution at nodes, and compute an error norm. 

       cout << ts;
   Vect<double> U(ms,"u_ex",theFinalTime);
   U.set("exp(-1)*exp(x+y)");
   cout << "Error in L2-Norm: " << (u-U).getWNorm2() << endl;
   return 0;
}