From 2155a8081a10b70e0cce6a20ac5b8711f974df98 Mon Sep 17 00:00:00 2001 From: Josh Hope-Collins Date: Wed, 18 Dec 2024 14:18:38 +0000 Subject: [PATCH] set some sensible default mass solve parameters --- gusto/solvers/__init__.py | 3 +- gusto/solvers/parameters.py | 62 +++++++++++++++++++ .../time_discretisation.py | 37 ++--------- .../semi_implicit_quasi_newton.py | 14 +---- 4 files changed, 70 insertions(+), 46 deletions(-) create mode 100644 gusto/solvers/parameters.py diff --git a/gusto/solvers/__init__.py b/gusto/solvers/__init__.py index efb4a5af4..0d3daf222 100644 --- a/gusto/solvers/__init__.py +++ b/gusto/solvers/__init__.py @@ -1,2 +1,3 @@ from gusto.solvers.linear_solvers import * # noqa -from gusto.solvers.preconditioners import * # noqa \ No newline at end of file +from gusto.solvers.preconditioners import * # noqa +from gusto.solvers.parameters import * # noqa diff --git a/gusto/solvers/parameters.py b/gusto/solvers/parameters.py new file mode 100644 index 000000000..9115c9288 --- /dev/null +++ b/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 diff --git a/gusto/time_discretisation/time_discretisation.py b/gusto/time_discretisation/time_discretisation.py index 6bd275507..808e74faf 100644 --- a/gusto/time_discretisation/time_discretisation.py +++ b/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).""" diff --git a/gusto/timestepping/semi_implicit_quasi_newton.py b/gusto/timestepping/semi_implicit_quasi_newton.py index f0f67a154..c5f8a7ece 100644 --- a/gusto/timestepping/semi_implicit_quasi_newton.py +++ b/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(