diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md index 99da871be..c1f2fa9fa 100644 --- a/docs/developers_guide/ocean/api.md +++ b/docs/developers_guide/ocean/api.md @@ -313,18 +313,21 @@ .. autosummary:: :toctree: generated/ - ConvergenceForward - ConvergenceForward.compute_cell_count - ConvergenceForward.dynamic_model_config - - ConvergenceAnalysis - ConvergenceAnalysis.compute_error - ConvergenceAnalysis.convergence_parameters - ConvergenceAnalysis.exact_solution - ConvergenceAnalysis.get_output_field - ConvergenceAnalysis.plot_convergence - ConvergenceAnalysis.run - ConvergenceAnalysis.setup + get_resolution_for_tasks + get_timestep_for_tasks + + forward.ConvergenceForward + forward.ConvergenceForward.compute_cell_count + forward.ConvergenceForward.dynamic_model_config + + analysis.ConvergenceAnalysis + analysis.ConvergenceAnalysis.compute_error + analysis.ConvergenceAnalysis.convergence_parameters + analysis.ConvergenceAnalysis.exact_solution + analysis.ConvergenceAnalysis.get_output_field + analysis.ConvergenceAnalysis.plot_convergence + analysis.ConvergenceAnalysis.run + analysis.ConvergenceAnalysis.setup ``` ### Spherical Convergence Tests diff --git a/docs/developers_guide/ocean/framework.md b/docs/developers_guide/ocean/framework.md index cca5f7986..d1c324f0a 100644 --- a/docs/developers_guide/ocean/framework.md +++ b/docs/developers_guide/ocean/framework.md @@ -219,16 +219,35 @@ tests on {ref}`dev-ocean-spherical-meshes` and planar meshes. The ocean framework includes shared config options and base classes for forward and analysis steps that are expected to be useful across these tests. +The key config options that control the convergence test are `base_resolution` +and `refinement_factors`. The `base_resolution` is multipled by the +`refinement_factors` to determine which resolutions to test when the +convergence is being tested in space (or space and time together). The +`base_resolution` is applied to all steps when convergence in time is tested. +`base_resolution` times `dt_per_km` determines the base timestep in that case +and is then multiplied by the `refinement_factors` to determine which time steps +to test. When spherical meshes are being tested, the values in the +`convergence` section are overridden by their values in the +`spherical_convergence` section with a prefix indicating the mesh type. + The shared config options are: ```cfg # config options for spherical convergence tests [spherical_convergence] +# The base resolution for the icosahedral mesh to which the refinement +# factors are applied +icos_base_resolution = 60. + # a list of icosahedral mesh resolutions (km) to test -icos_resolutions = 60, 120, 240, 480 +icos_refinement_factors = 8., 4., 2., 1. + +# The base resolution for the quasi-uniform mesh to which the refinement +# factors are applied +qu_base_resolution = 120. # a list of quasi-uniform mesh resolutions (km) to test -qu_resolutions = 60, 90, 120, 150, 180, 210, 240 +qu_refinement_factors = 0.5, 0.75, 1., 1.25, 1.5, 1.75, 2. [convergence] @@ -241,6 +260,15 @@ convergence_thresh = 1.0 # Type of error to compute error_type = l2 +# the base mesh resolution (km) to which refinement_factors +# are applied if refinement is 'space' or 'both' on a planar mesh +# base resolutions for spherical meshes are given in section spherical_convergence +base_resolution = 120 + +# refinement factors for a planar mesh applied to either space or time +# refinement factors for a spherical mesh given in section spherical_convergence +refinement_factors = 4., 2., 1., 0.5 + # config options for convergence forward steps [convergence_forward] @@ -325,7 +353,8 @@ class Forward(SphericalConvergenceForward): bell test case """ - def __init__(self, component, name, subdir, resolution, mesh, init): + def __init__(self, component, name, subdir, mesh, init, + refinement_factor, refinement='both'): """ Create a new step @@ -348,6 +377,9 @@ class Forward(SphericalConvergenceForward): init : polaris.Step The init step + + refinement_factor : float + The factor by which to scale space, time or both """ package = 'polaris.ocean.tasks.cosine_bell' validate_vars = ['normalVelocity', 'tracer1'] @@ -356,7 +388,10 @@ class Forward(SphericalConvergenceForward): init=init, package=package, yaml_filename='forward.yaml', output_filename='output.nc', - validate_vars=validate_vars) + validate_vars=validate_vars, + graph_target=f'{init.path}/graph.info', + refinement_factor=refinement_factor, + refinement=refinement) ``` Each convergence test must define a YAML file with model config options, called `forward.yaml` by default. The `package` parameter is the location of this @@ -419,7 +454,7 @@ class Analysis(ConvergenceAnalysis): """ A step for analyzing the output from the cosine bell test case """ - def __init__(self, component, resolutions, subdir, dependencies): + def __init__(self, component, subdir, dependencies, refinement='both'): """ Create the step @@ -436,6 +471,9 @@ class Analysis(ConvergenceAnalysis): dependencies : dict of dict of polaris.Steps The dependencies of this step + + refinement : str, optional + Whether to refine in space, time or both space and time """ convergence_vars = [{'name': 'tracer1', 'title': 'tracer1', @@ -443,7 +481,8 @@ class Analysis(ConvergenceAnalysis): super().__init__(component=component, subdir=subdir, resolutions=resolutions, dependencies=dependencies, - convergence_vars=convergence_vars) + convergence_vars=convergence_vars, + refinement=refinement) ``` Many tasks will also need to override the diff --git a/docs/users_guide/ocean/tasks/cosine_bell.md b/docs/users_guide/ocean/tasks/cosine_bell.md index f49afa97a..5230fa36b 100644 --- a/docs/users_guide/ocean/tasks/cosine_bell.md +++ b/docs/users_guide/ocean/tasks/cosine_bell.md @@ -33,34 +33,56 @@ These tasks support only MPAS-Ocean. (ocean-cosine-bell-mesh)= ## mesh -Two global mesh variants are tested, quasi-uniform (QU) and icosohydral. Thus, -there are 4 variants of the task: +Two global mesh variants are tested, quasi-uniform (QU) and icosohydral. There +are also variants to test convergence in space, time, or both space and time. +In addition, the tests can be set up with or without the viz step. Thus, there +are 12 variants of the task: ``` -ocean/spherical/icos/cosine_bell -ocean/spherical/icos/cosine_bell/with_viz -ocean/spherical/qu/cosine_bell -ocean/spherical/qu/cosine_bell/with_viz +ocean/spherical/icos/cosine_bell/convergence_both/ +ocean/spherical/icos/cosine_bell/convergence_both/with_viz +ocean/spherical/qu/cosine_bell/convergence_both/ +ocean/spherical/qu/cosine_bell/convergence_both/with_viz +ocean/spherical/icos/cosine_bell/convergence_space/ +ocean/spherical/icos/cosine_bell/convergence_space/with_viz +ocean/spherical/qu/cosine_bell/convergence_space/ +ocean/spherical/qu/cosine_bell/convergence_space/with_viz +ocean/spherical/icos/cosine_bell/convergence_time/ +ocean/spherical/icos/cosine_bell/convergence_time/with_viz +ocean/spherical/qu/cosine_bell/convergence_time/ +ocean/spherical/qu/cosine_bell/convergence_time/with_viz ``` The default resolutions used in the task depends on the mesh type. -For the `icos` mesh type, the defaults are: +For the `icos` mesh type, the defaults are 60, 120, 240, 480 km, as determined +by the following config options. See {ref}`dev-ocean-convergence` for more +details. ```cfg # config options for spherical convergence tests [spherical_convergence] +# The base resolution for the icosahedral mesh to which the refinement +# factors are applied +icos_base_resolution = 60. + +# a list of icosahedral mesh resolutions (km) to test +icos_refinement_factors = 8., 4., 2., 1. # a list of icosahedral mesh resolutions (km) to test -icos_resolutions = 60, 120, 240, 480 ``` -for the `qu` mesh type, they are: +For the `qu` mesh type, they are 60, 90, 120, 150, 180, 210, 240 km as +determined by the following config options: ```cfg # config options for spherical convergence tests [spherical_convergence] +# The base resolution for the quasi-uniform mesh to which the refinement +# factors are applied +qu_base_resolution = 120. + # a list of quasi-uniform mesh resolutions (km) to test -qu_resolutions = 60, 90, 120, 150, 180, 210, 240 +qu_refinement_factors = 0.5, 0.75, 1., 1.25, 1.5, 1.75, 2. ``` To alter the resolutions used in this task, you will need to create your own diff --git a/docs/users_guide/ocean/tasks/geostrophic.md b/docs/users_guide/ocean/tasks/geostrophic.md index 82301a154..614eeacce 100644 --- a/docs/users_guide/ocean/tasks/geostrophic.md +++ b/docs/users_guide/ocean/tasks/geostrophic.md @@ -4,20 +4,25 @@ ## description -The `geostrophic` and `geostrophic/with_viz` tasks implement the "Global Steady +The `geostrophic` tasks implement the "Global Steady State Nonlinear Zonal Geostrophic Flow" test case described in [Williamson et al. 1992]() -The task is a convergence test with time step varying proportionately to grid -size. The result of the `analysis` step of the task are plots like the -following showing convergence of water-column thickness and normal velocity as -functions of the mesh resolution: +The task `geostrophic/convergence_both` is a space-and-time convergence test +with time step varying proportionately to grid size. Convergence tests in +either space or time are also available as `convergence_space` or +`convergence_time`. The result of the `analysis` step of the task are plots +like the following showing convergence of water-column thickness and normal +velocity as functions of the mesh resolution: ```{image} images/geostrophic_convergence.png :align: center :width: 500 px ``` +Visualizations of the fields themselves can be added in viz steps by appending +`with_viz` to the task name at setup. + ## suppported models These tasks support only MPAS-Ocean. diff --git a/docs/users_guide/ocean/tasks/inertial_gravity_wave.md b/docs/users_guide/ocean/tasks/inertial_gravity_wave.md index f858f232e..709166e4c 100644 --- a/docs/users_guide/ocean/tasks/inertial_gravity_wave.md +++ b/docs/users_guide/ocean/tasks/inertial_gravity_wave.md @@ -11,6 +11,14 @@ The implementation is from ## description +There are 3 versions of the convergence test case, `convergence_space`, +`convergence_time`, and `convergence_both` corresponding to space, time, and +space and time convergence tests. Tests involving spatial convergence run the +inertial gravity wave simulation for 4 different resolutions: 200, 100, 50, +and 25 km. Tests involving temporal convergence use the parameter `dt_per_km` +at the `base_resolution` multiplied by `refinement_factors` (see +{ref}`dev-ocean-convergence` for more details on how to change resolutions or +time steps tested). The `inertial_gravity_wave` task runs the inertial gravity wave simulation for 4 different resolutions: 200, 100, 50, and 25 km. diff --git a/docs/users_guide/ocean/tasks/manufactured_solution.md b/docs/users_guide/ocean/tasks/manufactured_solution.md index f9e840e0f..866f2d85d 100644 --- a/docs/users_guide/ocean/tasks/manufactured_solution.md +++ b/docs/users_guide/ocean/tasks/manufactured_solution.md @@ -22,8 +22,14 @@ These tasks support both MPAS-Ocean and Omega. ### description -The `convergence` test case runs the manufactured solution simulation for 4 -different resolutions: 200, 100, 50, and 25 km. +There are 3 versions of the convergence test case, `convergence_space`, +`convergence_time`, and `convergence_both` corresponding to space, time, and +space and time convergence tests. Tests involving spatial convergence run the +manufactured solution simulation for 4 different resolutions: 200, 100, 50, +and 25 km. Tests involving temporal convergence use the parameter `dt_per_km` +at the `base_resolution` multiplied by `refinement_factors` (see +{ref}`dev-ocean-convergence` for more details on how to change resolutions or +time steps tested). The forward step for each resolution runs the simulation for 10 hours. The model is configured without vertical advection and mixing. No tracers are enabled diff --git a/polaris/ocean/convergence/__init__.py b/polaris/ocean/convergence/__init__.py index 7fb71d53f..0a04e50e7 100644 --- a/polaris/ocean/convergence/__init__.py +++ b/polaris/ocean/convergence/__init__.py @@ -1,2 +1,106 @@ -from polaris.ocean.convergence.analysis import ConvergenceAnalysis -from polaris.ocean.convergence.forward import ConvergenceForward +import numpy as np + + +def get_resolution_for_task(config, refinement_factor, + refinement='both'): + """ + Get the resolution for a step in a convergence task + + Parameters + ---------- + config : polaris.Config + The config options for this task + + refinement_factor : float + The factor by which either resolution or time is refined for this step + + refinement : str, optional + Whether to refine in space, time or both + + Returns + ------- + resolution : float + The resolution corresponding to the refinement_factor and convergence + test type + """ + if refinement == 'both': + option = 'refinement_factors_space' + else: + option = f'refinement_factors_{refinement}' + base_resolution = config.getfloat('convergence', 'base_resolution') + refinement_factors = config.getlist('convergence', option, dtype=float) + + if refinement_factor not in refinement_factors: + raise ValueError( + 'refinement_factor not found in config option refinement_factors') + + if refinement == 'time': + resolution = base_resolution + else: + resolution = refinement_factor * base_resolution + + return resolution + + +def get_timestep_for_task(config, refinement_factor, + refinement='both'): + """ + Get the time step for a forward step in a convergence task + + Parameters + ---------- + config : polaris.Config + The config options for this task + + refinement_factor : float + The factor by which either resolution or time is refined for this step + + refinement : str, optional + Whether to refine in space, time or both + + Returns + ------- + resolution : float + The resolution corresponding to the refinement_factor and convergence + test type + """ + + if refinement == 'both': + option = 'refinement_factors_space' + else: + option = f'refinement_factors_{refinement}' + base_resolution = config.getfloat('convergence', 'base_resolution') + refinement_factors = config.getlist('convergence', option, dtype=float) + + if refinement_factor not in refinement_factors: + raise ValueError( + 'refinement_factor not found in config option refinement_factors') + + resolution = get_resolution_for_task( + config, refinement_factor, refinement=refinement) + + section = config['convergence_forward'] + time_integrator = section.get('time_integrator') + # dt is proportional to resolution: default 30 seconds per km + if time_integrator == 'RK4': + dt_per_km = section.getfloat('rk4_dt_per_km') + btr_timestep = 0. + else: + dt_per_km = section.getfloat('split_dt_per_km') + btr_dt_per_km = section.getfloat('btr_dt_per_km') + for idx, refinement_factor in enumerate(refinement_factors): + if refinement == 'time': + btr_timestep = btr_dt_per_km * base_resolution * \ + refinement_factor + elif refinement == 'space': + btr_timestep = btr_dt_per_km * base_resolution + else: + btr_timestep = btr_dt_per_km * resolution + if refinement == 'time': + timestep = dt_per_km * refinement_factor * resolution + elif refinement == 'space': + timestep = dt_per_km * base_resolution + else: + timestep = dt_per_km * resolution + + return timestep, btr_timestep diff --git a/polaris/ocean/convergence/analysis.py b/polaris/ocean/convergence/analysis.py index f2bbf85a1..54d3b5ec4 100644 --- a/polaris/ocean/convergence/analysis.py +++ b/polaris/ocean/convergence/analysis.py @@ -5,8 +5,11 @@ import pandas as pd from polaris.mpas import area_for_field, time_index_from_xtime +from polaris.ocean.convergence import ( + get_resolution_for_task, + get_timestep_for_task, +) from polaris.ocean.model import OceanIOStep -from polaris.ocean.resolution import resolution_to_subdir from polaris.viz import use_mplstyle @@ -16,9 +19,6 @@ class ConvergenceAnalysis(OceanIOStep): Attributes ---------- - resolutions : list of float - The resolutions of the meshes that have been run - dependencies_dict : dict of dict of polaris.Steps The dependencies of this step must be given as separate keys in the dict: @@ -52,9 +52,13 @@ class ConvergenceAnalysis(OceanIOStep): zidx : int The z-index to use for variables that have an nVertLevels dimension, which should be None for variables that don't + + refinement : str + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ - def __init__(self, component, resolutions, subdir, dependencies, - convergence_vars): + def __init__(self, component, subdir, dependencies, convergence_vars, + refinement='both'): """ Create the step @@ -63,9 +67,6 @@ def __init__(self, component, resolutions, subdir, dependencies, component : polaris.Component The component the step belongs to - resolutions : list of float - The resolutions of the meshes that have been run - subdir : str The subdirectory that the step resides in @@ -74,17 +75,17 @@ def __init__(self, component, resolutions, subdir, dependencies, dict: mesh : dict of polaris.Steps - Keys of the dict correspond to `resolutions` + Keys of the dict correspond to `refinement_factors` Values of the dict are polaris.Steps, which must have the attribute `path`, the path to `base_mesh.nc` of that resolution init : dict of polaris.Steps - Keys of the dict correspond to `resolutions` + Keys of the dict correspond to `refinement_factors` Values of the dict are polaris.Steps, which must have the attribute `path`, the path to `initial_state.nc` of that resolution forward : dict of polaris.Steps - Keys of the dict correspond to `resolutions` + Keys of the dict correspond to `refinement_factors` Values of the dict are polaris.Steps, which must have the attribute `path`, the path to `forward.nc` of that resolution @@ -102,11 +103,15 @@ def __init__(self, component, resolutions, subdir, dependencies, zidx : int The z-index to use for variables that have an nVertLevels dimension, which should be None for variables that don't + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ super().__init__(component=component, name='analysis', subdir=subdir) - self.resolutions = resolutions self.dependencies_dict = dependencies self.convergence_vars = convergence_vars + self.refinement = refinement for var in convergence_vars: self.add_output_file(f'convergence_{var["name"]}.png') @@ -116,21 +121,26 @@ def setup(self): Add input files based on resolutions, which may have been changed by user config options """ + config = self.config dependencies = self.dependencies_dict - - for resolution in self.resolutions: - mesh_name = resolution_to_subdir(resolution) - base_mesh = dependencies['mesh'][resolution] - init = dependencies['init'][resolution] - forward = dependencies['forward'][resolution] + if self.refinement == 'time': + option = 'refinement_factors_time' + else: + option = 'refinement_factors_space' + refinement_factors = config.getlist('convergence', option, + dtype=float) + for refinement_factor in refinement_factors: + base_mesh = dependencies['mesh'][refinement_factor] + init = dependencies['init'][refinement_factor] + forward = dependencies['forward'][refinement_factor] self.add_input_file( - filename=f'{mesh_name}_mesh.nc', + filename=f'mesh_r{refinement_factor:02g}.nc', work_dir_target=f'{base_mesh.path}/base_mesh.nc') self.add_input_file( - filename=f'{mesh_name}_init.nc', + filename=f'init_r{refinement_factor:02g}.nc', work_dir_target=f'{init.path}/initial_state.nc') self.add_input_file( - filename=f'{mesh_name}_output.nc', + filename=f'output_r{refinement_factor:02g}.nc', work_dir_target=f'{forward.path}/output.nc') def run(self): @@ -162,37 +172,59 @@ def plot_convergence(self, variable_name, title, zidx): The z-index to use for variables that have an nVertLevels dimension, which should be None for variables that don't """ - resolutions = self.resolutions + config = self.config logger = self.logger conv_thresh, error_type = self.convergence_parameters( field_name=variable_name) error = [] - for resolution in resolutions: - mesh_name = resolution_to_subdir(resolution) + if self.refinement == 'time': + option = 'refinement_factors_time' + header = 'time step' + else: + option = 'refinement_factors_space' + header = 'resolution' + refinement_factors = config.getlist('convergence', option, + dtype=float) + resolutions = list() + timesteps = list() + for refinement_factor in refinement_factors: + timestep, _ = get_timestep_for_task( + config, refinement_factor, refinement=self.refinement) + timesteps.append(timestep) + resolution = get_resolution_for_task( + config, refinement_factor, refinement=self.refinement) + resolutions.append(resolution) error_res = self.compute_error( - mesh_name=mesh_name, + refinement_factor=refinement_factor, variable_name=variable_name, zidx=zidx, error_type=error_type) error.append(error_res) - res_array = np.array(resolutions) + if self.refinement == 'time': + refinement_array = np.array(timesteps) + x_label = 'Time (s)' + else: + refinement_array = np.array(resolutions) + x_label = 'Horizontal resolution (km)' error_array = np.array(error) filename = f'convergence_{variable_name}.csv' - data = np.stack((res_array, error_array), axis=1) - df = pd.DataFrame(data, columns=['resolution', error_type]) + data = np.stack((refinement_array, error_array), axis=1) + df = pd.DataFrame(data, columns=[header, error_type]) df.to_csv(f'convergence_{variable_name}.csv', index=False) convergence_failed = False - poly = np.polyfit(np.log10(res_array), np.log10(error_array), 1) + poly = np.polyfit(np.log10(refinement_array), np.log10(error_array), 1) convergence = poly[0] conv_round = np.round(convergence, 3) - fit = res_array ** poly[0] * 10 ** poly[1] + fit = refinement_array ** poly[0] * 10 ** poly[1] - order1 = 0.5 * error_array[-1] * (res_array / res_array[-1]) - order2 = 0.5 * error_array[-1] * (res_array / res_array[-1])**2 + order1 = 0.5 * error_array[-1] * \ + (refinement_array / refinement_array[-1]) + order2 = 0.5 * error_array[-1] * \ + (refinement_array / refinement_array[-1])**2 use_mplstyle() fig = plt.figure() @@ -205,9 +237,9 @@ def plot_convergence(self, variable_name, title, zidx): alpha=0.3) ax.loglog(resolutions, order2, 'k', label='second order', alpha=0.3) - ax.loglog(res_array, fit, 'k', + ax.loglog(refinement_array, fit, 'k', label=f'linear fit (order={conv_round})') - ax.loglog(res_array, error_array, 'o', label='numerical') + ax.loglog(refinement_array, error_array, 'o', label='numerical') if self.baseline_dir is not None: baseline_filename = os.path.join(self.baseline_dir, filename) @@ -217,20 +249,20 @@ def plot_convergence(self, variable_name, title, zidx): raise ValueError( f'{error_type} not available for baseline') else: - res_array = data.resolution.to_numpy() + refinement_array = data.resolution.to_numpy() error_array = data[error_type].to_numpy() poly = np.polyfit( - np.log10(res_array), np.log10(error_array), 1) + np.log10(refinement_array), np.log10(error_array), 1) base_convergence = poly[0] conv_round = np.round(base_convergence, 3) - fit = res_array ** poly[0] * 10 ** poly[1] - ax.loglog(res_array, error_array, 'o', color='#ff7f0e', - label='baseline') - ax.loglog(res_array, fit, color='#ff7f0e', + fit = refinement_array ** poly[0] * 10 ** poly[1] + ax.loglog(refinement_array, error_array, 'o', + color='#ff7f0e', label='baseline') + ax.loglog(refinement_array, fit, color='#ff7f0e', label=f'linear fit, baseline ' f'(order={conv_round})') - ax.set_xlabel('resolution (km)') + ax.set_xlabel(x_label) ax.set_ylabel(f'{error_title}') ax.set_title(f'Error Convergence of {title}') ax.legend(loc='lower left') @@ -251,15 +283,15 @@ def plot_convergence(self, variable_name, title, zidx): if convergence_failed: raise ValueError('Convergence rate below minimum tolerance.') - def compute_error(self, mesh_name, variable_name, zidx=None, + def compute_error(self, refinement_factor, variable_name, zidx=None, error_type='l2'): """ Compute the error for a given resolution Parameters ---------- - mesh_name : str - The name of the mesh + refinement_factor : float + The factor by which step is refined in space, time or both variable_name : str The variable name of which to evaluate error with respect to the @@ -277,16 +309,16 @@ def compute_error(self, mesh_name, variable_name, zidx=None, The error of the variable given by variable_name """ norm_type = {'l2': None, 'inf': np.inf} - ds_mesh = self.open_model_dataset(f'{mesh_name}_mesh.nc') + ds_mesh = self.open_model_dataset(f'mesh_r{refinement_factor:02g}.nc') config = self.config section = config['convergence'] eval_time = section.getfloat('convergence_eval_time') s_per_hour = 3600.0 - field_exact = self.exact_solution(mesh_name, variable_name, + field_exact = self.exact_solution(refinement_factor, variable_name, time=eval_time * s_per_hour, zidx=zidx) - field_mpas = self.get_output_field(mesh_name, variable_name, + field_mpas = self.get_output_field(refinement_factor, variable_name, time=eval_time * s_per_hour, zidx=zidx) diff = field_exact - field_mpas @@ -304,14 +336,14 @@ def compute_error(self, mesh_name, variable_name, zidx=None, return error - def exact_solution(self, mesh_name, field_name, time, zidx=None): + def exact_solution(self, refinement_factor, field_name, time, zidx=None): """ Get the exact solution Parameters ---------- - mesh_name : str - The mesh name which is the prefix for the initial condition file + refinement_factor : float + The factor by which step is refined in space, time or both field_name : str The name of the variable of which we evaluate convergence @@ -332,20 +364,20 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): The exact solution as derived from the initial condition """ - ds_init = self.open_model_dataset(f'{mesh_name}_init.nc') + ds_init = self.open_model_dataset(f'init_r{refinement_factor:02g}.nc') ds_init = ds_init.isel(Time=0) if zidx is not None: ds_init = ds_init.isel(nVertLevels=zidx) return ds_init[field_name] - def get_output_field(self, mesh_name, field_name, time, zidx=None): + def get_output_field(self, refinement_factor, field_name, time, zidx=None): """ Get the model output field at the given time and z index Parameters ---------- - mesh_name : str + e str The mesh name which is the prefix for the output file field_name : str @@ -362,7 +394,7 @@ def get_output_field(self, mesh_name, field_name, time, zidx=None): field_mpas : xarray.DataArray model output field """ - ds_out = self.open_model_dataset(f'{mesh_name}_output.nc') + ds_out = self.open_model_dataset(f'output_r{refinement_factor:02g}.nc') tidx = time_index_from_xtime(ds_out.xtime.values, time) ds_out = ds_out.isel(Time=tidx) diff --git a/polaris/ocean/convergence/convergence.cfg b/polaris/ocean/convergence/convergence.cfg index 5f8b0752d..9957bc0d2 100644 --- a/polaris/ocean/convergence/convergence.cfg +++ b/polaris/ocean/convergence/convergence.cfg @@ -4,11 +4,24 @@ convergence_eval_time = 24.0 # Convergence threshold below which a test fails -convergence_thresh = 1.0 +convergence_thresh_space = 1.0 +convergence_thresh_time = 1.0 # Type of error to compute error_type = l2 +# the base mesh resolution (km) to which refinement_factors are applied +# base resolutions for spherical meshes are given in section spherical_convergence +base_resolution = 120 + +# refinement factors for a planar mesh applied to either space or both space and time +# refinement factors for a spherical mesh given in section spherical_convergence +refinement_factors_space = 4., 2., 1., 0.5 + +# refinement factors for a planar mesh applied to time with the base timestep +# determined by base_resolution * dt_per_km +refinement_factors_time = 1., 0.5, 0.25 + # config options for convergence forward steps [convergence_forward] @@ -18,9 +31,11 @@ error_type = l2 time_integrator = RK4 # RK4 time step per resolution (s/km), since dt is proportional to resolution +# if using convergence in time only, this is used for the largest resolution rk4_dt_per_km = 3.0 # split time step per resolution (s/km), since dt is proportional to resolution +# if using convergence in time only, this is used for the largest resolution split_dt_per_km = 30.0 # the barotropic time step (s/km) for simulations using split time stepping, diff --git a/polaris/ocean/convergence/forward.py b/polaris/ocean/convergence/forward.py index 2866853e2..7a7ecbc2e 100644 --- a/polaris/ocean/convergence/forward.py +++ b/polaris/ocean/convergence/forward.py @@ -1,3 +1,4 @@ +from polaris.ocean.convergence import get_timestep_for_task from polaris.ocean.model import OceanModelStep, get_time_interval_string @@ -8,21 +9,25 @@ class ConvergenceForward(OceanModelStep): Attributes ---------- - resolution : float - The resolution of the (uniform) mesh in km - package : Package The package name or module object that contains the YAML file yaml_filename : str A YAML file that is a Jinja2 template with model config options + refinement : str + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time + + refinement_factor : float + Refinement factor use to scale space, time or both depending on + refinement option """ - def __init__(self, component, name, subdir, resolution, mesh, init, - package, yaml_filename='forward.yaml', options=None, - graph_target=None, output_filename='output.nc', - validate_vars=None): + def __init__(self, component, name, subdir, refinement_factor, mesh, init, + package, yaml_filename='forward.yaml', + options=None, graph_target=None, output_filename='output.nc', + validate_vars=None, refinement='both'): """ Create a new step @@ -37,8 +42,9 @@ def __init__(self, component, name, subdir, resolution, mesh, init, subdir : str The subdirectory for the step - resolution : float - The resolution of the (uniform) mesh in km + refinement_factor : float + Refinement factor use to scale space, time or both depending on + refinement option package : Package The package name or module object that contains the YAML file @@ -60,11 +66,16 @@ def __init__(self, component, name, subdir, resolution, mesh, init, validate_vars : list of str, optional A list of variables to validate against a baseline if requested + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ super().__init__(component=component, name=name, subdir=subdir, openmp_threads=1, graph_target=graph_target) - self.resolution = resolution + self.refinement = refinement + self.refinement_factor = refinement_factor self.package = package self.yaml_filename = yaml_filename @@ -109,27 +120,19 @@ def dynamic_model_config(self, at_setup): super().dynamic_model_config(at_setup=at_setup) config = self.config + refinement_factor = self.refinement_factor vert_levels = config.getfloat('vertical_grid', 'vert_levels') if not at_setup and vert_levels == 1: self.add_yaml_file('polaris.ocean.config', 'single_layer.yaml') - section = config['convergence_forward'] - - time_integrator = section.get('time_integrator') - # dt is proportional to resolution: default 30 seconds per km - if time_integrator == 'RK4': - dt_per_km = section.getfloat('rk4_dt_per_km') - else: - dt_per_km = section.getfloat('split_dt_per_km') - dt_str = get_time_interval_string(seconds=dt_per_km * self.resolution) - - # btr_dt is also proportional to resolution: default 1.5 seconds per km - btr_dt_per_km = section.getfloat('btr_dt_per_km') - btr_dt_str = get_time_interval_string( - seconds=btr_dt_per_km * self.resolution) + timestep, btr_timestep = get_timestep_for_task( + config, refinement_factor, refinement=self.refinement) + dt_str = get_time_interval_string(seconds=timestep) + btr_dt_str = get_time_interval_string(seconds=btr_timestep) s_per_hour = 3600. + section = config['convergence_forward'] run_duration = section.getfloat('run_duration') run_duration_str = get_time_interval_string( seconds=run_duration * s_per_hour) @@ -138,6 +141,7 @@ def dynamic_model_config(self, at_setup): output_interval_str = get_time_interval_string( seconds=output_interval * s_per_hour) + time_integrator = section.get('time_integrator') time_integrator_map = dict([('RK4', 'RungeKutta4')]) model = config.get('ocean', 'model') if model == 'omega': diff --git a/polaris/ocean/convergence/spherical/forward.py b/polaris/ocean/convergence/spherical/forward.py index 88be45217..26f61ac47 100644 --- a/polaris/ocean/convergence/spherical/forward.py +++ b/polaris/ocean/convergence/spherical/forward.py @@ -1,4 +1,5 @@ -from polaris.ocean.convergence import ConvergenceForward +from polaris.ocean.convergence import get_resolution_for_task +from polaris.ocean.convergence.forward import ConvergenceForward class SphericalConvergenceForward(ConvergenceForward): @@ -30,5 +31,8 @@ def compute_cell_count(self): The approximate number of cells in the mesh """ # use a heuristic based on QU30 (65275 cells) and QU240 (10383 cells) - cell_count = 6e8 / self.resolution**2 + resolution = get_resolution_for_task( + self.config, self.refinement_factor, + refinement=self.refinement) + cell_count = 6e8 / resolution**2 return cell_count diff --git a/polaris/ocean/convergence/spherical/spherical.cfg b/polaris/ocean/convergence/spherical/spherical.cfg index 69ea2dd22..cb4b8d4fd 100644 --- a/polaris/ocean/convergence/spherical/spherical.cfg +++ b/polaris/ocean/convergence/spherical/spherical.cfg @@ -1,8 +1,18 @@ # config options for spherical convergence tests [spherical_convergence] -# a list of icosahedral mesh resolutions (km) to test -icos_resolutions = 60, 120, 240, 480 +# The base resolution for the icosahedral mesh to which the refinement +# factors are applied +icos_base_resolution = 60. -# a list of quasi-uniform mesh resolutions (km) to test -qu_resolutions = 60, 90, 120, 150, 180, 210, 240 +# The factors by which to scale space or time based on icos_base_resolution +icos_refinement_factors_space = 8., 4., 2., 1. +icos_refinement_factors_time = 1., 0.5, 0.25 + +# The base resolution for the quasi-uniform mesh to which the refinement +# factors are applied +qu_base_resolution = 120. + +# The factors by which to scale space or time based on qu_base_resolution +qu_refinement_factors_space = 2., 1.75, 1.5, 1.25, 1., 0.75, 0.5 +qu_refinement_factors_time = 1., 0.5, 0.25 diff --git a/polaris/ocean/suites/convergence.txt b/polaris/ocean/suites/convergence.txt index 56a8e63c7..9a92ef47e 100644 --- a/polaris/ocean/suites/convergence.txt +++ b/polaris/ocean/suites/convergence.txt @@ -1,22 +1,22 @@ -ocean/planar/inertial_gravity_wave -ocean/planar/manufactured_solution +ocean/planar/inertial_gravity_wave/convergence_both +ocean/planar/manufactured_solution/convergence_both ocean/spherical/icos/correlated_tracers_2d cached: icos_base_mesh_60km icos_base_mesh_120km icos_base_mesh_240km icos_base_mesh_480km ocean/spherical/qu/correlated_tracers_2d cached: qu_base_mesh_60km qu_base_mesh_90km qu_base_mesh_120km qu_base_mesh_150km cached: qu_base_mesh_180km qu_base_mesh_210km qu_base_mesh_240km -ocean/spherical/icos/cosine_bell +ocean/spherical/icos/cosine_bell/convergence_both cached: icos_base_mesh_60km icos_init_60km icos_base_mesh_120km icos_init_120km cached: icos_base_mesh_240km icos_init_240km icos_base_mesh_480km icos_init_480km -ocean/spherical/qu/cosine_bell +ocean/spherical/qu/cosine_bell/convergence_both cached: qu_base_mesh_60km qu_init_60km qu_base_mesh_90km qu_init_90km cached: qu_base_mesh_120km qu_init_120km qu_base_mesh_150km qu_init_150km cached: qu_base_mesh_180km qu_init_180km qu_base_mesh_210km qu_init_210km cached: qu_base_mesh_240km qu_init_240km -ocean/spherical/icos/geostrophic +ocean/spherical/icos/geostrophic/convergence_both cached: icos_base_mesh_60km icos_init_60km icos_base_mesh_120km icos_init_120km cached: icos_base_mesh_240km icos_init_240km icos_base_mesh_480km icos_init_480km -ocean/spherical/qu/geostrophic +ocean/spherical/qu/geostrophic/convergence_both cached: qu_base_mesh_60km qu_init_60km qu_base_mesh_90km qu_init_90km cached: qu_base_mesh_120km qu_init_120km qu_base_mesh_150km qu_init_150km cached: qu_base_mesh_180km qu_init_180km qu_base_mesh_210km qu_init_210km diff --git a/polaris/ocean/suites/geostrophic.txt b/polaris/ocean/suites/geostrophic.txt index a8a573e24..7a88f6da9 100644 --- a/polaris/ocean/suites/geostrophic.txt +++ b/polaris/ocean/suites/geostrophic.txt @@ -1,2 +1,2 @@ -ocean/spherical/icos/geostrophic -ocean/spherical/qu/geostrophic +ocean/spherical/icos/geostrophic/convergence_both +ocean/spherical/qu/geostrophic/convergence_both diff --git a/polaris/ocean/suites/nightly.txt b/polaris/ocean/suites/nightly.txt index 725b504e1..e2b2ef4aa 100644 --- a/polaris/ocean/suites/nightly.txt +++ b/polaris/ocean/suites/nightly.txt @@ -3,5 +3,5 @@ ocean/planar/baroclinic_channel/10km/decomp ocean/planar/baroclinic_channel/10km/restart ocean/planar/ice_shelf_2d/5km/z-star/default/with_restart ocean/planar/ice_shelf_2d/5km/z-level/default/with_restart -ocean/planar/inertial_gravity_wave +ocean/planar/inertial_gravity_wave/convergence_both # ocean/planar/manufactured_solution diff --git a/polaris/ocean/suites/pr.txt b/polaris/ocean/suites/pr.txt index 3d63ada67..4bb3cf52a 100644 --- a/polaris/ocean/suites/pr.txt +++ b/polaris/ocean/suites/pr.txt @@ -3,7 +3,7 @@ ocean/planar/baroclinic_channel/10km/decomp ocean/planar/baroclinic_channel/10km/restart ocean/planar/ice_shelf_2d/5km/z-star/default/with_restart ocean/planar/ice_shelf_2d/5km/z-level/default/with_restart -ocean/planar/inertial_gravity_wave +ocean/planar/inertial_gravity_wave/convergence_both ocean/planar/internal_wave/standard/default ocean/planar/internal_wave/vlr/default # ocean/planar/manufactured_solution diff --git a/polaris/ocean/tasks/cosine_bell/__init__.py b/polaris/ocean/tasks/cosine_bell/__init__.py index 89124c644..07b4e20bc 100644 --- a/polaris/ocean/tasks/cosine_bell/__init__.py +++ b/polaris/ocean/tasks/cosine_bell/__init__.py @@ -1,7 +1,12 @@ +from math import ceil from typing import Dict from polaris import Step, Task from polaris.config import PolarisConfigParser +from polaris.ocean.convergence import ( + get_resolution_for_task, + get_timestep_for_task, +) from polaris.ocean.mesh.spherical import add_spherical_base_mesh_step from polaris.ocean.tasks.cosine_bell.analysis import Analysis from polaris.ocean.tasks.cosine_bell.forward import Forward @@ -28,11 +33,13 @@ def add_cosine_bell_tasks(component): config.add_from_package('polaris.ocean.tasks.cosine_bell', 'cosine_bell.cfg') - for include_viz in [False, True]: - component.add_task(CosineBell(component=component, - config=config, - icosahedral=icosahedral, - include_viz=include_viz)) + for refinement in ['space', 'time', 'both']: + for include_viz in [False, True]: + component.add_task(CosineBell(component=component, + config=config, + icosahedral=icosahedral, + include_viz=include_viz, + refinement=refinement)) class CosineBell(Task): @@ -41,8 +48,9 @@ class CosineBell(Task): Attributes ---------- - resolutions : list of float - A list of mesh resolutions + refinement : str + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time icosahedral : bool Whether to use icosahedral, as opposed to less regular, JIGSAW meshes @@ -50,7 +58,8 @@ class CosineBell(Task): include_viz : bool Include VizMap and Viz steps for each resolution """ - def __init__(self, component, config, icosahedral, include_viz): + def __init__(self, component, config, icosahedral, include_viz, + refinement='both'): """ Create the convergence test @@ -68,29 +77,30 @@ def __init__(self, component, config, icosahedral, include_viz): include_viz : bool Include VizMap and Viz steps for each resolution + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ if icosahedral: prefix = 'icos' else: prefix = 'qu' - subdir = f'spherical/{prefix}/cosine_bell' - name = f'{prefix}_cosine_bell' + subdir = f'spherical/{prefix}/cosine_bell/convergence_{refinement}' + name = f'{prefix}_cosine_bell_convergence_{refinement}' if include_viz: subdir = f'{subdir}/with_viz' name = f'{name}_with_viz' - link = 'cosine_bell.cfg' - else: - # config options live in the task already so no need for a symlink - link = None + link = 'cosine_bell.cfg' super().__init__(component=component, name=name, subdir=subdir) - self.resolutions = list() + self.refinement = refinement self.icosahedral = icosahedral self.include_viz = include_viz self.set_shared_config(config, link=link) - self._setup_steps() + self._setup_steps(refinement=refinement) def configure(self): """ @@ -99,10 +109,17 @@ def configure(self): super().configure() # set up the steps again in case a user has provided new resolutions - self._setup_steps() + self._setup_steps(self.refinement) + + def _setup_steps(self, refinement): + """ + setup steps given resolutions - def _setup_steps(self): - """ setup steps given resolutions """ + Parameters + ---------- + refinement : str + refinement type. One of 'space', 'time' or 'both' + """ icosahedral = self.icosahedral config = self.config config_filename = self.config_filename @@ -112,85 +129,107 @@ def _setup_steps(self): else: prefix = 'qu' - resolutions = config.getlist('spherical_convergence', - f'{prefix}_resolutions', dtype=float) - - if self.resolutions == resolutions: - return + if refinement == 'time': + option = 'refinement_factors_time' + else: + option = 'refinement_factors_space' + refinement_factors = config.getlist('spherical_convergence', + f'{prefix}_{option}', dtype=str) + refinement_factors = ', '.join(refinement_factors) + config.set('convergence', option, value=refinement_factors) + refinement_factors = config.getlist('convergence', + option, dtype=float) + + base_resolution = config.getfloat('spherical_convergence', + f'{prefix}_base_resolution') + config.set('convergence', 'base_resolution', + value=f'{base_resolution:03g}') # start fresh with no steps for step in list(self.steps.values()): self.remove_step(step) - self.resolutions = resolutions - component = self.component analysis_dependencies: Dict[str, Dict[str, Step]] = ( dict(mesh=dict(), init=dict(), forward=dict())) - for resolution in resolutions: + + resolutions = list() + timesteps = list() + + case_dir = f'spherical/{prefix}/cosine_bell' + + for refinement_factor in refinement_factors: + resolution = get_resolution_for_task( + config, refinement_factor, refinement=refinement) + base_mesh_step, mesh_name = add_spherical_base_mesh_step( component, resolution, icosahedral) - self.add_step(base_mesh_step, symlink=f'base_mesh/{mesh_name}') - analysis_dependencies['mesh'][resolution] = base_mesh_step - - cos_bell_dir = f'spherical/{prefix}/cosine_bell' + analysis_dependencies['mesh'][refinement_factor] = base_mesh_step name = f'{prefix}_init_{mesh_name}' - subdir = f'{cos_bell_dir}/init/{mesh_name}' - if self.include_viz: - symlink = f'init/{mesh_name}' - else: - symlink = None + subdir = f'{case_dir}/init/{mesh_name}' if subdir in component.steps: init_step = component.steps[subdir] else: - init_step = Init(component=component, name=name, subdir=subdir, - base_mesh=base_mesh_step) + init_step = Init(component=component, name=name, + subdir=subdir, base_mesh=base_mesh_step) init_step.set_shared_config(config, link=config_filename) - self.add_step(init_step, symlink=symlink) - analysis_dependencies['init'][resolution] = init_step + analysis_dependencies['init'][refinement_factor] = init_step - name = f'{prefix}_forward_{mesh_name}' - subdir = f'{cos_bell_dir}/forward/{mesh_name}' - if self.include_viz: - symlink = f'forward/{mesh_name}' - else: - symlink = None + if resolution not in resolutions: + self.add_step(base_mesh_step, symlink=f'base_mesh/{mesh_name}') + self.add_step(init_step, symlink=f'init/{mesh_name}') + resolutions.append(resolution) + + timestep, _ = get_timestep_for_task( + config, refinement_factor, refinement=refinement) + timestep = ceil(timestep) + timesteps.append(timestep) + + subdir = f'{case_dir}/forward/{mesh_name}_{timestep}s' + symlink = f'forward/{mesh_name}_{timestep}s' if subdir in component.steps: forward_step = component.steps[subdir] else: - forward_step = Forward(component=component, name=name, - subdir=subdir, resolution=resolution, - mesh=base_mesh_step, - init=init_step) - forward_step.set_shared_config(config, link=config_filename) + name = f'{prefix}_forward_{mesh_name}_{timestep}s' + forward_step = Forward( + component=component, name=name, + subdir=subdir, + refinement_factor=refinement_factor, + mesh=base_mesh_step, + init=init_step, refinement=refinement) + forward_step.set_shared_config( + config, link=config_filename) self.add_step(forward_step, symlink=symlink) - analysis_dependencies['forward'][resolution] = forward_step + analysis_dependencies['forward'][refinement_factor] = forward_step if self.include_viz: - with_viz_dir = f'spherical/{prefix}/cosine_bell/with_viz' - name = f'{prefix}_viz_{mesh_name}' - subdir = f'{with_viz_dir}/viz/{mesh_name}' - step = Viz(component=component, name=name, - subdir=subdir, base_mesh=base_mesh_step, - init=init_step, forward=forward_step, - mesh_name=mesh_name) - step.set_shared_config(config, link=config_filename) - self.add_step(step) - - subdir = f'spherical/{prefix}/cosine_bell/analysis' - if self.include_viz: - symlink = 'analysis' - else: - symlink = None + name = f'{prefix}_viz_{mesh_name}_{timestep}s' + subdir = f'{case_dir}/viz/{mesh_name}_{timestep}s' + if subdir in component.steps: + viz_step = component.steps[subdir] + else: + viz_step = Viz(component=component, name=name, + subdir=subdir, base_mesh=base_mesh_step, + init=init_step, forward=forward_step, + mesh_name=mesh_name) + viz_step.set_shared_config(config, link=config_filename) + self.add_step(viz_step) + + subdir = f'{case_dir}/convergence_{refinement}/analysis' if subdir in component.steps: step = component.steps[subdir] step.resolutions = resolutions step.dependencies_dict = analysis_dependencies else: - step = Analysis(component=component, resolutions=resolutions, + step = Analysis(component=component, subdir=subdir, - dependencies=analysis_dependencies) + dependencies=analysis_dependencies, + refinement=refinement) step.set_shared_config(config, link=config_filename) - self.add_step(step, symlink=symlink) + if self.include_viz: + symlink = f'analysis_{refinement}' + self.add_step(step, symlink=symlink) + else: + self.add_step(step) diff --git a/polaris/ocean/tasks/cosine_bell/analysis.py b/polaris/ocean/tasks/cosine_bell/analysis.py index 7ba98f8ae..28fd7576e 100644 --- a/polaris/ocean/tasks/cosine_bell/analysis.py +++ b/polaris/ocean/tasks/cosine_bell/analysis.py @@ -3,7 +3,7 @@ from mpas_tools.transects import lon_lat_to_cartesian from mpas_tools.vector import Vector -from polaris.ocean.convergence import ConvergenceAnalysis +from polaris.ocean.convergence.analysis import ConvergenceAnalysis from polaris.ocean.tasks.cosine_bell.init import cosine_bell @@ -11,7 +11,7 @@ class Analysis(ConvergenceAnalysis): """ A step for analyzing the output from the cosine bell test case """ - def __init__(self, component, resolutions, subdir, dependencies): + def __init__(self, component, subdir, dependencies, refinement='both'): """ Create the step @@ -20,31 +20,31 @@ def __init__(self, component, resolutions, subdir, dependencies): component : polaris.Component The component the step belongs to - resolutions : list of float - The resolutions of the meshes that have been run - subdir : str The subdirectory that the step resides in dependencies : dict of dict of polaris.Steps The dependencies of this step + + refinement : str, optional + Whether to refine in space, time or both space and time """ convergence_vars = [{'name': 'tracer1', 'title': 'tracer1', 'zidx': 0}] super().__init__(component=component, subdir=subdir, - resolutions=resolutions, dependencies=dependencies, - convergence_vars=convergence_vars) + convergence_vars=convergence_vars, + refinement=refinement) - def exact_solution(self, mesh_name, field_name, time, zidx=None): + def exact_solution(self, refinement_factor, field_name, time, zidx=None): """ Get the exact solution Parameters ---------- - mesh_name : str - The mesh name which is the prefix for the initial condition file + refinement_factor : float + The factor by which to scale space, time or both field_name : str The name of the variable of which we evaluate convergence @@ -61,7 +61,7 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): Returns ------- - solution: xarray.DataArray + solution : xarray.DataArray The exact solution with dimension nCells """ @@ -76,10 +76,10 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): psi0 = config.getfloat('cosine_bell', 'psi0') vel_pd = config.getfloat('cosine_bell', 'vel_pd') - ds_mesh = xr.open_dataset(f'{mesh_name}_mesh.nc') + ds_mesh = xr.open_dataset(f'mesh_r{refinement_factor:02g}.nc') sphere_radius = ds_mesh.sphere_radius - ds_init = xr.open_dataset(f'{mesh_name}_init.nc') + ds_init = xr.open_dataset(f'init_r{refinement_factor:02g}.nc') latCell = ds_init.latCell.values lonCell = ds_init.lonCell.values diff --git a/polaris/ocean/tasks/cosine_bell/forward.py b/polaris/ocean/tasks/cosine_bell/forward.py index 19ae04ee0..8c477d9e7 100644 --- a/polaris/ocean/tasks/cosine_bell/forward.py +++ b/polaris/ocean/tasks/cosine_bell/forward.py @@ -7,7 +7,8 @@ class Forward(SphericalConvergenceForward): bell test case """ - def __init__(self, component, name, subdir, resolution, mesh, init): + def __init__(self, component, name, subdir, mesh, init, + refinement_factor, refinement='both'): """ Create a new step @@ -22,21 +23,26 @@ def __init__(self, component, name, subdir, resolution, mesh, init): subdir : str The subdirectory for the step - resolution : float - The resolution of the (uniform) mesh in km - mesh : polaris.Step The base mesh step init : polaris.Step The init step + + refinement_factor : float + The factor by which to scale space, time or both + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ package = 'polaris.ocean.tasks.cosine_bell' validate_vars = ['normalVelocity', 'tracer1'] super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, mesh=mesh, - init=init, package=package, + mesh=mesh, init=init, package=package, yaml_filename='forward.yaml', output_filename='output.nc', validate_vars=validate_vars, - graph_target=f'{mesh.path}/graph.info') + graph_target=f'{mesh.path}/graph.info', + refinement_factor=refinement_factor, + refinement=refinement) diff --git a/polaris/ocean/tasks/geostrophic/__init__.py b/polaris/ocean/tasks/geostrophic/__init__.py index d29ad36e7..9d67179ce 100644 --- a/polaris/ocean/tasks/geostrophic/__init__.py +++ b/polaris/ocean/tasks/geostrophic/__init__.py @@ -1,7 +1,12 @@ +from math import ceil from typing import Dict from polaris import Step, Task from polaris.config import PolarisConfigParser +from polaris.ocean.convergence import ( + get_resolution_for_task, + get_timestep_for_task, +) from polaris.ocean.mesh.spherical import add_spherical_base_mesh_step from polaris.ocean.tasks.geostrophic.analysis import Analysis from polaris.ocean.tasks.geostrophic.forward import Forward @@ -28,11 +33,13 @@ def add_geostrophic_tasks(component): config.add_from_package('polaris.ocean.tasks.geostrophic', 'geostrophic.cfg') - for include_viz in [False, True]: - component.add_task(Geostrophic(component=component, - config=config, - icosahedral=icosahedral, - include_viz=include_viz)) + for refinement in ['space', 'time', 'both']: + for include_viz in [False, True]: + component.add_task(Geostrophic(component=component, + config=config, + icosahedral=icosahedral, + include_viz=include_viz, + refinement=refinement)) class Geostrophic(Task): @@ -41,8 +48,9 @@ class Geostrophic(Task): Attributes ---------- - resolutions : list of float - A list of mesh resolutions + refinement : str + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time icosahedral : bool Whether to use icosahedral, as opposed to less regular, JIGSAW meshes @@ -50,7 +58,8 @@ class Geostrophic(Task): include_viz : bool Include VizMap and Viz steps for each resolution """ - def __init__(self, component, config, icosahedral, include_viz): + def __init__(self, component, config, icosahedral, include_viz, + refinement='both'): """ Create the convergence test @@ -68,29 +77,31 @@ def __init__(self, component, config, icosahedral, include_viz): include_viz : bool Include VizMap and Viz steps for each resolution + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ if icosahedral: prefix = 'icos' else: prefix = 'qu' - subdir = f'spherical/{prefix}/geostrophic' - name = f'{prefix}_geostrophic' + subdir = f'spherical/{prefix}/geostrophic/convergence_{refinement}' + name = f'{prefix}_geostrophic_convergence_{refinement}' if include_viz: subdir = f'{subdir}/with_viz' name = f'{name}_with_viz' - link = 'geostrophic.cfg' - else: - # config options live in the task already so no need for a symlink - link = None + link = 'geostrophic.cfg' + super().__init__(component=component, name=name, subdir=subdir) - self.resolutions = list() + self.refinement = refinement self.icosahedral = icosahedral self.include_viz = include_viz self.set_shared_config(config, link=link) - self._setup_steps() + self._setup_steps(refinement=refinement) def configure(self): """ @@ -99,96 +110,127 @@ def configure(self): super().configure() # set up the steps again in case a user has provided new resolutions - self._setup_steps() + self._setup_steps(self.refinement) - def _setup_steps(self): - """ setup steps given resolutions """ + def _setup_steps(self, refinement): + """ + setup steps given resolutions + + Parameters + ---------- + refinement : str + refinement type. One of 'space', 'time' or 'both' + """ icosahedral = self.icosahedral config = self.config config_filename = self.config_filename + if icosahedral: prefix = 'icos' else: prefix = 'qu' - resolutions = config.getlist('spherical_convergence', - f'{prefix}_resolutions', dtype=float) - - if self.resolutions == resolutions: - return + if refinement == 'time': + option = 'refinement_factors_time' + else: + option = 'refinement_factors_space' + refinement_factors = config.getlist('spherical_convergence', + f'{prefix}_{option}', dtype=str) + refinement_factors = ', '.join(refinement_factors) + config.set('convergence', option, value=refinement_factors) + refinement_factors = config.getlist('convergence', + option, dtype=float) + + base_resolution = config.getfloat('spherical_convergence', + f'{prefix}_base_resolution') + config.set('convergence', 'base_resolution', + value=f'{base_resolution:03g}') # start fresh with no steps for step in list(self.steps.values()): self.remove_step(step) - self.resolutions = resolutions - component = self.component + resolutions = list() + timesteps = list() + case_dir = f'spherical/{prefix}/geostrophic' analysis_dependencies: Dict[str, Dict[float, Step]] = ( dict(mesh=dict(), init=dict(), forward=dict())) - for resolution in resolutions: + + for refinement_factor in refinement_factors: + resolution = get_resolution_for_task( + config, refinement_factor, refinement=refinement) + base_mesh_step, mesh_name = add_spherical_base_mesh_step( component, resolution, icosahedral) - self.add_step(base_mesh_step, symlink=f'base_mesh/{mesh_name}') - analysis_dependencies['mesh'][resolution] = base_mesh_step + analysis_dependencies['mesh'][refinement_factor] = base_mesh_step name = f'{prefix}_init_{mesh_name}' subdir = f'{case_dir}/init/{mesh_name}' - if self.include_viz: - symlink = f'init/{mesh_name}' - else: - symlink = None if subdir in component.steps: init_step = component.steps[subdir] else: - init_step = Init(component=component, name=name, subdir=subdir, - base_mesh=base_mesh_step) + init_step = Init(component=component, name=name, + subdir=subdir, base_mesh=base_mesh_step) init_step.set_shared_config(config, link=config_filename) - self.add_step(init_step, symlink=symlink) - analysis_dependencies['init'][resolution] = init_step + analysis_dependencies['init'][refinement_factor] = init_step - name = f'{prefix}_forward_{mesh_name}' - subdir = f'{case_dir}/forward/{mesh_name}' - if self.include_viz: - symlink = f'forward/{mesh_name}' - else: - symlink = None + if resolution not in resolutions: + self.add_step(base_mesh_step, symlink=f'base_mesh/{mesh_name}') + self.add_step(init_step, symlink=f'init/{mesh_name}') + resolutions.append(resolution) + + timestep, _ = get_timestep_for_task( + config, refinement_factor, refinement=refinement) + timestep = ceil(timestep) + timesteps.append(timestep) + + subdir = f'{case_dir}/forward/{mesh_name}_{timestep}s' + symlink = f'forward/{mesh_name}_{timestep}s' if subdir in component.steps: forward_step = component.steps[subdir] else: - forward_step = Forward(component=component, name=name, - subdir=subdir, resolution=resolution, - mesh=base_mesh_step, - init=init_step) - forward_step.set_shared_config(config, link=config_filename) + name = f'{prefix}_forward_{mesh_name}_{timestep}s' + forward_step = Forward( + component=component, name=name, + subdir=subdir, + refinement_factor=refinement_factor, + mesh=base_mesh_step, + init=init_step, refinement=refinement) + forward_step.set_shared_config( + config, link=config_filename) self.add_step(forward_step, symlink=symlink) - analysis_dependencies['forward'][resolution] = forward_step + analysis_dependencies['forward'][refinement_factor] = forward_step if self.include_viz: - with_viz_dir = f'{case_dir}/with_viz' - name = f'{prefix}_viz_{mesh_name}' - subdir = f'{with_viz_dir}/viz/{mesh_name}' - step = Viz(component=component, name=name, - subdir=subdir, base_mesh=base_mesh_step, - init=init_step, forward=forward_step, - mesh_name=mesh_name) - step.set_shared_config(config, link=config_filename) - self.add_step(step) - - subdir = f'{case_dir}/analysis' - if self.include_viz: - symlink = 'analysis' - else: - symlink = None + name = f'{prefix}_viz_{mesh_name}_{timestep}s' + subdir = f'{case_dir}/viz/{mesh_name}_{timestep}s' + if subdir in component.steps: + viz_step = component.steps[subdir] + else: + viz_step = Viz(component=component, name=name, + subdir=subdir, base_mesh=base_mesh_step, + init=init_step, forward=forward_step, + mesh_name=mesh_name) + viz_step.set_shared_config(config, link=config_filename) + self.add_step(viz_step) + + subdir = f'{case_dir}/convergence_{refinement}/analysis' if subdir in component.steps: step = component.steps[subdir] step.resolutions = resolutions step.dependencies_dict = analysis_dependencies else: - step = Analysis(component=component, resolutions=resolutions, - subdir=subdir, dependencies=analysis_dependencies) + step = Analysis(component=component, + subdir=subdir, + dependencies=analysis_dependencies, + refinement=refinement) step.set_shared_config(config, link=config_filename) - self.add_step(step, symlink=symlink) + if self.include_viz: + symlink = f'analysis_{refinement}' + self.add_step(step, symlink=symlink) + else: + self.add_step(step) diff --git a/polaris/ocean/tasks/geostrophic/analysis.py b/polaris/ocean/tasks/geostrophic/analysis.py index f56847cb8..a26f01fdb 100644 --- a/polaris/ocean/tasks/geostrophic/analysis.py +++ b/polaris/ocean/tasks/geostrophic/analysis.py @@ -1,6 +1,6 @@ import xarray as xr -from polaris.ocean.convergence import ConvergenceAnalysis +from polaris.ocean.convergence.analysis import ConvergenceAnalysis from polaris.ocean.tasks.geostrophic.exact_solution import ( compute_exact_solution, ) @@ -10,7 +10,7 @@ class Analysis(ConvergenceAnalysis): """ A step for analyzing the output from the geostrophic test case """ - def __init__(self, component, resolutions, subdir, dependencies): + def __init__(self, component, subdir, dependencies, refinement='both'): """ Create the step @@ -19,14 +19,14 @@ def __init__(self, component, resolutions, subdir, dependencies): component : polaris.Component The component the step belongs to - resolutions : list of float - The resolutions of the meshes that have been run - subdir : str The subdirectory that the step resides in dependencies : dict of dict of polaris.Steps The dependencies of this step + + refinement : str, optional + Whether to refine in space, time or both space and time """ convergence_vars = [{'name': 'h', 'title': 'water-column thickness', @@ -35,18 +35,18 @@ def __init__(self, component, resolutions, subdir, dependencies): 'title': 'normal velocity', 'zidx': 0}] super().__init__(component=component, subdir=subdir, - resolutions=resolutions, dependencies=dependencies, - convergence_vars=convergence_vars) + convergence_vars=convergence_vars, + refinement=refinement) - def exact_solution(self, mesh_name, field_name, time, zidx=None): + def exact_solution(self, refinement_factor, field_name, time, zidx=None): """ Get the exact solution Parameters ---------- - mesh_name : str - The mesh name which is the prefix for the initial condition file + refinement_factor : float + The factor by which to scale space, time or both field_name : str The name of the variable of which we evaluate convergence @@ -78,7 +78,7 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): vel_period = section.getfloat('vel_period') gh_0 = section.getfloat('gh_0') - mesh_filename = f'{mesh_name}_mesh.nc' + mesh_filename = f'mesh_r{refinement_factor:02g}.nc' h, _, _, normalVelocity = compute_exact_solution( alpha, vel_period, gh_0, mesh_filename) @@ -88,7 +88,7 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): else: return normalVelocity - def get_output_field(self, mesh_name, field_name, time, zidx=None): + def get_output_field(self, refinement_factor, field_name, time, zidx=None): """ Get the model output field at the given time and z index @@ -117,15 +117,15 @@ def get_output_field(self, mesh_name, field_name, time, zidx=None): f'geostrophic test case') if field_name == 'normalVelocity': - return super().get_output_field(mesh_name=mesh_name, - field_name=field_name, - time=time, zidx=zidx) + return super().get_output_field( + refinement_factor=refinement_factor, + field_name=field_name, time=time, zidx=zidx) else: - ds_init = xr.open_dataset(f'{mesh_name}_init.nc') + ds_init = xr.open_dataset(f'init_r{refinement_factor:02g}.nc') bottom_depth = ds_init.bottomDepth - ssh = super().get_output_field(mesh_name=mesh_name, - field_name='ssh', - time=time, zidx=None) + ssh = super().get_output_field( + refinement_factor=refinement_factor, + field_name='ssh', time=time, zidx=None) h = ssh + bottom_depth return h diff --git a/polaris/ocean/tasks/geostrophic/forward.py b/polaris/ocean/tasks/geostrophic/forward.py index 1e788b468..8e8188ea8 100644 --- a/polaris/ocean/tasks/geostrophic/forward.py +++ b/polaris/ocean/tasks/geostrophic/forward.py @@ -7,7 +7,8 @@ class Forward(SphericalConvergenceForward): bell test case """ - def __init__(self, component, name, subdir, resolution, mesh, init): + def __init__(self, component, name, subdir, mesh, init, + refinement_factor, refinement='both'): """ Create a new step @@ -22,14 +23,18 @@ def __init__(self, component, name, subdir, resolution, mesh, init): subdir : str The subdirectory for the step - resolution : float - The resolution of the (uniform) mesh in km - mesh : polaris.Step The base mesh step init : polaris.Step The init step + + refinement_factor : float + The factor by which to scale space, time or both + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ package = 'polaris.ocean.tasks.geostrophic' validate_vars = ['temperature', @@ -37,9 +42,10 @@ def __init__(self, component, name, subdir, resolution, mesh, init): 'layerThickness', 'normalVelocity'] super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, mesh=mesh, - init=init, package=package, + mesh=mesh, init=init, package=package, yaml_filename='forward.yaml', output_filename='output.nc', validate_vars=validate_vars, - graph_target=f'{mesh.path}/graph.info') + graph_target=f'{mesh.path}/graph.info', + refinement_factor=refinement_factor, + refinement=refinement) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py index a3023a74e..f9b2de3a0 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/__init__.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/__init__.py @@ -1,6 +1,12 @@ +from math import ceil from typing import Dict from polaris import Step, Task +from polaris.config import PolarisConfigParser +from polaris.ocean.convergence import ( + get_resolution_for_task, + get_timestep_for_task, +) from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.inertial_gravity_wave.analysis import Analysis from polaris.ocean.tasks.inertial_gravity_wave.forward import Forward @@ -15,7 +21,19 @@ def add_inertial_gravity_wave_tasks(component): component : polaris.ocean.Ocean the ocean component that the task will be added to """ - component.add_task(InertialGravityWave(component=component)) + basedir = 'planar/inertial_gravity_wave' + config_filename = 'inertial_gravity_wave.cfg' + filepath = f'{basedir}/{config_filename}' + config = PolarisConfigParser(filepath=filepath) + config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') + config.add_from_package( + 'polaris.ocean.tasks.inertial_gravity_wave', + config_filename) + for refinement in ['space', 'time', 'both']: + component.add_task(InertialGravityWave(component=component, + config=config, + refinement=refinement)) class InertialGravityWave(Task): @@ -23,7 +41,7 @@ class InertialGravityWave(Task): The convergence test case for inertial gravity waves """ - def __init__(self, component): + def __init__(self, component, config, refinement='both'): """ Create the test case @@ -31,39 +49,82 @@ def __init__(self, component): ---------- component : polaris.ocean.Ocean The ocean component that this task belongs to + + config : polaris.config.PolarisConfigParser + A shared config parser + + refinement : str, optional + Whether to refine in space, time or both space and time """ - name = 'inertial_gravity_wave' - subdir = f'planar/{name}' + name = f'inertial_gravity_wave_convergence_{refinement}' + basedir = 'planar/inertial_gravity_wave' + subdir = f'{basedir}/convergence_{refinement}' + config_filename = 'inertial_gravity_wave.cfg' super().__init__(component=component, name=name, subdir=subdir) + self.set_shared_config(config, link=config_filename) - self.resolutions = [200., 100., 50., 25.] analysis_dependencies: Dict[str, Dict[float, Step]] = ( dict(mesh=dict(), init=dict(), forward=dict())) - for resolution in self.resolutions: + + if refinement == 'time': + option = 'refinement_factors_time' + else: + option = 'refinement_factors_space' + refinement_factors = self.config.getlist( + 'convergence', option, dtype=float) + timesteps = list() + resolutions = list() + for refinement_factor in refinement_factors: + resolution = get_resolution_for_task( + self.config, refinement_factor, refinement=refinement) mesh_name = resolution_to_subdir(resolution) - init_step = Init(component=component, resolution=resolution, - taskdir=self.subdir) - self.add_step(init_step) - forward_step = Forward(component=component, - name=f'forward_{mesh_name}', - resolution=resolution, - subdir=f'{self.subdir}/forward/{mesh_name}', - init=init_step) - self.add_step(forward_step) - analysis_dependencies['mesh'][resolution] = init_step - analysis_dependencies['init'][resolution] = init_step - analysis_dependencies['forward'][resolution] = forward_step + + subdir = f'{basedir}/init/{mesh_name}' + symlink = f'init/{mesh_name}' + if subdir in component.steps: + init_step = component.steps[subdir] + else: + init_step = Init(component=component, resolution=resolution, + subdir=subdir) + init_step.set_shared_config(self.config, link=config_filename) + if resolution not in resolutions: + self.add_step(init_step, symlink=symlink) + resolutions.append(resolution) + + timestep, _ = get_timestep_for_task( + config, refinement_factor, refinement=refinement) + timestep = ceil(timestep) + timesteps.append(timestep) + + subdir = f'{basedir}/forward/{mesh_name}_{timestep}s' + symlink = f'forward/{mesh_name}_{timestep}s' + if subdir in component.steps: + forward_step = component.steps[subdir] + else: + forward_step = Forward( + component=component, + refinement=refinement, + refinement_factor=refinement_factor, + name=f'forward_{mesh_name}_{timestep}s', + subdir=subdir, + init=init_step) + forward_step.set_shared_config( + config, link=config_filename) + self.add_step(forward_step, symlink=symlink) + + analysis_dependencies['mesh'][refinement_factor] = init_step + analysis_dependencies['init'][refinement_factor] = init_step + analysis_dependencies['forward'][refinement_factor] = forward_step self.add_step(Analysis(component=component, - resolutions=self.resolutions, - subdir=f'{subdir}/analysis', + subdir=f'{self.subdir}/analysis', + refinement=refinement, dependencies=analysis_dependencies)) - self.add_step(Viz(component=component, resolutions=self.resolutions, + self.add_step(Viz(component=component, resolutions=resolutions, taskdir=self.subdir), run_by_default=False) - - self.config.add_from_package('polaris.ocean.convergence', - 'convergence.cfg') - self.config.add_from_package( + config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') + config.add_from_package( 'polaris.ocean.tasks.inertial_gravity_wave', - 'inertial_gravity_wave.cfg') + config_filename) diff --git a/polaris/ocean/tasks/inertial_gravity_wave/analysis.py b/polaris/ocean/tasks/inertial_gravity_wave/analysis.py index c9383358e..4cdb4e59a 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/analysis.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/analysis.py @@ -1,6 +1,6 @@ import xarray as xr -from polaris.ocean.convergence import ConvergenceAnalysis +from polaris.ocean.convergence.analysis import ConvergenceAnalysis from polaris.ocean.tasks.inertial_gravity_wave.exact_solution import ( ExactSolution, ) @@ -16,7 +16,7 @@ class Analysis(ConvergenceAnalysis): resolutions : list of float The resolutions of the meshes that have been run """ - def __init__(self, component, resolutions, subdir, dependencies): + def __init__(self, component, subdir, dependencies, refinement='both'): """ Create the step @@ -33,23 +33,26 @@ def __init__(self, component, resolutions, subdir, dependencies): dependencies : dict of dict of polaris.Steps The dependencies of this step + + refinement : str, optional + Whether to refine in space, time or both space and time """ convergence_vars = [{'name': 'ssh', 'title': 'SSH', 'zidx': None}] super().__init__(component=component, subdir=subdir, - resolutions=resolutions, dependencies=dependencies, - convergence_vars=convergence_vars) + convergence_vars=convergence_vars, + refinement=refinement) - def exact_solution(self, mesh_name, field_name, time, zidx=None): + def exact_solution(self, refinement_factor, field_name, time, zidx=None): """ Get the exact solution Parameters ---------- - mesh_name : str - The mesh name which is the prefix for the initial condition file + refinement_factor : float + The factor by which to scale space, time or both field_name : str The name of the variable of which we evaluate convergence @@ -69,7 +72,7 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): solution : xarray.DataArray The exact solution as derived from the initial condition """ - init = xr.open_dataset(f'{mesh_name}_init.nc') + init = xr.open_dataset(f'init_r{refinement_factor:02g}.nc') exact = ExactSolution(init, self.config) if field_name != 'ssh': raise ValueError(f'{field_name} is not currently supported') diff --git a/polaris/ocean/tasks/inertial_gravity_wave/forward.py b/polaris/ocean/tasks/inertial_gravity_wave/forward.py index abca61c8f..68930dde6 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/forward.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/forward.py @@ -1,7 +1,8 @@ import numpy as np from polaris.mesh.planar import compute_planar_hex_nx_ny -from polaris.ocean.convergence import ConvergenceForward +from polaris.ocean.convergence import get_resolution_for_task +from polaris.ocean.convergence.forward import ConvergenceForward class Forward(ConvergenceForward): @@ -11,10 +12,15 @@ class Forward(ConvergenceForward): Attributes ---------- - resolution : float - The resolution of the test case in km + refinement_factor : float + The factor by which to scale space, time or both + + refinement : str + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ - def __init__(self, component, name, resolution, subdir, init): + def __init__(self, component, name, refinement_factor, subdir, + init, refinement='both'): """ Create a new test case @@ -23,18 +29,26 @@ def __init__(self, component, name, resolution, subdir, init): component : polaris.Component The component the step belongs to - resolution : km - The resolution of the test case in km + name : str + The name of the step + refinement_factor : float + The factor by which to scale space, time or both subdir : str - The subdirectory that the step belongs to + The subdirectory that the task belongs to init : polaris.Step The step which generates the mesh and initial condition + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, mesh=init, init=init, + refinement=refinement, + refinement_factor=refinement_factor, + mesh=init, init=init, package='polaris.ocean.tasks.inertial_gravity_wave', yaml_filename='forward.yaml', graph_target=f'{init.path}/culled_graph.info', @@ -52,8 +66,10 @@ def compute_cell_count(self): The approximate number of cells in the mesh """ section = self.config['inertial_gravity_wave'] + resolution = get_resolution_for_task( + self.config, self.refinement_factor, refinement=self.refinement) lx = section.getfloat('lx') ly = np.sqrt(3.0) / 2.0 * lx - nx, ny = compute_planar_hex_nx_ny(lx, ly, self.resolution) + nx, ny = compute_planar_hex_nx_ny(lx, ly, resolution) cell_count = nx * ny return cell_count diff --git a/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg b/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg index c33e82230..aebb28655 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg +++ b/polaris/ocean/tasks/inertial_gravity_wave/inertial_gravity_wave.cfg @@ -57,6 +57,13 @@ convergence_thresh = ${inertial_gravity_wave:conv_thresh} # Type of error to compute error_type = l2 +# the base mesh resolution (km) to which refinement_factors +# are applied if refinement is 'space' or 'both' on a planar mesh +base_resolution = 100. + +# refinement factors for a planar mesh applied to either space or time +refinement_factors_space = 2., 1., 0.5, 0.25 +refinement_factors_time = 1., 0.5, 0.25 # config options for spherical convergence tests [convergence_forward] diff --git a/polaris/ocean/tasks/inertial_gravity_wave/init.py b/polaris/ocean/tasks/inertial_gravity_wave/init.py index 4935e912a..f19d89563 100644 --- a/polaris/ocean/tasks/inertial_gravity_wave/init.py +++ b/polaris/ocean/tasks/inertial_gravity_wave/init.py @@ -23,7 +23,7 @@ class Init(Step): resolution : float The resolution of the test case in km """ - def __init__(self, component, resolution, taskdir): + def __init__(self, component, resolution, subdir): """ Create the step @@ -35,13 +35,13 @@ def __init__(self, component, resolution, taskdir): resolution : float The resolution of the test case in km - taskdir : str + subdir : str The subdirectory that the task belongs to """ mesh_name = resolution_to_subdir(resolution) super().__init__(component=component, name=f'init_{mesh_name}', - subdir=f'{taskdir}/init/{mesh_name}') + subdir=subdir) self.resolution = resolution for filename in ['culled_mesh.nc', 'initial_state.nc', 'culled_graph.info']: diff --git a/polaris/ocean/tasks/manufactured_solution/__init__.py b/polaris/ocean/tasks/manufactured_solution/__init__.py index 2c95b634e..f85917698 100644 --- a/polaris/ocean/tasks/manufactured_solution/__init__.py +++ b/polaris/ocean/tasks/manufactured_solution/__init__.py @@ -1,6 +1,12 @@ +from math import ceil from typing import Dict from polaris import Step, Task +from polaris.config import PolarisConfigParser +from polaris.ocean.convergence import ( + get_resolution_for_task, + get_timestep_for_task, +) from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.manufactured_solution.analysis import Analysis from polaris.ocean.tasks.manufactured_solution.forward import Forward @@ -16,19 +22,27 @@ def add_manufactured_solution_tasks(component): component : polaris.ocean.Ocean the ocean component that the task will be added to """ - component.add_task(ManufacturedSolution(component=component)) + basedir = 'planar/manufactured_solution' + config_filename = 'manufactured_solution.cfg' + filepath = f'{basedir}/{config_filename}' + config = PolarisConfigParser(filepath=filepath) + config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') + config.add_from_package( + 'polaris.ocean.tasks.manufactured_solution', + config_filename) + for refinement in ['space', 'time', 'both']: + component.add_task(ManufacturedSolution(component=component, + config=config, + refinement=refinement)) class ManufacturedSolution(Task): """ The convergence test case using the method of manufactured solutions - - Attributes - ---------- - resolutions : list of floats - The resolutions of the test case in km """ - def __init__(self, component): + + def __init__(self, component, config, refinement='both'): """ Create the test case @@ -36,38 +50,82 @@ def __init__(self, component): ---------- component : polaris.ocean.Ocean The ocean component that this task belongs to + + config : polaris.config.PolarisConfigParser + A shared config parser + + refinement : str, optional + Whether to refine in space, time or both space and time """ - name = 'manufactured_solution' - subdir = f'planar/{name}' + name = f'manufactured_solution_convergence_{refinement}' + basedir = 'planar/manufactured_solution' + subdir = f'{basedir}/convergence_{refinement}' + config_filename = 'manufactured_solution.cfg' super().__init__(component=component, name=name, subdir=subdir) + self.set_shared_config(config, link=config_filename) - self.resolutions = [200., 100., 50., 25.] analysis_dependencies: Dict[str, Dict[float, Step]] = ( dict(mesh=dict(), init=dict(), forward=dict())) - for resolution in self.resolutions: + + if refinement == 'time': + option = 'refinement_factors_time' + else: + option = 'refinement_factors_space' + refinement_factors = self.config.getlist( + 'convergence', option, dtype=float) + timesteps = list() + resolutions = list() + for refinement_factor in refinement_factors: + resolution = get_resolution_for_task( + self.config, refinement_factor, refinement=refinement) mesh_name = resolution_to_subdir(resolution) - init_step = Init(component=component, resolution=resolution, - taskdir=self.subdir) - self.add_step(init_step) - forward_step = Forward(component=component, resolution=resolution, - name=f'forward_{mesh_name}', - subdir=f'{self.subdir}/forward/{mesh_name}', - init=init_step) - self.add_step(forward_step) - analysis_dependencies['mesh'][resolution] = init_step - analysis_dependencies['init'][resolution] = init_step - analysis_dependencies['forward'][resolution] = forward_step + + subdir = f'{basedir}/init/{mesh_name}' + symlink = f'init/{mesh_name}' + if subdir in component.steps: + init_step = component.steps[subdir] + else: + init_step = Init(component=component, resolution=resolution, + subdir=subdir) + init_step.set_shared_config(self.config, link=config_filename) + if resolution not in resolutions: + self.add_step(init_step, symlink=symlink) + resolutions.append(resolution) + + timestep, _ = get_timestep_for_task( + config, refinement_factor, refinement=refinement) + timestep = ceil(timestep) + timesteps.append(timestep) + + subdir = f'{basedir}/forward/{mesh_name}_{timestep}s' + symlink = f'forward/{mesh_name}_{timestep}s' + if subdir in component.steps: + forward_step = component.steps[subdir] + else: + forward_step = Forward( + component=component, + refinement=refinement, + refinement_factor=refinement_factor, + name=f'forward_{mesh_name}_{timestep}s', + subdir=subdir, + init=init_step) + forward_step.set_shared_config( + config, link=config_filename) + self.add_step(forward_step, symlink=symlink) + + analysis_dependencies['mesh'][refinement_factor] = init_step + analysis_dependencies['init'][refinement_factor] = init_step + analysis_dependencies['forward'][refinement_factor] = forward_step self.add_step(Analysis(component=component, - resolutions=self.resolutions, subdir=f'{self.subdir}/analysis', + refinement=refinement, dependencies=analysis_dependencies)) - self.add_step(Viz(component=component, resolutions=self.resolutions, + self.add_step(Viz(component=component, resolutions=resolutions, taskdir=self.subdir), run_by_default=False) - - self.config.add_from_package('polaris.ocean.convergence', - 'convergence.cfg') - self.config.add_from_package( + config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') + config.add_from_package( 'polaris.ocean.tasks.manufactured_solution', - 'manufactured_solution.cfg') + config_filename) diff --git a/polaris/ocean/tasks/manufactured_solution/analysis.py b/polaris/ocean/tasks/manufactured_solution/analysis.py index e7947192e..06ba74046 100644 --- a/polaris/ocean/tasks/manufactured_solution/analysis.py +++ b/polaris/ocean/tasks/manufactured_solution/analysis.py @@ -1,4 +1,4 @@ -from polaris.ocean.convergence import ConvergenceAnalysis +from polaris.ocean.convergence.analysis import ConvergenceAnalysis from polaris.ocean.tasks.manufactured_solution.exact_solution import ( ExactSolution, ) @@ -14,7 +14,7 @@ class Analysis(ConvergenceAnalysis): resolutions : list of float The resolutions of the meshes that have been run """ - def __init__(self, component, resolutions, subdir, dependencies): + def __init__(self, component, subdir, dependencies, refinement='both'): """ Create the step @@ -31,23 +31,26 @@ def __init__(self, component, resolutions, subdir, dependencies): dependencies : dict of dict of polaris.Steps The dependencies of this step + + refinement : str, optional + Whether to refine in space, time or both space and time """ convergence_vars = [{'name': 'ssh', 'title': 'SSH', 'zidx': None}] super().__init__(component=component, subdir=subdir, - resolutions=resolutions, dependencies=dependencies, - convergence_vars=convergence_vars) + convergence_vars=convergence_vars, + refinement=refinement) - def exact_solution(self, mesh_name, field_name, time, zidx=None): + def exact_solution(self, refinement_factor, field_name, time, zidx=None): """ Get the exact solution Parameters ---------- - mesh_name : str - The mesh name which is the prefix for the initial condition file + refinement_factor : float + The factor by which to scale space, time or both field_name : str The name of the variable of which we evaluate convergence @@ -67,7 +70,7 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None): solution : xarray.DataArray The exact solution as derived from the initial condition """ - init = self.open_model_dataset(f'{mesh_name}_init.nc') + init = self.open_model_dataset(f'init_r{refinement_factor:02g}.nc') exact = ExactSolution(self.config, init) if field_name != 'ssh': raise ValueError(f'{field_name} is not currently supported') diff --git a/polaris/ocean/tasks/manufactured_solution/forward.py b/polaris/ocean/tasks/manufactured_solution/forward.py index 62bd76a4d..012bea98f 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.py +++ b/polaris/ocean/tasks/manufactured_solution/forward.py @@ -1,7 +1,8 @@ import numpy as np from polaris.mesh.planar import compute_planar_hex_nx_ny -from polaris.ocean.convergence import ConvergenceForward +from polaris.ocean.convergence import get_resolution_for_task +from polaris.ocean.convergence.forward import ConvergenceForward from polaris.ocean.tasks.manufactured_solution.exact_solution import ( ExactSolution, ) @@ -14,10 +15,15 @@ class Forward(ConvergenceForward): Attributes ---------- - resolution : float - The resolution of the test case in km + refinement_factor : float + The factor by which to scale space, time or both + + refinement : str + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ - def __init__(self, component, name, resolution, subdir, init): + def __init__(self, component, name, refinement_factor, subdir, + init, refinement='both'): """ Create a new test case @@ -26,18 +32,26 @@ def __init__(self, component, name, resolution, subdir, init): component : polaris.Component The component the step belongs to - resolution : km - The resolution of the test case in km + name : str + The name of the step + + refinement_factor : float + The factor by which to scale space, time or both subdir : str The subdirectory that the task belongs to init : polaris.Step The step which generates the mesh and initial condition + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, mesh=init, init=init, + refinement_factor=refinement_factor, + mesh=init, init=init, refinement=refinement, package='polaris.ocean.tasks.manufactured_solution', yaml_filename='forward.yaml', graph_target=f'{init.path}/culled_graph.info', @@ -66,10 +80,12 @@ def compute_cell_count(self): The approximate number of cells in the mesh """ # no file to read from, so we'll compute it based on config options + resolution = get_resolution_for_task( + self.config, self.refinement_factor, refinement=self.refinement) section = self.config['manufactured_solution'] lx = section.getfloat('lx') ly = np.sqrt(3.0) / 2.0 * lx - nx, ny = compute_planar_hex_nx_ny(lx, ly, self.resolution) + nx, ny = compute_planar_hex_nx_ny(lx, ly, resolution) cell_count = nx * ny return cell_count diff --git a/polaris/ocean/tasks/manufactured_solution/init.py b/polaris/ocean/tasks/manufactured_solution/init.py index 63cf166f2..d08cf15fc 100644 --- a/polaris/ocean/tasks/manufactured_solution/init.py +++ b/polaris/ocean/tasks/manufactured_solution/init.py @@ -22,7 +22,7 @@ class Init(OceanIOStep): resolution : float The resolution of the test case in km """ - def __init__(self, component, resolution, taskdir): + def __init__(self, component, resolution, subdir): """ Create the step @@ -40,7 +40,7 @@ def __init__(self, component, resolution, taskdir): mesh_name = resolution_to_subdir(resolution) super().__init__(component=component, name=f'init_{mesh_name}', - subdir=f'{taskdir}/init/{mesh_name}') + subdir=subdir) self.resolution = resolution def setup(self): diff --git a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg index 079efdd92..a4edc7e49 100644 --- a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg +++ b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg @@ -65,6 +65,13 @@ convergence_thresh = ${manufactured_solution:conv_thresh} # Type of error to compute error_type = l2 +# the base mesh resolution (km) to which refinement_factors +# are applied if refinement is 'space' or 'both' on a planar mesh +base_resolution = 50. + +# refinement factors for a planar mesh applied to either space or time +refinement_factors_space = 4., 2., 1., 0.5 +refinement_factors_time = 1., 0.5, 0.25 # config options for spherical convergence tests [convergence_forward] diff --git a/polaris/ocean/tasks/sphere_transport/__init__.py b/polaris/ocean/tasks/sphere_transport/__init__.py index 972199fad..599bedf7b 100644 --- a/polaris/ocean/tasks/sphere_transport/__init__.py +++ b/polaris/ocean/tasks/sphere_transport/__init__.py @@ -1,7 +1,12 @@ +from math import ceil from typing import Dict from polaris import Step, Task from polaris.config import PolarisConfigParser +from polaris.ocean.convergence import ( + get_resolution_for_task, + get_timestep_for_task, +) from polaris.ocean.mesh.spherical import add_spherical_base_mesh_step from polaris.ocean.tasks.sphere_transport.analysis import Analysis from polaris.ocean.tasks.sphere_transport.filament_analysis import ( @@ -50,13 +55,18 @@ class SphereTransport(Task): resolutions : list of float A list of mesh resolutions + refinement : str + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time + icosahedral : bool Whether to use icosahedral, as opposed to less regular, JIGSAW meshes include_viz : bool Include VizMap and Viz steps for each resolution """ - def __init__(self, component, config, case_name, icosahedral, include_viz): + def __init__(self, component, config, case_name, icosahedral, include_viz, + refinement='both'): """ Create test case for creating a global MPAS-Ocean mesh @@ -78,6 +88,10 @@ def __init__(self, component, config, case_name, icosahedral, include_viz): include_viz : bool Include VizMap and Viz steps for each resolution + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ if icosahedral: prefix = 'icos' @@ -95,6 +109,7 @@ def __init__(self, component, config, case_name, icosahedral, include_viz): link = None super().__init__(component=component, name=name, subdir=subdir) self.resolutions = list() + self.refinement = refinement self.icosahedral = icosahedral self.case_name = case_name self.include_viz = include_viz @@ -102,7 +117,7 @@ def __init__(self, component, config, case_name, icosahedral, include_viz): self.set_shared_config(config, link=link) # add the steps with default resolutions so they can be listed - self._setup_steps() + self._setup_steps(refinement=refinement) def configure(self): """ @@ -111,9 +126,9 @@ def configure(self): super().configure() # set up the steps again in case a user has provided new resolutions - self._setup_steps() + self._setup_steps(self.refinement) - def _setup_steps(self): + def _setup_steps(self, refinement): # noqa: C901 """ setup steps given resolutions """ case_name = self.case_name icosahedral = self.icosahedral @@ -125,17 +140,27 @@ def _setup_steps(self): else: prefix = 'qu' - resolutions = config.getlist('spherical_convergence', - f'{prefix}_resolutions', dtype=float) - - if self.resolutions == resolutions: - return + if refinement == 'time': + option = 'refinement_factors_time' + else: + option = 'refinement_factors_space' + refinement_factors = config.getlist('spherical_convergence', + f'{prefix}_{option}', dtype=str) + refinement_factors = ', '.join(refinement_factors) + config.set('convergence', option, value=refinement_factors) + refinement_factors = config.getlist('convergence', + option, dtype=float) + + base_resolution = config.getfloat('spherical_convergence', + f'{prefix}_base_resolution') + config.set('convergence', 'base_resolution', + value=f'{base_resolution:03g}') # start fresh with no steps for step in list(self.steps.values()): self.remove_step(step) - self.resolutions = resolutions + resolutions = self.resolutions component = self.component @@ -143,11 +168,13 @@ def _setup_steps(self): analysis_dependencies: Dict[str, Dict[str, Step]] = ( dict(mesh=dict(), init=dict(), forward=dict())) - for resolution in resolutions: + timesteps = list() + for idx, refinement_factor in enumerate(refinement_factors): + resolution = get_resolution_for_task( + config, refinement_factor, refinement=refinement) base_mesh_step, mesh_name = add_spherical_base_mesh_step( component, resolution, icosahedral) - self.add_step(base_mesh_step, symlink=f'base_mesh/{mesh_name}') - analysis_dependencies['mesh'][resolution] = base_mesh_step + analysis_dependencies['mesh'][refinement_factor] = base_mesh_step name = f'{prefix}_init_{mesh_name}' subdir = f'{sph_trans_dir}/init/{mesh_name}' @@ -158,34 +185,45 @@ def _setup_steps(self): if subdir in component.steps: init_step = component.steps[subdir] else: - init_step = Init(component=component, name=name, subdir=subdir, - base_mesh=base_mesh_step, case_name=case_name) + init_step = Init(component=component, name=name, + subdir=subdir, base_mesh=base_mesh_step, + case_name=case_name) init_step.set_shared_config(config, link=config_filename) + analysis_dependencies['init'][refinement_factor] = init_step + self.add_step(base_mesh_step, symlink=f'base_mesh/{mesh_name}') self.add_step(init_step, symlink=symlink) - analysis_dependencies['init'][resolution] = init_step - - name = f'{prefix}_forward_{mesh_name}' - subdir = f'{sph_trans_dir}/forward/{mesh_name}' + if resolution not in resolutions: + resolutions.append(resolution) + + timestep, _ = get_timestep_for_task( + config, refinement_factor, refinement=refinement) + timestep = ceil(timestep) + timesteps.append(timestep) + name = f'{prefix}_forward_{mesh_name}_{timestep}s' + subdir = f'{sph_trans_dir}/forward/{mesh_name}_{timestep}s' if self.include_viz: - symlink = f'forward/{mesh_name}' + symlink = f'forward/{mesh_name}_{timestep}s' else: symlink = None if subdir in component.steps: forward_step = component.steps[subdir] else: - forward_step = Forward(component=component, name=name, - subdir=subdir, resolution=resolution, - base_mesh=base_mesh_step, - init=init_step, - case_name=case_name) + forward_step = Forward( + component=component, name=name, + subdir=subdir, + base_mesh=base_mesh_step, + init=init_step, + case_name=case_name, + refinement_factor=refinement_factor, + refinement=refinement) forward_step.set_shared_config(config, link=config_filename) self.add_step(forward_step, symlink=symlink) - analysis_dependencies['forward'][resolution] = forward_step + analysis_dependencies['forward'][refinement_factor] = forward_step if self.include_viz: with_viz_dir = f'{sph_trans_dir}/with_viz' - name = f'{prefix}_viz_{mesh_name}' - subdir = f'{with_viz_dir}/viz/{mesh_name}' + name = f'{prefix}_viz_{mesh_name}_{timestep}s' + subdir = f'{with_viz_dir}/viz/{mesh_name}_{timestep}s' step = Viz(component=component, name=name, subdir=subdir, base_mesh=base_mesh_step, init=init_step, forward=forward_step, @@ -193,9 +231,9 @@ def _setup_steps(self): step.set_shared_config(config, link=config_filename) self.add_step(step) - subdir = f'{sph_trans_dir}/analysis' + subdir = f'{sph_trans_dir}/analysis/{refinement}' if self.include_viz: - symlink = 'analysis' + symlink = f'analysis_{refinement}' else: symlink = None if subdir in component.steps: @@ -203,9 +241,10 @@ def _setup_steps(self): step.resolutions = resolutions step.dependencies_dict = analysis_dependencies else: - step = Analysis(component=component, resolutions=resolutions, + step = Analysis(component=component, subdir=subdir, case_name=case_name, - dependencies=analysis_dependencies) + dependencies=analysis_dependencies, + refinement=refinement) step.set_shared_config(config, link=config_filename) self.add_step(step, symlink=symlink) @@ -219,10 +258,11 @@ def _setup_steps(self): step = component.steps[subdir] else: step = MixingAnalysis(component=component, - resolutions=resolutions, icosahedral=icosahedral, subdir=subdir, + refinement_factors=refinement_factors, case_name=case_name, - dependencies=analysis_dependencies) + dependencies=analysis_dependencies, + refinement=refinement) step.set_shared_config(config, link=config_filename) self.add_step(step, symlink=symlink) @@ -236,9 +276,10 @@ def _setup_steps(self): step = component.steps[subdir] else: step = FilamentAnalysis(component=component, - resolutions=resolutions, + refinement_factors=refinement_factors, icosahedral=icosahedral, subdir=subdir, case_name=case_name, - dependencies=analysis_dependencies) + dependencies=analysis_dependencies, + refinement=refinement) step.set_shared_config(config, link=config_filename) self.add_step(step, symlink=symlink) diff --git a/polaris/ocean/tasks/sphere_transport/analysis.py b/polaris/ocean/tasks/sphere_transport/analysis.py index 620a80651..17983cce3 100644 --- a/polaris/ocean/tasks/sphere_transport/analysis.py +++ b/polaris/ocean/tasks/sphere_transport/analysis.py @@ -1,4 +1,4 @@ -from polaris.ocean.convergence import ConvergenceAnalysis +from polaris.ocean.convergence.analysis import ConvergenceAnalysis class Analysis(ConvergenceAnalysis): @@ -13,8 +13,8 @@ class Analysis(ConvergenceAnalysis): case_name : str The name of the test case """ - def __init__(self, component, resolutions, subdir, case_name, - dependencies): + def __init__(self, component, subdir, case_name, + dependencies, refinement='both'): """ Create the step @@ -23,9 +23,6 @@ def __init__(self, component, resolutions, subdir, case_name, component : polaris.Component The component the step belongs to - resolutions : list of float - The resolutions of the meshes that have been run - subdir : str The subdirectory that the step resides in @@ -34,6 +31,10 @@ def __init__(self, component, resolutions, subdir, case_name, dependencies : dict of dict of polaris.Steps The dependencies of this step + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ self.case_name = case_name convergence_vars = [{'name': 'tracer1', @@ -46,9 +47,9 @@ def __init__(self, component, resolutions, subdir, case_name, 'title': 'tracer3', 'zidx': 1}] super().__init__(component=component, subdir=subdir, - resolutions=resolutions, dependencies=dependencies, - convergence_vars=convergence_vars) + convergence_vars=convergence_vars, + refinement=refinement) # Note: there is no need to overwrite the default method exact_solution # which uses the initial condition diff --git a/polaris/ocean/tasks/sphere_transport/filament_analysis.py b/polaris/ocean/tasks/sphere_transport/filament_analysis.py index a3f196439..f5a3d6895 100644 --- a/polaris/ocean/tasks/sphere_transport/filament_analysis.py +++ b/polaris/ocean/tasks/sphere_transport/filament_analysis.py @@ -5,6 +5,7 @@ from polaris import Step from polaris.mpas import time_index_from_xtime +from polaris.ocean.convergence import get_resolution_for_task from polaris.ocean.resolution import resolution_to_subdir from polaris.viz import use_mplstyle @@ -25,8 +26,8 @@ class FilamentAnalysis(Step): case_name : str The name of the test case """ - def __init__(self, component, resolutions, icosahedral, subdir, - case_name, dependencies): + def __init__(self, component, refinement_factors, icosahedral, subdir, + case_name, dependencies, refinement='both'): """ Create the step @@ -53,22 +54,22 @@ def __init__(self, component, resolutions, icosahedral, subdir, """ super().__init__(component=component, name='filament_analysis', subdir=subdir) - self.resolutions = resolutions + self.refinement_factors = refinement_factors + self.refinement = refinement self.case_name = case_name - for resolution in resolutions: - mesh_name = resolution_to_subdir(resolution) - base_mesh = dependencies['mesh'][resolution] - init = dependencies['init'][resolution] - forward = dependencies['forward'][resolution] + for refinement_factor in refinement_factors: + base_mesh = dependencies['mesh'][refinement_factor] + init = dependencies['init'][refinement_factor] + forward = dependencies['forward'][refinement_factor] self.add_input_file( - filename=f'{mesh_name}_mesh.nc', + filename=f'mesh_r{refinement_factor:02g}.nc', work_dir_target=f'{base_mesh.path}/base_mesh.nc') self.add_input_file( - filename=f'{mesh_name}_init.nc', + filename=f'init_r{refinement_factor:02g}.nc', work_dir_target=f'{init.path}/initial_state.nc') self.add_input_file( - filename=f'{mesh_name}_output.nc', + filename=f'output_r{refinement_factor:02g}.nc', work_dir_target=f'{forward.path}/output.nc') self.add_output_file('filament.png') @@ -77,7 +78,11 @@ def run(self): Run this step of the test case """ plt.switch_backend('Agg') - resolutions = self.resolutions + resolutions = list() + for refinement_factor in self.refinement_factors: + resolution = get_resolution_for_task( + self.config, refinement_factor, self.refinement) + resolutions.append(resolution) config = self.config section = config[self.case_name] eval_time = section.getfloat('filament_evaluation_time') @@ -89,9 +94,9 @@ def run(self): filament_norm = np.zeros((len(resolutions), num_tau)) use_mplstyle() fig, ax = plt.subplots() - for i, resolution in enumerate(resolutions): - mesh_name = resolution_to_subdir(resolution) - ds = xr.open_dataset(f'{mesh_name}_output.nc') + for i, refinement_factor in enumerate(self.refinement_factors): + mesh_name = resolution_to_subdir(resolutions[i]) + ds = xr.open_dataset(f'output_r{refinement_factor:02g}.nc') tidx = time_index_from_xtime(ds.xtime.values, eval_time * s_per_day) tracer = ds[variable_name] diff --git a/polaris/ocean/tasks/sphere_transport/forward.py b/polaris/ocean/tasks/sphere_transport/forward.py index c795993af..6ce2195b8 100644 --- a/polaris/ocean/tasks/sphere_transport/forward.py +++ b/polaris/ocean/tasks/sphere_transport/forward.py @@ -7,8 +7,8 @@ class Forward(SphericalConvergenceForward): transport test case """ - def __init__(self, component, name, subdir, resolution, base_mesh, init, - case_name): + def __init__(self, component, name, subdir, base_mesh, init, + case_name, refinement_factor, refinement): """ Create a new step @@ -23,9 +23,6 @@ def __init__(self, component, name, subdir, resolution, base_mesh, init, subdir : str The subdirectory for the step - resolution : float - The resolution of the (uniform) mesh in km - base_mesh : polaris.Step The base mesh step @@ -34,6 +31,13 @@ def __init__(self, component, name, subdir, resolution, base_mesh, init, case_name: str The name of the test case + + refinement_factor : float + The factor by which to scale space, time or both + + refinement : str, optional + Refinement type. One of 'space', 'time' or 'both' indicating both + space and time """ package = 'polaris.ocean.tasks.sphere_transport' flow_id = {'rotation_2d': 1, @@ -47,10 +51,12 @@ def __init__(self, component, name, subdir, resolution, base_mesh, init, } validate_vars = ['normalVelocity', 'tracer1', 'tracer2', 'tracer3'] super().__init__(component=component, name=name, subdir=subdir, - resolution=resolution, mesh=base_mesh, + mesh=base_mesh, init=init, package=package, yaml_filename='forward.yaml', output_filename='output.nc', validate_vars=validate_vars, options=namelist_options, - graph_target=f'{base_mesh.path}/graph.info') + graph_target=f'{base_mesh.path}/graph.info', + refinement_factor=refinement_factor, + refinement=refinement) diff --git a/polaris/ocean/tasks/sphere_transport/forward.yaml b/polaris/ocean/tasks/sphere_transport/forward.yaml index e8811b654..be8a53ec0 100644 --- a/polaris/ocean/tasks/sphere_transport/forward.yaml +++ b/polaris/ocean/tasks/sphere_transport/forward.yaml @@ -49,6 +49,8 @@ mpas-ocean: - tracers - mesh - xtime + - velocityZonal + - velocityMeridional - normalVelocity - layerThickness - refZMid diff --git a/polaris/ocean/tasks/sphere_transport/mixing_analysis.py b/polaris/ocean/tasks/sphere_transport/mixing_analysis.py index 750518a7f..88ee06702 100644 --- a/polaris/ocean/tasks/sphere_transport/mixing_analysis.py +++ b/polaris/ocean/tasks/sphere_transport/mixing_analysis.py @@ -7,6 +7,7 @@ from polaris import Step from polaris.mpas import time_index_from_xtime +from polaris.ocean.convergence import get_resolution_for_task from polaris.ocean.resolution import resolution_to_subdir from polaris.viz import use_mplstyle @@ -27,8 +28,8 @@ class MixingAnalysis(Step): case_name : str The name of the test case """ - def __init__(self, component, resolutions, icosahedral, subdir, - case_name, dependencies): + def __init__(self, component, refinement_factors, icosahedral, subdir, + case_name, dependencies, refinement='both'): """ Create the step @@ -55,22 +56,22 @@ def __init__(self, component, resolutions, icosahedral, subdir, """ super().__init__(component=component, name='mixing_analysis', subdir=subdir) - self.resolutions = resolutions + self.refinement_factors = refinement_factors + self.refinement = refinement self.case_name = case_name - for resolution in resolutions: - mesh_name = resolution_to_subdir(resolution) - base_mesh = dependencies['mesh'][resolution] - init = dependencies['init'][resolution] - forward = dependencies['forward'][resolution] + for refinement_factor in refinement_factors: + base_mesh = dependencies['mesh'][refinement_factor] + init = dependencies['init'][refinement_factor] + forward = dependencies['forward'][refinement_factor] self.add_input_file( - filename=f'{mesh_name}_mesh.nc', + filename=f'mesh_r{refinement_factor:02g}.nc', work_dir_target=f'{base_mesh.path}/base_mesh.nc') self.add_input_file( - filename=f'{mesh_name}_init.nc', + filename=f'init_r{refinement_factor:02g}.nc', work_dir_target=f'{init.path}/initial_state.nc') self.add_input_file( - filename=f'{mesh_name}_output.nc', + filename=f'output_r{refinement_factor:02g}.nc', work_dir_target=f'{forward.path}/output.nc') self.add_output_file('triplots.png') @@ -79,7 +80,11 @@ def run(self): Run this step of the test case """ plt.switch_backend('Agg') - resolutions = self.resolutions + resolutions = list() + for refinement_factor in self.refinement_factors: + resolution = get_resolution_for_task( + self.config, refinement_factor, self.refinement) + resolutions.append(resolution) config = self.config section = config[self.case_name] eval_time = section.getfloat('mixing_evaluation_time') @@ -89,12 +94,12 @@ def run(self): use_mplstyle() fig, axes = plt.subplots(nrows=nrows, ncols=2, sharex=True, sharey=True, figsize=(5.5, 7)) - for i, resolution in enumerate(resolutions): + for i, refinement_factor in enumerate(self.refinement_factors): ax = axes[int(i / 2), i % 2] _init_triplot_axes(ax) - mesh_name = resolution_to_subdir(resolution) + mesh_name = resolution_to_subdir(resolutions[i]) ax.set(title=mesh_name) - ds = xr.open_dataset(f'{mesh_name}_output.nc') + ds = xr.open_dataset(f'output_r{refinement_factor:02g}.nc') if i % 2 == 0: ax.set_ylabel("tracer3") if int(i / 2) == nrows - 1: diff --git a/polaris/ocean/tasks/sphere_transport/resources/flow_types.py b/polaris/ocean/tasks/sphere_transport/resources/flow_types.py index 78a9adf95..3645d2d27 100644 --- a/polaris/ocean/tasks/sphere_transport/resources/flow_types.py +++ b/polaris/ocean/tasks/sphere_transport/resources/flow_types.py @@ -43,7 +43,7 @@ def flow_nondivergent(t, lon, lat, u_0, tau): def flow_divergent(t, lon, lat, u_0, tau): """ - Compute a nondivergent velocity field + Compute a divergent velocity field Parameters ---------- @@ -81,7 +81,7 @@ def flow_divergent(t, lon, lat, u_0, tau): def flow_rotation(lon, lat, omega, tau, sphere_radius): """ - Compute a nondivergent velocity field + Compute a rotational velocity field Parameters ----------