Merge branch 'main' into conservative_projection

.github/workflows/build.yml

@@ -26,23 +26,37 @@ jobs:
         run: |
           cd ..
           rm -rf build
-      - name: Setup python
-        uses: actions/setup-python@v4
-        with:
-          python-version: 3.8
       - name: Install Gusto
         run: |
           . /home/firedrake/firedrake/bin/activate
           python -m pip install -r requirements.txt
           python -m pip install -e .
           python -m pip install  \
-            pytest-cov pytest-timeout pytest-xdist
+            pytest-timeout pytest-xdist
       - name: Gusto tests
         run: |
           . /home/firedrake/firedrake/bin/activate
-          which firedrake-clean
+          firedrake-clean
+          export GUSTO_PARALLEL_LOG=FILE
           python -m pytest \
             -n 12 --dist worksteal \
             --durations=100 \
-            --cov gusto \
+            --timeout=3600 \
+            --timeout-method=thread \
+            -o faulthandler_timeout=3660 \
             -v unit-tests integration-tests examples
+        timeout-minutes: 120
+      - name: Prepare logs
+        if: always()
+        run: |
+          mkdir logs
+          cd /tmp/pytest-of-firedrake/pytest-0/
+          find . -name "*.log" -exec cp --parents {} /__w/gusto/gusto/logs/ \;
+      - name: Upload artifact
+        if: always()
+        uses: actions/upload-pages-artifact@v1
+        with:
+          name: log-files
+          path: /__w/gusto/gusto/logs
+          retention-days: 5
+

README.md

@@ -0,0 +1,47 @@
+# Gusto
+
+Gusto is a Python code library, providing a toolkit of finite element methods for modelling geophysical fluids, such as the atmosphere and ocean.
+The methods used by Gusto are underpinned by the [Firedrake](http://firedrakeproject.org) finite element code generation software.
+
+Gusto is particularly targeted at the numerical methods used by the dynamical cores used in numerical weather prediction and climate models.
+Gusto focuses on compatible finite element discretisations, in which variables lie in function spaces that preserve the underlying geometric structure of the equations.
+These compatible methods underpin the Met Office's next-generation model, [LFRic](https://www.metoffice.gov.uk/research/approach/modelling-systems/lfric).
+
+### Gusto is designed to provide:
+- a testbed for **rapid prototyping** of novel numerical methods
+- a **flexible framework** for exploring different modelling choices for geophysical fluid dynamics
+- a **simple environment** for setting up and running test cases
+
+## Download
+
+The best way to install Gusto is as an additional package when installing [Firedrake](http://firedrakeproject.org). Usually, for a Mac with Homebrew or an Ubuntu installation this is done by downloading the Firedrake install script and executing it:
+```
+curl -0 https://raw.githubusercontent/com/firedrakeproject/firedrake/master/scripts/firedrake-install
+python3 firedrake-install --install gusto
+```
+For an up-to-date installation guide, see the [firedrake installation instructions](http://firedrakeproject.org/download.html). Once installed, Gusto must be run from within the Firedrake virtual environment, which is activated via
+```
+source firedrake/bin/activate
+```
+## Getting Started
+
+To test your Gusto installation, run the test-suites:
+```
+cd firedrake/src/gusto
+make test
+```
+
+The `examples` directory contains several test cases, which you can play with to get started with Gusto.
+You can also see the [gusto case studies repository](https://github.com/firedrakeproject/gusto_case_studies), which contains a larger collection of test cases that use Gusto.
+
+Gusto is documented [here](https://www.firedrakeproject.org/gusto-docs/), which is generated from the doc-strings in the codebase.
+
+## Visualisation
+
+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.
+
+## 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/).

docs/Gemfile.lock

@@ -222,7 +222,7 @@ GEM
     rb-fsevent (0.11.2)
     rb-inotify (0.10.1)
       ffi (~> 1.0)
-    rexml (3.3.2)
+    rexml (3.3.6)
       strscan
     rouge (3.26.0)
     ruby2_keywords (0.0.5)

examples/compressible/dry_bryan_fritsch.py

@@ -9,10 +9,8 @@
 from firedrake import (IntervalMesh, ExtrudedMesh,
                        SpatialCoordinate, conditional, cos, pi, sqrt,
                        TestFunction, dx, TrialFunction, Constant, Function,
-                       LinearVariationalProblem, LinearVariationalSolver,
-                       FunctionSpace, VectorFunctionSpace)
+                       LinearVariationalProblem, LinearVariationalSolver)
 import sys
-
 # ---------------------------------------------------------------------------- #
 # Test case parameters
 # ---------------------------------------------------------------------------- #
@@ -61,19 +59,14 @@
 io = IO(domain, output, diagnostic_fields=diagnostic_fields)
 
 # Transport schemes -- set up options for using recovery wrapper
-VDG1 = domain.spaces("DG1_equispaced")
-VCG1 = FunctionSpace(mesh, "CG", 1)
-Vu_DG1 = VectorFunctionSpace(mesh, VDG1.ufl_element())
-Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1)
-
-u_opts = RecoveryOptions(embedding_space=Vu_DG1,
-                         recovered_space=Vu_CG1,
-                         boundary_method=BoundaryMethod.taylor)
-rho_opts = RecoveryOptions(embedding_space=VDG1,
-                           recovered_space=VCG1,
-                           boundary_method=BoundaryMethod.taylor)
-theta_opts = RecoveryOptions(embedding_space=VDG1,
-                             recovered_space=VCG1)
+boundary_methods = {'DG': BoundaryMethod.taylor,
+                    'HDiv': BoundaryMethod.taylor}
+
+recovery_spaces = RecoverySpaces(domain, boundary_methods, use_vector_spaces=True)
+
+u_opts = recovery_spaces.HDiv_options
+rho_opts = recovery_spaces.DG_options
+theta_opts = recovery_spaces.theta_options
 
 transported_fields = [SSPRK3(domain, "rho", options=rho_opts),
                       SSPRK3(domain, "theta", options=theta_opts),

examples/compressible/mountain_hydrostatic.py

@@ -142,34 +142,70 @@
 exner = Function(Vr)
 rho_b = Function(Vr)
 
-exner_params = {'ksp_type': 'gmres',
-                'ksp_monitor_true_residual': None,
-                'pc_type': 'python',
-                'mat_type': 'matfree',
-                'pc_python_type': 'gusto.VerticalHybridizationPC',
-                # Vertical trace system is only coupled vertically in columns
-                # block ILU is a direct solver!
-                'vert_hybridization': {'ksp_type': 'preonly',
-                                       'pc_type': 'bjacobi',
-                                       'sub_pc_type': 'ilu'}}
-
-compressible_hydrostatic_balance(eqns, theta_b, rho_b, exner,
-                                 top=True, exner_boundary=0.5,
-                                 params=exner_params)
-
-# Use kernel as a parallel-safe method of computing minimum
-min_kernel = kernels.MinKernel()
-
-p0 = min_kernel.apply(exner)
-compressible_hydrostatic_balance(eqns, theta_b, rho_b, exner,
-                                 top=True, params=exner_params)
-p1 = min_kernel.apply(exner)
-alpha = 2.*(p1-p0)
-beta = p1-alpha
-exner_top = (1.-beta)/alpha
-compressible_hydrostatic_balance(eqns, theta_b, rho_b, exner,
-                                 top=True, exner_boundary=exner_top, solve_for_rho=True,
-                                 params=exner_params)
+exner_surf = 1.0         # maximum value of Exner pressure at surface
+max_iterations = 10      # maximum number of hydrostatic balance iterations
+tolerance = 1e-7         # tolerance for hydrostatic balance iteration
+
+# Set up kernels to evaluate global minima and maxima of fields
+min_kernel = MinKernel()
+max_kernel = MaxKernel()
+
+# First solve hydrostatic balance that gives Exner = 1 at bottom boundary
+# This gives us a guess for the top boundary condition
+bottom_boundary = Constant(exner_surf, domain=mesh)
+logger.info(f'Solving hydrostatic with bottom Exner of {exner_surf}')
+compressible_hydrostatic_balance(
+    eqns, theta_b, rho_b, exner, top=False, exner_boundary=bottom_boundary
+)
+
+# Solve hydrostatic balance again, but now use minimum value from first
+# solve as the *top* boundary condition for Exner
+top_value = min_kernel.apply(exner)
+top_boundary = Constant(top_value, domain=mesh)
+logger.info(f'Solving hydrostatic with top Exner of {top_value}')
+compressible_hydrostatic_balance(
+    eqns, theta_b, rho_b, exner, top=True, exner_boundary=top_boundary
+)
+
+max_bottom_value = max_kernel.apply(exner)
+
+# Now we iterate, adjusting the top boundary condition, until this gives
+# a maximum value of 1.0 at the surface
+lower_top_guess = 0.9*top_value
+upper_top_guess = 1.2*top_value
+for i in range(max_iterations):
+    # If max bottom Exner value is equal to desired value, stop iteration
+    if abs(max_bottom_value - exner_surf) < tolerance:
+        break
+
+    # Make new guess by average of previous guesses
+    top_guess = 0.5*(lower_top_guess + upper_top_guess)
+    top_boundary.assign(top_guess)
+
+    logger.info(
+        f'Solving hydrostatic balance iteration {i}, with top Exner value '
+        + f'of {top_guess}'
+    )
+
+    compressible_hydrostatic_balance(
+        eqns, theta_b, rho_b, exner, top=True, exner_boundary=top_boundary
+    )
+
+    max_bottom_value = max_kernel.apply(exner)
+
+    # Adjust guesses based on new value
+    if max_bottom_value < exner_surf:
+        lower_top_guess = top_guess
+    else:
+        upper_top_guess = top_guess
+
+logger.info(f'Final max bottom Exner value of {max_bottom_value}')
+
+# Perform a final solve to obtain hydrostatically balanced rho
+compressible_hydrostatic_balance(
+    eqns, theta_b, rho_b, exner, top=True, exner_boundary=top_boundary,
+    solve_for_rho=True
+)
 
 theta0.assign(theta_b)
 rho0.assign(rho_b)

examples/compressible/skamarock_klemp_hydrostatic.py

@@ -39,10 +39,9 @@
 domain = Domain(mesh, dt, "RTCF", 1)
 
 # Equation
-parameters = CompressibleParameters()
-Omega = as_vector((0., 0., 0.5e-4))
+parameters = CompressibleParameters(Omega=0.5e-4)
 balanced_pg = as_vector((0., -1.0e-4*20, 0.))
-eqns = CompressibleEulerEquations(domain, parameters, Omega=Omega,
+eqns = CompressibleEulerEquations(domain, parameters,
                                   extra_terms=[("u", balanced_pg)])
 
 # I/O

examples/compressible/unsaturated_bubble.py

@@ -6,11 +6,12 @@
 Limiters are applied to the transport of the water species.
 """
 from gusto import *
+from gusto.equations import thermodynamics
 from firedrake import (PeriodicIntervalMesh, ExtrudedMesh,
                        SpatialCoordinate, conditional, cos, pi, sqrt, exp,
                        TestFunction, dx, TrialFunction, Constant, Function,
                        LinearVariationalProblem, LinearVariationalSolver,
-                       FunctionSpace, VectorFunctionSpace, errornorm)
+                       errornorm)
 from firedrake.slope_limiter.vertex_based_limiter import VertexBasedLimiter
 import sys
 
@@ -61,21 +62,17 @@
 diagnostic_fields = [RelativeHumidity(eqns), Perturbation('theta'),
                      Perturbation('water_vapour'), Perturbation('rho'), Perturbation('RelativeHumidity')]
 io = IO(domain, output, diagnostic_fields=diagnostic_fields)
-
 # Transport schemes -- specify options for using recovery wrapper
+boundary_methods = {'DG': BoundaryMethod.taylor,
+                    'HDiv': BoundaryMethod.taylor}
+
+recovery_spaces = RecoverySpaces(domain, boundary_method=boundary_methods, use_vector_spaces=True)
+
+u_opts = recovery_spaces.HDiv_options
+rho_opts = recovery_spaces.DG_options
+theta_opts = recovery_spaces.theta_options
+
 VDG1 = domain.spaces("DG1_equispaced")
-VCG1 = FunctionSpace(mesh, "CG", 1)
-Vu_DG1 = VectorFunctionSpace(mesh, VDG1.ufl_element())
-Vu_CG1 = VectorFunctionSpace(mesh, "CG", 1)
-Vt = domain.spaces("theta")
-
-u_opts = RecoveryOptions(embedding_space=Vu_DG1,
-                         recovered_space=Vu_CG1,
-                         boundary_method=BoundaryMethod.taylor)
-rho_opts = RecoveryOptions(embedding_space=VDG1,
-                           recovered_space=VCG1,
-                           boundary_method=BoundaryMethod.taylor)
-theta_opts = RecoveryOptions(embedding_space=VDG1, recovered_space=VCG1)
 limiter = VertexBasedLimiter(VDG1)
 
 transported_fields = [SSPRK3(domain, "u", options=u_opts),
@@ -92,6 +89,7 @@
 
 # Physics schemes
 # NB: can't yet use wrapper or limiter options with physics
+Vt = domain.spaces('theta')
 rainfall_method = DGUpwind(eqns, 'rain', outflow=True)
 zero_limiter = MixedFSLimiter(eqns, {'water_vapour': ZeroLimiter(Vt),
                                      'cloud_water': ZeroLimiter(Vt)})

examples/shallow_water/shallow_water_1d.py

@@ -1,4 +1,5 @@
 import numpy as np
+import sys
 
 from firedrake import *
 from gusto import *
@@ -9,6 +10,10 @@
 delta = L/n
 mesh = PeriodicIntervalMesh(128, L)
 dt = 0.0001
+if '--running-tests' in sys.argv:
+    T = 0.0005
+else:
+    T = 1
 
 domain = Domain(mesh, dt, 'CG', 1)
 
@@ -61,4 +66,4 @@
 
 D += parameters.H
 
-stepper.run(0, 1)
+stepper.run(0, T)

examples/test_examples_run.py

@@ -3,6 +3,7 @@
 import subprocess
 import glob
 import sys
+import os
 
 
 examples_dir = abspath(dirname(__file__))
@@ -19,4 +20,19 @@ def test_example_runs(example_file, tmpdir, monkeypatch):
     # This ensures that the test writes output in a temporary
     # directory, rather than where pytest was run from.
     monkeypatch.chdir(tmpdir)
-    subprocess.check_call([sys.executable, example_file, "--running-tests"])
+    subprocess.run(
+        [sys.executable, example_file, "--running-tests"],
+        check=True,
+        env=os.environ | {"PYOP2_CFLAGS": "-O0"}
+    )
+
+
+def test_example_runs_parallel(example_file, tmpdir, monkeypatch):
+    # This ensures that the test writes output in a temporary
+    # directory, rather than where pytest was run from.
+    monkeypatch.chdir(tmpdir)
+    subprocess.run(
+        ["mpiexec", "-n", "4", sys.executable, example_file, "--running-tests"],
+        check=True,
+        env=os.environ | {"PYOP2_CFLAGS": "-O0"}
+    )

gusto/core/configuration.py

@@ -115,6 +115,7 @@ class BoussinesqParameters(Configuration):
     g = 9.810616
     N = 0.01  # Brunt-Vaisala frequency (1/s)
     cs = 340  # sound speed (for compressible case) (m/s)
+    Omega = None
 
 
 class CompressibleParameters(Configuration):
@@ -137,6 +138,7 @@ class CompressibleParameters(Configuration):
     w_sat2 = -17.27  # second const. in Teten's formula (no units)
     w_sat3 = 35.86  # third const. in Teten's formula (K)
     w_sat4 = 610.9  # fourth const. in Teten's formula (Pa)
+    Omega = None    # Rotation rate
 
 
 class ShallowWaterParameters(Configuration):