{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Implementation\n",
"\n",
"Author: Jørgen S. Dokken\n",
"\n",
"## Test problem\n",
"To solve a test problem, we need to choose the right hand side $f$, the coefficient $q(u)$, and the boundary $u_D$. Previously, we have worked with manufactured solutions that can be reproduced without approximation errors. This is more difficult in nonlinear problems, and the algebra is more tedious. However, we will utilize the UFL differentiation capabilities to obtain a manufactured solution.\n",
"\n",
"For this problem, we will choose $q(u) = 1 + u^2$ and define a two dimensional manufactured solution that is linear in $x$ and $y$:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import ufl\n",
"import numpy\n",
"\n",
"from mpi4py import MPI\n",
"from petsc4py import PETSc\n",
"\n",
"from dolfinx import mesh, fem, io, nls, log\n",
"from dolfinx.fem.petsc import NonlinearProblem\n",
"from dolfinx.nls.petsc import NewtonSolver\n",
"\n",
"def q(u):\n",
" return 1 + u**2\n",
"\n",
"\n",
"domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)\n",
"x = ufl.SpatialCoordinate(domain)\n",
"u_ufl = 1 + x[0] + 2 * x[1]\n",
"f = - ufl.div(q(u_ufl) * ufl.grad(u_ufl))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that since `x` is a 2D vector, the first component (index 0) represents $x$, while the second component (index 1) represents $y$. The resulting function `f` can be directly used in variational formulations in DOLFINx.\n",
"\n",
"As we now have defined our source term and an exact solution, we can create the appropriate function space and boundary conditions.\n",
"Note that as we have already defined the exact solution, we only have to convert it to a Python function that can be evaluated in the interpolation function. We do this by employing the Python `eval` and `lambda`-functions."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"V = fem.functionspace(domain, (\"Lagrange\", 1))\n",
"def u_exact(x): return eval(str(u_ufl))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"u_D = fem.Function(V)\n",
"u_D.interpolate(u_exact)\n",
"fdim = domain.topology.dim - 1\n",
"boundary_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: numpy.full(x.shape[1], True, dtype=bool))\n",
"bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are now ready to define the variational formulation. Note that as the problem is nonlinear, we have to replace the `TrialFunction` with a `Function`, which serves as the unknown of our problem."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"uh = fem.Function(V)\n",
"v = ufl.TestFunction(V)\n",
"F = q(uh) * ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - f * v * ufl.dx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Newton's method\n",
"The next step is to define the non-linear problem. As it is non-linear we will use [Newtons method](https://en.wikipedia.org/wiki/Newton%27s_method).\n",
"For details about how to implement a Newton solver, see [Custom Newton solvers](../chapter4/newton-solver.ipynb).\n",
"Newton's method requires methods for evaluating the residual `F` (including application of boundary conditions), as well as a method for computing the Jacobian matrix. DOLFINx provides the function `NonlinearProblem` that implements these methods. In addition to the boundary conditions, you can supply the variational form for the Jacobian (computed if not supplied), and form and jit parameters, see the [JIT parameters section](../chapter4/compiler_parameters.ipynb)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"problem = NonlinearProblem(F, uh, bcs=[bc])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we use the DOLFINx Newton solver. We can set the convergence criteria for the solver by changing the absolute tolerance (`atol`), relative tolerance (`rtol`) or the convergence criterion (`residual` or `incremental`)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"solver = NewtonSolver(MPI.COMM_WORLD, problem)\n",
"solver.convergence_criterion = \"incremental\"\n",
"solver.rtol = 1e-6\n",
"solver.report = True"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can modify the linear solver in each Newton iteration by accessing the underlying `PETSc` object."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"ksp = solver.krylov_solver\n",
"opts = PETSc.Options()\n",
"option_prefix = ksp.getOptionsPrefix()\n",
"opts[f\"{option_prefix}ksp_type\"] = \"cg\"\n",
"opts[f\"{option_prefix}pc_type\"] = \"gamg\"\n",
"opts[f\"{option_prefix}pc_factor_mat_solver_type\"] = \"mumps\"\n",
"ksp.setFromOptions()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are now ready to solve the non-linear problem. We assert that the solver has converged and print the number of iterations."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of interations: 8\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"2023-10-28 17:59:12.824 ( 0.944s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.827 ( 0.947s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.002726, 0.000000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.827 ( 0.947s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.828 ( 0.948s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000803, 0.000000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.828 ( 0.948s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 2: r (abs) = 20.379 (tol = 1e-10) r (rel) = 0.922522(tol = 1e-06)\n",
"2023-10-28 17:59:12.828 ( 0.948s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.828 ( 0.948s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000200, 0.010000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.829 ( 0.949s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 3: r (abs) = 6.95259 (tol = 1e-10) r (rel) = 0.314732(tol = 1e-06)\n",
"2023-10-28 17:59:12.829 ( 0.949s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.830 ( 0.950s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000229, 0.000000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.830 ( 0.950s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 4: r (abs) = 2.93565 (tol = 1e-10) r (rel) = 0.132892(tol = 1e-06)\n",
"2023-10-28 17:59:12.830 ( 0.950s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.830 ( 0.950s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000203, 0.000000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.830 ( 0.950s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 5: r (abs) = 0.700567 (tol = 1e-10) r (rel) = 0.0317135(tol = 1e-06)\n",
"2023-10-28 17:59:12.830 ( 0.950s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.831 ( 0.951s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000232, 0.000000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.831 ( 0.951s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 6: r (abs) = 0.0490758 (tol = 1e-10) r (rel) = 0.00222158(tol = 1e-06)\n",
"2023-10-28 17:59:12.831 ( 0.951s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.831 ( 0.951s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000179, 0.000000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.831 ( 0.951s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 7: r (abs) = 0.000299251 (tol = 1e-10) r (rel) = 1.35466e-05(tol = 1e-06)\n",
"2023-10-28 17:59:12.831 ( 0.951s) [main ] petsc.cpp:698 INFO| PETSc Krylov solver starting to solve system.\n",
"2023-10-28 17:59:12.832 ( 0.952s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000265, 0.000000, 0.000000 (PETSc Krylov solver)\n",
"2023-10-28 17:59:12.832 ( 0.952s) [main ] NewtonSolver.cpp:36 INFO| Newton iteration 8: r (abs) = 1.62359e-08 (tol = 1e-10) r (rel) = 7.34971e-10(tol = 1e-06)\n",
"2023-10-28 17:59:12.832 ( 0.952s) [main ] NewtonSolver.cpp:255 INFO| Newton solver finished in 8 iterations and 109 linear solver iterations.\n"
]
}
],
"source": [
"log.set_log_level(log.LogLevel.INFO)\n",
"n, converged = solver.solve(uh)\n",
"assert (converged)\n",
"print(f\"Number of interations: {n:d}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We observe that the solver converges after $8$ iterations.\n",
"If we think of the problem in terms of finite differences on a uniform mesh, $\\mathcal{P}_1$ elements mimic standard second-order finite differences, which compute the derivative of a linear or quadratic funtion exactly. Here $\\nabla u$ is a constant vector, which is multiplied by $1+u^2$, giving a second order polynomial in $x$ and $y$, which the finite difference operator would compute exactly. We can therefore, even with $\\mathcal{P}_1$ elements, expect the manufactured solution to be reproduced by the numerical method. However, if we had chosen a nonlinearity, such as $1+u^4$, this would not be the case, and we would need to verify convergence rates."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"L2-error: 5.94e-15\n",
"Error_max: 1.60e-14\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"2023-10-28 17:59:12.848 ( 0.968s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000023, 0.000000, 0.000000 (Init dofmap from element dofmap)\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000006, 0.000000, 0.000000 (Compute dof reordering map)\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:166 INFO| Computing communication graph edges (using NBX algorithm). Number of input edges: 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] MPI.cpp:237 INFO| Finished graph edge discovery using NBX algorithm. Number of discovered edges 0\n",
"2023-10-28 17:59:12.848 ( 0.968s) [main ] TimeLogger.cpp:28 INFO| Elapsed wall, usr, sys time: 0.000236, 0.000000, 0.000000 (Build dofmap data)\n"
]
}
],
"source": [
"# Compute L2 error and error at nodes\n",
"V_ex = fem.functionspace(domain, (\"Lagrange\", 2))\n",
"u_ex = fem.Function(V_ex)\n",
"u_ex.interpolate(u_exact)\n",
"error_local = fem.assemble_scalar(fem.form((uh - u_ex)**2 * ufl.dx))\n",
"error_L2 = numpy.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))\n",
"if domain.comm.rank == 0:\n",
" print(f\"L2-error: {error_L2:.2e}\")\n",
"\n",
"# Compute values at mesh vertices\n",
"error_max = domain.comm.allreduce(numpy.max(numpy.abs(uh.x.array - u_D.x.array)), op=MPI.MAX)\n",
"if domain.comm.rank == 0:\n",
" print(f\"Error_max: {error_max:.2e}\")"
]
}
],
"metadata": {
"jupytext": {
"formats": "ipynb,py:light"
},
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}