Change from vector to petsc_vec (#206)

chapter2/diffusion_code.ipynb

@@ -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",
@@ -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",

chapter2/diffusion_code.py

@@ -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.
@@ -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)

chapter2/heat_code.ipynb

@@ -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",

chapter2/heat_code.py

@@ -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)

chapter2/linearelasticity_code.ipynb

@@ -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",

chapter2/linearelasticity_code.py

@@ -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)

chapter2/ns_code1.ipynb

@@ -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",
@@ -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",
@@ -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",

chapter2/ns_code1.py

@@ -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
@@ -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[:]
@@ -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}")

chapter2/ns_code2.ipynb

@@ -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",
@@ -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",

chapter2/ns_code2.py

@@ -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
@@ -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)
 

chapter4/newton-solver.ipynb

@@ -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",
@@ -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",

chapter4/newton-solver.py

@@ -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
@@ -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)))