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 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 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 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 @*/ 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