<jupyter_start><jupyter_text>solcing the advection diffusion equation $v\cdot\nabla z_{gl} + \epsilon \Delta z_{gl} = 0$ in domain, $z_{gl}=$ on grounding lineboundary conditions- Neumann (or second-type) boundary condition specifies derivatives- Dirichlet boundary condition specifies the values of the solution itself-> we want Dirichlet BC at the grounding line, Neumann BC at the ice shelf front, or maybe it is possible to set Dirichlet there too, as velocities are converging there only anyways steps: 1. grid- need to interpolate velocity and draft data onto common grid- figure out where the data need to live on the grid<jupyter_code>import sys
import numpy as np
import xarray as xr
import matplotlib
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

sys.path.append("..")
from tqdm.notebook import tqdm
from advect import advect_grl
# %matplotlib widget
%matplotlib inline
%config InlineBackend.print_figure_kwargs={'bbox_inches':None}
%load_ext autoreload
%autoreload 2
ds = xr.open_dataset('../../results/PICOP/Dotson.nc')
ds
plt.figure(figsize=(10,10))
ds.v.plot()
ds
ds = xr.open_dataset('../../results/PICOP/Totten.nc')  # created on Cartesius in CESM env in `regrid_melt.ipynb`
ds['u'] = ds.u.fillna(0)
ds['v'] = ds.v.fillna(0)
# ds['u'] = xr.where(np.isnan(ds.u), 0, ds.u)
# ds['v'] = xr.where(np.isnan(ds.v), 0, ds.v)
evo = advect_grl(ds, eps=.04, T=50)
grl_advect = evo.isel(time=-1).drop('time')
grl_advect.name = 'adv_grl'
ds.draft.plot()
plt.figure(figsize=(10,10))
((ds.draft.rolling(x=5).mean()+ds.draft.rolling(y=5).mean())/2).plot()
# (ds.draft-(ds.draft.rolling(x=5).mean()+ds.draft.rolling(y=5).mean())/2).plot()

xr.merge([ds, grl_advect.isel(dict(x=slice(1,-1), y=slice(1,-1)))])
xr.merge([ds, grl_advect])
f = plt.figure(figsize=(15,10))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=116.33, central_latitude=-67))
# ax.set_extent([-179,179,-85,-60])
im = ax.pcolormesh(evo.x, evo.y, evo.isel(time=-1),#.isel({'time':-1,'x':slice(1,-1),'y':slice(1,-1)}),
              transform=ccrs.SouthPolarStereo(true_scale_latitude=-71))
# ax.add_feature(cartopy.feature.LAND)
ax.coastlines()
# xx, yy = np.meshgrid(ds.x, ds.y)
# ax.streamplot(xx, yy, ds.u, ds.v, color='w', transform=ccrs.SouthPolarStereo(true_scale_latitude=-71))#, 
#               linewidth=np.sqrt(ds.u**2+ds.v**2).fillna(0).values/5e2)
ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
plt.colorbar(im, ax=ax)
plt.savefig('../../results/advect_grl_Totten')<jupyter_output><empty_output><jupyter_text>$\frac{1}{N}\Sigma\sqrt{(z_t-z_{t-1})^2}$'  test advection-diffusion algorithm<jupyter_code># test DataArray
nx, ny = 11, 5
mask = np.ones((ny,nx))
kwargs = {'dims':['y','x'], 'coords':{'x':np.linspace(0,1e5,nx), 'y':np.linspace(0,1e5,ny)}}
da3 = xr.DataArray(data=mask , name='mask' , **kwargs)

# test 1: constant slope, constant u, v=0
draft, _ = np.meshgrid(np.linspace(-1000,-500,nx), np.ones((ny)))
grl = np.zeros((ny,nx))
isf = np.zeros((ny,nx))
u = np.ones_like(draft)*1000
v = np.zeros_like(draft)
grl[:,0] = 1
isf[:,-1] = 1
da0 = xr.DataArray(data=draft, name='draft', **kwargs)
da1 = xr.DataArray(data=isf  , name='isf'  , **kwargs)
da2 = xr.DataArray(data=grl  , name='grl'  , **kwargs)
da4 = xr.DataArray(data=u    , name='u'    , **kwargs)
da5 = xr.DataArray(data=v    , name='v'    , **kwargs)
test1 = xr.merge([da0, da1, da2, da3, da4, da5])

# test 2: same as test 1 but in y direction
_, draft = np.meshgrid( np.ones((nx)), np.linspace(-1000,-500,ny))
grl = np.zeros((ny,nx))
isf = np.zeros((ny,nx))
u = np.ones_like(draft)*500
v = np.ones_like(draft)*1000
grl[0,:] = 1
isf[-1,:] = 1
kwargs = {'dims':['y','x'], 'coords':{'x':np.linspace(0,1e5,nx), 'y':np.linspace(0,1e5,ny)}}
da0 = xr.DataArray(data=draft, name='draft', **kwargs)
da1 = xr.DataArray(data=isf  , name='isf'  , **kwargs)
da2 = xr.DataArray(data=grl  , name='grl'  , **kwargs)
da4 = xr.DataArray(data=u    , name='u'    , **kwargs)
da5 = xr.DataArray(data=v    , name='v'    , **kwargs)
test2 = xr.merge([da0, da1, da2, da3, da4, da5])

# test 3: curving grounding line
xx, yy = np.meshgrid(np.linspace(1,0,nx), np.linspace(0,np.pi,ny))
draft = -500-500*np.sin(yy)*xx
grl = np.zeros((ny,nx))
isf = np.zeros((ny,nx))
u = np.ones_like(draft)*1000
v = np.ones_like(draft)*500
grl[:,0] = 1
isf[:,-1] = 1
da0 = xr.DataArray(data=draft, name='draft', **kwargs)
da1 = xr.DataArray(data=isf  , name='isf'  , **kwargs)
da2 = xr.DataArray(data=grl  , name='grl'  , **kwargs)
da4 = xr.DataArray(data=u    , name='u'    , **kwargs)
da5 = xr.DataArray(data=v    , name='v'    , **kwargs)
test3 = xr.merge([da0, da1, da2, da3, da4, da5])
evo = advect_grl(test1, eps=0, T=100)
plt.figure()
evo.isel(time=-1).plot()
evo = advect_grl(test2, eps=0, T=100)
plt.figure()
evo.isel(time=-1).plot()
evo = advect_grl(test3, eps=100, T=500)
plt.figure()
evo.isel(time=-1).plot()
f, ax = plt.subplots(10, 5, sharey=True, sharex=True, figsize=(10,12))
for i in range(50):
    evo.isel(time=i).plot(ax=ax[int(i/5),i%5])#vmin=-1000, vmax=-500)
# test 4: diffusion only, curving grounding line
xx, yy = np.meshgrid(np.linspace(1,0,nx), np.linspace(0,np.pi,ny))
draft = -500-500*np.sin(yy)*xx
grl = np.zeros((ny,nx))
isf = np.zeros((ny,nx))
u = np.zeros_like(draft)
v = np.zeros_like(draft)
grl[:,0] = 1
isf[:,-1] = 1
da0 = xr.DataArray(data=draft, name='draft', **kwargs)
da1 = xr.DataArray(data=isf  , name='isf'  , **kwargs)
da2 = xr.DataArray(data=grl  , name='grl'  , **kwargs)
da4 = xr.DataArray(data=u    , name='u'    , **kwargs)
da5 = xr.DataArray(data=v    , name='v'    , **kwargs)
test4 = xr.merge([da0, da1, da2, da3, da4, da5])

evo = advect_grl(test4, eps=1e7, T=20)
plt.figure()
evo.isel(time=-1).plot()<jupyter_output><empty_output><jupyter_text>test: realistic geometry, uniform velocity<jupyter_code># real Totten IS geometry
daT = xr.open_dataset(f'../../results/Totten.nc')
daT['u'] = xr.full_like(daT.draft,  500)
daT['v'] = xr.full_like(daT.draft, -500)
evo = advect_grl(daT, eps=.04, T=200)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20685657502484267, 'fraction_numerical': 0.04306061609804571, 'mean_word_length': 5.168539325842697, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '393349', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2754, 'n_tokens_neox': 2613, 'n_words': 489, 'path': 'src/notebooks/Advection.ipynb', 'repo_name': 'AJueling/melt', 'size': '6038'}