schaer mountain working and settings

examples/compressible_euler/schaer_mountain.py

@@ -57,8 +57,8 @@ def schaer_mountain(
     cp = 1004.5              # specific heat capacity at constant pressure
     mu_dt = 1.2              # strength of sponge layer, no units
     exner_surf = 1.0         # maximum value of Exner pressure at surface
-    max_iterations = 100     # maximum number of hydrostatic balance iterations
-    tolerance = 1e-12        # tolerance for hydrostatic balance iteration
+    max_iterations = 20      # maximum number of hydrostatic balance iterations
+    tolerance = 1e-8         # tolerance for hydrostatic balance iteration
 
     # ------------------------------------------------------------------------ #
     # Our settings for this set up

gusto/solvers/parameters.py

@@ -121,4 +121,4 @@ def mass_parameters(V, spaces=None, ignore_vertical=True):
             'sub_pc_type': 'ilu'
         }
     }
-}
+}

gusto/time_discretisation/time_discretisation.py

@@ -419,11 +419,12 @@ def update_subcycling(self):
 
             # Cap number of subcycles
             if max_subcycles is not None:
-                self.ncycles = min(self.ncycles, max_subcycles)
-                logger.warning(
-                    'Adaptive subcycling: capping number of subcycles at '
-                    f'{max_subcycles}'
-                )
+                if self.ncycles > max_subcycles:
+                    logger.warning(
+                        'Adaptive subcycling: capping number of subcycles at '
+                        f'{max_subcycles}'
+                    )
+                    self.ncycles = max_subcycles
 
             logger.debug(f'Performing {self.ncycles} subcycles')
             self.dt.assign(self.original_dt/self.ncycles)

plotting/compressible_euler/plot_schaer_mountain.py

@@ -8,14 +8,13 @@
 import matplotlib.pyplot as plt
 from netCDF4 import Dataset
 import numpy as np
-import pandas as pd
 from tomplot import (
     set_tomplot_style, tomplot_cmap, plot_contoured_field,
     add_colorbar_ax, tomplot_field_title, tomplot_contours,
     extract_gusto_coords, extract_gusto_field, reshape_gusto_data
 )
 
-test = 'schaer_alpha_0p51'
+test = 'schaer_mountain'
 
 # ---------------------------------------------------------------------------- #
 # Directory for results and plots
@@ -35,8 +34,8 @@
     r'$w$ (m s$^{-1}$)', r'$\Delta\theta$ (K)'
 ]
 final_contours = [
-    np.linspace(-1.0, 1.0, 21), np.linspace(-1.0, 1.0, 21),
-    np.linspace(-1.0, 1.0, 21), np.linspace(-1.0, 1.0, 21)
+    np.linspace(-1.1, 1.1, 23), np.linspace(-1.4, 1.4, 15),
+    np.linspace(-1.1, 1.1, 23), np.linspace(-1.4, 1.4, 31)
 ]
 
 # ---------------------------------------------------------------------------- #
@@ -129,25 +128,6 @@
     coords_X, coords_Y = extract_gusto_coords(data_file, field_name)
     time = data_file['time'][time_idx]
 
-    # # Filter data for panels that are zoomed in mountain region
-    # if i in [2, 3]:
-    #     data_dict = {
-    #         'X': coords_X,
-    #         'Y': coords_Y,
-    #         'field': field_data
-    #     }
-    #     data_frame = pd.DataFrame(data_dict)
-
-    #     data_frame = data_frame[
-    #         (data_frame['X'] >= xlims_zoom[0])
-    #         & (data_frame['X'] <= xlims_zoom[1])
-    #         & (data_frame['Y'] >= ylims_zoom[0])
-    #         & (data_frame['Y'] <= ylims_zoom[1])
-    #     ]
-    #     field_data = data_frame['field'].values[:]
-    #     coords_X = data_frame['X'].values[:]
-    #     coords_Y = data_frame['Y'].values[:]
-
     field_data, coords_X, coords_Y = \
         reshape_gusto_data(field_data, coords_X, coords_Y)
 
@@ -160,9 +140,12 @@
     )
 
     add_colorbar_ax(ax, cf, field_label, location='bottom')
-    tomplot_field_title(
-        ax, None, minmax=True, field_data=field_data, minmax_format='.3f'
-    )
+
+    if i in [0, 1]:
+        # Only print min/max for top plots
+        tomplot_field_title(
+            ax, None, minmax=True, field_data=field_data, minmax_format='.3f'
+        )
 
     # Labels -------------------------------------------------------------------
     ax.set_xlabel(r'$x$ (km)', labelpad=-10)