Error running LayeredModel

[worthamc]

I've run into two errors when running the LayeredModel, and both appear to be linked to array sizes with zonal jet feature #159. @mbueti, can you look into this?

The first error is related to the size of Vbg in the APEgenspec diagnostic in layered_model. Looks like Vbg has the wrong size. This is triggered when a run passes tavestart.

The second error is layered_model _calc_cfl. It looks like Ubg has an extra np.newaxis. This is triggered when a run passes twrite.

These errors can be generated by the 3-layer baroclinic instability code in the documentation, but I'll copy it here anyway. Setting twrite=1 should show the second error.

import numpy as np
import pyqg
from matplotlib import pyplot as plt

L =  1000.e3     # length scale of box    [m]
Ld = 15.e3       # deformation scale      [m]
kd = 1./Ld       # deformation wavenumber [m^-1]
Nx = 64          # number of grid points

H1 = 500.        # layer 1 thickness  [m]
H2 = 1750.       # layer 2
H3 = 1750.       # layer 3

U1 = 0.05          # layer 1 zonal velocity [m/s]
U2 = 0.025         # layer 2
U3 = 0.00          # layer 3

rho1 = 1025.
rho2 = 1025.275
rho3 = 1025.640

rek = 1.e-7       # linear bottom drag coeff.  [s^-1]
f0  = 0.0001236812857687059 # coriolis param [s^-1]
beta = 1.2130692965249345e-11 # planetary vorticity gradient [m^-1 s^-1]

Ti = Ld/(abs(U1))  # estimate of most unstable e-folding time scale [s]
dt = Ti/200.   # time-step [s]
tmax = 20*Ti      # simulation time [s]

m = pyqg.LayeredModel(nx=Nx, nz=3, U = [U1,U2,U3],V = [0.,0.,0.],L=L,f=f0,beta=beta,
                         H = [H1,H2,H3], rho=[rho1,rho2,rho3],rek=rek,
                        dt=dt,tmax=tmax, twrite=5000, tavestart=Ti*10, ntd=1)

sig = 1.e-7
qi = sig*np.vstack([np.random.randn(m.nx,m.ny)[np.newaxis,],
                    np.random.randn(m.nx,m.ny)[np.newaxis,],
                    np.random.randn(m.nx,m.ny)[np.newaxis,]])
m.set_q(qi)

m.run()

[rabernat]

Thanks a lot for the bug report @worthamc! It looks like our test suite missed this scenario. Really appreciate your feedback.

[mbueti]

I’ll try giving this a look at some point next week, but unfortunately won’t be able to sit down with the code before then. Sorry for the trouble, and thanks for both your patience and catching this!

Mike Bueti

On Jun 6, 2019, at 9:43 PM, Cimarron Wortham notifications@github.com wrote:

I've run into two errors when running the LayeredModel, and both appear to be linked to array sizes with zonal jet feature #159 https://github.com/pyqg/pyqg/pull/159. @mbueti https://github.com/mbueti, can you look into this?

The first error is related to the size of Vbg in the APEgenspec diagnostic in layered_model. Looks like Vbg has the wrong size. This is triggered when a run passes tavestart.

The second error is layered_model _calc_cfl. It looks like Ubg has an extra np.newaxis. This is triggered when a run passes twrite.

These errors can be generated by the 3-layer baroclinic instability code in the documentation, but I'll copy it here anyway. Setting twrite=1 should show the second error.

########################

import numpy as np import pyqg from matplotlib import pyplot as plt

L = 1000.e3 # length scale of box [m] Ld = 15.e3 # deformation scale [m] kd = 1./Ld # deformation wavenumber [m^-1] Nx = 64 # number of grid points

H1 = 500. # layer 1 thickness [m] H2 = 1750. # layer 2 H3 = 1750. # layer 3

U1 = 0.05 # layer 1 zonal velocity [m/s] U2 = 0.025 # layer 2 U3 = 0.00 # layer 3

rho1 = 1025. rho2 = 1025.275 rho3 = 1025.640

rek = 1.e-7 # linear bottom drag coeff. [s^-1] f0 = 0.0001236812857687059 # coriolis param [s^-1] beta = 1.2130692965249345e-11 # planetary vorticity gradient [m^-1 s^-1]

Ti = Ld/(abs(U1)) # estimate of most unstable e-folding time scale [s] dt = Ti/200. # time-step [s] tmax = 20*Ti # simulation time [s]

m = pyqg.LayeredModel(nx=Nx, nz=3, U = [U1,U2,U3],V = [0.,0.,0.],L=L,f=f0,beta=beta, H = [H1,H2,H3], rho=[rho1,rho2,rho3],rek=rek, dt=dt,tmax=tmax, twrite=5000, tavestart=Ti*10, ntd=1)

sig = 1.e-7 qi = sig*np.vstack([np.random.randn(m.nx,m.ny)[np.newaxis,], np.random.randn(m.nx,m.ny)[np.newaxis,], np.random.randn(m.nx,m.ny)[np.newaxis,]]) m.set_q(qi)

m.run()

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/pyqg/pyqg/issues/192?email_source=notifications&email_token=AAACAZWWEIGPV6XKNSH72ULPZG4NLA5CNFSM4HVOD7LKYY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4GYE4SYQ, or mute the thread https://github.com/notifications/unsubscribe-auth/AAACAZSWD52BYIKJN5PY4LTPZG4NLANCNFSM4HVOD7LA.

[mbueti]

@rabernat okay, have an update (https://github.com/pyqg/pyqg/pull/193) that should address this issue! I'll take a look about adding some more test coverage here, but in the meantime this should unblock work for @worthamc

[MFJansen]

This issue should be resolved now (zonal jet feature has been reverted). @worthamc please feel free to re-open if issue is not resolved.