source: CIVL/examples/mpi-omp/AMG2013/struct_mv/new_assemble.c

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 34.2 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14
15#include "headers.h"
16
17
18
19int hypre_FillResponseStructAssembleAP(void *, int, int, void *,
20 MPI_Comm, void **, int *);
21
22
23
24
25#define DEBUG 0
26#define AP_STAT 0
27#define NEIGH_PRINT 0
28
29#if DEBUG
30char filename[255];
31FILE *file;
32int my_rank;
33#endif
34
35
36
37/*--------------------------------------------------------------------------
38 * new version of
39 * hypre_StructGridAssemble (for use with "no global partition" option)
40 * AHB 6/05
41 *------------------------------------------------------------------------*/
42
43
44int hypre_StructGridAssembleWithAP( hypre_StructGrid *grid )
45{
46
47
48
49 int ierr = 0;
50 int tmp_i;
51
52 int size, global_num_boxes, num_local_boxes;
53 int i, j, d, k, index;
54 int num_procs, myid;
55 int *sendbuf8, *recvbuf8, *sendbuf2, *recvbuf2;
56 int min_box_size, max_box_size;
57 int global_min_box_size, global_max_box_size;
58 int *ids;
59 int max_regions, max_refinements, ologp;
60 double gamma;
61 hypre_Index min_index, max_index;
62
63
64 int prune;
65
66 hypre_Box *box;
67
68
69 MPI_Comm comm = hypre_StructGridComm(grid);
70 hypre_Box *bounding_box = hypre_StructGridBoundingBox(grid);
71 hypre_BoxArray *local_boxes = hypre_StructGridBoxes(grid);
72 int dim = hypre_StructGridDim(grid);
73 hypre_BoxNeighbors *neighbors = hypre_StructGridNeighbors(grid);
74 int max_distance = hypre_StructGridMaxDistance(grid);
75 hypre_IndexRef periodic = hypre_StructGridPeriodic(grid);
76
77 int *local_boxnums;
78
79 double dbl_global_size, tmp_dbl;
80
81 hypre_BoxArray *my_partition;
82 int *part_ids, *part_boxnums;
83
84 int *proc_array, proc_count, proc_alloc, count;
85 int *tmp_proc_ids = NULL;
86
87 int max_response_size;
88 int *ap_proc_ids, *send_buf, *send_buf_starts;
89 int *response_buf, *response_buf_starts;
90
91 hypre_BoxArray *neighbor_boxes, *n_boxes_copy;
92 int *neighbor_proc_ids, *neighbor_boxnums;
93
94 int *order_index, *delete_array;
95 int tmp_id, start, first_local;
96
97 int grow, grow_array[6];
98 hypre_Box *grow_box;
99
100
101 int *numghost;
102 int ghostsize;
103 hypre_Box *ghostbox;
104
105 hypre_StructAssumedPart *assumed_part;
106 hypre_DataExchangeResponse response_obj;
107
108 int px = hypre_IndexX(periodic);
109 int py = hypre_IndexY(periodic);
110 int pz = hypre_IndexZ(periodic);
111
112 int i_periodic = px ? 1 : 0;
113 int j_periodic = py ? 1 : 0;
114 int k_periodic = pz ? 1 : 0;
115
116 int num_periods, multiple_ap, p;
117 hypre_Box *result_box, *period_box;
118 hypre_Index *pshifts;
119
120 hypre_IndexRef pshift;
121
122#if NEIGH_PRINT
123 double start_time, end_time;
124
125#endif
126
127
128
129/*---------------------------------------------
130 Step 1: Initializations
131 -----------------------------------------------*/
132
133 prune = 1; /* default is to prune */
134
135 num_local_boxes = hypre_BoxArraySize(local_boxes);
136
137 num_periods = (1+2*i_periodic) * (1+2*j_periodic) * (1+2*k_periodic);
138
139
140 MPI_Comm_size(comm, &num_procs);
141 MPI_Comm_rank(comm, &myid);
142
143
144
145/*---------------------------------------------
146 Step 2: Determine the global size, total number of boxes,
147 and global bounding box.
148 Also get the min and max box sizes
149 since it is convenient to do so.
150 -----------------------------------------------*/
151
152 if (neighbors == NULL)
153 {
154
155 /*these may not be needed - check later */
156 ids = hypre_TAlloc(int, num_local_boxes);
157
158 /* for the vol and number of boxes */
159 sendbuf2 = hypre_CTAlloc(int, 2);
160 recvbuf2 = hypre_CTAlloc(int, 2);
161 size = 0;
162
163 bounding_box = hypre_BoxCreate();
164 grow_box = hypre_BoxCreate();
165
166
167 if (num_local_boxes)
168 {
169
170 min_box_size = hypre_BoxVolume( hypre_BoxArrayBox(local_boxes, 0));
171 max_box_size = hypre_BoxVolume( hypre_BoxArrayBox(local_boxes, 0));
172
173
174 /* initialize min and max */
175 for (d=0; d<3; d++)
176 {
177 hypre_IndexD(min_index, d) = pow(2,30);
178 hypre_IndexD(max_index, d) = -pow(2,30);
179 }
180
181
182 hypre_ForBoxI(i, local_boxes)
183 {
184 box = hypre_BoxArrayBox(local_boxes, i);
185 /* get global size and number of boxes */
186 tmp_i = hypre_BoxVolume(box);
187 size += tmp_i;
188 min_box_size = hypre_min(min_box_size, tmp_i);
189 max_box_size = hypre_max(max_box_size, tmp_i);
190
191
192 /* set id */
193 ids[i] = i;
194
195
196 /* 1/3/05 we need this for the case of holes in the domain. (I had
197 commented
198 it out on 12/04 - as I thought this was not necessary. */
199
200
201 /* zero volume boxes - still look at for getting the bounding box */
202 if (hypre_BoxVolume(box) == 0) /* zero volume boxes - still count */
203 {
204 hypre_CopyBox(box, grow_box);
205 for (d = 0; d < 3; d++)
206 {
207 if(!hypre_BoxSizeD(box, d))
208 {
209 grow = (hypre_BoxIMinD(box, d) - hypre_BoxIMaxD(box, d) + 1)/2;
210 grow_array[2*d] = grow;
211 grow_array[2*d+1] = grow;
212 }
213 else
214 {
215 grow_array[2*d] = 0;
216 grow_array[2*d+1] = 0;
217 }
218 }
219 /* expand the box */
220 hypre_BoxExpand(grow_box, grow_array);
221 box = grow_box; /*pointer copy*/
222 }
223 /*now we have a vol > 0 box */
224
225
226 for (d = 0; d < dim; d++) /* for each dimension */
227 {
228 hypre_IndexD(min_index, d) = hypre_min( hypre_IndexD(min_index, d),
229 hypre_BoxIMinD(box, d));
230 hypre_IndexD(max_index, d) = hypre_max( hypre_IndexD(max_index, d),
231 hypre_BoxIMaxD(box, d));
232 }
233
234 }/*end for each box loop */
235
236 /* bounding box extents */
237 hypre_BoxSetExtents(bounding_box, min_index, max_index);
238
239 }
240 else /* no boxes owned*/
241 {
242 for (d=0; d<3; d++)
243 {
244 hypre_BoxIMinD(bounding_box, d) = pow(2,30) ;
245 hypre_BoxIMaxD(bounding_box, d) = -pow(2,30);
246 }
247
248 min_box_size = 0;
249 max_box_size = 0;
250 size = 0;
251 }
252
253
254 /* if dim < 3 then set the extra dimensions to zero */
255 for (d = dim; d < 3; d++)
256 {
257 hypre_BoxIMinD(bounding_box, d) = 0;
258 hypre_BoxIMaxD(bounding_box, d) = 0;
259 }
260
261
262 /* populate the vol and number of boxes buffer */
263 sendbuf2[0] = size;
264 sendbuf2[1] = hypre_BoxArraySize(local_boxes);
265
266
267 /* set local size (volume) */
268 hypre_StructGridLocalSize(grid) = size;
269
270 MPI_Allreduce(sendbuf2, recvbuf2, 2, MPI_INT, MPI_SUM, comm);
271
272 /* now set global size */
273 hypre_StructGridGlobalSize(grid) = recvbuf2[0]; /*this could easily overflow!
274 need to change to a double
275 in struct grid data
276 structure */
277 global_num_boxes = recvbuf2[1];
278
279 /* in case of overflow (until we change the datatype in the struct grid
280 object) - we use this global vol for our calculations */
281 tmp_dbl = (double) size;
282 MPI_Allreduce(&tmp_dbl, &dbl_global_size, 1, MPI_DOUBLE, MPI_SUM, comm);
283 hypre_TFree(sendbuf2);
284 hypre_TFree(recvbuf2);
285
286#if AP_STAT
287 if (myid ==0)
288 {
289 printf("myid = %d, GLOBAL number of boxes = %d\n", myid, global_num_boxes);
290 }
291#endif
292
293 /* don't need the grow_box */
294 hypre_BoxDestroy(grow_box);
295
296
297 /* now get the bounding box and min and max box sizes */
298
299 sendbuf8 = hypre_CTAlloc(int, 8);
300 recvbuf8 = hypre_CTAlloc(int, 8);
301
302
303 /*get the min global lower extents and max upper extents */
304 /* note: min(x)= -max(-x) */
305
306 for (d = 0; d < 3; d++)
307 {
308 sendbuf8[d] = hypre_BoxIMinD(bounding_box, d);
309 sendbuf8[d+3] = -hypre_BoxIMaxD(bounding_box, d);
310 }
311
312 /*also collect the min and max box sizes*/
313 sendbuf8[6] = min_box_size;
314 sendbuf8[7] = -max_box_size;
315
316
317 MPI_Allreduce(sendbuf8, recvbuf8, 8, MPI_INT, MPI_MIN, comm);
318
319
320 for (d = 0; d < 3; d++)
321 {
322 hypre_BoxIMinD(bounding_box, d) = recvbuf8[d];
323 hypre_BoxIMaxD(bounding_box, d) = -recvbuf8[d+3];
324 }
325
326 global_min_box_size = recvbuf8[6];
327 global_max_box_size = -recvbuf8[7];
328
329
330 hypre_TFree(sendbuf8);
331 hypre_TFree(recvbuf8);
332
333
334
335#if 0
336 printf("myid = %d, GLOBAL min box size = %d\n", myid, global_min_box_size);
337 printf("myid = %d, GLOBAL max box size = %d\n", myid, global_max_box_size);
338#endif
339
340
341 /* in theory there should be at least one box with vol > 0, but
342 just in case: */
343 if (global_num_boxes == 0)
344 {
345 for (d = 0; d < 3; d++)
346 {
347 hypre_BoxIMinD(bounding_box, d) = 0;
348 hypre_BoxIMaxD(bounding_box, d) = 0;
349 }
350 global_min_box_size = 1;
351 global_max_box_size = 1;
352
353 /* do we need to set the global grid size ? */
354 hypre_StructGridGlobalSize(grid) = 1;
355
356 }
357
358 if (hypre_BoxVolume(bounding_box) == 0)
359 {
360 if (myid ==0) printf("Error: bounding box has zero volume - "
361 "this shouldn't happen!\n");
362 return -1;
363 }
364
365 hypre_StructGridBoundingBox(grid) = bounding_box;
366
367
368
369
370 /* ids are in order r: 0 - num(local-boxes) - 1*/
371 /* only set if they have not already been set! */
372 if ( hypre_StructGridIDs(grid) == NULL)
373 {
374 hypre_StructGridIDs(grid) = ids;
375 }
376 else
377 {
378 /* not needed! */
379 hypre_TFree(ids);
380 }
381
382
383
384 }
385
386
387
388
389/*---------------------------------------------
390 Step 3: Create an assumed partition
391 -----------------------------------------------*/
392
393/* want #regions < #procs */
394/* want #regions bounded */
395
396
397 if (neighbors == NULL)
398 {
399 /* Initializations - these should probably be
400 adjusted depending on the problem */
401
402 /* what is the max number of regions we can store ? */
403 /* what is the max number of refinements we can do? */
404 /* what fraction of the region do we want to be full? */
405
406 /* estimate of log(num_procs) */
407
408 d = num_procs/2;
409 ologp = 0;
410 while ( d > 0)
411 {
412 d = d/2; /* note - d is an int - so this is floored */
413 ologp++;
414 }
415
416 max_regions = hypre_min(pow(2, ologp+1), 10*ologp);
417 /* max_regions = 100; */
418
419 max_refinements = ologp;
420 /* max_refinements = 1;*/
421
422 /* TEMP - lots of refinement */
423 /*max_regions = num_procs;
424 max_refinements = ologp + 10;*/
425
426 /* new 1/10/04*/
427 /* max_refinements = ologp/2;
428 max_regions = 64;*/
429
430
431
432 gamma = .6; /* percentage a region must be full to
433 avoid refinement */
434
435
436 /* assign boxnums */
437 local_boxnums = hypre_CTAlloc(int, num_local_boxes);
438 for (i = 0; i< num_local_boxes; i++)
439 {
440 local_boxnums[i] = i;
441 }
442
443
444 ierr = hypre_StructAssumedPartitionCreate(dim, bounding_box, dbl_global_size,
445 global_num_boxes,
446 local_boxes, local_boxnums,
447 max_regions,
448 max_refinements, gamma,
449 comm, &assumed_part);
450
451
452 hypre_TFree(local_boxnums);
453
454
455
456 /*Now we have the assumed partition built */
457
458#if AP_STAT
459 my_partition = hypre_StructAssumedPartMyPartition(assumed_part);
460 hypre_ForBoxI(i, my_partition)
461 {
462 box = hypre_BoxArrayBox(my_partition, i);
463 printf("*****myid = %d: MY ASSUMED Partitions (%d): (%d, %d, %d) "
464 "x (%d, %d, %d)\n",
465 myid, i,
466 hypre_BoxIMinX(box),
467 hypre_BoxIMinY(box),
468 hypre_BoxIMinZ(box),
469 hypre_BoxIMaxX(box),
470 hypre_BoxIMaxY(box),
471 hypre_BoxIMaxZ(box));
472 }
473#endif
474
475 my_partition = hypre_StructAssumedPartMyPartitionBoxes(assumed_part);
476 part_ids = hypre_StructAssumedPartMyPartitionProcIds(assumed_part);
477 part_boxnums = hypre_StructAssumedPartMyPartitionBoxnums(assumed_part);
478
479#if AP_STAT
480 printf("*****myid = %d: I have %d boxes in my AP from %d procs\n", myid,
481 hypre_BoxArraySize(my_partition),
482 hypre_StructAssumedPartMyPartitionNumDistinctProcs(assumed_part));
483
484#endif
485
486#if AP_STAT
487
488 hypre_ForBoxI(i, my_partition)
489 {
490 box = hypre_BoxArrayBox(my_partition, i);
491 printf("*****myid = %d: BOXES in MY AP (number %d): (%d, %d, %d) "
492 "x (%d, %d, %d), boxnum = %d, owned by proc %d\n",
493 myid, i,
494 hypre_BoxIMinX(box),
495 hypre_BoxIMinY(box),
496 hypre_BoxIMinZ(box),
497 hypre_BoxIMaxX(box),
498 hypre_BoxIMaxY(box),
499 hypre_BoxIMaxZ(box), part_boxnums[i], part_ids[i]);
500 }
501#endif
502
503 }
504
505
506
507/*---------------------------------------------
508 Step 5: Use the assumed partition to find a
509 shortened list of *potential* neighbors
510 -----------------------------------------------*/
511
512
513/* set neighbors */
514 if (neighbors == NULL)
515 {
516
517 /* need to pass in a list of boxes (to BoxNeighborsAssemble)-
518 not all the boxes - that
519 contains potential nearest neighbors. Also a corresponding list of
520 processor numbers (same length). Also provide an id
521 number for each of these and
522 an index that indicates which one is the first of the local boxes.
523 Previously the list of boxes contained
524 ALL of the boxes in the entire domain! */
525
526 /*first determine the shifting needed for periodic boxes */
527 pshifts = hypre_CTAlloc(hypre_Index, num_periods); /* this is deallocated
528 in boxneighbordestroy */
529 pshift = pshifts[0];
530 hypre_ClearIndex(pshift);
531
532 if( num_periods > 1 )
533 {
534 int ip, jp, kp;
535 p = 1;
536 for (ip = -i_periodic; ip <= i_periodic; ip++)
537 {
538 for (jp = -j_periodic; jp <= j_periodic; jp++)
539 {
540 for (kp = -k_periodic; kp <= k_periodic; kp++)
541 {
542 if( !(ip == 0 && jp == 0 && kp == 0) )
543 {
544 pshift = pshifts[p];
545 hypre_SetIndex(pshift, ip*px, jp*py, kp*pz);
546 p++;
547 }
548 }
549 }
550 }
551 }
552
553 /* now we go through all our boxes and
554 we grow them in each dimension before checking to see whose AP
555 they lie in - in case we have any near boundaries. We need to
556 store the list of processors form the AP that our boxes intersect*/
557
558
559 /*to store info from one box */
560 proc_count = 0;
561 proc_alloc = 8;
562 proc_array = hypre_CTAlloc(int, proc_alloc);
563
564 /* probably there will mostly be one proc per box */
565 size = 2*hypre_BoxArraySize(local_boxes);
566 tmp_proc_ids = hypre_CTAlloc(int, size);
567 count = 0;
568
569 box = hypre_BoxCreate();
570 result_box = hypre_BoxCreate();
571 period_box = hypre_BoxCreate();
572
573 /* loop through all boxes */
574 hypre_ForBoxI(i, local_boxes)
575 {
576 multiple_ap = 0;
577 hypre_CopyBox(hypre_BoxArrayBox(local_boxes, i) ,box);
578 hypre_BoxExpandConstant(box, max_distance);
579
580 if (hypre_BoxVolume(box) == 0) /* zero volume boxes - still
581 grow more*/
582 {
583 for (d = 0; d < 3; d++)
584 {
585 if(!hypre_BoxSizeD(box, d))
586 {
587 grow = (hypre_BoxIMinD(box, d) - hypre_BoxIMaxD(box, d) + 1)/2;
588 grow_array[2*d] = grow;
589 grow_array[2*d+1] = grow;
590 }
591 else
592 {
593 grow_array[2*d] = 0;
594 grow_array[2*d+1] = 0;
595 }
596 }
597 /* expand the box */
598 hypre_BoxExpand(box, grow_array);
599 }
600
601 if (num_periods > 1)
602 {
603 /*treat periodic separately if it's at the edge of the domain -
604 then we need to look for periodic neighbors*/
605 hypre_IntersectBoxes( box, bounding_box, result_box);
606 if ( hypre_BoxVolume(result_box) < hypre_BoxVolume(box) )
607 /* on the edge of the bounding box */
608 {
609 multiple_ap = 1;
610 }
611 }
612
613 if (!multiple_ap) /* only one assumed partition (AP) call */
614 {
615 hypre_StructAssumedPartitionGetProcsFromBox(assumed_part, box, &proc_count,
616 &proc_alloc, &proc_array);
617 if ((count + proc_count) > size)
618 {
619 size = size + proc_count + 2*(hypre_BoxArraySize(local_boxes)-i);
620 tmp_proc_ids = hypre_TReAlloc(tmp_proc_ids, int, size);
621 }
622 for (j = 0; j< proc_count; j++)
623 {
624 tmp_proc_ids[count] = proc_array[j];
625 count++;
626 }
627
628 }
629
630 else /* this periodic box needs to be shifted as it is near a boundary -
631 so we'll have multiple AP calls*/
632 {
633 for (k=0; k < num_periods; k++)
634 {
635 /* get the periodic box (k=0 is the actual box) */
636 hypre_CopyBox(box, period_box);
637 pshift = pshifts[k];
638 hypre_BoxShiftPos(period_box, pshift);
639 /* see if the shifted box intersects the domain */
640 hypre_IntersectBoxes(period_box, bounding_box, result_box);
641 if (hypre_BoxVolume(result_box) > 0)
642 {
643 hypre_StructAssumedPartitionGetProcsFromBox(assumed_part, period_box,
644 &proc_count, &proc_alloc,
645 &proc_array);
646 if ((count + proc_count) > size)
647 {
648 size = size + proc_count +
649 2*(hypre_BoxArraySize(local_boxes)-i);
650 tmp_proc_ids = hypre_TReAlloc(tmp_proc_ids, int, size);
651 }
652 for (j = 0; j< proc_count; j++)
653 {
654 tmp_proc_ids[count] = proc_array[j];
655 count++;
656 }
657 }
658 }
659 } /* end of if multiple AP calls */
660
661 } /* end of loop through boxes */
662
663 hypre_TFree(proc_array);
664 hypre_BoxDestroy(box);
665 hypre_BoxDestroy(result_box);
666 hypre_BoxDestroy(period_box);
667
668
669 /* now get rid of redundencies in tmp_proc_ids - put in ap_proc_ids*/
670 qsort0(tmp_proc_ids, 0, count-1);
671 proc_count = 0;
672 ap_proc_ids = hypre_CTAlloc(int, count);
673
674 if (count)
675 {
676 ap_proc_ids[0] = tmp_proc_ids[0];
677 proc_count++;
678 }
679 for (i = 1; i < count; i++)
680 {
681 if (tmp_proc_ids[i] != ap_proc_ids[proc_count-1])
682 {
683 ap_proc_ids[proc_count] = tmp_proc_ids[i];
684 proc_count++;
685 }
686 }
687
688 hypre_TFree(tmp_proc_ids);
689
690#if NEIGH_PRINT
691 printf("*****myid = %d: I have %d procs to contact for their AP boxes \n", myid, proc_count);
692#endif
693
694
695
696 /* now we have a sorted list with no duplicates in ap_proc_ids */
697 /* for each of these processor ids, we need to get the boxes in their
698 assumed partition as these are potential neighbors */
699
700
701 /* need to do an exchange data for this! (use flag 2 since flag 1 is in create
702 assumed partition)*/
703 /* contact message will be empty. the response message should be all the
704 boxes + corresp. proc id + boxnum in that processor's assumed partition
705 (so 8 ints for each)*/
706
707 /*response object*/
708 response_obj.fill_response = hypre_FillResponseStructAssembleAP;
709 response_obj.data1 = assumed_part; /* needed to fill responses*/
710 response_obj.data2 = NULL;
711
712 send_buf = NULL;
713 send_buf_starts = hypre_CTAlloc(int, proc_count + 1);
714 for (i=0; i< proc_count+1; i++)
715 {
716 send_buf_starts[i] = 0;
717 }
718
719
720 response_buf = NULL; /*this and the next are allocated in exchange data */
721 response_buf_starts = NULL;
722
723
724 /*we expect back an array of items consisting of 8 integers*/
725 size = 8*sizeof(int);
726
727 /* this needs to be the same for all processors */
728 max_response_size = global_num_boxes/num_procs + 10;
729 max_response_size = (global_num_boxes/num_procs)*2;
730
731
732#if NEIGH_PRINT
733 start_time = MPI_Wtime();
734#endif
735
736
737
738 hypre_DataExchangeList(proc_count, ap_proc_ids,
739 send_buf, send_buf_starts,
740 0, size, &response_obj, max_response_size, 2,
741 comm, (void**) &response_buf, &response_buf_starts);
742
743
744
745
746#if NEIGH_PRINT
747 end_time = MPI_Wtime();
748 end_time = end_time - start_time;
749 printf("myid = %d EXCHANGE 2 time = %f sec.\n", myid, end_time);
750#endif
751
752 /*clean right away */
753 hypre_TFree(send_buf_starts);
754 hypre_TFree(ap_proc_ids);
755 /*we no longer need the assumed partition */
756 hypre_StructAssumedPartitionDestroy(assumed_part);
757
758#if NEIGH_PRINT
759 start_time = MPI_Wtime();
760#endif
761
762 /*now we have the potential list of neighbor boxes and corresponding
763 processors in response_buf*/
764
765 /*how many neighbor boxes do we have?*/
766 size = response_buf_starts[proc_count];
767 hypre_TFree(response_buf_starts);
768
769 neighbor_proc_ids = hypre_CTAlloc(int, size);
770 neighbor_boxnums = hypre_CTAlloc(int, size);
771 neighbor_boxes = hypre_BoxArrayCreate(size);
772 box = hypre_BoxCreate();
773 order_index = hypre_CTAlloc(int, size);
774
775 index = 0;
776 for (i=0; i< size; i++) /* for each neighbor box */
777 {
778 neighbor_proc_ids[i] = response_buf[index++];
779 neighbor_boxnums[i] = response_buf[index++];
780 for (d=0; d< 3; d++)
781 {
782 hypre_BoxIMinD(box, d) = response_buf[index++];
783 hypre_BoxIMaxD(box, d) = response_buf[index++];
784 }
785 hypre_CopyBox(box, hypre_BoxArrayBox(neighbor_boxes, i));
786 order_index[i] = i;
787 }
788
789 hypre_BoxDestroy(box);
790 hypre_TFree(response_buf);
791
792
793
794 /* now we have an array of boxes, boxnums, and proc_ids -
795 not in any order and there
796 may be duplicate pairs of proc id and boxes - this list includes all of MY
797 boxes as well*/
798#if 0
799 printf("*****myid = %d: I have %d potential neighbor boxes \n", myid, size);
800
801 hypre_ForBoxI(i, neighbor_boxes)
802 {
803 box = hypre_BoxArrayBox(neighbor_boxes, i);
804 printf("*****myid = %d: BOXES in my list of potential neighbors "
805 "(number %d): (%d, %d, %d) x (%d, %d, %d) , boxnum = %d, "
806 "proc %d\n",
807 myid, i,
808 hypre_BoxIMinX(box),
809 hypre_BoxIMinY(box),
810 hypre_BoxIMinZ(box),
811 hypre_BoxIMaxX(box),
812 hypre_BoxIMaxY(box),
813 hypre_BoxIMaxZ(box), neighbor_boxnums[i], neighbor_proc_ids[i]);
814 }
815
816#endif
817
818
819 /* to set neighbors - pass in reduced list of potential neighbor boxes.
820 We'll put processors in ascending order. Also we need to remove
821 duplicate boxes first */
822
823 size = hypre_BoxArraySize(neighbor_boxes);
824
825 /*sort on proc_id - move boxnums and order_index*/
826 hypre_qsort3i(neighbor_proc_ids, neighbor_boxnums, order_index, 0, size-1);
827
828 delete_array = hypre_CTAlloc(int, size);
829 index = 0;
830 first_local = 0;
831
832 /*now within each proc_id, we need to sort the boxnums and order index*/
833 if (size)
834 {
835 tmp_id = neighbor_proc_ids[0];
836 }
837 start = 0;
838 for (i=1; i< size; i++)
839 {
840 if (neighbor_proc_ids[i] != tmp_id)
841 {
842 hypre_qsort2i(neighbor_boxnums, order_index, start, i-1);
843 /*now find duplicate boxnums */
844 for (j=start+1; j< i; j++)
845 {
846 if (neighbor_boxnums[j] == neighbor_boxnums[j-1])
847 {
848 delete_array[index++] = j;
849 }
850 }
851 /* update start and tmp_id */
852 start = i;
853 tmp_id = neighbor_proc_ids[i];
854 if(tmp_id == myid )
855 {
856 /* subtract the value of index as this is how many previous we will
857 delete */
858 first_local = i - index;
859 }
860 }
861 }
862
863 /* final sort and purge (the last group doesn't
864 get caught in the above loop) */
865 if (size > 1)
866 {
867 hypre_qsort2i(neighbor_boxnums, order_index, start, size-1);
868 /*now find duplicate boxnums */
869 for (j=start+1; j<size; j++)
870 {
871 if (neighbor_boxnums[j] == neighbor_boxnums[j-1])
872 {
873 delete_array[index++] = j;
874 }
875 }
876 }
877
878 /* now index = the number in delete_array */
879
880 /*now sort the boxes according to index_order */
881 n_boxes_copy = hypre_BoxArrayDuplicate(neighbor_boxes);
882 for (i=0; i< size; i++)
883 {
884 hypre_CopyBox(hypre_BoxArrayBox(n_boxes_copy, order_index[i]),
885 hypre_BoxArrayBox(neighbor_boxes,i ));
886 }
887 hypre_TFree(order_index);
888 hypre_BoxArrayDestroy(n_boxes_copy);
889
890 /*now delete (index) duplicates from boxes */
891 hypre_DeleteMultipleBoxes( neighbor_boxes, delete_array, index);
892 /*now delete (index) duplicates from ids and boxnums */
893 /* size = num of neighbor boxes */
894 if (index)
895 {
896 start = delete_array[0];
897 j = 0;
898 for (i = start; (i + j) < size; i++)
899 {
900 if (j < index)
901 {
902 while ((i+j) == delete_array[j]) /* see if deleting consec. items */
903 {
904 j++; /*increase the shift*/
905 if (j == index) break;
906 }
907 }
908 if ((i+j) < size) /* if deleteing the last item then no moving */
909 {
910 neighbor_boxnums[i] = neighbor_boxnums[i+j];
911 neighbor_proc_ids[i] = neighbor_proc_ids[i+j];
912 }
913
914 }
915 }
916
917 /* finished removing duplicates! */
918 /* new number of boxes is (size - index) */
919
920
921
922 hypre_TFree(delete_array);
923
924
925#if AP_STAT
926
927 printf("!!!!!*myid = %d: I have %d potential "
928 "neighbor boxes \n", myid, size-index);
929
930#endif
931
932#if 0
933
934/* now let's check the shortened list */
935
936 printf("!!!!!*myid = %d: I have %d potential neighbor boxes \n",
937 myid, size-index);
938 hypre_ForBoxI(i, neighbor_boxes)
939 {
940 box = hypre_BoxArrayBox(neighbor_boxes, i);
941 printf("!!!!!myid = %d: BOXES in my list of potential neighbors "
942 "(number %d): (%d, %d, %d) x (%d, %d, %d) , boxnum = %d,"
943 "proc %d\n",
944 myid, i,
945 hypre_BoxIMinX(box),
946 hypre_BoxIMinY(box),
947 hypre_BoxIMinZ(box),
948 hypre_BoxIMaxX(box),
949 hypre_BoxIMaxY(box),
950 hypre_BoxIMaxZ(box), neighbor_boxnums[i], neighbor_proc_ids[i]);
951 }
952#endif
953
954
955 /* create the neighbors structure! */
956 hypre_BoxNeighborsCreateWithAP(neighbor_boxes, neighbor_proc_ids,
957 neighbor_boxnums,
958 first_local, num_local_boxes, pshifts, &neighbors);
959
960
961
962
963
964
965 hypre_StructGridNeighbors(grid) = neighbors;
966
967
968 }
969
970
971#if NEIGH_PRINT
972 end_time = MPI_Wtime();
973 end_time = end_time - start_time;
974 printf("myid = %d Getting ready for neighbor assemble stuff time = %f sec.\n", myid, end_time);
975#endif
976
977
978/*---------------------------------------------
979 Step 6: Find your neighbors
980 -----------------------------------------------*/
981
982
983#if NEIGH_PRINT
984 start_time = MPI_Wtime();
985#endif
986
987 hypre_BoxNeighborsAssembleWithAP(neighbors, periodic, max_distance, prune);
988
989
990#if NEIGH_PRINT
991 end_time = MPI_Wtime();
992 end_time = end_time - start_time;
993 printf("myid = %d Box Neighbor assemble time = %f sec.\n", myid, end_time);
994#endif
995
996/*---------------------------------------------
997 Step 7: Expand to include ghosts
998-----------------------------------------------*/
999
1000
1001 numghost = hypre_StructGridNumGhost(grid) ;
1002 ghostsize = 0;
1003 ghostbox = hypre_BoxCreate();
1004 hypre_ForBoxI(i, local_boxes)
1005 {
1006 box = hypre_BoxArrayBox(local_boxes, i);
1007
1008 hypre_CopyBox(box, ghostbox);
1009 hypre_BoxExpand(ghostbox, numghost);
1010
1011 ghostsize += hypre_BoxVolume(ghostbox);
1012
1013 }
1014
1015 hypre_StructGridGhlocalSize(grid) = ghostsize;
1016 hypre_BoxDestroy(ghostbox);
1017
1018
1019
1020
1021 /* clean up */
1022
1023 /*neighbor_ids or neighbor_boxes, neighbor_boxnums, and
1024 neighbor_proc_ids
1025 get aliased and destroyed elsewhere (BoxNeighborDestroy) */
1026
1027
1028
1029
1030/*---------------------------------------------
1031 DEBUG
1032-----------------------------------------------*/
1033
1034#if DEBUG
1035{
1036 hypre_BoxNeighbors *neighbors;
1037 int *procs, *boxnums;
1038 int num_neighbors, i;
1039
1040 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
1041
1042 sprintf(filename, "zgrid.%05d", my_rank);
1043
1044 if ((file = fopen(filename, "a")) == NULL)
1045 {
1046 printf("Error: can't open output file %s\n", filename);
1047 exit(1);
1048 }
1049
1050 fprintf(file, "\n\n============================\n\n");
1051
1052 hypre_StructGridPrint(file, grid);
1053
1054 neighbors = hypre_StructGridNeighbors(grid);
1055 num_neighbors = hypre_BoxArraySize(hypre_BoxNeighborsBoxes(neighbors));
1056 procs = hypre_BoxNeighborsProcs(neighbors);
1057 boxnums = hypre_BoxNeighborsBoxnums(neighbors);
1058 fprintf(file, "num_neighbors = %d\n", num_neighbors);
1059 for (i = 0; i < num_neighbors; i++)
1060 {
1061 fprintf(file, "%d: (%d, %d)\n", i, procs[i], boxnums[i]);
1062 }
1063
1064 fflush(file);
1065 fclose(file);
1066}
1067#endif
1068
1069
1070 return ierr;
1071}
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081/******************************************************************************
1082 *
1083 * fillResponseStructAssembleAP
1084 *
1085 *****************************************************************************/
1086
1087int
1088hypre_FillResponseStructAssembleAP(void *p_recv_contact_buf,
1089 int contact_size, int contact_proc, void *ro,
1090 MPI_Comm comm, void **p_send_response_buf,
1091 int *response_message_size )
1092{
1093
1094
1095 int myid, i, index, d;
1096 int size, num_objects;
1097 int *send_response_buf = (int *) *p_send_response_buf;
1098 int *proc_ids;
1099 int *boxnums;
1100
1101
1102 hypre_BoxArray *box_array;
1103 hypre_Box *box;
1104
1105
1106 hypre_DataExchangeResponse *response_obj = ro;
1107 hypre_StructAssumedPart *assumed_part = response_obj->data1;
1108
1109 int overhead = response_obj->send_response_overhead;
1110
1111 /*initialize stuff */
1112 MPI_Comm_rank(comm, &myid );
1113
1114
1115 proc_ids = hypre_StructAssumedPartMyPartitionProcIds(assumed_part);
1116 box_array = hypre_StructAssumedPartMyPartitionBoxes(assumed_part);
1117 boxnums = hypre_StructAssumedPartMyPartitionBoxnums(assumed_part);
1118
1119 /*we need to send back the list of all of the boxes in our
1120 assumed partition along with the corresponding processor id */
1121
1122 /*how many boxes and ids do we have?*/
1123 num_objects = hypre_StructAssumedPartMyPartitionIdsSize(assumed_part);
1124 /* num_objects is then how much we need to send*/
1125
1126
1127 /*check storage in send_buf for adding the information */
1128 /* note: we are returning objects that are 8 ints in size */
1129
1130 if ( response_obj->send_response_storage < num_objects )
1131 {
1132 response_obj->send_response_storage = hypre_max(num_objects, 10);
1133 size = 8*(response_obj->send_response_storage + overhead);
1134 send_response_buf = hypre_TReAlloc( send_response_buf, int,
1135 size);
1136 *p_send_response_buf = send_response_buf; /* needed when using ReAlloc */
1137 }
1138
1139 /* populate send_response_buf */
1140 index = 0;
1141 for (i = 0; i< num_objects; i++)
1142 {
1143 /* processor id */
1144 send_response_buf[index++] = proc_ids[i];
1145 send_response_buf[index++] = boxnums[i];
1146 box = hypre_BoxArrayBox(box_array, i);
1147
1148 for (d=0; d< 3; d++)
1149 {
1150 send_response_buf[index++] = hypre_BoxIMinD(box, d);
1151 send_response_buf[index++] = hypre_BoxIMaxD(box, d);
1152 }
1153 }
1154
1155 /* return variable */
1156 *response_message_size = num_objects;
1157 *p_send_response_buf = send_response_buf;
1158
1159 return(0);
1160
1161
1162}
1163
1164
1165/*--------------------------------------------------------------------------
1166 * hypre_StructGridSetIDs
1167 *--------------------------------------------------------------------------*/
1168
1169int
1170hypre_StructGridSetIDs( hypre_StructGrid *grid,
1171 int *ids )
1172{
1173 int ierr = 0;
1174
1175 hypre_TFree(hypre_StructGridIDs(grid));
1176 hypre_StructGridIDs(grid) = ids;
1177
1178 return ierr;
1179}
Note: See TracBrowser for help on using the repository browser.