latest stand of the code before starting on the revision

the instability structure after an hour resembles the analytical solution, and the initial oscillation in the energy is reduced significantly.

Author: ray-chew

RKLM_Python/inputs/boundary.py

@@ -433,14 +433,18 @@ def get_bottom_tau_y(ud, elem, node, alpha, cutoff=0.5):
     tauc_y = dd * tauc_y / np.abs(tauc_y).max()
     taun_y = dd * taun_y / np.abs(taun_y).max()
 
+    # tauc_y = tauc_y[2:-2]
+    # taun_y = taun_y[2:-2]
+
     return tauc_y, taun_y
 
 
 
-def rayleigh_damping(Sol, mpv, ud, forcing=None):
-    u = Sol.rhou / Sol.rho
-    v = Sol.rhov / Sol.rho
-    Y = Sol.rhoY / Sol.rho
+def rayleigh_damping(Sol, mpv, ud, elem, node, forcing=None):
+    u = (Sol.rhou / Sol.rho)#[elem.i2]
+    v = (Sol.rhov / Sol.rho)#[elem.i2]
+    Y = (Sol.rhoY / Sol.rho)#[elem.i2]
+    rho = Sol.rho#[elem.i2]
 
     if ud.bdry_type[1] == BdryType.RAYLEIGH:
         tcy, tny = ud.tcy, ud.tny
@@ -479,11 +483,13 @@ def rayleigh_damping(Sol, mpv, ud, forcing=None):
     # assuming 2D vertical slice - not dimension agnostic
     u += tcy * (u - ud.u_wind_speed) + c_f * (tcy_f * (u - ud.u_wind_speed) + np.abs(tcy_f) * mfac * u_f)
     v += tcy * (v - ud.v_wind_speed) + c_f * (tcy_f * (v - ud.v_wind_speed) + np.abs(tcy_f) * mfac * v_f)
-    Y += tcy * (Y - mpv.HydroState.Y0.reshape(1,-1)) + c_f * (tcy_f * (Y - mpv.HydroState.Y0.reshape(1,-1)) + np.abs(tcy_f) * mfac * Y_f)
 
-    Sol.rhou[...] = Sol.rho * u
-    Sol.rhov[...] = Sol.rho * v
-    Sol.rhoY[...] = Sol.rho * Y
+    Ybar = mpv.HydroState.Y0.reshape(1,-1)
+    Y += tcy * (Y - Ybar) + c_f * (tcy_f * (Y - Ybar) + np.abs(tcy_f) * mfac * Y_f)
+
+    Sol.rhou[...] = rho * u
+    Sol.rhov[...] = rho * v
+    Sol.rhoY[...] = rho * Y
 
 
 

RKLM_Python/inputs/lamb_wave_perturb.py

@@ -123,9 +123,9 @@ def __init__(self):
 
 
         self.tout = [36.0]
-        # self.tout = np.arange(0,36.1,0.1)
+        # self.tout = np.arange(0,361,1.0)
         # self.tout = np.append(self.tout, [720.0])
-        self.stepmax = 10000
+        self.stepmax = 31
         self.output_timesteps = True
 
         self.output_base_name = "_mark_wave"
@@ -138,7 +138,7 @@ def __init__(self):
         self.init_forcing = self.forcing
 
         self.rayleigh_forcing = True
-        self.rayleigh_forcing_type = 'file' # func or file
+        self.rayleigh_forcing_type = 'func' # func or file
         self.rayleigh_forcing_fn = 'output_mark_wave_ensemble=1_601_240_bottom_forcing_S16.h5'
         self.rayleigh_forcing_path = './output_mark_wave'
 
@@ -218,12 +218,12 @@ def eigenfunction(self, t, s, grid='c'):
             eigval, eigvec = np.linalg.eig( self.T_matrix )
 
             # Find index of eigenvalue 
-            # with greatest real part aka the insrtability growth rate
+            # with greatest real part aka the instability growth rate
             ind = np.argmax( np.real( eigval ) )
 
             # construct solution according to eq. 2.27 and 2.19
             exponentials = np.exp( 1j * self.k * x + self.mu * z 
-                                + ( eigval[ind] ) * t )
+                                + ( eigval[ind] ) * (t) + 1j * s * t )
             chi_u  = self.ampl * np.real( eigvec[0,ind] * exponentials )
             chi_w  = self.ampl * np.real( eigvec[1,ind] * exponentials )
             chi_th = self.ampl * np.real( eigvec[2,ind] * exponentials )
@@ -252,7 +252,8 @@ def dehatter(self, th, grid='c'):
 def sol_init(Sol, mpv, elem, node, th, ud, seeds=None):
 
     if hasattr(ud, 'rayleigh_bdry_switch'):
-        ud.bdry_type[1] = BdryType.RAYLEIGH
+        if ud.rayleigh_bdry_switch:
+            ud.bdry_type[1] = BdryType.RAYLEIGH
 
     if ud.bdry_type[1] == BdryType.RAYLEIGH:
         ud.tcy, ud.tny = get_tau_y(ud, elem, node, 0.5)
@@ -264,7 +265,7 @@ def sol_init(Sol, mpv, elem, node, th, ud, seeds=None):
 
     A0 = 1.0e-3
     Msq = ud.Msq
-    g = ud.gravity_strength[1]
+    g = ud.gravity_strength[1] * ud.Rg
 
     x = elem.x.reshape(-1,1)
     y = elem.y.reshape(1,-1)
@@ -299,7 +300,11 @@ def sol_init(Sol, mpv, elem, node, th, ud, seeds=None):
     ud.Cs = Cs
     ud.Ns = N
     # dimensionless Coriolis strength
-    F = ud.coriolis_strength[2] #+ 1e-30
+    if ud.coriolis_strength[2] == 0.0:
+        ud.coriolis_strength[2] += 1e-15
+    F = ud.coriolis_strength[2]
+    # if F == 0.0:
+    #     F += 1e-15
 
     G =  np.sqrt( 9. / 40. )
     Gamma = G * N / Cs
@@ -342,6 +347,7 @@ def sol_init(Sol, mpv, elem, node, th, ud, seeds=None):
     set_explicit_boundary_data(Sol,elem,ud,th,mpv)
 
     if hasattr(ud, 'mixed_run'):
-        ud.coriolis_strength[2] = 2.0 * 7.292 * 1e-5 * ud.t_ref
+        if ud.mixed_run:
+            ud.coriolis_strength[2] = 2.0 * 7.292 * 1e-5 * ud.t_ref
 
     return Sol

RKLM_Python/management/data.py

@@ -184,7 +184,7 @@ def time_update(Sol,flux,mpv,t,tout,ud,elem,node,steps,th,bld=None,writer=None,d
 
         if ud.bdry_type[1] == BdryType.RAYLEIGH: 
             # top rayleight damping
-            boundary.rayleigh_damping(Sol, mpv, ud)
+            boundary.rayleigh_damping(Sol, mpv, ud, elem, node)
 
         # bottom rayleigh forcing
         if hasattr(ud, 'rayleigh_forcing'):
@@ -210,14 +210,16 @@ def time_update(Sol,flux,mpv,t,tout,ud,elem,node,steps,th,bld=None,writer=None,d
                     boundary.rayleigh_damping(Sol, mpv, ud, [up, vp, Yp, pi, t+0.5*dt])
 
                 elif ud.rayleigh_forcing_type == 'func':
-                    boundary.set_explicit_boundary_data(Sol, elem, ud, th, mpv)
-                    ud.rf_bot.eigenfunction((t+0.5*dt), 1)
+                    # boundary.set_explicit_boundary_data(Sol, elem, ud, th, mpv)
+                    
+                    s = 5.0e-3+1e-4+0e-5
+                    ud.rf_bot.eigenfunction((t+0.5*dt), s)
                     up, vp, Yp, pi = ud.rf_bot.dehatter(th)
 
-                    ud.rf_bot.eigenfunction((t+0.5*dt), 1, grid='n')
+                    ud.rf_bot.eigenfunction((t+0.5*dt), s, grid='n')
                     _, _, _, pi_n = ud.rf_bot.dehatter(th, grid='n')
 
-                    boundary.rayleigh_damping(Sol, mpv, ud, [up, vp, Yp, pi_n, t+0.5*dt])
+                    boundary.rayleigh_damping(Sol, mpv, ud, elem, node, [up, vp, Yp, pi_n, t+0.5*dt])
                 boundary.set_explicit_boundary_data(Sol, elem, ud, th, mpv)
 
 
@@ -325,7 +327,7 @@ def time_update(Sol,flux,mpv,t,tout,ud,elem,node,steps,th,bld=None,writer=None,d
 
         if ud.bdry_type[1] == BdryType.RAYLEIGH: 
             # top rayleight damping
-            boundary.rayleigh_damping(Sol, mpv, ud)
+            boundary.rayleigh_damping(Sol, mpv, ud, elem, node)
 
         # bottom rayleigh forcing
         if hasattr(ud, 'rayleigh_forcing'):
@@ -350,16 +352,15 @@ def time_update(Sol,flux,mpv,t,tout,ud,elem,node,steps,th,bld=None,writer=None,d
                     boundary.rayleigh_damping(Sol, mpv, ud, [up, vp, Yp, pi, t+dt])
 
                 elif ud.rayleigh_forcing_type == 'func':
-                    ud.rf_bot.eigenfunction((t+dt), 1)
+
+                    s = 5.0e-3+1e-4+0e-5
+                    ud.rf_bot.eigenfunction((t+dt), s)
                     up, vp, Yp, pi = ud.rf_bot.dehatter(th)
 
-                    ud.rf_bot.eigenfunction((t+dt), 1, grid='n')
+                    ud.rf_bot.eigenfunction((t+dt), s, grid='n')
                     _, _, _, pi_n = ud.rf_bot.dehatter(th, grid='n')
-                    # kernel = np.ones([2]*mpv.p2_nodes.ndim)
-                    # pi_n = np.zeros_like(mpv.p2_nodes)
-                    # pi_n[node.i1] = signal.fftconvolve(pi, kernel, mode='valid') / kernel.sum()
 
-                    boundary.rayleigh_damping(Sol, mpv, ud, [up, vp, Yp, pi_n, t+dt])
+                    boundary.rayleigh_damping(Sol, mpv, ud, elem, node, [up, vp, Yp, pi_n, t+dt])
                 boundary.set_explicit_boundary_data(Sol, elem, ud, th, mpv)
 
 

RKLM_Python/queue_run_lamb_forcing.py

@@ -19,14 +19,15 @@
 resol_x = [301]
 resol_y = [120]
 resol_t = [1,2,4,8,10,12,14,16]
-resol_t = [1]
+resol_t = [10.0]
 # resol_t = [10,12,14,16.0]
 # omegas = [0.0, 2.0 * omega * t_ref]
 # omegas = [2.0 * omega * t_ref]
 omegas = [0.0]
 
 tsteps = [1800, 900, 450, 360, 300, 258, 225]
-tsteps = [3600, 300, 258, 225]
+tsteps = [360, 300, 258, 225]
+# tsteps = [21]
 
 ud = {}
 dap = {

RKLM_Python/queue_run_lamb_mixed.py

@@ -9,18 +9,18 @@
 
 rp.N = 1
 
-# rp.tc = 'lw_p'
-rp.tc = 'mark'
+rp.tc = 'lw_p'
+# rp.tc = 'mark'
 
 t_ref = 100.0
 omega = 7.292 * 1e-5
 # resol_x = [151,301,601,1201]
 # resol_y = [15,30,60,120]
-resol_x = [151]
-resol_y = [60]
+resol_x = [301]
+resol_y = [120]
 resol_t = [1,2,4,8,10,12,14,16]
 # resol_t = [200,400,600,800]
-resol_t = [1.0]
+resol_t = [10.0]
 # omegas = [0.0, 2.0 * omega * t_ref]
 # omegas = [2.0 * omega * t_ref]
 omegas = [0.0]
@@ -37,9 +37,9 @@
 
     ud['ymax'] = 8.0
     ud['rayleigh_forcing'] = False
-    ud['mixed_run'] = False
-    ud['do_advection'] = False
-    # ud['output_timesteps'] = False
+    ud['mixed_run'] = True
+    # ud['do_advection'] = False
+    ud['output_timesteps'] = False
 
 
     for t in resol_t:
@@ -49,12 +49,12 @@
         ud['rayleigh_forcing_fn'] = 'output_mark_wave_ensemble=1_%i_%i_36.000000_bottom_forcing_S%i.h5' %(x,int(4*y),t)
         ud['rayleigh_forcing_type'] = 'func'
 
-        for om in omegas:
+    for om in omegas:
             ud['coriolis_strength'] = [0.0, 0.0, om]
             if om > 0:
                 ud['aux'] = 'test_run_S%i_a05' %t
             else:
-                ud['aux'] = 'test_run_S%i' %t
+                ud['aux'] = 'test_run_S%i_mix' %t
 
             print(ud)
             # run simulation

RKLM_Python/queue_run_lamb_unstable.py

@@ -20,16 +20,16 @@
 resol_y = [30]
 # resol_t = [10,12,14,16]
 # resol_t = [200,400,600,800,1000,1200,1400,1600]
-resol_t = [1,2,4,8]
-# resol_t = [10.0]
+# resol_t = [1,2,4,8]
+resol_t = [10.0]
 # omegas = [0.0, 2.0 * omega * t_ref]
 omegas = [2.0 * omega * t_ref]
 # omegas = [0.0]
 
 
 # tsteps = [3600, 1800, 900, 450, 360, 300, 258, 225]
-tsteps = [3600, 1800, 900, 450]
-# tsteps = [361]
+# tsteps = [3600, 1800, 900, 450]
+tsteps = [360]
 
 ud = {}
 dap = {
@@ -40,6 +40,8 @@
     ud['inx'] = x+1
     ud['iny'] = y+1
     ud['rayleigh_bdry_switch'] = True
+    ud['rayleigh_forcing'] = True
+    ud['do_advection'] = True
 
     for t_idx, t in enumerate(resol_t):
         ud['dtfixed0'] = t / t_ref