diff --git a/docs/developers_guide/ocean/api.md b/docs/developers_guide/ocean/api.md index 812efda9b..185998b97 100644 --- a/docs/developers_guide/ocean/api.md +++ b/docs/developers_guide/ocean/api.md @@ -104,8 +104,7 @@ forward.Forward.dynamic_model_config analysis.Analysis - analysis.Analysis.run - analysis.Analysis.rmse + analysis.Analysis.exact_solution viz.VizMap viz.VizMap.run @@ -178,6 +177,12 @@ :toctree: generated/ SphericalConvergenceForward + SphericalConvergenceAnalysis + SphericalConvergenceAnalysis.compute_error + SphericalConvergenceAnalysis.convergence_parameters + SphericalConvergenceAnalysis.exact_solution + SphericalConvergenceAnalysis.run + SphericalConvergenceAnalysis.plot_convergence ``` ### Ocean Model diff --git a/docs/developers_guide/ocean/framework.md b/docs/developers_guide/ocean/framework.md index 44d33b4d6..f2aea1908 100644 --- a/docs/developers_guide/ocean/framework.md +++ b/docs/developers_guide/ocean/framework.md @@ -159,6 +159,12 @@ qu_resolutions = 60, 90, 120, 150, 180, 210, 240 # Evaluation time for convergence analysis (in days) convergence_eval_time = 1.0 +# Convergence threshold below which a test fails +convergence_thresh = 1.0 + +# Type of error to compute +error_type = l2 + # config options for spherical convergence forward steps [spherical_convergence_forward] @@ -183,14 +189,19 @@ run_duration = ${spherical_convergence:convergence_eval_time} output_interval = ${run_duration} ``` The first 2 are the default resolutions for icosahedral and quasi-uniform -base meshes, respectively. The `time_integrator` will typically be overridden -by the specific convergence task's config options, and indicates which time -integrator to use for the forward run. Depending on the time integrator, -either `rk4_dt_per_km` or `split_dt_per_km` will be used to determine an -appropriate time step for each mesh resolution (proportional to the cell size). -For split time integrators, `btr_dt_per_km` will be used to compute the -barotropic time step in a similar way. The `run_duration` and -`output_interval` are typically the same, and they are given in days. +base meshes, respectively. The `convergence_eval_time` will generally be +modified by each test case. The `convergence_thresh` will also be modified by +each test case, and will depend on the numerical methods being tested. The +`error_type` is the L2 norm by default. + +`time_integrator` will typically be overridden by the specific convergence +task's config options, and indicates which time integrator to use for the +forward run. Depending on the time integrator, either `rk4_dt_per_km` or +`split_dt_per_km` will be used to determine an appropriate time step for each +mesh resolution (proportional to the cell size). For split time integrators, +`btr_dt_per_km` will be used to compute the barotropic time step in a similar +way. The `run_duration` and `output_interval` are typically the same, and +they are given in days. Each convergence test can override these defaults with its own defaults by defining them in its own config file. Convergence tests should bring in this @@ -311,6 +322,57 @@ omega: on the associated config options (first at setup and again at runtime in case the config options have changed). +In addition, the {py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceAnalysis` +step can serve as a parent class for analysis steps in convergence tests. This +parent class computes the error norm for the output from each resolution's +forward step. It also produces the convergence plot. + +This is an example of how a task's analysis step can descend from the parent +class: + +```python +class Analysis(SphericalConvergenceAnalysis): + """ + A step for analyzing the output from the cosine bell test case + """ + def __init__(self, component, resolutions, icosahedral, subdir, + dependencies): + """ + Create the step + + Parameters + ---------- + component : polaris.Component + The component the step belongs to + + resolutions : list of float + The resolutions of the meshes that have been run + + icosahedral : bool + Whether to use icosahedral, as opposed to less regular, JIGSAW + meshes + + subdir : str + The subdirectory that the step resides in + + dependencies : dict of dict of polaris.Steps + The dependencies of this step + """ + convergence_vars = [{'name': 'tracer1', + 'title': 'tracer1', + 'zidx': 0}] + super().__init__(component=component, subdir=subdir, + resolutions=resolutions, + icosahedral=icosahedral, + dependencies=dependencies, + convergence_vars=convergence_vars) +``` + +Many tasks will also need to override the `exact_solution` method given in +{py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceAnalysis`. +If not overridden, the analysis step will compute the difference of the output +from the initial state. + (dev-ocean-framework-vertical)= ## Vertical coordinate diff --git a/docs/developers_guide/ocean/tasks/cosine_bell.md b/docs/developers_guide/ocean/tasks/cosine_bell.md index 180012744..e825e39f4 100644 --- a/docs/developers_guide/ocean/tasks/cosine_bell.md +++ b/docs/developers_guide/ocean/tasks/cosine_bell.md @@ -44,8 +44,11 @@ section. Other model config options are taken from `forward.yaml`. ### analysis The class {py:class}`polaris.ocean.tasks.cosine_bell.analysis.Analysis` -defines a step for computing the RMSE (root-mean-squared error) for the results -at each resolution and plotting them in `convergence.png`. +descends from +{py:class}`polaris.ocean.convergence.spherical.SphericalConvergenceAnalysis`, +and defines a step for computing the error norm (L2) for the results +at each resolution, saving them in `convergence_tracer1.csv` and plotting them +in `convergence_tracer1.png`. ### viz diff --git a/docs/users_guide/ocean/tasks/cosine_bell.md b/docs/users_guide/ocean/tasks/cosine_bell.md index fd1c44b11..b344fae14 100644 --- a/docs/users_guide/ocean/tasks/cosine_bell.md +++ b/docs/users_guide/ocean/tasks/cosine_bell.md @@ -176,12 +176,6 @@ The convergence_eval_time, run duration and output interval are the period for advection to make a full rotation around the globe, 24 days: ```cfg -# config options for spherical convergence tests -[spherical_convergence] - -# Evaluation time for convergence analysis (in days) -convergence_eval_time = ${cosine_bell:vel_pd} - # config options for spherical convergence tests [spherical_convergence_forward] @@ -192,7 +186,7 @@ run_duration = ${cosine_bell:vel_pd} output_interval = ${cosine_bell:vel_pd} ``` -Her, `${cosine_bell:vel_pd}` means that the same value is used as in the +Here, `${cosine_bell:vel_pd}` means that the same value is used as in the option `vel_pd` in section `[cosine_bell]`, see below. ## config options @@ -224,17 +218,8 @@ psi0 = 1.0 # time (days) for bell to transit equator once vel_pd = 24.0 -# convergence threshold below which the test fails for QU meshes -qu_conv_thresh = 1.8 - -# Convergence rate above which a warning is issued for QU meshes -qu_conv_max = 2.2 - -# convergence threshold below which the test fails for icosahedral meshes -icos_conv_thresh = 1.8 - -# Convergence rate above which a warning is issued for icosahedral meshes -icos_conv_max = 2.2 +# convergence threshold below which the test fails +convergence_thresh = 1.8 # options for visualization for the cosine bell convergence test case @@ -270,6 +255,21 @@ when the convergence rates are not within the expected range. The options in the `cosine_bell_viz` section are used in visualizing the initial and final states on a lon-lat grid for `cosine_bell/with_viz` tasks. +By default, the convergence analysis step analyzes convergence after the +cosine bell has circulated the globe once. It also computes the L2 norm. Both +of these config options can be changed here: + +```cfg +# config options for spherical convergence tests +[spherical_convergence] + +# Evaluation time for convergence analysis (in days) +convergence_eval_time = ${cosine_bell:vel_pd} + +# Type of error to compute +error_type = l2 +``` + ## cores The target and minimum number of cores are determined by `goal_cells_per_core` diff --git a/docs/users_guide/ocean/tasks/images/cosine_bell_convergence.png b/docs/users_guide/ocean/tasks/images/cosine_bell_convergence.png index 0cfeae17d..f2c26b969 100644 Binary files a/docs/users_guide/ocean/tasks/images/cosine_bell_convergence.png and b/docs/users_guide/ocean/tasks/images/cosine_bell_convergence.png differ diff --git a/polaris/ocean/convergence/spherical/__init__.py b/polaris/ocean/convergence/spherical/__init__.py index 702a1352f..418252ded 100644 --- a/polaris/ocean/convergence/spherical/__init__.py +++ b/polaris/ocean/convergence/spherical/__init__.py @@ -1,3 +1,6 @@ +from polaris.ocean.convergence.spherical.analysis import ( + SphericalConvergenceAnalysis, +) from polaris.ocean.convergence.spherical.forward import ( SphericalConvergenceForward, ) diff --git a/polaris/ocean/convergence/spherical/analysis.py b/polaris/ocean/convergence/spherical/analysis.py new file mode 100644 index 000000000..2ccc7791f --- /dev/null +++ b/polaris/ocean/convergence/spherical/analysis.py @@ -0,0 +1,403 @@ +import datetime +import os + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import xarray as xr + +from polaris import Step +from polaris.ocean.resolution import resolution_to_subdir + + +class SphericalConvergenceAnalysis(Step): + """ + A step for analyzing the output from the geostrophic convergence test + + Attributes + ---------- + resolutions : list of float + The resolutions of the meshes that have been run + + icosahedral : bool + Whether to use icosahedral, as opposed to less regular, JIGSAW + meshes + + dependencies_dict : dict of dict of polaris.Steps + The dependencies of this step must be given as separate keys in the + dict: + + mesh : dict of polaris.Steps + Keys of the dict correspond to `resolutions` + 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` + 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` + Values of the dict are polaris.Steps, which must have the + attribute `path`, the path to `forward.nc` of that + resolution + + convergence_vars: list of dict + The attributes for each variable for which to analyze the convergence + rate. Each dict must contain the following keys: + + name : str + The name of the variable as given in the output netcdf file + + title : str + The name of the variable to use in the plot title + + zidx : int + The z-index to use for variables that have an nVertLevels + dimension, which should be None for variables that don't + """ + def __init__(self, component, resolutions, icosahedral, subdir, + dependencies, convergence_vars): + """ + Create the step + + Parameters + ---------- + component : polaris.Component + The component the step belongs to + + resolutions : list of float + The resolutions of the meshes that have been run + + icosahedral : bool + Whether to use icosahedral, as opposed to less regular, JIGSAW + meshes + + subdir : str + The subdirectory that the step resides in + + dependencies : dict of dict of polaris.Steps + The dependencies of this step must be given as separate keys in the + dict: + + mesh : dict of polaris.Steps + Keys of the dict correspond to `resolutions` + 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` + 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` + Values of the dict are polaris.Steps, which must have the + attribute `path`, the path to `forward.nc` of that + resolution + + convergence_vars: list of dict + The convergence attributes for each variable. Each dict must + contain the following keys: + + name : str + The name of the variable as given in the output netcdf file + + title : str + The name of the variable to use in the plot title + + zidx : int + The z-index to use for variables that have an nVertLevels + dimension, which should be None for variables that don't + """ + super().__init__(component=component, name='analysis', subdir=subdir) + self.resolutions = resolutions + self.icosahedral = icosahedral + self.dependencies_dict = dependencies + self.convergence_vars = convergence_vars + + for var in convergence_vars: + self.add_output_file(f'convergence_{var["name"]}.png') + + def setup(self): + """ + Add input files based on resolutions, which may have been changed by + user config options + """ + 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] + self.add_input_file( + filename=f'{mesh_name}_mesh.nc', + work_dir_target=f'{base_mesh.path}/base_mesh.nc') + self.add_input_file( + filename=f'{mesh_name}_init.nc', + work_dir_target=f'{init.path}/initial_state.nc') + self.add_input_file( + filename=f'{mesh_name}_output.nc', + work_dir_target=f'{forward.path}/output.nc') + + def run(self): + """ + Run this step of the test case + """ + plt.switch_backend('Agg') + convergence_vars = self.convergence_vars + for var in convergence_vars: + self.plot_convergence( + variable_name=var["name"], + title=var["title"], + zidx=var["zidx"]) + + def plot_convergence(self, variable_name, title, zidx): + """ + Compute the error norm for each resolution and produce a convergence + plot + + Parameters + ---------- + variable_name : str + The name of the variable as given in the output netcdf file + + title : str + The name of the variable to use in the plot title + + zidx : int + The z-index to use for variables that have an nVertLevels + dimension, which should be None for variables that don't + """ + resolutions = self.resolutions + 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) + error_res = self.compute_error( + mesh_name=mesh_name, + variable_name=variable_name, + zidx=zidx, + error_type=error_type) + error.append(error_res) + + res_array = np.array(resolutions) + 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]) + 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) + convergence = poly[0] + conv_round = np.round(convergence, 3) + + fit = res_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 + + fig = plt.figure() + + error_dict = {'l2': 'L2 norm', 'inf': 'L-infinity norm'} + error_title = error_dict[error_type] + + ax = fig.add_subplot(111) + ax.loglog(resolutions, order1, '--k', label='first order', + alpha=0.3) + ax.loglog(resolutions, order2, 'k', label='second order', + alpha=0.3) + ax.loglog(res_array, fit, 'k', + label=f'linear fit (order={conv_round})') + ax.loglog(res_array, error_array, 'o', label='numerical') + + if self.baseline_dir is not None: + baseline_filename = os.path.join(self.baseline_dir, filename) + if os.path.exists(baseline_filename): + data = pd.read_csv(baseline_filename) + if error_type not in data.keys(): + raise ValueError( + f'{error_type} not available for baseline') + else: + res_array = data.resolution.to_numpy() + error_array = data[error_type].to_numpy() + poly = np.polyfit( + np.log10(res_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', + label=f'linear fit, baseline ' + f'(order={conv_round})') + ax.set_xlabel('resolution (km)') + ax.set_ylabel(f'{error_title}') + ax.set_title(f'Error Convergence of {title}') + ax.legend(loc='lower left') + ax.invert_xaxis() + fig.savefig(f'convergence_{variable_name}.png', + bbox_inches='tight', pad_inches=0.1) + plt.close() + + logger.info(f'Order of convergence for {title}: ' + f'{conv_round}') + + if convergence < conv_thresh: + logger.error(f'Error: order of convergence for {title}\n' + f' {conv_round} < min tolerance ' + f'{conv_thresh}') + convergence_failed = True + + if convergence_failed: + raise ValueError('Convergence rate below minimum tolerance.') + + def compute_error(self, mesh_name, variable_name, zidx=None, + error_type='l2'): + """ + Compute the error for a given resolution + + Parameters + ---------- + mesh_name : str + The name of the mesh + + variable_name : str + The variable name of which to evaluate error with respect to the + exact solution + + zidx : int, optional + The z index to use to slice the field given by variable name + + error_type: {'l2', 'inf'}, optional + The type of error to compute + + Returns + ------- + error : float + The error of the variable given by variable_name + """ + ds_mesh = xr.open_dataset(f'{mesh_name}_mesh.nc') + ds_out = xr.open_dataset(f'{mesh_name}_output.nc') + config = self.config + section = config['spherical_convergence'] + eval_time = section.getfloat('convergence_eval_time') + s_per_day = 86400.0 + tidx = _time_index_from_xtime(ds_out.xtime.values, + eval_time * s_per_day) + ds_out = ds_out.isel(Time=tidx) + + if zidx is not None: + ds_out = ds_out.isel(nVertLevels=zidx) + field_exact = self.exact_solution(mesh_name, variable_name, + time=eval_time * s_per_day, + zidx=zidx) + field_mpas = ds_out[variable_name].values + diff = field_exact - field_mpas + + if error_type == 'l2': + area_cell = ds_mesh.areaCell.values + total_area = np.sum(area_cell) + den_l2 = np.sum(field_exact**2 * area_cell) / total_area + num_l2 = np.sum(diff**2 * area_cell) / total_area + error = np.sqrt(num_l2) / np.sqrt(den_l2) + elif error_type == 'inf': + error = np.amax(diff) / np.amax(np.abs(field_exact)) + else: + raise ValueError(f'Unsupported error type {error_type}') + + return error + + def exact_solution(self, mesh_name, 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 + + 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 + """ + + ds_init = xr.open_dataset(f'{mesh_name}_init.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 convergence_parameters(self, field_name=None): + """ + Get convergence parameters + + Parameters + ---------- + 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 + + Returns + ------- + conv_thresh : float + The minimum convergence rate + + error_type : {'l2', 'inf'}, str + The error norm to compute + """ + config = self.config + section = config['spherical_convergence'] + conv_thresh = section.getfloat('convergence_thresh') + error_type = section.get('error_type') + return conv_thresh, error_type + + +def _time_index_from_xtime(xtime, dt_target): + """ + Determine the time index at which to evaluate convergence + + Parameters + ---------- + xtime : list of str + Times in the dataset + + dt_target : float + Time in seconds at which to evaluate convergence + + Returns + ------- + tidx : int + Index in xtime that is closest to dt_target + """ + t0 = datetime.datetime.strptime(xtime[0].decode(), + '%Y-%m-%d_%H:%M:%S') + dt = np.zeros((len(xtime))) + for idx, xt in enumerate(xtime): + t = datetime.datetime.strptime(xt.decode(), + '%Y-%m-%d_%H:%M:%S') + dt[idx] = (t - t0).total_seconds() + return np.argmin(np.abs(np.subtract(dt, dt_target))) diff --git a/polaris/ocean/convergence/spherical/spherical.cfg b/polaris/ocean/convergence/spherical/spherical.cfg index 5b2e00269..630250da4 100644 --- a/polaris/ocean/convergence/spherical/spherical.cfg +++ b/polaris/ocean/convergence/spherical/spherical.cfg @@ -10,6 +10,12 @@ qu_resolutions = 60, 90, 120, 150, 180, 210, 240 # Evaluation time for convergence analysis (in days) convergence_eval_time = 1.0 +# Convergence threshold below which a test fails +convergence_thresh = 1.0 + +# Type of error to compute +error_type = l2 + # config options for spherical convergence forward steps [spherical_convergence_forward] diff --git a/polaris/ocean/tasks/cosine_bell/analysis.py b/polaris/ocean/tasks/cosine_bell/analysis.py index e9f888bbf..e749e8f81 100644 --- a/polaris/ocean/tasks/cosine_bell/analysis.py +++ b/polaris/ocean/tasks/cosine_bell/analysis.py @@ -1,28 +1,12 @@ -import warnings - -import matplotlib.pyplot as plt import numpy as np import xarray as xr -from polaris import Step -from polaris.ocean.resolution import resolution_to_subdir +from polaris.ocean.convergence.spherical import SphericalConvergenceAnalysis -class Analysis(Step): +class Analysis(SphericalConvergenceAnalysis): """ A step for analyzing the output from the cosine bell test case - - Attributes - ---------- - resolutions : list of float - The resolutions of the meshes that have been run - - icosahedral : bool - Whether to use icosahedral, as opposed to less regular, JIGSAW - meshes - - dependencies_dict : dict of dict of polaris.Steps - The dependencies of this step """ def __init__(self, component, resolutions, icosahedral, subdir, dependencies): @@ -47,143 +31,72 @@ def __init__(self, component, resolutions, icosahedral, subdir, dependencies : dict of dict of polaris.Steps The dependencies of this step """ - super().__init__(component=component, name='analysis', subdir=subdir) - self.resolutions = resolutions - self.icosahedral = icosahedral - self.dependencies_dict = dependencies - - self.add_output_file('convergence.png') - - def setup(self): - """ - Add input files based on resolutions, which may have been changed by - user config options - """ - 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] - self.add_input_file( - filename=f'{mesh_name}_mesh.nc', - work_dir_target=f'{base_mesh.path}/base_mesh.nc') - self.add_input_file( - filename=f'{mesh_name}_init.nc', - work_dir_target=f'{init.path}/initial_state.nc') - self.add_input_file( - filename=f'{mesh_name}_output.nc', - work_dir_target=f'{forward.path}/output.nc') - - def run(self): - """ - Run this step of the test case + convergence_vars = [{'name': 'tracer1', + 'title': 'tracer1', + 'zidx': 0}] + super().__init__(component=component, subdir=subdir, + resolutions=resolutions, + icosahedral=icosahedral, + dependencies=dependencies, + convergence_vars=convergence_vars) + + def exact_solution(self, mesh_name, field_name, time, zidx=None): """ - plt.switch_backend('Agg') - resolutions = self.resolutions - xdata = list() - ydata = list() - for resolution in resolutions: - mesh_name = resolution_to_subdir(resolution) - rmseValue, nCells = self.rmse(mesh_name) - xdata.append(nCells) - ydata.append(rmseValue) - xdata = np.asarray(xdata) - ydata = np.asarray(ydata) - - p = np.polyfit(np.log10(xdata), np.log10(ydata), 1) - conv = abs(p[0]) * 2.0 - - yfit = xdata ** p[0] * 10 ** p[1] - - plt.loglog(xdata, yfit, 'k') - plt.loglog(xdata, ydata, 'or') - plt.annotate(f'Order of Convergence = {np.round(conv, 3)}', - xycoords='axes fraction', xy=(0.3, 0.95), fontsize=14) - plt.xlabel('Number of Grid Cells', fontsize=14) - plt.ylabel('L2 Norm', fontsize=14) - plt.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1) - - section = self.config['cosine_bell'] - if self.icosahedral: - conv_thresh = section.getfloat('icos_conv_thresh') - conv_max = section.getfloat('icos_conv_max') - else: - conv_thresh = section.getfloat('qu_conv_thresh') - conv_max = section.getfloat('qu_conv_max') - - 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}') - - def rmse(self, mesh_name): - """ - Compute the RMSE for a given resolution + Get the exact solution Parameters ---------- mesh_name : str - The name of the mesh + 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 ------- - rmseValue : float - The root-mean-squared error - - nCells : int - The number of cells in the mesh + solution: xarray.DataArray + The exact solution with dimension nCells """ + if field_name != 'tracer1': + print(f'Variable {field_name} not available as an analytic ' + 'solution for the cosine_bell test case') + config = self.config latCent = config.getfloat('cosine_bell', 'lat_center') lonCent = config.getfloat('cosine_bell', 'lon_center') radius = config.getfloat('cosine_bell', 'radius') psi0 = config.getfloat('cosine_bell', 'psi0') - convergence_eval_time = config.getfloat('spherical_convergence', - 'convergence_eval_time') + vel_pd = config.getfloat('cosine_bell', 'vel_pd') ds_mesh = xr.open_dataset(f'{mesh_name}_mesh.nc') ds_init = xr.open_dataset(f'{mesh_name}_init.nc') - # find time since the beginning of run - ds = xr.open_dataset(f'{mesh_name}_output.nc') - for j in range(len(ds.xtime)): - tt = str(ds.xtime[j].values) - tt.rfind('_') - DY = float(tt[10:12]) - 1 - if DY == convergence_eval_time: - sliceTime = j - break - HR = float(tt[13:15]) - MN = float(tt[16:18]) - t = 86400.0 * DY + HR * 3600. + MN # find new location of blob center # center is based on equatorial velocity R = ds_mesh.sphere_radius - distTrav = 2.0 * 3.14159265 * R / (86400.0 * convergence_eval_time) * t # distance in radians is - distRad = distTrav / R + distRad = 2.0 * np.pi * time / (86400.0 * vel_pd) newLon = lonCent + distRad if newLon > 2.0 * np.pi: newLon -= 2.0 * np.pi - - # construct analytic tracer tracer = np.zeros_like(ds_init.tracer1[0, :, 0].values) - latC = ds_mesh.latCell.values - lonC = ds_mesh.lonCell.values + latC = ds_init.latCell.values + lonC = ds_init.lonCell.values + # TODO replace this with great circle distance temp = R * np.arccos(np.sin(latCent) * np.sin(latC) + np.cos(latCent) * np.cos(latC) * np.cos( lonC - newLon)) mask = temp < radius tracer[mask] = (psi0 / 2.0 * (1.0 + np.cos(3.1415926 * temp[mask] / radius))) - - # oad forward mode data - tracerF = ds.tracer1[sliceTime, :, 0].values - rmseValue = np.sqrt(np.mean((tracerF - tracer)**2)) - - return rmseValue, ds_init.dims['nCells'] + return xr.DataArray(data=tracer, dims=('nCells',)) diff --git a/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg b/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg index 5909e8fa6..c2b8b5ff7 100644 --- a/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg +++ b/polaris/ocean/tasks/cosine_bell/cosine_bell.cfg @@ -26,6 +26,12 @@ min_pc_fraction = 0.1 # Evaluation time for convergence analysis (in days) convergence_eval_time = ${cosine_bell:vel_pd} +# Convergence threshold below which a test fails +convergence_thresh = ${cosine_bell:convergence_thresh} + +# Type of error to compute +error_type = l2 + # config options for spherical convergence tests [spherical_convergence_forward] @@ -67,17 +73,8 @@ psi0 = 1.0 # time (days) for bell to transit equator once vel_pd = 24.0 -# convergence threshold below which the test fails for QU meshes -qu_conv_thresh = 1.8 - -# Convergence rate above which a warning is issued for QU meshes -qu_conv_max = 2.2 - -# convergence threshold below which the test fails for icosahedral meshes -icos_conv_thresh = 1.8 - -# Convergence rate above which a warning is issued for icosahedral meshes -icos_conv_max = 2.2 +# convergence threshold below which the test fails +convergence_thresh = 1.8 # options for visualization for the cosine bell convergence test case