Source code for compass.ocean.tests.merry_go_round.convergence_test.analysis

import math
import numpy
import xarray
import matplotlib
import matplotlib.pyplot as plt
import cmocean

from compass.step import Step

[docs] class Analysis(Step): """ A step for plotting the convergence of the solution with resolution and time step in the merry-go-round test group Attributes ---------- resolution : str The resolution of the test case """
[docs] def __init__(self, test_case, resolutions, name='analysis'): """ Create the step Parameters ---------- test_case : compass.TestCase The test case this step belongs to resolutions : list of str The resolutions of the test case name: str The name of the step """ super().__init__(test_case=test_case, name=name) self.resolutions = resolutions self.add_output_file(filename='convergence_plot.png')
[docs] def run(self): """ Run this step of the test case """ _plot(self.outputs[0], self.resolutions)
def _plot(filename, resolutions): """ Plot section of the merry-go-round TODO Parameters ---------- filename : str The output file name resolutions : list of str The resolutions of the test case """ plt.switch_backend('Agg') fig = plt.gcf() dt = [3, 6, 12] order2 = [0.01, 0.04, 0.16] fields = ['tracer1'] L2 = numpy.zeros((len(resolutions))) for k, field in enumerate(fields): for i, resolution in enumerate(resolutions): ds = xarray.open_dataset(f'../forward_{resolution}/') areaCell = ds.areaCell.values final_field = ds[field].isel(Time=1, nVertLevels=0).values initial_field = ds[field].isel(Time=0, nVertLevels=0).values diff = abs(final_field - initial_field) multDen = (initial_field**2)*areaCell multNum = (diff**2)*areaCell denL2 = numpy.sum(multDen)/numpy.sum(areaCell) numL2 = numpy.sum(multNum)/numpy.sum(areaCell) L2[i] = numpy.sqrt(numL2)/numpy.sqrt(denL2) print(f'Order of convergence from dt 6 min to 3 min: ', f'{math.log2(L2[0]/L2[1])}') print(f'Order of convergence from dt 12 min to 6 min: ', f'{math.log2(L2[1]/L2[2])}') plt.loglog(dt, L2[:], '-x', label=f'Simulated {field}') plt.loglog(dt, order2, 'k', label='Order 2 convergence') plt.title('Convergence to the exact solution') plt.ylabel('l_2 error norm') plt.legend() plt.grid() plt.xticks(dt, dt) plt.xlabel('time step (min)') plt.savefig(filename)