Skip to content

Commit 4086e05

Browse files
matthewhoffmanmatthewhoffman
committed
Add thickness advection version of horizontal_advection test
This is identical to the horizontal_advection test of a passive tracer, but Gaussian bump is applied to thickness field and tracer advection is not considered. This isolates the analysis to a single advection call; if there are issues with the thickness advection, we presumably cannot expect tracer advection to converge properly. It might have been possible to create an argument to toggle between versions, but it was easier (and arguably less confusing) to just make a copy. This uses a Gaussian bump of 1000 m sitting on top of a 1000 m base layer in a doubly periodic mesh. This ensures a margin-less domain, which would avoid the issue of advection degrading to 1st order at margins.
1 parent 3ee1139 commit 4086e05

File tree

7 files changed

+315
-0
lines changed

7 files changed

+315
-0
lines changed

‎compass/landice/tests/mesh_convergence/__init__.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,9 @@
22
from compass.landice.tests.mesh_convergence.horizontal_advection import (
33
HorizontalAdvection,
44
)
5+
from compass.landice.tests.mesh_convergence.horizontal_advection_thickness import ( # noqa
6+
HorizontalAdvectionThickness,
7+
)
58
from compass.testgroup import TestGroup
69

710

@@ -17,4 +20,5 @@ def __init__(self, mpas_core):
1720
super().__init__(mpas_core=mpas_core, name='mesh_convergence')
1821

1922
self.add_test_case(HorizontalAdvection(test_group=self))
23+
self.add_test_case(HorizontalAdvectionThickness(test_group=self))
2024
self.add_test_case(Halfar(test_group=self))

‎compass/landice/tests/mesh_convergence/horizontal_advection_thickness/__init__.py‎

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
from compass.landice.tests.mesh_convergence.conv_test_case import ConvTestCase
2+
from compass.landice.tests.mesh_convergence.horizontal_advection_thickness.analysis import ( # noqa
3+
Analysis,
4+
)
5+
from compass.landice.tests.mesh_convergence.horizontal_advection_thickness.init import ( # noqa
6+
Init,
7+
)
8+
9+
10+
class HorizontalAdvectionThickness(ConvTestCase):
11+
"""
12+
A test case for testing horizontal advection of thickness in MALI with
13+
planar, doubly periodic meshes
14+
"""
15+
def __init__(self, test_group):
16+
"""
17+
Create test case for creating a MALI mesh
18+
19+
Parameters
20+
----------
21+
test_group : compass.landice.tests.mesh_convergence.MeshConvergence
22+
The landice test group that this test case belongs to
23+
"""
24+
super().__init__(test_group=test_group,
25+
name='horizontal_advection_thickness')
26+
27+
self.add_step(Analysis(test_case=self, resolutions=self.resolutions))
28+
29+
def create_init(self, resolution):
30+
"""
31+
Child class must override this to return an instance of a
32+
ConvInit step
33+
34+
Parameters
35+
----------
36+
resolution : int
37+
The resolution of the test case
38+
39+
Returns
40+
-------
41+
init : compass.landice.tests.mesh_convergence.conv_init.ConvInit
42+
The init step object
43+
"""
44+
45+
return Init(test_case=self, resolution=resolution)
46+
47+
def create_analysis(self, resolutions):
48+
"""
49+
50+
Child class must override this to return an instance of a
51+
ConvergenceInit step
52+
53+
Parameters
54+
----------
55+
resolutions : list of int
56+
The resolutions of the other steps in the test case
57+
58+
Returns
59+
-------
60+
analysis : compass.landice.tests.mesh_convergence.conv_analysis.ConvAnalysis # noqa
61+
The init step object
62+
"""
63+
64+
return Analysis(test_case=self, resolutions=resolutions)

‎compass/landice/tests/mesh_convergence/horizontal_advection_thickness/analysis.py‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
import warnings
2+
3+
import matplotlib.pyplot as plt
4+
import numpy as np
5+
import xarray as xr
6+
7+
from compass.landice.tests.mesh_convergence.conv_analysis import ConvAnalysis
8+
9+
10+
class Analysis(ConvAnalysis):
11+
"""
12+
A step for visualizing the output from the advection convergence test case
13+
"""
14+
def __init__(self, test_case, resolutions):
15+
"""
16+
Create the step
17+
18+
Parameters
19+
----------
20+
test_case : compass.TestCase
21+
The test case this step belongs to
22+
23+
resolutions : list of int
24+
The resolutions of the meshes that have been run
25+
"""
26+
super().__init__(test_case=test_case, resolutions=resolutions)
27+
self.resolutions = resolutions
28+
self.add_output_file('convergence.png')
29+
30+
def run(self):
31+
"""
32+
Run this step of the test case
33+
"""
34+
plt.switch_backend('Agg')
35+
resolutions = self.resolutions
36+
ncells_list = list()
37+
errors = list()
38+
for res in resolutions:
39+
rms_error, ncells = self.rmse(res, variable='thickness')
40+
ncells_list.append(ncells)
41+
errors.append(rms_error)
42+
43+
ncells = np.array(ncells_list)
44+
errors = np.array(errors)
45+
46+
p = np.polyfit(np.log10(ncells), np.log10(errors), 1)
47+
conv = abs(p[0]) * 2.0
48+
49+
error_fit = ncells**p[0] * 10**p[1]
50+
51+
plt.loglog(ncells, error_fit, 'k')
52+
plt.loglog(ncells, errors, 'or')
53+
plt.annotate('Order of Convergence = {}'.format(np.round(conv, 3)),
54+
xycoords='axes fraction', xy=(0.3, 0.95), fontsize=14)
55+
plt.xlabel('Number of Grid Cells', fontsize=14)
56+
plt.ylabel('L2 Norm', fontsize=14)
57+
section = self.config['mesh_convergence']
58+
duration = section.getfloat('duration')
59+
plt.title(f'Thickness horizontal advection convergence test,'
60+
f'{duration} yrs')
61+
plt.savefig('convergence.png', bbox_inches='tight', pad_inches=0.1)
62+
63+
section = self.config['horizontal_advection']
64+
conv_thresh = section.getfloat('conv_thresh')
65+
conv_max = section.getfloat('conv_max')
66+
67+
if conv < conv_thresh:
68+
raise ValueError(f'order of convergence '
69+
f' {conv} < min tolerence {conv_thresh}')
70+
71+
if conv > conv_max:
72+
warnings.warn(f'order of convergence '
73+
f'{conv} > max tolerence {conv_max}')
74+
75+
def rmse(self, resolution, variable):
76+
"""
77+
Compute the RMSE for a given resolution
78+
79+
Parameters
80+
----------
81+
resolution : int
82+
The resolution of the (uniform) mesh in km
83+
84+
variable : str
85+
The name of a variable in the output file to analyze.
86+
87+
Returns
88+
-------
89+
rms_error : float
90+
The root-mean-squared error
91+
92+
ncells : int
93+
The number of cells in the mesh
94+
"""
95+
res_tag = '{}km'.format(resolution)
96+
97+
ds = xr.open_dataset('{}_output.nc'.format(res_tag))
98+
init = ds[variable].isel(Time=0)
99+
final = ds[variable].isel(Time=-1)
100+
diff = final - init
101+
rms_error = np.sqrt((diff**2).mean()).values
102+
ncells = ds.sizes['nCells']
103+
return rms_error, ncells

‎compass/landice/tests/mesh_convergence/horizontal_advection_thickness/horizontal_advection_thickness.cfg‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
# options for planar horizontal advection test case
2+
[horizontal_advection]
3+
4+
# Number of vertical levels
5+
vert_levels = 3
6+
7+
# ice thickness (m)
8+
ice_thickness = 1000.0
9+
10+
# bed elevation (m)
11+
bed_elevation = 0.0
12+
13+
# center of the tracer gaussian (km)
14+
x_center = 0.
15+
y_center = 0.
16+
17+
# width of gaussian tracer "blob" (km)
18+
gaussian_width = 50
19+
20+
# whether to advect in x, y, or both
21+
advect_x = True
22+
advect_y = True
23+
24+
# convergence threshold below which the test fails
25+
conv_thresh = 0.6
26+
27+
# Convergence rate above which a warning is issued
28+
conv_max = 3.0
29+
30+
# options for mesh convergence test cases
31+
[mesh_convergence]
32+
33+
# a list of resolutions (km) to test
34+
resolutions = 16, 8, 4, 2

‎compass/landice/tests/mesh_convergence/horizontal_advection_thickness/init.py‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
import numpy
2+
import xarray
3+
from mpas_tools.io import write_netcdf
4+
5+
from compass.landice.tests.mesh_convergence.conv_init import ConvInit
6+
7+
8+
class Init(ConvInit):
9+
"""
10+
A step for creating an initial_condition for advection convergence test
11+
case
12+
"""
13+
def __init__(self, test_case, resolution):
14+
"""
15+
Create the step
16+
17+
Parameters
18+
----------
19+
test_case : compass.TestCase
20+
The test case this step belongs to
21+
22+
resolution : int
23+
The resolution of the test case
24+
"""
25+
super().__init__(test_case=test_case, resolution=resolution)
26+
self.add_output_file('initial_state.nc')
27+
28+
def run(self):
29+
"""
30+
Run this step of the test case
31+
"""
32+
# create the mesh and graph.info
33+
super().run()
34+
35+
config = self.config
36+
37+
section = config['horizontal_advection']
38+
ice_thickness = section.getfloat('ice_thickness')
39+
bed_elevation = section.getfloat('bed_elevation')
40+
x_center = 1e3 * section.getfloat('x_center')
41+
y_center = 1e3 * section.getfloat('y_center')
42+
advect_x = section.getboolean('advect_x')
43+
advect_y = section.getboolean('advect_y')
44+
gaussian_width = 1e3 * section.getfloat('gaussian_width')
45+
nVertLevels = section.getint('vert_levels')
46+
47+
section = config['mesh_convergence']
48+
duration = section.getfloat('duration') * 3600.0 * 24.0 * 365.0
49+
50+
ds = xarray.open_dataset('mesh.nc')
51+
xCell = ds.xCell
52+
yCell = ds.yCell
53+
54+
if advect_x:
55+
x_vel = ds.attrs['x_period'] / duration
56+
else:
57+
x_vel = 0.
58+
59+
if advect_y:
60+
y_vel = ds.attrs['y_period'] / duration
61+
else:
62+
y_vel = 0.
63+
64+
layerThicknessFractions = xarray.DataArray(
65+
data=1.0 / nVertLevels * numpy.ones((nVertLevels,)),
66+
dims=['nVertLevels'])
67+
68+
# create base layer of uniform thickness over entire domain
69+
# to ensure no margins
70+
thickness = ice_thickness * xarray.ones_like(xCell)
71+
# now add Gaussian bump on top
72+
# using bump height equal to base height
73+
dist_sq = (xCell - x_center)**2 + (yCell - y_center)**2
74+
thickness += ice_thickness * numpy.exp(-0.5 * dist_sq /
75+
gaussian_width**2)
76+
thickness = thickness.expand_dims(dim='Time', axis=0)
77+
78+
bedTopography = bed_elevation * xarray.ones_like(thickness)
79+
80+
uReconstructX = x_vel * xarray.ones_like(xCell)
81+
uReconstructX = uReconstructX.expand_dims(dim={"nVertInterfaces":
82+
nVertLevels + 1},
83+
axis=1)
84+
uReconstructX = uReconstructX.expand_dims(dim='Time', axis=0)
85+
86+
uReconstructY = y_vel * xarray.ones_like(xCell)
87+
uReconstructY = uReconstructY.expand_dims(dim={"nVertInterfaces":
88+
nVertLevels + 1},
89+
axis=1)
90+
uReconstructY = uReconstructY.expand_dims(dim='Time', axis=0)
91+
92+
ds['layerThicknessFractions'] = layerThicknessFractions
93+
ds['thickness'] = thickness
94+
ds['bedTopography'] = bedTopography
95+
ds['uReconstructX'] = uReconstructX
96+
ds['uReconstructY'] = uReconstructY
97+
98+
write_netcdf(ds, 'initial_state.nc')

‎compass/landice/tests/mesh_convergence/horizontal_advection_thickness/namelist.forward‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
config_velocity_solver = 'none'
2+
config_unrealistic_velocity = 1.0e36
3+
config_thickness_advection = 'fo'
4+
config_tracer_advection = 'none'
5+
config_start_time = '0001-01-01_00:00:00'

‎compass/landice/tests/mesh_convergence/horizontal_advection_thickness/streams.forward‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
<streams>
2+
<stream name="output">
3+
<var name="uReconstructX"/>
4+
<var name="uReconstructY"/>
5+
<var name="thickness"/>
6+
</stream>
7+
</streams>

0 commit comments

Comments
 (0)