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 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 PetscValidPointer(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(0); 61 } 62 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 84 vtx = *iter; 85 /* Get adjacency information for current vertex - i.e., all elements of dimension (dim) that connects 86 to the current vertex. We can then decipher if a vertex is ghosted or not and compute the 87 non-zero pattern accordingly. */ 88 adjs.clear(); 89 if (dmmoab->hlevel && (dmmoab->pcomm->size() == 1)) { 90 merr = dmmoab->hierarchy->get_adjacencies(vtx, dmmoab->dim, adjs); MBERRNM(merr); 91 } else { 92 merr = dmmoab->mbiface->get_adjacencies(&vtx, 1, dmmoab->dim, true, adjs, moab::Interface::UNION); MBERRNM(merr); 93 } 94 95 /* reset counters */ 96 n_nnz = n_onz = 0; 97 found.clear(); 98 99 /* loop over vertices and update the number of connectivity */ 100 for (unsigned jter = 0; jter < adjs.size(); ++jter) { 101 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); MBERRNM(merr); 106 107 /* loop over each element connected to the adjacent vertex and update as needed */ 108 for (i = 0; i < vpere; ++i) { 109 /* find the truly user-expected layer of ghosted entities to decipher NNZ pattern */ 110 if (connect[i] == vtx || found.find(connect[i]) != found.end()) continue; /* make sure we don't double count shared vertices */ 111 if (dmmoab->vghost->find(connect[i]) != dmmoab->vghost->end()) n_onz++; /* update out-of-proc onz */ 112 else n_nnz++; /* else local vertex */ 113 found.insert(connect[i]); 114 } 115 } 116 storage.clear(); 117 118 if (isbaij) { 119 nnz[ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 120 if (onz) { 121 onz[ivtx] = n_onz; /* add ghost non-owned nodes */ 122 } 123 } else { /* AIJ matrices */ 124 if (!isinterlaced) { 125 for (f = 0; f < nfields; f++) { 126 nnz[f * nloc + ivtx] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 127 if (onz) 128 onz[f * nloc + ivtx] = n_onz; /* add ghost non-owned nodes */ 129 } 130 } else { 131 for (f = 0; f < nfields; f++) { 132 nnz[nfields * ivtx + f] = n_nnz; /* leave out self to avoid repeats -> node shared by multiple elements */ 133 if (onz) 134 onz[nfields * ivtx + f] = n_onz; /* add ghost non-owned nodes */ 135 } 136 } 137 } 138 } 139 140 for (i = 0; i < nlsiz; i++) 141 nnz[i] += 1; /* self count the node */ 142 143 for (ivtx = 0; ivtx < nloc; ivtx++) { 144 if (!isbaij) { 145 for (ibs = 0; ibs < nfields; ibs++) { 146 if (dmmoab->dfill) { /* first address the diagonal block */ 147 /* just add up the ints -- easier/faster rather than branching based on "1" */ 148 for (jbs = 0, inbsize = 0; jbs < nfields; jbs++) 149 inbsize += dmmoab->dfill[ibs * nfields + jbs]; 150 } 151 else inbsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 152 if (isinterlaced) nnz[ivtx * nfields + ibs] *= inbsize; 153 else nnz[ibs * nloc + ivtx] *= inbsize; 154 155 if (onz) { 156 if (dmmoab->ofill) { /* next address the off-diagonal block */ 157 /* just add up the ints -- easier/faster rather than branching based on "1" */ 158 for (jbs = 0, iobsize = 0; jbs < nfields; jbs++) 159 iobsize += dmmoab->dfill[ibs * nfields + jbs]; 160 } 161 else iobsize = nfields; /* dense coupling since user didn't specify the component fill explicitly */ 162 if (isinterlaced) onz[ivtx * nfields + ibs] *= iobsize; 163 else onz[ibs * nloc + ivtx] *= iobsize; 164 } 165 } 166 } else { 167 /* check if we got overzealous in our nnz and onz computations */ 168 nnz[ivtx] = (nnz[ivtx] > dmmoab->nloc ? dmmoab->nloc : nnz[ivtx]); 169 if (onz) onz[ivtx] = (onz[ivtx] > dmmoab->nloc ? dmmoab->nloc : onz[ivtx]); 170 } 171 } 172 /* update innz and ionz based on local maxima */ 173 if (innz || ionz) { 174 if (innz) *innz = 0; 175 if (ionz) *ionz = 0; 176 for (i = 0; i < nlsiz; i++) { 177 if (innz && (nnz[i] > *innz)) *innz = nnz[i]; 178 if ((ionz && onz) && (onz[i] > *ionz)) *ionz = onz[i]; 179 } 180 } 181 PetscFunctionReturn(0); 182 } 183 184 static PetscErrorCode DMMoabSetBlockFills_Private(PetscInt w, const PetscInt *fill, PetscInt **rfill) 185 { 186 PetscInt i, j, *ifill; 187 188 PetscFunctionBegin; 189 if (!fill) PetscFunctionReturn(0); 190 PetscCall(PetscMalloc1(w * w, &ifill)); 191 for (i = 0; i < w; i++) { 192 for (j = 0; j < w; j++) 193 ifill[i * w + j] = fill[i * w + j]; 194 } 195 196 *rfill = ifill; 197 PetscFunctionReturn(0); 198 } 199 200 /*@C 201 DMMoabSetBlockFills - Sets the fill pattern in each block for a multi-component problem 202 of the matrix returned by DMCreateMatrix(). 203 204 Logically Collective on da 205 206 Input Parameters: 207 + dm - the DMMoab object 208 . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block) 209 - ofill - the fill pattern in the off-diagonal blocks 210 211 Level: developer 212 213 Notes: 214 This only makes sense when you are doing multicomponent problems but using the 215 MPIAIJ matrix format 216 217 The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 218 representing coupling and 0 entries for missing coupling. For example 219 $ dfill[9] = {1, 0, 0, 220 $ 1, 1, 0, 221 $ 0, 1, 1} 222 means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 223 itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 224 diagonal block). 225 226 DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 227 can be represented in the dfill, ofill format 228 229 Contributed by Glenn Hammond 230 231 .seealso `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()` 232 233 @*/ 234 PetscErrorCode DMMoabSetBlockFills(DM dm, const PetscInt *dfill, const PetscInt *ofill) 235 { 236 DM_Moab *dmmoab = (DM_Moab*)dm->data; 237 238 PetscFunctionBegin; 239 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 240 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, dfill, &dmmoab->dfill)); 241 PetscCall(DMMoabSetBlockFills_Private(dmmoab->numFields, ofill, &dmmoab->ofill)); 242 PetscFunctionReturn(0); 243 } 244