|
84 | 84 | "cell_type": "markdown", |
85 | 85 | "metadata": {}, |
86 | 86 | "source": [ |
87 | | - "The next step is to solve the non-linear problem. To do so, we use [Newtons method](https://en.wikipedia.org/wiki/Newton%27s_method). \n", |
88 | | - "We start by creating a class containing the core functions that we require to solve the problem." |
| 87 | + "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", |
| 88 | + "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)." |
89 | 89 | ] |
90 | 90 | }, |
91 | 91 | { |
|
94 | 94 | "metadata": {}, |
95 | 95 | "outputs": [], |
96 | 96 | "source": [ |
97 | | - "from typing import List\n", |
98 | | - "class NonlinearPDEProblem:\n", |
99 | | - " \"\"\"Nonlinear problem class for solving the non-linear problem\n", |
100 | | - " F(u, v) = 0 for all v in V\n", |
101 | | - " \"\"\"\n", |
102 | | - " def __init__(self, F: ufl.form.Form, u: dolfinx.Function, bcs: List[dolfinx.DirichletBC]):\n", |
103 | | - " \"\"\"\n", |
104 | | - " Input:\n", |
105 | | - " - F: The PDE residual F(u, v)\n", |
106 | | - " - u: The unknown\n", |
107 | | - " - bcs: List of Dirichlet boundary conditions\n", |
108 | | - " This class set up structures for solving the non-linear problem using Newton's method, \n", |
109 | | - " dF/du(u) du = -F(u)\n", |
110 | | - " \"\"\"\n", |
111 | | - " V = u.function_space\n", |
112 | | - " du = ufl.TrialFunction(V)\n", |
113 | | - " self.L = F\n", |
114 | | - " # Create the Jacobian matrix, dF/du\n", |
115 | | - " self.a = ufl.derivative(F, u, du)\n", |
116 | | - " self.bcs = bcs\n", |
117 | | - "\n", |
118 | | - " # Create matrix and vector to be used for assembly\n", |
119 | | - " # of the non-linear problem\n", |
120 | | - " self.matrix = dolfinx.fem.create_matrix(self.a)\n", |
121 | | - " self.vector = dolfinx.fem.create_vector(self.L)\n", |
122 | | - "\n", |
123 | | - " def form(self, x: PETSc.Vec):\n", |
124 | | - " \"\"\"\n", |
125 | | - " This function is called before the residual or Jacobian is computed. This is usually used to update ghost values.\n", |
126 | | - " Input: \n", |
127 | | - " x: The vector containing the latest solution\n", |
128 | | - " \"\"\"\n", |
129 | | - " x.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)\n", |
130 | | - "\n", |
131 | | - " def F(self, x: PETSc.Vec, b: PETSc.Vec):\n", |
132 | | - " \"\"\"Assemble the residual F into the vector b. \n", |
133 | | - " Input:\n", |
134 | | - " x: The vector containing the latest solution\n", |
135 | | - " b: Vector to assemble the residual into\n", |
136 | | - " \"\"\"\n", |
137 | | - " # Reset the residual vector\n", |
138 | | - " with b.localForm() as b_local:\n", |
139 | | - " b_local.set(0.0)\n", |
140 | | - " dolfinx.fem.assemble_vector(b, self.L)\n", |
141 | | - " # Apply boundary condition\n", |
142 | | - " dolfinx.fem.apply_lifting(b, [self.a], [self.bcs], [x], -1.0)\n", |
143 | | - " b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)\n", |
144 | | - " dolfinx.fem.set_bc(b, self.bcs, x, -1.0)\n", |
145 | | - "\n", |
146 | | - " def J(self, x: PETSc.Vec, A: PETSc.Mat):\n", |
147 | | - " \"\"\"Assemble the Jacobian matrix.\n", |
148 | | - " Input:\n", |
149 | | - " - x: The vector containing the latest solution\n", |
150 | | - " - A: The matrix to assemble the Jacobian into\n", |
151 | | - " \"\"\"\n", |
152 | | - " A.zeroEntries()\n", |
153 | | - " dolfinx.fem.assemble_matrix(A, self.a, self.bcs)\n", |
154 | | - " A.assemble()" |
| 97 | + "problem = dolfinx.fem.NonlinearProblem(F, uh, bcs=[bc])" |
155 | 98 | ] |
156 | 99 | }, |
157 | 100 | { |
158 | 101 | "cell_type": "markdown", |
159 | 102 | "metadata": {}, |
160 | 103 | "source": [ |
161 | | - "As we have now defined the the solution strategy, we can use the built-in Newton-solver in dolfin-X to solve the variational problem. We can set the convergence criterions for the solver by changing the absolute tolerance (`atol`), relative tolerance (`rtol`) or the convergence criterion (`residual` or `incremental`)." |
| 104 | + "Next, we use the dolfinx Newton solver. We can set the convergence criterions for the solver by changing the absolute tolerance (`atol`), relative tolerance (`rtol`) or the convergence criterion (`residual` or `incremental`)." |
162 | 105 | ] |
163 | 106 | }, |
164 | 107 | { |
165 | 108 | "cell_type": "code", |
166 | 109 | "execution_count": 5, |
167 | 110 | "metadata": {}, |
| 111 | + "outputs": [], |
| 112 | + "source": [ |
| 113 | + "solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)\n", |
| 114 | + "solver.convergence_criterion = \"incremental\"\n", |
| 115 | + "solver.rtol = 1e-6\n", |
| 116 | + "solver.report = True" |
| 117 | + ] |
| 118 | + }, |
| 119 | + { |
| 120 | + "cell_type": "markdown", |
| 121 | + "metadata": {}, |
| 122 | + "source": [ |
| 123 | + "We can modify the linear solver in each Newton iteration by accessing the underlying `PETSc` object." |
| 124 | + ] |
| 125 | + }, |
| 126 | + { |
| 127 | + "cell_type": "code", |
| 128 | + "execution_count": 6, |
| 129 | + "metadata": {}, |
| 130 | + "outputs": [], |
| 131 | + "source": [ |
| 132 | + "ksp = solver.krylov_solver\n", |
| 133 | + "opts = PETSc.Options()\n", |
| 134 | + "option_prefix = ksp.getOptionsPrefix()\n", |
| 135 | + "opts[f\"{option_prefix}ksp_type\"] = \"cg\"\n", |
| 136 | + "opts[f\"{option_prefix}pc_type\"] = \"gamg\"\n", |
| 137 | + "opts[f\"{option_prefix}pc_factor_mat_solver_type\"] = \"mumps\"\n", |
| 138 | + "ksp.setFromOptions()" |
| 139 | + ] |
| 140 | + }, |
| 141 | + { |
| 142 | + "cell_type": "markdown", |
| 143 | + "metadata": {}, |
| 144 | + "source": [ |
| 145 | + "We are now ready to solve the non-linear problem. We assert that the solver has converged and print the number of iterations." |
| 146 | + ] |
| 147 | + }, |
| 148 | + { |
| 149 | + "cell_type": "code", |
| 150 | + "execution_count": 7, |
| 151 | + "metadata": {}, |
168 | 152 | "outputs": [ |
169 | 153 | { |
170 | 154 | "name": "stdout", |
|
175 | 159 | } |
176 | 160 | ], |
177 | 161 | "source": [ |
178 | | - "# Create nonlinear problem\n", |
179 | | - "problem = NonlinearPDEProblem(F, uh, [bc])\n", |
180 | | - "\n", |
181 | | - "# Create Newton solver\n", |
182 | | - "solver = dolfinx.cpp.nls.NewtonSolver(MPI.COMM_WORLD)\n", |
183 | | - "\n", |
184 | | - "# Set Newton solver options\n", |
185 | | - "solver.atol = 1e-8\n", |
186 | | - "solver.rtol = 1e-8\n", |
187 | | - "solver.convergence_criterion = \"residual\"\n", |
188 | | - "\n", |
189 | | - "# Set non-linear problem for Newton solver\n", |
190 | | - "solver.setF(problem.F, problem.vector)\n", |
191 | | - "solver.setJ(problem.J, problem.matrix)\n", |
192 | | - "solver.set_form(problem.form)\n", |
193 | | - "\n", |
194 | | - "# Solve non-linear problem\n", |
195 | 162 | "dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)\n", |
196 | | - "n, converged = solver.solve(uh.vector)\n", |
| 163 | + "n, converged = solver.solve(uh)\n", |
197 | 164 | "dolfinx.cpp.la.scatter_forward(uh.x)\n", |
198 | 165 | "assert(converged)\n", |
199 | 166 | "print(f\"Number of interations: {n:d}\")" |
|
209 | 176 | }, |
210 | 177 | { |
211 | 178 | "cell_type": "code", |
212 | | - "execution_count": 1, |
| 179 | + "execution_count": 8, |
213 | 180 | "metadata": {}, |
214 | 181 | "outputs": [ |
215 | 182 | { |
216 | | - "ename": "NameError", |
217 | | - "evalue": "name 'dolfinx' is not defined", |
218 | | - "output_type": "error", |
219 | | - "traceback": [ |
220 | | - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", |
221 | | - "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", |
222 | | - "\u001b[0;32m<ipython-input-1-1940f870bef7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m# Compute L2 error and error at nodes\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mV_ex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdolfinx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFunctionSpace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmesh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m\"CG\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0mu_ex\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdolfinx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mFunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mV_ex\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mu_ex\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu_exact\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mdolfinx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcpp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mla\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter_forward\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu_ex\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", |
223 | | - "\u001b[0;31mNameError\u001b[0m: name 'dolfinx' is not defined" |
| 183 | + "name": "stdout", |
| 184 | + "output_type": "stream", |
| 185 | + "text": [ |
| 186 | + "L2-error: 3.65e-15\n", |
| 187 | + "Error_max: 1.42e-14\n" |
224 | 188 | ] |
225 | 189 | } |
226 | 190 | ], |
|
0 commit comments