[761] | 1 | from dynamico import unstructured as unst |
---|
| 2 | from dynamico import dyn |
---|
| 3 | from dynamico import time_step |
---|
| 4 | from dynamico import DCMIP |
---|
| 5 | from dynamico import meshes |
---|
| 6 | from dynamico import xios |
---|
| 7 | from dynamico import precision as prec |
---|
| 8 | from dynamico.meshes import Cartesian_mesh as Mesh |
---|
| 9 | |
---|
| 10 | import math as math |
---|
| 11 | import numpy as np |
---|
| 12 | import time |
---|
[764] | 13 | import argparse |
---|
[774] | 14 | import logging |
---|
[761] | 15 | from numpy import pi, log, exp, sin, cos |
---|
| 16 | |
---|
| 17 | # Baroclinic instability test based on Ullrich et al. 2015, QJRMS |
---|
| 18 | |
---|
[764] | 19 | parser = argparse.ArgumentParser() |
---|
| 20 | |
---|
[771] | 21 | parser.add_argument("--mpi_ni", type=int, default=1, |
---|
[764] | 22 | help="number of x processors") |
---|
[771] | 23 | parser.add_argument("--mpi_nj", type=int, default=1, |
---|
[764] | 24 | help="number of y processors") |
---|
[774] | 25 | parser.add_argument("--T", type=float, default=3600., |
---|
[769] | 26 | help="Length of time slice in seconds") |
---|
[774] | 27 | parser.add_argument("--N", type=int, default=24, |
---|
| 28 | help="Number of time slices") |
---|
| 29 | parser.add_argument("--log",type=str,default='DEBUG',choices=('DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'), |
---|
| 30 | help="Set verbosity level of logging") |
---|
| 31 | parser.add_argument("--Davies_N1", type=int, default=3) |
---|
| 32 | parser.add_argument("--Davies_N2", type=int, default=3) |
---|
[769] | 33 | |
---|
[771] | 34 | |
---|
[764] | 35 | args = parser.parse_args() |
---|
| 36 | |
---|
| 37 | |
---|
[761] | 38 | def baroclinic_3D(Lx,nx,Ly,ny,llm,ztop=25000.): |
---|
| 39 | Rd = 287.0 # Gas constant for dryy air (j kg^-1 K^-1) |
---|
| 40 | T0 = 288.0 # Reference temperature (K) |
---|
| 41 | lap = 0.005 # Lapse rate (K m^-1) |
---|
| 42 | b = 2. # Non dimensional vertical width parameter |
---|
| 43 | u0 = 35. # Reference zonal wind speed (m s^-1) |
---|
| 44 | a = 6.371229e6 # Radius of the Earth (m) |
---|
| 45 | ptop = 2000. |
---|
| 46 | y0 = Ly*0.5 |
---|
| 47 | Cpd = 1004.5 |
---|
| 48 | p0 = 1e5 |
---|
[774] | 49 | up = 1. # amplitude (m/s) |
---|
| 50 | xc,yc,lp = 0.,Ly/2.,600000. |
---|
[761] | 51 | |
---|
| 52 | omega = 7.292e-5 # Angular velocity of the Earth (s^-1) |
---|
[771] | 53 | phi0 = 90.*np.pi/180.0 # Reference latitude North pi/4 (deg) |
---|
[761] | 54 | f0 = 2*omega*np.sin(phi0) |
---|
| 55 | beta0 = 2*omega*np.cos(phi0)/a |
---|
| 56 | fb = 2*omega*np.sin(phi0) - y0*2*omega*np.cos(phi0)/a |
---|
| 57 | |
---|
[769] | 58 | def Phi_xy(y): |
---|
[761] | 59 | fc = y*y - (Ly*y/pi)*sin(2*pi*y/Ly) |
---|
| 60 | fd = Ly*Ly/(2*pi*pi)*cos(2*pi*y/Ly) |
---|
| 61 | return .5*u0*( fb*(y-y0-Ly/(2*pi)*sin(2*pi*y/Ly)) + .5*beta0*(fc-fd-(Ly*Ly/3.)- Ly*Ly/(2*pi*pi)) ) |
---|
| 62 | |
---|
[769] | 63 | def Phi_xyeta(y,eta): return T0*g/lap*(1-eta**(Rd*lap/g)) + Phi_xy(y)*log(eta)*exp(-((log(eta)/b)**2)) |
---|
[774] | 64 | def ulon(x,y,eta): |
---|
| 65 | u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b))) |
---|
| 66 | # u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2)) |
---|
| 67 | return u |
---|
| 68 | |
---|
[761] | 69 | def tmean(eta) : return T0*eta**(Rd*lap/g) |
---|
[769] | 70 | def T(y,eta) : return tmean(eta)+(Phi_xy(y)/Rd)*(((2/(b*b))*(log(eta))**2)-1)*exp(-((0.5*log(eta))**2)) |
---|
[761] | 71 | def p(eta): return p0*eta # eta = p/p_s |
---|
| 72 | |
---|
| 73 | def eta(alpha) : return (1-(lap*ztop*alpha/(T0)))**(g/(Rd*lap)) # roughly equispaced levels |
---|
[774] | 74 | def coriolis(x,y): return f0+beta0*y |
---|
[761] | 75 | |
---|
| 76 | filename = 'cart_%03d_%03d.nc'%(nx,ny) |
---|
[774] | 77 | logging.info('Reading Cartesian mesh ...') |
---|
[761] | 78 | meshfile = meshes.DYNAMICO_Format(filename) |
---|
| 79 | pmesh = meshes.Unstructured_PMesh(comm,meshfile) |
---|
[763] | 80 | pmesh.partition_curvilinear(args.mpi_ni,args.mpi_nj) |
---|
[769] | 81 | nqdyn, radius = 1, None |
---|
[761] | 82 | mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis) |
---|
| 83 | |
---|
[769] | 84 | |
---|
[761] | 85 | alpha_k = (np.arange(llm) +.5)/llm |
---|
| 86 | alpha_l = (np.arange(llm+1)+ 0.)/llm |
---|
[774] | 87 | x_ik, alpha_ik = np.meshgrid(mesh.lon_i, alpha_k, indexing='ij') |
---|
| 88 | y_ik, alpha_ik = np.meshgrid(mesh.lat_i+.5*Ly, alpha_k, indexing='ij') |
---|
| 89 | x_il, alpha_il = np.meshgrid(mesh.lon_i, alpha_l, indexing='ij') |
---|
| 90 | y_il, alpha_il = np.meshgrid(mesh.lat_i+.5*Ly, alpha_l, indexing='ij') |
---|
| 91 | x_ek, alpha_ek = np.meshgrid(mesh.lon_e, alpha_k, indexing='ij') |
---|
| 92 | y_ek, alpha_ek = np.meshgrid(mesh.lat_e+.5*Ly, alpha_k, indexing='ij') |
---|
[761] | 93 | |
---|
[774] | 94 | print('ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g))) |
---|
[761] | 95 | print(np.shape(alpha_k),np.shape(alpha_l)) |
---|
| 96 | thermo = dyn.Ideal_perfect(Cpd, Rd, p0, T0) |
---|
| 97 | |
---|
| 98 | eta_il = eta(alpha_il) |
---|
| 99 | eta_ik = eta(alpha_ik) |
---|
| 100 | eta_ek = eta(alpha_ek) |
---|
| 101 | print('min max eta_il', np.min(eta_il),np.max(eta_il)) |
---|
| 102 | |
---|
[769] | 103 | Phi_il = Phi_xyeta(y_il, eta_il) |
---|
| 104 | Phi_ik = Phi_xyeta(y_ik, eta_ik) |
---|
[761] | 105 | p_ik = p(eta_ik) |
---|
[769] | 106 | T_ik = T(y_ik, eta_ik) #ik full level(40), il(41) |
---|
[761] | 107 | |
---|
| 108 | gas = thermo.set_pT(p_ik,T_ik) |
---|
| 109 | mass_ik = mesh.field_mass() |
---|
| 110 | for l in range(llm): |
---|
| 111 | mass_ik[:,l]=(Phi_il[:,l+1]-Phi_il[:,l])/(g*gas.v[:,l]) |
---|
[774] | 112 | Sik, Wil = gas.s*mass_ik, mesh.field_w() |
---|
[761] | 113 | |
---|
| 114 | u_ek = mesh.ucov3D(ulon(x_ek, y_ek, eta_ek), 0.*eta_ek) |
---|
| 115 | |
---|
[774] | 116 | print('ztop (m) = ', Phi_il[0,-1]/g, ztop) |
---|
[761] | 117 | ptop = p(eta(1.)) |
---|
[774] | 118 | print( 'ptop (Pa) = ', gas.p[0,-1], ptop) |
---|
[761] | 119 | |
---|
| 120 | params=dyn.Struct() |
---|
| 121 | params.ptop=ptop |
---|
| 122 | params.dx=dx |
---|
| 123 | params.dx_g0=dx/g |
---|
| 124 | params.g = g |
---|
| 125 | |
---|
| 126 | # define parameters for lower BC |
---|
| 127 | pbot = p(eta_il[:,0]) |
---|
[774] | 128 | print( 'min p, T :', pbot.min(), tmean(pbot/p0)) |
---|
[761] | 129 | gas_bot = thermo.set_pT(pbot, tmean(pbot/p0)) |
---|
| 130 | params.pbot = gas_bot.p |
---|
| 131 | params.rho_bot = 1e6/gas_bot.v |
---|
| 132 | |
---|
[774] | 133 | return thermo, mesh, params, prec.asnum([mass_ik,Sik,u_ek,Phi_il,Wil]), gas |
---|
[761] | 134 | |
---|
[769] | 135 | def diagnose(Phi,S,m,W): |
---|
| 136 | s=S/m ; s=.5*(s+abs(s)) |
---|
| 137 | for l in range(llm): |
---|
| 138 | v[:,l]=(Phi[:,l+1]-Phi[:,l])/(g*m[:,l]) |
---|
| 139 | w[:,l]=.5*params.g*(W[:,l+1]+W[:,l])/m[:,l] |
---|
| 140 | z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g |
---|
| 141 | gas = thermo.set_vs(v,s) |
---|
| 142 | return gas, w, z |
---|
| 143 | |
---|
[771] | 144 | class Davies: |
---|
| 145 | def __init__(self,N1,N2,x_i,y_i,x_e,y_e): |
---|
| 146 | self.N1, self.N2 = N1, N2 |
---|
| 147 | self.beta_i = self.mask(x_i,y_i) |
---|
| 148 | self.beta_e = self.mask(x_e,y_e) |
---|
[774] | 149 | self.x_e,self.y_e = x_e,y_e |
---|
[771] | 150 | def mask0(self,c,c0,delta): # 1D building block for Davies relaxation |
---|
| 151 | N1, N2 = self.N1, self.N2 |
---|
| 152 | m = np.zeros(c.size) |
---|
[774] | 153 | c1,c2 = c0-N1*delta, c0-(N1+N2)*delta |
---|
[771] | 154 | for i in range(c.size): |
---|
| 155 | ci=c[i] |
---|
[774] | 156 | m[i] = (1.+np.cos((ci-c2)*np.pi/(N2*delta)))/2.0 |
---|
| 157 | if ci < c2 : m[i]=1. |
---|
| 158 | if ci > c1 : m[i]=0. |
---|
[771] | 159 | return m |
---|
| 160 | def relax(self, llm, step, flow): |
---|
| 161 | beta_i, beta_e = self.beta_i, self.beta_e |
---|
| 162 | m,S,u,Phi,W=flow |
---|
| 163 | for l in range(llm): |
---|
| 164 | step.mass[:,l] = beta_i*step.mass[:,l] + (1.-beta_i)*m[:,l] |
---|
| 165 | step.theta_rhodz[:,l] = beta_i*step.theta_rhodz[:,l] + (1.-beta_i)*S[:,l] |
---|
| 166 | step.u[:,l] = beta_e*step.u[:,l] + (1.-beta_e)*u[:,l] |
---|
| 167 | for l in range(llm+1): |
---|
| 168 | step.geopot[:,l] = beta_i*step.geopot[:,l] + (1.-beta_i)*Phi[:,l] |
---|
| 169 | step.W[:,l] = beta_i*step.W[:,l] + (1.-beta_i)*W[:,l] |
---|
| 170 | |
---|
| 171 | class myDavies(Davies): |
---|
| 172 | def mask(self,x,y): |
---|
[774] | 173 | logging.debug('davies dy = %f' % dy) |
---|
| 174 | return self.mask0(y,.5*Ly,dy)*self.mask0(-y,.5*Ly,dy) |
---|
| 175 | |
---|
| 176 | |
---|
[771] | 177 | |
---|
[763] | 178 | with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator |
---|
| 179 | comm = client.comm |
---|
| 180 | mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() |
---|
[761] | 181 | |
---|
[774] | 182 | #-------------Logging verbosity configuration--------------------------- |
---|
| 183 | myloglevel = args.log |
---|
| 184 | myloglevel = getattr(logging, myloglevel.upper()) |
---|
| 185 | logging.basicConfig(filename='out.log',filemode='w',level=myloglevel, |
---|
| 186 | format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' ) |
---|
| 187 | |
---|
| 188 | logging.debug('%d/%d starting'%(mpi_rank,mpi_size)) |
---|
| 189 | |
---|
[763] | 190 | g, Lx, Ly = 9.81, 4e7, 6e6 |
---|
[769] | 191 | nx, ny, llm = 200, 30, 25 |
---|
[774] | 192 | dx,dy = Lx/nx,Ly/ny |
---|
[761] | 193 | |
---|
[763] | 194 | unst.setvar('g',g) |
---|
[761] | 195 | |
---|
[763] | 196 | thermo, mesh, params, flow0, gas0 = baroclinic_3D(Lx,nx,Ly,ny,llm) |
---|
[761] | 197 | |
---|
[769] | 198 | mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) |
---|
[774] | 199 | print( 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk]) |
---|
[769] | 200 | unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) |
---|
[761] | 201 | |
---|
[774] | 202 | T = args.T |
---|
| 203 | dt = 360. |
---|
[769] | 204 | dz = flow0[3].max()/(params.g*llm) |
---|
| 205 | nt = int(math.ceil(T/dt)) |
---|
| 206 | dt = T/nt |
---|
[774] | 207 | logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)) |
---|
[769] | 208 | |
---|
| 209 | caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass |
---|
| 210 | |
---|
| 211 | if False: # time stepping in Python |
---|
| 212 | caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,thermo,params,params.g) |
---|
| 213 | scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) |
---|
| 214 | def next_flow(m,S,u,Phi,W): |
---|
| 215 | return scheme.advance((m,S,u,Phi,W),nt) |
---|
| 216 | |
---|
| 217 | else: # time stepping in Fortran |
---|
| 218 | scheme = time_step.ARK2(None, dt) |
---|
[771] | 219 | caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g) |
---|
| 220 | davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e) |
---|
[774] | 221 | print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() ) |
---|
| 222 | logging.debug('mask min/max :') |
---|
| 223 | logging.debug(davies.beta_i.min()) |
---|
| 224 | logging.debug(davies.beta_i.max()) |
---|
[769] | 225 | def next_flow(m,S,u,Phi,W): |
---|
| 226 | # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) |
---|
| 227 | caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u |
---|
| 228 | caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W |
---|
[771] | 229 | for i in range(nt): |
---|
| 230 | caldyn_step.next() |
---|
| 231 | davies.relax(llm, caldyn_step, flow0) |
---|
[769] | 232 | return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(), |
---|
[771] | 233 | caldyn_step.geopot.copy(), caldyn_step.W.copy()) |
---|
[769] | 234 | |
---|
[774] | 235 | m,S,u,Phi,W = flow0 |
---|
[769] | 236 | if caldyn_thermo == unst.thermo_theta: |
---|
| 237 | s=S/m |
---|
| 238 | theta = thermo.T0*np.exp(s/thermo.Cpd) |
---|
| 239 | S=m*theta |
---|
| 240 | title_format = 'Potential temperature at t=%g s (K)' |
---|
| 241 | else: |
---|
| 242 | title_format = 'Specific entropy at t=%g s (J/K/kg)' |
---|
| 243 | |
---|
| 244 | w=mesh.field_mass() |
---|
| 245 | z=mesh.field_mass() |
---|
| 246 | |
---|
| 247 | |
---|
| 248 | #T, nslice, dt = 3600., 1, 3600. |
---|
[774] | 249 | Nslice = args.N |
---|
[769] | 250 | |
---|
| 251 | with xios.Context_Curvilinear(mesh,1, 24*3600) as context: |
---|
[763] | 252 | # now XIOS knows about the mesh and we can write to disk |
---|
[769] | 253 | v = mesh.field_mass() # specific volume (diagnosed) |
---|
| 254 | |
---|
| 255 | for i in range(Nslice): |
---|
[763] | 256 | context.update_calendar(i) |
---|
[769] | 257 | |
---|
| 258 | # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W |
---|
| 259 | gas, w, z = diagnose(Phi,S,m,W) |
---|
| 260 | |
---|
| 261 | # write to disk |
---|
[774] | 262 | context.send_field_primal('ps', davies.beta_i) |
---|
[763] | 263 | context.send_field_primal('temp', gas0.T) |
---|
[769] | 264 | context.send_field_primal('p', gas.p) |
---|
| 265 | context.send_field_primal('theta', gas.s) |
---|
| 266 | context.send_field_primal('uz', w) |
---|
| 267 | |
---|
[774] | 268 | print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')) |
---|
[769] | 269 | |
---|
| 270 | time1, elapsed1 =time.time(), unst.getvar('elapsed') |
---|
| 271 | m,S,u,Phi,W = next_flow(m,S,u,Phi,W) |
---|
| 272 | time2, elapsed2 =time.time(), unst.getvar('elapsed') |
---|
| 273 | factor = 1000./nt |
---|
[774] | 274 | print( 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)) |
---|
[769] | 275 | factor = 1e9/(4*nt*nx*ny*llm) |
---|
[774] | 276 | print( 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)) |
---|
[769] | 277 | |
---|
| 278 | context.update_calendar(Nslice+1) # make sure XIOS writes last iteration |
---|
[774] | 279 | logging.info('************DONE************') |
---|