source: codes/icosagcm/devel/Python/test/py/bubble.py @ 617

import math as math
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as manimation

import dynamico.dyn as dyn
import dynamico.time_step as time_step
import dynamico.DCMIP as DCMIP
import dynamico.advect as advect
from dynamico import unstructured as unst

def reldiff(a,b):
    return (a-b)/(np.absolute(a).max()+np.absolute(b).max()+1e-20)

def go(model,courant,scheme,start):
    def plotter(metric, flow, scalars, t):
        mik, gas_ik, ujk, Phi_il, Wil = model.diagnose(*flow)
        Q1, Q2 = scalars
        
        mil = metric.ml_ext(mik)
        wil = Wil/mil*metric.g0
        wik = metric.mk_ext(Wil)/mik*metric.g0
        xik = metric.xik/1000
        xjl = metric.xjl/1000
        xil = metric.xil/1000
        zil = Phi_il/metric.g0/1000
        zik = metric.mk_int(zil)
        zjl = metric.mj(zil)
        print 's ',np.max(gas_ik.s),np.min(gas_ik.s)
        print 'ujk ',np.max(ujk),np.min(ujk)
        print 'wil ',np.max(wil),np.min(wil)
        sik = gas_ik.s
        sik = .5*(abs(sik)+sik)

        plt.clf()
        plt.contourf(xik, zik, sik, 20 )
#        plt.contourf(xik, zik, wik, 20 )
        plt.xlim((-1,1)), plt.xlabel('x (km)')
        plt.ylim((0,ztop/1000.)), plt.ylabel('z (km)')
        plt.colorbar()
        plt.title("t = %f s" % (t))
        plt.savefig('fig_bubble/bubble.png',bbox_inches='tight')
        writer.grab_frame()
        
    transport = advect.SplitTransport(model.metric, advect.HorizVanLeer, advect.VertVanLeer)
    replace = lambda fd,fv : fv
    return DCMIP.run_with_scalars(model,scheme,transport,replace, plotter, courant, start,scalars, T1, N1, N2)

Lx, ztop, nx, llm = 2000.,1000., 100, 50
nqdyn, ny, dx = 1, 1, Lx/nx
T1, N1, N2 = 5., 1, 200
#T1, N1, N2 = 5., 2, 5

metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.thermal_bubble(Lx,nx,llm, ztop=ztop, zc=350.)
Sik = m0ik*gas0.s
start = (m0ik, Sik, u0_jk, Phi0_il, 0*Phi0_il)
scalars = (Sik, m0ik*metric.mk_int(Phi0_il))

model = dyn.Nonhydro(thermo,metric,BC)
solver, caldyn_eta = dyn.MassBasedNH(model), unst.eta_mass
#solver, caldyn_eta = dyn.LagrangianNH(model), unst.eta_lag

rho_bot=BC.rho_bot
if False: # explicit time-stepping with rigid lid
    def rigid(flow,tau): # impose rigid BCs at top and bottom
        flow, fast, slow = solver.bwd_fast_slow(flow,0.)
        fast[4][:,0]=0.
        fast[4][:,-1]=0.
        return flow,fast,slow
    BC.rho_bot = 100*rho_bot # moderate stiffness
    courant, time_scheme = 2.0, lambda dt : time_step.RK4(rigid, dt)
else: # HEVI time-stepping with pressure BC
    # Note : parameters BC.rigid_top and BC.rigid_bottom are True by default, but are ignored by the Python code...
    BC.rigid_top=False
    BC.rigid_bottom=False
    BC.rho_bot = 1e6*rho_bot # large stiffness
#    BC.rho_bot = 100*rho_bot # moderate stiffness
    courant, time_scheme = 3.0, lambda dt : time_step.ARK2(solver.bwd_fast_slow, dt, a32=0.7)

FFMpegWriter = manimation.writers['ffmpeg']
metadata = dict(title='Aquaplanet', artist='DYNAMICO_LMDZ5',
                comment='Movie support!')
writer = FFMpegWriter(fps=15, metadata=metadata, bitrate=4000, codec="h263p")
fig = plt.figure(figsize=(12,5))
with writer.saving(fig, "fig_bubble/bubble.avi", 100):
    m0ik, gas0_ik, u0jk, Phi0_il, W0il = go(model, courant, time_scheme, start)

#-------------------------------------------------------------------------------------------------------
# Now use final state to check that Fortran and Python implementations produce the exact same tendencies
#-------------------------------------------------------------------------------------------------------

#W0il[:,0]=0.
#W0il[:,-1]=0.

unst.setvar('nb_threads', 1)
mesh = unst.Cartesian_mesh(nx,ny,llm,nqdyn,Lx,ny*dx,0.)
xx_ik, xx_il, ll = mesh.xx[:,0,:]/1000, mesh.xxp1[:,0,:]/1000, mesh.ll[:,0,:]

mass_bl,mass_dak,mass_dbk=unst.compute_hybrid_coefs(m0ik)
unst.init_hybrid(mass_bl,mass_dak,mass_dbk)

caldyn_thermo = unst.thermo_entropy
#caldyn_thermo = unst.thermo_theta
g = mesh.dx/metric.dx_g0
caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g)
unst.setvar('rigid',False)

s0ik = gas0_ik.s
if caldyn_thermo == unst.thermo_theta:
    s0ik = thermo.T0*np.exp(s0ik/thermo.Cpd)

# 2 momentum components in DYNAMICO vs 1 component in slice
u0jk = u0jk + 100*dx # shift horizontal wind by 100 m/s
u0_ker = mesh.field_u()
mesh.set_ucomp(u0_ker,u0jk)

tau=0.1
# impose hybdrid distribution of mass (Fortran)
flow = (m0ik/dx,s0ik*m0ik/dx,u0_ker,Phi0_il,W0il/dx) # NB : different definitions of mass field in DYNAMICO / Python
flow,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
# shift Phi from target value Phi_bot
Phi0_il = Phi0_il+g*100.*np.exp(-10.*xx_il*xx_il)
m,S,u,Phi,W = flow
flow = m,S,u,Phi0_il,W
# compute tendencies (0=Python)
flow0 = (m*dx,S*dx,u0jk,Phi0_il,W0il, 0., 0.) # NB : different definitions of mass field in DYNAMICO / Python
flow0, fast0, slow0 = solver.bwd_fast_slow(flow0,tau,True) 
# compute tendencies (Fortran)
flow,fast,slow = caldyn.bwd_fast_slow(flow, tau) 


zz_il = Phi0_il/metric.g0/1000.
zz_ik=.5*(zz_il[:,0:-1]+zz_il[:,1:])
print
print 'ptop, model top (km) :', unst.getvar('ptop'), zz_il.max()
        
for data0,data,dname in ((flow0,flow,'flow'),(slow0,slow,'slow'),(fast0,fast,'fast')):
    def plot(xx,zz,var,vname):
        title = "%s_%s" % (vname,dname)
        filename = "bubble/bubble_%s.png" % title
        var = 1.*var
        if var.size>1:
            plt.figure(figsize=(12,5))
            plt.contourf(xx,zz,var,20)
#                    plt.pcolor(xx,zz,var)
            plt.xlim((-1.,1.)), plt.xlabel('x (km)')
            plt.ylim((0.,1.)), plt.ylabel('z (km)')        
            plt.colorbar()
            plt.title(title)
            plt.savefig(filename,bbox_inches='tight')
            print title, var.min(), var.max()
            if var.max()>1e100:
                raise OverflowError
                    
    dm,dS,du,dPhi,dW = data    
    du=mesh.ucomp(du)/dx
    du=du[0,:,:]

    dm0, dS0, du0, dPhi0, dW0, junk, junk = data0
    dm0, dS0, du0, dW0 = dm0/dx, dS0/dx, du0/dx, dW0/dx

    if dname == 'flow':
        hflux=mesh.ucomp(caldyn.hflux)
        plot(xx_ik, zz_ik, dm,'m')
        plot(xx_ik, zz_ik, du,'u')
        plot(xx_ik, zz_ik, du-du0,'uu')
        plot(xx_ik, zz_ik, dS/m,'s')
        plot(xx_il, zz_il, dW,'w')
        plot(xx_il, zz_il, dW-dW0,'ww')
    else:
        dS=reldiff(dS,dS0)
        dW=dW-dW0
        dPhi=dPhi-dPhi0
        du=reldiff(du,du0)
        plot(xx_ik, zz_ik, dS,'dS')
        plot(xx_il, zz_il, dPhi,'dPhi')
        plot(xx_il, zz_il, dW,'dW')
        plot(xx_ik, zz_ik, du,'du')

if unst.getvar('rigid'):
    print 'Fortran BCs : rigid'
else:
    print 'Fortran BCs : pressure/soft', unst.getvars('ptop','rho_bot','Phi_bot','pbot')

if BC.rigid_top:
    print 'Top BC : rigid'
else:
    print 'Top BC : pressure', BC.ptop

if BC.rigid_bottom:
    print 'Bottom BC : rigid'
else:
    print 'Bottom BC : soft', BC.rho_bot[0], BC.Phi_bot[0], BC.pbot[0]