You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
The purpose of this demo is to show how to use Firedrake's multigrid solver on a hierarchy of adaptively refined Netgen meshes.
5
-
We will first have a look at how to use the *AdaptiveMeshHierarchy* to construct the mesh hierarchy with Netgen meshes, then we will consider a solution to the Poisson problem on an L-shaped domain.
6
-
Finally, we will show how to use the *AdaptiveMeshHierarchy* and *AdaptiveTransferManager* as arguments to Firedrake solvers. The *AdaptiveMeshHierarchy* contains information of the mesh hierarchy and the parent child relations between the meshes.
7
-
The *AdaptiveTransferManager* deals with the transfer operator logic across any given levels in the hierarchy.
8
+
We will first have a look at how to use the :class:`.AdaptiveMeshHierarchy` to construct the mesh hierarchy with Netgen meshes, then we will consider a solution to the Poisson problem on an L-shaped domain.
9
+
Finally, we will show how to use the :class:`.AdaptiveMeshHierarchy` and :class:`.AdaptiveTransferManager` to construct a scalable solver. The :class:`.AdaptiveMeshHierarchy` contains information of the mesh hierarchy and the parent child relations between the meshes.
10
+
The :class:`.AdaptiveTransferManager` deals with the transfer operator logic across any given levels in the hierarchy.
8
11
We begin by importing the necessary libraries ::
9
12
10
13
from firedrake import *
11
14
from netgen.occ import *
12
-
from ngsPETSc import AdaptiveMeshHierarchy, AdaptiveTransferManager
15
+
import numpy
13
16
14
17
Constructing the Mesh Hierarchy
15
-
---------------------------
16
-
We first must construct the domain over which we will solve the problem. For a more comprehensive demo on how to use Open Cascade Technology (OCC) and Constructive Solid Geometry (CSG), see `Netgen integration in Firedrake <https://www.firedrakeproject.org/demos/netgen_mesh.py.html>`_ .
18
+
-------------------------------
19
+
We first must construct the domain over which we will solve the problem. For a more comprehensive demo on how to use Open Cascade Technology (OCC) and Constructive Solid Geometry (CSG),
20
+
see `Netgen integration in Firedrake <netgen_mesh.py>`_.
17
21
We begin with the L-shaped domain, which we build as the union of two rectangles: ::
It is important to convert the initial Netgen mesh into a Firedrake mesh before constructing the *AdaptiveMeshHierarchy*. To call the constructor to the hierarchy, we must pass the initial mesh. Our initial mesh looks like this:
31
+
It is important to convert the initial Netgen mesh into a Firedrake mesh before constructing the :class:`.AdaptiveMeshHierarchy`. To call the constructor to the hierarchy, we must pass the initial mesh. Our initial mesh looks like this:
27
32
28
33
.. figure:: initial_mesh.png
29
34
:align:center
30
35
:alt:Initial mesh.
31
36
32
-
We will also initialize the *AdaptiveTransferManager* here: ::
37
+
We will also initialize the :class:`.AdaptiveTransferManager` here: ::
33
38
34
-
amh = AdaptiveMeshHierarchy(mesh)
35
-
atm = AdaptiveTransferManager()
39
+
amh = AdaptiveMeshHierarchy(mesh)
40
+
atm = AdaptiveTransferManager()
36
41
37
42
Poisson Problem
38
-
-------------------------
43
+
---------------
39
44
Now we can define a simple Poisson problem
40
45
41
46
.. math::
42
47
43
48
- \nabla^2 u = f \text{ in } \Omega, \quad u = 0\text{ on } \partial\Omega
44
49
45
-
Our approach strongly follows the similar problem in this `lecture course <https://github.com/pefarrell/icerm2024>`_. We define the function ``solve_poisson``. The first lines correspond to finding a solution in the CG1 space. The variational problem is formulated with F, where f is the constant function equal to 1. Since we want Dirichlet boundary conditions, we construct the *DirichletBC* object and apply it to the entire boundary: ::
50
+
Our approach strongly follows the similar problem in this `lecture course <https://github.com/pefarrell/icerm2024>`_. We define the function ``solve_poisson``. The first lines correspond to finding a solution in the CG1 space. The variational problem is formulated with F, where f is the constant function equal to 1. Since we want Dirichlet boundary conditions, we construct the :class:`.DirichletBC` object and apply it to the entire boundary: ::
46
51
47
52
def solve_poisson(mesh, params):
48
53
V = FunctionSpace(mesh, "CG", 1)
49
-
uh = Function(V, name="Solution")
54
+
uh = Function(V, name="solution")
50
55
v = TestFunction(V)
51
-
bc = DirichletBC(V, Constant(0), "on_boundary")
56
+
bc = DirichletBC(V, 0, "on_boundary")
52
57
f = Constant(1)
53
58
F = inner(grad(uh), grad(v))*dx - inner(f, v)*dx
54
59
@@ -58,57 +63,47 @@ Our approach strongly follows the similar problem in this `lecture course <https
58
63
solver.solve()
59
64
return uh
60
65
61
-
Note the code after the construction of the *NonlinearVariationalProblem()*. To use the *AdaptiveMeshHierarchy* with the existing Firedrake solver, we have to set the *AdaptiveTransferManager* as the transfer manager of the multigrid solver.
62
-
For the parameters of the multigrid solver, we will be using patch relaxation, which we define with ::
66
+
Note the code after the construction of the :class:`.NonlinearVariationalProblem`. To use the :class:`.AdaptiveMeshHierarchy` with the existing Firedrake solver, we have to set the :class:`.AdaptiveTransferManager` as the transfer manager of the multigrid solver.
67
+
Since we are using linear CG elements, we will employ Jacobi as the multigrid relaxation, which we define with ::
For more information about patch relaxation, see `Using patch relaxation for multigrid <https://www.firedrakeproject.org/demos/poisson_mg_patches.py.html>`_. The initial solution is shown below.
70
+
"mat_type": "matfree",
71
+
"ksp_type": "cg",
72
+
"pc_type": "mg",
73
+
"mg_levels": {
74
+
"ksp_type": "chebyshev",
75
+
"ksp_max_it": 1,
76
+
"pc_type": "jacobi",
77
+
},
78
+
"mg_coarse": {
79
+
"mat_type": "aij",
80
+
"pc_type": "lu",
81
+
},
82
+
}
83
+
84
+
Alternatively for high-order CG elements, it is recommended to use patch relaxation
85
+
to achieve degree-independent multigrid convergence.
86
+
For more information
87
+
see :doc:`Using patch relaxation for multigrid <poisson_mg_patches.py>`.
88
+
The initial solution is shown below.
94
89
95
90
.. figure:: solution_l1.png
96
91
:align:center
97
92
:alt:Initial Solution from multigrid with initial mesh.
98
93
99
94
100
95
Adaptive Mesh Refinement
101
-
-------------------------
102
-
In this section we will discuss how to adaptively refine select elements and add the newly refined mesh into the *AdaptiveMeshHierarchy*.
96
+
------------------------
97
+
In this section we will discuss how to adaptively refine select elements and add the newly refined mesh into the :class:`.AdaptiveMeshHierarchy`.
103
98
For this problem, we will be using the Babuška-Rheinbolt a-posteriori estimate for an element:
104
99
105
100
.. math::
106
-
\eta_K^2 = h_K^2\int_K | f + \nabla^2 u_h |^2\mathrm{d}x + \frac{h_K}{2} \int_{\partial K \setminus\partial\Omega} \llbracket\nabla u_h \cdot n \rrbracket^2\mathrm{d}s,
101
+
\eta_K^2 = h_K^2\int_K \| f + \nabla^2 u_h \|^2\mathrm{d}x + \frac{h_K}{2} \int_{\partial K \setminus\partial\Omega} ⟦\nabla u_h \cdot n ⟧^2\mathrm{d}s,
107
102
108
-
where :math:`K` is the element, :math:`h_K` is the diameter of the element, :math:`n` is the normal, and :math:`\llbracket\cdot\rrbracket` is the jump operator. The a-posteriori estimator is computed using the solution at the current level :math:`h`. Integrating over the domain and using the fact that the components of the estimator are piecewise constant on each cell, we can transform the above estimator into the variational problem
103
+
where :math:`K` is the element, :math:`h_K` is the diameter of the element, :math:`n` is the normal, and :math:`⟦\cdot⟧` is the jump operator. The a-posteriori estimator is computed using the solution at the current level :math:`h`. Integrating over the domain and using the fact that the components of the estimator are piecewise constant on each cell, we can transform the above estimator into the variational problem
109
104
110
105
.. math::
111
-
\int_\Omega\eta_K^2 w \mathrm{d}x = \int_\Omega\sum_K h_K^2\int_K (f + \text{div} (\text{grad} u_h) )^2\mathrm{d}x w \mathrm{d}x + \int_\Omega\sum_K \frac{h_K}{2} \int_{\partial K \setminus\partial\Omega} \llbracket\nabla u_h \cdot n \rrbracket^2\mathrm{d}s w \mathrm{d}x
106
+
\int_\Omega\eta_K^2 w \mathrm{d}x = \int_\Omega\sum_K h_K^2\int_K (f + \text{div} (\text{grad} u_h) )^2\mathrm{d}x w \mathrm{d}x + \int_\Omega\sum_K \frac{h_K}{2} \int_{\partial K \setminus\partial\Omega} ⟦\nabla u_h \cdot n ⟧^2\mathrm{d}s w \mathrm{d}x
112
107
113
108
Our approach will be to compute the estimator over all elements and selectively choose to refine only those that contribute most to the error. To compute the error estimator, we use the function below to solve the variational formulation of the error estimator. Since our estimator is a constant per element, we use a DG0 function space. ::
114
109
@@ -117,54 +112,73 @@ Our approach will be to compute the estimator over all elements and selectively
117
112
eta_sq = Function(W)
118
113
w = TestFunction(W)
119
114
f = Constant(1)
120
-
h = CellDiameter(mesh) # symbols for mesh quantities
eta = Function(W).interpolate(sqrt(eta_sq)) # compute eta from eta^2
134
131
135
132
with eta.dat.vec_ro as eta_: # compute estimate for error in energy norm
136
-
error_est = sqrt(eta_.dot(eta_))
137
-
return (eta, error_est)
133
+
error_est = eta_.norm()
134
+
return eta, error_est
138
135
139
-
The next step is to choose which elements to refine. For this we use Dörfler marking [Dörfler1996]:
136
+
The next step is to choose which elements to refine. For this we use Dörfler marking :cite:`Dorfler1996`:
140
137
141
138
.. math::
142
139
\eta_K \geq\theta\text{max}_L \eta_L
143
140
144
-
The logic is to select an element :math:`K` to refine if the estimator is greater than some factor :math:`\theta` of the maximum error estimate of the mesh, where :math:`\theta` ranges from 0 to 1. In our code we choose :math:`theta=0.5`.
141
+
The logic is to select an element :math:`K` to refine if the estimator is greater than some factor :math:`\theta` of the maximum error estimate of the mesh, where :math:`\theta` ranges from 0 to 1. In our code we choose :math:`\theta=0.5`.
145
142
With these helper functions complete, we can solve the system iteratively. In the max_iterations is the number of total levels we want to perform multigrid on. We will solve for 15 levels. At every level :math:`l`, we first compute the solution using multigrid with patch relaxation up till level :math:`l`. We then use the current approximation of the solution to estimate the error across the mesh. Finally, we refine the mesh and repeat. ::
To perform Dörfler marking, refine the current mesh, and add the mesh to the *AdaptiveMeshHierarchy*, we us the ``amh.adapt(eta, theta)`` method. In this method the input is the recently computed error estimator :math:`eta` and the Dörfler marking parameter :math:`theta`. The method always performs this on the current fine mesh in the hierarchy. There is another method for adding a mesh to the hierarchy: ``amh.add_mesh(mesh)``. In this method, refinement on the mesh is performed externally by some custom procedure and the resulting mesh directly gets added to the hierarchy.
167
-
The meshes now refine according to the error estimator. The error estimators at levels 3,5, and 15 are shown below. Zooming into the vertex of the L at level 15 shows the error indicator remains strongest there. Further refinements will focus on that area.
To perform Dörfler marking, refine the current mesh, and add the mesh to the :class:`.AdaptiveMeshHierarchy`, we use the ``amh.adapt(eta, theta)`` method. In this method the input is the recently computed error estimator ``eta`` and the Dörfler marking parameter ``theta``. The method always performs this on the current fine mesh in the hierarchy. There is another method for adding a mesh to the hierarchy: ``amh.add_mesh(mesh)``. In this method, refinement on the mesh is performed externally by some custom procedure and the resulting mesh directly gets added to the hierarchy.
181
+
The meshes now refine according to the error estimator. The error estimators at levels 3,5, and 15 are shown below. Zooming into the vertex of the L-shape at level 15 shows the error indicator remains strongest there. Further refinements will focus on that area.
0 commit comments