[663] | 1 | #print 'import dynamico' |
---|
| 2 | #import dynamico |
---|
[622] | 3 | print 'import dynamico.unstructured' |
---|
| 4 | import dynamico.unstructured as unst |
---|
| 5 | print 'import dynamico.xios' |
---|
| 6 | import dynamico.xios as xios |
---|
| 7 | print 'Done.' |
---|
[631] | 8 | from dynamico.meshes import MPAS_Mesh as Mesh, radian |
---|
[622] | 9 | |
---|
[663] | 10 | from mpi4py import MPI |
---|
| 11 | comm = MPI.COMM_WORLD |
---|
| 12 | mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size() |
---|
| 13 | print '%d/%d starting'%(mpi_rank,mpi_size) |
---|
| 14 | |
---|
| 15 | import numpy as np |
---|
| 16 | |
---|
[622] | 17 | def boundaries(degree,points,lon,lat): |
---|
| 18 | n, nvertex = len(degree), degree.max() |
---|
| 19 | bnd_lon, bnd_lat = np.zeros((n,nvertex)),np.zeros((n,nvertex)) |
---|
| 20 | for ij in range(n): |
---|
| 21 | nb=degree[ij] |
---|
| 22 | for k in range(nb): |
---|
| 23 | vertex = points[ij,k] |
---|
| 24 | bnd_lon[ij,k]=lon[vertex] |
---|
| 25 | bnd_lat[ij,k]=lat[vertex] |
---|
| 26 | for k in range(nb,nvertex): # repeat last vertex if nb<nvertex |
---|
| 27 | bnd_lon[ij,k]=bnd_lon[ij,nb-1] |
---|
| 28 | bnd_lat[ij,k]=bnd_lat[ij,nb-1] |
---|
| 29 | return bnd_lon, bnd_lat |
---|
| 30 | |
---|
| 31 | #----------------------------- read MPAS mesh -------------------------------- |
---|
| 32 | |
---|
| 33 | grid, llm, nqdyn, nqtot = 10242, 1,1,1 # 2562, 10242, 40962 |
---|
| 34 | Omega, radius, g, gh0 = 2.*np.pi/86400., 6.4e6, 1., 2.94e4 |
---|
| 35 | N, T, courant = 40, 10800., 1.2 # simulation length = N*T |
---|
| 36 | print 'Omega, planetary PV', Omega, 2*Omega/gh0 |
---|
| 37 | |
---|
| 38 | def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter |
---|
| 39 | print 'Reading MPAS mesh ...' |
---|
[631] | 40 | mesh = Mesh('grids/x1.%d.grid.nc'%grid, llm, nqdyn, radius, f) |
---|
[622] | 41 | print '...Done' |
---|
| 42 | |
---|
| 43 | lon_i, lat_i, lon_v, lat_v = [x*radian for x in mesh.lon_i, mesh.lat_i, mesh.lon_v, mesh.lat_v] |
---|
| 44 | #lon_i, lat_i, lon_v, lat_v = mesh.lon_i, mesh.lat_i, mesh.lon_v, mesh.lat_v |
---|
| 45 | |
---|
| 46 | #--------------------------------- initialize XIOS -------------------------- |
---|
| 47 | |
---|
| 48 | def setup_domain(cat,id, degree,vertex,lon,lat,lon_v,lat_v): # cat in 'domain','domaingroup' |
---|
| 49 | bnd_lon, bnd_lat = boundaries(degree, vertex, lon_v, lat_v) |
---|
| 50 | ncell_tot, nvertex = bnd_lon.shape |
---|
| 51 | mesh = xios.Handle(cat,id) |
---|
| 52 | mesh.set_attr(type='unstructured') |
---|
| 53 | mesh.set_attr(ni_glo=ncell_tot, ibegin=0, ni=ncell_tot) |
---|
| 54 | mesh.set_attr(i_index=np.arange(ncell_tot,dtype=np.int32)) |
---|
| 55 | mesh.set_attr(lonvalue_1d=lon, latvalue_1d=lat) |
---|
| 56 | mesh.set_attr(nvertex=nvertex,bounds_lon_1d=bnd_lon, bounds_lat_1d=bnd_lat) |
---|
| 57 | |
---|
| 58 | unst.ker.dynamico_setup_xios() |
---|
| 59 | |
---|
| 60 | if True: |
---|
| 61 | lev, levp1, nq = [xios.Handle('axis',name) for name in 'lev', 'levp1', 'nq'] |
---|
| 62 | lev.set_attr(n_glo=llm, value=np.arange(llm)) |
---|
| 63 | levp1.set_attr(n_glo=llm+1, value=np.arange(llm+1)) |
---|
| 64 | nq.set_attr(n_glo=nqtot, value=np.arange(nqtot)) |
---|
| 65 | mesh_i=setup_domain('domaingroup','i', mesh.primal_deg, mesh.primal_vertex-1, lon_i, lat_i, lon_v, lat_v) # primal mesh |
---|
| 66 | # mesh_v=setup_domain('domain','v', mesh.dual_deg, mesh.dual_vertex-1, lon_v, lat_v, lon_i, lat_i) # dual mesh |
---|
| 67 | |
---|
| 68 | # CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ, ni=ncell) |
---|
| 69 | # CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo) |
---|
| 70 | # CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat) |
---|
| 71 | #mesh_e = xios.Handle('domain','u') |
---|
| 72 | |
---|
| 73 | # works up to that point ! |
---|
| 74 | |
---|
| 75 | dtime = xios.Duration(second=3600.) |
---|
| 76 | |
---|
| 77 | # CALL xios_set_timestep(dtime) |
---|
| 78 | if True: |
---|
| 79 | calendar=xios.Handle('calendar_wrapper') |
---|
| 80 | print (dtime.year, dtime.month, dtime.day, dtime.hour, dtime.second) |
---|
| 81 | calendar.set_attr(timestep=dtime) |
---|
| 82 | calendar.update_timestep() |
---|
| 83 | else: |
---|
| 84 | unst.ker.dynamico_xios_set_timestep(3600.) |
---|
| 85 | |
---|
| 86 | # CALL xios_set_fieldgroup_attr("standard_output", freq_op=itau_out*xios_timestep, freq_offset=(itau_out-1)*xios_timestep) |
---|
| 87 | if False: |
---|
| 88 | std=xios.Handle('fieldgroup','standard_output') |
---|
| 89 | std.set_attr(freq_op=dtime) |
---|
| 90 | freq_offset = xios.Duration(second=0) |
---|
| 91 | std.set_attr(freq_offset=freq_offset) |
---|
| 92 | |
---|
| 93 | if True: |
---|
| 94 | print 'xios.context_close_definition()' |
---|
| 95 | xios.context_close_definition() |
---|
| 96 | |
---|
| 97 | no_error=True |
---|
| 98 | try: |
---|
| 99 | for i in range(100): |
---|
| 100 | if True: |
---|
| 101 | unst.ker.dynamico_xios_update_calendar(i) |
---|
| 102 | else: |
---|
| 103 | xios.update_calendar(i) |
---|
| 104 | print 'send_field', i |
---|
| 105 | xios.send_field('ps', lat_i) |
---|
| 106 | except: |
---|
| 107 | print 'An error has occured, trying to close XIOS properly' |
---|
| 108 | no_error=False |
---|
| 109 | # CALL xios_recv_field(name,field) |
---|
| 110 | # CALL xios_send_field(name,field_tmp) |
---|
| 111 | |
---|
| 112 | print 'xios.context_finalize()' |
---|
| 113 | xios.context_finalize() |
---|
| 114 | |
---|
| 115 | print 'xios.finalize()' |
---|
| 116 | xios.finalize() |
---|
| 117 | |
---|
| 118 | print 'Done' |
---|
| 119 | |
---|
| 120 | if not no_error : raise |
---|