Poisson Equation (03-poisson)

Watch the Video

Git reference: Tutorial example 03-poisson.

This example shows how to solve a simple PDE that describes stationary heat transfer in an object that is heated by constant volumetric heat sources (such as with a DC current). The temperature equals to a prescribed constant on the boundary. We will learn how to:

  • Define a weak formulation.
  • Initialize matrix solver.
  • Solve the discrete problem.
  • Output the solution and element orders in VTK format (to be visualized, e.g., using Paraview).
  • Visualize the solution using Hermes’ native OpenGL-based functionality.

Model problem

Consider the Poisson equation

(1)-\mbox{div}(\lambda \nabla u) = C_{src}

in the L-shaped domain \Omega from the previous examples. The equation is equipped with constant Dirichlet boundary conditions

(2)u = u_0\ \ \ \mbox{on}\ \partial \Omega.

Here u is an unknown temperature distribution, C_{src} a real number representing volumetric heat sources/losses, and \lambda > 0 is thermal conductivity of the material.

The weak formulation is derived as usual, first by multiplying equation (1) with a test function v, then integrating over the domain \Omega, and then applying the Green’s theorem (integration by parts) to the second derivatives. Because of the Dirichlet condition (2), there are no surface integrals. Since the product of the two gradients in the volumetric weak form needs to be integrable for all u and v in V, the proper space for the solution is V = H^1_0(\Omega). The weak formulation reads: Find u \in V such that

(3)\int_\Omega \lambda \nabla u \cdot \nabla v \;\mbox{d\bfx} = \int_\Omega C_{src} v \;\mbox{d\bfx}\ \ \ \mbox{for all}\ v \in V.

Hermes, however, needs the equation in the form

(4)\int_\Omega \lambda \nabla u \cdot \nabla v \;\mbox{d\bfx} - \int_\Omega C_{src} v \;\mbox{d\bfx} = 0\ \ \ \mbox{for all}\ v \in V.

Let us explain why.

Jacobian-residual formulation

Hermes always assumes that the problem is nonlinear, and by default it uses the Newton’s method to solve it. Other methods for the solution of nonlinear problems are available as well (to be discussed later).

For linear problems, the Newton’s method converges in one step. Well, at least that’s what the theory says.

In practice the Newton’s method may take more than one step for a linear problem if the matrix solver does not do a good job. This is not unusual at all, in particular when an iterative solver is used. By checking the residual of the equation, the Newton’s method always makes sure that the problem is solved correctly, or it fails in a clearly visible way. This is the reason No. 1 why the Newton’s method should be used even for problems that are linear.

Another reason is that a consistent approach to linear and nonlinear problems allows Hermes’ users to first formulate and solve a simplified linear version of the problem, and then extend it to a nonlinear version effortlessly. Let’s explain how this works.

Consistent approach to linear and nonlinear problems

First assume that \lambda = \lambda_{al} in \Omega_{al} and \lambda = \lambda_{cu} in \Omega_{cu} where both \lambda_{al} and \lambda_{cu} are constants. Then the problem is linear and the weak form for the Jacobian is

\int_{\Omega_{al}} \lambda_{al} \nabla u \cdot \nabla v \, \mbox{d}x \mbox{d}y + \int_{\Omega_{cu}} \lambda_{cu} \nabla u \cdot \nabla v \, \mbox{d}x \mbox{d}y,

where u stands for a basis function and v for a test function. The reader does not have to worry about the word “Jacobian” since for linear problems this is the same as “stiffness matrix”. Simply forget from the left-hand side of the weak formulation (4) all expressions that do not contain u. A detailed explanation of the Newton’s method for nonlinear problems will be provided at the beginning of the tutorial part P02.

The residual weak form is the entire left-hand side of (4) where u is now the approximate solution (not a basis function as above):

\int_{\Omega_{al}} \lambda_{al} \nabla u \cdot \nabla v \, \mbox{d}x \mbox{d}y + \int_{\Omega_{cu}} \lambda_{cu} \nabla u \cdot \nabla v \, \mbox{d}x \mbox{d}y - \int_{\Omega} C_{src} v \, \mbox{d}x \mbox{d}y.

This is the constructor of the corresponding weak formulation in Hermes:

CustomWeakFormPoisson::CustomWeakFormPoisson(std::string marker_al, HermesFunction* lambda_al,
                                             std::string marker_cu, HermesFunction* lambda_cu,
                                             HermesFunction* src_term) : WeakForm(1)
{
  // Jacobian forms.
  add_matrix_form(new WeakFormsH1::DefaultJacobianDiffusion(0, 0, marker_al, lambda_al));
  add_matrix_form(new WeakFormsH1::DefaultJacobianDiffusion(0, 0, marker_cu, lambda_cu));

  // Residual forms.
  add_vector_form(new WeakFormsH1::DefaultResidualDiffusion(0, marker_al, lambda_al));
  add_vector_form(new WeakFormsH1::DefaultResidualDiffusion(0, marker_cu, lambda_cu));
  add_vector_form(new WeakFormsH1::DefaultVectorFormVol(0, HERMES_ANY, src_term));
};

Here HERMES_ANY means that the volumetric vector form will be assigned to all material markers.

For constant LAMBDA_AL and LAMBDA_CU, the form is instantiated as follows:

CustomWeakFormPoisson wf("Aluminum", new HermesFunction(LAMBDA_AL), "Copper",
                         new HermesFunction(LAMBDA_CU), new HermesFunction(-VOLUME_HEAT_SRC));

To replace the constants with cubic splines, one just needs to do

CubicSpline lambda_al(...);
CubicSpline lambda_cu(...);
CustomWeakFormPoisson wf("Aluminum", &lambda_al, "Copper",
                         &lambda_cu, new HermesFunction(-VOLUME_HEAT_SRC));

If the reader guesses that CubicSpline is a descendant of HermesFunction, then this is a good guess:

class CubicSpline : public HermesFunction
...

The constant VOLUME_HEAT_SRC can be replaced with an arbitrary function of x and y by subclassing HermesFunction:

class CustomNonConstSrc : public HermesFunction
...

If cubic splines are not enough, then one can subclass HermesFunction to define arbitrary nonlinearities:

class CustomLanbdaAl : public HermesFunction
...
class CustomLanbdaCu : public HermesFunction
...

Note that we are able to extend a linear problem to a nonlinear one without touching the CustomWeakFormPoisson class. This is also how all default weak forms in Hermes work.

In the rest of part P01 we will focus on linear problems.

Default Jacobian for the diffusion operator

Hermes provides default weak forms for many common PDE operators. The above default weak forms DefaultJacobianDiffusion, DefaultResidualDiffusion and DefaultVectorFormVol can be found in the file weakforms_h1.h. To begin with, the line

add_matrix_form(new DefaultJacobianDiffusion(0, 0, marker_al, lambda_al));

adds to the Jacobian weak form the integral

\int_{\Omega_{al}} \lambda_{al} \nabla u \cdot \nabla v \, \mbox{d}x \mbox{d}y

where u is a basis function and v a test function.

It has the following constructor:

DefaultJacobianDiffusion(int i = 0, int j = 0, std::string area = HERMES_ANY,
                         HermesFunction* coeff = HERMES_ONE,
                         SymFlag sym = HERMES_NONSYM, GeomType gt = HERMES_PLANAR);

The pair of indices ‘i’ and ‘j’ identifies a block in the Jacobian matrix (for systems of equations). For a single equation it is i = j = 0.

The parameter ‘area’ identifies the material marker of elements to which the weak form will be assigned. HERMES_ANY means to any material marker.

The parameter ‘coeff’ can be a constant, cubic spline, or a general nonlinear function of the solution u. HERMES_ONE means constant 1.0.

SymFlag is the symmetry flag. If SymFlag sym == HERMES_NONSYM, then Hermes evaluates the form at both symmetric positions r, s and s, r in the stiffness matrix. If sym == HERMES_SYM, only the integral at the position r, s is evaluated, and its value is copied to the symmetric position s, r. If sym == HERMES_ANTISYM, the value is copied with a minus sign.

The GeomType parameter tells Hermes whether the form is planar (HERMES_PLANAR), axisymmetrix with respect to the x-axis (HERMES_AXISYM_X), or axisymmetrix with respect to the y-axis (HERMES_AXISYM_Y).

The form can be linked to multiple material markers:

DefaultJacobianDiffusion(int i, int j, Hermes::vector<std::string> areas,
                         HermesFunction* coeff = HERMES_ONE,
                         SymFlag sym = HERMES_NONSYM, GeomType gt = HERMES_PLANAR);

Here, Hermes::vector is just a std::vector equipped with additional constructors for comfort. Sample usage:

Hermes::vector<std::string> areas("marker_1", "marker_2", "marker_3");

Default residual for the diffusion operator

Similarly, the line

add_vector_form(new DefaultResidualDiffusion(0, marker_al, lambda_al));

adds to the residual weak form the integral

\int_{\Omega_{al}} \lambda_{al} \nabla u \cdot \nabla v \, \mbox{d}x \mbox{d}y,

where u is the approximate solution and v a test function.

It has the following constructor:

DefaultResidualDiffusion(int i = 0, std::string area = HERMES_ANY,
                         HermesFunction* coeff = HERMES_ONE,
                         GeomType gt = HERMES_PLANAR);

The index ‘i’ identifies a block in the residual vector (for systems of equations). For a single equation it is i = 0. Again the form can be assigned to multiple material markers:

DefaultResidualDiffusion(int i, Hermes::vector<std::string> areas,
                         HermesFunction* coeff = HERMES_ONE,
                         GeomType gt = HERMES_PLANAR);

Default volumetric vector form

The last default weak form used in the CustomWeakFormPoisson class above is

add_vector_form(new DefaultVectorFormVol(0, HERMES_ANY, c));

It adds to the residual weak form the integral

\int_{\Omega} c v \, \mbox{d}x \mbox{d}y

and thus it completes (4). The constructors of this class are:

DefaultVectorFormVol(int i = 0, std::string area = HERMES_ANY,
                     HermesFunction* coeff = HERMES_ONE,
                     GeomType gt = HERMES_PLANAR);

and

DefaultVectorFormVol(int i, Hermes::vector<std::string> areas,
                     HermesFunction* coeff = HERMES_ONE,
                     GeomType gt = HERMES_PLANAR);

Loading the mesh

The main.cpp file typically begins with loading the mesh:

// Load the mesh.
Mesh mesh;
H2DReader mloader;
mloader.load("domain.mesh", &mesh);

Performing initial mesh refinements

A number of initial refinement operations can be done as explained above. In this case we just perform optional uniform mesh refinements:

// Perform initial mesh refinements (optional).
for (int i=0; i < INIT_REF_NUM; i++) mesh.refine_all_elements();

Initializing the weak formulation

Next, an instance of the corresponding weak form class is created:

// Initialize the weak formulation.
CustomWeakFormPoisson wf("Aluminum", new HermesFunction(LAMBDA_AL), "Copper",
                         new HermesFunction(LAMBDA_CU), new HermesFunction(-VOLUME_HEAT_SRC));

Setting constant Dirichlet boundary conditions

Constant Dirichlet boundary conditions are assigned to the boundary markers “Bottom”, “Inner”, “Outer”, and “Left” as follows:

// Initialize essential boundary conditions.
DefaultEssentialBCConst bc_essential(Hermes::vector<std::string>("Bottom", "Inner", "Outer", "Left"), FIXED_BDY_TEMP);
EssentialBCs bcs(&bc_essential);

Do not worry about the complicated-looking Hermes::vector, this is just std::vector enhanced with a few extra constructors. It is used to avoid using variable-length arrays.

The treatment of nonzero Dirichlet and other boundary conditions will be explained in more detail, and illustrated on examples, in the following examples. For the moment, let’s proceed to the finite element space.

Initializing finite element space

As a next step, we initialize the FE space in the same way as in the previous tutorial example 02-space:

// Create an H1 space with default shapeset.
H1Space space(&mesh, &bcs, P_INIT);
int ndof = space.get_num_dofs();
info("ndof = %d", ndof);

Here P_INIT is a uniform polynomial degree of mesh elements (an integer number between 1 and 10).

Initializing discrete problem

The weak formulation and finite element space(s) constitute a finite element problem. To define it, one needs to create an instance of the DiscreteProblem class:

// Initialize the FE problem.
DiscreteProblem dp(&wf, &space);

Initializing matrix solver

Next one needs to choose a matrix solver:

MatrixSolverType matrix_solver = SOLVER_UMFPACK;

Besides UMFPACK, one can use SOLVER_AMESOS, SOLVER_MUMPS, SOLVER_PETSC, and SOLVER_SUPERLU (and matrix-free SOLVER_NOX for nonlinear problems - to be discussed later).

After that one needs to create instances of a matrix, vector, and matrix solver as follows:

// Set up the solver, matrix, and rhs according to the solver selection.
SparseMatrix* matrix = create_matrix(matrix_solver);
Vector* rhs = create_vector(matrix_solver);
Solver* solver = create_linear_solver(matrix_solver, matrix, rhs);

Solving the discrete problem

Before solving the discrete problem, one has to create a coefficient vector:

// Initial coefficient vector for the Newton's method.
scalar* coeff_vec = new scalar[ndof];
memset(coeff_vec, 0, ndof*sizeof(scalar));

The discrete problem is solved via the Newton’s method:

// Perform Newton's iteration.
if (!hermes2d.solve_newton(coeff_vec, &dp, solver, matrix, rhs)) error("Newton's iteration failed.");

This function comes with a number of optional parameters, see the file hermes2d/src/h2d_common.h for more details.

Translating the coefficient vector into a solution

The coefficient vector can be converted into a piecewise-polynomial Solution via the function Solution::vector_to_solution():

// Translate the resulting coefficient vector into a Solution.
Solution sln;
Solution::vector_to_solution(coeff_vec, &space, &sln);

Saving solution in VTK format

The solution can be saved in the VTK format to be visualized, for example, using Paraview. To do this, one uses the Linearizer class that has the ability to approximate adaptively a higher-order polynomial solution using linear triangles:

// Output solution in VTK format.
Linearizer lin;
bool mode_3D = true;
lin.save_solution_vtk(&sln, "sln.vtk", "Temperature", mode_3D);
info("Solution in VTK format saved to file %s.", "sln.vtk");

The function save_solution_vtk() can be found in hermes2d/src/linearizer/ and its complete header is:

// Saves a MeshFunction (Solution, Filter) in VTK format.
virtual void save_solution_vtk(MeshFunction* meshfn, const char* file_name, const char* quantity_name,
                               bool mode_3D = true, int item = H2D_FN_VAL_0,
                               double eps = HERMES_EPS_NORMAL, double max_abs = -1.0,
                               MeshFunction* xdisp = NULL, MeshFunction* ydisp = NULL,
                               double dmult = 1.0);

Only the first three arguments are mandatory, the remaining ones are optional. Their meaning is as follows:

  • mode_3D ... select either 2D or 3D rendering (default is 3D).
  • item: H2D_FN_VAL_0 ... show function values, H2D_FN_DX_0 ... show x-derivative, H2D_FN_DY_0 ... show y-derivative, H2D_FN_DXX_0 ... show xx-derivative, H2D_FN_DXY_0 ... show xy-derivative, H2D_FN_DYY_0 ... show yy-derivative,
  • eps: HERMES_EPS_LOW ... low resolution (small output file), HERMES_EPS_NORMAL ... normal resolution (medium output file), HERMES_EPS_HIGH ... high resolution (large output file), HERMES_EPS_VERYHIGH ... high resolution (very large output file).
  • max_abs: technical parameter, see file src/linearizer/linear.h.
  • xdisp, ydisp, dmult: Can be used to deform the domain. Typical applications are elasticity, plasticity, etc.

The following figure shows the corresponding Paraview visualization:

Solution of the Poisson equation.

Visualizing the solution using OpenGL

The solution can also be visualized via the ScalarView class:

// Visualize the solution.
ScalarView view("Solution", new WinGeom(0, 0, 440, 350));
view.show(&sln);
View::wait();

Hermes’ built-in OpenGL visualization looks as follows:

Solution of the Poisson equation.

Visualization quality

The method show() has an optional second parameter – the visualization accuracy. It can have the values HERMES_EPS_LOW, HERMES_EPS_NORMAL (default), HERMES_EPS_HIGH and HERMES_EPS_VERYHIGH. This parameter influences the number of linear triangles that Hermes uses to approximate higher-order polynomial solutions with linear triangles for OpenGL. In fact, the EPS value is a stopping criterion for automatic adaptivity that Hermes uses to keep the number of the linear triangles as low as possible.

IMPORTANT: If you notice in the image white points or even discontinuities where the approximation is continuous, try to move from HERMES_EPS_NORMAL to HERMES_EPS_HIGH. If the interval of solution values is very small compared to the solution magnitude, such as if the solution values lie in the interval (50, 50.5), then you may need HERMES_EPS_VERYHIGH.

Before pressing ‘s’ to save the image, make sure to press ‘h’ to render high-quality image.

Visualization of derivatives

The method show() also has an optional third parameter to indicate whether function values or partial derivatives should be displayed. For example, HERMES_FN_VAL_0 stands for the function value of solution component 0 (first solution component which in this case is the VonMises stress). HERMES_FN_VAL_1 would mean the function value of the second solution component (relevant for vector-valued Hcurl or Hdiv elements only), HERMES_FN_DX_0 means the x-derivative of the first solution component, etc.