changes to petsc vec access for dlx version 0.8.0

example/poisson_dirichlet_example.py

@@ -168,7 +168,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
     m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")
 
-    V_P1 = dlx.fem.FunctionSpace(msh, ("Lagrange", 1))
+    V_P1 = dlx.fem.functionspace(msh, ("Lagrange", 1))
 
     u_true_fun = dlx.fem.Function(V_P1, name="u_true")
     u_true_fun.interpolate(hpx.vector2Function(u_true, Vh[hpx.STATE]))
@@ -209,9 +209,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
             )
         )
 
-    temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
-    Omega = hpx.MultiVector(temp_para_vec, k + p)
-    temp_para_vec.destroy()
+    Omega = hpx.MultiVector(x[hpx.PARAMETER].petsc_vec, k + p)
 
     hpx.parRandom.normal(1.0, Omega)
 

example/poisson_dirichlet_example_reg.py

@@ -176,7 +176,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
     m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")
 
-    V_P1 = dlx.fem.FunctionSpace(msh, ("Lagrange", 1))
+    V_P1 = dlx.fem.functionspace(msh, ("Lagrange", 1))
 
     u_true_fun = dlx.fem.Function(V_P1, name="u_true")
     u_true_fun.interpolate(hpx.vector2Function(u_true, Vh[hpx.STATE]))
@@ -219,9 +219,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
             )
         )
 
-    temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
-    Omega = hpx.MultiVector(temp_para_vec, k + p)
-    temp_para_vec.destroy()
+    Omega = hpx.MultiVector(x[hpx.PARAMETER].petsc_vec, k + p)
 
     hpx.parRandom.normal(1.0, Omega)
 

example/poisson_example.py

@@ -148,7 +148,7 @@ def run_inversion(
     m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
     m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")
 
-    V_P1 = dlx.fem.FunctionSpace(msh, ("Lagrange", 1))
+    V_P1 = dlx.fem.functionspace(msh, ("Lagrange", 1))
 
     u_true_fun = dlx.fem.Function(V_P1, name="u_true")
     u_true_fun.interpolate(hpx.vector2Function(u_true, Vh[hpx.STATE]))
@@ -189,9 +189,7 @@ def run_inversion(
             )
         )
 
-    temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
-    Omega = hpx.MultiVector(temp_para_vec, k + p)
-    temp_para_vec.destroy()
+    Omega = hpx.MultiVector(x[hpx.PARAMETER].petsc_vec, k + p)
 
     hpx.parRandom.normal(1.0, Omega)
 

example/poisson_example_reg.py

@@ -151,7 +151,7 @@ def run_inversion(
     m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
     m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")
 
-    V_P1 = dlx.fem.FunctionSpace(msh, ("Lagrange", 1))
+    V_P1 = dlx.fem.functionspace(msh, ("Lagrange", 1))
 
     u_true_fun = dlx.fem.Function(V_P1, name="u_true")
     u_true_fun.interpolate(hpx.vector2Function(u_true, Vh[hpx.STATE]))
@@ -192,9 +192,7 @@ def run_inversion(
             )
         )
 
-    temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
-    Omega = hpx.MultiVector(temp_para_vec, k + p)
-    temp_para_vec.destroy()
+    Omega = hpx.MultiVector(x[hpx.PARAMETER].petsc_vec, k + p)
 
     hpx.parRandom.normal(1.0, Omega)
 

example/sfsi_toy_gaussian.py

@@ -176,7 +176,7 @@ def run_inversion(
     m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
     m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")
 
-    V_P1 = dlx.fem.FunctionSpace(msh, ("Lagrange", 1))
+    V_P1 = dlx.fem.functionspace(msh, ("Lagrange", 1))
 
     u_true_fun = dlx.fem.Function(V_P1, name="u_true")
     u_true_fun.interpolate(hpx.vector2Function(u_true, Vh[hpx.STATE]))
@@ -217,9 +217,7 @@ def run_inversion(
             )
         )
 
-    temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
-    Omega = hpx.MultiVector(temp_para_vec, k + p)
-    temp_para_vec.destroy()
+    Omega = hpx.MultiVector(x[hpx.PARAMETER].petsc_vec, k + p)
 
     hpx.parRandom.normal(1.0, Omega)
 

example/sfsi_toy_gaussian_reg.py

@@ -182,7 +182,7 @@ def run_inversion(
     m_fun = hpx.vector2Function(x[hpx.PARAMETER], Vh[hpx.PARAMETER], name="m_map")
     m_true_fun = hpx.vector2Function(m_true, Vh[hpx.PARAMETER], name="m_true")
 
-    V_P1 = dlx.fem.FunctionSpace(msh, ("Lagrange", 1))
+    V_P1 = dlx.fem.functionspace(msh, ("Lagrange", 1))
 
     u_true_fun = dlx.fem.Function(V_P1, name="u_true")
     u_true_fun.interpolate(hpx.vector2Function(u_true, Vh[hpx.STATE]))
@@ -223,9 +223,7 @@ def run_inversion(
             )
         )
 
-    temp_para_vec = dlx.la.create_petsc_vector_wrap(x[hpx.PARAMETER])
-    Omega = hpx.MultiVector(temp_para_vec, k + p)
-    temp_para_vec.destroy()
+    Omega = hpx.MultiVector(x[hpx.PARAMETER].petsc_vec, k + p)
 
     hpx.parRandom.normal(1.0, Omega)
 

hippylibX/algorithms/cgsolverSteihaug.py

@@ -161,11 +161,12 @@ def update_x_with_TR(self, x, alpha, d):
             return True
 
     def solve(self, b: dlx.la.Vector, x: dlx.la.Vector) -> None:
-        temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x)
-        temp_petsc_vec_b = dlx.la.create_petsc_vector_wrap(b)
-        self.solve_petsc(temp_petsc_vec_b, temp_petsc_vec_x)
-        temp_petsc_vec_x.destroy()
-        temp_petsc_vec_b.destroy()
+        # temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x)
+        # temp_petsc_vec_b = dlx.la.create_petsc_vector_wrap(b)
+        # self.solve_petsc(temp_petsc_vec_b, temp_petsc_vec_x)
+        # temp_petsc_vec_x.destroy()
+        # temp_petsc_vec_b.destroy()
+        self.solve_petsc(b.petsc_vec, x.petsc_vec)
 
     def solve_petsc(self, b: petsc4py.PETSc.Vec, x: petsc4py.PETSc.Vec) -> None:
         """

hippylibX/modeling/PDEProblem.py

@@ -130,11 +130,7 @@ def solveFwd(self, state: dlx.la.Vector, x: list) -> None:
                 petsc4py.PETSc.InsertMode.ADD_VALUES, petsc4py.PETSc.ScatterMode.REVERSE
             )
             dlx.fem.petsc.set_bc(b, self.bc)
-            state_vec = dlx.la.create_petsc_vector_wrap(state)
-
-            self.solver.solve(b, state_vec)
-
-            state_vec.destroy()
+            self.solver.solve(b, state.petsc_vec)
 
             A.destroy()
             b.destroy()
@@ -169,12 +165,8 @@ def solveAdj(
 
         self.solver.setOperators(Aadj)
 
-        temp_petsc_vec_adj = dlx.la.create_petsc_vector_wrap(adj)
-        temp_petsc_vec_adj_rhs = dlx.la.create_petsc_vector_wrap(adj_rhs)
+        self.solver.solve(adj_rhs.petsc_vec, adj.petsc_vec)
 
-        self.solver.solve(temp_petsc_vec_adj_rhs, temp_petsc_vec_adj)
-        temp_petsc_vec_adj.destroy()
-        temp_petsc_vec_adj_rhs.destroy()
         Aadj.destroy()
 
     def evalGradientParameter(self, x: list, out: dlx.la.Vector) -> None:
@@ -193,14 +185,13 @@ def evalGradientParameter(self, x: list, out: dlx.la.Vector) -> None:
         res_form = self.varf_handler(u, m, p)
 
         out.array[:] = 0.0
-        tmp_out = dlx.la.create_petsc_vector_wrap(out)
+
         dlx.fem.petsc.assemble_vector(
-            tmp_out, dlx.fem.form(ufl.derivative(res_form, m, dm))
+            out.petsc_vec, dlx.fem.form(ufl.derivative(res_form, m, dm))
         )
-        tmp_out.ghostUpdate(
+        out.petsc_vec.ghostUpdate(
             petsc4py.PETSc.InsertMode.ADD_VALUES, petsc4py.PETSc.ScatterMode.REVERSE
         )
-        tmp_out.destroy()
 
     def _createLUSolver(self) -> petsc4py.PETSc.KSP:
         ksp = petsc4py.PETSc.KSP().create(self.Vh[0].mesh.comm)
@@ -399,23 +390,18 @@ def apply_ij(self, i: int, j: int, dir: dlx.la.Vector, out: dlx.la.Vector) -> No
         KKT[ADJOINT, STATE] = self.A
         KKT[ADJOINT, PARAMETER] = self.C
 
-        temp_petsc_vec_out = dlx.la.create_petsc_vector_wrap(out)
-        temp_petsc_vec_dir = dlx.la.create_petsc_vector_wrap(dir)
-
         if i >= j:
             if KKT[i, j] is None:
-                temp_petsc_vec_out.scale(0.0)
+                out.petsc_vec.scale(0.0)
             else:
-                KKT[i, j].mult(temp_petsc_vec_dir, temp_petsc_vec_out)
+                KKT[i, j].mult(dir.petsc_vec, out.petsc_vec)
 
         else:
             if KKT[j, i] is None:
-                temp_petsc_vec_out.scale(0.0)
-            else:
-                KKT[j, i].multTranspose(temp_petsc_vec_dir, temp_petsc_vec_out)
+                out.petsc_vec.scale(0.0)
 
-        temp_petsc_vec_out.destroy()
-        temp_petsc_vec_dir.destroy()
+            else:
+                KKT[j, i].multTranspose(dir.petsc_vec, out.petsc_vec)
 
     def solveIncremental(
         self, out: dlx.la.Vector, rhs: dlx.la.Vector, is_adj: bool
@@ -434,16 +420,10 @@ def solveIncremental(
 
             .. math:: \\delta_{up} F(u, m, p; \\hat{u}, \\tilde{p}) = \\mbox{rhs},\\quad \\forall \\hat{u}.
         """
-        temp_petsc_vec_rhs = dlx.la.create_petsc_vector_wrap(rhs)
-        temp_petsc_vec_out = dlx.la.create_petsc_vector_wrap(out)
-
         if is_adj:
             self.n_calls["incremental_adjoint"] += 1
-            self.solver_adj_inc.solve(temp_petsc_vec_rhs, temp_petsc_vec_out)
+            self.solver_adj_inc.solve(rhs.petsc_vec, out.petsc_vec)
 
         else:
             self.n_calls["incremental_forward"] += 1
-            self.solver_fwd_inc.solve(temp_petsc_vec_rhs, temp_petsc_vec_out)
-
-        temp_petsc_vec_out.destroy()
-        temp_petsc_vec_rhs.destroy()
+            self.solver_fwd_inc.solve(rhs.petsc_vec, out.petsc_vec)

hippylibX/modeling/Regularization.py

@@ -106,19 +106,19 @@ def cost(self, m: dlx.la.Vector) -> float:
     def grad(self, m: dlx.la.Vector, out: dlx.la.Vector) -> dlx.la.Vector:
         hpx.updateFromVector(self.mfun, m)
         out.array[:] = 0.0
-        tmp_out = dlx.la.create_petsc_vector_wrap(out)
+
         dlx.fem.petsc.assemble_vector(
-            tmp_out,
+            out.petsc_vec,
             dlx.fem.form(
                 ufl.derivative(
                     self.functional_handler(self.mfun), self.mfun, self.mtest
                 )
             ),
         )
-        tmp_out.ghostUpdate(
+
+        out.petsc_vec.ghostUpdate(
             petsc4py.PETSc.InsertMode.ADD_VALUES, petsc4py.PETSc.ScatterMode.REVERSE
         )
-        tmp_out.destroy()
 
     def setLinearizationPoint(
         self, m: dlx.la.Vector, gauss_newton_approx=False

hippylibX/modeling/laplaceApproximation.py

@@ -102,14 +102,9 @@ def createVecLeft(self) -> petsc4py.PETSc.Vec:
         return self.prior.R.createVecLeft()
 
     def mult(self, noise: dlx.la.Vector, s: dlx.la.Vector):
-        temp_petsc_vec_noise = dlx.la.create_petsc_vector_wrap(noise)
-        temp_petsc_vec_s = dlx.la.create_petsc_vector_wrap(s)
-
-        self.prior.R.mult(temp_petsc_vec_noise, self.help)
-        self.lrsqrt.mult(self.help, temp_petsc_vec_s)
-        temp_petsc_vec_s.axpby(-1.0, 1.0, temp_petsc_vec_noise)
-        temp_petsc_vec_noise.destroy()
-        temp_petsc_vec_s.destroy()
+        self.prior.R.mult(noise.petsc_vec, self.help)
+        self.lrsqrt.mult(self.help, s.petsc_vec)
+        s.petsc_vec.axpby(-1.0, 1.0, noise.petsc_vec)
 
 
 class LaplaceApproximator:
@@ -152,11 +147,9 @@ def cost(self, m: dlx.la.Vector) -> float:
             dm = m
         else:
             dm = m - self.mean
-        temp_petsc_vec_dm = dlx.la.create_petsc_vector_wrap(dm)
-        self.Hlr.mult(temp_petsc_vec_dm, self.Hlr.help1)
-        return_value = 0.5 * self.Hlr.help1.dot(temp_petsc_vec_dm)
-        temp_petsc_vec_dm.destroy()
-        return return_value
+
+        self.Hlr.mult(dm.petsc_vec, self.Hlr.help1)
+        return 0.5 * self.Hlr.help1.dot(dm.petsc_vec)
 
     def sample(self, *args, **kwargs):
         """

hippylibX/modeling/misfit.py

@@ -41,16 +41,15 @@ def grad(self, i: int, x: list, out: dlx.la.Vector) -> None:
         x_fun = [u_fun, m_fun]
 
         out.array[:] = 0.0
-        tmp_out = dlx.la.create_petsc_vector_wrap(out)
+
         dlx.fem.petsc.assemble_vector(
-            tmp_out,
+            out.petsc_vec,
             dlx.fem.form(ufl.derivative(self.form(*x_fun), x_fun[i], self.x_test[i])),
         )
-        tmp_out.ghostUpdate(
+        out.petsc_vec.ghostUpdate(
             petsc4py.PETSc.InsertMode.ADD_VALUES, petsc4py.PETSc.ScatterMode.REVERSE
         )
-        dlx.fem.petsc.set_bc(tmp_out, self.bc0)
-        tmp_out.destroy()
+        dlx.fem.petsc.set_bc(out.petsc_vec, self.bc0)
 
     def setLinearizationPoint(self, x: list, gauss_newton_approx=False) -> None:
         hpx.updateFromVector(self.xfun[hpx.STATE], x[hpx.STATE])
@@ -74,11 +73,9 @@ def apply_ij(self, i: int, j: int, dir: dlx.la.Vector, out: dlx.la.Vector) -> No
             )
         )
         out.array[:] = 0.0
-        tmp_out = dlx.la.create_petsc_vector_wrap(out)
-        dlx.fem.petsc.assemble_vector(tmp_out, action)
-        tmp_out.ghostUpdate(
+        dlx.fem.petsc.assemble_vector(out.petsc_vec, action)
+        out.petsc_vec.ghostUpdate(
             petsc4py.PETSc.InsertMode.ADD_VALUES, petsc4py.PETSc.ScatterMode.REVERSE
         )
         if i == hpx.STATE:
-            dlx.fem.petsc.set_bc(tmp_out, self.bc0)
-        tmp_out.destroy()
+            dlx.fem.petsc.set_bc(out.petsc_vec, self.bc0)

hippylibX/modeling/model.py

@@ -160,24 +160,16 @@ def evalGradientParameter(
 
         self.misfit.grad(PARAMETER, x, tmp)
 
-        mg_petsc = dlx.la.create_petsc_vector_wrap(mg)
-        tmp_petsc = dlx.la.create_petsc_vector_wrap(tmp)
-
-        mg_petsc.axpy(1.0, tmp_petsc)
+        mg.petsc_vec.axpy(1.0, tmp.petsc_vec)
 
         if not misfit_only:
             self.prior.grad(x[PARAMETER], tmp)
 
-            mg_petsc.axpy(1.0, tmp_petsc)
-
-        self.prior.Msolver.solve(mg_petsc, tmp_petsc)
-
-        return_value = math.sqrt(mg_petsc.dot(tmp_petsc))
+            mg.petsc_vec.axpy(1.0, tmp.petsc_vec)
 
-        mg_petsc.destroy()
-        tmp_petsc.destroy()
+        self.prior.Msolver.solve(mg.petsc_vec, tmp.petsc_vec)
 
-        return return_value
+        return math.sqrt(mg.petsc_vec.dot(tmp.petsc_vec))
 
     def setPointForHessianEvaluations(self, x: list, gauss_newton_approx=False) -> None:
         """
@@ -320,11 +312,7 @@ def applyR(self, dm: dlx.la.Vector, out: dlx.la.Vector) -> None:
 
         .. note:: This routine assumes that :code:`out` has the correct shape.
         """
-        temp_petsc_vec_dm = dlx.la.create_petsc_vector_wrap(dm)
-        temp_petsc_vec_out = dlx.la.create_petsc_vector_wrap(out)
-        self.prior.R.mult(temp_petsc_vec_dm, temp_petsc_vec_out)
-        temp_petsc_vec_dm.destroy()
-        temp_petsc_vec_out.destroy()
+        self.prior.R.mult(dm.petsc_vec, out.petsc_vec)
 
     def Rsolver(self) -> Any:
         """