Loading

Solution of the time dependent heat equation

We present a program to solve the time dependent 2-D heat equation using
the P_{1} 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; }

Copyright © 1998-2018 Rachid Touzani