diff --git a/polaris/ocean/tasks/manufactured_solution/__init__.py b/polaris/ocean/tasks/manufactured_solution/__init__.py index 056e89047..2c95b634e 100644 --- a/polaris/ocean/tasks/manufactured_solution/__init__.py +++ b/polaris/ocean/tasks/manufactured_solution/__init__.py @@ -1,4 +1,7 @@ -from polaris import Task +from typing import Dict + +from polaris import Step, 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 from polaris.ocean.tasks.manufactured_solution.init import Init @@ -39,19 +42,32 @@ def __init__(self, component): super().__init__(component=component, name=name, subdir=subdir) self.resolutions = [200., 100., 50., 25.] - for res in self.resolutions: - self.add_step(Init(component=component, resolution=res, - taskdir=self.subdir)) - self.add_step(Forward(component=component, resolution=res, - taskdir=self.subdir)) + analysis_dependencies: Dict[str, Dict[float, Step]] = ( + dict(mesh=dict(), init=dict(), forward=dict())) + for resolution in self.resolutions: + 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 self.add_step(Analysis(component=component, resolutions=self.resolutions, - taskdir=self.subdir)) + subdir=f'{self.subdir}/analysis', + dependencies=analysis_dependencies)) self.add_step(Viz(component=component, resolutions=self.resolutions, taskdir=self.subdir), run_by_default=False) + self.config.add_from_package('polaris.ocean.convergence', + 'convergence.cfg') self.config.add_from_package( 'polaris.ocean.tasks.manufactured_solution', 'manufactured_solution.cfg') diff --git a/polaris/ocean/tasks/manufactured_solution/analysis.py b/polaris/ocean/tasks/manufactured_solution/analysis.py index a1748049e..25a17e717 100644 --- a/polaris/ocean/tasks/manufactured_solution/analysis.py +++ b/polaris/ocean/tasks/manufactured_solution/analysis.py @@ -1,18 +1,12 @@ -import datetime -import warnings - -import cmocean # noqa: F401 -import numpy as np import xarray as xr -from polaris import Step -from polaris.ocean.resolution import resolution_to_subdir +from polaris.ocean.convergence import ConvergenceAnalysis from polaris.ocean.tasks.manufactured_solution.exact_solution import ( ExactSolution, ) -class Analysis(Step): +class Analysis(ConvergenceAnalysis): """ A step for analysing the output from the manufactured solution test case @@ -22,7 +16,7 @@ class Analysis(Step): resolutions : list of float The resolutions of the meshes that have been run """ - def __init__(self, component, resolutions, taskdir): + def __init__(self, component, resolutions, subdir, dependencies): """ Create the step @@ -34,54 +28,49 @@ def __init__(self, component, resolutions, taskdir): resolutions : list of float The resolutions of the meshes that have been run - taskdir : str - The subdirectory that the task belongs to - """ - super().__init__(component=component, name='analysis', indir=taskdir) - self.resolutions = resolutions - - for resolution in resolutions: - mesh_name = resolution_to_subdir(resolution) - self.add_input_file( - filename=f'init_{mesh_name}.nc', - target=f'../init/{mesh_name}/initial_state.nc') - self.add_input_file( - filename=f'output_{mesh_name}.nc', - target=f'../forward/{mesh_name}/output.nc') + subdir : str + The subdirectory that the step resides in - def run(self): + dependencies : dict of dict of polaris.Steps + The dependencies of this step """ - Run this step of the test case + convergence_vars = [{'name': 'ssh', + 'title': 'SSH', + 'zidx': None}] + super().__init__(component=component, subdir=subdir, + resolutions=resolutions, + dependencies=dependencies, + convergence_vars=convergence_vars) + + def exact_solution(self, mesh_name, field_name, time, zidx=None): """ - config = self.config - resolutions = self.resolutions - - section = config['manufactured_solution'] - conv_thresh = section.getfloat('conv_thresh') - conv_max = section.getfloat('conv_max') - - rmse = [] - for i, res in enumerate(resolutions): - mesh_name = f'{res:g}km' - init = xr.open_dataset(f'init_{mesh_name}.nc') - ds = xr.open_dataset(f'output_{mesh_name}.nc') - exact = ExactSolution(config, init) - - t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(), - '%Y-%m-%d_%H:%M:%S') - tf = datetime.datetime.strptime(ds.xtime.values[-1].decode(), - '%Y-%m-%d_%H:%M:%S') - t = (tf - t0).total_seconds() - ssh_model = ds.ssh.values[-1, :] - rmse.append(np.sqrt(np.mean((ssh_model - exact.ssh(t).values)**2))) + Get the exact solution - p = np.polyfit(np.log10(resolutions), np.log10(rmse), 1) - conv = p[0] - - if conv < conv_thresh: - raise ValueError(f'order of convergence ' - f' {conv} < min tolerence {conv_thresh}') - - if conv > conv_max: - warnings.warn(f'order of convergence ' - f'{conv} > max tolerence {conv_max}') + Parameters + ---------- + mesh_name : str + The mesh name which is the prefix for the initial condition file + + field_name : str + The name of the variable of which we evaluate convergence + For the default method, we use the same convergence rate for all + fields + + time : float + The time at which to evaluate the exact solution in seconds. + For the default method, we always use the initial state. + + zidx : int, optional + The z-index for the vertical level at which to evaluate the exact + solution + + Returns + ------- + solution : xarray.DataArray + The exact solution as derived from the initial condition + """ + init = xr.open_dataset(f'{mesh_name}_init.nc') + exact = ExactSolution(self.config, init) + if field_name != 'ssh': + raise ValueError(f'{field_name} is not currently supported') + return exact.ssh(time) diff --git a/polaris/ocean/tasks/manufactured_solution/forward.py b/polaris/ocean/tasks/manufactured_solution/forward.py index d8f56c92b..198d38c1b 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.py +++ b/polaris/ocean/tasks/manufactured_solution/forward.py @@ -1,16 +1,13 @@ -import time - import numpy as np from polaris.mesh.planar import compute_planar_hex_nx_ny -from polaris.ocean.model import OceanModelStep -from polaris.ocean.resolution import resolution_to_subdir +from polaris.ocean.convergence import ConvergenceForward from polaris.ocean.tasks.manufactured_solution.exact_solution import ( ExactSolution, ) -class Forward(OceanModelStep): +class Forward(ConvergenceForward): """ A step for performing forward ocean component runs as part of manufactured solution test cases. @@ -20,8 +17,7 @@ class Forward(OceanModelStep): resolution : float The resolution of the test case in km """ - def __init__(self, component, resolution, taskdir, - ntasks=None, min_tasks=None, openmp_threads=1): + def __init__(self, component, name, resolution, subdir, init): """ Create a new test case @@ -35,40 +31,14 @@ def __init__(self, component, resolution, taskdir, taskdir : str The subdirectory that the task belongs to - - ntasks : int, optional - the number of tasks the step would ideally use. If fewer tasks - are available on the system, the step will run on all available - tasks as long as this is not below ``min_tasks`` - - min_tasks : int, optional - the number of tasks the step requires. If the system has fewer - than this number of tasks, the step will fail - - openmp_threads : int, optional - the number of OpenMP threads the step will use """ - self.resolution = resolution - mesh_name = resolution_to_subdir(resolution) super().__init__(component=component, - name=f'forward_{mesh_name}', - subdir=f'{taskdir}/forward/{mesh_name}', - ntasks=ntasks, min_tasks=min_tasks, - openmp_threads=openmp_threads) - - self.add_input_file(filename='initial_state.nc', - target=f'../../init/{mesh_name}/initial_state.nc') - self.add_input_file(filename='graph.info', - target=f'../../init/{mesh_name}/culled_graph.info') - - self.add_output_file( - filename='output.nc', - validate_vars=['layerThickness', 'normalVelocity']) - - self.add_yaml_file('polaris.ocean.config', - 'single_layer.yaml') - self.add_yaml_file('polaris.ocean.tasks.manufactured_solution', - 'forward.yaml') + name=name, subdir=subdir, + resolution=resolution, base_mesh=init, init=init, + package='polaris.ocean.tasks.manufactured_solution', + yaml_filename='forward.yaml', + output_filename='output.nc', + validate_vars=['layerThickness', 'normalVelocity']) def compute_cell_count(self): """ @@ -100,16 +70,8 @@ def dynamic_model_config(self, at_setup): """ super().dynamic_model_config(at_setup=at_setup) - # dt is proportional to resolution - config = self.config - section = config['manufactured_solution'] - dt_per_km = section.getfloat('dt_per_km') - dt = dt_per_km * self.resolution - # https://stackoverflow.com/a/1384565/7728169 - dt_str = time.strftime('%H:%M:%S', time.gmtime(dt)) - exact_solution = ExactSolution(config) - options = {'config_dt': dt_str, - 'config_manufactured_solution_amplitude': + exact_solution = ExactSolution(self.config) + options = {'config_manufactured_solution_amplitude': exact_solution.eta0, 'config_manufactured_solution_wavelength_x': exact_solution.lambda_x, diff --git a/polaris/ocean/tasks/manufactured_solution/forward.yaml b/polaris/ocean/tasks/manufactured_solution/forward.yaml index 6b3502fce..aa40ed4b1 100644 --- a/polaris/ocean/tasks/manufactured_solution/forward.yaml +++ b/polaris/ocean/tasks/manufactured_solution/forward.yaml @@ -1,8 +1,9 @@ omega: time_management: - config_run_duration: 10:00:00 + config_run_duration: {{ run_duration }} time_integration: - config_time_integrator: RK4 + config_dt: {{ dt }} + config_time_integrator: {{ time_integrator }} bottom_drag: config_bottom_drag_mode: implicit config_implicit_bottom_drag_type: constant @@ -13,14 +14,14 @@ omega: config_disable_vel_hmix: true streams: mesh: - filename_template: initial_state.nc + filename_template: init.nc input: - filename_template: initial_state.nc + filename_template: init.nc restart: {} output: type: output filename_template: output.nc - output_interval: 0000_00:20:00 + output_interval: {{ output_interval }} clobber_mode: truncate contents: - xtime diff --git a/polaris/ocean/tasks/manufactured_solution/init.py b/polaris/ocean/tasks/manufactured_solution/init.py index 73c997d48..7942b10b9 100644 --- a/polaris/ocean/tasks/manufactured_solution/init.py +++ b/polaris/ocean/tasks/manufactured_solution/init.py @@ -5,6 +5,7 @@ from mpas_tools.planar_hex import make_planar_hex_mesh from polaris import Step +from polaris.io import symlink from polaris.mesh.planar import compute_planar_hex_nx_ny from polaris.ocean.resolution import resolution_to_subdir from polaris.ocean.tasks.manufactured_solution.exact_solution import ( @@ -43,6 +44,9 @@ def __init__(self, component, resolution, taskdir): name=f'init_{mesh_name}', subdir=f'{taskdir}/init/{mesh_name}') self.resolution = resolution + for filename in ['culled_mesh.nc', 'initial_state.nc', + 'culled_graph.info']: + self.add_output_file(filename=filename) def run(self): """ @@ -69,6 +73,7 @@ def run(self): ds_mesh = cull(ds_mesh, logger=logger) ds_mesh = convert(ds_mesh, graphInfoFileName='culled_graph.info', logger=logger) + symlink('culled_graph.info', 'graph.info') write_netcdf(ds_mesh, 'culled_mesh.nc') bottom_depth = config.getfloat('vertical_grid', 'bottom_depth') diff --git a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg index 14d42ca80..f241fcb8e 100644 --- a/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg +++ b/polaris/ocean/tasks/manufactured_solution/manufactured_solution.cfg @@ -30,8 +30,8 @@ dt_per_km = 3.0 # Convergence threshold below which the test fails conv_thresh = 1.8 -# Convergence rate above which a warning is issued -conv_max = 2.2 +# Run duration in days +run_duration = 0.4166667 [vertical_grid] @@ -52,3 +52,31 @@ partial_cell_type = None # The minimum fraction of a layer for partial cells min_pc_fraction = 0.1 + +# config options for spherical convergence tests +[convergence] + +# Evaluation time for convergence analysis (in days) +convergence_eval_time = ${manufactured_solution:run_duration} + +# Convergence threshold below which a test fails +convergence_thresh = ${manufactured_solution:conv_thresh} + +# Type of error to compute +error_type = l2 + + +# config options for spherical convergence tests +[convergence_forward] + +# time integrator: {'split_explicit', 'RK4'} +time_integrator = RK4 + +# RK4 time step per resolution (s/km), since dt is proportional to resolution +rk4_dt_per_km = ${manufactured_solution:dt_per_km} + +# Run duration in days +run_duration = ${manufactured_solution:run_duration} + +# Output interval in days +output_interval = ${manufactured_solution:run_duration}