Demo: A 2-D transient diffusion code

This code is developed to solve 2-D transient (time dependent) heat transfer problems. We use here an alternative to the presentation of the previous demo code. More specifically, we show how to perform all important phases of a time dependent finite element code.

  • As usual, the program starts by including the header file OFELI.h that itself includes all kernel class definitions, then we include thermal dedicated classes described in Therm.h. We also declare the namespace to simplify the code writing:

    #include "OFELI.h"
#include "Therm.h"
using namespace OFELI;

  • We construct an instance of class IPF to manage project data. We extract useful information from this instance: The maximal time, the time step, a flag for saving results and another for controlling output. 

       IPF data("ttd2 - 1.1",argv[1]);
   theFinalTime = data.getMaxTime();
   theTimeStep = data.getTimeStep();
   int save_flag = data.getSave();
   int verbose = data.getVerbose();

  • We instantiate class IOField to store the solution for plotting: 

       IOField pf(data.getPlotFile(),IOField::OUT);

  • We now define the Mesh instance. We next instantiate vectors b and u that will contain respectively the right-hand side and the solution. For commodity we give a name to the solution vector: 

       Mesh ms(data.getMeshFile());
   Vect<double> b(ms), u(ms);
   u.setName("Temperature");

  • The initial solution is given through the class Prescription. 

       Prescription pr(ms,data.getDataFile());
   pr.get(INITIAL_FIELD,u);

  • We define vectors that contain boundary conditions, body and boundary forces. These are only intantiated here but will be assigned at each time step in the time loop: 

       Vect<double> bc(ms), body_f(ms), bound_f(ms);

  • We can now define the equation to solve: We choose the diffusion-convection equation in 2-D, solved by triangular finite elements and (class DC2DT3) and transmit initial condition to the equation instance. 

       DC2DT3 eq(ms);
   eq.setInput(INITIAL,u);

  • We start the time step loop using the macro TimeLoop: 

       TimeLoop {

  • For the current time step, we get (Dirichlet) boundary conditions, sources and fluxes (Neumann boundary conditions): 

          pr.get(BOUNDARY_CONDITION,bc,theTime);
      pr.get(SOURCE,body_f,theTime);
      pr.get(FLUX,bound_f,theTime);
    Note the use of the global variable theTime which is automatically updated each time step.

  • We next trasmit these data to the equation by using the member function setInput: 

          eq.setInput(BOUNDARY_CONDITION,bc);
      eq.setInput(SOURCE,body_f);
      eq.setInput(FLUX,bound_f);

  • We can now run the time step by using the member function run with the appropriate parameter 

          eq.run(TRANSIENT_ONE_STEP);
    Note that we didn't define a time integration scheme. This means that the Backward Euler scheme (Default value) is used.

  • It remains now to store the solution in a file if the parameter save_flag enables it. The procedure it to store the solution each multiple of save_flag step. 

          u.setTime(theTime);
      if (theStep%save_flag == 0)
         pf.put(u);
   }

    Note that we have stored in the solution vector the time value, so that this one can be retrieved for plotting.

      

    An example

    Let us run this program with the data presented in the following project file: 

    <?xml version="1.0" encoding="ISO-8859-1" ?>
<OFELI_File>
   <info>
      <title>Finite Element Mesh of a beam</title>
      <date>January 1, 2010</date>
      <author>R. Touzani</author>
   </info>
   <Project name="proj">
      <mesh_file>proj-5x10.m</mesh_file>
      <plot_file>proj.pl</plot_file>
      <time_step>0.01</time_step>
      <max_time>1.0</max_time>
      <verbose>0</verbose>
      <output>0</output>
      <save>1</save>
   </Project>
   <Prescription>
      <BoundaryCondition code="1">-tanh(10)*(exp(t)-1)</BoundaryCondition>
      <BoundaryCondition code="2"> tanh(10)*(exp(t)-1)</BoundaryCondition>
      <Source>tanh(10*y)*(exp(t)+200*(exp(t)-1)/(cosh(10*y)*cosh(10*y)))</Source>
   </Prescription>
</OFELI_File>

    In summary, this file looks mainly like the one in the previous example. We show here in addition how to prescribe a boundary condition by an algebraic expression. The library is equipped with the expression parser fparser. We use here the variables x, y, z, and t for space and time variables.