Merge changes from main into future

.github/workflows/schedule_future.yml

@@ -0,0 +1,13 @@
+name: Weekly future tests
+
+on:
+  schedule:
+    # Scheduled for 0430 UTC on Monday mornings to detect bitrot
+    - cron:  '30 4 * * 1'
+jobs:
+  test_future:
+    uses: ./.github/workflows/test.yml
+    with:
+      source_ref: future
+      firedrake_docker_version: dev-main
+    secrets: inherit

.github/workflows/schedule_main.yml

@@ -0,0 +1,13 @@
+name: Weekly main tests
+
+on:
+  schedule:
+    # Scheduled for 0130 UTC on Monday mornings to detect bitrot
+    - cron:  '30 1 * * 1'
+jobs:
+  test_main:
+    uses: ./.github/workflows/test.yml
+    with:
+      source_ref: main
+      firedrake_docker_version: dev-release
+    secrets: inherit

.github/workflows/test.yml

@@ -87,6 +87,8 @@ jobs:
           pip install --no-binary netCDF4 --no-build-isolation netCDF4
 
       - name: Run tests (nprocs = 1)
+        # Run even if earlier tests failed
+        if: success() || steps.install-two.conclusion == 'success'
         run: |
           . venv-gusto/bin/activate
           : # Use pytest-xdist here so we can have a single collated output (not possible
@@ -95,7 +97,6 @@ jobs:
         timeout-minutes: 60
 
       - name: Run tests (nprocs = 2)
-        # Run even if earlier tests failed
         if: success() || steps.install-two.conclusion == 'success'
         run: |
           . venv-gusto/bin/activate
@@ -120,15 +121,15 @@ jobs:
         uses: actions/upload-artifact@v4
         if: success() || steps.install-two.conclusion == 'success'
         with:
-          name: pytest-logs
+          name: pytest-logs-${{ inputs.firedrake_docker_version }}
           path: pytest_*.log
           retention-days: 5
 
       - name: Upload Gusto log files
         uses: actions/upload-artifact@v4
         if: success() || steps.install-two.conclusion == 'success'
         with:
-          name: gusto-logs
+          name: gusto-logs-${{ inputs.firedrake_docker_version }}
           path: gusto*.log
           retention-days: 5
 

README.md

@@ -44,6 +44,19 @@ Gusto can produce output in two formats:
 - VTU files, which can be viewed with the [Paraview](https://www.paraview.org/) software
 - netCDF files, which has data that can be plotted using standard python packages such as matplotlib. We suggest using the [tomplot](https://github.com/tommbendall/tomplot) Python library, which contains several routines to simplify the plotting of Gusto output.
 
+## Working Practices
+
+Gusto's default development happens on the `main` branch, which is fixed to the latest release of Firedrake.
+We also maintain a `future` branch which builds from the head of Firedrake's `main` branch.
+Contributions to Gusto's `main` branch must pass our Continuous Integration tests.
+
+The relationship between `future` and `main` is as follows:
+- whenever there is a new release of Firedrake, `future` should be merged into `main` *and* `main` should be merged into `future` so that the two are in-sync
+- otherwise, `future` should not be merged into `main`
+- we *may* choose to merge `main` into `future` out-of-cycle of Firedrake releases, which might be desirable after we fix bugs or make significant API changes
+
+We also use some tests with "KGOs" (Known Good Output), which are used to detect changes in science. These can be regenerated using the `integration-tests/data/update_kgos.py` script.
+
 ## Website
 
 For more information, please see our [website](https://www.firedrakeproject.org/gusto/), and please do get in touch via the Gusto channel on the Firedrake project [Slack workspace](https://firedrakeproject.slack.com/).

examples/boussinesq/skamarock_klemp_compressible.py

@@ -15,8 +15,7 @@
 from gusto import (
     Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind,
     TrapeziumRule, SUPGOptions, Divergence, Perturbation, CourantNumber,
-    BoussinesqParameters, BoussinesqEquations, BoussinesqSolver,
-    boussinesq_hydrostatic_balance
+    BoussinesqParameters, BoussinesqEquations, boussinesq_hydrostatic_balance
 )
 
 skamarock_klemp_compressible_bouss_defaults = {
@@ -90,13 +89,9 @@ def skamarock_klemp_compressible_bouss(
         DGUpwind(eqns, "b", ibp=b_opts.ibp)
     ]
 
-    # Linear solver
-    linear_solver = BoussinesqSolver(eqns)
-
     # Time stepper
     stepper = SemiImplicitQuasiNewton(
-        eqns, io, transported_fields, transport_methods,
-        linear_solver=linear_solver
+        eqns, io, transported_fields, transport_methods
     )
 
     # ------------------------------------------------------------------------ #

examples/boussinesq/skamarock_klemp_incompressible.py

@@ -15,8 +15,7 @@
 from gusto import (
     Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind,
     TrapeziumRule, SUPGOptions, Divergence, Perturbation, CourantNumber,
-    BoussinesqParameters, BoussinesqEquations, BoussinesqSolver,
-    boussinesq_hydrostatic_balance
+    BoussinesqParameters, BoussinesqEquations, boussinesq_hydrostatic_balance
 )
 
 skamarock_klemp_incompressible_bouss_defaults = {
@@ -87,13 +86,9 @@ def skamarock_klemp_incompressible_bouss(
         DGUpwind(eqns, "b", ibp=b_opts.ibp)
     ]
 
-    # Linear solver
-    linear_solver = BoussinesqSolver(eqns)
-
-    # Time stepper
+    # Timestepper - pressure p is diagnostic
     stepper = SemiImplicitQuasiNewton(
-        eqns, io, transported_fields, transport_methods,
-        linear_solver=linear_solver
+        eqns, io, transported_fields, transport_methods, solver_prognostics=["u", "b"]
     )
 
     # ------------------------------------------------------------------------ #

examples/compressible_euler/dcmip_3_1_gravity_wave.py

@@ -15,7 +15,7 @@
 from gusto import (
     Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind,
     TrapeziumRule, SUPGOptions, lonlatr_from_xyz, CompressibleParameters,
-    CompressibleEulerEquations, CompressibleSolver, ZonalComponent,
+    CompressibleEulerEquations, ZonalComponent,
     compressible_hydrostatic_balance, Perturbation, GeneralCubedSphereMesh,
     SubcyclingOptions
 )
@@ -105,13 +105,9 @@ def dcmip_3_1_gravity_wave(
         DGUpwind(eqns, field) for field in ["u", "rho", "theta"]
     ]
 
-    # Linear solver
-    linear_solver = CompressibleSolver(eqns)
-
     # Time stepper
     stepper = SemiImplicitQuasiNewton(
-        eqns, io, transported_fields, transport_methods,
-        linear_solver=linear_solver
+        eqns, io, transported_fields, transport_methods
     )
 
     # ------------------------------------------------------------------------ #

examples/compressible_euler/dry_bryan_fritsch.py

@@ -15,8 +15,7 @@
 from gusto import (
     Domain, IO, OutputParameters, SemiImplicitQuasiNewton, SSPRK3, DGUpwind,
     RecoverySpaces, BoundaryMethod, Perturbation, CompressibleParameters,
-    CompressibleEulerEquations, CompressibleSolver,
-    compressible_hydrostatic_balance
+    CompressibleEulerEquations, compressible_hydrostatic_balance
 )
 
 dry_bryan_fritsch_defaults = {
@@ -81,8 +80,7 @@ def dry_bryan_fritsch(
     io = IO(domain, output, diagnostic_fields=diagnostic_fields)
 
     # Transport schemes -- set up options for using recovery wrapper
-    boundary_methods = {'DG': BoundaryMethod.taylor,
-                        'HDiv': BoundaryMethod.taylor}
+    boundary_methods = {'DG': BoundaryMethod.taylor}
 
     recovery_spaces = RecoverySpaces(
         domain, boundary_methods, use_vector_spaces=True
@@ -102,13 +100,10 @@ def dry_bryan_fritsch(
         DGUpwind(eqns, field) for field in ["u", "rho", "theta"]
     ]
 
-    # Linear solver
-    linear_solver = CompressibleSolver(eqns)
-
     # Time stepper
     stepper = SemiImplicitQuasiNewton(
         eqns, io, transported_fields, transport_methods,
-        linear_solver=linear_solver
+        tau_values={'rho': 1.0, 'theta': 1.0}
     )
 
     # ------------------------------------------------------------------------ #

examples/compressible_euler/schaer_mountain.py

@@ -15,7 +15,7 @@
     SpatialCoordinate, exp, pi, cos, Function, Mesh, Constant
 )
 from gusto import (
-    Domain, CompressibleParameters, CompressibleSolver, logger,
+    Domain, CompressibleParameters, logger,
     OutputParameters, IO, SSPRK3, DGUpwind, SemiImplicitQuasiNewton,
     compressible_hydrostatic_balance, SpongeLayerParameters, Exner, ZComponent,
     Perturbation, SUPGOptions, TrapeziumRule, MaxKernel, MinKernel,
@@ -97,11 +97,11 @@ def schaer_mountain(
 
     # Equation
     parameters = CompressibleParameters(mesh, g=g, cp=cp)
-    sponge = SpongeLayerParameters(
+    sponge_params = SpongeLayerParameters(
         mesh, H=domain_height, z_level=domain_height-sponge_depth, mubar=mu_dt/dt
     )
     eqns = CompressibleEulerEquations(
-        domain, parameters, sponge_options=sponge, u_transport_option=u_eqn_type
+        domain, parameters, sponge_options=sponge_params, u_transport_option=u_eqn_type
     )
 
     # I/O
@@ -134,14 +134,11 @@ def schaer_mountain(
         DGUpwind(eqns, "theta", ibp=theta_opts.ibp)
     ]
 
-    # Linear solver
-    tau_values = {'rho': 1.0, 'theta': 1.0}
-    linear_solver = CompressibleSolver(eqns, alpha, tau_values=tau_values)
-
     # Time stepper
+    tau_values = {'rho': 1.0, 'theta': 1.0}
     stepper = SemiImplicitQuasiNewton(
         eqns, io, transported_fields, transport_methods,
-        linear_solver=linear_solver, alpha=alpha, spinup_steps=spinup_steps
+        alpha=alpha, tau_values=tau_values, spinup_steps=spinup_steps
     )
 
     # ------------------------------------------------------------------------ #

examples/compressible_euler/skamarock_klemp_nonhydrostatic.py

@@ -13,21 +13,19 @@
 """
 from argparse import ArgumentParser, ArgumentDefaultsHelpFormatter
 
-from petsc4py import PETSc
-PETSc.Sys.popErrorHandler()
 import itertools
 from firedrake import (
     as_vector, SpatialCoordinate, PeriodicIntervalMesh, ExtrudedMesh, exp, sin,
-    Function, pi, COMM_WORLD, sqrt
+    PETSc, Function, pi, COMM_WORLD, sqrt
 )
 import numpy as np
 from gusto import (
-    Domain, IO, OutputParameters, TRBDF2QuasiNewton, SemiImplicitQuasiNewton, SSPRK3,
-    DGUpwind, logger, SUPGOptions, Perturbation, CompressibleParameters,
+    Domain, IO, OutputParameters, TRBDF2QuasiNewton, SemiImplicitQuasiNewton,
+    DGUpwind, logger, EmbeddedDGOptions, Perturbation, CompressibleParameters,
     CompressibleEulerEquations, HydrostaticCompressibleEulerEquations,
-    compressible_hydrostatic_balance, RungeKuttaFormulation, CompressibleSolver,
-    hydrostatic_parameters, SubcyclingOptions,
+    compressible_hydrostatic_balance, RungeKuttaFormulation, SubcyclingOptions, SSPRK3
 )
+PETSc.Sys.popErrorHandler()
 
 skamarock_klemp_nonhydrostatic_defaults = {
     'ncolumns': 150,
@@ -70,6 +68,11 @@ def skamarock_klemp_nonhydrostatic(
 
     element_order = 1
 
+    if timestepper == 'TR-BDF2':
+        gamma = (1-sqrt(2)/2)
+    else:
+        alpha = 0.5
+
     # ------------------------------------------------------------------------ #
     # Set up model objects
     # ------------------------------------------------------------------------ #
@@ -94,6 +97,11 @@ def skamarock_klemp_nonhydrostatic(
     # Adjust default directory name
     if hydrostatic and dirname == skamarock_klemp_nonhydrostatic_defaults['dirname']:
         dirname = f'hyd_switch_{dirname}'
+    if timestepper == 'TR-BDF2' and (
+        dirname == skamarock_klemp_nonhydrostatic_defaults['dirname']
+        or dirname == f'hyd_switch_{skamarock_klemp_nonhydrostatic_defaults["dirname"]}'
+    ):
+        dirname = f'{dirname}_trbdf2'
 
     # Dumping point data using legacy PointDataOutput is not supported in parallel
     if COMM_WORLD.size == 1:
@@ -116,11 +124,10 @@ def skamarock_klemp_nonhydrostatic(
     io = IO(domain, output, diagnostic_fields=diagnostic_fields)
 
     # Transport schemes
-    theta_opts = SUPGOptions()
-    if timestepper == 'SIQN':
-        subcycling_options = SubcyclingOptions(subcycle_by_courant=0.25)
-    else:
-        subcycling_options = None
+    # We include subcycling here for test coverage,
+    # it is not necessary to subcycle for this test case!
+    subcycling_options = SubcyclingOptions(subcycle_by_courant=0.25)
+    theta_opts = EmbeddedDGOptions()
     transported_fields = [
         SSPRK3(domain, "u", subcycling_options=subcycling_options),
         SSPRK3(
@@ -135,38 +142,26 @@ def skamarock_klemp_nonhydrostatic(
     transport_methods = [
         DGUpwind(eqns, "u"),
         DGUpwind(eqns, "rho", advective_then_flux=True),
-        DGUpwind(eqns, "theta", ibp=theta_opts.ibp)
+        DGUpwind(eqns, "theta")
     ]
 
     # Linear solver
-    if hydrostatic:
-        if timestepper == 'TR-BDF2':
-            raise ValueError('Hydrostatic equations not implmented for TR-BDF2')
-
-        linear_solver = CompressibleSolver(
-            eqns, solver_parameters=hydrostatic_parameters,
-            overwrite_solver_parameters=True
-        )
-    else:
-        linear_solver = CompressibleSolver(eqns)
+    # The use of advective-then-flux formulation and 2x2 Quasi-Newton iterations
+    # requires tau values to take implicit values for rho and theta
+    tau_values = {'rho': 1.0, 'theta': 1.0}
+    if hydrostatic and timestepper == 'TR-BDF2':
+        raise ValueError('Hydrostatic equations not implmented for TR-BDF2')
 
     # Time stepper
     if timestepper == 'TR-BDF2':
-        gamma = (1-sqrt(2)/2)
-        gamma2 = (1 - 2*float(gamma))/(2 - 2*float(gamma))
-        tr_solver = CompressibleSolver(eqns, alpha=gamma)
-        bdf_solver = CompressibleSolver(eqns, alpha=gamma2)
         stepper = TRBDF2QuasiNewton(
             eqns, io, transported_fields, transport_methods,
-            gamma=gamma,
-            tr_solver=tr_solver,
-            bdf_solver=bdf_solver
+            gamma=gamma, tau_values_tr=tau_values, tau_values_bdf=tau_values
         )
 
     elif timestepper == 'SIQN':
         stepper = SemiImplicitQuasiNewton(
-            eqns, io, transported_fields, transport_methods,
-            linear_solver=linear_solver
+            eqns, io, transported_fields, transport_methods, alpha=alpha, tau_values=tau_values
         )
     # ------------------------------------------------------------------------ #
     # Initial conditions
@@ -269,7 +264,7 @@ def skamarock_klemp_nonhydrostatic(
         default=skamarock_klemp_nonhydrostatic_defaults['hydrostatic']
     )
     parser.add_argument(
-        'timestepper',
+        '--timestepper',
         help='Which time stepper to use, takes SIQN or TR-BDF2',
         type=str,
         choices=['SIQN', 'TR-BDF2'],