PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
scatters.cc
Go to the documentation of this file.
1 /*!
2  * This code sets up scatters of a distributed Vec to a given rank to process data using
3  * serial code.
4  *
5  * It is *not* used at the moment.
6  */
7 #if 0
8 struct VecAndScatter {
9  VecScatter scatter;
10  Vec v;
11 };
12 
13 /*!
14  * Allocate the scatter from a part of a parallel Vec to a target rank.
15  *
16  * The caller is responsible for de-allocating both the scatter and the target Vec.
17  */
18 static VecAndScatter scatter_part(Vec v_in, int start, int length, int target_rank) {
19  PetscErrorCode ierr;
20  int rank;
21  VecAndScatter result;
22  IS is;
23 
24  MPI_Comm_rank(PetscObjectComm((PetscObject)v_in), &rank);
25 
26  if (rank != target_rank) {
27  length = 0;
28  }
29 
30  ierr = VecCreateSeq(PETSC_COMM_SELF, length, &result.v);
31  PISM_CHK(ierr, "VecCreateSeq");
32 
33  ierr = ISCreateStride(PETSC_COMM_SELF, length, start, 1, &is);
34  PISM_CHK(ierr, "ISCreateStride");
35 
36  ierr = VecScatterCreate(v_in, is, result.v, NULL, &result.scatter);
37  PISM_CHK(ierr, "VecScatterCreate");
38 
39  ierr = ISDestroy(&is);
40  PISM_CHK(ierr, "ISDestroy");
41 
42  return result;
43 }
44 
45 /*!
46  * Allocate a natural Vec for a given DM.
47  *
48  * The caller is responsible for de-allocating the Vec returned by this function.
49  */
50 static Vec get_natural_work(DM dm) {
51  PetscErrorCode ierr;
52  Vec result;
53 
54  ierr = PetscObjectQuery((PetscObject)dm, "natural_work", (PetscObject*)&result);
55  PISM_CHK(ierr, "PetscObjectQuery");
56 
57  if (result == NULL) {
58  Vec v = NULL;
59  ierr = DMDACreateNaturalVector(dm, &v);
60  PISM_CHK(ierr, "DMDACreateNaturalVector");
61 
62  ierr = PetscObjectCompose((PetscObject)dm, "natural_work", (PetscObject)(v));
63  PISM_CHK(ierr, "PetscObjectCompose");
64 
65  result = v;
66 
67  ierr = VecDestroy(&v);
68  PISM_CHK(ierr, "VecDestroy");
69  }
70 
71  return result;
72 }
73 
74 /*!
75  * Given a DM, allocate a rank 0 target Vec that can be used to gather a part of a
76  * "global" Vec on rank 0. Arguments "start" and "length" define the part in question.
77  *
78  * The caller is responsible for de-allocating the Vec returned by this function.
79  */
80 static Vec proc0_copy(DM dm, int start, int length) {
81  Vec v_proc0 = NULL;
82  PetscErrorCode ierr = 0;
83 
84  ierr = PetscObjectQuery((PetscObject)dm, "v_proc0", (PetscObject*)&v_proc0);
85  PISM_CHK(ierr, "PetscObjectQuery");
86  ;
87  if (v_proc0 == NULL) {
88 
89  // natural_work will be destroyed at the end of scope, but it will
90  // only decrement the reference counter incremented by
91  // PetscObjectCompose below.
92  auto natural_work = get_natural_work(dm);
93 
94  // scatter_to_zero will be destroyed at the end of scope, but it
95  // will only decrement the reference counter incremented by
96  // PetscObjectCompose below.
97  auto vs = scatter_part(natural_work, start, length, 0);
98 
99  // this increments the reference counter of scatter_to_zero
100  ierr = PetscObjectCompose((PetscObject)dm, "scatter_to_zero",
101  (PetscObject)(vs.scatter));
102  PISM_CHK(ierr, "PetscObjectCompose");
103 
104  // this increments the reference counter of v_proc0
105  ierr = PetscObjectCompose((PetscObject)dm, "v_proc0",
106  (PetscObject)vs.v);
107  PISM_CHK(ierr, "PetscObjectCompose");
108 
109  VecScatterDestroy(&vs.scatter);
110 
111  // We DO NOT call VecDestroy(v_proc0): the petsc::Vec wrapper will
112  // take care of this.
113  return vs.v;
114  }
115  return v_proc0;
116 }
117 #endif
#define PISM_CHK(errcode, name)