Skip to content

Commit

Permalink
Use shared convergence for manufactured solution
Browse files Browse the repository at this point in the history
  • Loading branch information
cbegeman committed Oct 6, 2023
1 parent 870ee9f commit 1d766ca
Show file tree
Hide file tree
Showing 6 changed files with 120 additions and 119 deletions.
30 changes: 23 additions & 7 deletions polaris/ocean/tasks/manufactured_solution/__init__.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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')
101 changes: 45 additions & 56 deletions polaris/ocean/tasks/manufactured_solution/analysis.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand All @@ -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)
60 changes: 11 additions & 49 deletions polaris/ocean/tasks/manufactured_solution/forward.py
Original file line number Diff line number Diff line change
@@ -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.
Expand All @@ -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
Expand All @@ -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):
"""
Expand Down Expand Up @@ -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,
Expand Down
11 changes: 6 additions & 5 deletions polaris/ocean/tasks/manufactured_solution/forward.yaml
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand Down
5 changes: 5 additions & 0 deletions polaris/ocean/tasks/manufactured_solution/init.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand Down Expand Up @@ -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):
"""
Expand All @@ -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')
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand All @@ -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}

0 comments on commit 1d766ca

Please sign in to comment.