Applying Filtering and Coarse-Graining to Generate Low-Resolution Datasets#
Now that we’ve run our high-resolution simulations, we can apply the different filtering and coarse-graining operators examined in the paper to generate low-resolution datasets from these simulations.
import xarray as xr
import fsspec
import numpy as np
import matplotlib.pyplot as plt
import pyqg
import pyqg_parameterization_benchmarks.coarsening_ops as coarsening
%matplotlib inline
We must first initialize a high-resolution model. It is from the simulations run on the high-resolution model that we will generate filtered and coarsened low-resolution simulations that approximate or generalize to simulations made on low-resolution initialized models. Below we are using the same eddy-configured high resolution model from the preceding section and running it through without stopping.
high_res_eddy_model = pyqg.QGModel(nx=256, dt=3600.0, tmax=311040000.0, tavestart=155520000.0, log_level=0)
high_res_eddy_model.run()
We now initialize our filtering and coarse-graining operators. Each of the following operators are defined in more detail within the paper. The operator methods take in two parameters: the high resolution model to be filtered and coarsened and the resolution at which to filter and coarsen down to. We will be taking our high resolution simulations with a grid of \(256\times256\) and coarsening them down to low resolution simulations with a grid of \(64\times64\).
op1 = coarsening.Operator1(high_res_eddy_model, 64) # spectral truncation + sharp filter
op2 = coarsening.Operator2(high_res_eddy_model, 64) # spectral truncation + softer Gaussian filter
op3 = coarsening.Operator3(high_res_eddy_model, 64) # GCM-Filters + real-space coarsening
ops = [op1, op2, op3]
INFO: Logger initialized
INFO: Logger initialized
INFO: Logger initialized
We can retrieve the filtered and coarsened, low-resolution simulations as pyqg.QGModel instances from each of the above operators and convert them into xarray.Dataset objects.
op1.m2.to_dataset()
<xarray.Dataset>
Dimensions: (time: 1, lev: 2, y: 64, x: 64, l: 64, k: 33, lev_mid: 1)
Coordinates:
* time (time) float64 0.0
* lev (lev) int64 1 2
* lev_mid (lev_mid) float64 1.5
* x (x) float64 7.812e+03 2.344e+04 3.906e+04 ... 9.766e+05 9.922e+05
* y (y) float64 7.812e+03 2.344e+04 3.906e+04 ... 9.766e+05 9.922e+05
* l (l) float64 0.0 6.283e-06 1.257e-05 ... -1.257e-05 -6.283e-06
* k (k) float64 0.0 6.283e-06 1.257e-05 ... 0.0001948 0.0002011
Data variables: (12/14)
q (time, lev, y, x) float64 1.576e-05 5.876e-06 ... -1.64e-06
u (time, lev, y, x) float64 0.008275 -0.02175 ... 0.008176 0.01083
v (time, lev, y, x) float64 0.05972 0.1077 ... -0.01099 -0.01272
ufull (time, lev, y, x) float64 0.03327 0.003253 ... 0.008176 0.01083
vfull (time, lev, y, x) float64 0.05972 0.1077 ... -0.01099 -0.01272
qh (time, lev, l, k) complex128 (0.002316570289261702+0j) ... (-8.8...
... ...
ph (time, lev, l, k) complex128 0j ... (-1.2670009605120546e-14+1.5...
dqhdt (time, lev, l, k) complex128 0j 0j 0j 0j 0j 0j ... 0j 0j 0j 0j 0j
Ubg (lev) float64 0.025 0.0
Qy (lev) float64 1.039e-10 -7.222e-12
dqdt (time, lev, y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
p (time, lev, y, x) float64 -1.952e+03 -560.5 999.5 ... 304.7 109.0
Attributes: (12/23)
pyqg:beta: 1.5e-11
pyqg:delta: 0.25
pyqg:del2: 0.8
pyqg:dt: 3600.0
pyqg:filterfac: 23.6
pyqg:L: 1000000.0
... ...
pyqg:tc: 0
pyqg:tmax: 1576800000.0
pyqg:twrite: 1000.0
pyqg:W: 1000000.0
title: pyqg: Python Quasigeostrophic Model
reference: https://pyqg.readthedocs.io/en/latest/index.htmlVisualizing Effects#
We can visualize the effects on various state variables and resulting subgrid forcing terms from applying the three different methods of filtering and coarse-graining respectively to high resolution eddy configurations.
# Helper function to display plots
def imshow(arr, vlim=3e-5):
plt.xticks([]); plt.yticks([])
return plt.imshow(arr, vmin=-vlim, vmax=vlim, cmap='bwr', interpolation='none')
# Visualizing effects on upper-level PV
fig = plt.figure(figsize=(15.5,3))
plt.subplot(1,4,1, title='High-res PV')
imshow(high_res_eddy_model.q[0])
for j, op in enumerate(ops):
plt.subplot(1,4,j+2, title=op.__class__.__name__)
im = imshow(op.m2.q[0])
fig.colorbar(im, ax=fig.axes, pad=0.02).set_label('Upper PV [$s^{-1}$]')
From each of the resulting coarsened PV snapshots, it is difficult to see the differences that distinguish them. We can turn to looking at the subgrid forcing to better depict the contrast. We can attain the forcing by calling subgrid_forcing() on the operator instance which takes in a state variable ( q, u, or v) as a string, to compute the subgrid forcing of.
# Visualizing effects on upper PV forcing term
fig = plt.figure(figsize=(14.5,3))
plt.subplot(1,4,1, title='Operator1 $S_{q_{1,tot}}$')
imshow(op1.q_forcing_total[0], 3e-11)
for j, op in enumerate(ops):
plt.subplot(1,4,j+2, title=op.__class__.__name__ + " $S_{q_1}$")
im = imshow(op.subgrid_forcing('q')[0], 5e-11)
cb = fig.colorbar(im, ax=fig.axes, pad=0.02).set_label('$S_{q_1}$ [$s^{-2}$]')
We can also compute subgrid flux terms with respect to \(u\) and \(v\), the \(x\)- and \(y\)- velocities relative to the background flow, by calling subgrid_fluxes() which similarily also takes in a state variable, as a string, to compute the subgrid flux of.
# Visualizing x-component of the lower-layer PV subgrid flux
fig = plt.figure(figsize=(14.5,3))
for j, op in enumerate(ops):
plt.subplot(1,4,j+1, title=op.__class__.__name__ + " $\phi_{q_2}$")
uq, vq = op.subgrid_fluxes('q')
im = imshow(uq[1], 1.5e-8)
fig.colorbar(im, ax=fig.axes, pad=0.02).set_label('$\phi_{q_2}$ [$m\,s^{-2}$]', fontsize=14)
Utilizing Preexisting Datasets#
As promised, we will discuss the structure of the remaining sub-directories housing the datasets.
eddy/
high_res.zarr
low_res.zarr
forcing1.zarr
forcing2.zarr
forcing3.zarr
jet/
high_res.zarr
low_res.zarr
forcing1.zarr
forcing2.zarr
forcing3.zarr
low_res.zarr contains snapshots and diagnostics for low resolution eddy- and jet-configured models (where nx=64). We can run low resolution simulations in the same manner in which we’ve previously ran our high resolution simulations with the modification of changing the argument value for the nx parameter defined in the pyqg.QGModel class when initializing our model.
forcing{1,2,3}.zarr contain high resolution snapshots that have been filtered and coarsened down to low resolution, along with relevant information on subgrid forcing variables. There are three of these datasets due to the fact three different filtering and coarse-graining operations on high resolution simulations were explored in the paper. Each dataset corresponds respectively to the three filtering and coarse-graining operators we initialized.
# Load a small subset of the subgrid forcing data
eddy_forcing1 = xr.open_zarr('gs://leap-persistent/pbluc/eddy/forcing1').isel(run=0).load() # Computed with Operator 1
eddy_forcing2 = xr.open_zarr('gs://leap-persistent/pbluc/eddy/forcing2').isel(run=0).load() # Computed with Operator 2
Subgrid Forcings#
Each of the forcing datasets, include additional variables that detail the underlying deviations between the high resolution model and its coarsened complement (where \(\overline{(\,)}\) denotes the coarse-graining operator):
q_subgrid_forcing: \(S_q \equiv \overline{(\mathbf{u} \cdot \nabla)q} - (\overline{\mathbf{u}} \cdot \overline{\nabla})\overline{q}\)u_subgrid_forcing: \(S_u \equiv \overline{(\mathbf{u} \cdot \nabla)u} - (\overline{\mathbf{u}} \cdot \overline{\nabla})\overline{u}\)v_subgrid_forcing: \(S_v \equiv \overline{(\mathbf{u} \cdot \nabla)v} - (\overline{\mathbf{u}} \cdot \overline{\nabla})\overline{v}\)uq_subgrid_flux: \(\phi_{uq} \equiv \overline{uq} - \bar{u}\bar{q}\)vq_subgrid_flux: \(\phi_{vq} \equiv \overline{vq} - \bar{v}\bar{q}\).uu_subgrid_flux: \(\phi_{uu} \equiv \overline{u^2} - \bar{u}^2\)vv_subgrid_flux: \(\phi_{vv} \equiv \overline{v^2} - \bar{v}^2\)uv_subgrid_flux: \(\phi_{uv} \equiv \overline{uv} - \bar{u}\bar{v}\).dqdt_bar: PV tendency from the high-resolution model, filtered and coarsened to low resolutiondqbar_dt: PV tendency from the low-resolution model, initialized at \(\overline{q}\).
plt.figure(figsize=(12,5)).suptitle("Subgrid forcing with filtering and coarse-graining schemes")
plt.subplot(121); eddy_forcing1.q_subgrid_forcing.isel(lev=1, time=-1).plot(
cbar_kwargs=dict(label="Upper PV Forcing [$s^{-2}$]"))
plt.title("eddy/forcing1")
plt.subplot(122); eddy_forcing2.q_subgrid_forcing.isel(lev=1, time=-1).plot(
cbar_kwargs=dict(label="Upper PV Forcing [$s^{-2}$]"))
plt.title("eddy/forcing2")
plt.tight_layout()
We can construct our own forcing datasets with these additional variables included as follows:
def append_subgrid_variables(snapshot, op):
subgrid_forcing_terms = [('q_subgrid_forcing', 'q'), ('u_subgrid_forcing', 'u'), ('v_subgrid_forcing', 'v')]
subgrid_flux_terms = [('uq_subgrid_flux', 'q', 0), ('vq_subgrid_flux', 'q', 1), ('uu_subgrid_flux', 'u', 0), ('uv_subgrid_flux', 'v', 0), ('vv_subgrid_flux', 'v', 1)]
snapshot['q_forcing_total'] = (['lev', 'y', 'x'], op.q_forcing_total)
for term in subgrid_forcing_terms:
snapshot[term[0]] = (['lev', 'y', 'x'], op.subgrid_forcing(term[1]))
for term in subgrid_flux_terms:
snapshot[term[0]] = (['lev', 'y', 'x'], op.subgrid_fluxes(term[1])[term[2]])
return snapshot
hires_model = pyqg.QGModel(nx=256, dt=3600.0, tmax=311040000.0, tavestart=155520000.0)
op1 = coarsening.Operator1(hires_model, 64)
initial_snapshot = op1.m2.to_dataset()
snapshots = [append_subgrid_variables(initial_snapshot, op1)]
for _ in hires_model.run_with_snapshots(tsnapint=1000*hires_model.dt):
op1 = coarsening.Operator1(hires_model, 64)
snapshot = op1.m2.to_dataset()
snapshots.append(append_subgrid_variables(snapshot, op1))
eddy_forcing1 = xr.concat(snapshots, dim='time')