diff --git a/chapter4/newton-solver.ipynb b/chapter4/newton-solver.ipynb index f979e173..750832f7 100644 --- a/chapter4/newton-solver.ipynb +++ b/chapter4/newton-solver.ipynb @@ -184,11 +184,24 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "0d5ce361", "metadata": {}, "outputs": [], "source": [ "solver = PETSc.KSP().create(mesh.comm)\n", + "solver.setType(\"preonly\")\n", + "solver.getPC().setType(\"lu\")\n", + "solver.getPC().setFactorSolverType(\"mumps\")\n", + "solver.setErrorIfNotConverged(True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ "solver.setOperators(A)\n", "du = dolfinx.fem.Function(V)" ] @@ -257,7 +270,7 @@ "\n", " # Compute norm of update\n", " correction_norm = du.x.petsc_vec.norm(0)\n", - " print(f\"Iteration {i}: Correction norm {correction_norm}\")\n", + " PETSc.Sys.Print(f\"Iteration {i}: Correction norm {correction_norm}\")\n", " if correction_norm < 1e-10:\n", " break\n", " solutions[i, :] = uh.x.array[sort_order]" @@ -279,7 +292,7 @@ "outputs": [], "source": [ "dolfinx.fem.petsc.assemble_vector(L, residual)\n", - "print(f\"Final residual {L.norm(0)}\")\n", + "PETSc.Sys.Print(f\"Final residual {L.norm(0)}\")\n", "A.destroy()\n", "L.destroy()\n", "solver.destroy()" @@ -314,7 +327,7 @@ " u_ex = root(x)\n", " L2_error = dolfinx.fem.form(ufl.inner(uh - u_ex, uh - u_ex) * ufl.dx)\n", " global_L2 = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_error), op=MPI.SUM)\n", - " print(f\"L2-error (root {j}) {np.sqrt(global_L2)}\")\n", + " PETSc.Sys.Print(f\"L2-error (root {j}) {np.sqrt(global_L2)}\")\n", "\n", " kwargs = {} if j == 0 else {\"label\": \"u_exact\"}\n", " plt.plot(x_spacing, root(x_spacing.reshape(1, -1)), *args, **kwargs)\n", @@ -410,7 +423,11 @@ "A = dolfinx.fem.petsc.create_matrix(jacobian)\n", "L = dolfinx.fem.petsc.create_vector(dolfinx.fem.extract_function_spaces(residual))\n", "solver = PETSc.KSP().create(mesh.comm)\n", - "solver.setOperators(A)" + "solver.setOperators(A)\n", + "solver.setType(\"preonly\")\n", + "solver.getPC().setType(\"lu\")\n", + "solver.getPC().setFactorSolverType(\"mumps\")\n", + "solver.setErrorIfNotConverged(True)" ] }, { @@ -461,15 +478,20 @@ " dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=[bc])\n", " A.assemble()\n", " dolfinx.fem.petsc.assemble_vector(L, residual)\n", + "\n", + " # Compute b - alpha * J(u_D-u_(i-1))\n", + " dolfinx.fem.petsc.apply_lifting(\n", + " L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=-1.0\n", + " )\n", " L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)\n", - " L.scale(-1)\n", "\n", - " # Compute b - J(u_D-u_(i-1))\n", - " dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=1)\n", - " # Set du|_bc = u_{i-1}-u_D\n", - " dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, 1.0)\n", + " # Set du|_bc = - (u_{i-1}-u_D)\n", + " dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, -1.0)\n", " L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)\n", "\n", + " # Compute negative residual\n", + " L.scale(-1)\n", + "\n", " # Solve linear problem\n", " solver.solve(L, du.x.petsc_vec)\n", " du.x.scatter_forward()\n", @@ -486,8 +508,10 @@ " np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM))\n", " )\n", " du_norm.append(correction_norm)\n", - "\n", - " print(f\"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}\")\n", + " PETSc.Sys.Print(\n", + " f\"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}\",\n", + " flush=True,\n", + " )\n", " if correction_norm < 1e-10:\n", " break" ] @@ -542,8 +566,7 @@ "outputs": [], "source": [ "error_max = domain.comm.allreduce(np.max(np.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}\")" + "PETSc.Sys.Print(f\"Error_max: {error_max:.2e}\")" ] }, { diff --git a/chapter4/newton-solver.py b/chapter4/newton-solver.py index 882ddae9..80b460af 100644 --- a/chapter4/newton-solver.py +++ b/chapter4/newton-solver.py @@ -6,7 +6,7 @@ # extension: .py # format_name: light # format_version: '1.5' -# jupytext_version: 1.18.1 +# jupytext_version: 1.19.1 # kernelspec: # display_name: Python 3 (ipykernel) # language: python @@ -99,6 +99,11 @@ def root_1(x): # Next, we create the linear solver and the vector to hold `du`. solver = PETSc.KSP().create(mesh.comm) +solver.setType("preonly") +solver.getPC().setType("lu") +solver.getPC().setFactorSolverType("mumps") +solver.setErrorIfNotConverged(True) + solver.setOperators(A) du = dolfinx.fem.Function(V) @@ -139,7 +144,7 @@ def root_1(x): # Compute norm of update correction_norm = du.x.petsc_vec.norm(0) - print(f"Iteration {i}: Correction norm {correction_norm}") + PETSc.Sys.Print(f"Iteration {i}: Correction norm {correction_norm}") if correction_norm < 1e-10: break solutions[i, :] = uh.x.array[sort_order] @@ -147,7 +152,7 @@ def root_1(x): # We now compute the magnitude of the residual. dolfinx.fem.petsc.assemble_vector(L, residual) -print(f"Final residual {L.norm(0)}") +PETSc.Sys.Print(f"Final residual {L.norm(0)}") A.destroy() L.destroy() solver.destroy() @@ -167,7 +172,7 @@ def root_1(x): u_ex = root(x) L2_error = dolfinx.fem.form(ufl.inner(uh - u_ex, uh - u_ex) * ufl.dx) global_L2 = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_error), op=MPI.SUM) - print(f"L2-error (root {j}) {np.sqrt(global_L2)}") + PETSc.Sys.Print(f"L2-error (root {j}) {np.sqrt(global_L2)}") kwargs = {} if j == 0 else {"label": "u_exact"} plt.plot(x_spacing, root(x_spacing.reshape(1, -1)), *args, **kwargs) @@ -227,6 +232,10 @@ def u_exact(x): L = dolfinx.fem.petsc.create_vector(dolfinx.fem.extract_function_spaces(residual)) solver = PETSc.KSP().create(mesh.comm) solver.setOperators(A) +solver.setType("preonly") +solver.getPC().setType("lu") +solver.getPC().setFactorSolverType("mumps") +solver.setErrorIfNotConverged(True) # Since this problem has strong Dirichlet conditions, we need to apply lifting to the right hand side of our Newton problem. # We previously had that we wanted to solve the system: @@ -262,15 +271,20 @@ def u_exact(x): dolfinx.fem.petsc.assemble_matrix(A, jacobian, bcs=[bc]) A.assemble() dolfinx.fem.petsc.assemble_vector(L, residual) + + # Compute b - alpha * J(u_D-u_(i-1)) + dolfinx.fem.petsc.apply_lifting( + L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=-1.0 + ) L.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) - L.scale(-1) - # Compute b - J(u_D-u_(i-1)) - dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], alpha=1) - # Set du|_bc = u_{i-1}-u_D - dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, 1.0) + # Set du|_bc = - (u_{i-1}-u_D) + dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, -1.0) L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD) + # Compute negative residual + L.scale(-1) + # Solve linear problem solver.solve(L, du.x.petsc_vec) du.x.scatter_forward() @@ -287,8 +301,10 @@ def u_exact(x): np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM)) ) du_norm.append(correction_norm) - - print(f"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}") + PETSc.Sys.Print( + f"Iteration {i}: Correction norm {correction_norm}, L2 error: {L2_error[-1]}", + flush=True, + ) if correction_norm < 1e-10: break @@ -315,8 +331,7 @@ def u_exact(x): # We compute the max error and plot the solution error_max = domain.comm.allreduce(np.max(np.abs(uh.x.array - u_D.x.array)), op=MPI.MAX) -if domain.comm.rank == 0: - print(f"Error_max: {error_max:.2e}") +PETSc.Sys.Print(f"Error_max: {error_max:.2e}") u_topology, u_cell_types, u_geometry = dolfinx.plot.vtk_mesh(V) u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)