- Timestamp:
- 11/21/18 16:17:24 (6 years ago)
- Location:
- codes/icosagcm/devel/Python/test
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py
r780 r785 25 25 parser.add_argument("--T", type=float, default=3600., 26 26 help="Length of time slice in seconds") 27 parser.add_argument("--N", type=int, default= 24,27 parser.add_argument("--N", type=int, default=48, 28 28 help="Number of time slices") 29 29 parser.add_argument("--log",type=str,default='DEBUG',choices=('DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'), … … 141 141 gas = thermo.set_vs(v,s) 142 142 return gas, w, z 143 144 def grad(mesh, s_ik): 145 llm=s_ik.shape[1] 146 edge_num, left, right = mesh.edge_num, mesh.left, mesh.right 147 grad_ek = prec.zeros((edge_num, llm)) 148 for e in range(edge_num): 149 left_e = left[e]-1 # Fortran => Python indexing 150 right_e = right[e]-1 # Fortran => Python indexing 151 grad_ek[e,:] = s_ik[right_e,:]-s_ik[left_e,:] 152 return grad_ek 153 154 def div(mesh, u_ek): 155 llm=u_ek.shape[1] 156 primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge 157 primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai 158 div_ik = prec.zeros((primal_num, llm)) 159 for i in range(primal_num): 160 div_i = 0. 161 for iedge in range(primal_deg[i]): 162 e = primal_edge[i,iedge]-1 # Fortran => Python indexing 163 div_i = div_i + u_ek[e,:]*primal_ne[i,iedge]*le_de[e] 164 div_ik[i,:] = div_i / Ai[i] 165 return div_ik 166 167 def curl(mesh, u_ek): 168 llm=u_ek.shape[1] 169 dual_num, dual_deg, dual_edge, dual_ne, Av = mesh.dual_num, mesh.dual_deg, mesh.dual_edge, mesh.dual_ne, mesh.Av 170 curl_vk = prec.zeros((dual_num, llm)) 171 for v in range(dual_num): 172 curl_v = 0. 173 for iedge in range(dual_deg[v]): 174 e = dual_edge[v,iedge]-1 # Fortran => Python indexing 175 curl_v = curl_v + u_ek[e,:]*dual_ne[v,iedge] 176 curl_vk[v,:] = curl_v / Av[v] 177 return curl_vk 143 178 144 179 class Davies: … … 189 224 190 225 g, Lx, Ly = 9.81, 4e7, 6e6 191 nx, ny, llm = 200, 30, 25226 nx, ny, llm = 200, 30, 30 192 227 dx,dy = Lx/nx,Ly/ny 193 228 … … 227 262 caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u 228 263 caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W 264 kappa = dx**4/43200. 229 265 for i in range(nt): 230 266 caldyn_step.next() 231 267 davies.relax(llm, caldyn_step, flow0) 268 m,S = caldyn_step.mass, caldyn_step.theta_rhodz 269 s = S/m 270 laps = mesh.field_mass() 271 bilaps = mesh.field_mass() 272 unst.ker.dynamico_scalar_laplacian(s,laps) 273 unst.ker.dynamico_scalar_laplacian(laps,bilaps) 274 # grads = grad(mesh,s) 275 # laps = div(mesh,grads) 276 # gradlaps = grad(mesh,laps) # bilaplacian 277 # bilaps = div(mesh,gradlaps) 278 caldyn_step.theta_rhodz[:] = S - dt*kappa*bilaps*m # Euler step 232 279 return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(), 233 280 caldyn_step.geopot.copy(), caldyn_step.W.copy()) … … 254 301 255 302 for i in range(Nslice): 256 context.update_calendar(i )303 context.update_calendar(i+1) 257 304 258 305 # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W 259 306 gas, w, z = diagnose(Phi,S,m,W) 260 307 curl_vk = curl(mesh, u) 261 308 # write to disk 262 309 context.send_field_primal('ps', davies.beta_i) 263 context.send_field_primal('temp', gas 0.T)310 context.send_field_primal('temp', gas.T) 264 311 context.send_field_primal('p', gas.p) 265 312 context.send_field_primal('theta', gas.s) 266 313 context.send_field_primal('uz', w) 267 314 context.send_field_primal('z', z) 315 context.send_field_primal('curl', curl_vk) 268 316 269 317 print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')) … … 277 325 print( 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)) 278 326 279 context.update_calendar(Nslice+1) # make sure XIOS writes last iteration280 327 logging.info('************DONE************') -
codes/icosagcm/devel/Python/test/xml/field_def_simple.xml
r782 r785 23 23 <field id="z" /> 24 24 <field id="theta" /> 25 <field id="curl" /> 25 26 </field_group> 26 27 -
codes/icosagcm/devel/Python/test/xml/filedef_simple.xml
r782 r785 20 20 <field field_ref="z" /> 21 21 <field field_ref="theta" /> 22 <field field_ref="curl" /> 22 23 <field field_ref="phis" operation="once" freq_offset="0ts" name="PHIS" standard_name="surface_geopotential" long_name="Surface geopotential" unit="m2/m2"/> 23 24 </field_group>
Note: See TracChangeset
for help on using the changeset viewer.