Add in riesz map and remove pc update

gusto/solvers/preconditioners.py

@@ -2,7 +2,7 @@
 
 from firedrake import (dot, jump, dS_h, ds_b, ds_t, ds,
                        FacetNormal, Tensor, AssembledVector,
-                       AuxiliaryOperatorPC, PETSc)
+                       AuxiliaryOperatorPC, PETSc, RieszMap)
 
 from firedrake.preconditioners import PCBase
 from firedrake.matrix_free.operators import ImplicitMatrixContext
@@ -548,6 +548,8 @@ def initialize(self, pc):
         self.xrhs = Function(self.W)
         self.y = Function(self.W)
 
+        self.riesz_map = RieszMap(self.W.dual(), solver_parameters=self.cg_ilu_parameters, constant_jacobian=True)
+
         self.y_hybrid = Function(self.W_hyb)
 
         # Split input functions
@@ -753,7 +755,7 @@ def apply(self, pc, x, y):
         with self.xstar.dat.vec_wo as xv:
             x.copy(xv)
 
-        self.xrhs.assign(self.xstar.riesz_representation())
+        self.xrhs.assign(self.riesz_map.__call__(self.xstar))
         # Solve hybridized system
         self.hybridized_solver.solve()
 

gusto/timestepping/semi_implicit_quasi_newton.py

@@ -20,7 +20,7 @@
 from gusto.time_discretisation.time_discretisation import ExplicitTimeDiscretisation
 from gusto.timestepping.timestepper import BaseTimestepper
 from gusto.solvers.solver_presets import hybridised_solver_parameters
 from gusto.equations.compressible_euler_equations import CompressibleEulerEquations
 
 
 __all__ = ["SemiImplicitQuasiNewton", "Forcing", "QuasiNewtonLinearSolver"]
@@ -937,12 +937,6 @@
             self.equation.update_reference_profiles()
             self.solver.invalidate_jacobian()
 
-            # TODO: Issue #686 is to address this reference profile update bug (pythonPC update not called)
-            # this line forces it to update for now
-            pc = self.solver.snes.getKSP().getPC()
-            if (isinstance(self.equation, CompressibleEulerEquations) and pc.getType() == "python"):
-                pc.getPythonContext().update(pc)
-
     def zero_non_prognostics(self, equation, xrhs, field_names, prognostic_names):
         """
         Zero rhs term F(x) for non-prognostics.