1 #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/
2 #include <petsc/private/vecimpl.h>
3
4 #include <petscdmmoab.h>
5 #include <MBTagConventions.hpp>
6 #include <moab/NestedRefine.hpp>
7
8 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscBool);
9
DMCreateMatrix_Moab(DM dm,Mat * J)10 PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J)
11 {
12 PetscInt innz = 0, ionz = 0, nlsiz;
13 DM_Moab *dmmoab = (DM_Moab *)dm->data;
14 PetscInt *nnz = 0, *onz = 0;
15 char *tmp = 0;
16 Mat A;
17 MatType mtype;
18
19 PetscFunctionBegin;
20 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
21 PetscAssertPointer(J, 3);
22
23 /* next, need to allocate the non-zero arrays to enable pre-allocation */
24 mtype = dm->mattype;
25 PetscCall(PetscStrstr(mtype, MATBAIJ, &tmp));
26 nlsiz = (tmp ? dmmoab->nloc : dmmoab->nloc * dmmoab->numFields);
27
28 /* allocate the nnz, onz arrays based on block size and local nodes */
29 PetscCall(PetscCalloc2(nlsiz, &nnz, nlsiz, &onz));
30
31 /* compute the nonzero pattern based on MOAB connectivity data for local elements */
32 PetscCall(DMMoab_Compute_NNZ_From_Connectivity(dm, &innz, nnz, &ionz, onz, tmp ? PETSC_TRUE : PETSC_FALSE));
33
34 /* create the Matrix and set its type as specified by user */
35 PetscCall(MatCreate(((PetscObject)dm)->comm, &A));
36 PetscCall(MatSetSizes(A, dmmoab->nloc * dmmoab->numFields, dmmoab->nloc * dmmoab->numFields, PETSC_DETERMINE, PETSC_DETERMINE));
37 PetscCall(MatSetType(A, mtype));
38 PetscCall(MatSetBlockSize(A, dmmoab->bs));
39 PetscCall(MatSetDM(A, dm)); /* set DM reference */
40 PetscCall(MatSetFromOptions(A));
41
42 PetscCheck(dmmoab->ltog_map, ((PetscObject)dm)->comm, PETSC_ERR_ORDER, "Cannot create a DMMoab Mat without calling DMSetUp first.");
43 PetscCall(MatSetLocalToGlobalMapping(A, dmmoab->ltog_map, dmmoab->ltog_map));
44
45 /* set preallocation based on different supported Mat types */
46 PetscCall(MatSeqAIJSetPreallocation(A, innz, nnz));
47 PetscCall(MatMPIAIJSetPreallocation(A, innz, nnz, ionz, onz));
48 PetscCall(MatSeqBAIJSetPreallocation(A, dmmoab->bs, innz, nnz));
49 PetscCall(MatMPIBAIJSetPreallocation(A, dmmoab->bs, innz, nnz, ionz, onz));
50
51 /* clean up temporary memory */
52 PetscCall(PetscFree2(nnz, onz));
53
54 /* set up internal matrix data-structures */
55 PetscCall(MatSetUp(A));
56
57 /* MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); */
58
59 *J = A;
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
DMMoab_Compute_NNZ_From_Connectivity(DM dm,PetscInt * innz,PetscInt * nnz,PetscInt * ionz,PetscInt * onz,PetscBool isbaij)63 PETSC_EXTERN PetscErrorCode DMMoab_Compute_NNZ_From_Connectivity(DM dm, PetscInt *innz, PetscInt *nnz, PetscInt *ionz, PetscInt *onz, PetscBool isbaij)
64 {
65 PetscInt i, f, nloc, vpere, bs, n_nnz, n_onz, ivtx = 0;
66 PetscInt ibs, jbs, inbsize, iobsize, nfields, nlsiz;
67 DM_Moab *dmmoab = (DM_Moab *)dm->data;
68 moab::Range found;
69 std::vector<moab::EntityHandle> adjs, storage;
70 PetscBool isinterlaced;
71 moab::EntityHandle vtx;
72 moab::ErrorCode merr;
73
74 PetscFunctionBegin;
75 bs = dmmoab->bs;
76 nloc = dmmoab->nloc;
77 nfields = dmmoab->numFields;
78 isinterlaced = (isbaij || bs == nfields ? PETSC_TRUE : PETSC_FALSE);
79 nlsiz = (isinterlaced ? nloc : nloc * nfields);
80
81 /* loop over the locally owned vertices and figure out the NNZ pattern using connectivity information */
82 for (moab::Range::const_iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, ivtx++) {
83 vtx = *iter;
84 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects
85 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the
86 non-zero pattern accordingly. */
87 adjs.clear();
88 if (dmmoab->hlevel && (dmmoab->pcomm->size() == 1)) {
89 merr = dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs);
90 MBERRNM(merr);
91 } else {
92 merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::UNION);
93 MBERRNM(merr);
94 }
95
96 /* reset counters */
97 n_nnz = n_onz = 0;
98 found.clear();
99
100 /* loop over vertices and update the number of connectivity */
101 for (unsigned jter = 0; jter < adjs.size(); ++jter) {
102 /* Get connectivity information in canonical ordering for the local element */
103 const moab::EntityHandle *connect;
104 std::vector<moab::EntityHandle> cconnect;
105 merr = dmmoab->mbiface->get_connectivity(adjs[jter], connect, vpere, false, &storage);
106 MBERRNM(merr);
107
108 /* loop over each element connected to the adjacent vertex and update as needed */
109 for (i = 0; i < vpere; ++i) {
110 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */
111 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */
112 if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++; /* update out-of-proc onz */
113 else n_nnz++; /* else local vertex */
114 found.insert(connect[i]);
115 }
116 }
117 storage.clear();
118
119 if (isbaij) {
120 nnz[ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
121 if (onz) onz[ivtx] = n_onz; /* add ghost non-owned nodes */
122 } else { /* AIJ matrices */
123 if (!isinterlaced) {
124 for (f = 0; f < nfields; f++) {
125 nnz[f * nloc + ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
126 if (onz) onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */
127 }
128 } else {
129 for (f = 0; f < nfields; f++) {
130 nnz[nfields * ivtx + f] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */
131 if (onz) onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */
132 }
133 }
134 }
135 }
136
137 for (i = 0; i < nlsiz; i++) nnz[i] += 1; /* self count the node */
138
139 for (ivtx = 0; ivtx < nloc; ivtx++) {
140 if (!isbaij) {
141 for (ibs = 0; ibs < nfields; ibs++) {
142 if (dmmoab->dfill) { /* first address the diagonal block */
143 /* just add up the ints -- easier/faster rather than branching based on "1" */
144 for (jbs = 0, inbsize = 0; jbs < nfields; jbs++) inbsize += dmmoab->dfill[ibs * nfields + jbs];
145 } else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
146 if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize;
147 else nnz[ibs * nloc + ivtx] *= inbsize;
148
149 if (onz) {
150 if (dmmoab->ofill) { /* next address the off-diagonal block */
151 /* just add up the ints -- easier/faster rather than branching based on "1" */
152 for (jbs = 0, iobsize = 0; jbs < nfields; jbs++) iobsize += dmmoab->dfill[ibs * nfields + jbs];
153 } else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */
154 if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize;
155 else onz[ibs * nloc + ivtx] *= iobsize;
156 }
157 }
158 } else {
159 /* check if we got overzealous in our nnz and onz computations */
160 nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]);
161 if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]);
162 }
163 }
164 /* update innz and ionz based on local maxima */
165 if (innz || ionz) {
166 if (innz) *innz = 0;
167 if (ionz) *ionz = 0;
168 for (i = 0; i < nlsiz; i++) {
169 if (innz && nnz[i] > *innz) *innz = nnz[i];
170 if (ionz && onz && onz[i] > *ionz) *ionz = onz[i];
171 }
172 }
173 PetscFunctionReturn(PETSC_SUCCESS);
174 }
175
DMMoabSetBlockFills_Private(PetscInt w,const PetscInt * fill,PetscInt ** rfill)176 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill)
177 {
178 PetscInt i, j, *ifill;
179
180 PetscFunctionBegin;
181 if (!fill) PetscFunctionReturn(PETSC_SUCCESS);
182 PetscCall(PetscMalloc1(w * w, &ifill));
183 for (i = 0; i < w; i++) {
184 for (j = 0; j < w; j++) ifill[i * w + j] = fill[i * w + j];
185 }
186
187 *rfill = ifill;
188 PetscFunctionReturn(PETSC_SUCCESS);
189 }
190
191 /*@C
192 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem
193 of the matrix returned by `DMCreateMatrix()`.
194
195 Logically Collective
196
197 Input Parameters:
198 + dm - the `DMMOAB` object
199 . dfill - the fill pattern in the diagonal block (may be `NULL`, means use dense block)
200 - ofill - the fill pattern in the off-diagonal blocks
201
202 Level: developer
203
204 Notes:
205 This only makes sense when you are doing multicomponent problems but using the
206 MPIAIJ matrix format
207
208 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
209 representing coupling and 0 entries for missing coupling. For example
210 .vb
211 dfill[9] = {1, 0, 0,
212 1, 1, 0,
213 0, 1, 1}
214 .ve
215 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
216 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
217 diagonal block).
218
219 `DMDASetGetMatrix()` allows you to provide general code for those more complicated nonzero patterns then
220 can be represented in the dfill, ofill format
221
222 Contributed by\:
223 Glenn Hammond
224
225 .seealso: `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()`
226 @*/
DMMoabSetBlockFills(DM dm,const PetscInt * dfill,const PetscInt * ofill)227 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill)
228 {
229 DM_Moab *dmmoab = (DM_Moab *)dm->data;
230
231 PetscFunctionBegin;
232 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
233 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill));
234 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill));
235 PetscFunctionReturn(PETSC_SUCCESS);
236 }
237