 v0.9.0
A nonlinear Poisson equation

In this tutorial, we show how to solve nonlinear Poisson equation. Note that structure of the main program remains almost unchanged compared to linear analysis as shown in the first Poisson tutorial, i.e. Solving the Poisson equation.

PDE problem

As a model problem for the solution of nonlinear PDEs, we take the following nonlinear Poisson equation:

$-\nabla \cdot (A(u)\nabla u) = f$

in domain $$\Omega$$ with $$u=\overline{u}$$ on boundary $$\partial\Omega$$. The source of nonlinearity is the coefficient $$A(u)$$ which is dependent on the solution $$u$$, here we restrict ourself to the case

$A(u) = 1+u^2$

Note, code implementation is general, and you can look at the alternative cases, at the beginning of analytical_nonlinear_poisson.cpp example you can find two classes. Those functions can be modified to investigate alternative solutions

struct FunA {
double operator()(const double u) { return 1+u*u; }
};
struct DiffFunA {
double operator()(const double u) { return 2*u; }
};

Note that second class has to be the exact derivative of $$A(u)$$.

Linearization

With that at hand we can express equation in the weak form, where residuals are

$\mathbf{r}= \left[ \begin{array}{c} r_F\\ r_g \end{array} \right]= \left\{ \begin{array}{l} (\nabla v,(1+u^2)\nabla u)_\Omega-(v,f)_\Omega+(v,\lambda)_{\partial\Omega} = 0\\ (\tau,u-\overline{u})_{\partial\Omega} = 0 \end{array} \right.,\quad \forall v,\tau,\; u,v \in H^1(\Omega),\;\tau,\lambda \in H^1(\partial\Omega).$

where $$\overline{u}$$ is given function on boundary, $$f$$ is a source term, $$u,\lambda$$ are unknown functions and $$v,\tau$$ are test functions. This notation, suitable for problems with many terms in the variational formulations, requires some explanation. First, we use the shorthand notation

$(v,u)_\Omega = \int_\Omega vu \textrm{d}\Omega,\quad (\lambda,u)_{\partial\Omega} = \int_{\partial \Omega} \lambda u \textrm{d}\partial\Omega$

With that at hand, we expand residuals in Taylor series and truncate after first order term, as a result we have

$\mathbf{r}^i + \frac{\partial \mathbf{r}^i}{\partial \{ u, \lambda \}} \left\{ \begin{array}{c} \delta u^{i+1}\\ \delta \lambda^{i+1} \end{array} \right\} = \mathbf{0}$

where $$(\cdot)^i$$ is quantity at Newton iteration $$i$$, $$u^{i+1} = u^i + \delta u^{i+1}$$ and $$\lambda^{i+1} = \lambda^i + \delta \lambda^{i+1}$$. The linearized system of equations takes form

$\left\{ \begin{array}{l} r^i_F + (\nabla v,(1+(u^i)^2)\nabla \delta u)_\Omega+(\nabla v,(2u^i) \nabla u^i \delta u^{i+1})_\Omega+(v,\delta\lambda^{i+1}) = 0 \\ r^i_g+ (\tau,\delta u^{i+1})_{\partial\Omega}=0 \end{array} \right.$

and applying discretization

$u^h = \sum_i^{n_u} \phi^i(\mathbf{x}) U_i = \boldsymbol\phi\cdot\mathbf{U} , \quad \mathbf{x} \in \Omega,\\ \lambda^h = \sum_i^{n_\lambda} \psi^i(\mathbf{x}) L_i = \boldsymbol\psi\cdot\mathbf{L} , \quad \mathbf{x} \in \partial\Omega$

we get

$\left[ \begin{array}{c} \mathbf{r}_F\\ \mathbf{r}_g \end{array} \right]+ \left[ \begin{array}{cc} \mathbf{K}_t & \mathbf{C}^T \\ \mathbf{C} & \mathbf{0} \end{array} \right] \left\{ \begin{array}{c} \delta\mathbf{U}^{i+1}\\ \delta\mathbf{L}^{i+1} \end{array} \right\} = \left[ \begin{array}{c} \mathbf{0}\\ \mathbf{0} \end{array} \right]$

where

$\mathbf{K}_t = \int_\Omega (\nabla \boldsymbol\phi)^\textrm{T} (1+(u^i)^2)\boldsymbol\phi + (\nabla \boldsymbol\phi)^\textrm{T} (2(u^i))\boldsymbol u^i \boldsymbol\phi \textrm{d}\Omega,\\ \mathbf{C} = \int_{\partial\Omega} \boldsymbol\psi^\textrm{T} \boldsymbol\phi \textrm{d}\partial\Omega,\\ \mathbf{r}_F = \mathbf{r}_{F_\Omega} + \mathbf{r}_{F_{\partial\Omega}} = \int_\Omega (\nabla \boldsymbol\phi)^\textrm{T} (1+(u^i)^2)\boldsymbol u^i - \boldsymbol\phi^\textrm{T}f \textrm{d}\Omega+ \int_{\partial\Omega} \boldsymbol\phi^\textrm{T} \lambda^i \textrm{d}\partial\Omega\\ \mathbf{r}_g = \int_{\partial\Omega} \boldsymbol\psi^\textrm{T} (u^i-\overline{u}) \textrm{d}\partial\Omega$

Analytical solution

To validate correctness of implementation, we first presume exact solution to be polynomial, such that it can be precisely represented by finite element approximation base, see function

struct ExactFunction {
double operator()(const double x,const double y,const double z) const {
return 1+x+y+pow(z,3);
}
};

where the gradient of that function is implemented here

FTensor::Tensor1<double,3> operator()(const double x,const double y,const double z) const {
}
};

and source term

double operator()(const double x,const double y,const double z) const {
return 0.4e1 + (double) (4 * x) + (double) (4 * y) + 0.4e1 * pow(z, 0.3e1) +
0.3e1 * (0.6e1 * z * z + 0.6e1 * (double) x * z * z +
0.6e1 * (double) y * z * z + 0.6e1 * pow(z, 0.5e1)) * z * z +
0.6e1 * (0.2e1 + (double) (2 * x) + (double) (2 * y) +
0.2e1 * pow(z, 0.3e1) + (double) (x * x) + (double) (2 * x * y) +
0.2e1 * (double) x * pow(z, 0.3e1) +
(double) (y * y) + 0.2e1 * (double) y * pow(z, 0.3e1) + pow(z, 0.6e1)) * z;
}
};

Note that source term has a complex form and is obtained from the following formula

$f = \nabla \cdot \left( A(u)\nabla u \right)$

For 3rd polynomial order, as above, results can be obtained after tedious calculations, to avoid that you can use symbolic derivations from Mathematica, Maple, Matlab or free and open alternative Maxima.

Code dissection

Implementation of operators

All operators from this example are implemented in PoissonOperators.cpp. In particular integration is implemented in the following methods

Finite element instances and setup of finite element operators, following Figure 1, is implanted in PoissonExample::CreateFiniteElements::createFEToAssembleMatrixAndVectorForNonlinearProblem, in particular we find there the following code

boost::shared_ptr<VectorDouble> values_at_integation_ptr = boost::make_shared<VectorDouble>();
boost::shared_ptr<MatrixDouble> grad_at_integation_ptr = boost::make_shared<MatrixDouble>();
domain_lhs_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues("U",values_at_integation_ptr));

You can note that instances of vector and matrix are created to store values of function $$u^h$$ and function gradients at integration points, respectively. Next operators are added to evaluate function and function gradients at integration points. Finally, operator to calculate tangent matrix is added. Operator constructor takes arguments of pointers of previously created instances of vector and matrix. The same procedure is applied to evaluate residuals by integration over domain and boundary. Figure 1. Operators of volume finite element

Nonlinear Poisson code

/**
* \file analytical_nonlinear_poisson.cpp
* \ingroup mofem_simple_interface
* \example analytical_nonlinear_poisson.cpp
*
* For more information and detailed explain of this
* example see \ref poisson_tut4
*
*
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
static char help[] = "...\n\n";
/**
* \brief Function
*
* This is prescribed exact function. If this function is given by polynomial
* order of "p" and order of approximation is "p" or higher, solution of
* finite element method is exact (with machine precision).
*
* \f[
* u = 1+x+2y+3z
* \f]
*
*/
struct ExactFunction {
double operator()(const double x, const double y, const double z) const {
return 1 + x + y + pow(z, 3);
}
};
/**
* \brief Exact gradient
*/
FTensor::Tensor1<double, 3> operator()(const double x, const double y,
const double z) const {
grad(2) = 3 * z * z;
}
};
/**
* \brief Laplacian of function.
*
*/
double operator()(const double x, const double y, const double z) const {
return 0.4e1 + (double)(4 * x) + (double)(4 * y) + 0.4e1 * pow(z, 0.3e1) +
0.3e1 *
(0.6e1 * z * z + 0.6e1 * (double)x * z * z +
0.6e1 * (double)y * z * z + 0.6e1 * pow(z, 0.5e1)) *
z * z +
0.6e1 *
(0.2e1 + (double)(2 * x) + (double)(2 * y) +
0.2e1 * pow(z, 0.3e1) + (double)(x * x) + (double)(2 * x * y) +
0.2e1 * (double)x * pow(z, 0.3e1) + (double)(y * y) +
0.2e1 * (double)y * pow(z, 0.3e1) + pow(z, 0.6e1)) *
z;
}
};
struct FunA {
double operator()(const double u) { return 1 + u * u; }
};
struct DiffFunA {
double operator()(const double u) { return 2 * u; }
};
int operator()(int, int, int p) const { return 2 * (p + 1); }
};
int main(int argc, char *argv[]) {
// Initialize PETSc
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
// Create MoAB database
moab::Core moab_core; // create database
moab::Interface &moab = moab_core; // create interface to database
try {
// Get command line options
int order = 3; // default approximation order
PetscBool flg_test = PETSC_FALSE; // true check if error is numerical error
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Poisson's problem options",
"none");
// Set approximation order
CHKERR PetscOptionsInt("-order", "approximation order", "", order, &order,
PETSC_NULL);
// Set testing (used by CTest)
CHKERR PetscOptionsBool("-test", "if true is ctest", "", flg_test,
&flg_test, PETSC_NULL);
ierr = PetscOptionsEnd();
// Create MoFEM database and link it to MoAB
MoFEM::Core mofem_core(moab); // create database
MoFEM::Interface &m_field = mofem_core; // create interface to database
// Register DM Manager
CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
// Create vector to store approximation global error
Vec global_error;
// First we crate elements, implementation of elements is problem
// independent, we do not know yet what fields are present in the problem,
// or even we do not decided yet what approximation base or spaces we are
// going to use. Implementation of element is free from those constrains and
// can be used in different context.
// Elements used by KSP & DM to assemble system of equations
boost::shared_ptr<ForcesAndSourcesCore>
domain_lhs_fe; ///< Volume element for the matrix
boost::shared_ptr<ForcesAndSourcesCore>
boundary_lhs_fe; ///< Boundary element for the matrix
boost::shared_ptr<ForcesAndSourcesCore>
domain_rhs_fe; ///< Volume element to assemble vector
boost::shared_ptr<ForcesAndSourcesCore>
boundary_rhs_fe; ///< Volume element to assemble vector
boost::shared_ptr<ForcesAndSourcesCore>
domain_error; ///< Volume element evaluate error
boost::shared_ptr<ForcesAndSourcesCore>
post_proc_volume; ///< Volume element to Post-process results
boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element do nothing
{
// Add problem specific operators the generic finite elements to calculate
// matrices and vectors.
domain_lhs_fe, boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe,
// Add problem specific operators the generic finite elements to calculate
// error on elements and global error in H1 norm
global_error, domain_error);
// Post-process results
.creatFEToPostProcessResults(post_proc_volume);
}
// Get simple interface is simplified version enabling quick and
// easy construction of problem.
Simple *simple_interface;
// Query interface and get pointer to Simple interface
CHKERR m_field.getInterface(simple_interface);
// Build problem with simple interface
{
// Get options for simple interface from command line
CHKERR simple_interface->getOptions();
// Load mesh file to database
// Add field on domain and boundary. Field is declared by space and base
// and rank. space can be H1. Hcurl, Hdiv and L2, base can be
// AINSWORTH_LEGENDRE_BASE, DEMKOWICZ_JACOBI_BASE and more, where rank is
// number of coefficients for dof.
//
// Simple interface assumes that there are four types of field; 1) defined
// on body domain, 2) fields defined on body boundary, 3) skeleton field
// defined on finite element skeleton. Finally data field is defined on
// body domain. Data field is not solved but used for post-process or to
// keep material parameters, etc. Here we using it to calculate
// approximation error on elements.
// Add domain filed "U" in space H1 and Legendre base, Ainsworth recipe is
// used to construct base functions.
CHKERR simple_interface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE,
1);
// Add Lagrange multiplier field on body boundary
// Add error (data) field, we need only L2 norm. Later order is set to 0,
// so this is piecewise discontinuous constant approx., i.e. 1 DOF for
// element. You can use more DOFs and collate moments of error to drive
// anisotropic h/p-adaptivity, however this is beyond this example.
// Set fields order domain and boundary fields.
CHKERR simple_interface->setFieldOrder("U",
order); // to approximate function
CHKERR simple_interface->setFieldOrder("L",
order); // to Lagrange multipliers
CHKERR simple_interface->setFieldOrder(
"ERROR", 0); // approximation order for error
// Setup problem. At that point database is constructed, i.e. fields,
// finite elements and problem data structures with local and global
// indexing.
CHKERR simple_interface->setUp();
}
// Get access to PETSC-MoFEM DM manager.
// or more derails see
// <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DM/index.html>
// Form that point internal MoFEM data structures are linked with PETSc
// interface. If DM functions contains string MoFEM is is MoFEM specific DM
// interface function, otherwise DM function part of standard PETSc
// interface.
DM dm;
// Get dm
CHKERR simple_interface->getDM(&dm);
// Set KSP context for DM. At that point only elements are added to DM
// operators. Calculations of matrices and vectors is executed by KSP
// solver. This part of the code makes connection between implementation of
// finite elements and data operators with finite element declarations in DM
// manager. From that point DM takes responsibility for executing elements,
// calculations of matrices and vectors, and solution of the problem.
{
// Set operators for SNES solver
CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getDomainFEName(),
domain_lhs_fe, null, null);
CHKERR DMMoFEMSNESSetJacobian(dm, simple_interface->getBoundaryFEName(),
boundary_lhs_fe, null, null);
// Set calculation of the right hand side vector for KSP solver
CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getDomainFEName(),
domain_rhs_fe, null, null);
CHKERR DMMoFEMSNESSetFunction(dm, simple_interface->getBoundaryFEName(),
boundary_rhs_fe, null, null);
}
// Solve problem, only PETEc interface here.
{
// Create the right hand side vector and vector of unknowns
Vec F, D;
CHKERR DMCreateGlobalVector(dm, &F);
// Create unknown vector by creating duplicate copy of F vector. only
// structure is duplicated no values.
CHKERR VecDuplicate(F, &D);
// Create solver and link it to DM
SNES solver;
CHKERR SNESCreate(PETSC_COMM_WORLD, &solver);
CHKERR SNESSetFromOptions(solver);
CHKERR SNESSetDM(solver, dm);
// Set-up solver, is type of solver and pre-conditioners
CHKERR SNESSetUp(solver);
// At solution process, KSP solver using DM creates matrices, Calculate
// values of the left hand side and the right hand side vector. then
// solves system of equations. Results are stored in vector D.
CHKERR SNESSolve(solver, F, D);
// Scatter solution on the mesh. Stores unknown vector on field on the
// mesh.
CHKERR DMoFEMMeshToGlobalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
// Clean data. Solver and vector are not needed any more.
CHKERR SNESDestroy(&solver);
CHKERR VecDestroy(&D);
CHKERR VecDestroy(&F);
}
// Calculate error
{
// Loop over all elements in mesh, and run users operators on each
// element.
CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
domain_error);
global_error);
if (flg_test == PETSC_TRUE) {
}
}
{
// Loop over all elements in the mesh and for each execute
// post_proc_volume element and operators on it.
CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
post_proc_volume);
// Write results
CHKERR boost::static_pointer_cast<PostProcVolumeOnRefinedMesh>(
post_proc_volume)
->writeFile("out_vol.h5m");
}
// Destroy DM, no longer needed.
CHKERR DMDestroy(&dm);
// Destroy ghost vector
CHKERR VecDestroy(&global_error);
}
// finish work cleaning memory, getting statistics, etc.
return 0;
}

Following code is practically the same which we used to linear Poisson example. In nonlinear case, instead of KSP solver, we set-up SNES solver, first by adding finite elements instances to MoFEM SNES solver context

CHKERR DMMoFEMSNESSetJacobian(dm,simple_interface->getDomainFEName(),domain_lhs_fe,null,null);
CHKERR DMMoFEMSNESSetJacobian(dm,simple_interface->getBoundaryFEName(),boundary_lhs_fe,null,null);
CHKERR DMMoFEMSNESSetFunction(dm,simple_interface->getDomainFEName(),domain_rhs_fe,null,null);
CHKERR DMMoFEMSNESSetFunction(dm,simple_interface->getBoundaryFEName(),boundary_rhs_fe,null,null);

Once we have done this, we create solver instance and kick-start calculations

SNES solver;
CHKERR SNESCreate(PETSC_COMM_WORLD,&solver);
CHKERR SNESSetFromOptions(solver);
CHKERR SNESSetDM(solver,dm);
CHKERR SNESSetUp(solver);
CHKERR SNESSolve(solver,F,D);

Note that SNES solver controls DM, and implicitly query of DM to create necessary matrices and assembly of the tangent matrix and residual vectors. SNES solver can be configured from a command line, see for details.

Running code

In order to run the program, one should first go to the directory where the problem is located, compile the code and then run the executable file. Typically, this can be done as follows

cd mofem_install/um/build/basic_finite_elements/poisson
make -j2
mpirun -np 2 ./analytical_nonlinear_poisson -file_name cube_2part.h5m -order 3 \
-snes_monitor -snes_converged_reason -snes_type newtonls \
-snes_linesearch_type cp -snes_linesearch_monitor

where SNES options are explained here.

As a result of above command, we get

0 SNES Function norm 6.585598470492e-01
0 KSP Residual norm 6.585598470492e-01
1 KSP Residual norm 5.914941314228e-15
Line search: lambdas = [1., 0.], ftys = [-3.52513, -12.2419]
Line search terminated: lambda = 0.595591, fnorms = 0.34441
1 SNES Function norm 3.444100959814e-01
0 KSP Residual norm 3.444100959814e-01
1 KSP Residual norm 4.038504742986e-15
Line search: lambdas = [1., 0.], ftys = [-0.443111, -2.89313]
Line search terminated: lambda = 0.81914, fnorms = 0.160021
2 SNES Function norm 1.600212778450e-01
0 KSP Residual norm 1.600212778450e-01
1 KSP Residual norm 1.468413895487e-15
Line search: lambdas = [1., 0.], ftys = [-0.00392301, -0.184194]
Line search terminated: lambda = 0.978238, fnorms = 0.0159661
3 SNES Function norm 1.596606319710e-02
0 KSP Residual norm 1.596606319710e-02
1 KSP Residual norm 4.135944054060e-17
Line search: lambdas = [1., 0.], ftys = [-2.21812e-07, -0.000140473]
Line search terminated: lambda = 0.998418, fnorms = 2.23868e-05
4 SNES Function norm 2.238684200528e-05
0 KSP Residual norm 2.238684200528e-05
1 KSP Residual norm 7.619868739966e-20
Line search: lambdas = [1., 0.], ftys = [-4.35284e-16, -4.97174e-10]
Line search terminated: lambda = 0.999999, fnorms = 4.34694e-11
5 SNES Function norm 4.346943842203e-11
Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 5
Approximation error 4.109e-10

Note fast convergence, after first two interactions, once iterative approximation approaches the exact solution. To obtain this solution, we are using Newton method with line-searcher, in this case, we use critical point secant line search, see details Line-search method has several alternative types and can be setup independently from solver, they can be found here.