Skip to content

Commit

Permalink
Assume mesh variables only in mesh file
Browse files Browse the repository at this point in the history
  • Loading branch information
cbegeman committed Jul 21, 2023
1 parent 2202e23 commit 7dbf393
Show file tree
Hide file tree
Showing 11 changed files with 65 additions and 52 deletions.
6 changes: 3 additions & 3 deletions polaris/ocean/tests/baroclinic_channel/init.py
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ def run(self):

temperature = temperature.expand_dims(dim='Time', axis=0)

normal_velocity = xr.zeros_like(ds.xEdge)
normal_velocity = xr.zeros_like(ds_mesh.xEdge)
normal_velocity, _ = xr.broadcast(normal_velocity, ds.refBottomDepth)
normal_velocity = normal_velocity.transpose('nEdges', 'nVertLevels')
normal_velocity = normal_velocity.expand_dims(dim='Time', axis=0)
Expand All @@ -149,8 +149,8 @@ def run(self):
ds['salinity'] = salinity * xr.ones_like(temperature)
ds['normalVelocity'] = normal_velocity
ds['fCell'] = coriolis_parameter * xr.ones_like(x_cell)
ds['fEdge'] = coriolis_parameter * xr.ones_like(ds.xEdge)
ds['fVertex'] = coriolis_parameter * xr.ones_like(ds.xVertex)
ds['fEdge'] = coriolis_parameter * xr.ones_like(ds_mesh.xEdge)
ds['fVertex'] = coriolis_parameter * xr.ones_like(ds_mesh.xVertex)

ds.attrs['nx'] = nx
ds.attrs['ny'] = ny
Expand Down
1 change: 0 additions & 1 deletion polaris/ocean/tests/baroclinic_channel/rpe/analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -149,5 +149,4 @@ def _plot(nx, ny, lx, ly, filename, nus, rpe):
ax.set_ylabel('y, km')
if iCol == num_files - 1:
fig.colorbar(dis, ax=axs[num_files - 1], aspect=40)

plt.savefig(filename)
6 changes: 3 additions & 3 deletions polaris/ocean/tests/baroclinic_channel/viz.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@ def __init__(self, test_case):
"""
super().__init__(test_case=test_case, name='viz')
self.add_input_file(
filename='initial_state.nc',
target='../init/initial_state.nc')
filename='mesh.nc',
target='../init/culled_mesh.nc')
self.add_input_file(
filename='output.nc',
target='../forward/output.nc')
Expand All @@ -32,7 +32,7 @@ def run(self):
"""
Run this step of the test case
"""
ds_mesh = xr.load_dataset('initial_state.nc')
ds_mesh = xr.load_dataset('mesh.nc')
ds = xr.load_dataset('output.nc')
t_index = ds.sizes['Time'] - 1
plot_horiz_field(ds, ds_mesh, 'temperature',
Expand Down
18 changes: 10 additions & 8 deletions polaris/ocean/tests/global_convergence/cosine_bell/analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,9 @@ def __init__(self, test_case, resolutions, icosahedral):
mesh_name = f'Icos{resolution}'
else:
mesh_name = f'QU{resolution}'
self.add_input_file(
filename=f'{mesh_name}_mesh.nc',
target=f'../{mesh_name}/mesh/mesh.nc')
self.add_input_file(
filename=f'{mesh_name}_init.nc',
target=f'../{mesh_name}/init/initial_state.nc')
Expand Down Expand Up @@ -127,7 +130,8 @@ def rmse(self, mesh_name):
psi0 = config.getfloat('cosine_bell', 'psi0')
pd = config.getfloat('cosine_bell', 'vel_pd')

init = xr.open_dataset(f'{mesh_name}_init.nc')
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)):
Expand All @@ -142,7 +146,7 @@ def rmse(self, mesh_name):
t = 86400.0 * DY + HR * 3600. + MN
# find new location of blob center
# center is based on equatorial velocity
R = init.sphere_radius
R = ds_mesh.sphere_radius
distTrav = 2.0 * 3.14159265 * R / (86400.0 * pd) * t
# distance in radians is
distRad = distTrav / R
Expand All @@ -151,9 +155,9 @@ def rmse(self, mesh_name):
newLon -= 2.0 * np.pi

# construct analytic tracer
tracer = np.zeros_like(init.tracer1[0, :, 0].values)
latC = init.latCell.values
lonC = init.lonCell.values
tracer = np.zeros_like(ds_init.tracer1[0, :, 0].values)
latC = ds_mesh.latCell.values
lonC = ds_mesh.lonCell.values
temp = R * np.arccos(np.sin(latCent) * np.sin(latC) +
np.cos(latCent) * np.cos(latC) * np.cos(
lonC - newLon))
Expand All @@ -165,6 +169,4 @@ def rmse(self, mesh_name):
tracerF = ds.tracer1[sliceTime, :, 0].values
rmseValue = np.sqrt(np.mean((tracerF - tracer)**2))

init.close()
ds.close()
return rmseValue, init.dims['nCells']
return rmseValue, ds_init.dims['nCells']
22 changes: 11 additions & 11 deletions polaris/ocean/tests/global_convergence/cosine_bell/init.py
Original file line number Diff line number Diff line change
Expand Up @@ -51,21 +51,21 @@ def run(self):
section = config['vertical_grid']
bottom_depth = section.getfloat('bottom_depth')

dsMesh = xr.open_dataset('mesh.nc')
angleEdge = dsMesh.angleEdge
latCell = dsMesh.latCell
latEdge = dsMesh.latEdge
lonCell = dsMesh.lonCell
sphere_radius = dsMesh.sphere_radius
ds_mesh = xr.open_dataset('mesh.nc')
angleEdge = ds_mesh.angleEdge
latCell = ds_mesh.latCell
latEdge = ds_mesh.latEdge
lonCell = ds_mesh.lonCell
sphere_radius = ds_mesh.sphere_radius

ds = dsMesh.copy()
ds = ds_mesh.copy()

ds['bottomDepth'] = bottom_depth * xr.ones_like(latCell)
ds['ssh'] = xr.zeros_like(latCell)

init_vertical_coord(config, ds)

temperature_array = temperature * xr.ones_like(ds.latCell)
temperature_array = temperature * xr.ones_like(ds_mesh.latCell)
temperature_array, _ = xr.broadcast(temperature_array, ds.refZMid)
ds['temperature'] = temperature_array.expand_dims(dim='Time', axis=0)
ds['salinity'] = salinity * xr.ones_like(ds.temperature)
Expand All @@ -92,8 +92,8 @@ def run(self):
velocity_array, _ = xr.broadcast(velocity, ds.refZMid)
ds['normalVelocity'] = velocity_array.expand_dims(dim='Time', axis=0)

ds['fCell'] = xr.zeros_like(ds.xCell)
ds['fEdge'] = xr.zeros_like(ds.xEdge)
ds['fVertex'] = xr.zeros_like(ds.xVertex)
ds['fCell'] = xr.zeros_like(ds_mesh.xCell)
ds['fEdge'] = xr.zeros_like(ds_mesh.xEdge)
ds['fVertex'] = xr.zeros_like(ds_mesh.xVertex)

write_netcdf(ds, 'initial_state.nc')
8 changes: 6 additions & 2 deletions polaris/ocean/tests/global_convergence/cosine_bell/viz.py
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,9 @@ def __init__(self, test_case, name, subdir, viz_map, mesh_name):
The name of the mesh
"""
super().__init__(test_case=test_case, name=name, subdir=subdir)
self.add_input_file(
filename='mesh.nc',
target='../init/mesh.nc')
self.add_input_file(
filename='initial_state.nc',
target='../init/initial_state.nc')
Expand All @@ -109,12 +112,13 @@ def run(self):

remapper = viz_map.get_remapper()

ds_mesh = xr.open_dataset('mesh.nc')
ds_init = xr.open_dataset('initial_state.nc')
ds_init = ds_init[['tracer1', ]].isel(Time=0, nVertLevels=0)
ds_init = remapper.remap(ds_init)
ds_init.to_netcdf('remapped_init.nc')

plot_global(ds_init.lon.values, ds_init.lat.values,
plot_global(ds_mesh.lon.values, ds_mesh.lat.values,
ds_init.tracer1.values,
out_filename='init.png', config=config,
colormap_section='cosine_bell_viz',
Expand All @@ -125,7 +129,7 @@ def run(self):
ds_out = remapper.remap(ds_out)
ds_out.to_netcdf('remapped_final.nc')

plot_global(ds_out.lon.values, ds_out.lat.values,
plot_global(ds_mesh.lon.values, ds_mesh.lat.values,
ds_out.tracer1.values,
out_filename='final.png', config=config,
colormap_section='cosine_bell_viz',
Expand Down
10 changes: 5 additions & 5 deletions polaris/ocean/tests/inertial_gravity_wave/init.py
Original file line number Diff line number Diff line change
Expand Up @@ -70,14 +70,14 @@ def run(self):

ds = ds_mesh.copy()

ds['ssh'] = xr.zeros_like(ds.xCell)
ds['bottomDepth'] = bottom_depth * xr.ones_like(ds.xCell)
ds['ssh'] = xr.zeros_like(ds_mesh.xCell)
ds['bottomDepth'] = bottom_depth * xr.ones_like(ds_mesh.xCell)

init_vertical_coord(config, ds)

ds['fCell'] = f0 * xr.ones_like(ds.xCell)
ds['fEdge'] = f0 * xr.ones_like(ds.xEdge)
ds['fVertex'] = f0 * xr.ones_like(ds.xVertex)
ds['fCell'] = f0 * xr.ones_like(ds_mesh.xCell)
ds['fEdge'] = f0 * xr.ones_like(ds_mesh.xEdge)
ds['fVertex'] = f0 * xr.ones_like(ds_mesh.xVertex)

ds_mesh['maxLevelCell'] = ds.maxLevelCell
exact_solution = ExactSolution(ds, config)
Expand Down
14 changes: 9 additions & 5 deletions polaris/ocean/tests/inertial_gravity_wave/viz.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,9 @@ def __init__(self, test_case, resolutions):
self.resolutions = resolutions

for resolution in resolutions:
self.add_input_file(
filename=f'mesh_{resolution}km.nc',
target=f'../{resolution}km/initial_state/culled_mesh.nc')
self.add_input_file(
filename=f'init_{resolution}km.nc',
target=f'../{resolution}km/init/initial_state.nc')
Expand All @@ -63,9 +66,10 @@ def run(self):
fig, axes = plt.subplots(nrows=nres, ncols=3, figsize=(12, 2 * nres))
rmse = []
for i, res in enumerate(resolutions):
init = xr.open_dataset(f'init_{res}km.nc')
ds_mesh = xr.open_dataset(f'mesh_{res}km.nc')
ds_init = xr.open_dataset(f'init_{res}km.nc')
ds = xr.open_dataset(f'output_{res}km.nc')
exact = ExactSolution(init, config)
exact = ExactSolution(ds_init, config)

t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(),
'%Y-%m-%d_%H:%M:%S')
Expand All @@ -81,13 +85,13 @@ def run(self):
if i == 0:
error_range = np.max(np.abs(ds.ssh_error.values))

plot_horiz_field(ds, init, 'ssh', ax=axes[i, 0],
plot_horiz_field(ds, ds_mesh, 'ssh', ax=axes[i, 0],
cmap='cmo.balance', t_index=ds.sizes["Time"] - 1,
vmin=-eta0, vmax=eta0, cmap_title="SSH (m)")
plot_horiz_field(ds, init, 'ssh_exact', ax=axes[i, 1],
plot_horiz_field(ds, ds_mesh, 'ssh_exact', ax=axes[i, 1],
cmap='cmo.balance',
vmin=-eta0, vmax=eta0, cmap_title="SSH (m)")
plot_horiz_field(ds, init, 'ssh_error', ax=axes[i, 2],
plot_horiz_field(ds, ds_mesh, 'ssh_error', ax=axes[i, 2],
cmap='cmo.balance',
cmap_title=r"$\Delta$ SSH (m)",
vmin=-error_range, vmax=error_range)
Expand Down
10 changes: 5 additions & 5 deletions polaris/ocean/tests/manufactured_solution/init.py
Original file line number Diff line number Diff line change
Expand Up @@ -70,14 +70,14 @@ def run(self):

ds = ds_mesh.copy()

ds['ssh'] = xr.zeros_like(ds.xCell)
ds['bottomDepth'] = bottom_depth * xr.ones_like(ds.xCell)
ds['ssh'] = xr.zeros_like(ds_mesh.xCell)
ds['bottomDepth'] = bottom_depth * xr.ones_like(ds_mesh.xCell)

init_vertical_coord(config, ds)

ds['fCell'] = coriolis_parameter * xr.ones_like(ds.xCell)
ds['fEdge'] = coriolis_parameter * xr.ones_like(ds.xEdge)
ds['fVertex'] = coriolis_parameter * xr.ones_like(ds.xVertex)
ds['fCell'] = coriolis_parameter * xr.ones_like(ds_mesh.xCell)
ds['fEdge'] = coriolis_parameter * xr.ones_like(ds_mesh.xEdge)
ds['fVertex'] = coriolis_parameter * xr.ones_like(ds_mesh.xVertex)

ds_mesh['maxLevelCell'] = ds.maxLevelCell

Expand Down
14 changes: 9 additions & 5 deletions polaris/ocean/tests/manufactured_solution/viz.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,9 @@ def __init__(self, test_case, resolutions):
self.resolutions = resolutions

for resolution in resolutions:
self.add_input_file(
filename=f'mesh_{resolution}km.nc',
target=f'../{resolution}km/initial_state/culled_mesh.nc')
self.add_input_file(
filename=f'init_{resolution}km.nc',
target=f'../{resolution}km/init/initial_state.nc')
Expand All @@ -62,9 +65,10 @@ def run(self):
fig, axes = plt.subplots(nrows=nres, ncols=3, figsize=(12, 2 * nres))
rmse = []
for i, res in enumerate(resolutions):
init = xr.open_dataset(f'init_{res}km.nc')
ds_mesh = xr.open_dataset(f'mesh_{res}km.nc')
ds_init = xr.open_dataset(f'init_{res}km.nc')
ds = xr.open_dataset(f'output_{res}km.nc')
exact = ExactSolution(config, init)
exact = ExactSolution(config, ds_init)

t0 = datetime.datetime.strptime(ds.xtime.values[0].decode(),
'%Y-%m-%d_%H:%M:%S')
Expand All @@ -80,13 +84,13 @@ def run(self):
if i == 0:
error_range = np.max(np.abs(ds.ssh_error.values))

plot_horiz_field(ds, init, 'ssh', ax=axes[i, 0],
plot_horiz_field(ds, ds_mesh, 'ssh', ax=axes[i, 0],
cmap='cmo.balance', t_index=ds.sizes["Time"] - 1,
vmin=-eta0, vmax=eta0, cmap_title="SSH")
plot_horiz_field(ds, init, 'ssh_exact', ax=axes[i, 1],
plot_horiz_field(ds, ds_mesh, 'ssh_exact', ax=axes[i, 1],
cmap='cmo.balance',
vmin=-eta0, vmax=eta0, cmap_title="SSH")
plot_horiz_field(ds, init, 'ssh_error', ax=axes[i, 2],
plot_horiz_field(ds, ds_mesh, 'ssh_error', ax=axes[i, 2],
cmap='cmo.balance', cmap_title="dSSH",
vmin=-error_range, vmax=error_range)

Expand Down
8 changes: 4 additions & 4 deletions polaris/ocean/tests/single_column/init.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ def run(self):
write_netcdf(ds_mesh, 'culled_mesh.nc')

ds = ds_mesh.copy()
x_cell = ds.xCell
x_cell = ds_mesh.xCell
bottom_depth = config.getfloat('vertical_grid', 'bottom_depth')
ds['bottomDepth'] = bottom_depth * xr.ones_like(x_cell)
ds['ssh'] = xr.zeros_like(x_cell)
Expand Down Expand Up @@ -116,7 +116,7 @@ def run(self):
salinity = salinity.expand_dims(dim='Time', axis=0)

normal_velocity, _ = xr.broadcast(
xr.zeros_like(ds.xEdge), ds.refBottomDepth)
xr.zeros_like(ds_mesh.xEdge), ds.refBottomDepth)
normal_velocity = normal_velocity.transpose('nEdges', 'nVertLevels')
normal_velocity = normal_velocity.expand_dims(dim='Time', axis=0)

Expand All @@ -127,8 +127,8 @@ def run(self):
ds['salinity'] = salinity
ds['normalVelocity'] = normal_velocity
ds['fCell'] = coriolis_parameter * xr.ones_like(x_cell)
ds['fEdge'] = coriolis_parameter * xr.ones_like(ds.xEdge)
ds['fVertex'] = coriolis_parameter * xr.ones_like(ds.xVertex)
ds['fEdge'] = coriolis_parameter * xr.ones_like(ds_mesh.xEdge)
ds['fVertex'] = coriolis_parameter * xr.ones_like(ds_mesh.xVertex)

ds.attrs['nx'] = nx
ds.attrs['ny'] = ny
Expand Down

0 comments on commit 7dbf393

Please sign in to comment.