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