Skip to content

Change from vector to petsc_vec #206

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.

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
merged 1 commit into from
Sep 26, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
4 changes: 2 additions & 2 deletions chapter2/diffusion_code.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@
"source": [
"## Updating the solution and right hand side per time step\n",
"To be able to solve the variation problem at each time step, we have to assemble the right hand side and apply the boundary condition before calling\n",
"`solver.solve(b, uh.vector)`. We start by resetting the values in `b` as we are reusing the vector at every time step.\n",
"`solver.solve(b, uh.x.petsc_vec)`. We start by resetting the values in `b` as we are reusing the vector at every time step.\n",
"The next step is to assemble the vector calling `dolfinx.fem.petsc.assemble_vector(b, L)`, which means that we are assembling the linear form `L(v)` into the vector `b`. Note that we do not supply the boundary conditions for assembly, as opposed to the left hand side.\n",
"This is because we want to use lifting to apply the boundary condition, which preserves symmetry of the matrix $A$ in the bilinear form $a(u,v)=a(v,u)$ without Dirichlet boundary conditions.\n",
"When we have applied the boundary condition, we can solve the linear system and update values that are potentially shared between processors.\n",
Expand All @@ -249,7 +249,7 @@
" set_bc(b, [bc])\n",
"\n",
" # Solve linear problem\n",
" solver.solve(b, uh.vector)\n",
" solver.solve(b, uh.x.petsc_vec)\n",
" uh.x.scatter_forward()\n",
"\n",
" # Update solution at previous time step (u_n)\n",
Expand Down
4 changes: 2 additions & 2 deletions chapter2/diffusion_code.py
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ def initial_condition(x, a=5):

# ## Updating the solution and right hand side per time step
# To be able to solve the variation problem at each time step, we have to assemble the right hand side and apply the boundary condition before calling
# `solver.solve(b, uh.vector)`. We start by resetting the values in `b` as we are reusing the vector at every time step.
# `solver.solve(b, uh.x.petsc_vec)`. We start by resetting the values in `b` as we are reusing the vector at every time step.
# The next step is to assemble the vector calling `dolfinx.fem.petsc.assemble_vector(b, L)`, which means that we are assembling the linear form `L(v)` into the vector `b`. Note that we do not supply the boundary conditions for assembly, as opposed to the left hand side.
# This is because we want to use lifting to apply the boundary condition, which preserves symmetry of the matrix $A$ in the bilinear form $a(u,v)=a(v,u)$ without Dirichlet boundary conditions.
# When we have applied the boundary condition, we can solve the linear system and update values that are potentially shared between processors.
Expand All @@ -162,7 +162,7 @@ def initial_condition(x, a=5):
set_bc(b, [bc])

# Solve linear problem
solver.solve(b, uh.vector)
solver.solve(b, uh.x.petsc_vec)
uh.x.scatter_forward()

# Update solution at previous time step (u_n)
Expand Down
2 changes: 1 addition & 1 deletion chapter2/heat_code.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -242,7 +242,7 @@
" set_bc(b, [bc])\n",
"\n",
" # Solve linear problem\n",
" solver.solve(b, uh.vector)\n",
" solver.solve(b, uh.x.petsc_vec)\n",
" uh.x.scatter_forward()\n",
"\n",
" # Update solution at previous time step (u_n)\n",
Expand Down
2 changes: 1 addition & 1 deletion chapter2/heat_code.py
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@ def __call__(self, x):
set_bc(b, [bc])

# Solve linear problem
solver.solve(b, uh.vector)
solver.solve(b, uh.x.petsc_vec)
uh.x.scatter_forward()

# Update solution at previous time step (u_n)
Expand Down
2 changes: 1 addition & 1 deletion chapter2/linearelasticity_code.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -314,7 +314,7 @@
}
],
"source": [
"warped.cell_data[\"VonMises\"] = stresses.vector.array\n",
"warped.cell_data[\"VonMises\"] = stresses.x.petsc_vec.array\n",
"warped.set_active_scalars(\"VonMises\")\n",
"p = pyvista.Plotter()\n",
"p.add_mesh(warped)\n",
Expand Down
2 changes: 1 addition & 1 deletion chapter2/linearelasticity_code.py
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ def sigma(u):

# In the previous sections, we have only visualized first order Lagrangian functions. However, the Von Mises stresses are piecewise constant on each cell. Therefore, we modify our plotting routine slightly. The first thing we notice is that we now set values for each cell, which has a one to one correspondence with the degrees of freedom in the function space.

warped.cell_data["VonMises"] = stresses.vector.array
warped.cell_data["VonMises"] = stresses.x.petsc_vec.array
warped.set_active_scalars("VonMises")
p = pyvista.Plotter()
p.add_mesh(warped)
Expand Down
8 changes: 4 additions & 4 deletions chapter2/ns_code1.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -494,7 +494,7 @@
" apply_lifting(b1, [a1], [bcu])\n",
" b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
" set_bc(b1, bcu)\n",
" solver1.solve(b1, u_.vector)\n",
" solver1.solve(b1, u_.x.petsc_vec)\n",
" u_.x.scatter_forward()\n",
"\n",
" # Step 2: Pressure corrrection step\n",
Expand All @@ -504,15 +504,15 @@
" apply_lifting(b2, [a2], [bcp])\n",
" b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
" set_bc(b2, bcp)\n",
" solver2.solve(b2, p_.vector)\n",
" solver2.solve(b2, p_.x.petsc_vec)\n",
" p_.x.scatter_forward()\n",
"\n",
" # Step 3: Velocity correction step\n",
" with b3.localForm() as loc_3:\n",
" loc_3.set(0)\n",
" assemble_vector(b3, L3)\n",
" b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
" solver3.solve(b3, u_.vector)\n",
" solver3.solve(b3, u_.x.petsc_vec)\n",
" u_.x.scatter_forward()\n",
" # Update variable with solution form this time step\n",
" u_n.x.array[:] = u_.x.array[:]\n",
Expand All @@ -524,7 +524,7 @@
"\n",
" # Compute error at current time-step\n",
" error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))\n",
" error_max = mesh.comm.allreduce(np.max(u_.vector.array - u_ex.vector.array), op=MPI.MAX)\n",
" error_max = mesh.comm.allreduce(np.max(u_.x.petsc_vec.array - u_ex.x.petsc_vec.array), op=MPI.MAX)\n",
" # Print error only every 20th step and at the last step\n",
" if (i % 20 == 0) or (i == num_steps - 1):\n",
" print(f\"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}\")\n",
Expand Down
8 changes: 4 additions & 4 deletions chapter2/ns_code1.py
Original file line number Diff line number Diff line change
Expand Up @@ -283,7 +283,7 @@ def u_exact(x):
apply_lifting(b1, [a1], [bcu])
b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b1, bcu)
solver1.solve(b1, u_.vector)
solver1.solve(b1, u_.x.petsc_vec)
u_.x.scatter_forward()

# Step 2: Pressure corrrection step
Expand All @@ -293,15 +293,15 @@ def u_exact(x):
apply_lifting(b2, [a2], [bcp])
b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b2, bcp)
solver2.solve(b2, p_.vector)
solver2.solve(b2, p_.x.petsc_vec)
p_.x.scatter_forward()

# Step 3: Velocity correction step
with b3.localForm() as loc_3:
loc_3.set(0)
assemble_vector(b3, L3)
b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver3.solve(b3, u_.vector)
solver3.solve(b3, u_.x.petsc_vec)
u_.x.scatter_forward()
# Update variable with solution form this time step
u_n.x.array[:] = u_.x.array[:]
Expand All @@ -313,7 +313,7 @@ def u_exact(x):

# Compute error at current time-step
error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))
error_max = mesh.comm.allreduce(np.max(u_.vector.array - u_ex.vector.array), op=MPI.MAX)
error_max = mesh.comm.allreduce(np.max(u_.x.petsc_vec.array - u_ex.x.petsc_vec.array), op=MPI.MAX)
# Print error only every 20th step and at the last step
if (i % 20 == 0) or (i == num_steps - 1):
print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}")
Expand Down
10 changes: 5 additions & 5 deletions chapter2/ns_code2.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -655,7 +655,7 @@
" apply_lifting(b1, [a1], [bcu])\n",
" b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
" set_bc(b1, bcu)\n",
" solver1.solve(b1, u_s.vector)\n",
" solver1.solve(b1, u_s.x.petsc_vec)\n",
" u_s.x.scatter_forward()\n",
"\n",
" # Step 2: Pressure corrrection step\n",
Expand All @@ -665,26 +665,26 @@
" apply_lifting(b2, [a2], [bcp])\n",
" b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
" set_bc(b2, bcp)\n",
" solver2.solve(b2, phi.vector)\n",
" solver2.solve(b2, phi.x.petsc_vec)\n",
" phi.x.scatter_forward()\n",
"\n",
" p_.vector.axpy(1, phi.vector)\n",
" p_.x.petsc_vec.axpy(1, phi.x.petsc_vec)\n",
" p_.x.scatter_forward()\n",
"\n",
" # Step 3: Velocity correction step\n",
" with b3.localForm() as loc:\n",
" loc.set(0)\n",
" assemble_vector(b3, L3)\n",
" b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)\n",
" solver3.solve(b3, u_.vector)\n",
" solver3.solve(b3, u_.x.petsc_vec)\n",
" u_.x.scatter_forward()\n",
"\n",
" # Write solutions to file\n",
" vtx_u.write(t)\n",
" vtx_p.write(t)\n",
"\n",
" # Update variable with solution form this time step\n",
" with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:\n",
" with u_.x.petsc_vec.localForm() as loc_, u_n.x.petsc_vec.localForm() as loc_n, u_n1.x.petsc_vec.localForm() as loc_n1:\n",
" loc_n.copy(loc_n1)\n",
" loc_.copy(loc_n)\n",
"\n",
Expand Down
10 changes: 5 additions & 5 deletions chapter2/ns_code2.py
Original file line number Diff line number Diff line change
Expand Up @@ -418,7 +418,7 @@ def __call__(self, x):
apply_lifting(b1, [a1], [bcu])
b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b1, bcu)
solver1.solve(b1, u_s.vector)
solver1.solve(b1, u_s.x.petsc_vec)
u_s.x.scatter_forward()

# Step 2: Pressure corrrection step
Expand All @@ -428,26 +428,26 @@ def __call__(self, x):
apply_lifting(b2, [a2], [bcp])
b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b2, bcp)
solver2.solve(b2, phi.vector)
solver2.solve(b2, phi.x.petsc_vec)
phi.x.scatter_forward()

p_.vector.axpy(1, phi.vector)
p_.x.petsc_vec.axpy(1, phi.x.petsc_vec)
p_.x.scatter_forward()

# Step 3: Velocity correction step
with b3.localForm() as loc:
loc.set(0)
assemble_vector(b3, L3)
b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
solver3.solve(b3, u_.vector)
solver3.solve(b3, u_.x.petsc_vec)
u_.x.scatter_forward()

# Write solutions to file
vtx_u.write(t)
vtx_p.write(t)

# Update variable with solution form this time step
with u_.vector.localForm() as loc_, u_n.vector.localForm() as loc_n, u_n1.vector.localForm() as loc_n1:
with u_.x.petsc_vec.localForm() as loc_, u_n.x.petsc_vec.localForm() as loc_n, u_n1.x.petsc_vec.localForm() as loc_n1:
loc_n.copy(loc_n1)
loc_.copy(loc_n)

Expand Down
12 changes: 6 additions & 6 deletions chapter4/newton-solver.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -258,14 +258,14 @@
" L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)\n",
"\n",
" # Solve linear problem\n",
" solver.solve(L, du.vector)\n",
" solver.solve(L, du.x.petsc_vec)\n",
" du.x.scatter_forward()\n",
" # Update u_{i+1} = u_i + delta u_i\n",
" uh.x.array[:] += du.x.array\n",
" i += 1\n",
"\n",
" # Compute norm of update\n",
" correction_norm = du.vector.norm(0)\n",
" correction_norm = du.x.petsc_vec.norm(0)\n",
" print(f\"Iteration {i}: Correction norm {correction_norm}\")\n",
" if correction_norm < 1e-10:\n",
" break\n",
Expand Down Expand Up @@ -516,21 +516,21 @@
" L.scale(-1)\n",
"\n",
" # Compute b - J(u_D-u_(i-1))\n",
" dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.vector], scale=1)\n",
" dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], scale=1)\n",
" # Set du|_bc = u_{i-1}-u_D\n",
" dolfinx.fem.petsc.set_bc(L, [bc], uh.vector, 1.0)\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",
" # Solve linear problem\n",
" solver.solve(L, du.vector)\n",
" solver.solve(L, du.x.petsc_vec)\n",
" du.x.scatter_forward()\n",
"\n",
" # Update u_{i+1} = u_i + delta u_i\n",
" uh.x.array[:] += du.x.array\n",
" i += 1\n",
"\n",
" # Compute norm of update\n",
" correction_norm = du.vector.norm(0)\n",
" correction_norm = du.x.petsc_vec.norm(0)\n",
"\n",
" # Compute L2 error comparing to the analytical solution\n",
" L2_error.append(np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM)))\n",
Expand Down
12 changes: 6 additions & 6 deletions chapter4/newton-solver.py
Original file line number Diff line number Diff line change
Expand Up @@ -129,14 +129,14 @@ def root_1(x):
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

# Solve linear problem
solver.solve(L, du.vector)
solver.solve(L, du.x.petsc_vec)
du.x.scatter_forward()
# Update u_{i+1} = u_i + delta u_i
uh.x.array[:] += du.x.array
i += 1

# Compute norm of update
correction_norm = du.vector.norm(0)
correction_norm = du.x.petsc_vec.norm(0)
print(f"Iteration {i}: Correction norm {correction_norm}")
if correction_norm < 1e-10:
break
Expand Down Expand Up @@ -257,21 +257,21 @@ def u_exact(x):
L.scale(-1)

# Compute b - J(u_D-u_(i-1))
dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.vector], scale=1)
dolfinx.fem.petsc.apply_lifting(L, [jacobian], [[bc]], x0=[uh.x.petsc_vec], scale=1)
# Set du|_bc = u_{i-1}-u_D
dolfinx.fem.petsc.set_bc(L, [bc], uh.vector, 1.0)
dolfinx.fem.petsc.set_bc(L, [bc], uh.x.petsc_vec, 1.0)
L.ghostUpdate(addv=PETSc.InsertMode.INSERT_VALUES, mode=PETSc.ScatterMode.FORWARD)

# Solve linear problem
solver.solve(L, du.vector)
solver.solve(L, du.x.petsc_vec)
du.x.scatter_forward()

# Update u_{i+1} = u_i + delta u_i
uh.x.array[:] += du.x.array
i += 1

# Compute norm of update
correction_norm = du.vector.norm(0)
correction_norm = du.x.petsc_vec.norm(0)

# Compute L2 error comparing to the analytical solution
L2_error.append(np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(error), op=MPI.SUM)))
Expand Down