-
Notifications
You must be signed in to change notification settings - Fork 37
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #703 from cbegeman/enhance-drying-slope-with-conve…
…rgence Add drying_slope convergence tests: This PR adds a convergence test case for both wetting-and-drying methods, `standard` and `ramp`. The RMSE is computed at all locations and times for which an analytic solution is available (the analytic solution was provided prior to this PR). It produces a convergence plot and compares both methods when both have been run. The analytic solution is only comparable to the `default` cases, which use Rayleigh damping.
- Loading branch information
Showing
20 changed files
with
591 additions
and
374 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,168 @@ | ||
| import os | ||
|
|
||
| import matplotlib.pyplot as plt | ||
| import numpy as np | ||
| import pandas as pd | ||
| import xarray as xr | ||
| from scipy.interpolate import interp1d | ||
|
|
||
| from compass.step import Step | ||
|
|
||
|
|
||
| class Analysis(Step): | ||
| """ | ||
| A step for analyzing the convergence of drying slope results and producing | ||
| a convergence plot. | ||
| Attributes | ||
| ---------- | ||
| damping_coeff : float | ||
| The Rayleigh damping coefficient used for the forward runs | ||
| resolutions : float | ||
| The resolution of the test case | ||
| times : list of float | ||
| The times at which to compare to the analytical solution | ||
| """ | ||
| def __init__(self, test_case, resolutions, damping_coeff): | ||
| """ | ||
| Create a forward step | ||
| Parameters | ||
| ---------- | ||
| test_case : compass.TestCase | ||
| The test case this step belongs to | ||
| resolutions : list of floats | ||
| The resolution of the test case | ||
| damping_coeff: float | ||
| the value of the rayleigh damping coefficient | ||
| """ | ||
|
|
||
| super().__init__(test_case=test_case, name='analysis') | ||
|
|
||
| self.damping_coeff = damping_coeff | ||
| self.resolutions = resolutions | ||
| self.times = ['0.05', '0.15', '0.25', '0.30', '0.40', '0.50'] | ||
|
|
||
| for resolution in resolutions: | ||
| if resolution < 1.: | ||
| res_name = f'{int(resolution*1e3)}m' | ||
| else: | ||
| res_name = f'{int(resolution)}km' | ||
| self.add_input_file(filename=f'output_{res_name}.nc', | ||
| target=f'../forward_{res_name}/output.nc') | ||
| for time in self.times: | ||
| filename = f'r{damping_coeff}d{time}-analytical'\ | ||
| '.csv' | ||
| self.add_input_file(filename=filename, target=filename, | ||
| database='drying_slope') | ||
|
|
||
| self.add_output_file(filename='convergence.png') | ||
|
|
||
| def run(self): | ||
| """ | ||
| Run this step of the test case | ||
| """ | ||
| self._plot_convergence() | ||
|
|
||
| def _compute_rmse(self, ds, t): | ||
| """ | ||
| Get the time step | ||
| Parameters | ||
| ---------- | ||
| ds : xarray.Dataset | ||
| the MPAS dataset containing output from a forward step | ||
| t : float | ||
| the time to evaluate the RMSE at | ||
| Returns | ||
| ------- | ||
| rmse : float | ||
| the root-mean-squared-error of the MPAS output relative to the | ||
| analytical solution at time t | ||
| """ | ||
|
|
||
| x_exact, ssh_exact = self._exact_solution(t) | ||
| tidx = int((float(t) / 0.2 + 1e-16) * 24.0) | ||
| ds = ds.drop_vars(np.setdiff1d([j for j in ds.variables], | ||
| ['yCell', 'ssh'])) | ||
| ds = ds.isel(Time=tidx) | ||
| x_mpas = ds.yCell.values / 1000.0 | ||
| ssh_mpas = ds.ssh.values | ||
| # Interpolate mpas solution to the points at which we have an exact | ||
| # solution | ||
| idx_min = np.argwhere(x_exact - x_mpas[0] >= 0.).item(0) | ||
| idx_max = np.argwhere(x_exact - x_mpas[-1] <= 0.).item(-1) | ||
| f = interp1d(x_mpas, ssh_mpas) | ||
| ssh_mpas_interp = f(x_exact[idx_min:idx_max]) | ||
| rmse = np.sqrt(np.mean(np.square(ssh_mpas_interp - | ||
| ssh_exact[idx_min:idx_max]))) | ||
| return rmse | ||
|
|
||
| def _plot_convergence(self): | ||
| """ | ||
| Plot convergence curves | ||
| """ | ||
| comparisons = [] | ||
| cases = {'standard': '../../../standard/convergence/analysis', | ||
| 'ramp': '../../../ramp/convergence/analysis'} | ||
| for case in cases: | ||
| include = True | ||
| for resolution in self.resolutions: | ||
| if resolution < 1.: | ||
| res_name = f'{int(resolution*1e3)}m' | ||
| else: | ||
| res_name = f'{int(resolution)}km' | ||
| if not os.path.exists(f'{cases[case]}/output_{res_name}.nc'): | ||
| include = False | ||
| if include: | ||
| comparisons.append(case) | ||
|
|
||
| fig, ax = plt.subplots(nrows=1, ncols=1) | ||
|
|
||
| max_rmse = 0 | ||
| for k, comp in enumerate(comparisons): | ||
| rmse = np.zeros((len(self.resolutions), len(self.times))) | ||
| for i, resolution in enumerate(self.resolutions): | ||
| if resolution < 1.: | ||
| res_name = f'{int(resolution*1e3)}m' | ||
| else: | ||
| res_name = f'{int(resolution)}km' | ||
| ds = xr.open_dataset( | ||
| f'{cases[comp]}/output_{res_name}.nc') | ||
| for j, time in enumerate(self.times): | ||
| rmse[i, j] = self._compute_rmse(ds, time) | ||
| rmse_tav = np.mean(rmse, axis=1) | ||
| if np.max(rmse_tav) > max_rmse: | ||
| max_rmse = np.max(rmse_tav) | ||
| ax.loglog(self.resolutions, rmse_tav, | ||
| linestyle='-', marker='o', label=comp) | ||
|
|
||
| rmse_1st_order = np.zeros(len(self.resolutions)) | ||
| rmse_1st_order[0] = max_rmse | ||
| for i in range(len(self.resolutions) - 1): | ||
| rmse_1st_order[i + 1] = rmse_1st_order[i] / 2.0 | ||
|
|
||
| ax.loglog(self.resolutions, np.flip(rmse_1st_order), | ||
| linestyle='-', color='k', alpha=.25, label='1st order') | ||
|
|
||
| ax.set_xlabel('Cell size (km)') | ||
| ax.set_ylabel('RMSE (m)') | ||
|
|
||
| ax.legend(loc='lower right') | ||
| ax.set_title('SSH convergence') | ||
| fig.tight_layout() | ||
| fig.savefig('convergence.png') | ||
|
|
||
| def _exact_solution(self, time): | ||
| """ | ||
| Returns distance, ssh of the analytic solution | ||
| """ | ||
| datafile = f'./r{self.damping_coeff}d{time}-'\ | ||
| f'analytical.csv' | ||
| data = pd.read_csv(datafile, header=None) | ||
| return data[0], -data[1] |
126 changes: 126 additions & 0 deletions
126
compass/ocean/tests/drying_slope/convergence/__init__.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,126 @@ | ||
| from numpy import ceil | ||
|
|
||
| from compass.config import CompassConfigParser | ||
| from compass.ocean.tests.drying_slope.analysis import Analysis | ||
| from compass.ocean.tests.drying_slope.forward import Forward | ||
| from compass.ocean.tests.drying_slope.initial_state import InitialState | ||
| from compass.ocean.tests.drying_slope.viz import Viz | ||
| from compass.testcase import TestCase | ||
| from compass.validate import compare_variables | ||
|
|
||
|
|
||
| class Convergence(TestCase): | ||
| """ | ||
| The convergence drying_slope test case | ||
| Attributes | ||
| ---------- | ||
| resolution : float | ||
| The resolution of the test case in km | ||
| coord_type : str | ||
| The type of vertical coordinate (``sigma``, ``single_layer``, etc.) | ||
| damping_coeffs: list of float | ||
| The damping coefficients at which to evaluate convergence. Must be of | ||
| length 1. | ||
| """ | ||
|
|
||
| def __init__(self, test_group, coord_type, method): | ||
| """ | ||
| Create the test case | ||
| Parameters | ||
| ---------- | ||
| test_group : compass.ocean.tests.drying_slope.DryingSlope | ||
| The test group that this test case belongs to | ||
| coord_type : str | ||
| The type of vertical coordinate (``sigma``, ``single_layer``) | ||
| method: str | ||
| The wetting-and-drying method (``standard``, ``ramp``) | ||
| """ | ||
| name = 'convergence' | ||
|
|
||
| self.coord_type = coord_type | ||
| damping_coeffs = [0.01] | ||
| self.damping_coeffs = damping_coeffs | ||
| subdir = f'{coord_type}/{method}/{name}' | ||
| super().__init__(test_group=test_group, name=name, | ||
| subdir=subdir) | ||
| self.resolutions = None | ||
| # add the steps with default resolutions so they can be listed | ||
| config = CompassConfigParser() | ||
| config.add_from_package('compass.ocean.tests.drying_slope', | ||
| 'drying_slope.cfg') | ||
| self._setup_steps(config, subdir=subdir, method=method) | ||
|
|
||
| def _setup_steps(self, config, subdir, method): | ||
| """ setup steps given resolutions """ | ||
|
|
||
| default_resolutions = '0.25, 0.5, 1, 2' | ||
|
|
||
| # set the default values that a user may change before setup | ||
| config.set('drying_slope_convergence', 'resolutions', | ||
| default_resolutions, | ||
| comment='a list of resolutions (km) to test') | ||
|
|
||
| # get the resolutions back, perhaps with values set in the user's | ||
| # config file | ||
| resolutions = config.getlist('drying_slope_convergence', | ||
| 'resolutions', dtype=float) | ||
|
|
||
| if self.resolutions is not None and self.resolutions == resolutions: | ||
| return | ||
|
|
||
| # start fresh with no steps | ||
| self.steps = dict() | ||
| self.steps_to_run = list() | ||
|
|
||
| self.resolutions = resolutions | ||
| section = config['drying_slope'] | ||
| ntasks_baseline = section.getint('ntasks_baseline') | ||
| min_tasks = section.getint('min_tasks') | ||
|
|
||
| for resolution in self.resolutions: | ||
|
|
||
| if resolution < 1.: | ||
| res_name = f'{int(resolution*1e3)}m' | ||
| else: | ||
| res_name = f'{int(resolution)}km' | ||
| init_name = f'initial_state_{res_name}' | ||
| self.add_step(InitialState(test_case=self, | ||
| name=init_name, | ||
| resolution=resolution, | ||
| coord_type=self.coord_type)) | ||
| ntasks = max(min_tasks, | ||
| int(ceil(ntasks_baseline / resolution**2.))) | ||
| forward_step = Forward(test_case=self, resolution=resolution, | ||
| name=f'forward_{res_name}', | ||
| input_path=f'../{init_name}', | ||
| ntasks=ntasks, min_tasks=min_tasks, | ||
| openmp_threads=1, | ||
| damping_coeff=self.damping_coeffs[0], | ||
| coord_type=self.coord_type) | ||
| if method == 'ramp': | ||
| forward_step.add_namelist_options( | ||
| {'config_zero_drying_velocity_ramp': ".true."}) | ||
| self.add_step(forward_step) | ||
| self.add_step(Analysis(test_case=self, | ||
| resolutions=resolutions, | ||
| damping_coeff=self.damping_coeffs[0])) | ||
|
|
||
| def validate(self): | ||
| """ | ||
| Validate variables against a baseline | ||
| """ | ||
| super().validate() | ||
| variables = ['layerThickness', 'normalVelocity'] | ||
| for resolution in self.resolutions: | ||
| if resolution < 1.: | ||
| res_name = f'{int(resolution*1e3)}m' | ||
| else: | ||
| res_name = f'{int(resolution)}km' | ||
| compare_variables(test_case=self, variables=variables, | ||
| filename1=f'forward_{res_name}/output.nc') |
Oops, something went wrong.