merge to main, conflicts fixed

examples/compressible/dcmip_3_1_meanflow_quads.py

@@ -87,8 +87,8 @@
 
 # Transport schemes
 transported_fields = [TrapeziumRule(domain, "u"),
-                      SSPRK3(domain, "rho", subcycles=2),
-                      SSPRK3(domain, "theta", options=SUPGOptions(), subcycles=2)]
+                      SSPRK3(domain, "rho", fixed_subcycles=2),
+                      SSPRK3(domain, "theta", options=SUPGOptions(), fixed_subcycles=2)]
 transport_methods = [DGUpwind(eqns, field) for field in ["u", "rho", "theta"]]
 
 # Linear solver

examples/compressible/mountain_hydrostatic.py

@@ -66,7 +66,7 @@
     dumpfreq=dumpfreq,
     dumplist=['u'],
 )
-diagnostic_fields = [CourantNumber(), VelocityZ(), HydrostaticImbalance(eqns),
+diagnostic_fields = [CourantNumber(), ZComponent('u'), HydrostaticImbalance(eqns),
                      Perturbation('theta'), Perturbation('rho')]
 io = IO(domain, output, diagnostic_fields=diagnostic_fields)
 

examples/compressible/skamarock_klemp_nonhydrostatic.py

@@ -9,7 +9,7 @@
 from gusto import *
 import itertools
 from firedrake import (as_vector, SpatialCoordinate, PeriodicIntervalMesh,
-                       ExtrudedMesh, exp, sin, Function, pi)
+                       ExtrudedMesh, exp, sin, Function, pi, COMM_WORLD)
 import numpy as np
 import sys
 
@@ -51,13 +51,28 @@
 points_z = [H/2.]
 points = np.array([p for p in itertools.product(points_x, points_z)])
 dirname = 'skamarock_klemp_nonlinear'
-output = OutputParameters(
-    dirname=dirname,
-    dumpfreq=dumpfreq,
-    pddumpfreq=dumpfreq,
-    dumplist=['u'],
-    point_data=[('theta_perturbation', points)],
-)
+
+# Dumping point data using legacy PointDataOutput is not supported in parallel
+if COMM_WORLD.size == 1:
+    output = OutputParameters(
+        dirname=dirname,
+        dumpfreq=dumpfreq,
+        pddumpfreq=dumpfreq,
+        dumplist=['u'],
+        point_data=[('theta_perturbation', points)],
+    )
+else:
+    logger.warning(
+        'Dumping point data using legacy PointDataOutput is not'
+        ' supported in parallel\nDisabling PointDataOutput'
+    )
+    output = OutputParameters(
+        dirname=dirname,
+        dumpfreq=dumpfreq,
+        pddumpfreq=dumpfreq,
+        dumplist=['u'],
+    )
+
 diagnostic_fields = [CourantNumber(), Gradient("u"), Perturbation('theta'),
                      Gradient("theta_perturbation"), Perturbation('rho'),
                      RichardsonNumber("theta", parameters.g/Tsurf), Gradient("theta")]

examples/shallow_water/moist_thermal_williamson5.py

@@ -60,7 +60,7 @@
 fexpr = 2*Omega*x[2]/R
 
 # Topography
-phi, lamda = latlon_coords(mesh)
+lamda, phi, _ = lonlatr_from_xyz(x[0], x[1], x[2])
 lsq = (lamda - lamda_c)**2
 thsq = (phi - phi_c)**2
 rsq = min_value(R0sq, lsq+thsq)

examples/shallow_water/thermal_williamson2.py

@@ -1,5 +1,5 @@
 from gusto import *
-from firedrake import (IcosahedralSphereMesh, SpatialCoordinate, sin, cos)
+from firedrake import IcosahedralSphereMesh, SpatialCoordinate, sin, cos
 import sys
 
 # ----------------------------------------------------------------- #
@@ -71,9 +71,9 @@
 D0 = stepper.fields("D")
 b0 = stepper.fields("b")
 
-phi, lamda = latlon_coords(mesh)
+lamda, phi, _ = lonlatr_from_xyz(x[0], x[1], x[2])
 
-uexpr = sphere_to_cartesian(mesh, u_max*cos(phi), 0)
+uexpr = xyz_vector_from_lonlatr(u_max*cos(phi), 0, 0, x)
 g = params.g
 w = Omega*R*u_max + (u_max**2)/2
 sigma = w/10

examples/shallow_water/williamson_2.py

@@ -7,7 +7,7 @@
 """
 
 from gusto import *
-from firedrake import IcosahedralSphereMesh, SpatialCoordinate, as_vector, pi
+from firedrake import IcosahedralSphereMesh, SpatialCoordinate, sin, cos, pi
 import sys
 
 # ---------------------------------------------------------------------------- #
@@ -28,6 +28,7 @@
 # setup shallow water parameters
 R = 6371220.
 H = 5960.
+rotated_pole = (0.0, pi/3)
 
 # setup input that doesn't change with ref level or dt
 parameters = ShallowWaterParameters(H=H)
@@ -42,11 +43,13 @@
     mesh = IcosahedralSphereMesh(radius=R,
                                  refinement_level=ref_level, degree=2)
     x = SpatialCoordinate(mesh)
-    domain = Domain(mesh, dt, 'BDM', 1)
+    domain = Domain(mesh, dt, 'BDM', 1, rotated_pole=rotated_pole)
 
     # Equation
     Omega = parameters.Omega
-    fexpr = 2*Omega*x[2]/R
+    _, lat, _ = rotated_lonlatr_coords(x, rotated_pole)
+    e_lon, _, _ = rotated_lonlatr_vectors(x, rotated_pole)
+    fexpr = 2*Omega*sin(lat)
     eqns = ShallowWaterEquations(domain, parameters, fexpr=fexpr)
 
     # I/O
@@ -63,12 +66,14 @@
                          ShallowWaterKineticEnergy(),
                          ShallowWaterPotentialEnergy(parameters),
                          ShallowWaterPotentialEnstrophy(),
-                         SteadyStateError('u'), SteadyStateError('D')]
+                         SteadyStateError('u'), SteadyStateError('D'),
+                         MeridionalComponent('u', rotated_pole),
+                         ZonalComponent('u', rotated_pole)]
     io = IO(domain, output, diagnostic_fields=diagnostic_fields)
 
     # Transport schemes
     transported_fields = [TrapeziumRule(domain, "u"),
-                          SSPRK3(domain, "D", subcycles=2)]
+                          SSPRK3(domain, "D", fixed_subcycles=2)]
     transport_methods = [DGUpwind(eqns, "u"), DGUpwind(eqns, "D")]
 
     # Time stepper
@@ -82,9 +87,9 @@
     D0 = stepper.fields("D")
     x = SpatialCoordinate(mesh)
     u_max = 2*pi*R/(12*day)  # Maximum amplitude of the zonal wind (m/s)
-    uexpr = as_vector([-u_max*x[1]/R, u_max*x[0]/R, 0.0])
+    uexpr = u_max*cos(lat)*e_lon
     g = parameters.g
-    Dexpr = H - ((R * Omega * u_max + u_max*u_max/2.0)*(x[2]*x[2]/(R*R)))/g
+    Dexpr = H - (R * Omega * u_max + u_max*u_max/2.0)*(sin(lat))**2/g
 
     u0.project(uexpr)
     D0.interpolate(Dexpr)

examples/shallow_water/williamson_5.py

@@ -47,7 +47,7 @@
     # Equation
     Omega = parameters.Omega
     fexpr = 2*Omega*x[2]/R
-    theta, lamda = latlon_coords(mesh)
+    lamda, theta, _ = lonlatr_from_xyz(x[0], x[1], x[2])
     R0 = pi/9.
     R0sq = R0**2
     lamda_c = -pi/2.

gusto/__init__.py

@@ -6,6 +6,7 @@
 from gusto.active_tracers import *                   # noqa
 from gusto.common_forms import *                     # noqa
 from gusto.configuration import *                    # noqa
+from gusto.coord_transforms import *                 # noqa
 from gusto.domain import *                           # noqa
 from gusto.diagnostics import *                      # noqa
 from gusto.diffusion_methods import *                # noqa

gusto/configuration.py

@@ -88,7 +88,8 @@ class OutputParameters(Configuration):
     dumplist = None
     dumplist_latlon = []
     dump_diagnostics = True
-    checkpoint = True
+    diagfreq = 1
+    checkpoint = False
     checkpoint_method = 'checkpointfile'
     checkpoint_pickup_filename = None
     chkptfreq = 1