This directory shows how to perform verification of the 2D thermal solver using the Method of Manufactured Solutions.

Make sure you follow up with the 3D thermal cases and the mechanical cases as well.

The chosen domain is a square of size one with the lower left corner at $x=0$ and $y=0$ as in

SetFactory("OpenCASCADE");
Rectangle(1) = {0, 0, 0, 1, 1, 0};

Physical Surface("bulk", 1) = {1};
Physical Curve("left", 2) = {4};
Physical Curve("right", 3) = {2};
Physical Curve("bottom", 4) = {1};
Physical Curve("top", 5) = {3};

The manufactured solution $T(x,y)$ and thermal conductivity $k(x,y)$ are given in the FeenoX input file thermal-square.fee as

# MMS data, set T_mms(x) and k_mms(x) as desired
T_mms(x,y) = 1 + sin(2*x)^2 * cos(3*y)^2
k_mms(x,y) = 1 + x - 0.5*y

As explained below, a script run.sh will

1. perform a parametric sweep over different mesh sizes
2. compute the numerical error with respect to the manufactured solution
3. fit the slope of the logarithm of the error with respect to the logarithm of the mesh size
4. create a Markdown table with the slopes
5. create a Pyxplot file to plot the logarithmic errors

The binaries feenox, gmsh and maxima need to be available to be called from the script.

• FeenoX
• Gmsh
• Maxima
• Pyxplot (optional)

Run the script run.sh (and clean.sh if you want to re-run):

./clean.sh
./run.sh

By default this will run over a few combinations of

• BCs
  • dirichlet
  • neumann
• element types
  • tri3
  • tri6
  • quad4
  • quad8
  • quad9
• meshing algorithms
  • struct
  • frontal
  • dealaunay
• refinement factors $c$ such that the mesh size $h$ is proportional to $1/c$,
  • 4
  • 8
  • 12
  • 16
  • 20
  • 24

Feel free to edit run.sh and modify these combinations:

bcs="dirichlet neumann"
elems="tri3 tri6 quad4 quad8 quad9"
algos="struct frontal delaunay"
cs="4 8 12 16 20 24"

You can also edit thermal-square.fee and change the expressions for the temperature and the conductivity. The run.sh script calls Maxima to compute the appropriate volumetric source $q(x,y)$ and the Neumann boundary conditions.

While running you'll see 8 columns:

1. logarithm of the mesh size $\log(h)$
2. logarithm of the $L-2$ error $\log(e_2)$
3. logarithm of the $L-\infty$ error $\log(e_\infty)$
4. $c$ factor for mesh size $h \propto 1/c$
5. number of cells (volumetric elements)
6. number of nodes (degrees of freedom)
7. total memory in Mb
8. wall time in seconds

After running, along the output for the individual combinations, you'll get aggregated results:

• thermal-square-results.md a Markdown table with the slopes of each combination of BC, element type and meshing algorithm
• thermal-square.ppl a Pyxplot file to create plots of the $L_2$ errors vs. element size

If you want to keep on VTK output for each case, run with vtk as the argument:

./run.sh vtk

Each combination of BC, element type and meshing algorithm has a data file with the errors ($e_2$ and $e_\infty$) vs. the element size\ $h$, which is fitted with FeenoX and plotted with Pyxplot:

pyxplot thermal-square.ppl

There's also a markdown file with a table listing the order of convergence of each combination at thermal-square-results.md.