Original WRF subgrid support version from John Michalakes without fire
[wrffire.git] / wrfv2_fire / external / RSL / RSL / rsl_hemiforce.c
blobb9ea57d7065ffa7653c1f784d31de7da42033502
1 #define KLUDGE_20000821
3 /***********************************************************************
5 COPYRIGHT
7 The following is a notice of limited availability of the code and
8 Government license and disclaimer which must be included in the
9 prologue of the code and in all source listings of the code.
11 Copyright notice
12 (c) 1977 University of Chicago
14 Permission is hereby granted to use, reproduce, prepare
15 derivative works, and to redistribute to others at no charge. If
16 you distribute a copy or copies of the Software, or you modify a
17 copy or copies of the Software or any portion of it, thus forming
18 a work based on the Software and make and/or distribute copies of
19 such work, you must meet the following conditions:
21 a) If you make a copy of the Software (modified or verbatim)
22 it must include the copyright notice and Government
23 license and disclaimer.
25 b) You must cause the modified Software to carry prominent
26 notices stating that you changed specified portions of
27 the Software.
29 This software was authored by:
31 Argonne National Laboratory
32 J. Michalakes: (630) 252-6646; email: michalak@mcs.anl.gov
33 Mathematics and Computer Science Division
34 Argonne National Laboratory, Argonne, IL 60439
36 ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES
37 OF ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT,
38 AND OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A
39 CONTRACT WITH THE DEPARTMENT OF ENERGY.
41 GOVERNMENT LICENSE AND DISCLAIMER
43 This computer code material was prepared, in part, as an account
44 of work sponsored by an agency of the United States Government.
45 The Government is granted for itself and others acting on its
46 behalf a paid-up, nonexclusive, irrevocable worldwide license in
47 this data to reproduce, prepare derivative works, distribute
48 copies to the public, perform publicly and display publicly, and
49 to permit others to do so. NEITHER THE UNITED STATES GOVERNMENT
50 NOR ANY AGENCY THEREOF, NOR THE UNIVERSITY OF CHICAGO, NOR ANY OF
51 THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
52 ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
53 COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS,
54 PRODUCT, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD
55 NOT INFRINGE PRIVATELY OWNED RIGHTS.
57 ***************************************************************************/
59 #include <stdio.h>
60 #include <stdlib.h>
61 #include "rsl.h"
63 /*@
64 RSL_TO_OH_INFO -- Get the next cell in a packing sequence for forcing.
66 Notes:
68 See also:
70 @*/
72 static rsl_domain_info_t *s_tinfo, *s_oinfo ;
73 static int s_oig, s_ojg ;
74 static int s_p, s_t, s_o ;
75 static int s_msize ;
76 static struct rsl_hemi_rec * s_q, * s_q1 ;
77 static int s_nlen_o ;
78 static int s_mlen_o ;
79 static rsl_point_t *s_tdomain, *s_odomain ;
80 static char * s_pointbuf = NULL ;
82 RSL_TO_OH_INFO ( t_p, o_p, msize_p, seed_p,
83 oig_p, ojg_p, retval_p )
84 int_p
85 t_p /* (I) RSL domain descriptor of this hemi. */
86 ,o_p /* (I) RSL domain descriptor of other hemi. */
87 ,msize_p /* (I) Message size in bytes. */
88 ,seed_p /* (I) =1 start the traversal; =0 (zero) continue traversal */
89 ,oig_p /* (O) Global M index of other domain point. */
90 ,ojg_p /* (O) Global N index of other domain point. */
91 ,retval_p ; /* (O) =1 if a valid point returned; =0 (zero) otherwise. */
93 int kiddex ;
94 int P ;
95 rsl_hemi_rec_t * q ;
96 #ifdef KLUDGE_20000821
97 rsl_hemi_rec_t * qnuke ;
98 rsl_hemi_rec_t * prev ;
99 #endif
100 int p, p1 ;
101 int globalhemiPlist[RSL_MAXPROC][RSL_MAXPROC], work[RSL_MAXPROC][RSL_MAXPROC] ;
103 #ifndef STUBS
104 s_msize = *msize_p ;
105 s_t = *t_p ;
106 s_o = *o_p ;
107 RSL_TEST_ERR( s_t < 0 || s_t > RSL_MAXDOMAINS,
108 "rsl_ready_bcast: bad 'this hemi' descriptor" ) ;
109 RSL_TEST_ERR( s_o < 0 || s_o > RSL_MAXDOMAINS,
110 "rsl_ready_bcast: bad 'other hemi' descriptor" ) ;
111 RSL_TEST_ERR( s_t == s_o,
112 "rsl_ready_bcast: hemispere cannot force itself" ) ;
114 s_tinfo = &( domain_info[s_t]) ;
115 s_oinfo = &( domain_info[s_o]) ;
116 s_mlen_o = s_oinfo->len_m ;
117 s_nlen_o = s_oinfo->len_n ;
118 s_odomain = s_oinfo->domain ;
120 if ( ! s_tinfo->other_hemi_proclist_built )
122 if ( *seed_p )
124 s_oig = 0 ; s_ojg = 0 ;
125 for ( p = 0 ; p < RSL_MAXPROC ; p++ )
127 #ifdef KLUDGE_20000821
128 for ( q = s_tinfo->other_hemi_procbufs[p], prev = NULL ; q ; )
130 if ( q->data ) RSL_FREE( q->data ) ;
131 qnuke = q ;
132 q = q->next ;
133 RSL_FREE( qnuke ) ;
135 #endif
136 s_tinfo->other_hemi_procbufs[p] = NULL ;
137 s_tinfo->hemi_sendPlist[p] = 0 ;
138 for ( p1 = 0 ; p1 < RSL_MAXPROC ; p1++ )
140 globalhemiPlist[p][p1] = 0 ;
144 else
146 s_oig++ ;
147 if ( s_oig >= s_oinfo->len_m )
149 s_oig = 0 ;
150 s_ojg++ ;
151 if ( s_ojg >= s_oinfo->len_n )
153 *retval_p = 0 ;
154 #ifndef KLUDGE_20000821
155 s_tinfo->other_hemi_proclist_built = 1 ; /* FIX 20000818 JM */
156 #endif
158 /* collapse the list and keep only entries that have data associated */
159 /* also fill entries for processors I must send to, indicating the number
160 of columns that go to each processors, or zero for processors I don't
161 send to */
162 for ( p = 0 ; p < RSL_MAXPROC ; p++ )
164 rsl_hemi_rec_t * prev ;
165 for ( q = s_tinfo->other_hemi_procbufs[p], prev = NULL ; q ; )
167 if ( q->data == NULL )
169 if ( prev == NULL )
171 s_tinfo->other_hemi_procbufs[p] = q->next ;
172 RSL_FREE(q) ;
173 q = s_tinfo->other_hemi_procbufs[p] ;
175 else if ( prev->next == q )
177 prev->next = q->next ;
178 RSL_FREE(q) ;
179 q = prev->next ;
181 else
182 RSL_TEST_ERR(1,"internal error") ;
184 else
186 s_tinfo->hemi_sendPlist[p]++ ;
187 prev = q ;
188 q = q->next ;
193 /* mpi all to all communication to share matrix of senders/receivers */
194 MPI_Gather( s_tinfo->hemi_sendPlist, RSL_MAXPROC, MPI_INT,
195 globalhemiPlist, RSL_MAXPROC, MPI_INT,
196 0, rsl_mpi_communicator ) ;
197 /* transpose */
198 for ( p = 0 ; p < RSL_MAXPROC ; p++ )
199 for ( p1 = 0 ; p1 < RSL_MAXPROC ; p1++ )
200 work[p][p1] = globalhemiPlist[p1][p] ;
201 for ( p = 0 ; p < RSL_MAXPROC ; p++ )
202 for ( p1 = 0 ; p1 < RSL_MAXPROC ; p1++ )
203 globalhemiPlist[p][p1] = work[p][p1] ;
204 MPI_Scatter( globalhemiPlist, RSL_MAXPROC, MPI_INT,
205 s_tinfo->hemi_recvPlist, RSL_MAXPROC, MPI_INT,
206 0, rsl_mpi_communicator ) ;
208 return ; /* EARLY RETURN */
212 kiddex = INDEX_2(s_ojg,s_oig,s_mlen_o) ;
213 P = s_odomain[ kiddex ].P ;
214 if ( s_tinfo->other_hemi_procbufs[P] == NULL )
216 q = RSL_MALLOC( rsl_hemi_rec_t, 1 ) ;
218 else
220 q = RSL_MALLOC( rsl_hemi_rec_t, 1 ) ;
221 q->next = s_tinfo->other_hemi_procbufs[P] ;
223 q->oig = s_oig ;
224 q->ojg = s_ojg ;
225 q->data = NULL ;
226 s_tinfo->other_hemi_procbufs[P] = q ;
227 s_q1 = q ;
229 else
231 int * x ;
232 if ( *seed_p )
234 s_p = -1 ;
235 s_q = NULL ;
237 if ( s_q == NULL )
239 s_p++ ;
240 while ( s_tinfo->other_hemi_procbufs[s_p] == NULL ) s_p++ ;
241 if ( s_p >= rsl_nproc_all )
243 *retval_p = 0 ;
244 return ; /* EARLY RETURN */
246 s_q = s_tinfo->other_hemi_procbufs[s_p] ;
248 s_oig = s_q->oig ;
249 s_ojg = s_q->ojg ;
250 s_q1 = s_q ;
251 s_q = s_q->next ;
253 *oig_p = s_oig + 1 ; /* C to Fortran */
254 *ojg_p = s_ojg + 1 ; /* C to Fortran */
255 *retval_p = 1 ;
256 #else
257 RSL_TEST_ERR( 1, "RSL_TO_OH_INFO STUBBED" ) ;
258 #endif
259 return ;
263 RSL_TO_OH_MSG -- Pack force data into a message for a nest point.
265 Notes:
266 See also:
270 RSL_TO_OH_MSG ( nbuf_p, buf )
271 int_p
272 nbuf_p ; /* (I) Number of bytes to be packed. */
273 char *
274 buf ; /* (I) Buffer containing the data to be packed. */
276 int kiddex ;
277 int nbuf ;
278 int P ;
280 RSL_TEST_ERR(buf==NULL,"2nd argument is NULL. Field allocated?") ;
281 nbuf = *nbuf_p ;
282 if ( s_q1->data == NULL )
284 s_q1->data = RSL_MALLOC( char, s_msize ) ;
285 s_q1->curs = 0 ;
287 if ( s_q1->curs+nbuf >= s_msize )
289 sprintf(mess,"RSL_TO_OH_MSG: store of %d bytes would overflow %d sized buffer.\n",nbuf,s_msize ) ;
290 RSL_TEST_ERR(1,mess) ;
292 bcopy( buf, &(s_q1->data[s_q1->curs]), nbuf ) ;
293 s_q1->curs += nbuf ;
297 RSL_FORCE_HEMI -- Convey forcing data from this hemi to other hemi
299 Notes:
301 See also:
304 RSL_FORCE_HEMI ()
306 int P ;
307 int msglen, mdest, mtag ;
308 int ii ;
309 int ig, jg ;
310 char * recvbuf, * sendbuf ;
311 rsl_hemi_rec_t * q ;
313 /* post receives */
315 for ( P = 0 ; P < rsl_nproc_all ; P++ )
317 if ( s_tinfo->hemi_recvPlist[P] > 0 )
319 msglen = s_msize * s_tinfo->hemi_recvPlist[P] + 3*sizeof(int) ;
320 recvbuf = buffer_for_proc( P, msglen, RSL_RECVBUF ) ;
321 mtag = MTYPE_FROMTO( MSG_FROM_PARENT, P, rsl_myproc ) ;
322 #ifdef DEBUGGAL
323 fprintf(stderr,"Posting receive on tag %d\n",mtag ) ;
324 #endif
325 RSL_RECVBEGIN( recvbuf, msglen, mtag ) ;
326 s_tinfo->hemi_recv_tags[P] = mtag ; /* store tag */
330 /* do sends */
332 for ( P = 0 ; P < rsl_nproc_all ; P++ )
334 #ifdef DEBUGGAL
335 fprintf(stderr,"s_tinfo->hemi_sendPlist[P] %d\n",s_tinfo->hemi_sendPlist[P]) ;
336 #endif
337 if ( s_tinfo->hemi_sendPlist[P] > 0 )
339 int curs ;
340 int endofdata ;
342 /* oig,ojg,nbytes,buffer * # of points + end of data */
343 msglen = (3*sizeof(int)+s_msize)*s_tinfo->hemi_sendPlist[P]+1*sizeof(int) ;
344 sendbuf = buffer_for_proc( P, msglen, RSL_SENDBUF ) ;
345 curs = 0 ;
346 for ( q = s_tinfo->other_hemi_procbufs[P] ; q ; q = q->next )
348 #ifdef DEBUGGAL
350 int *dp ;
351 dp = (int *) q->data ;
352 fprintf(stderr,"> curs %d, msglen %d msize %d (%d %d) data %d\n", curs, msglen, s_msize, q->oig, q->ojg, *dp ) ;
354 #endif
355 bcopy( &(q->oig), &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
356 bcopy( &(q->ojg), &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
357 bcopy( &(q->curs), &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
358 bcopy( q->data, &(sendbuf[curs]), q->curs) ; curs += q->curs ;
360 endofdata = RSL_INVALID ;
361 bcopy( &endofdata, &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
362 mtag = MTYPE_FROMTO( MSG_FROM_PARENT, rsl_myproc, P ) ;
363 #ifdef DEBUGGAL
364 fprintf(stderr,"sending sendbuf to %d, curs = %d\n",P,curs) ;
365 #endif
366 RSL_SEND( sendbuf, curs, mtag, P ) ;
372 RSL_FROM_TH_INFO -- Get the next cell in a unpacking sequence for forcing.
374 Notes:
376 See also:
379 static int s_endofdata, s_remaining, s_ndata, s_curs ;
380 static char * s_recvbuf ;
382 RSL_FROM_TH_INFO ( seed_p, oig_p, ojg_p, retval_p )
383 int_p
384 seed_p /* (I) =1 if first call; =0 otherwise */
385 ,oig_p /* (O) Global index in M dimension of nest. */
386 ,ojg_p /* (O) Global index in N dimension of nest. */
387 ,retval_p ; /* (O) Return value; =1 valid point, =0 done. */
389 int mtag ;
391 #ifdef DEBUGGAL
392 fprintf(stderr,"RSL_FROM_TH_INFO seed = %d, s_endofdata %d\n",*seed_p,s_endofdata) ;
393 #endif
395 if ( *seed_p == 1 )
397 if ( s_pointbuf != NULL ) RSL_FREE(s_pointbuf) ;
398 s_pointbuf = RSL_MALLOC( char, 2*s_msize ) ; /* 2 times for safety */
399 s_p = 0 ;
400 s_endofdata = 1 ;
403 nextproc:
404 if ( s_endofdata )
406 while ( s_tinfo->hemi_recvPlist[s_p] <= 0 ) s_p++ ;
407 if ( s_p >= rsl_nproc_all )
409 *retval_p = 0 ;
410 #ifdef DEBUGGAL
411 fprintf(stderr,"EARLY RETURN retval = 0\n") ;
412 #endif
413 return ; /* EARLY RETURN */
415 mtag = s_tinfo->hemi_recv_tags[s_p] ;
416 #ifdef DEBUGGAL
417 fprintf(stderr,"Waiting for receive on tag %d\n",mtag ) ;
418 #endif
419 RSL_RECVEND ( mtag ) ;
420 #ifdef DEBUGGAL
421 fprintf(stderr,"got receive\n") ;
422 #endif
423 s_recvbuf = buffer_for_proc( s_p, 0, RSL_RECVBUF ) ;
424 s_p++ ;
425 s_curs = 0 ;
426 s_endofdata = 0 ;
429 #ifdef DEBUGGAL
430 fprintf(stderr,"before bcopy %d, s_recvbuf %08x\n",s_curs, s_recvbuf) ;
431 #endif
433 bcopy ( &(s_recvbuf[s_curs]), oig_p, sizeof(int) ) ; s_curs += sizeof(int) ;
434 if ( *oig_p == RSL_INVALID )
436 #ifdef DEBUGGAL
437 fprintf(stderr,"hit end of data for s_p %d, %d\n", s_p, *oig_p ) ;
438 #endif
439 s_endofdata = 1 ;
440 goto nextproc ;
442 bcopy ( &(s_recvbuf[s_curs]), ojg_p, sizeof(int) ) ; s_curs += sizeof(int) ;
443 bcopy ( &(s_recvbuf[s_curs]), &s_ndata, sizeof(int) ) ; s_curs += sizeof(int) ;
444 bcopy ( &(s_recvbuf[s_curs]), s_pointbuf, s_ndata ) ; s_curs += s_ndata ;
445 s_remaining = s_ndata ;
446 #ifdef DEBUGGAL
447 fprintf(stderr,"s_remaining = %d\n",s_remaining) ;
448 #endif
450 (*oig_p) ++ ;
451 (*ojg_p) ++ ;
452 #ifdef DEBUGGAL
453 fprintf(stderr,"RETURN oig ojg %d %d\n", *oig_p, *ojg_p ) ;
454 #endif
455 *retval_p = 1 ;
456 return ;
460 RSL_FROM_TH_MSG -- Unpack feedback data into a nest point.
462 Notes:
465 RSL_FROM_TH_MSG ( len_p, buf )
466 int_p
467 len_p ; /* (I) Number of bytes to unpack. */
468 char *
469 buf ; /* (O) Destination buffer. */
471 if ( *len_p <= 0 ) return ;
472 if ( *len_p > s_remaining )
474 sprintf(mess,
475 "RSL_FROM_TH_MSG:\n Requested number of bytes (%d) exceeds %d, the number remaining for this point.\n", *len_p, s_remaining) ;
476 RSL_TEST_WRN(1,mess) ;
478 bcopy( &(s_pointbuf[s_ndata-s_remaining]),
479 buf,
480 *len_p ) ;
482 s_remaining -= *len_p ;
485 /* retval =1 if point is local, =0 otherwise */
486 RSL_POINT_ON_PROC ( d_p, ig_p, jg_p, retval_p )
487 int_p d_p, ig_p, jg_p, retval_p ;
489 int d ;
490 int kiddex ;
491 int P ;
492 int ig, jg ;
494 rsl_domain_info_t * info ;
495 rsl_point_t * domain ;
496 ig = *ig_p - 1 ;
497 jg = *jg_p - 1 ;
499 d = *d_p ;
500 RSL_TEST_ERR( d < 0 || d > RSL_MAXDOMAINS,
501 "rsl_ready_bcast: bad 'this hemi' descriptor" ) ;
502 info = &( domain_info[d]) ;
503 /* added 12/27/01 -- JM */
504 if ( ig < 0 || ig >= info->len_m ||
505 jg < 0 || jg >= info->len_n ) { *retval_p = 0 ; return ; }
506 domain = info->domain ;
507 kiddex = INDEX_2(jg,ig,info->len_m ) ;
508 P = domain[ kiddex ].P ;
509 *retval_p = 0 ;
510 if ( P == rsl_myproc ) *retval_p = 1 ;
511 return ;
514 /* given a global point, return the processor number */
515 RSL_PROC_FOR_POINT ( d_p, ig_p, jg_p, retval_p )
516 int_p d_p, ig_p, jg_p, retval_p ;
518 int d ;
519 int kiddex ;
520 int P ;
521 int ig, jg ;
523 rsl_domain_info_t * info ;
524 rsl_point_t * domain ;
525 ig = *ig_p - 1 ;
526 jg = *jg_p - 1 ;
528 d = *d_p ;
529 RSL_TEST_ERR( d < 0 || d > RSL_MAXDOMAINS,
530 "rsl_point_on_proc: bad descriptor" ) ;
531 info = &( domain_info[d]) ;
532 domain = info->domain ;
533 kiddex = INDEX_2(jg,ig,info->len_m ) ;
534 *retval_p = domain[ kiddex ].P ;
535 return ;