diff --git a/demos/netgen/netgen_mesh.py.rst b/demos/netgen/netgen_mesh.py.rst index 2c18c1d74c..7fc649654a 100755 --- a/demos/netgen/netgen_mesh.py.rst +++ b/demos/netgen/netgen_mesh.py.rst @@ -382,7 +382,8 @@ We will now show how to solve the Poisson problem on a high-order mesh, of order bc = DirichletBC(V, 0.0, [1]) A = assemble(a, bcs=bc) b = assemble(l) - bc.apply(b) + b_riesz = b.riesz_representation(riesz_map="l2") + bc.apply(b_riesz) solve(A, sol, b, solver_parameters={"ksp_type": "cg", "pc_type": "lu"}) VTKFile("output/Sphere.pvd").write(sol) diff --git a/docs/source/solving-interface.rst b/docs/source/solving-interface.rst index 61000af7b2..2dd547c818 100644 --- a/docs/source/solving-interface.rst +++ b/docs/source/solving-interface.rst @@ -144,8 +144,8 @@ pass in. In the pre-assembled case, we are solving a linear system: Where :math:`A` is a known matrix, :math:`\vec{b}` is a known right hand side vector and :math:`\vec{x}` is the unknown solution vector. In Firedrake, :math:`A` is represented as a -:py:class:`~.Matrix`, while :math:`\vec{x}` is a :py:class:`~.Function`, and -:math:`\vec{b}` a :py:class:`~.Cofunction`. +:py:class:`~.MatrixBase`, while :math:`\vec{b}` and +:math:`\vec{x}` can be :py:class:`~.Function`\s or :py:class:`~.Cofunction`\s. We build these values by calling ``assemble`` on the UFL forms that define our problem, which, as before are denoted ``a`` and ``L``. Similarly to the linear variational case, we first need a function in diff --git a/firedrake/adjoint_utils/blocks/dirichlet_bc.py b/firedrake/adjoint_utils/blocks/dirichlet_bc.py index 7a4784812c..146bb71a59 100644 --- a/firedrake/adjoint_utils/blocks/dirichlet_bc.py +++ b/firedrake/adjoint_utils/blocks/dirichlet_bc.py @@ -51,7 +51,7 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, adj_output = None for adj_input in adj_inputs: if isconstant(c): - adj_value = firedrake.Function(self.parent_space.dual()) + adj_value = firedrake.Function(self.parent_space) adj_input.apply(adj_value) if self.function_space != self.parent_space: vec = extract_bc_subvector( @@ -87,11 +87,12 @@ def evaluate_adj_component(self, inputs, adj_inputs, block_variable, idx, # you can even use the Function outside its domain. # For now we will just assume the FunctionSpace is the same for # the BC and the Function. - adj_value = firedrake.Function(self.parent_space.dual()) + adj_value = firedrake.Function(self.parent_space) adj_input.apply(adj_value) - r = extract_bc_subvector( + output = extract_bc_subvector( adj_value, c.function_space(), bc ) + r = output.riesz_representation(riesz_map="l2") if adj_output is None: adj_output = r else: diff --git a/firedrake/adjoint_utils/blocks/solving.py b/firedrake/adjoint_utils/blocks/solving.py index 5eb99dccd2..93931d3e8d 100644 --- a/firedrake/adjoint_utils/blocks/solving.py +++ b/firedrake/adjoint_utils/blocks/solving.py @@ -190,8 +190,10 @@ def _assemble_and_solve_adj_eq(self, dFdu_adj_form, dJdu, compute_bdy): kwargs["bcs"] = bcs dFdu = firedrake.assemble(dFdu_adj_form, **kwargs) + dJdu_func = dJdu.riesz_representation(riesz_map="l2") for bc in bcs: - bc.apply(dJdu) + bc.apply(dJdu_func) + dJdu.assign(dJdu_func.riesz_representation(riesz_map="l2")) adj_sol = firedrake.Function(self.function_space) firedrake.solve( @@ -201,7 +203,7 @@ def _assemble_and_solve_adj_eq(self, dFdu_adj_form, dJdu, compute_bdy): adj_sol_bdy = None if compute_bdy: adj_sol_bdy = firedrake.Function( - self.function_space.dual(), + self.function_space, dJdu_copy.dat - firedrake.assemble( firedrake.action(dFdu_adj_form, adj_sol) ).dat @@ -811,7 +813,7 @@ def __init__(self, source, target_space, target, bcs=[], **kwargs): self.add_dependency(bc, no_duplicates=True) def apply_mixedmass(self, a): - b = firedrake.Function(self.target_space) + b = self.backend.Cofunction(self.target_space.dual()) with a.dat.vec_ro as vsrc, b.dat.vec_wo as vrhs: self.mixed_mass.mult(vsrc, vrhs) return b diff --git a/firedrake/assemble.py b/firedrake/assemble.py index 9090377dae..ad22130dee 100644 --- a/firedrake/assemble.py +++ b/firedrake/assemble.py @@ -1201,7 +1201,7 @@ def _check_tensor(self, tensor): rank = len(self._form.arguments()) if rank == 1: test, = self._form.arguments() - if tensor is not None and test.function_space() != tensor.function_space(): + if tensor is not None and test.function_space() != tensor.function_space().dual(): raise ValueError("Form's argument does not match provided result tensor") @staticmethod diff --git a/firedrake/functionspaceimpl.py b/firedrake/functionspaceimpl.py index a2ece362d8..a9d4293d8c 100644 --- a/firedrake/functionspaceimpl.py +++ b/firedrake/functionspaceimpl.py @@ -297,7 +297,8 @@ def restore_work_function(self, function): def __eq__(self, other): try: - return self.topological == other.topological and \ + return type(self) == type(other) and \ + self.topological == other.topological and \ self.mesh() is other.mesh() except AttributeError: return False diff --git a/firedrake/linear_solver.py b/firedrake/linear_solver.py index c1dfbcc07e..7a9f7d807d 100644 --- a/firedrake/linear_solver.py +++ b/firedrake/linear_solver.py @@ -147,6 +147,11 @@ def solve(self, x, b): if not isinstance(b, (function.Function, cofunction.Cofunction)): raise TypeError("Provided RHS is a '%s', not a Function or Cofunction" % type(b).__name__) + if x.function_space() != self.trial_space or b.function_space() != self.test_space.dual(): + # When solving `Ax = b`, with A: V x U -> R, or equivalently A: V -> U*, + # we need to make sure that x and b belong to V and U*, respectively. + raise ValueError("Mismatching function spaces.") + if len(self.trial_space) > 1 and self.nullspace is not None: self.nullspace._apply(self.trial_space.dof_dset.field_ises) if len(self.test_space) > 1 and self.transpose_nullspace is not None: diff --git a/tests/regression/test_netgen.py b/tests/regression/test_netgen.py index bc92cb5b16..461c272ea9 100644 --- a/tests/regression/test_netgen.py +++ b/tests/regression/test_netgen.py @@ -52,7 +52,8 @@ def poisson(h, degree=2): # Assembling matrix A = assemble(a, bcs=bc) b = assemble(l) - bc.apply(b) + b_riesz = b.riesz_representation(riesz_map="l2") + bc.apply(b_riesz) # Solving the problem solve(A, u, b, solver_parameters={"ksp_type": "preonly", "pc_type": "lu"}) @@ -96,7 +97,8 @@ def poisson3D(h, degree=2): # Assembling matrix A = assemble(a, bcs=bc) b = assemble(l) - bc.apply(b) + b_riesz = b.riesz_representation(riesz_map="l2") + bc.apply(b_riesz) # Solving the problem solve(A, u, b, solver_parameters={"ksp_type": "preonly", "pc_type": "lu"}) diff --git a/tests/regression/test_solving_interface.py b/tests/regression/test_solving_interface.py index 0493de99d5..dfbe4823ad 100644 --- a/tests/regression/test_solving_interface.py +++ b/tests/regression/test_solving_interface.py @@ -238,3 +238,24 @@ def test_solve_cofunction_rhs(): Aw = assemble(action(a, w)) assert isinstance(Aw, Cofunction) assert np.allclose(Aw.dat.data_ro, L.dat.data_ro) + + +def test_linear_solver_check_spaces(): + mesh = UnitSquareMesh(10, 10) + V = FunctionSpace(mesh, "CG", 1) + + u = TrialFunction(V) + v = TestFunction(V) + a = inner(u, v) * dx + A = assemble(a) + + L = Cofunction(V.dual()) + L.vector()[:] = 1. + + Lf = L.riesz_representation(riesz_map="l2") + + w = Function(V) + solve(A, w, L) + + with pytest.raises(ValueError): + solve(A, w, Lf)