set some sensible default mass solve parameters

gusto/solvers/__init__.py

@@ -1,2 +1,3 @@
 from gusto.solvers.linear_solvers import *         # noqa
-from gusto.solvers.preconditioners import *        # noqa
+from gusto.solvers.preconditioners import *        # noqa
+from gusto.solvers.parameters import *             # noqa

gusto/solvers/parameters.py

@@ -0,0 +1,62 @@
+"""
+This module provides some parameters sets that are good defaults
+for particular kinds of system.
+"""
+from gusto.time_discretisation.wrappers import is_cg
+
+__all__ = ['mass_parameters']
+
+
+def mass_parameters(V, spaces=None, ignore_vertical=True):
+    """
+    PETSc solver parameters for mass matrices.
+
+    Args:
+        spaces: Optional `Spaces` object. If present, any subspace
+            of V that came from the `Spaces` object will use the
+            continuity information from `spaces`.
+            If not present, continuity is checked with `is_cg`.
+
+        ignore_vertical: whether to include the vertical direction when checking
+            field continuity on extruded meshes. If True, only the horizontal
+            continuity will be considered, e.g. the standard theta space will
+            be treated as discontinuous.
+    """
+
+    extruded = hasattr(V.mesh, "_base_mesh")
+
+    continuous_fields = set()
+    for i, Vsub in V.subfunctions:
+        field = Vsub.name or str(i)
+
+        if spaces is not None:
+            continuous = spaces.continuity.get(field, is_cg(Vsub))
+        else:
+            continuous = is_cg(Vsub)
+
+        # For extruded meshes the continuity is recorded
+        # separately for the horizontal and vertical directions.
+        if extruded:
+            if ignore_vertical:
+                continuous = continuous['horizontal']
+            else:
+                continuous = (continuous['horizontal']
+                              and continuous['vertical'])
+
+        if continuous:
+            continuous_fields.add(field)
+
+    parameters = {
+        'ksp_type': 'preonly',
+        'pc_type': 'fieldsplit',
+        'pc_fieldsplit_type': 'additive',
+        'pc_fieldsplit_0_fields': ','.join(continuous_fields),
+        'fieldsplit': {
+            'ksp_type': 'preonly',
+            'pc_type': 'bjacobi',
+            'sub_pc_type': 'ilu'
+        },
+        'fieldsplit_0_ksp_type': 'cg',
+    }
+
+    return parameters

gusto/time_discretisation/time_discretisation.py

@@ -21,6 +21,7 @@
                                mass_weighted, nonlinear_time_derivative)
 from gusto.core.logging import logger, DEBUG, logging_ksp_monitor_true_residual
 from gusto.time_discretisation.wrappers import *
+from gusto.solvers import mass_parameters
 
 __all__ = ["TimeDiscretisation", "ExplicitTimeDiscretisation", "BackwardEuler",
            "ThetaMethod", "TrapeziumRule", "TR_BDF2"]
@@ -433,10 +434,9 @@ def __init__(self, domain, field_name=None, subcycling_options=None,
 
         # get default solver options if none passed in
         if solver_parameters is None:
-            self.solver_parameters = {'snes_type': 'ksponly',
-                                      'ksp_type': 'cg',
-                                      'pc_type': 'bjacobi',
-                                      'sub_pc_type': 'ilu'}
+            self.solver_parameters = mass_parameters(
+                equation.function_space, equation.domain.spaces)
+            self.solver_parameters['snes_type'] = 'ksponly'
         else:
             self.solver_parameters = solver_parameters
 
@@ -479,41 +479,12 @@ def setup(self, equation, apply_bcs=True, *active_labels):
                        + ' as the time derivative term is nonlinear')
             logger.warning(message)
             self.solver_parameters['snes_type'] = 'newtonls'
-
-            # TODO: move this outside the if-else and just set the fieldsplit_type here
-            if len(self.x0.subfunctions) > 1:
-                self.solver_parameters.update({
-                    'ksp_type': 'preonly',
-                    'pc_type': 'fieldsplit',
-                    'pc_fieldsplit_type': 'multiplicative',
-                    'fieldsplit': {
-                        'ksp_type': 'preonly',
-                        'pc_type': 'bjacobi',
-                        'sub_pc_type': 'ilu'
-                    },
-                    'fieldsplit_HDiv_ksp_type': 'cg'
-                })
         else:
             self.solver_parameters.setdefault('snes_lag_jacobian', -2)
             self.solver_parameters.setdefault('snes_lag_jacobian_persists', None)
             self.solver_parameters.setdefault('snes_lag_preconditioner', -2)
             self.solver_parameters.setdefault('snes_lag_preconditioner_persists', None)
 
-            # The time derivatives for each field are independent if all time
-            # time derivatives are linear, so we solve each one independently.
-            if len(self.x0.subfunctions) > 1:
-                self.solver_parameters.update({
-                    'ksp_type': 'preonly',
-                    'pc_type': 'fieldsplit',
-                    'pc_fieldsplit_type': 'additive',
-                    'fieldsplit': {
-                        'ksp_type': 'preonly',
-                        'pc_type': 'bjacobi',
-                        'sub_pc_type': 'ilu'
-                    },
-                    'fieldsplit_HDiv_ksp_type': 'cg'
-                })
-
     @cached_property
     def lhs(self):
         """Set up the discretisation's left hand side (the time derivative)."""

gusto/timestepping/semi_implicit_quasi_newton.py

@@ -13,7 +13,7 @@
 from gusto.core.labels import (transport, diffusion, time_derivative,
                                linearisation, prognostic, hydrostatic,
                                physics_label, sponge, incompressible)
-from gusto.solvers import LinearTimesteppingSolver
+from gusto.solvers import LinearTimesteppingSolver, mass_parameters
 from gusto.core.logging import logger, DEBUG, logging_ksp_monitor_true_residual
 from gusto.time_discretisation.time_discretisation import ExplicitTimeDiscretisation
 from gusto.timestepping.timestepper import BaseTimestepper
@@ -571,17 +571,7 @@ def __init__(self, equation, alpha):
             constant_jacobian=True
         )
 
-        self.solver_parameters = {
-            'ksp_type': 'preonly',
-            'pc_type': 'fieldsplit',
-            'pc_fieldsplit_type': 'additive',
-            'fieldsplit': {
-                'ksp_type': 'preonly',
-                'pc_type': 'bjacobi',
-                'sub_pc_type': 'ilu'
-            },
-            'fieldsplit_HDiv_ksp_type': 'cg'
-        }
+        self.solver_parameters = mass_parameters(W, equation.domain.spaces)
 
         self.solvers = {}
         self.solvers["explicit"] = LinearVariationalSolver(