<jupyter_start><jupyter_text>Standalone Convergence Checker for the numerical vKdV solverCopied from Standalone Convergence Checker for the numerical KdV solver - just add bathyDoes not save or require any input data<jupyter_code>import xarray as xr
from iwaves.kdv.kdvimex import KdVImEx#from_netcdf
from iwaves.kdv.vkdv import vKdV 
from iwaves.kdv.solve import solve_kdv
from iwaves.utils.plot import vKdV_plot
import iwaves.utils.initial_conditions as ics

import numpy as np
from scipy.interpolate import PchipInterpolator as pchip
import matplotlib.pyplot as plt

%matplotlib inline

from matplotlib import rcParams

# Set font sizes
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Bitstream Vera Sans']
rcParams['font.serif'] = ['Bitstream Vera Sans']
rcParams["font.size"] = "14"
rcParams['axes.labelsize']='large'
# CONSTANTS FOR WHOLE NOTEBOOK
d = 252.5
L_d = 4.0e5
Nz = 100

# Functions

def run_kdv(args):
    """
    Main function for generating different soliton scenarios
    """
    rho_params, bathy_params, a0, L_d, mode, nu_H, dx, runtime, dt, Lw = args
    ####################################################
    # Inputs
    mode = 0
    
    Nz = 100

    ntout = 1800.0

    z = np.linspace(0, -d, Nz)

    dz = np.abs(z[1]-z[0])
    
    x = np.arange(-2*dx,L_d+dx,dx)
        
    h = ics.depth_tanh2(bathy_params, x)

    kdvargs = dict(\
        verbose=False,\
        a0=a0,\
        Lw=Lw,\
        mode=mode,
        dt=dt,\
        nu_H=nu_H,\
        ekdv=False,\
        wavefunc=ics.eta_fullsine,\
        #L_d = L_d,
        x=x,\
        Nsubset=10,
        nonlinear=False,
        nonhydrostatic=False
    )

    ###
#     THIS WAS COPIED FROM THE KdV VERSION. IT INITIALISES EACH vKdV 3 TIMES - QUITE SLOW. 
    ###
    
    ii=0
    #rhoz = single_tanh_rho(
    #        z, pp['rho0'][ii], pp['drho1'][ii], pp['z1'][ii], pp['h1'][ii])
    rhoz = ics.rho_double_tanh_rayson(rho_params,z)
    ######

    ## Call the vKdV run function
    mykdv, Bda = solve_kdv(rhoz, z, runtime,\
            solver='vkdv', h=h, ntout=ntout, outfile=None, **kdvargs)

    print('Done with dx={} and dt={}'.format(dx, dt))
    
    return mykdv, Bda
dx = 10
x = np.arange(-2*dx,L_d+dx,dx)
bathy_params = [L_d*0.5, 50000, d+50, d-50]       

h = ics.depth_tanh2(bathy_params, x)
    
plt.figure(figsize=(9,5))
plt.plot(x, h, 'k')
plt.ylabel('h (m)')
plt.xlabel('x (m)')
plt.title('vKdV bathy')
#betas = [1023.7, 1.12, 105, 52, 155, 43] # ~April 5data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAkAAAAFdCAYAAADrOtLuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3gVVf7H8fc3JISEJhEwu1Q7CKhI7GBFsdd1VViwYkFxVZAisiDCAqIgP0XWhrAWXGyoWFAQhAVRKUsTsABSpEVKBCIlOb8/ZoKXaxISUuaWz+t55rk3M+fO/R6G8uHMnBlzziEiIiISTxKCLkBERESkvCkAiYiISNxRABIREZG4owAkIiIicUcBSEREROKOApCIiIjEHQUgEYk4ZtbXzEr1Hh1m1tDMnJn1KM39FvJ9o83st/L4LhEpPgUgESlVZnazHzROC1ufamaTzWyvmbUtxv5O8vc3qJA2h5rZbjN7uyS1F5eZ1fXD2onl+b0iUnIKQCJS5swsFfgAOBvo4Jx7vaifdc7NBb4F2pqZFdDsr0AS8GpJay2mukAfQAFIJMooAIlImQoLP+2LE35CvAbUA1oVsL0dsAX48KCKFJG4owAkIgdkZn/xT0OdX8xt4eFnbD5tLjWz+Wb2m5l9Z2a35VPCa4DDCzrhn28InAG86ZzbXcT+dDazFWaWbWazzOzUsO0NzOxZM1tqZjvNbKuZfWBmTULanAN86f/4sv9r4Mysb9i+DjOz/5hZlpltMbPnzaxSyPYZZja/gDq/MbPZRemTiBSPApCIFMUE4Ffghny23QBsAKaGrQ8NP38rIPycB7wHJAO9gVeAfwJXhbZzzv0ETAf+YmYVw3bTFjCKfvqrLdANeB7oD9QHJpnZMSFtTvbrfge4H3jKXzfNzNL9NkuAvv7754H2/vJOyH4SgE+B3f53jgc6Av8IaTMaON7Mjg8t0q8nA+/XRERKm3NOixYtWg644P1D/AuQFLKuCrATeDpk3c14ozUrgb3ADYXscw6wCUgLWdfY/5wLa9vR3+8VYesX+99lB6i/of/5XUDDkPXHAHuA10PWpebz+aOA34BeIetO8/d5cz7tR/vb+oetHw9sCvm5OpANPB7W7jG/rtpBH3stWmJx0QiQiBTVG0AacEHIuiuBFH9buMPwwsbK/Hbmj6ScBPzbObc5b71zbgkwMZ+PvOnvb99pMH/21XHAa865ok6bf985t68m59x3/vddHLJuZ8h3pJrZocBW4DugRRG/J8+zYT9/AdQ0s6r+d23DC0VtzSzB/07D6+enzrmNxfw+ESkCBSARKapPgc3sfxrsBmA1MDOf9ncB24GPwk/v+Br6r8vy2faHdc65rXin4i7PCw/8HoaKM/vr+3zWfQccYmbVAcyskpk9bmY/AzuATLyRqmbAIcX4rj3OuZ/D1m3xX9NC1o0B6gDn+j+fCRyOTn+JlBkFIBEpEufcHuBt4Eo/IBwCXAj8p4DRl2VAG/zrYMzs6LDteVPa8/tsQdPdX8UbcbraHy25EZjrjxoVuStF+L7hQBdgHN4U+zZ4I1+LKd7fm7mFbAv9zs+AdcDf/J//hnfN1XvF+C4RKQYFIBEpjjeAasAlwLVARfI//<KEY>///<KEY>
#betas = [1023.5, 1.22, 67, 55, 157, 52] # ~March 1

betas_w = [1023.8229810318612,
     0.9865506702797462,
     143.5428700089361,
     46.1265812512485,
     136.66278860120943,
     41.57014327398592] # 15 July 2016

betas_s =[1023.6834358117951,
     1.2249066117658955,
     156.78804559089772,
     53.66835548728355,
     73.14183287436342,
     40.21031777315428] # 1st April 2017

a0 = 20.
mode =0
nu_H = 0
runtime = 2.5*86400.

# Going to make Lw an input for the vKdV as it will really speed things up. 
dx = 100
dt = 10

z = np.linspace(0, -d, Nz)
rhoz_w = ics.rho_double_tanh_rayson(betas_w, z)
rhoz_s = ics.rho_double_tanh_rayson(betas_s, z)
Lw_w = ics.get_Lw(rhoz_w, z, z0=max(h), mode=0)
Lw_s = ics.get_Lw(rhoz_s, z, z0=max(h), mode=0)

print(Lw_w)
print(Lw_s)
dxs =[1600,800,400,200,100,75,50,37.5,25]
dxs =[800,400,200,100,75,50,35]
dxs =[800,200,100,50]
dt = 8.

all_kdv_dx_w = []
all_kdv_dx_s = []

for dx in dxs:
    
    print(' ')
    print('Running dx={}'.format(dx))
    print(' ')
    
    mykdv, B = run_kdv( (betas_w, bathy_params, a0, L_d, mode, nu_H, dx, runtime, dt, Lw_w))
    all_kdv_dx_w.append(mykdv)
    
    mykdv, B = run_kdv( (betas_s, bathy_params, a0, L_d, mode, nu_H, dx, runtime, dt, Lw_s))
    all_kdv_dx_s.append(mykdv)
    
    print(' ')
    print('Completed dx={}'.format(dx))
    print(' ')
plt.figure(figsize=(9,5))
for mykdv in all_kdv_dx_s:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dx_s)
# plt.xlim((162200, 163600))
plt.legend()
plt.show()

plt.figure(figsize=(9,5))
for mykdv in all_kdv_dx_s:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dx_s)
# plt.xlim((162200, 163600))
plt.ylim((-65, 40))
plt.xlim((250000, 300000))
plt.legend()

plt.figure(figsize=(9,5))
for mykdv in all_kdv_dx_w:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dx_s)
plt.legend()
plt.show()

plt.figure(figsize=(9,5))
for mykdv in all_kdv_dx_w:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dx_s)
    
plt.legend()
plt.ylim((-40, 10))
plt.xlim((250000, 300000))
# Compute the errors
X = np.arange(0,L_d, 10.)
nx = X.shape[0]
ndx = len(dxs)

solns = np.zeros((ndx, nx))
for ii, mykdv in enumerate(all_kdv_dx_w):
    Fx = pchip(mykdv.x, mykdv.B)
    solns[ii,:] = Fx(X)

# Compute the error between each solution
#err = np.diff(solns, axis=0)
err = solns - solns[-1,:]

err_rms_w = np.linalg.norm(err, ord=2, axis=1) # L2-norm
#err_rms_w = np.sqrt(np.mean(err**2,axis=1))

solns = np.zeros((ndx, nx))
for ii, mykdv in enumerate(all_kdv_dx_s):
    Fx = pchip(mykdv.x, mykdv.B)
    solns[ii,:] = Fx(X)

# Compute the error between each solution
#err = np.diff(solns, axis=0)
err = solns - solns[-1,:]

err_rms_s = np.linalg.norm(err, ord=2, axis=1) # L2-norm
#err_rms_s = np.sqrt(np.mean(err**2,axis=1))


plt.figure(figsize=(9,8))
plt.loglog(dxs[:-1],err_rms_s[:-1],'ko')
plt.loglog(dxs[:-1],err_rms_w[:-1],'s', color='0.5')
plt.xlim(2e1,2e3)
plt.ylim(1e-2,2e3)
plt.grid(b=True)
x0 = np.array([50,100.])

plt.plot(x0, 100/x0[0]**2*x0**2, 'k--')
plt.plot(x0, 100/x0[0]**1*x0**1, 'k:')

plt.ylabel('L2-norm Error [m]')
plt.xlabel('$\Delta x$ [m]')

alpha_s = -2*all_kdv_dx_s[0].c1*all_kdv_dx_s[0].r10 
beta_s = -1*all_kdv_dx_s[0].r01
alpha_w = -2*all_kdv_dx_w[0].c1*all_kdv_dx_w[0].r10 
beta_w = -1*all_kdv_dx_w[0].r01

plt.legend((r'~\Deltax$^1$',
            r'$\alpha$ = (%3.4f,%3.4f), $\beta$ = (%3.4f,%3.4f)'%(min(alpha_s), max(alpha_s), min(beta_s), max(beta_s)),
           r'$\alpha$ = (%3.4f,%3.4f), $\beta$ = (%3.4f,%3.4f)'%(min(alpha_w), max(alpha_w), min(beta_w), max(beta_w))), loc='lower right')
# Delta t comparison
dts = [20,10.,5,2.5,1.25,0.6,0.3]
dx = 50.

all_kdv_dt_w = []
all_kdv_dt_s = []

for dt in dts:
    
    print(' ')
    print('Running dt={}'.format(dt))
    print(' ')
    
    mykdv, B = run_kdv( (betas_w, bathy_params, a0, L_d, mode, nu_H, dx, runtime, dt, Lw_w))
    all_kdv_dt_w.append(mykdv)
    
    mykdv, B = run_kdv( (betas_s, bathy_params, a0, L_d, mode, nu_H, dx, runtime, dt, Lw_s))
    all_kdv_dt_s.append(mykdv)
    
    print(' ')
    print('Completed dt={}'.format(dt))
    print(' ')
    
plt.figure(figsize=(9,5))
for mykdv in all_kdv_dt_s:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dt_s)
plt.legend()
plt.show()

plt.figure(figsize=(9,5))
for mykdv in all_kdv_dt_s:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dt_s)
    
plt.legend()
plt.ylim((-50, 30))
plt.xlim((195000, 210000))
plt.figure(figsize=(9,5))
for mykdv in all_kdv_dt_w:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dt_s)
plt.legend()
plt.show()

plt.figure(figsize=(9,5))
for mykdv in all_kdv_dt_w:
    plt.plot(mykdv.x, mykdv.B, label=mykdv.dt_s)
    
plt.legend()
plt.ylim((-30, 1))
plt.xlim((175000, 205000))
# Compute the errors
X = np.arange(0,L_d, 10.)
nx = X.shape[0]
ndx = len(dts)

solns = np.zeros((ndx, nx))
for ii, mykdv in enumerate(all_kdv_dt_w):
    print(ii)
    Fx = pchip(mykdv.x, mykdv.B)
    solns[ii,:] = Fx(X)

# Compute the error between each solution
#err = np.diff(solns, axis=0)
err = solns - solns[-1,:]

err_rms_w_t = np.linalg.norm(err, ord=2, axis=1) # L2-norm
#err_rms_w = np.sqrt(np.mean(err**2,axis=1))

solns = np.zeros((ndx, nx))
for ii, mykdv in enumerate(all_kdv_dt_s):
    print(ii)
    Fx = pchip(mykdv.x, mykdv.B)
    solns[ii,:] = Fx(X)

# Compute the error between each solution
#err = np.diff(solns, axis=0)
err = solns - solns[-1,:]

err_rms_s_t = np.linalg.norm(err, ord=2, axis=1) # L2-norm
#err_rms_s = np.sqrt(np.mean(err**2,axis=1))
plt.figure(figsize=(12,8))
ax=plt.subplot(121)
plt.loglog(dxs[:-1],err_rms_s[:-1],'ko', markersize=6)
plt.loglog(dxs[:-1],err_rms_w[:-1],'s', color='0.5', markersize=4)
plt.xlim(2e1,2e3)
plt.ylim(1e0,2e3)
plt.grid(b=True)
x0 = np.array([50,100.])
plt.plot(x0, 100/x0[0]**2*x0**2, 'k--')
plt.plot(x0, 100/x0[0]**1*x0**1, 'k:')
plt.ylabel('L2-norm Error [m]')
plt.xlabel('$\Delta x$ [m]')

alpha_s = -2*all_kdv_dx_s[0].c1*all_kdv_dx_s[0].r10 
beta_s = -1*all_kdv_dx_s[0].r01
alpha_w = -2*all_kdv_dx_w[0].c1*all_kdv_dx_w[0].r10 
beta_w = -1*all_kdv_dx_w[0].r01

plt.legend((r'$\alpha$ = (%3.3f, %3.3f), $\beta$ = (%3.0f, %3.0f)'%(min(alpha_s), max(alpha_s), min(beta_s), max(beta_s)),
           r'$\alpha$ = (%3.3f, %3.3f), $\beta$ = (%3.0f, %3.0f)'%(min(alpha_w), max(alpha_w), min(beta_w), max(beta_w))), loc='lower right')

plt.text(0.05,0.95,'(a)',transform=ax.transAxes)


ax=plt.subplot(122)
plt.loglog(dts[:-1],err_rms_s_t[:-1],'kd', markersize=6)
plt.loglog(dts[:-1],err_rms_w_t[:-1],'s', color='0.5', markersize=4)
plt.xlim(0,0.5e2)
plt.ylim(1e-2,1e3)
plt.grid(b=True)
x0 = np.array([5,20])
plt.plot(x0, 10/x0[0]**2*x0**2, 'k--')
plt.plot(x0, 10/x0[0]**1*x0**1, 'k:')
#plt.ylabel('L2-norm Error [m]')
plt.xlabel('$\Delta t$ [s]')

plt.text(0.05,0.95,'(b)',transform=ax.transAxes)
alpha_s = -2*all_kdv_dt_s[0].c1*all_kdv_dt_s[0].r10 
beta_s = -1*all_kdv_dt_s[0].r01
alpha_w = -2*all_kdv_dt_w[0].c1*all_kdv_dt_w[0].r10 
beta_w = -1*all_kdv_dt_w[0].r01

plt.legend((r'$\alpha$ = (%3.3f, %3.3f), $\beta$ = (%3.0f, %3.0f)'%(min(alpha_s), max(alpha_s), min(beta_s), max(beta_s)),
           r'$\alpha$ = (%3.3f, %3.3f), $\beta$ = (%3.0f, %3.0f)'%(min(alpha_w), max(alpha_w), min(beta_w), max(beta_w))), loc='lower right')

plt.savefig('../FIGURES/vkdv_convergence_dxdt.png',dpi=150)
plt.savefig('../FIGURES/vkdv_convergence_dxdt.pdf',dpi=150)<jupyter_output>C:\ProgramData\Anaconda3\lib\site-packages\ipykernel_launcher.py:28: UserWarning: Attempted to set non-positive left xlim on a log-scaled axis.
Invalid limit will be ignored.
{'fraction_non_alphanumeric': 0.1784361005856011, 'fraction_numerical': 0.08887357905614882, 'mean_word_length': 4.653846153846154, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '502962', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 6678, 'n_tokens_neox': 6043, 'n_words': 824, 'path': 'tests/DELETE standalone_vkdv_convergence linearHydrostatic.ipynb', 'repo_name': 'iosonobert/iwaves', 'size': '31206'}