1 | MODULE data_unstructured_mod |
---|
2 | USE ISO_C_BINDING |
---|
3 | USE OMP_LIB |
---|
4 | IMPLICIT NONE |
---|
5 | SAVE |
---|
6 | |
---|
7 | #include "unstructured.h90" |
---|
8 | |
---|
9 | INTEGER, PARAMETER :: eta_mass=1, eta_lag=2, & |
---|
10 | thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, & |
---|
11 | caldyn_vert_cons=1, max_nb_stage=5 |
---|
12 | INDEX, BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, & |
---|
13 | caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0 |
---|
14 | LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_hevi_solver=.TRUE. |
---|
15 | |
---|
16 | INDEX, BIND(C) :: llm, nqdyn, edge_num, primal_num, dual_num, & |
---|
17 | max_primal_deg, max_dual_deg, max_trisk_deg |
---|
18 | INDEX, ALLOCATABLE :: & ! deg(ij) = nb of vertices = nb of edges of primal/dual cell ij |
---|
19 | primal_deg(:), primal_edge(:,:), primal_vertex(:,:), primal_ne(:,:), & |
---|
20 | dual_deg(:), dual_edge(:,:), dual_vertex(:,:), dual_ne(:,:), & |
---|
21 | trisk_deg(:), trisk(:,:), & |
---|
22 | left(:), right(:), up(:), down(:) |
---|
23 | ! left and right are adjacent primal cells |
---|
24 | ! flux is positive when going from left to right |
---|
25 | ! up and down are adjacent dual cells |
---|
26 | ! circulation is positive when going from down to up |
---|
27 | |
---|
28 | TIME, BIND(C) :: elapsed |
---|
29 | NUM, BIND(C) :: g, ptop, cpp, cppv, Rd, Rv, preff, Treff, pbot, Phi_bot, rho_bot |
---|
30 | NUM :: kappa |
---|
31 | NUM1(max_nb_stage), BIND(C) :: tauj ! diagonal of fast Butcher tableau |
---|
32 | NUM2(max_nb_stage,max_nb_stage), BIND(C) :: cslj, cflj ! slow and fast modified Butcher tableaus |
---|
33 | NUM1(:), ALLOCATABLE :: le_de, fv, Av, Ai |
---|
34 | NUM2(:,:), ALLOCATABLE :: Riv2, wee, ap,bp, mass_bl, mass_dak, mass_dbk |
---|
35 | |
---|
36 | INTEGER(C_INT), BIND(C) :: comm_icosa, dynamico_mpi_rank=0 |
---|
37 | |
---|
38 | INTEGER, PARAMETER :: id_dev1=1, id_dev2=2, & |
---|
39 | id_pvort_only=3, id_slow_hydro=4, id_fast=5, id_coriolis=6, id_theta=7, id_geopot=8, id_vert=9, & |
---|
40 | id_solver=10, id_slow_NH=11, id_NH_geopot=12, id_vert_NH=13, id_update=14, nb_routines=14 |
---|
41 | TIME, PRIVATE :: start_time, time_spent(nb_routines) ! time spent in each kernel |
---|
42 | INTEGER, PRIVATE :: current_id, nb_calls(nb_routines) |
---|
43 | INTEGER(KIND=8), PRIVATE :: bytes(nb_routines) ! bytes read or written by each kernel |
---|
44 | CHARACTER(len = 10) :: id_name(nb_routines) = & |
---|
45 | (/'dev1 ', 'dev2 ', & |
---|
46 | 'pvort_only', 'slow_hydro', 'fast ', 'coriolis ', 'theta ', 'geopot ', 'vert ', & |
---|
47 | 'solver ', 'slow_NH ', 'NH_geopot ', 'vert_NH ', 'update ' /) |
---|
48 | |
---|
49 | INTEGER, PARAMETER ::transfer_primal=1, transfer_edge=2, transfer_dual=3, transfer_max=3 |
---|
50 | TYPE Halo_transfer |
---|
51 | INTEGER :: ranks ! size of arrays rank, len |
---|
52 | INTEGER, ALLOCATABLE :: rank(:), & ! MPI ranks to communicate with |
---|
53 | num(:), & ! number of cells to send to / receive from other MPI ranks |
---|
54 | cells(:) ! local indices of cells to send/receive |
---|
55 | NUM, ALLOCATABLE :: buf2(:,:) |
---|
56 | END TYPE Halo_transfer |
---|
57 | TYPE(Halo_transfer), TARGET :: send_info(transfer_max), recv_info(transfer_max) |
---|
58 | |
---|
59 | CONTAINS |
---|
60 | |
---|
61 | !---------------------------- PROFILING -------------------------- |
---|
62 | |
---|
63 | SUBROUTINE init_trace() |
---|
64 | !$OMP MASTER |
---|
65 | time_spent(:)=0. |
---|
66 | bytes(:)=0 |
---|
67 | nb_calls(:)=0 |
---|
68 | !$OMP END MASTER |
---|
69 | END SUBROUTINE init_trace |
---|
70 | |
---|
71 | SUBROUTINE print_trace() |
---|
72 | INTEGER :: id |
---|
73 | TIME :: total_spent |
---|
74 | !$OMP MASTER |
---|
75 | IF(dynamico_mpi_rank==0) THEN |
---|
76 | total_spent=SUM(time_spent) |
---|
77 | IF(total_spent>1.) THEN |
---|
78 | PRINT *, '========================= Performance metrics =========================' |
---|
79 | PRINT *, 'Total time spent in instrumented code (seconds) :', total_spent |
---|
80 | PRINT *, 'Name, #calls, %time, microsec/call, MB/sec' |
---|
81 | DO id=1,nb_routines |
---|
82 | IF(nb_calls(id)>0) PRINT *, id_name(id), nb_calls(id), INT(100.*time_spent(id)/total_spent), & |
---|
83 | INT(1e6*time_spent(id)/nb_calls(id)), INT(1e-6*bytes(id)/time_spent(id)) |
---|
84 | END DO |
---|
85 | CALL init_trace() |
---|
86 | END IF |
---|
87 | END IF |
---|
88 | !$OMP END MASTER |
---|
89 | END SUBROUTINE print_trace |
---|
90 | |
---|
91 | SUBROUTINE enter_trace(id, nbytes) |
---|
92 | INTEGER :: id, nbytes |
---|
93 | !$OMP MASTER |
---|
94 | current_id = id |
---|
95 | bytes(id) = bytes(id) + nbytes |
---|
96 | nb_calls(id)=nb_calls(id)+1 |
---|
97 | start_time = OMP_GET_WTIME() |
---|
98 | !$OMP END MASTER |
---|
99 | END SUBROUTINE enter_trace |
---|
100 | |
---|
101 | SUBROUTINE exit_trace() |
---|
102 | TIME :: elapsed |
---|
103 | !$OMP MASTER |
---|
104 | elapsed = OMP_GET_WTIME()-start_time |
---|
105 | IF(elapsed<0.) elapsed=0. |
---|
106 | time_spent(current_id) = time_spent(current_id) + elapsed |
---|
107 | !$OMP END MASTER |
---|
108 | END SUBROUTINE exit_trace |
---|
109 | |
---|
110 | !---------------------------- CONTEXT INITIALIZATION -------------------------- |
---|
111 | |
---|
112 | #define ALLOC1(v,n1) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1)) |
---|
113 | #define ALLOC2(v,n1,n2) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2)) |
---|
114 | |
---|
115 | SUBROUTINE init_mesh( & |
---|
116 | primal_deg_, primal_edge_, primal_ne_, & |
---|
117 | dual_deg_, dual_edge_, dual_ne_, dual_vertex_, & |
---|
118 | left_, right_, up_, down_ ,& |
---|
119 | trisk_deg_, trisk_) BINDC(init_mesh) |
---|
120 | INDEX :: primal_deg_(primal_num), primal_edge_(max_primal_deg,primal_num), & |
---|
121 | primal_ne_(max_primal_deg,primal_num), & |
---|
122 | dual_deg_(dual_num), dual_edge_(max_dual_deg,dual_num), & |
---|
123 | dual_ne_(max_dual_deg,dual_num), & |
---|
124 | dual_vertex_(max_dual_deg,dual_num), & |
---|
125 | trisk_deg_(edge_num), trisk_(max_trisk_deg, edge_num) |
---|
126 | INDEX, DIMENSION(edge_num) :: left_, right_, down_, up_ |
---|
127 | |
---|
128 | PRINT *, 'init_mesh ...' |
---|
129 | PRINT *, 'Primal mesh : ', primal_num, max_primal_deg |
---|
130 | PRINT *, 'Dual mesh : ', dual_num, max_dual_deg |
---|
131 | PRINT *, 'Edge mesh : ', edge_num, max_trisk_deg |
---|
132 | PRINT *, 'Vertical levels :', llm |
---|
133 | ALLOC1(primal_deg, primal_num) |
---|
134 | ALLOC2(primal_edge, max_primal_deg,primal_num) |
---|
135 | ALLOC2(primal_ne, max_primal_deg,primal_num) |
---|
136 | ALLOC1(dual_deg,dual_num) |
---|
137 | ALLOC2(dual_edge, max_dual_deg,dual_num) |
---|
138 | ALLOC2(dual_ne, max_dual_deg,dual_num) |
---|
139 | ALLOC2(dual_vertex, max_dual_deg,dual_num) |
---|
140 | ALLOC1(trisk_deg, edge_num) |
---|
141 | ALLOC2(trisk, max_trisk_deg, edge_num) |
---|
142 | PRINT *, SHAPE(trisk), edge_num |
---|
143 | ALLOC1(left, edge_num) |
---|
144 | ALLOC1(right, edge_num) |
---|
145 | ALLOC1(up, edge_num) |
---|
146 | ALLOC1(down, edge_num) |
---|
147 | primal_deg(:) = primal_deg_(:) |
---|
148 | primal_edge(:,:) = primal_edge_(:,:) |
---|
149 | primal_ne(:,:) = primal_ne_(:,:) |
---|
150 | dual_deg(:) = dual_deg_(:) |
---|
151 | dual_edge(:,:) = dual_edge_(:,:) |
---|
152 | dual_ne(:,:) = dual_ne_(:,:) |
---|
153 | dual_vertex(:,:) = dual_vertex_(:,:) |
---|
154 | IF(MINVAL(dual_deg)<3) THEN |
---|
155 | STOP 'At least one dual cell has less than 3 vertices' |
---|
156 | END IF |
---|
157 | IF(MINVAL(primal_deg)<3) THEN |
---|
158 | STOP 'At least one primal cell has less than 3 vertices' |
---|
159 | END IF |
---|
160 | left(:)=left_(:) |
---|
161 | right(:)=right_(:) |
---|
162 | down(:)=down_(:) |
---|
163 | up=up_(:) |
---|
164 | trisk_deg(:)=trisk_deg_(:) |
---|
165 | trisk(:,:)=trisk_(:,:) |
---|
166 | PRINT *, MAXVAL(primal_edge), edge_num |
---|
167 | PRINT *, MAXVAL(dual_edge), edge_num |
---|
168 | PRINT *, MAXVAL(dual_vertex), dual_num |
---|
169 | PRINT *, MAXVAL(trisk), edge_num |
---|
170 | PRINT *, MAX(MAXVAL(left),MAXVAL(right)), primal_num |
---|
171 | PRINT *, MAX(MAXVAL(up),MAXVAL(down)), dual_num |
---|
172 | PRINT *, SHAPE(trisk), edge_num |
---|
173 | PRINT *,' ... Done.' |
---|
174 | END SUBROUTINE init_mesh |
---|
175 | |
---|
176 | ! Input arrays to init_metric and init_hybrid are declared DBL |
---|
177 | ! => always float64 on the Python side |
---|
178 | ! They are copied to Fortran arrays of type NUM (float or double) |
---|
179 | |
---|
180 | SUBROUTINE init_metric(Ai_, Av_, fv_, le_de_, Riv2_, wee_) BINDC(init_metric) |
---|
181 | DBL :: Ai_(primal_num), Av_(dual_num), fv_(dual_num), le_de_(edge_num), & |
---|
182 | Riv2_(max_dual_deg,dual_num), wee_(max_trisk_deg,edge_num) |
---|
183 | PRINT *, 'init_metric ...' |
---|
184 | ALLOC1(Ai,primal_num) |
---|
185 | ALLOC1(Av,dual_num) |
---|
186 | ALLOC1(fv,dual_num) |
---|
187 | ALLOC1(le_de,edge_num) |
---|
188 | ALLOC2(Riv2, max_dual_deg, dual_num) |
---|
189 | ALLOC2(wee, max_trisk_deg, edge_num) |
---|
190 | Ai(:) = Ai_(:) |
---|
191 | Av(:) = Av_(:) |
---|
192 | fv(:) = fv_(:) |
---|
193 | le_de(:) = le_de_(:) |
---|
194 | Riv2(:,:)=Riv2_(:,:) |
---|
195 | wee(:,:) = wee_(:,:) |
---|
196 | PRINT *, MAXVAL(ABS(Ai)) |
---|
197 | PRINT *, MAXVAL(ABS(Av)) |
---|
198 | PRINT *, MAXVAL(ABS(fv)) |
---|
199 | PRINT *, MAXVAL(ABS(le_de)) |
---|
200 | PRINT *, MAXVAL(ABS(Riv2)) |
---|
201 | PRINT *, MAXVAL(ABS(wee)) |
---|
202 | PRINT *, MINVAL(right), MAXVAL(right) |
---|
203 | PRINT *, MINVAL(right), MAXVAL(left) |
---|
204 | PRINT *,' ... Done.' |
---|
205 | IF(nb_threads==0) nb_threads=OMP_GET_MAX_THREADS() |
---|
206 | PRINT *,'OpenMP : max_threads, num_procs, nb_threads', OMP_GET_MAX_THREADS(), OMP_GET_NUM_PROCS(), nb_threads |
---|
207 | END SUBROUTINE init_metric |
---|
208 | ! |
---|
209 | SUBROUTINE show_openmp() BINDC(show_openmp) |
---|
210 | PRINT *,'OpenMP : max_threads, num_procs', OMP_GET_MAX_THREADS(), OMP_GET_NUM_PROCS() |
---|
211 | END SUBROUTINE show_openmp |
---|
212 | ! |
---|
213 | SUBROUTINE init_params() BINDC(init_params) |
---|
214 | PRINT *, 'Setting physical parameters ...' |
---|
215 | IF(hydrostatic) THEN |
---|
216 | PRINT *, 'Hydrostatic dynamics (HPE)' |
---|
217 | ELSE |
---|
218 | PRINT *, 'Non-hydrostatic dynamics (Euler)' |
---|
219 | END IF |
---|
220 | kappa = Rd/cpp |
---|
221 | PRINT *, 'g = ',g |
---|
222 | PRINT *, 'preff = ',preff |
---|
223 | PRINT *, 'Treff = ',Treff |
---|
224 | PRINT *, 'Rd = ',Rd |
---|
225 | PRINT *, 'cpp = ',cpp |
---|
226 | PRINT *, 'kappa = ',kappa |
---|
227 | PRINT *, '... Done' |
---|
228 | CALL init_trace |
---|
229 | END SUBROUTINE init_params |
---|
230 | ! |
---|
231 | SUBROUTINE init_hybrid(bl,dak,dbk) BINDC(init_hybrid) |
---|
232 | DBL :: bl(llm+1, primal_num), & |
---|
233 | dak(llm, primal_num), dbk(llm, primal_num) |
---|
234 | PRINT *, 'Setting hybrid coefficients ...' |
---|
235 | ALLOC2(mass_bl, llm+1, primal_num) |
---|
236 | ALLOC2(mass_dak, llm, primal_num) |
---|
237 | ALLOC2(mass_dbk, llm, primal_num) |
---|
238 | mass_bl(:,:) = bl(:,:) |
---|
239 | mass_dak(:,:) = dak(:,:) |
---|
240 | mass_dbk(:,:) = dbk(:,:) |
---|
241 | PRINT *, '... Done, llm = ', llm |
---|
242 | END SUBROUTINE Init_hybrid |
---|
243 | |
---|
244 | END MODULE data_unstructured_mod |
---|