Skip to content

Update newton solver to be safe in parallel #305

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jorgensd merged 1 commit into from
Feb 26, 2026
Merged

Update newton solver to be safe in parallel #305

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 37 additions & 14 deletions chapter4/newton-solver.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)"
]
Expand Down Expand Up @@ -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]"
Expand All @@ -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()"
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -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"
]
Expand Down Expand Up @@ -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}\")"
]
},
{
Expand Down
41 changes: 28 additions & 13 deletions chapter4/newton-solver.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -139,15 +144,15 @@ 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]

# 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()
Expand All @@ -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)
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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()
Expand All @@ -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

Expand All @@ -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)
Expand Down
Loading