Tip

Before we start the exercise 1, please download (test_laplacianFoam) the modified Laplacian solver first. Then put it in the shared folder and compile (wmake) it in the docker container. This modified solver writes out detailed information on matrix coefficients into the log file, so that it is possible to check the derived values against the computed values.

OpenFoam implementation

Ok, now let’s do some practical tings, (0) generate mesh; (1) read mesh and do some useful calculation, e.g. cell volume, face area, …, etc; (2) discretize Laplacian term and get coefficients; (3) discretize transient term and get additional coefficients; (4) construct the final coefficient matrix and RHS; (5) solve the system of algebraic equations; (6) write solution to file.

Goal

Deeply look into OpenFOAM and understand how it works!

In order to better understand OpenFOAM’s logic and its work flow, we have to look at the basic structure of the source code of a basic solver, laplacianFoam.

Listing 27 Source code of laplacianFoam
 1#include "fvCFD.H"               // Basic head file of OF
 2#include "simpleControl.H"       // Basic head file of OF
 3int main(int argc, char *argv[]) // Typical c++ main control function
 4{
 5   #include "setRootCaseLists.H" // Do some case/file path-related thing
 6   #include "createTime.H"       // Create a time object: read controlDict, ...
 7   #include "createMesh.H"       // (2) Create mesh object: read mesh and do some useful calculation
 8   simpleControl simple(mesh);   // Time loop control object
 9   #include "createFields.H"     // (3) Read input data: T and D in the PDE
10
11   Info<< "\nCalculating temperature distribution\n" << endl;
12   while (simple.loop(runTime)) // do time loop
13   {
14      Info<< "Time = " << runTime.timeName() << nl << endl;
15      while (simple.correctNonOrthogonal())  // Non-orthogonal correction loop for the unstructured/non-orthogonal mesh
16      {
17         fvScalarMatrix TEqn                 // (5) construct the final coefficient matrix, include RHS (.source)
18            (
19               fvm::ddt(T)                   // (4) Discretization of the transient term, return a fvMatrix object
20               -
21               fvm::laplacian(DT, T)         // (3) Discretization of the Laplacian, return a fvMatrix object
22            );
23            TEqn.solve();                    // (6) Solve the system of algebraic equations
24      }
25      #include "write.H"                     // (7) Save solution to file.
26      Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
27            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
28            << nl << endl;
29   }
30   Info<< "End\n" << endl;
31   return 0;
32}
Table 1 Comparison between equation and code

Term

Equation

Code

Coefficients(ldu:b)

Reference

Transient

\(\frac{\partial T}{\partial t}\)

fvm::ddt(T)

d(\(FluxC\)), b(\(FluxC^o\))

equation (78)

Laplacian

\(\nabla \cdot D \nabla T\)

fvm::laplacian(DT, T)

d(\(\sum\limits_{F\sim NB(C)}a_F\)), u(\(a_F\) only internal faces)

equation (58), equation (65), equation (67)

Listing 28 createFields.H
 1volScalarField T //Read input data of T, include ICs and BCs
 2(
 3   IOobject
 4   (
 5      "T",
 6      runTime.timeName(),
 7      mesh,
 8      IOobject::MUST_READ,
 9      IOobject::AUTO_WRITE
10   ),
11   mesh
12);
13
14IOdictionary transportProperties // Read dictionary file of transportProperties
15(
16   IOobject
17   (
18      "transportProperties",
19      runTime.constant(),
20      mesh,
21      IOobject::MUST_READ_IF_MODIFIED,
22      IOobject::NO_WRITE
23   )
24);
25dimensionedScalar DT // Read the diffusivity D (constant) from transportProperties object
26(
27   transportProperties.lookup("DT")
28);
Listing 29 constant/transportProperties
1FoamFile
2{
3   version     2.0;
4   format      ascii;
5   class       dictionary;
6   location    "constant";
7   object      transportProperties;
8}
9DT              DT [0 2 -1 0 0 0 0] 4e-05;
Listing 30 0/T
 1FoamFile
 2{
 3   version     2.0;
 4   format      ascii;
 5   class       volScalarField;
 6   object      T;
 7}
 8dimensions      [0 0 0 1 0 0 0];
 9internalField   uniform 273;
10boundaryField
11{
12   left
13   {
14      type            fixedValue;
15      value           uniform 273;
16   }
17   right
18   {
19      type            fixedValue;
20      value           uniform 573;
21   }
22   "(top|bottom)"
23   {
24      type            zeroGradient;
25   }
26}

Step0, Generate mesh

Nothing special, just a meshing process. Here we use the OpenFOAM utility blockMesh to generate a regular mesh (Regular box case) same as Fig. 40.

Listing 31 blockMeshDict of a regular box mesh we shown above
 1/*--------------------------------*- C++ -*----------------------------------*\
 2| =========                 |                                                 |
 3| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
 4|  \\    /   O peration     | Version:  5                                     |
 5|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 6|    \\/     M anipulation  |                                                 |
 7\*---------------------------------------------------------------------------*/
 8FoamFile
 9{
10   version     2.0;
11   format      ascii;
12   class       dictionary;
13   object      blockMeshDict;
14}
15// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
16
17convertToMeters 1;
18xmin -0.05;    //variable definition
19xmax 0.05;
20ymin -0.015;
21ymax 0.015;
22Lz 0.01;
23vertices    //vertices definition
24(
25   ($xmin      $ymin   0)  //coordinate of vertex 0
26   ($xmax    $ymin   0)  //coordinate of vertex 1
27   ($xmax    $ymax   0)  //coordinate of vertex 2
28   ($xmin      $ymax   0)  //coordinate of vertex 3
29   ($xmin      $ymin   $Lz)//coordinate of vertex 4
30   ($xmax    $ymin   $Lz)//coordinate of vertex 5
31   ($xmax    $ymax   $Lz)//coordinate of vertex 6
32   ($xmin      $ymax   $Lz)//coordinate of vertex 7
33);
34blocks
35(
36   hex (0 1 2 3 4 5 6 7) (10 3 1) simpleGrading (1 1 1)
37);
38boundary
39(
40   left    //patch name
41   {
42      type patch ;
43      faces   //face list
44      (
45            (0 4 7 3)
46      );
47   }
48   right
49   {
50      type patch;
51      faces
52      (
53            (2 6 5 1)
54      );
55   }
56   top
57   {
58      type patch;
59      faces
60      (
61            (3 7 6 2)
62      );
63   }
64   bottom
65   {
66      type patch;
67      faces
68      (
69            (1 5 4 0)
70      );
71   }
72   frontAndBack    //patch name
73   {
74      type empty;
75      faces //face list
76      (
77            (0 3 2 1)   //back face
78            (4 5 6 7)   //front face
79      );
80   }
81);
82// ************************************************************************* //

Ok, let’s exploring the main steps of the laplacianFoam. Please download (test_laplacianFoam) and do the following practice/debug steps in test_laplacianFoam.C.

Step 1, Read mesh and input field

The basic properties of cells and faces, e.g. area, volume, face normal vector, will be evaluated after mesh reading, all these processes are happened in the mesh object. It means that after calling creatFields.H all these properties are evaluated and stored in the mesh object. Of course the temperature field object T (volScalarField) with BCs and ICs, and thermal diffusivity DT are also initialized from input data after calling creatFields.H. It should be noted that the part of Laplacian discretization coefficients are also calculated in this step. If the mesh is not changed during simulation time, the mesh related coefficients just need to be calculated one time.

Table 2 Mesh- and Field-related coefficients

Object

Equation

Code

Faces

Reference

mesh

\(1/\delta_{C\leftrightarrow F}\)

mesh.deltaCoeffs()

All faces

equation (58), equation (65), equation (67)

T

Dependents on BCs type

gradientInternalCoeffs (diagonal), gradientBoundaryCoeffs (source)

Boundary patch/faces

equation (65), equation (67)

Listing 32 Access mesh properties and field boundary properties
 1// (1). read data
 2#include "createMesh.H"
 3#include "createFields.H"
 4// 1.1 access internal faces
 5Info<<"\n\nAccess internal mesh"<<endl;
 6surfaceVectorField Cf = mesh.Cf();
 7surfaceVectorField Sf = mesh.Sf();
 8surfaceScalarField S = mesh.magSf();
 9forAll(Cf, iface)
10{
11   Info<<iface<<": face center "<<Cf[iface]<<endl;
12   Info<<iface<<": face area vector "<<Sf[iface]<<endl;
13   Info<<iface<<": face area "<<S[iface]<<endl;
14   Info<<iface<<": face delta coeff "<<mesh.deltaCoeffs()[iface]<<endl;
15   Info<<iface<<": coeff(D*magSf*deltacoeff) "<<mesh.deltaCoeffs()[iface]*DT*S[iface]<<endl;
16}
17// 1.2. access boundary mesh
18Info<<"\n\nAccess boundary mesh"<<endl;
19const fvBoundaryMesh& boundaryMesh = mesh.boundary();
20forAll(boundaryMesh, patchI)
21{
22   const fvPatch& patch = boundaryMesh[patchI];
23   forAll(patch, faceI)
24   {
25      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face center "<<patch.Cf()[faceI]<<endl;
26      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face area vector "<<patch.Sf()[faceI]<<endl;
27      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face area "<<patch.magSf()[faceI]<<endl;
28      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": face delta coeff "<<patch.deltaCoeffs()[faceI]<<endl;
29      Info<<"Patch "<<patch.name()<<" face "<<faceI<<": owner cell "<<patch.faceCells()[faceI]<<endl;
30   }
31}
32// 1.3. access boundary field, boundary field coefficients,
33forAll(T.boundaryField(), patchI)
34{
35   Info<<"Boundary patch: "<<mesh.boundary()[patchI].name()<<endl;
36   Info<<"Is coupled ? "<<mesh.boundary()[patchI].coupled()<<endl;
37   Info<<"gradientInternalCoeffs of field T "<<endl;
38   Info<<T.boundaryField()[patchI].gradientInternalCoeffs()<<endl; //Diagonal coeff [A]
39   Info<<"gradientBoundaryCoeffs of field T "<<endl;
40   Info<<T.boundaryField()[patchI].gradientBoundaryCoeffs()<<"\n"<<endl; //source coeff, [B]
41}
Listing 33 Internal mesh properties of the regular mesh shown in Fig. 40.
  146: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07Access internal mesh
  20: face center (-0.04 -0.01 0.005)
  30: face area vector (0.0001 0 0)
  40: face area 0.0001
  50: face delta coeff 100
  60: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
  71: face center (-0.045 -0.005 0.005)
  81: face area vector (0 0.0001 0)
  91: face area 0.0001
 101: face delta coeff 100
 111: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 122: face center (-0.03 -0.01 0.005)
 132: face area vector (0.0001 0 0)
 142: face area 0.0001
 152: face delta coeff 100
 162: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 173: face center (-0.035 -0.005 0.005)
 183: face area vector (0 0.0001 0)
 193: face area 0.0001
 203: face delta coeff 100
 213: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 224: face center (-0.02 -0.01 0.005)
 234: face area vector (0.0001 0 0)
 244: face area 0.0001
 254: face delta coeff 100
 264: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 275: face center (-0.025 -0.005 0.005)
 285: face area vector (0 0.0001 0)
 295: face area 0.0001
 305: face delta coeff 100
 315: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 326: face center (-0.01 -0.01 0.005)
 336: face area vector (0.0001 0 0)
 346: face area 0.0001
 356: face delta coeff 100
 366: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 377: face center (-0.015 -0.005 0.005)
 387: face area vector (0 0.0001 0)
 397: face area 0.0001
 407: face delta coeff 100
 417: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 428: face center (0 -0.01 0.005)
 438: face area vector (0.0001 0 0)
 448: face area 0.0001
 458: face delta coeff 100
 468: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 479: face center (-0.005 -0.005 0.005)
 489: face area vector (0 0.0001 0)
 499: face area 0.0001
 509: face delta coeff 100
 519: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 5210: face center (0.01 -0.01 0.005)
 5310: face area vector (0.0001 0 0)
 5410: face area 0.0001
 5510: face delta coeff 100
 5610: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 5711: face center (0.005 -0.005 0.005)
 5811: face area vector (0 0.0001 0)
 5911: face area 0.0001
 6011: face delta coeff 100
 6111: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 6212: face center (0.02 -0.01 0.005)
 6312: face area vector (0.0001 0 0)
 6412: face area 0.0001
 6512: face delta coeff 100
 6612: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 6713: face center (0.015 -0.005 0.005)
 6813: face area vector (0 0.0001 0)
 6913: face area 0.0001
 7013: face delta coeff 100
 7113: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 7214: face center (0.03 -0.01 0.005)
 7314: face area vector (0.0001 0 0)
 7414: face area 0.0001
 7514: face delta coeff 100
 7614: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 7715: face center (0.025 -0.005 0.005)
 7815: face area vector (0 0.0001 0)
 7915: face area 0.0001
 8015: face delta coeff 100
 8115: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 8216: face center (0.04 -0.01 0.005)
 8316: face area vector (0.0001 0 0)
 8416: face area 0.0001
 8516: face delta coeff 100
 8616: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 8717: face center (0.035 -0.005 0.005)
 8817: face area vector (0 0.0001 0)
 8917: face area 0.0001
 9017: face delta coeff 100
 9117: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 9218: face center (0.045 -0.005 0.005)
 9318: face area vector (0 0.0001 0)
 9418: face area 0.0001
 9518: face delta coeff 100
 9618: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
 9719: face center (-0.04 0 0.005)
 9819: face area vector (0.0001 0 0)
 9919: face area 0.0001
10019: face delta coeff 100
10119: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
10220: face center (-0.045 0.005 0.005)
10320: face area vector (0 0.0001 0)
10420: face area 0.0001
10520: face delta coeff 100
10620: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
10721: face center (-0.03 0 0.005)
10821: face area vector (0.0001 0 0)
10921: face area 0.0001
11021: face delta coeff 100
11121: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
11222: face center (-0.035 0.005 0.005)
11322: face area vector (0 0.0001 0)
11422: face area 0.0001
11522: face delta coeff 100
11622: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
11723: face center (-0.02 0 0.005)
11823: face area vector (0.0001 0 0)
11923: face area 0.0001
12023: face delta coeff 100
12123: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
12224: face center (-0.025 0.005 0.005)
12324: face area vector (0 0.0001 0)
12424: face area 0.0001
12524: face delta coeff 100
12624: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
12725: face center (-0.01 0 0.005)
12825: face area vector (0.0001 0 0)
12925: face area 0.0001
13025: face delta coeff 100
13125: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
13226: face center (-0.015 0.005 0.005)
13326: face area vector (0 0.0001 0)
13426: face area 0.0001
13526: face delta coeff 100
13626: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
13727: face center (0 0 0.005)
13827: face area vector (0.0001 0 0)
13927: face area 0.0001
14027: face delta coeff 100
14127: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
14228: face center (-0.005 0.005 0.005)
14328: face area vector (0 0.0001 0)
14428: face area 0.0001
14528: face delta coeff 100
14628: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
14729: face center (0.01 0 0.005)
14829: face area vector (0.0001 0 0)
14929: face area 0.0001
15029: face delta coeff 100
15129: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
15230: face center (0.005 0.005 0.005)
15330: face area vector (0 0.0001 0)
15430: face area 0.0001
15530: face delta coeff 100
15630: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
15731: face center (0.02 0 0.005)
15831: face area vector (0.0001 0 0)
15931: face area 0.0001
16031: face delta coeff 100
16131: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
16232: face center (0.015 0.005 0.005)
16332: face area vector (0 0.0001 0)
16432: face area 0.0001
16532: face delta coeff 100
16632: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
16733: face center (0.03 0 0.005)
16833: face area vector (0.0001 0 0)
16933: face area 0.0001
17033: face delta coeff 100
17133: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
17234: face center (0.025 0.005 0.005)
17334: face area vector (0 0.0001 0)
17434: face area 0.0001
17534: face delta coeff 100
17634: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
17735: face center (0.04 0 0.005)
17835: face area vector (0.0001 0 0)
17935: face area 0.0001
18035: face delta coeff 100
18135: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
18236: face center (0.035 0.005 0.005)
18336: face area vector (0 0.0001 0)
18436: face area 0.0001
18536: face delta coeff 100
18636: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
18737: face center (0.045 0.005 0.005)
18837: face area vector (0 0.0001 0)
18937: face area 0.0001
19037: face delta coeff 100
19137: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
19238: face center (-0.04 0.01 0.005)
19338: face area vector (0.0001 0 0)
19438: face area 0.0001
19538: face delta coeff 100
19638: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
19739: face center (-0.03 0.01 0.005)
19839: face area vector (0.0001 0 0)
19939: face area 0.0001
20039: face delta coeff 100
20139: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
20240: face center (-0.02 0.01 0.005)
20340: face area vector (0.0001 0 0)
20440: face area 0.0001
20540: face delta coeff 100
20640: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
20741: face center (-0.01 0.01 0.005)
20841: face area vector (0.0001 0 0)
20941: face area 0.0001
21041: face delta coeff 100
21141: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
21242: face center (0 0.01 0.005)
21342: face area vector (0.0001 0 0)
21442: face area 0.0001
21542: face delta coeff 100
21642: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
21743: face center (0.01 0.01 0.005)
21843: face area vector (0.0001 0 0)
21943: face area 0.0001
22043: face delta coeff 100
22143: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
22244: face center (0.02 0.01 0.005)
22344: face area vector (0.0001 0 0)
22444: face area 0.0001
22544: face delta coeff 100
22644: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
22745: face center (0.03 0.01 0.005)
22845: face area vector (0.0001 0 0)
22945: face area 0.0001
23045: face delta coeff 100
23145: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
23246: face center (0.04 0.01 0.005)
23346: face area vector (0.0001 0 0)
23446: face area 0.0001
23546: face delta coeff 100
23646: coeff(D*magSf*deltacoeff) ((100*DT)*0.0001) [0 2 -1 0 0 0 0] 4e-07
Listing 34 Boundary mesh properties of the regular mesh shown in Fig. 40.
  1Patch bottom face 9: owner cell 9Access boundary mesh
  2Patch left face 0: face center (-0.05 -0.01 0.005)
  3Patch left face 0: face area vector (-0.0001 0 0)
  4Patch left face 0: face area 0.0001
  5Patch left face 0: face delta coeff 200
  6Patch left face 0: owner cell 0
  7Patch left face 1: face center (-0.05 0 0.005)
  8Patch left face 1: face area vector (-0.0001 0 0)
  9Patch left face 1: face area 0.0001
 10Patch left face 1: face delta coeff 200
 11Patch left face 1: owner cell 10
 12Patch left face 2: face center (-0.05 0.01 0.005)
 13Patch left face 2: face area vector (-0.0001 0 0)
 14Patch left face 2: face area 0.0001
 15Patch left face 2: face delta coeff 200
 16Patch left face 2: owner cell 20
 17Patch right face 0: face center (0.05 -0.01 0.005)
 18Patch right face 0: face area vector (0.0001 0 0)
 19Patch right face 0: face area 0.0001
 20Patch right face 0: face delta coeff 200
 21Patch right face 0: owner cell 9
 22Patch right face 1: face center (0.05 0 0.005)
 23Patch right face 1: face area vector (0.0001 0 0)
 24Patch right face 1: face area 0.0001
 25Patch right face 1: face delta coeff 200
 26Patch right face 1: owner cell 19
 27Patch right face 2: face center (0.05 0.01 0.005)
 28Patch right face 2: face area vector (0.0001 0 0)
 29Patch right face 2: face area 0.0001
 30Patch right face 2: face delta coeff 200
 31Patch right face 2: owner cell 29
 32Patch top face 0: face center (-0.045 0.015 0.005)
 33Patch top face 0: face area vector (0 0.0001 0)
 34Patch top face 0: face area 0.0001
 35Patch top face 0: face delta coeff 200
 36Patch top face 0: owner cell 20
 37Patch top face 1: face center (-0.035 0.015 0.005)
 38Patch top face 1: face area vector (0 0.0001 0)
 39Patch top face 1: face area 0.0001
 40Patch top face 1: face delta coeff 200
 41Patch top face 1: owner cell 21
 42Patch top face 2: face center (-0.025 0.015 0.005)
 43Patch top face 2: face area vector (0 0.0001 0)
 44Patch top face 2: face area 0.0001
 45Patch top face 2: face delta coeff 200
 46Patch top face 2: owner cell 22
 47Patch top face 3: face center (-0.015 0.015 0.005)
 48Patch top face 3: face area vector (0 0.0001 0)
 49Patch top face 3: face area 0.0001
 50Patch top face 3: face delta coeff 200
 51Patch top face 3: owner cell 23
 52Patch top face 4: face center (-0.005 0.015 0.005)
 53Patch top face 4: face area vector (0 0.0001 0)
 54Patch top face 4: face area 0.0001
 55Patch top face 4: face delta coeff 200
 56Patch top face 4: owner cell 24
 57Patch top face 5: face center (0.005 0.015 0.005)
 58Patch top face 5: face area vector (0 0.0001 0)
 59Patch top face 5: face area 0.0001
 60Patch top face 5: face delta coeff 200
 61Patch top face 5: owner cell 25
 62Patch top face 6: face center (0.015 0.015 0.005)
 63Patch top face 6: face area vector (0 0.0001 0)
 64Patch top face 6: face area 0.0001
 65Patch top face 6: face delta coeff 200
 66Patch top face 6: owner cell 26
 67Patch top face 7: face center (0.025 0.015 0.005)
 68Patch top face 7: face area vector (0 0.0001 0)
 69Patch top face 7: face area 0.0001
 70Patch top face 7: face delta coeff 200
 71Patch top face 7: owner cell 27
 72Patch top face 8: face center (0.035 0.015 0.005)
 73Patch top face 8: face area vector (0 0.0001 0)
 74Patch top face 8: face area 0.0001
 75Patch top face 8: face delta coeff 200
 76Patch top face 8: owner cell 28
 77Patch top face 9: face center (0.045 0.015 0.005)
 78Patch top face 9: face area vector (0 0.0001 0)
 79Patch top face 9: face area 0.0001
 80Patch top face 9: face delta coeff 200
 81Patch top face 9: owner cell 29
 82Patch bottom face 0: face center (-0.045 -0.015 0.005)
 83Patch bottom face 0: face area vector (0 -0.0001 0)
 84Patch bottom face 0: face area 0.0001
 85Patch bottom face 0: face delta coeff 200
 86Patch bottom face 0: owner cell 0
 87Patch bottom face 1: face center (-0.035 -0.015 0.005)
 88Patch bottom face 1: face area vector (0 -0.0001 0)
 89Patch bottom face 1: face area 0.0001
 90Patch bottom face 1: face delta coeff 200
 91Patch bottom face 1: owner cell 1
 92Patch bottom face 2: face center (-0.025 -0.015 0.005)
 93Patch bottom face 2: face area vector (0 -0.0001 0)
 94Patch bottom face 2: face area 0.0001
 95Patch bottom face 2: face delta coeff 200
 96Patch bottom face 2: owner cell 2
 97Patch bottom face 3: face center (-0.015 -0.015 0.005)
 98Patch bottom face 3: face area vector (0 -0.0001 0)
 99Patch bottom face 3: face area 0.0001
100Patch bottom face 3: face delta coeff 200
101Patch bottom face 3: owner cell 3
102Patch bottom face 4: face center (-0.005 -0.015 0.005)
103Patch bottom face 4: face area vector (0 -0.0001 0)
104Patch bottom face 4: face area 0.0001
105Patch bottom face 4: face delta coeff 200
106Patch bottom face 4: owner cell 4
107Patch bottom face 5: face center (0.005 -0.015 0.005)
108Patch bottom face 5: face area vector (0 -0.0001 0)
109Patch bottom face 5: face area 0.0001
110Patch bottom face 5: face delta coeff 200
111Patch bottom face 5: owner cell 5
112Patch bottom face 6: face center (0.015 -0.015 0.005)
113Patch bottom face 6: face area vector (0 -0.0001 0)
114Patch bottom face 6: face area 0.0001
115Patch bottom face 6: face delta coeff 200
116Patch bottom face 6: owner cell 6
117Patch bottom face 7: face center (0.025 -0.015 0.005)
118Patch bottom face 7: face area vector (0 -0.0001 0)
119Patch bottom face 7: face area 0.0001
120Patch bottom face 7: face delta coeff 200
121Patch bottom face 7: owner cell 7
122Patch bottom face 8: face center (0.035 -0.015 0.005)
123Patch bottom face 8: face area vector (0 -0.0001 0)
124Patch bottom face 8: face area 0.0001
125Patch bottom face 8: face delta coeff 200
126Patch bottom face 8: owner cell 8
127Patch bottom face 9: face center (0.045 -0.015 0.005)
128Patch bottom face 9: face area vector (0 -0.0001 0)
129Patch bottom face 9: face area 0.0001
130Patch bottom face 9: face delta coeff 200
131Patch bottom face 9: owner cell 9
Listing 35 Boundary properties of field T of the regular mesh shown in Fig. 40.
 10()Boundary patch: left
 2Is coupled ? 0
 3gradientInternalCoeffs of field T 
 43(-200 -200 -200)
 5gradientBoundaryCoeffs of field T 
 63(54600 54600 54600)
 7
 8Boundary patch: right
 9Is coupled ? 0
10gradientInternalCoeffs of field T 
113(-200 -200 -200)
12gradientBoundaryCoeffs of field T 
133(114600 114600 114600)
14
15Boundary patch: top
16Is coupled ? 0
17gradientInternalCoeffs of field T 
1810{0}
19gradientBoundaryCoeffs of field T 
2010{0}
21
22Boundary patch: bottom
23Is coupled ? 0
24gradientInternalCoeffs of field T 
2510{0}
26gradientBoundaryCoeffs of field T 
2710{0}
28
29Boundary patch: frontAndBack
30Is coupled ? 0
31gradientInternalCoeffs of field T 
320()
33gradientBoundaryCoeffs of field T 
340()
../../_images/Coordinate_delta_internalcell_regularBox.svg

Fig. 45 Information of internal cell (\(C_{12}\))

../../_images/Coordinate_delta_boundary_regularBox.svg

Fig. 46 Information of boundary cell (\(C_{19}\))

Step 2, discretize Laplacian term

Because discretization coefficients matrix of Laplacian term is a symmetric matrix, so fvm::Laplacian(DT, T) will return a fvMatrix object only has diagonal and upper. What fvm::Laplacian actually did is evaluate (1) \(a_F\) (see equation (58)) for each internal faces, (2) \(a_C\) for each cells, which is the negative summation of \(a_F\), (3) store the field BCs-related coefficients as internalCoeffs and boundaryCoeffs, respectively. All of these are implemented in gaussLaplacianScheme.C . Note that the Gaussian scheme is the only choice for Laplacian discretization in OF.

Listing 36 Access fvm::Laplacian discretization
1fvScalarMatrix Laplacian(fvm::laplacian(DT, T));
2Info<<"fvm::laplacian(DT, T): "<<"\n"
3   <<"\tLower"<<Laplacian.lower()<<"\n"
4   <<"\tDiagonal"<<Laplacian.diag()<<"\n"
5   <<"\tUpper"<<Laplacian.upper()<<"\n"
6   <<"\tinternalCoeffs"<<Laplacian.internalCoeffs()<<"\n"
7   <<"\tboundaryCoeffs"<<Laplacian.boundaryCoeffs()<<"\n"
8   <<"\tSource"<<Laplacian.source()<<"\n"
9   <<endl;

Tip

The source coefficients come from BCs are not stored in the .source(), but in boundaryCoeffs. So if you print Laplacian.source(), it will display zero. The BCs-related source will be added into the TEqn.source when TEqn.solve() is calling. There is a protected member function named addBoundarySource will be called in solve() function.

Listing 37 Key source code of fvm::Laplacian (gaussLaplacianScheme.C before nonOrthogonal correcting)
 1template<class Type, class GType>
 2tmp<fvMatrix<Type>>
 3gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
 4(
 5   const surfaceScalarField& gammaMagSf,
 6   const surfaceScalarField& deltaCoeffs,
 7   const GeometricField<Type, fvPatchField, volMesh>& vf
 8)
 9{
10   tmp<fvMatrix<Type>> tfvm
11   (
12      new fvMatrix<Type>
13      (
14            vf,
15            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
16      )
17   );
18   fvMatrix<Type>& fvm = tfvm.ref();
19
20   fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
21   fvm.negSumDiag();
22
23   forAll(vf.boundaryField(), patchi)
24   {
25      const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
26      const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
27      const fvsPatchScalarField& pDeltaCoeffs =
28            deltaCoeffs.boundaryField()[patchi];
29
30      if (pvf.coupled())
31      {
32            fvm.internalCoeffs()[patchi] =
33               pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
34            fvm.boundaryCoeffs()[patchi] =
35               -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
36      }
37      else
38      {
39            fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
40            fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
41      }
42   }
43
44   return tfvm;
45}
Listing 38 fvm::Laplacian Coefficients of the regular mesh shown in Fig. 40.
  1)
  2fvm::laplacian(DT, T): 
  3	Lower
  447
  5(
  64e-07
  74e-07
  84e-07
  94e-07
 104e-07
 114e-07
 124e-07
 134e-07
 144e-07
 154e-07
 164e-07
 174e-07
 184e-07
 194e-07
 204e-07
 214e-07
 224e-07
 234e-07
 244e-07
 254e-07
 264e-07
 274e-07
 284e-07
 294e-07
 304e-07
 314e-07
 324e-07
 334e-07
 344e-07
 354e-07
 364e-07
 374e-07
 384e-07
 394e-07
 404e-07
 414e-07
 424e-07
 434e-07
 444e-07
 454e-07
 464e-07
 474e-07
 484e-07
 494e-07
 504e-07
 514e-07
 524e-07
 53)
 54
 55	Diagonal
 5630
 57(
 58-8e-07			//0
 59-1.2e-06
 60-1.2e-06
 61-1.2e-06
 62-1.2e-06
 63-1.2e-06		//5
 64-1.2e-06
 65-1.2e-06
 66-1.2e-06
 67-8e-07			//9
 68-1.2e-06		//10
 69-1.6e-06
 70-1.6e-06		//12
 71-1.6e-06
 72-1.6e-06
 73-1.6e-06
 74-1.6e-06
 75-1.6e-06
 76-1.6e-06
 77-1.2e-06		//19
 78-8e-07
 79-1.2e-06
 80-1.2e-06
 81-1.2e-06
 82-1.2e-06
 83-1.2e-06
 84-1.2e-06
 85-1.2e-06
 86-1.2e-06
 87-8e-07
 88)
 89
 90	Upper
 9147
 92(
 934e-07
 944e-07
 954e-07
 964e-07
 974e-07
 984e-07
 994e-07
1004e-07
1014e-07
1024e-07
1034e-07
1044e-07
1054e-07
1064e-07
1074e-07
1084e-07
1094e-07
1104e-07
1114e-07
1124e-07
1134e-07
1144e-07
1154e-07
1164e-07
1174e-07
1184e-07
1194e-07
1204e-07
1214e-07
1224e-07
1234e-07
1244e-07
1254e-07
1264e-07
1274e-07
1284e-07
1294e-07
1304e-07
1314e-07
1324e-07
1334e-07
1344e-07
1354e-07
1364e-07
1374e-07
1384e-07
1394e-07
140)
141
142	internalCoeffs
1435
144(
1453(-8e-07 -8e-07 -8e-07)
1463(-8e-07 -8e-07 -8e-07)
14710{0}
14810{0}
1490()
150)
151
152	boundaryCoeffs
1535
154(
1553(-0.0002184 -0.0002184 -0.0002184)
1563(-0.0004584 -0.0004584 -0.0004584)
15710{-0}
15810{-0}
1590()
160)
161
162	Source
16330
164(
1650
1661.20371e-36
167-1.20371e-36
168-1.20371e-36
1691.20371e-36
1701.20371e-36
1711.67048e-52
172-1.20371e-36
1732.40741e-36
174-1.66533e-19
1750
1760
1770
1780
1790
1800
1810
1820
1830
1841.66533e-19
1850
1860
1870
1880
1890
1900
1910
1920
1930
1940
195)
../../_images/Coordinate_Laplacian_internalcell_regularBox.svg

Fig. 47 Information of internal cell (\(C_{12}\))

../../_images/Coordinate_Laplacian_boundary_C19_regularBox.svg

Fig. 48 Information of boundary cell (\(C_{19}\))

../../_images/Coordinate_Laplacian_boundary_C9_regularBox.svg

Fig. 49 Information of boundary cell (\(C_{9}\))

../../_images/Coordinate_Laplacian_boundary_C0_regularBox.svg

Fig. 50 Information of boundary cell (\(C_{0}\))

../../_images/Coordinate_Laplacian_boundary_C10_regularBox.svg

Fig. 51 Information of boundary cell (\(C_{10}\))

../../_images/Coordinate_Laplacian_boundary_C5_regularBox.svg

Fig. 52 Information of boundary cell (\(C_{5}\))

Step 3, discretize transient term

For implicit discretization, fvm::ddt(T) will return a fvMatrix object contains diagonal coefficients and source. The coefficients depend on discretization scheme. For example Euler scheme, the diagonal coefficients are calculated from equation (74) and equation (76). All these are implemented in EulerDdtScheme.C.

Tip

There are 8 schemes for transient discretization in OpenFOAM

  1. CoEuler

  2. CrankNicolson

  3. Euler

  4. SLTS

  5. backward

  6. bounded

  7. localEuler

  8. steadyState

Listing 39 Access fvm::ddt discretization
1Info<<"fvm::ddt(T): "<<"\n"
2   <<"\tLower"<<ddt.lower()<<"\n"
3   <<"\tDiagonal"<<ddt.diag()<<"\n"
4   <<"\tUpper"<<ddt.upper()<<"\n"
5   <<"\tinternalCoeffs"<<ddt.internalCoeffs()<<"\n" //actually this is not necessary for fvm::ddt, this is definitely equal to zero
6   <<"\tboundaryCoeffs"<<ddt.boundaryCoeffs()<<"\n"
7   <<"\tSource"<<ddt.source()<<"\n"
8   <<endl;
Listing 40 Key source code of fvm::ddt (EulerDdtScheme.C )
 1template<class Type>
 2tmp<fvMatrix<Type>>
 3EulerDdtScheme<Type>::fvmDdt
 4(
 5   const GeometricField<Type, fvPatchField, volMesh>& vf
 6)
 7{
 8   tmp<fvMatrix<Type>> tfvm
 9   (
10      new fvMatrix<Type>
11      (
12            vf,
13            vf.dimensions()*dimVol/dimTime
14      )
15   );
16
17   fvMatrix<Type>& fvm = tfvm.ref();
18
19   scalar rDeltaT = 1.0/mesh().time().deltaTValue(); // 1/dt
20
21   fvm.diag() = rDeltaT*mesh().Vsc(); // Vc/dt (FluxC)
22
23   if (mesh().moving())
24   {
25      fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc0();
26   }
27   else
28   {
29      fvm.source() = rDeltaT*vf.oldTime().primitiveField()*mesh().Vsc(); // -T_old*Vc/dt
30   }
31
32   return tfvm;
33}

\(\Delta t = 0.05\ s\)

Listing 41 fvm::ddt Coefficients of the regular mesh shown in Fig. 40.
 1)fvm::ddt(T): 
 2	Lower47{0}
 3	Diagonal
 430
 5(
 62e-05
 72e-05
 82e-05
 92e-05
102e-05
112e-05
122e-05
132e-05
142e-05
152e-05
162e-05
172e-05
182e-05
192e-05
202e-05
212e-05
222e-05
232e-05
242e-05
252e-05
262e-05
272e-05
282e-05
292e-05
302e-05
312e-05
322e-05
332e-05
342e-05
352e-05
36)
37
38	Upper47{0}
39	internalCoeffs
405
41(
423{0}
433{0}
4410{0}
4510{0}
460()
47)
48
49	boundaryCoeffs
505
51(
523{0}
533{0}
5410{0}
5510{0}
560()
57)
58
59	Source
6030
61(
620.00546
630.00546
640.00546
650.00546
660.00546
670.00546
680.00546
690.00546
700.00546
710.00546
720.00546
730.00546
740.00546
750.00546
760.00546
770.00546
780.00546
790.00546
800.00546
810.00546
820.00546
830.00546
840.00546
850.00546
860.00546
870.00546
880.00546
890.00546
900.00546
910.00546
92)

Step 4, construct the final coefficient matrix and RHS

The final coefficient matrix is constructed by simply adding the matrix of Step 2, discretize Laplacian term and Step 3, discretize transient term.

  • The diagonal coefficients come from \(a_C\) of Laplacian term and \(FluxC\) of transient term

  • The off-diagonal coefficients only come from \(a_F (internal\ face)\) of Laplacian term

  • The RHS comes from \(c_F\) of Laplacian term when at boundary faces and \(FluxC^oT_C^o\) of transient term.

Listing 42 Final matrix coefficients of the regular mesh shown in Fig. 40 at the first time step.
  1)TEqn: 
  2	Lower
  347
  4(
  5-4e-07
  6-4e-07
  7-4e-07
  8-4e-07
  9-4e-07
 10-4e-07
 11-4e-07
 12-4e-07
 13-4e-07
 14-4e-07
 15-4e-07
 16-4e-07
 17-4e-07
 18-4e-07
 19-4e-07
 20-4e-07
 21-4e-07
 22-4e-07
 23-4e-07
 24-4e-07
 25-4e-07
 26-4e-07
 27-4e-07
 28-4e-07
 29-4e-07
 30-4e-07
 31-4e-07
 32-4e-07
 33-4e-07
 34-4e-07
 35-4e-07
 36-4e-07
 37-4e-07
 38-4e-07
 39-4e-07
 40-4e-07
 41-4e-07
 42-4e-07
 43-4e-07
 44-4e-07
 45-4e-07
 46-4e-07
 47-4e-07
 48-4e-07
 49-4e-07
 50-4e-07
 51-4e-07
 52)
 53
 54	Diagonal
 5530
 56(
 572.08e-05        //0
 582.12e-05
 592.12e-05
 602.12e-05
 612.12e-05
 622.12e-05        //5
 632.12e-05
 642.12e-05
 652.12e-05
 662.08e-05        //9
 672.12e-05        //10
 682.16e-05
 692.16e-05        //12
 702.16e-05
 712.16e-05
 722.16e-05
 732.16e-05
 742.16e-05
 752.16e-05
 762.12e-05        //19
 772.08e-05
 782.12e-05
 792.12e-05
 802.12e-05
 812.12e-05
 822.12e-05
 832.12e-05
 842.12e-05
 852.12e-05
 862.08e-05
 87)
 88
 89	Upper
 9047
 91(
 92-4e-07
 93-4e-07
 94-4e-07
 95-4e-07
 96-4e-07
 97-4e-07
 98-4e-07
 99-4e-07
100-4e-07
101-4e-07
102-4e-07
103-4e-07
104-4e-07
105-4e-07
106-4e-07
107-4e-07
108-4e-07
109-4e-07
110-4e-07
111-4e-07
112-4e-07
113-4e-07
114-4e-07
115-4e-07
116-4e-07
117-4e-07
118-4e-07
119-4e-07
120-4e-07
121-4e-07
122-4e-07
123-4e-07
124-4e-07
125-4e-07
126-4e-07
127-4e-07
128-4e-07
129-4e-07
130-4e-07
131-4e-07
132-4e-07
133-4e-07
134-4e-07
135-4e-07
136-4e-07
137-4e-07
138-4e-07
139)
140
141	internalCoeffs
1425
143(
1443(8e-07 8e-07 8e-07)
1453(8e-07 8e-07 8e-07)
14610{0}
14710{0}
1480()
149)
150
151	boundaryCoeffs
1525
153(
1543(0.0002184 0.0002184 0.0002184)
1553(0.0004584 0.0004584 0.0004584)
15610{0}
15710{0}
1580()
159)
160
161	Source
16230
163(
1640.00546
1650.00546
1660.00546
1670.00546
1680.00546
1690.00546
1700.00546
1710.00546
1720.00546
1730.00546
1740.00546
1750.00546
1760.00546
1770.00546
1780.00546
1790.00546
1800.00546
1810.00546
1820.00546
1830.00546
1840.00546
1850.00546
1860.00546
1870.00546
1880.00546
1890.00546
1900.00546
1910.00546
1920.00546
1930.00546
194)

Step 5, solve

For the Regular box case case, we can use PBiCG solver and DILU preconditioner.

Available preconditioner in OpenFOAM

  • diagonal : for symmetric & nonsymmetric matrices (not very effective)

  • DIC : Diagonal Incomplete Cholesky preconditioner for symmetric matrices

  • DILU : Diagonal Incomplete LU preconditioner for nonsymmetric matrices

  • FDIC : Fast Diagonal Incomplete Cholesky preconditioner

  • GAMG : Geometric Agglomerated algebraic MultiGrid preconditioner

Available solver in OpenFOAM

  • BICCG: Diagonal incomplete LU preconditioned BiCG solver

  • diagonalSolver: Solver for symmetric and nonsymmetric matrices

  • GAMG: Geometric Agglomerated algebraic Multi-Grid solver

  • ICC: Incomplete Cholesky Conjugate Gradient solver

  • PBiCG: Bi-Conjugate Gradient solver with preconditioner

  • PCG: Conjugate Gradient solver with preconditioner

  • smoothSolver: Iterative solver with run-time selectable smoother

Krylov-subspace solvers

  • CG: The Conjugate Gradient algorithm applies to systems where A is symmetric positive definite (SPD)

  • GMRES: The Generalized Minimal RESidual algorithm is the first method to try if A is not SPD.

  • BiCG: The BiConjugate Gradient algorithm applies to general linear systems, but the convergence can be quite erratic.

  • BiCGstab: The stabilized version of the BiConjugate Gradient algorithm.

Step 6, write

Jupyter notebook