Demo: A 2-D steady state diffusion-convection code

This code is developed to solve 2-D steady state diffusion-convection problems. We use here a user friendly implementation. That is, rather than expliciting vector and matrix construction, matrix and right-hand side assembly, ... we provide functions that perform all-in-one operations.

The Program

  • 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. The constructor enables giving the program name and release and the program argument is obtained here: 

       IPF data("stdc2 - 2.0",argv[1]);

  • The member function getVerbose() returns the parameter that controls message output. 

       if (data.getVerbose())
      ...

  • We now read mesh data and construct an instance of class Mesh and an instance of class Prescription to prescribe various data that we will describe later. Note that these data prescription are provided in the file data.getPrescriptionFile(): 

       Mesh ms(data.getMeshFile(),true);
   Prescription p(ms,data.getPrescriptionFile());

  • Other types of data are the vectors that will contain problem data. Here we choose a sparse storage in order to use an iterative solution method. Vector data are the solution u, Dirichlet boundary condition contained in bc and boundary fluxes (Neumann boundary conditions) in bound_f. Note that the constructors of theses vectors use the Mesh instance. The vectors contain then mesh information. 

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

  • Vectors u, bc and bound_f are initialized from member function get of the instance of class Prescription: 

       p.get(BOUNDARY_CONDITION,bc);
   p.get(SOURCE,body_f);
   p.get(FLUX,bound_f);

    Note that BOUNDARY_CONDITION, SOURCE and FLUX are integers defined in an enum variable. They enable selecting the type of condition to enforce.

  • Since we deal with a diffusion-convection problem, we have to provide a velocity vector. We first instantiate the velocity vector, check if the project file contains the flag for velocity. In this case, we read by using class IOField the vector v: 

       Vect<double> v(ms,2);
   if (data.getInteger("v_flag")) {
      IOField ff(data.getMeshFile(),data.getString("v_file"),ms,IOField::IN);
      ff.get(v);
   }

  • We can now declare an equation. We instantiate class DC2DT3 that solves diffusion-convection equations in 2-D using triangular P1 (3-node) elements. After instantiation, we transfer to the class various informations: the mesh, the vector u that will contain solution at nodes, boundary condition vector, source vector, and velocity vector. 

       DC2DT3 eq;
   eq.setMesh(ms);
   eq.setInput(SOLUTION,u);
   eq.setInput(BOUNDARY_CONDITION,bc);
   eq.setInput(SOURCE,body_f);
   eq.setInput(FLUX,bound_f);
   eq.setInput(VELOCITY_FIELD,v);

  • With all these data, we can choose some options for the solver. Namely, we define the terms in the equation (Diffusion and convection) and then choose an iterative solver (Gmres) coupled with the ILU preconditioner. We then solve the problem by invoking the member function run: 

       eq.setTerms(DIFFUSION|CONVECTION);
   eq.setSolver(GMRES_SOLVER,ILU_PREC);
   int it = eq.run();

  • We end the program by outputting the solution and saving it for plotting purposes: 

       if (data.getOutput() > 0)
      cout << u;
   if (data.getPlot()) {
      IOField pf(data.getPlotFile(),IOField::OUT);
      pf.put(u);
   }

  

An example

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

  <Project name="proj">
    <mesh_file>proj.m</mesh_file>
    <parameter label="v_flag" value="1"/>
    <parameter label="v_file">proj.v</parameter>
    <plot_file>proj.t</plot_file>
    <verbose>1</verbose>
    <output>1</output>
    <plot>1</plot>
  </Project>
  <Prescription>
    <BoundaryCondition code="1">1.0</BoundaryCondition>
  </Prescription>

Let us give a brief description of this file contents:

  • The project is called proj. The name of a project yields some defaults for parameters and file names.
  • The mesh file is called proj.m
  • We have a flag parameter that says that a velocity is present for the convection contribution. This parameter can be retrieved by the label v_flag
  • The vector containing the given velocity field is present in the file proj.v. This file name can be retrieved by the label proj.v
  • The tag plot_file refers to a string that will be accessed in the IPF class in order to define a file for storing the solution
  • The next three lines describe parameters for message outputting, solution output and the opportunity of storing the solution in a plot file
  • Once the Project tag is closed, we define a Prescription tag to define prescribed boundaty conditions. Here, the tag BoundaryCondition chooses the type of prescription. The option code says that nodes with code equal to 1 are concerned with the prescription. Then, the value 1.0 is prescribed for this code.