EMPIRE DA  v1.9.1
Data assimilation codes using EMPIRE communication
 All Classes Files Functions Variables Pages
minimal_model_v3.f90
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! minimal_model_comms just sets up communicators for empire for
3 ! testing purposes
4 !
5 !The MIT License (MIT)
6 !
7 !Copyright (c) 2015 Philip A. Browne
8 !
9 !Permission is hereby granted, free of charge, to any person obtaining a copy
10 !of this software and associated documentation files (the "Software"), to deal
11 !in the Software without restriction, including without limitation the rights
12 !to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13 !copies of the Software, and to permit persons to whom the Software is
14 !furnished to do so, subject to the following conditions:
15 !
16 !The above copyright notice and this permission notice shall be included in all
17 !copies or substantial portions of the Software.
18 !
19 !THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20 !IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 !FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22 !AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 !LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24 !OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25 !SOFTWARE.
26 !
27 !Email: p.browne@reading.ac.uk
28 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
29 
31  implicit none
32 
33  include 'mpif.h'
34  integer :: mpi_err
35  integer :: world_rank
36  integer :: world_size
37  integer, parameter :: mdl_num_proc=5
38  integer :: mdl_mpi_comm
39  integer :: temp_mdls_comm
40  integer :: temp_mdls_rank
41  integer :: mdl_rank
42  integer :: da
43  integer, parameter :: rk = kind(1.0d0)
44  real(kind=rk), allocatable, dimension(:) :: state_vector
45  integer :: state_dim
46  integer :: i
47  integer :: cpl_root
48  integer :: particle_rank
49  integer :: nda
50  integer :: nens
51  integer :: temp_mdls_size
52  integer :: temp_cpl_comm
53  integer :: first_ptcl
54  integer :: final_ptcl
55  integer :: cpl_mpi_comm
56  integer :: null_mpi_comm
57  integer :: total_timesteps
58  integer :: tmp_colour_2
59  integer :: tmp_cpl_comm2
60  integer :: tmp_cpl_rank
61  integer :: status(mpi_status_size)
62 
63  print*,'RUNNING MINIMAL_MODEL_V3'
64 
65  call mpi_init(mpi_err)
66  if(mpi_err .eq. 0) then
67  print*,'mpi_init successful'
68  else
69  print*,'mpi_init unsuccessful'
70  end if
71 
72 
73  call mpi_comm_rank(mpi_comm_world,world_rank,mpi_err)
74  if(mpi_err .eq. 0) then
75  print*,'mpi_comm_rank successful'
76  print*,'world_rank = ',world_rank
77  else
78  print*,'mpi_comm_rank unsuccessful'
79  end if
80 
81 
82  call mpi_comm_size(mpi_comm_world,world_size,mpi_err)
83  if(mpi_err .eq. 0) then
84  print*,'mpi_comm_size successful'
85  print*,'world_size = ',world_size
86  else
87  print*,'mpi_comm_size unsuccessful'
88  end if
89 
90 
91  cpl_root = world_size-1
92  print*,'rank = ',world_rank,' on mpi_comm_world which has size '&
93  &,world_size
94  da = 0
95 
96 
97  call mpi_allreduce(mdl_num_proc,i,1,mpi_integer,mpi_max&
98  &,mpi_comm_world,mpi_err)
99  if(mpi_err .eq. 0) then
100  print*,'mpi_allreduce successful'
101  print*,'i = ',i
102  else
103  print*,'mpi_allreduce unsuccessful'
104  end if
105 !!!==================================================!!!
106 !!!================ BRUCE LEE =======================!!!
107 !!!==================================================!!!
108 
109 
110  call mpi_comm_split(mpi_comm_world,da,world_rank,temp_mdls_comm&
111  &,mpi_err)
112  if(mpi_err .eq. 0) then
113  print*,'mpi_comm_split successful: temp_mdls_comm created'
114  else
115  print*,'mpi_comm_split unsuccessful: temp_mdls_comm not created'
116  end if
117 
118  call mpi_comm_size(temp_mdls_comm,temp_mdls_size,mpi_err)
119  if(mpi_err .eq. 0) then
120  print*,'mpi_comm_size successful'
121  print*,'temp_mdls_size = ',temp_mdls_size
122  else
123  print*,'mpi_comm_size unsuccessful'
124  end if
125 
126 
127  if(mod(temp_mdls_size,mdl_num_proc) .ne. 0) then
128  print*,'MINIMAL MODEL LAUNCH ERROR.'
129  print*,'MUST LAUNCH A MULTIPLE OF ',mdl_num_proc,' copies of the &
130  &model'
131  stop
132  end if
133 
134 
135  nda = (world_size-temp_mdls_size)/mdl_num_proc
136  if(nda .lt. 1) then
137  print*,'MINIMAL MODEL COMMS v3 ERROR: nda is less than 1.'
138  print*,'Make sure you launch with a DA CODE'
139  stop
140  end if
141 
142 !!!==================================================!!!
143 !!!================ JEAN CLAUDE VAN DAMME ===========!!!
144 !!!==================================================!!!
145 
146 
147  nens = temp_mdls_size/mdl_num_proc
148  call mpi_comm_rank(temp_mdls_comm,temp_mdls_rank,mpi_err)
149  if(mpi_err .eq. 0) then
150  print*,'mpi_comm_rank successful'
151  print*,'temp_mdls_rank = ',temp_mdls_rank
152  else
153  print*,'mpi_comm_rank unsuccessful'
154  end if
155 
156 
157  particle_rank = temp_mdls_rank/mdl_num_proc
158 
159  call mpi_comm_split(temp_mdls_comm,particle_rank,temp_mdls_rank&
160  &,mdl_mpi_comm,mpi_err)
161  if(mpi_err .eq. 0) then
162  print*,'mpi_comm_split successful: mdl_mpi_comm created'
163  else
164  print*,'mpi_comm_split unsuccessful: mdl_mpi_comm not created'
165  end if
166 
167 
168 
169  call mpi_comm_rank(mdl_mpi_comm,mdl_rank,mpi_err)
170  if(mpi_err .eq. 0) then
171  print*,'mpi_comm_rank successful'
172  print*,'mdl_rank = ',mdl_rank
173  else
174  print*,'mpi_comm_rank unsuccessful'
175  end if
176 
177 !!!==================================================!!!
178 !!!============== STEVEN SEGAL ======================!!!
179 !!!==================================================!!!
180 
181 
182 
183  cpl_root = nda*particle_rank/nens
184  print*,'cpl_root = ',cpl_root
185 
186  if(cpl_root .lt. 0) then
187  print*,'MINIMAL MODEL COMMS v3 ERROR: cpl_root is less than 0.'
188  print*,'Make sure you launch with a DA CODE'
189  stop
190  end if
191 
192  call mpi_comm_split(mpi_comm_world,cpl_root,temp_mdls_rank&
193  &,temp_cpl_comm,mpi_err)
194  if(mpi_err .eq. 0) then
195  print*,'mpi_comm_split successful: temp_cpl_comm created'
196  else
197  print*,'mpi_comm_split unsuccessful: temp_cpl_comm not created'
198  end if
199 
200 !!!==================================================!!!
201 !!!============== CHUCK NORRIS ======================!!!
202 !!!==================================================!!!
203 
204  call mpi_comm_rank(temp_cpl_comm,tmp_cpl_rank,mpi_err)
205  if(mpi_err .eq. 0) then
206  print*,'mpi_comm_rank successful: tmp_cpl_rank = ',tmp_cpl_rank
207  else
208  print*,'mpl_comm_rank unsuccessful: tmp_cpl_rank not detected'
209  end if
210 
211  tmp_colour_2 = mod(tmp_cpl_rank,mdl_num_proc)
212  print*,'tmp_colour_2 = ',tmp_colour_2
213  call mpi_comm_split(temp_cpl_comm,tmp_colour_2,tmp_cpl_rank&
214  &,tmp_cpl_comm2,mpi_err)
215  if(mpi_err .eq. 0) then
216  print*,'mpi_comm_split successful: tmp_cpl_comm2 created'
217  else
218  print*,'mpi_comm_split unsuccessful: tmp_cpl_comm2 not created'
219  end if
220 
221 !!!==================================================!!!
222 !!!=============== JACKIE CHAN ======================!!!
223 !!!==================================================!!!
224 
225 
226  first_ptcl = -((-cpl_root)*nens/nda)
227  final_ptcl = -((-cpl_root-1)*nens/nda)-1
228 
229  first_ptcl = ceiling(real(cpl_root)*real(nens)/real(nda))
230  final_ptcl = ceiling(real(cpl_root+1)*real(nens)/real(nda))-1
231 
232 
233  print*,'range of particles = ',first_ptcl,final_ptcl
234 
235 
236 
237  do i = first_ptcl,final_ptcl
238  print*,'i = ',i,' particle_rank = ',particle_rank
239  if(i .eq. particle_rank) then
240  call mpi_comm_split(tmp_cpl_comm2,1,temp_mdls_rank&
241  &,cpl_mpi_comm,mpi_err)
242  print*,'created cpl_mpi_comm',cpl_mpi_comm
243  else
244  print*,'doing null splitting'
245  call mpi_comm_split(tmp_cpl_comm2,0,temp_mdls_rank&
246  &,null_mpi_comm,mpi_err)
247  print*,'created mpi_comm_null'
248  call mpi_comm_free(null_mpi_comm,mpi_err)
249  print*,'freed up null_mpi_comm'
250  end if
251 
252 
253  end do
254 
255  cpl_root = 1
256 
257 
258 !!!==================================================!!!
259 !!!=============== MICHELLE YEOH ====================!!!
260 !!!==================================================!!!
261 
262  !free up the temporary communicators
263  call mpi_comm_free(temp_cpl_comm ,mpi_err)
264  call mpi_comm_free(tmp_cpl_comm2,mpi_err)
265 
266 
267 
268  select case(mdl_rank)
269  case(0)
270  state_dim = 1
271  case(1)
272  state_dim = 3
273  case(2)
274  state_dim = 2
275  case(3)
276  state_dim = 5
277  case(4)
278  state_dim = 1
279  case default
280  print*,'it was at this point, model realised, he fucked up'
281  stop
282  end select
283  allocate(state_vector(state_dim))
284 
285 
286 
287 
288  state_vector = 10*mdl_rank + (/ (real(i,rk), i = 1,state_dim) /)
289 
290  print*,'state_vector = '
291  print*,state_vector
292 
293  print*,'doing a send on cpl_mpi_comm'
294 
295  call mpi_send(state_dim,1,mpi_integer,cpl_root,1,cpl_mpi_comm&
296  &,mpi_err)
297  print*,'finished the send on cpl_mpi_comm'
298 
299 
300 
301 
302 
303 
304 
305 
306 
307 
308 
309 
310 
311 
312 
313 
314 
315 
316 
317  print*,'Reading total_timesteps from file timesteps: '
318  open(11,file='timesteps',action='read',status='old')
319  read(11,*) total_timesteps
320  close(11)
321 
322  !send the state to da code with mpi_send
323  call mpi_send(state_vector,state_dim,mpi_double_precision,1,1&
324  &,cpl_mpi_comm,mpi_err)
325 
326 
327  !get the state back from da code with mpi_recv
328  call mpi_recv(state_vector,state_dim,mpi_double_precision,1&
329  &,mpi_any_tag,cpl_mpi_comm,status,mpi_err)
330  print*,'Received tag = ',status(mpi_tag)
331 
332 
333  do i = 1,total_timesteps
334  print*,'Timestep = ',i
335 
336  !send the state to da code with mpi_send
337  call mpi_send(state_vector,state_dim,mpi_double_precision,1,1&
338  &,cpl_mpi_comm,mpi_err)
339 
340  !get the state back from da code with mpi_recv
341  call mpi_recv(state_vector,state_dim,mpi_double_precision,1&
342  &,mpi_any_tag,cpl_mpi_comm,status,mpi_err)
343  print*,'Received tag = ',status(mpi_tag)
344  end do
345 
346 
347 
348 
349 
350 
351 
352 
353 
354 
355 
356 
357 
358 
359 
360 
361 
362 
363 
364  call mpi_finalize(mpi_err)
365  print*,'MINIMAL_MODEL_COMMS_v3 finished nicely'
366 end program minimal_model_v3
program minimal_model_v3