Fix failing Skamarock Klemp test and add PointDataOutput warning/error

firedrakeproject · Aug 2, 2023 · 4882454 · 4882454
1 parent d0042e4
commit 4882454
Show file tree

Hide file tree

Showing 2 changed files with 31 additions and 7 deletions.
diff --git a/examples/compressible/skamarock_klemp_nonhydrostatic.py b/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,12 +51,30 @@
 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)],
-                          log_level='INFO')
+
+# 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)],
+        log_level='INFO'
+    )
+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'],
+        log_level='INFO'
+    )
+
 diagnostic_fields = [CourantNumber(), Gradient("u"), Perturbation('theta'),
                      Gradient("theta_perturbation"), Perturbation('rho'),
                      RichardsonNumber("theta", parameters.g/Tsurf), Gradient("theta")]

diff --git a/gusto/io.py b/gusto/io.py
@@ -15,6 +15,10 @@
 __all__ = ["pick_up_mesh", "IO"]
 
 
+class GustoIOError(IOError):
+    pass
+
+
 def pick_up_mesh(output, mesh_name):
     """
     Picks up a checkpointed mesh. This must be the first step of any model being
@@ -68,6 +72,8 @@ def __init__(self, filename, field_points, description,
         self.field_points = field_points
         self.tolerance = tolerance
         self.comm = comm
+        if self.comm.size > 1:
+            raise GustoIOError("PointDataOutput does not work in parallel")
         if not create:
             return
         if self.comm.rank == 0: