Dualspace update

firedrake/adjoint_utils/blocks/dirichlet_bc.py

@@ -40,17 +40,15 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
         for adj_input in adj_inputs:
             if self.compat.isconstant(c):
                 adj_value = self.backend.Function(self.parent_space.dual())
-                adj_input.apply(adj_value.vector())
+                adj_input.apply(adj_value)
                 if self.function_space != self.parent_space:
                     vec = self.compat.extract_bc_subvector(
                         adj_value, self.collapsed_space, bc
                     )
-                    adj_value = self.compat.function_from_vector(
-                        self.collapsed_space, vec
-                    )
+                    adj_value = self.backend.Function(self.collapsed_space, vec.dat)
 
                 if adj_value.ufl_shape == () or adj_value.ufl_shape[0] <= 1:
-                    r = adj_value.vector().sum()
+                    r = adj_value.dat.data_ro.sum()
                 else:
                     output = []
                     subindices = _extract_subindices(self.function_space)
@@ -64,7 +62,7 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
                         output.append(
                             current_subfunc.sub(
                                 prev_idx, deepcopy=True
-                            ).vector().sum()
+                            ).dat.data_ro.sum()
                         )
 
                     r = self.backend.cpp.la.Vector(self.backend.MPI.comm_world,
@@ -78,13 +76,13 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
                 # For now we will just assume the FunctionSpace is the same for
                 # the BC and the Function.
                 adj_value = self.backend.Function(self.parent_space.dual())
-                adj_input.apply(adj_value.vector())
+                adj_input.apply(adj_value)
                 r = self.compat.extract_bc_subvector(
                     adj_value, c.function_space(), bc
                 )
             elif isinstance(c, self.compat.Expression):
                 adj_value = self.backend.Function(self.parent_space.dual())
-                adj_input.apply(adj_value.vector())
+                adj_input.apply(adj_value)
                 output = self.compat.extract_bc_subvector(
                     adj_value, self.collapsed_space, bc
                 )

firedrake/adjoint_utils/blocks/function.py

@@ -54,7 +54,7 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
             if isinstance(block_variable.output, AdjFloat):
                 try:
                     # Adjoint of a broadcast is just a sum
-                    return adj_inputs[0].vector().sum()
+                    return adj_inputs[0].dat.data_ro.sum()
                 except AttributeError:
                     # Catch the case where adj_inputs[0] is just a float
                     return adj_inputs[0]
@@ -67,8 +67,7 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
                 adj_output = self.backend.Function(
                     block_variable.output.function_space())
                 adj_output.assign(prepared)
-                adj_output = self.backend.Cofunction(adj_output.function_space().dual(),
-                                                     val=adj_output.vector())
+                adj_output = adj_output.riesz_representation(riesz_map="l2")
                 return adj_output
         else:
             # Linear combination
@@ -82,8 +81,7 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
                 )
                 # Firedrake does not support assignment of conjugate functions
                 adj_output.interpolate(ufl.conj(diff_expr))
-                adj_output = self.backend.Cofunction(adj_output.function_space().dual(),
-                                                     val=adj_output.vector())
+                adj_output = adj_output.riesz_representation(riesz_map="l2")
             else:
                 mesh = adj_output.function_space().mesh()
                 diff_expr = ufl.algorithms.expand_derivatives(
@@ -94,7 +92,7 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx,
                     )
                 )
                 adj_output.assign(diff_expr)
-                return adj_output.vector().inner(adj_input_func.vector())
+                return adj_output.dat.inner(adj_input_func.dat)
 
             if self.compat.isconstant(block_variable.output):
                 R = block_variable.output._ad_function_space(
@@ -109,14 +107,14 @@ def _adj_assign_constant(self, adj_output, constant_fs):
         shape = r.ufl_shape
         if shape == () or shape[0] == 1:
             # Scalar Constant
-            r.vector()[:] = adj_output.vector().sum()
+            r.dat.data[:] = adj_output.dat.data_ro.sum()
         else:
             # We assume the shape of the constant == shape of the output
             # function if not scalar. This assumption is due to FEniCS not
             # supporting products with non-scalar constants in assign.
             values = []
             for i in range(shape[0]):
-                values.append(adj_output.sub(i, deepcopy=True).vector().sum())
+                values.append(adj_output.sub(i, deepcopy=True).dat.data_ro.sum())
             r.assign(self.backend.Constant(values))
         return r
 
@@ -139,7 +137,7 @@ def evaluate_tlm_component(self, inputs, tlm_inputs, block_variable, idx,
                 dudmi.assign(ufl.algorithms.expand_derivatives(
                     ufl.derivative(expr, dep.saved_output,
                                    dep.tlm_value)))
-                dudm.vector().axpy(1.0, dudmi.vector())
+                dudm.dat += 1.0 * dudmi.dat
 
         return dudm
 

firedrake/adjoint_utils/blocks/solving.py

@@ -151,7 +151,7 @@ def prepare_evaluate_adj(self, inputs, adj_inputs, relevant_dependencies):
             )
         )
         dFdu_form = self.backend.adjoint(dFdu)
-        dJdu = dJdu.vector().copy()
+        dJdu = dJdu.copy()
 
         compute_bdy = self._should_compute_boundary_adjoint(
             relevant_dependencies
@@ -184,14 +184,14 @@ def _assemble_and_solve_adj_eq(self, dFdu_adj_form, dJdu, compute_bdy):
 
         adj_sol = self.compat.create_function(self.function_space)
         self.compat.linalg_solve(
-            dFdu, adj_sol.vector(), dJdu, *self.adj_args, **self.adj_kwargs
+            dFdu, adj_sol, dJdu, *self.adj_args, **self.adj_kwargs
         )
 
         adj_sol_bdy = None
         if compute_bdy:
-            adj_sol_bdy = self.compat.function_from_vector(
+            adj_sol_bdy = self.backend.Function(
                 self.function_space.dual(),
-                dJdu_copy - self.compat.assemble_adjoint_value(self.backend.action(dFdu_adj_form, adj_sol)).vector()
+                dJdu_copy.dat - self.compat.assemble_adjoint_value(self.backend.action(dFdu_adj_form, adj_sol)).dat
             )
 
         return adj_sol, adj_sol_bdy
@@ -359,8 +359,7 @@ def _assemble_and_solve_soa_eq(self, dFdu_form, adj_sol, hessian_input,
         b = self._assemble_soa_eq_rhs(dFdu_form, adj_sol, hessian_input,
                                       d2Fdu2)
         dFdu_form = self.backend.adjoint(dFdu_form)
-        adj_sol2, adj_sol2_bdy = self._assemble_and_solve_adj_eq(dFdu_form, b.vector(),
-                                                                 compute_bdy)
+        adj_sol2, adj_sol2_bdy = self._assemble_and_solve_adj_eq(dFdu_form, b, compute_bdy)
         if self.adj2_cb is not None:
             self.adj2_cb(adj_sol2)
         if self.adj2_bdy_cb is not None and compute_bdy:
@@ -523,7 +522,7 @@ def _assembled_solve(self, lhs, rhs, func, bcs, **kwargs):
         for bc in bcs:
             bc.apply(rhs_func)
         rhs.assign(rhs_func.riesz_representation(riesz_map="l2"))
-        self.backend.solve(lhs, func.vector(), rhs, **kwargs)
+        self.backend.solve(lhs, func, rhs, **kwargs)
         return func
 
     def recompute_component(self, inputs, block_variable, idx, prepared):
@@ -662,7 +661,7 @@ def prepare_evaluate_adj(self, inputs, adj_inputs, relevant_dependencies):
         F_form = self._create_F_form()
 
         dFdu_form = self.adj_F
-        dJdu = dJdu.vector().copy()
+        dJdu = dJdu.copy()
 
         # Replace the form coefficients with checkpointed values.
         replace_map = self._replace_map(dFdu_form)

firedrake/adjoint_utils/function.py

@@ -234,7 +234,7 @@ def _ad_convert_riesz(self, value, options=None):
             raise TypeError("Expected a Cofunction, Function or a float")
 
         if riesz_representation == "l2":
-            value = value.vector() if isinstance(value, (Cofunction, Function)) else value
+            value = value.dat if isinstance(value, (Cofunction, Function)) else value
             return Function(V, val=value)
 
         elif riesz_representation in ("L2", "H1"):
@@ -308,7 +308,7 @@ def _ad_dot(self, other, options=None):
         options = {} if options is None else options
         riesz_representation = options.get("riesz_representation", "l2")
         if riesz_representation == "l2":
-            return self.vector().inner(other.vector())
+            return self.dat.inner(other.dat)
         elif riesz_representation == "L2":
             return assemble(firedrake.inner(self, other)*firedrake.dx)
         elif riesz_representation == "H1":

firedrake/cofunction.py

@@ -47,7 +47,7 @@ def __init__(self, function_space, val=None, name=None, dtype=ScalarType):
             V = V.function_space()
             # Deep copy prevents modifications to Vector copies.
             # Also, this discard the value of `val` if it was specified (consistent with Function)
-            val = function_space.copy(deepcopy=True).vector()
+            val = function_space.copy(deepcopy=True).dat
         elif not isinstance(V, functionspaceimpl.FiredrakeDualSpace):
             raise NotImplementedError("Can't make a Cofunction defined on a "
                                       + str(type(function_space)))
@@ -79,12 +79,12 @@ def __del__(self):
         if hasattr(self, "_comm"):
             mpi.decref(self._comm)
 
-    def copy(self, deepcopy=False):
+    def copy(self, deepcopy=True):
         r"""Return a copy of this :class:`firedrake.function.CoordinatelessFunction`.
 
-        :kwarg deepcopy: If ``True``, the new
+        :kwarg deepcopy: If ``True``, the default, the new
             :class:`firedrake.function.CoordinatelessFunction` will allocate new space
-            and copy values.  If ``False``, the default, then the new
+            and copy values.  If ``False``, then the new
             :class:`firedrake.function.CoordinatelessFunction` will share the dof values.
         """
         if deepcopy:

firedrake/function.py

@@ -487,7 +487,7 @@ def riesz_representation(self, riesz_map='L2'):
 
         V = self.function_space()
         if riesz_map == "l2":
-            return Cofunction(V.dual(), val=self.vector())
+            return Cofunction(V.dual(), val=self.dat)
 
         elif riesz_map in ("L2", "H1"):
             a = self._define_riesz_map_form(riesz_map, V)

firedrake/matrix_free/operators.py

@@ -188,9 +188,11 @@ def _assemble_diagonal(self):
 
     def getDiagonal(self, mat, vec):
         self._assemble_diagonal()
+        diagonal_func = self._diagonal.riesz_representation(riesz_map="l2")
         for bc in self.bcs:
             # Operator is identity on boundary nodes
-            bc.set(self._diagonal, 1)
+            bc.set(diagonal_func, 1)
+        self._diagonal.assign(diagonal_func.riesz_representation(riesz_map="l2"))
         with self._diagonal.dat.vec_ro as v:
             v.copy(vec)
 

firedrake/ufl_expr.py

@@ -109,12 +109,17 @@ def function_space(self):
     def make_dat(self):
         return self.function_space().make_dat()
 
-    def _analyze_form_arguments(self):
+    def _analyze_form_arguments(self, outer_form=None):
         # Returns the argument found in the Coargument object
+        self._coefficients = ()
         # Coarguments map from V* to V*, i.e. V* -> V*, or equivalently V* x V -> R.
         # So they have one argument in the primal space and one in the dual space.
-        self._arguments = (Argument(self.function_space().dual(), 0), self)
-        self._coefficients = ()
+        # However, when they are composed with linear forms with dual arguments, such as BaseFormOperators,
+        # the primal argument is discarded when analysing the argument as Coarguments.
+        if not outer_form:
+            self._arguments = (Argument(self.function_space().dual(), 0), self)
+        else:
+            self._arguments = (self,)
 
     def reconstruct(self, function_space=None,
                     number=None, part=None):

tests/multigrid/test_p_multigrid.py

@@ -326,7 +326,7 @@ def test_p_multigrid_mixed(mat_type):
 
     # Make the Function spanning the nullspace
     c_basis = assemble(TestFunction(Z.sub(1))*dx)
-    f_basis = Function(c_basis.function_space().dual(), val=c_basis.vector())
+    f_basis = Function(c_basis.function_space().dual(), val=c_basis.dat)
 
     basis = VectorSpaceBasis([f_basis])
     basis.orthonormalize()

tests/regression/test_adjoint_operators.py

@@ -316,7 +316,7 @@ def test_interpolate_hessian_linear_expr():
     dJdm = J.block_variable.tlm_value
     assert isinstance(f.block_variable.adj_value, Cofunction)
     assert isinstance(f.block_variable.hessian_value, Cofunction)
-    Hm = f.block_variable.hessian_value.vector().inner(h.vector())
+    Hm = f.block_variable.hessian_value.dat.inner(h.dat)
     # If the new interpolate block has the right hessian, taylor test
     # convergence rate should be as for the unmodified test.
     assert taylor_test(Jhat, g, h, dJdm=dJdm, Hm=Hm) > 2.9
@@ -374,7 +374,7 @@ def test_interpolate_hessian_nonlinear_expr():
     dJdm = J.block_variable.tlm_value
     assert isinstance(f.block_variable.adj_value, Cofunction)
     assert isinstance(f.block_variable.hessian_value, Cofunction)
-    Hm = f.block_variable.hessian_value.vector().inner(h.vector())
+    Hm = f.block_variable.hessian_value.dat.inner(h.dat)
     # If the new interpolate block has the right hessian, taylor test
     # convergence rate should be as for the unmodified test.
     assert taylor_test(Jhat, g, h, dJdm=dJdm, Hm=Hm) > 2.9
@@ -436,7 +436,7 @@ def test_interpolate_hessian_nonlinear_expr_multi():
     dJdm = J.block_variable.tlm_value
     assert isinstance(f.block_variable.adj_value, Cofunction)
     assert isinstance(f.block_variable.hessian_value, Cofunction)
-    Hm = f.block_variable.hessian_value.vector().inner(h.vector())
+    Hm = f.block_variable.hessian_value.dat.inner(h.dat)
     # If the new interpolate block has the right hessian, taylor test
     # convergence rate should be as for the unmodified test.
     assert taylor_test(Jhat, g, h, dJdm=dJdm, Hm=Hm) > 2.9
@@ -640,7 +640,7 @@ def test_supermesh_project_hessian(vector):
     dJdm = J.block_variable.tlm_value
     assert isinstance(source.block_variable.adj_value, Cofunction)
     assert isinstance(source.block_variable.hessian_value, Cofunction)
-    Hm = source.block_variable.hessian_value.vector().inner(h.vector())
+    Hm = source.block_variable.hessian_value.dat.inner(h.dat)
     assert taylor_test(rf, source, h, dJdm=dJdm, Hm=Hm) > 2.9
 
 

tests/regression/test_assemble_baseform.py

@@ -228,7 +228,7 @@ def test_cofunction_riesz_representation(a):
                     M.petscmat.mult(v_vec, res_vec)
         else:
             # l2 mass matrix is identity
-            Mr = Function(V, val=r.vector())
+            Mr = Function(V, val=r.dat)
 
         # Check residual
         for a, b in zip(Mr.subfunctions, c.subfunctions):
@@ -265,7 +265,7 @@ def test_function_riesz_representation(f):
                     M.petscmat.mult(v_vec, res_vec)
         else:
             # l2 mass matrix is identity
-            Mf = Cofunction(V.dual(), val=f.vector())
+            Mf = Cofunction(V.dual(), val=f.dat)
 
         # Check residual
         for a, b in zip(Mf.subfunctions, r.subfunctions):