source: codes/icosagcm/devel/src/unstructured/data_unstructured.F90 @ 813

Last change on this file since 813 was 813, checked in by jisesh, 5 years ago

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