1 2 /* 3 Code for interpolating between grids represented by DMDAs 4 */ 5 6 /* 7 For linear elements there are two branches of code to compute the interpolation. They should compute the same results but may not. The "new version" does 8 not work for periodic domains, the old does. Change NEWVERSION to 1 to compile in the new version. Eventually when we are sure the two produce identical results 9 we will remove/merge the new version. Based on current tests, these both produce the same results. We are leaving NEWVERSION for now in the code since some 10 consider it cleaner, but old version is turned on since it handles periodic case. 11 */ 12 #define NEWVERSION 0 13 14 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 15 16 /* 17 Since the interpolation uses MATMAIJ for dof > 0 we convert request for non-MATAIJ baseded matrices to MATAIJ. 18 This is a bit of a hack, the reason for it is partially because -dm_mat_type defines the 19 matrix type for both the operator matrices and the interpolation matrices so that users 20 can select matrix types of base MATAIJ for accelerators 21 */ 22 static PetscErrorCode ConvertToAIJ(MatType intype, MatType *outtype) { 23 PetscInt i; 24 char const *types[3] = {MATAIJ, MATSEQAIJ, MATMPIAIJ}; 25 PetscBool flg; 26 27 PetscFunctionBegin; 28 *outtype = MATAIJ; 29 for (i = 0; i < 3; i++) { 30 PetscCall(PetscStrbeginswith(intype, types[i], &flg)); 31 if (flg) { 32 *outtype = intype; 33 PetscFunctionReturn(0); 34 } 35 } 36 PetscFunctionReturn(0); 37 } 38 39 PetscErrorCode DMCreateInterpolation_DA_1D_Q1(DM dac, DM daf, Mat *A) { 40 PetscInt i, i_start, m_f, Mx; 41 const PetscInt *idx_f, *idx_c; 42 PetscInt m_ghost, m_ghost_c; 43 PetscInt row, col, i_start_ghost, mx, m_c, nc, ratio; 44 PetscInt i_c, i_start_c, i_start_ghost_c, cols[2], dof; 45 PetscScalar v[2], x; 46 Mat mat; 47 DMBoundaryType bx; 48 ISLocalToGlobalMapping ltog_f, ltog_c; 49 MatType mattype; 50 51 PetscFunctionBegin; 52 PetscCall(DMDAGetInfo(dac, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bx, NULL, NULL, NULL)); 53 PetscCall(DMDAGetInfo(daf, NULL, &mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 54 if (bx == DM_BOUNDARY_PERIODIC) { 55 ratio = mx / Mx; 56 PetscCheck(ratio * Mx == mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mx/Mx must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 57 } else { 58 ratio = (mx - 1) / (Mx - 1); 59 PetscCheck(ratio * (Mx - 1) == mx - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 60 } 61 62 PetscCall(DMDAGetCorners(daf, &i_start, NULL, NULL, &m_f, NULL, NULL)); 63 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, NULL, NULL, &m_ghost, NULL, NULL)); 64 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 65 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 66 67 PetscCall(DMDAGetCorners(dac, &i_start_c, NULL, NULL, &m_c, NULL, NULL)); 68 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, NULL, NULL, &m_ghost_c, NULL, NULL)); 69 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 70 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 71 72 /* create interpolation matrix */ 73 PetscCall(MatCreate(PetscObjectComm((PetscObject)dac), &mat)); 74 #if defined(PETSC_HAVE_CUDA) 75 /* 76 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU 77 we don't want the original unconverted matrix copied to the GPU 78 */ 79 if (dof > 1) { PetscCall(MatBindToCPU(mat, PETSC_TRUE)); } 80 #endif 81 PetscCall(MatSetSizes(mat, m_f, m_c, mx, Mx)); 82 PetscCall(ConvertToAIJ(dac->mattype, &mattype)); 83 PetscCall(MatSetType(mat, mattype)); 84 PetscCall(MatSeqAIJSetPreallocation(mat, 2, NULL)); 85 PetscCall(MatMPIAIJSetPreallocation(mat, 2, NULL, 1, NULL)); 86 87 /* loop over local fine grid nodes setting interpolation for those*/ 88 if (!NEWVERSION) { 89 for (i = i_start; i < i_start + m_f; i++) { 90 /* convert to local "natural" numbering and then to PETSc global numbering */ 91 row = idx_f[i - i_start_ghost]; 92 93 i_c = (i / ratio); /* coarse grid node to left of fine grid node */ 94 PetscCheck(i_c >= i_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 95 i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT, 96 i_start, i_c, i_start_ghost_c); 97 98 /* 99 Only include those interpolation points that are truly 100 nonzero. Note this is very important for final grid lines 101 in x direction; since they have no right neighbor 102 */ 103 x = ((PetscReal)(i - i_c * ratio)) / ((PetscReal)ratio); 104 nc = 0; 105 /* one left and below; or we are right on it */ 106 col = (i_c - i_start_ghost_c); 107 cols[nc] = idx_c[col]; 108 v[nc++] = -x + 1.0; 109 /* one right? */ 110 if (i_c * ratio != i) { 111 cols[nc] = idx_c[col + 1]; 112 v[nc++] = x; 113 } 114 PetscCall(MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES)); 115 } 116 117 } else { 118 PetscScalar *xi; 119 PetscInt li, nxi, n; 120 PetscScalar Ni[2]; 121 122 /* compute local coordinate arrays */ 123 nxi = ratio + 1; 124 PetscCall(PetscMalloc1(nxi, &xi)); 125 for (li = 0; li < nxi; li++) { xi[li] = -1.0 + (PetscScalar)li * (2.0 / (PetscScalar)(nxi - 1)); } 126 127 for (i = i_start; i < i_start + m_f; i++) { 128 /* convert to local "natural" numbering and then to PETSc global numbering */ 129 row = idx_f[(i - i_start_ghost)]; 130 131 i_c = (i / ratio); /* coarse grid node to left of fine grid node */ 132 PetscCheck(i_c >= i_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 133 i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT, 134 i_start, i_c, i_start_ghost_c); 135 136 /* remainders */ 137 li = i - ratio * (i / ratio); 138 if (i == mx - 1) li = nxi - 1; 139 140 /* corners */ 141 col = (i_c - i_start_ghost_c); 142 cols[0] = idx_c[col]; 143 Ni[0] = 1.0; 144 if ((li == 0) || (li == nxi - 1)) { 145 PetscCall(MatSetValue(mat, row, cols[0], Ni[0], INSERT_VALUES)); 146 continue; 147 } 148 149 /* edges + interior */ 150 /* remainders */ 151 if (i == mx - 1) i_c--; 152 153 col = (i_c - i_start_ghost_c); 154 cols[0] = idx_c[col]; /* one left and below; or we are right on it */ 155 cols[1] = idx_c[col + 1]; 156 157 Ni[0] = 0.5 * (1.0 - xi[li]); 158 Ni[1] = 0.5 * (1.0 + xi[li]); 159 for (n = 0; n < 2; n++) { 160 if (PetscAbsScalar(Ni[n]) < 1.0e-32) cols[n] = -1; 161 } 162 PetscCall(MatSetValues(mat, 1, &row, 2, cols, Ni, INSERT_VALUES)); 163 } 164 PetscCall(PetscFree(xi)); 165 } 166 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 167 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 168 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 169 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 170 PetscCall(MatCreateMAIJ(mat, dof, A)); 171 PetscCall(MatDestroy(&mat)); 172 PetscFunctionReturn(0); 173 } 174 175 PetscErrorCode DMCreateInterpolation_DA_1D_Q0(DM dac, DM daf, Mat *A) { 176 PetscInt i, i_start, m_f, Mx; 177 const PetscInt *idx_f, *idx_c; 178 ISLocalToGlobalMapping ltog_f, ltog_c; 179 PetscInt m_ghost, m_ghost_c; 180 PetscInt row, col, i_start_ghost, mx, m_c, nc, ratio; 181 PetscInt i_c, i_start_c, i_start_ghost_c, cols[2], dof; 182 PetscScalar v[2], x; 183 Mat mat; 184 DMBoundaryType bx; 185 MatType mattype; 186 187 PetscFunctionBegin; 188 PetscCall(DMDAGetInfo(dac, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bx, NULL, NULL, NULL)); 189 PetscCall(DMDAGetInfo(daf, NULL, &mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 190 if (bx == DM_BOUNDARY_PERIODIC) { 191 PetscCheck(Mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be positive", Mx); 192 ratio = mx / Mx; 193 PetscCheck(ratio * Mx == mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mx/Mx must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 194 } else { 195 PetscCheck(Mx >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be greater than 1", Mx); 196 ratio = (mx - 1) / (Mx - 1); 197 PetscCheck(ratio * (Mx - 1) == mx - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 198 } 199 200 PetscCall(DMDAGetCorners(daf, &i_start, NULL, NULL, &m_f, NULL, NULL)); 201 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, NULL, NULL, &m_ghost, NULL, NULL)); 202 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 203 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 204 205 PetscCall(DMDAGetCorners(dac, &i_start_c, NULL, NULL, &m_c, NULL, NULL)); 206 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, NULL, NULL, &m_ghost_c, NULL, NULL)); 207 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 208 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 209 210 /* create interpolation matrix */ 211 PetscCall(MatCreate(PetscObjectComm((PetscObject)dac), &mat)); 212 #if defined(PETSC_HAVE_CUDA) 213 /* 214 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU 215 we don't want the original unconverted matrix copied to the GPU 216 */ 217 if (dof > 1) { PetscCall(MatBindToCPU(mat, PETSC_TRUE)); } 218 #endif 219 PetscCall(MatSetSizes(mat, m_f, m_c, mx, Mx)); 220 PetscCall(ConvertToAIJ(dac->mattype, &mattype)); 221 PetscCall(MatSetType(mat, mattype)); 222 PetscCall(MatSeqAIJSetPreallocation(mat, 2, NULL)); 223 PetscCall(MatMPIAIJSetPreallocation(mat, 2, NULL, 0, NULL)); 224 225 /* loop over local fine grid nodes setting interpolation for those*/ 226 for (i = i_start; i < i_start + m_f; i++) { 227 /* convert to local "natural" numbering and then to PETSc global numbering */ 228 row = idx_f[(i - i_start_ghost)]; 229 230 i_c = (i / ratio); /* coarse grid node to left of fine grid node */ 231 232 /* 233 Only include those interpolation points that are truly 234 nonzero. Note this is very important for final grid lines 235 in x direction; since they have no right neighbor 236 */ 237 x = ((PetscReal)(i - i_c * ratio)) / ((PetscReal)ratio); 238 nc = 0; 239 /* one left and below; or we are right on it */ 240 col = (i_c - i_start_ghost_c); 241 cols[nc] = idx_c[col]; 242 v[nc++] = -x + 1.0; 243 /* one right? */ 244 if (i_c * ratio != i) { 245 cols[nc] = idx_c[col + 1]; 246 v[nc++] = x; 247 } 248 PetscCall(MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES)); 249 } 250 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 251 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 252 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 253 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 254 PetscCall(MatCreateMAIJ(mat, dof, A)); 255 PetscCall(MatDestroy(&mat)); 256 PetscCall(PetscLogFlops(5.0 * m_f)); 257 PetscFunctionReturn(0); 258 } 259 260 PetscErrorCode DMCreateInterpolation_DA_2D_Q1(DM dac, DM daf, Mat *A) { 261 PetscInt i, j, i_start, j_start, m_f, n_f, Mx, My, dof; 262 const PetscInt *idx_c, *idx_f; 263 ISLocalToGlobalMapping ltog_f, ltog_c; 264 PetscInt m_ghost, n_ghost, m_ghost_c, n_ghost_c, *dnz, *onz; 265 PetscInt row, col, i_start_ghost, j_start_ghost, cols[4], mx, m_c, my, nc, ratioi, ratioj; 266 PetscInt i_c, j_c, i_start_c, j_start_c, n_c, i_start_ghost_c, j_start_ghost_c, col_shift, col_scale; 267 PetscMPIInt size_c, size_f, rank_f; 268 PetscScalar v[4], x, y; 269 Mat mat; 270 DMBoundaryType bx, by; 271 MatType mattype; 272 273 PetscFunctionBegin; 274 PetscCall(DMDAGetInfo(dac, NULL, &Mx, &My, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL)); 275 PetscCall(DMDAGetInfo(daf, NULL, &mx, &my, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 276 if (bx == DM_BOUNDARY_PERIODIC) { 277 PetscCheck(Mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be positive", Mx); 278 ratioi = mx / Mx; 279 PetscCheck(ratioi * Mx == mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mx/Mx must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 280 } else { 281 PetscCheck(Mx >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be greater than 1", Mx); 282 ratioi = (mx - 1) / (Mx - 1); 283 PetscCheck(ratioi * (Mx - 1) == mx - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 284 } 285 if (by == DM_BOUNDARY_PERIODIC) { 286 PetscCheck(My, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of y coarse grid points %" PetscInt_FMT " must be positive", My); 287 ratioj = my / My; 288 PetscCheck(ratioj * My == my, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: my/My must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 289 } else { 290 PetscCheck(My >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of y coarse grid points %" PetscInt_FMT " must be greater than 1", My); 291 ratioj = (my - 1) / (My - 1); 292 PetscCheck(ratioj * (My - 1) == my - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (my - 1)/(My - 1) must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 293 } 294 295 PetscCall(DMDAGetCorners(daf, &i_start, &j_start, NULL, &m_f, &n_f, NULL)); 296 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, NULL, &m_ghost, &n_ghost, NULL)); 297 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 298 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 299 300 PetscCall(DMDAGetCorners(dac, &i_start_c, &j_start_c, NULL, &m_c, &n_c, NULL)); 301 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, NULL, &m_ghost_c, &n_ghost_c, NULL)); 302 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 303 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 304 305 /* 306 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 307 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 308 processors). It's effective length is hence 4 times its normal length, this is 309 why the col_scale is multiplied by the interpolation matrix column sizes. 310 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 311 copy of the coarse vector. A bit of a hack but you do better. 312 313 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 314 */ 315 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dac), &size_c)); 316 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)daf), &size_f)); 317 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)daf), &rank_f)); 318 col_scale = size_f / size_c; 319 col_shift = Mx * My * (rank_f / size_c); 320 321 MatPreallocateBegin(PetscObjectComm((PetscObject)daf), m_f * n_f, col_scale * m_c * n_c, dnz, onz); 322 for (j = j_start; j < j_start + n_f; j++) { 323 for (i = i_start; i < i_start + m_f; i++) { 324 /* convert to local "natural" numbering and then to PETSc global numbering */ 325 row = idx_f[(m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 326 327 i_c = (i / ratioi); /* coarse grid node to left of fine grid node */ 328 j_c = (j / ratioj); /* coarse grid node below fine grid node */ 329 330 PetscCheck(j_c >= j_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 331 j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT, 332 j_start, j_c, j_start_ghost_c); 333 PetscCheck(i_c >= i_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 334 i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT, 335 i_start, i_c, i_start_ghost_c); 336 337 /* 338 Only include those interpolation points that are truly 339 nonzero. Note this is very important for final grid lines 340 in x and y directions; since they have no right/top neighbors 341 */ 342 nc = 0; 343 /* one left and below; or we are right on it */ 344 col = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 345 cols[nc++] = col_shift + idx_c[col]; 346 /* one right and below */ 347 if (i_c * ratioi != i) cols[nc++] = col_shift + idx_c[col + 1]; 348 /* one left and above */ 349 if (j_c * ratioj != j) cols[nc++] = col_shift + idx_c[col + m_ghost_c]; 350 /* one right and above */ 351 if (i_c * ratioi != i && j_c * ratioj != j) cols[nc++] = col_shift + idx_c[col + (m_ghost_c + 1)]; 352 PetscCall(MatPreallocateSet(row, nc, cols, dnz, onz)); 353 } 354 } 355 PetscCall(MatCreate(PetscObjectComm((PetscObject)daf), &mat)); 356 #if defined(PETSC_HAVE_CUDA) 357 /* 358 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU 359 we don't want the original unconverted matrix copied to the GPU 360 */ 361 if (dof > 1) { PetscCall(MatBindToCPU(mat, PETSC_TRUE)); } 362 #endif 363 PetscCall(MatSetSizes(mat, m_f * n_f, col_scale * m_c * n_c, mx * my, col_scale * Mx * My)); 364 PetscCall(ConvertToAIJ(dac->mattype, &mattype)); 365 PetscCall(MatSetType(mat, mattype)); 366 PetscCall(MatSeqAIJSetPreallocation(mat, 0, dnz)); 367 PetscCall(MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz)); 368 MatPreallocateEnd(dnz, onz); 369 370 /* loop over local fine grid nodes setting interpolation for those*/ 371 if (!NEWVERSION) { 372 for (j = j_start; j < j_start + n_f; j++) { 373 for (i = i_start; i < i_start + m_f; i++) { 374 /* convert to local "natural" numbering and then to PETSc global numbering */ 375 row = idx_f[(m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 376 377 i_c = (i / ratioi); /* coarse grid node to left of fine grid node */ 378 j_c = (j / ratioj); /* coarse grid node below fine grid node */ 379 380 /* 381 Only include those interpolation points that are truly 382 nonzero. Note this is very important for final grid lines 383 in x and y directions; since they have no right/top neighbors 384 */ 385 x = ((PetscReal)(i - i_c * ratioi)) / ((PetscReal)ratioi); 386 y = ((PetscReal)(j - j_c * ratioj)) / ((PetscReal)ratioj); 387 388 nc = 0; 389 /* one left and below; or we are right on it */ 390 col = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 391 cols[nc] = col_shift + idx_c[col]; 392 v[nc++] = x * y - x - y + 1.0; 393 /* one right and below */ 394 if (i_c * ratioi != i) { 395 cols[nc] = col_shift + idx_c[col + 1]; 396 v[nc++] = -x * y + x; 397 } 398 /* one left and above */ 399 if (j_c * ratioj != j) { 400 cols[nc] = col_shift + idx_c[col + m_ghost_c]; 401 v[nc++] = -x * y + y; 402 } 403 /* one right and above */ 404 if (j_c * ratioj != j && i_c * ratioi != i) { 405 cols[nc] = col_shift + idx_c[col + (m_ghost_c + 1)]; 406 v[nc++] = x * y; 407 } 408 PetscCall(MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES)); 409 } 410 } 411 412 } else { 413 PetscScalar Ni[4]; 414 PetscScalar *xi, *eta; 415 PetscInt li, nxi, lj, neta; 416 417 /* compute local coordinate arrays */ 418 nxi = ratioi + 1; 419 neta = ratioj + 1; 420 PetscCall(PetscMalloc1(nxi, &xi)); 421 PetscCall(PetscMalloc1(neta, &eta)); 422 for (li = 0; li < nxi; li++) { xi[li] = -1.0 + (PetscScalar)li * (2.0 / (PetscScalar)(nxi - 1)); } 423 for (lj = 0; lj < neta; lj++) { eta[lj] = -1.0 + (PetscScalar)lj * (2.0 / (PetscScalar)(neta - 1)); } 424 425 /* loop over local fine grid nodes setting interpolation for those*/ 426 for (j = j_start; j < j_start + n_f; j++) { 427 for (i = i_start; i < i_start + m_f; i++) { 428 /* convert to local "natural" numbering and then to PETSc global numbering */ 429 row = idx_f[(m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 430 431 i_c = (i / ratioi); /* coarse grid node to left of fine grid node */ 432 j_c = (j / ratioj); /* coarse grid node below fine grid node */ 433 434 /* remainders */ 435 li = i - ratioi * (i / ratioi); 436 if (i == mx - 1) li = nxi - 1; 437 lj = j - ratioj * (j / ratioj); 438 if (j == my - 1) lj = neta - 1; 439 440 /* corners */ 441 col = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 442 cols[0] = col_shift + idx_c[col]; /* left, below */ 443 Ni[0] = 1.0; 444 if ((li == 0) || (li == nxi - 1)) { 445 if ((lj == 0) || (lj == neta - 1)) { 446 PetscCall(MatSetValue(mat, row, cols[0], Ni[0], INSERT_VALUES)); 447 continue; 448 } 449 } 450 451 /* edges + interior */ 452 /* remainders */ 453 if (i == mx - 1) i_c--; 454 if (j == my - 1) j_c--; 455 456 col = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 457 cols[0] = col_shift + idx_c[col]; /* left, below */ 458 cols[1] = col_shift + idx_c[col + 1]; /* right, below */ 459 cols[2] = col_shift + idx_c[col + m_ghost_c]; /* left, above */ 460 cols[3] = col_shift + idx_c[col + (m_ghost_c + 1)]; /* right, above */ 461 462 Ni[0] = 0.25 * (1.0 - xi[li]) * (1.0 - eta[lj]); 463 Ni[1] = 0.25 * (1.0 + xi[li]) * (1.0 - eta[lj]); 464 Ni[2] = 0.25 * (1.0 - xi[li]) * (1.0 + eta[lj]); 465 Ni[3] = 0.25 * (1.0 + xi[li]) * (1.0 + eta[lj]); 466 467 nc = 0; 468 if (PetscAbsScalar(Ni[0]) < 1.0e-32) cols[0] = -1; 469 if (PetscAbsScalar(Ni[1]) < 1.0e-32) cols[1] = -1; 470 if (PetscAbsScalar(Ni[2]) < 1.0e-32) cols[2] = -1; 471 if (PetscAbsScalar(Ni[3]) < 1.0e-32) cols[3] = -1; 472 473 PetscCall(MatSetValues(mat, 1, &row, 4, cols, Ni, INSERT_VALUES)); 474 } 475 } 476 PetscCall(PetscFree(xi)); 477 PetscCall(PetscFree(eta)); 478 } 479 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 480 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 481 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 482 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 483 PetscCall(MatCreateMAIJ(mat, dof, A)); 484 PetscCall(MatDestroy(&mat)); 485 PetscFunctionReturn(0); 486 } 487 488 /* 489 Contributed by Andrei Draganescu <aidraga@sandia.gov> 490 */ 491 PetscErrorCode DMCreateInterpolation_DA_2D_Q0(DM dac, DM daf, Mat *A) { 492 PetscInt i, j, i_start, j_start, m_f, n_f, Mx, My, dof; 493 const PetscInt *idx_c, *idx_f; 494 ISLocalToGlobalMapping ltog_f, ltog_c; 495 PetscInt m_ghost, n_ghost, m_ghost_c, n_ghost_c, *dnz, *onz; 496 PetscInt row, col, i_start_ghost, j_start_ghost, cols[4], mx, m_c, my, nc, ratioi, ratioj; 497 PetscInt i_c, j_c, i_start_c, j_start_c, n_c, i_start_ghost_c, j_start_ghost_c, col_shift, col_scale; 498 PetscMPIInt size_c, size_f, rank_f; 499 PetscScalar v[4]; 500 Mat mat; 501 DMBoundaryType bx, by; 502 MatType mattype; 503 504 PetscFunctionBegin; 505 PetscCall(DMDAGetInfo(dac, NULL, &Mx, &My, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL)); 506 PetscCall(DMDAGetInfo(daf, NULL, &mx, &my, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 507 PetscCheck(Mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be positive", Mx); 508 PetscCheck(My, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of y coarse grid points %" PetscInt_FMT " must be positive", My); 509 ratioi = mx / Mx; 510 ratioj = my / My; 511 PetscCheck(ratioi * Mx == mx, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Fine grid points must be multiple of coarse grid points in x"); 512 PetscCheck(ratioj * My == my, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Fine grid points must be multiple of coarse grid points in y"); 513 PetscCheck(ratioi == 2, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Coarsening factor in x must be 2"); 514 PetscCheck(ratioj == 2, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Coarsening factor in y must be 2"); 515 516 PetscCall(DMDAGetCorners(daf, &i_start, &j_start, NULL, &m_f, &n_f, NULL)); 517 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, NULL, &m_ghost, &n_ghost, NULL)); 518 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 519 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 520 521 PetscCall(DMDAGetCorners(dac, &i_start_c, &j_start_c, NULL, &m_c, &n_c, NULL)); 522 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, NULL, &m_ghost_c, &n_ghost_c, NULL)); 523 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 524 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 525 526 /* 527 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 528 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 529 processors). It's effective length is hence 4 times its normal length, this is 530 why the col_scale is multiplied by the interpolation matrix column sizes. 531 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 532 copy of the coarse vector. A bit of a hack but you do better. 533 534 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 535 */ 536 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dac), &size_c)); 537 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)daf), &size_f)); 538 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)daf), &rank_f)); 539 col_scale = size_f / size_c; 540 col_shift = Mx * My * (rank_f / size_c); 541 542 MatPreallocateBegin(PetscObjectComm((PetscObject)daf), m_f * n_f, col_scale * m_c * n_c, dnz, onz); 543 for (j = j_start; j < j_start + n_f; j++) { 544 for (i = i_start; i < i_start + m_f; i++) { 545 /* convert to local "natural" numbering and then to PETSc global numbering */ 546 row = idx_f[(m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 547 548 i_c = (i / ratioi); /* coarse grid node to left of fine grid node */ 549 j_c = (j / ratioj); /* coarse grid node below fine grid node */ 550 551 PetscCheck(j_c >= j_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 552 j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT, 553 j_start, j_c, j_start_ghost_c); 554 PetscCheck(i_c >= i_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 555 i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT, 556 i_start, i_c, i_start_ghost_c); 557 558 /* 559 Only include those interpolation points that are truly 560 nonzero. Note this is very important for final grid lines 561 in x and y directions; since they have no right/top neighbors 562 */ 563 nc = 0; 564 /* one left and below; or we are right on it */ 565 col = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 566 cols[nc++] = col_shift + idx_c[col]; 567 PetscCall(MatPreallocateSet(row, nc, cols, dnz, onz)); 568 } 569 } 570 PetscCall(MatCreate(PetscObjectComm((PetscObject)daf), &mat)); 571 #if defined(PETSC_HAVE_CUDA) 572 /* 573 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU 574 we don't want the original unconverted matrix copied to the GPU 575 */ 576 if (dof > 1) { PetscCall(MatBindToCPU(mat, PETSC_TRUE)); } 577 #endif 578 PetscCall(MatSetSizes(mat, m_f * n_f, col_scale * m_c * n_c, mx * my, col_scale * Mx * My)); 579 PetscCall(ConvertToAIJ(dac->mattype, &mattype)); 580 PetscCall(MatSetType(mat, mattype)); 581 PetscCall(MatSeqAIJSetPreallocation(mat, 0, dnz)); 582 PetscCall(MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz)); 583 MatPreallocateEnd(dnz, onz); 584 585 /* loop over local fine grid nodes setting interpolation for those*/ 586 for (j = j_start; j < j_start + n_f; j++) { 587 for (i = i_start; i < i_start + m_f; i++) { 588 /* convert to local "natural" numbering and then to PETSc global numbering */ 589 row = idx_f[(m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 590 591 i_c = (i / ratioi); /* coarse grid node to left of fine grid node */ 592 j_c = (j / ratioj); /* coarse grid node below fine grid node */ 593 nc = 0; 594 /* one left and below; or we are right on it */ 595 col = (m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 596 cols[nc] = col_shift + idx_c[col]; 597 v[nc++] = 1.0; 598 599 PetscCall(MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES)); 600 } 601 } 602 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 603 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 604 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 605 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 606 PetscCall(MatCreateMAIJ(mat, dof, A)); 607 PetscCall(MatDestroy(&mat)); 608 PetscCall(PetscLogFlops(13.0 * m_f * n_f)); 609 PetscFunctionReturn(0); 610 } 611 612 /* 613 Contributed by Jianming Yang <jianming-yang@uiowa.edu> 614 */ 615 PetscErrorCode DMCreateInterpolation_DA_3D_Q0(DM dac, DM daf, Mat *A) { 616 PetscInt i, j, l, i_start, j_start, l_start, m_f, n_f, p_f, Mx, My, Mz, dof; 617 const PetscInt *idx_c, *idx_f; 618 ISLocalToGlobalMapping ltog_f, ltog_c; 619 PetscInt m_ghost, n_ghost, p_ghost, m_ghost_c, n_ghost_c, p_ghost_c, nc, *dnz, *onz; 620 PetscInt row, col, i_start_ghost, j_start_ghost, l_start_ghost, cols[8], mx, m_c, my, n_c, mz, p_c, ratioi, ratioj, ratiol; 621 PetscInt i_c, j_c, l_c, i_start_c, j_start_c, l_start_c, i_start_ghost_c, j_start_ghost_c, l_start_ghost_c, col_shift, col_scale; 622 PetscMPIInt size_c, size_f, rank_f; 623 PetscScalar v[8]; 624 Mat mat; 625 DMBoundaryType bx, by, bz; 626 MatType mattype; 627 628 PetscFunctionBegin; 629 PetscCall(DMDAGetInfo(dac, NULL, &Mx, &My, &Mz, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL)); 630 PetscCheck(Mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be positive", Mx); 631 PetscCheck(My, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of y coarse grid points %" PetscInt_FMT " must be positive", My); 632 PetscCheck(Mz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of z coarse grid points %" PetscInt_FMT " must be positive", Mz); 633 PetscCall(DMDAGetInfo(daf, NULL, &mx, &my, &mz, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 634 ratioi = mx / Mx; 635 ratioj = my / My; 636 ratiol = mz / Mz; 637 PetscCheck(ratioi * Mx == mx, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Fine grid points must be multiple of coarse grid points in x"); 638 PetscCheck(ratioj * My == my, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Fine grid points must be multiple of coarse grid points in y"); 639 PetscCheck(ratiol * Mz == mz, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Fine grid points must be multiple of coarse grid points in z"); 640 PetscCheck(ratioi == 2 || ratioi == 1, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Coarsening factor in x must be 1 or 2"); 641 PetscCheck(ratioj == 2 || ratioj == 1, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Coarsening factor in y must be 1 or 2"); 642 PetscCheck(ratiol == 2 || ratiol == 1, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_WRONG, "Coarsening factor in z must be 1 or 2"); 643 644 PetscCall(DMDAGetCorners(daf, &i_start, &j_start, &l_start, &m_f, &n_f, &p_f)); 645 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, &l_start_ghost, &m_ghost, &n_ghost, &p_ghost)); 646 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 647 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 648 649 PetscCall(DMDAGetCorners(dac, &i_start_c, &j_start_c, &l_start_c, &m_c, &n_c, &p_c)); 650 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, &l_start_ghost_c, &m_ghost_c, &n_ghost_c, &p_ghost_c)); 651 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 652 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 653 654 /* 655 Used for handling a coarse DMDA that lives on 1/4 the processors of the fine DMDA. 656 The coarse vector is then duplicated 4 times (each time it lives on 1/4 of the 657 processors). It's effective length is hence 4 times its normal length, this is 658 why the col_scale is multiplied by the interpolation matrix column sizes. 659 sol_shift allows each set of 1/4 processors do its own interpolation using ITS 660 copy of the coarse vector. A bit of a hack but you do better. 661 662 In the standard case when size_f == size_c col_scale == 1 and col_shift == 0 663 */ 664 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dac), &size_c)); 665 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)daf), &size_f)); 666 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)daf), &rank_f)); 667 col_scale = size_f / size_c; 668 col_shift = Mx * My * Mz * (rank_f / size_c); 669 670 MatPreallocateBegin(PetscObjectComm((PetscObject)daf), m_f * n_f * p_f, col_scale * m_c * n_c * p_c, dnz, onz); 671 for (l = l_start; l < l_start + p_f; l++) { 672 for (j = j_start; j < j_start + n_f; j++) { 673 for (i = i_start; i < i_start + m_f; i++) { 674 /* convert to local "natural" numbering and then to PETSc global numbering */ 675 row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 676 677 i_c = (i / ratioi); /* coarse grid node to left of fine grid node */ 678 j_c = (j / ratioj); /* coarse grid node below fine grid node */ 679 l_c = (l / ratiol); 680 681 PetscCheck(l_c >= l_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 682 l_start %" PetscInt_FMT " l_c %" PetscInt_FMT " l_start_ghost_c %" PetscInt_FMT, 683 l_start, l_c, l_start_ghost_c); 684 PetscCheck(j_c >= j_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 685 j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT, 686 j_start, j_c, j_start_ghost_c); 687 PetscCheck(i_c >= i_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 688 i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT, 689 i_start, i_c, i_start_ghost_c); 690 691 /* 692 Only include those interpolation points that are truly 693 nonzero. Note this is very important for final grid lines 694 in x and y directions; since they have no right/top neighbors 695 */ 696 nc = 0; 697 /* one left and below; or we are right on it */ 698 col = (m_ghost_c * n_ghost_c * (l_c - l_start_ghost_c) + m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 699 cols[nc++] = col_shift + idx_c[col]; 700 PetscCall(MatPreallocateSet(row, nc, cols, dnz, onz)); 701 } 702 } 703 } 704 PetscCall(MatCreate(PetscObjectComm((PetscObject)daf), &mat)); 705 #if defined(PETSC_HAVE_CUDA) 706 /* 707 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU 708 we don't want the original unconverted matrix copied to the GPU 709 */ 710 if (dof > 1) { PetscCall(MatBindToCPU(mat, PETSC_TRUE)); } 711 #endif 712 PetscCall(MatSetSizes(mat, m_f * n_f * p_f, col_scale * m_c * n_c * p_c, mx * my * mz, col_scale * Mx * My * Mz)); 713 PetscCall(ConvertToAIJ(dac->mattype, &mattype)); 714 PetscCall(MatSetType(mat, mattype)); 715 PetscCall(MatSeqAIJSetPreallocation(mat, 0, dnz)); 716 PetscCall(MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz)); 717 MatPreallocateEnd(dnz, onz); 718 719 /* loop over local fine grid nodes setting interpolation for those*/ 720 for (l = l_start; l < l_start + p_f; l++) { 721 for (j = j_start; j < j_start + n_f; j++) { 722 for (i = i_start; i < i_start + m_f; i++) { 723 /* convert to local "natural" numbering and then to PETSc global numbering */ 724 row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 725 726 i_c = (i / ratioi); /* coarse grid node to left of fine grid node */ 727 j_c = (j / ratioj); /* coarse grid node below fine grid node */ 728 l_c = (l / ratiol); 729 nc = 0; 730 /* one left and below; or we are right on it */ 731 col = (m_ghost_c * n_ghost_c * (l_c - l_start_ghost_c) + m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 732 cols[nc] = col_shift + idx_c[col]; 733 v[nc++] = 1.0; 734 735 PetscCall(MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES)); 736 } 737 } 738 } 739 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 740 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 741 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 742 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 743 PetscCall(MatCreateMAIJ(mat, dof, A)); 744 PetscCall(MatDestroy(&mat)); 745 PetscCall(PetscLogFlops(13.0 * m_f * n_f * p_f)); 746 PetscFunctionReturn(0); 747 } 748 749 PetscErrorCode DMCreateInterpolation_DA_3D_Q1(DM dac, DM daf, Mat *A) { 750 PetscInt i, j, i_start, j_start, m_f, n_f, Mx, My, dof, l; 751 const PetscInt *idx_c, *idx_f; 752 ISLocalToGlobalMapping ltog_f, ltog_c; 753 PetscInt m_ghost, n_ghost, m_ghost_c, n_ghost_c, Mz, mz; 754 PetscInt row, col, i_start_ghost, j_start_ghost, cols[8], mx, m_c, my, nc, ratioi, ratioj, ratiok; 755 PetscInt i_c, j_c, i_start_c, j_start_c, n_c, i_start_ghost_c, j_start_ghost_c; 756 PetscInt l_start, p_f, l_start_ghost, p_ghost, l_start_c, p_c; 757 PetscInt l_start_ghost_c, p_ghost_c, l_c, *dnz, *onz; 758 PetscScalar v[8], x, y, z; 759 Mat mat; 760 DMBoundaryType bx, by, bz; 761 MatType mattype; 762 763 PetscFunctionBegin; 764 PetscCall(DMDAGetInfo(dac, NULL, &Mx, &My, &Mz, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL)); 765 PetscCall(DMDAGetInfo(daf, NULL, &mx, &my, &mz, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 766 if (mx == Mx) { 767 ratioi = 1; 768 } else if (bx == DM_BOUNDARY_PERIODIC) { 769 PetscCheck(Mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be positive", Mx); 770 ratioi = mx / Mx; 771 PetscCheck(ratioi * Mx == mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mx/Mx must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 772 } else { 773 PetscCheck(Mx >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of x coarse grid points %" PetscInt_FMT " must be greater than 1", Mx); 774 ratioi = (mx - 1) / (Mx - 1); 775 PetscCheck(ratioi * (Mx - 1) == mx - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 776 } 777 if (my == My) { 778 ratioj = 1; 779 } else if (by == DM_BOUNDARY_PERIODIC) { 780 PetscCheck(My, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of y coarse grid points %" PetscInt_FMT " must be positive", My); 781 ratioj = my / My; 782 PetscCheck(ratioj * My == my, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: my/My must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 783 } else { 784 PetscCheck(My >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of y coarse grid points %" PetscInt_FMT " must be greater than 1", My); 785 ratioj = (my - 1) / (My - 1); 786 PetscCheck(ratioj * (My - 1) == my - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (my - 1)/(My - 1) must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 787 } 788 if (mz == Mz) { 789 ratiok = 1; 790 } else if (bz == DM_BOUNDARY_PERIODIC) { 791 PetscCheck(Mz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of z coarse grid points %" PetscInt_FMT " must be positive", Mz); 792 ratiok = mz / Mz; 793 PetscCheck(ratiok * Mz == mz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mz/Mz must be integer: mz %" PetscInt_FMT " Mz %" PetscInt_FMT, mz, Mz); 794 } else { 795 PetscCheck(Mz >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of z coarse grid points %" PetscInt_FMT " must be greater than 1", Mz); 796 ratiok = (mz - 1) / (Mz - 1); 797 PetscCheck(ratiok * (Mz - 1) == mz - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %" PetscInt_FMT " Mz %" PetscInt_FMT, mz, Mz); 798 } 799 800 PetscCall(DMDAGetCorners(daf, &i_start, &j_start, &l_start, &m_f, &n_f, &p_f)); 801 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, &l_start_ghost, &m_ghost, &n_ghost, &p_ghost)); 802 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 803 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 804 805 PetscCall(DMDAGetCorners(dac, &i_start_c, &j_start_c, &l_start_c, &m_c, &n_c, &p_c)); 806 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, &l_start_ghost_c, &m_ghost_c, &n_ghost_c, &p_ghost_c)); 807 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 808 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 809 810 /* create interpolation matrix, determining exact preallocation */ 811 MatPreallocateBegin(PetscObjectComm((PetscObject)dac), m_f * n_f * p_f, m_c * n_c * p_c, dnz, onz); 812 /* loop over local fine grid nodes counting interpolating points */ 813 for (l = l_start; l < l_start + p_f; l++) { 814 for (j = j_start; j < j_start + n_f; j++) { 815 for (i = i_start; i < i_start + m_f; i++) { 816 /* convert to local "natural" numbering and then to PETSc global numbering */ 817 row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 818 i_c = (i / ratioi); 819 j_c = (j / ratioj); 820 l_c = (l / ratiok); 821 PetscCheck(l_c >= l_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 822 l_start %" PetscInt_FMT " l_c %" PetscInt_FMT " l_start_ghost_c %" PetscInt_FMT, 823 l_start, l_c, l_start_ghost_c); 824 PetscCheck(j_c >= j_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 825 j_start %" PetscInt_FMT " j_c %" PetscInt_FMT " j_start_ghost_c %" PetscInt_FMT, 826 j_start, j_c, j_start_ghost_c); 827 PetscCheck(i_c >= i_start_ghost_c, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 828 i_start %" PetscInt_FMT " i_c %" PetscInt_FMT " i_start_ghost_c %" PetscInt_FMT, 829 i_start, i_c, i_start_ghost_c); 830 831 /* 832 Only include those interpolation points that are truly 833 nonzero. Note this is very important for final grid lines 834 in x and y directions; since they have no right/top neighbors 835 */ 836 nc = 0; 837 col = (m_ghost_c * n_ghost_c * (l_c - l_start_ghost_c) + m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 838 cols[nc++] = idx_c[col]; 839 if (i_c * ratioi != i) { cols[nc++] = idx_c[col + 1]; } 840 if (j_c * ratioj != j) { cols[nc++] = idx_c[col + m_ghost_c]; } 841 if (l_c * ratiok != l) { cols[nc++] = idx_c[col + m_ghost_c * n_ghost_c]; } 842 if (j_c * ratioj != j && i_c * ratioi != i) { cols[nc++] = idx_c[col + (m_ghost_c + 1)]; } 843 if (j_c * ratioj != j && l_c * ratiok != l) { cols[nc++] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c)]; } 844 if (i_c * ratioi != i && l_c * ratiok != l) { cols[nc++] = idx_c[col + (m_ghost_c * n_ghost_c + 1)]; } 845 if (i_c * ratioi != i && l_c * ratiok != l && j_c * ratioj != j) { cols[nc++] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c + 1)]; } 846 PetscCall(MatPreallocateSet(row, nc, cols, dnz, onz)); 847 } 848 } 849 } 850 PetscCall(MatCreate(PetscObjectComm((PetscObject)dac), &mat)); 851 #if defined(PETSC_HAVE_CUDA) 852 /* 853 Temporary hack: Since the MAIJ matrix must be converted to AIJ before being used by the GPU 854 we don't want the original unconverted matrix copied to the GPU 855 */ 856 if (dof > 1) { PetscCall(MatBindToCPU(mat, PETSC_TRUE)); } 857 #endif 858 PetscCall(MatSetSizes(mat, m_f * n_f * p_f, m_c * n_c * p_c, mx * my * mz, Mx * My * Mz)); 859 PetscCall(ConvertToAIJ(dac->mattype, &mattype)); 860 PetscCall(MatSetType(mat, mattype)); 861 PetscCall(MatSeqAIJSetPreallocation(mat, 0, dnz)); 862 PetscCall(MatMPIAIJSetPreallocation(mat, 0, dnz, 0, onz)); 863 MatPreallocateEnd(dnz, onz); 864 865 /* loop over local fine grid nodes setting interpolation for those*/ 866 if (!NEWVERSION) { 867 for (l = l_start; l < l_start + p_f; l++) { 868 for (j = j_start; j < j_start + n_f; j++) { 869 for (i = i_start; i < i_start + m_f; i++) { 870 /* convert to local "natural" numbering and then to PETSc global numbering */ 871 row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 872 873 i_c = (i / ratioi); 874 j_c = (j / ratioj); 875 l_c = (l / ratiok); 876 877 /* 878 Only include those interpolation points that are truly 879 nonzero. Note this is very important for final grid lines 880 in x and y directions; since they have no right/top neighbors 881 */ 882 x = ((PetscReal)(i - i_c * ratioi)) / ((PetscReal)ratioi); 883 y = ((PetscReal)(j - j_c * ratioj)) / ((PetscReal)ratioj); 884 z = ((PetscReal)(l - l_c * ratiok)) / ((PetscReal)ratiok); 885 886 nc = 0; 887 /* one left and below; or we are right on it */ 888 col = (m_ghost_c * n_ghost_c * (l_c - l_start_ghost_c) + m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 889 890 cols[nc] = idx_c[col]; 891 v[nc++] = .125 * (1. - (2.0 * x - 1.)) * (1. - (2.0 * y - 1.)) * (1. - (2.0 * z - 1.)); 892 893 if (i_c * ratioi != i) { 894 cols[nc] = idx_c[col + 1]; 895 v[nc++] = .125 * (1. + (2.0 * x - 1.)) * (1. - (2.0 * y - 1.)) * (1. - (2.0 * z - 1.)); 896 } 897 898 if (j_c * ratioj != j) { 899 cols[nc] = idx_c[col + m_ghost_c]; 900 v[nc++] = .125 * (1. - (2.0 * x - 1.)) * (1. + (2.0 * y - 1.)) * (1. - (2.0 * z - 1.)); 901 } 902 903 if (l_c * ratiok != l) { 904 cols[nc] = idx_c[col + m_ghost_c * n_ghost_c]; 905 v[nc++] = .125 * (1. - (2.0 * x - 1.)) * (1. - (2.0 * y - 1.)) * (1. + (2.0 * z - 1.)); 906 } 907 908 if (j_c * ratioj != j && i_c * ratioi != i) { 909 cols[nc] = idx_c[col + (m_ghost_c + 1)]; 910 v[nc++] = .125 * (1. + (2.0 * x - 1.)) * (1. + (2.0 * y - 1.)) * (1. - (2.0 * z - 1.)); 911 } 912 913 if (j_c * ratioj != j && l_c * ratiok != l) { 914 cols[nc] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c)]; 915 v[nc++] = .125 * (1. - (2.0 * x - 1.)) * (1. + (2.0 * y - 1.)) * (1. + (2.0 * z - 1.)); 916 } 917 918 if (i_c * ratioi != i && l_c * ratiok != l) { 919 cols[nc] = idx_c[col + (m_ghost_c * n_ghost_c + 1)]; 920 v[nc++] = .125 * (1. + (2.0 * x - 1.)) * (1. - (2.0 * y - 1.)) * (1. + (2.0 * z - 1.)); 921 } 922 923 if (i_c * ratioi != i && l_c * ratiok != l && j_c * ratioj != j) { 924 cols[nc] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c + 1)]; 925 v[nc++] = .125 * (1. + (2.0 * x - 1.)) * (1. + (2.0 * y - 1.)) * (1. + (2.0 * z - 1.)); 926 } 927 PetscCall(MatSetValues(mat, 1, &row, nc, cols, v, INSERT_VALUES)); 928 } 929 } 930 } 931 932 } else { 933 PetscScalar *xi, *eta, *zeta; 934 PetscInt li, nxi, lj, neta, lk, nzeta, n; 935 PetscScalar Ni[8]; 936 937 /* compute local coordinate arrays */ 938 nxi = ratioi + 1; 939 neta = ratioj + 1; 940 nzeta = ratiok + 1; 941 PetscCall(PetscMalloc1(nxi, &xi)); 942 PetscCall(PetscMalloc1(neta, &eta)); 943 PetscCall(PetscMalloc1(nzeta, &zeta)); 944 for (li = 0; li < nxi; li++) xi[li] = -1.0 + (PetscScalar)li * (2.0 / (PetscScalar)(nxi - 1)); 945 for (lj = 0; lj < neta; lj++) eta[lj] = -1.0 + (PetscScalar)lj * (2.0 / (PetscScalar)(neta - 1)); 946 for (lk = 0; lk < nzeta; lk++) zeta[lk] = -1.0 + (PetscScalar)lk * (2.0 / (PetscScalar)(nzeta - 1)); 947 948 for (l = l_start; l < l_start + p_f; l++) { 949 for (j = j_start; j < j_start + n_f; j++) { 950 for (i = i_start; i < i_start + m_f; i++) { 951 /* convert to local "natural" numbering and then to PETSc global numbering */ 952 row = idx_f[(m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))]; 953 954 i_c = (i / ratioi); 955 j_c = (j / ratioj); 956 l_c = (l / ratiok); 957 958 /* remainders */ 959 li = i - ratioi * (i / ratioi); 960 if (i == mx - 1) li = nxi - 1; 961 lj = j - ratioj * (j / ratioj); 962 if (j == my - 1) lj = neta - 1; 963 lk = l - ratiok * (l / ratiok); 964 if (l == mz - 1) lk = nzeta - 1; 965 966 /* corners */ 967 col = (m_ghost_c * n_ghost_c * (l_c - l_start_ghost_c) + m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 968 cols[0] = idx_c[col]; 969 Ni[0] = 1.0; 970 if ((li == 0) || (li == nxi - 1)) { 971 if ((lj == 0) || (lj == neta - 1)) { 972 if ((lk == 0) || (lk == nzeta - 1)) { 973 PetscCall(MatSetValue(mat, row, cols[0], Ni[0], INSERT_VALUES)); 974 continue; 975 } 976 } 977 } 978 979 /* edges + interior */ 980 /* remainders */ 981 if (i == mx - 1) i_c--; 982 if (j == my - 1) j_c--; 983 if (l == mz - 1) l_c--; 984 985 col = (m_ghost_c * n_ghost_c * (l_c - l_start_ghost_c) + m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c)); 986 cols[0] = idx_c[col]; /* one left and below; or we are right on it */ 987 cols[1] = idx_c[col + 1]; /* one right and below */ 988 cols[2] = idx_c[col + m_ghost_c]; /* one left and above */ 989 cols[3] = idx_c[col + (m_ghost_c + 1)]; /* one right and above */ 990 991 cols[4] = idx_c[col + m_ghost_c * n_ghost_c]; /* one left and below and front; or we are right on it */ 992 cols[5] = idx_c[col + (m_ghost_c * n_ghost_c + 1)]; /* one right and below, and front */ 993 cols[6] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c)]; /* one left and above and front*/ 994 cols[7] = idx_c[col + (m_ghost_c * n_ghost_c + m_ghost_c + 1)]; /* one right and above and front */ 995 996 Ni[0] = 0.125 * (1.0 - xi[li]) * (1.0 - eta[lj]) * (1.0 - zeta[lk]); 997 Ni[1] = 0.125 * (1.0 + xi[li]) * (1.0 - eta[lj]) * (1.0 - zeta[lk]); 998 Ni[2] = 0.125 * (1.0 - xi[li]) * (1.0 + eta[lj]) * (1.0 - zeta[lk]); 999 Ni[3] = 0.125 * (1.0 + xi[li]) * (1.0 + eta[lj]) * (1.0 - zeta[lk]); 1000 1001 Ni[4] = 0.125 * (1.0 - xi[li]) * (1.0 - eta[lj]) * (1.0 + zeta[lk]); 1002 Ni[5] = 0.125 * (1.0 + xi[li]) * (1.0 - eta[lj]) * (1.0 + zeta[lk]); 1003 Ni[6] = 0.125 * (1.0 - xi[li]) * (1.0 + eta[lj]) * (1.0 + zeta[lk]); 1004 Ni[7] = 0.125 * (1.0 + xi[li]) * (1.0 + eta[lj]) * (1.0 + zeta[lk]); 1005 1006 for (n = 0; n < 8; n++) { 1007 if (PetscAbsScalar(Ni[n]) < 1.0e-32) cols[n] = -1; 1008 } 1009 PetscCall(MatSetValues(mat, 1, &row, 8, cols, Ni, INSERT_VALUES)); 1010 } 1011 } 1012 } 1013 PetscCall(PetscFree(xi)); 1014 PetscCall(PetscFree(eta)); 1015 PetscCall(PetscFree(zeta)); 1016 } 1017 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 1018 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 1019 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY)); 1020 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); 1021 1022 PetscCall(MatCreateMAIJ(mat, dof, A)); 1023 PetscCall(MatDestroy(&mat)); 1024 PetscFunctionReturn(0); 1025 } 1026 1027 PetscErrorCode DMCreateInterpolation_DA(DM dac, DM daf, Mat *A, Vec *scale) { 1028 PetscInt dimc, Mc, Nc, Pc, mc, nc, pc, dofc, sc, dimf, Mf, Nf, Pf, mf, nf, pf, doff, sf; 1029 DMBoundaryType bxc, byc, bzc, bxf, byf, bzf; 1030 DMDAStencilType stc, stf; 1031 DM_DA *ddc = (DM_DA *)dac->data; 1032 1033 PetscFunctionBegin; 1034 PetscValidHeaderSpecific(dac, DM_CLASSID, 1); 1035 PetscValidHeaderSpecific(daf, DM_CLASSID, 2); 1036 PetscValidPointer(A, 3); 1037 if (scale) PetscValidPointer(scale, 4); 1038 1039 PetscCall(DMDAGetInfo(dac, &dimc, &Mc, &Nc, &Pc, &mc, &nc, &pc, &dofc, &sc, &bxc, &byc, &bzc, &stc)); 1040 PetscCall(DMDAGetInfo(daf, &dimf, &Mf, &Nf, &Pf, &mf, &nf, &pf, &doff, &sf, &bxf, &byf, &bzf, &stf)); 1041 PetscCheck(dimc == dimf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Dimensions of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, dimc, dimf); 1042 PetscCheck(dofc == doff, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "DOF of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, dofc, doff); 1043 PetscCheck(sc == sf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Stencil width of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, sc, sf); 1044 PetscCheck(bxc == bxf && byc == byf && bzc == bzf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Boundary type different in two DMDAs"); 1045 PetscCheck(stc == stf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Stencil type different in two DMDAs"); 1046 PetscCheck(Mc >= 2 || Mf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coarse grid requires at least 2 points in x direction"); 1047 PetscCheck(dimc <= 1 || Nc >= 2 || Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coarse grid requires at least 2 points in y direction"); 1048 PetscCheck(dimc <= 2 || Pc >= 2 || Pf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coarse grid requires at least 2 points in z direction"); 1049 1050 if (ddc->interptype == DMDA_Q1) { 1051 if (dimc == 1) { 1052 PetscCall(DMCreateInterpolation_DA_1D_Q1(dac, daf, A)); 1053 } else if (dimc == 2) { 1054 PetscCall(DMCreateInterpolation_DA_2D_Q1(dac, daf, A)); 1055 } else if (dimc == 3) { 1056 PetscCall(DMCreateInterpolation_DA_3D_Q1(dac, daf, A)); 1057 } else SETERRQ(PetscObjectComm((PetscObject)daf), PETSC_ERR_SUP, "No support for this DMDA dimension %" PetscInt_FMT " for interpolation type %d", dimc, (int)ddc->interptype); 1058 } else if (ddc->interptype == DMDA_Q0) { 1059 if (dimc == 1) { 1060 PetscCall(DMCreateInterpolation_DA_1D_Q0(dac, daf, A)); 1061 } else if (dimc == 2) { 1062 PetscCall(DMCreateInterpolation_DA_2D_Q0(dac, daf, A)); 1063 } else if (dimc == 3) { 1064 PetscCall(DMCreateInterpolation_DA_3D_Q0(dac, daf, A)); 1065 } else SETERRQ(PetscObjectComm((PetscObject)daf), PETSC_ERR_SUP, "No support for this DMDA dimension %" PetscInt_FMT " for interpolation type %d", dimc, (int)ddc->interptype); 1066 } 1067 if (scale) PetscCall(DMCreateInterpolationScale((DM)dac, (DM)daf, *A, scale)); 1068 PetscFunctionReturn(0); 1069 } 1070 1071 PetscErrorCode DMCreateInjection_DA_1D(DM dac, DM daf, VecScatter *inject) { 1072 PetscInt i, i_start, m_f, Mx, dof; 1073 const PetscInt *idx_f; 1074 ISLocalToGlobalMapping ltog_f; 1075 PetscInt m_ghost, m_ghost_c; 1076 PetscInt row, i_start_ghost, mx, m_c, nc, ratioi; 1077 PetscInt i_start_c, i_start_ghost_c; 1078 PetscInt *cols; 1079 DMBoundaryType bx; 1080 Vec vecf, vecc; 1081 IS isf; 1082 1083 PetscFunctionBegin; 1084 PetscCall(DMDAGetInfo(dac, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bx, NULL, NULL, NULL)); 1085 PetscCall(DMDAGetInfo(daf, NULL, &mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 1086 if (bx == DM_BOUNDARY_PERIODIC) { 1087 ratioi = mx / Mx; 1088 PetscCheck(ratioi * Mx == mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mx/Mx must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 1089 } else { 1090 ratioi = (mx - 1) / (Mx - 1); 1091 PetscCheck(ratioi * (Mx - 1) == mx - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 1092 } 1093 1094 PetscCall(DMDAGetCorners(daf, &i_start, NULL, NULL, &m_f, NULL, NULL)); 1095 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, NULL, NULL, &m_ghost, NULL, NULL)); 1096 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 1097 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 1098 1099 PetscCall(DMDAGetCorners(dac, &i_start_c, NULL, NULL, &m_c, NULL, NULL)); 1100 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, NULL, NULL, &m_ghost_c, NULL, NULL)); 1101 1102 /* loop over local fine grid nodes setting interpolation for those*/ 1103 nc = 0; 1104 PetscCall(PetscMalloc1(m_f, &cols)); 1105 1106 for (i = i_start_c; i < i_start_c + m_c; i++) { 1107 PetscInt i_f = i * ratioi; 1108 1109 PetscCheck(i_f >= i_start_ghost && i_f < i_start_ghost + m_ghost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\ni_c %" PetscInt_FMT " i_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]", i, i_f, i_start_ghost, i_start_ghost + m_ghost); 1110 1111 row = idx_f[(i_f - i_start_ghost)]; 1112 cols[nc++] = row; 1113 } 1114 1115 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 1116 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)daf), dof, nc, cols, PETSC_OWN_POINTER, &isf)); 1117 PetscCall(DMGetGlobalVector(dac, &vecc)); 1118 PetscCall(DMGetGlobalVector(daf, &vecf)); 1119 PetscCall(VecScatterCreate(vecf, isf, vecc, NULL, inject)); 1120 PetscCall(DMRestoreGlobalVector(dac, &vecc)); 1121 PetscCall(DMRestoreGlobalVector(daf, &vecf)); 1122 PetscCall(ISDestroy(&isf)); 1123 PetscFunctionReturn(0); 1124 } 1125 1126 PetscErrorCode DMCreateInjection_DA_2D(DM dac, DM daf, VecScatter *inject) { 1127 PetscInt i, j, i_start, j_start, m_f, n_f, Mx, My, dof; 1128 const PetscInt *idx_c, *idx_f; 1129 ISLocalToGlobalMapping ltog_f, ltog_c; 1130 PetscInt m_ghost, n_ghost, m_ghost_c, n_ghost_c; 1131 PetscInt row, i_start_ghost, j_start_ghost, mx, m_c, my, nc, ratioi, ratioj; 1132 PetscInt i_start_c, j_start_c, n_c, i_start_ghost_c, j_start_ghost_c; 1133 PetscInt *cols; 1134 DMBoundaryType bx, by; 1135 Vec vecf, vecc; 1136 IS isf; 1137 1138 PetscFunctionBegin; 1139 PetscCall(DMDAGetInfo(dac, NULL, &Mx, &My, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL)); 1140 PetscCall(DMDAGetInfo(daf, NULL, &mx, &my, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 1141 if (bx == DM_BOUNDARY_PERIODIC) { 1142 ratioi = mx / Mx; 1143 PetscCheck(ratioi * Mx == mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mx/Mx must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 1144 } else { 1145 ratioi = (mx - 1) / (Mx - 1); 1146 PetscCheck(ratioi * (Mx - 1) == mx - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 1147 } 1148 if (by == DM_BOUNDARY_PERIODIC) { 1149 ratioj = my / My; 1150 PetscCheck(ratioj * My == my, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: my/My must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 1151 } else { 1152 ratioj = (my - 1) / (My - 1); 1153 PetscCheck(ratioj * (My - 1) == my - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (my - 1)/(My - 1) must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 1154 } 1155 1156 PetscCall(DMDAGetCorners(daf, &i_start, &j_start, NULL, &m_f, &n_f, NULL)); 1157 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, NULL, &m_ghost, &n_ghost, NULL)); 1158 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 1159 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 1160 1161 PetscCall(DMDAGetCorners(dac, &i_start_c, &j_start_c, NULL, &m_c, &n_c, NULL)); 1162 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, NULL, &m_ghost_c, &n_ghost_c, NULL)); 1163 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 1164 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 1165 1166 /* loop over local fine grid nodes setting interpolation for those*/ 1167 nc = 0; 1168 PetscCall(PetscMalloc1(n_f * m_f, &cols)); 1169 for (j = j_start_c; j < j_start_c + n_c; j++) { 1170 for (i = i_start_c; i < i_start_c + m_c; i++) { 1171 PetscInt i_f = i * ratioi, j_f = j * ratioj; 1172 PetscCheck(j_f >= j_start_ghost && j_f < j_start_ghost + n_ghost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 1173 j_c %" PetscInt_FMT " j_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]", 1174 j, j_f, j_start_ghost, j_start_ghost + n_ghost); 1175 PetscCheck(i_f >= i_start_ghost && i_f < i_start_ghost + m_ghost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Processor's coarse DMDA must lie over fine DMDA\n\ 1176 i_c %" PetscInt_FMT " i_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]", 1177 i, i_f, i_start_ghost, i_start_ghost + m_ghost); 1178 row = idx_f[(m_ghost * (j_f - j_start_ghost) + (i_f - i_start_ghost))]; 1179 cols[nc++] = row; 1180 } 1181 } 1182 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 1183 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 1184 1185 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)daf), dof, nc, cols, PETSC_OWN_POINTER, &isf)); 1186 PetscCall(DMGetGlobalVector(dac, &vecc)); 1187 PetscCall(DMGetGlobalVector(daf, &vecf)); 1188 PetscCall(VecScatterCreate(vecf, isf, vecc, NULL, inject)); 1189 PetscCall(DMRestoreGlobalVector(dac, &vecc)); 1190 PetscCall(DMRestoreGlobalVector(daf, &vecf)); 1191 PetscCall(ISDestroy(&isf)); 1192 PetscFunctionReturn(0); 1193 } 1194 1195 PetscErrorCode DMCreateInjection_DA_3D(DM dac, DM daf, VecScatter *inject) { 1196 PetscInt i, j, k, i_start, j_start, k_start, m_f, n_f, p_f, Mx, My, Mz; 1197 PetscInt m_ghost, n_ghost, p_ghost, m_ghost_c, n_ghost_c, p_ghost_c; 1198 PetscInt i_start_ghost, j_start_ghost, k_start_ghost; 1199 PetscInt mx, my, mz, ratioi, ratioj, ratiok; 1200 PetscInt i_start_c, j_start_c, k_start_c; 1201 PetscInt m_c, n_c, p_c; 1202 PetscInt i_start_ghost_c, j_start_ghost_c, k_start_ghost_c; 1203 PetscInt row, nc, dof; 1204 const PetscInt *idx_c, *idx_f; 1205 ISLocalToGlobalMapping ltog_f, ltog_c; 1206 PetscInt *cols; 1207 DMBoundaryType bx, by, bz; 1208 Vec vecf, vecc; 1209 IS isf; 1210 1211 PetscFunctionBegin; 1212 PetscCall(DMDAGetInfo(dac, NULL, &Mx, &My, &Mz, NULL, NULL, NULL, NULL, NULL, &bx, &by, &bz, NULL)); 1213 PetscCall(DMDAGetInfo(daf, NULL, &mx, &my, &mz, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); 1214 1215 if (bx == DM_BOUNDARY_PERIODIC) { 1216 ratioi = mx / Mx; 1217 PetscCheck(ratioi * Mx == mx, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mx/Mx must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 1218 } else { 1219 ratioi = (mx - 1) / (Mx - 1); 1220 PetscCheck(ratioi * (Mx - 1) == mx - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mx - 1)/(Mx - 1) must be integer: mx %" PetscInt_FMT " Mx %" PetscInt_FMT, mx, Mx); 1221 } 1222 if (by == DM_BOUNDARY_PERIODIC) { 1223 ratioj = my / My; 1224 PetscCheck(ratioj * My == my, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: my/My must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 1225 } else { 1226 ratioj = (my - 1) / (My - 1); 1227 PetscCheck(ratioj * (My - 1) == my - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (my - 1)/(My - 1) must be integer: my %" PetscInt_FMT " My %" PetscInt_FMT, my, My); 1228 } 1229 if (bz == DM_BOUNDARY_PERIODIC) { 1230 ratiok = mz / Mz; 1231 PetscCheck(ratiok * Mz == mz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: mz/Mz must be integer: mz %" PetscInt_FMT " My %" PetscInt_FMT, mz, Mz); 1232 } else { 1233 ratiok = (mz - 1) / (Mz - 1); 1234 PetscCheck(ratiok * (Mz - 1) == mz - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Ratio between levels: (mz - 1)/(Mz - 1) must be integer: mz %" PetscInt_FMT " Mz %" PetscInt_FMT, mz, Mz); 1235 } 1236 1237 PetscCall(DMDAGetCorners(daf, &i_start, &j_start, &k_start, &m_f, &n_f, &p_f)); 1238 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, &k_start_ghost, &m_ghost, &n_ghost, &p_ghost)); 1239 PetscCall(DMGetLocalToGlobalMapping(daf, <og_f)); 1240 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_f, &idx_f)); 1241 1242 PetscCall(DMDAGetCorners(dac, &i_start_c, &j_start_c, &k_start_c, &m_c, &n_c, &p_c)); 1243 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, &k_start_ghost_c, &m_ghost_c, &n_ghost_c, &p_ghost_c)); 1244 PetscCall(DMGetLocalToGlobalMapping(dac, <og_c)); 1245 PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog_c, &idx_c)); 1246 1247 /* loop over local fine grid nodes setting interpolation for those*/ 1248 nc = 0; 1249 PetscCall(PetscMalloc1(n_f * m_f * p_f, &cols)); 1250 for (k = k_start_c; k < k_start_c + p_c; k++) { 1251 for (j = j_start_c; j < j_start_c + n_c; j++) { 1252 for (i = i_start_c; i < i_start_c + m_c; i++) { 1253 PetscInt i_f = i * ratioi, j_f = j * ratioj, k_f = k * ratiok; 1254 PetscCheck(k_f >= k_start_ghost && k_f < k_start_ghost + p_ghost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 1255 "Processor's coarse DMDA must lie over fine DMDA " 1256 "k_c %" PetscInt_FMT " k_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]", 1257 k, k_f, k_start_ghost, k_start_ghost + p_ghost); 1258 PetscCheck(j_f >= j_start_ghost && j_f < j_start_ghost + n_ghost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 1259 "Processor's coarse DMDA must lie over fine DMDA " 1260 "j_c %" PetscInt_FMT " j_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]", 1261 j, j_f, j_start_ghost, j_start_ghost + n_ghost); 1262 PetscCheck(i_f >= i_start_ghost && i_f < i_start_ghost + m_ghost, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, 1263 "Processor's coarse DMDA must lie over fine DMDA " 1264 "i_c %" PetscInt_FMT " i_f %" PetscInt_FMT " fine ghost range [%" PetscInt_FMT ",%" PetscInt_FMT "]", 1265 i, i_f, i_start_ghost, i_start_ghost + m_ghost); 1266 row = idx_f[(m_ghost * n_ghost * (k_f - k_start_ghost) + m_ghost * (j_f - j_start_ghost) + (i_f - i_start_ghost))]; 1267 cols[nc++] = row; 1268 } 1269 } 1270 } 1271 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_f, &idx_f)); 1272 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog_c, &idx_c)); 1273 1274 PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)daf), dof, nc, cols, PETSC_OWN_POINTER, &isf)); 1275 PetscCall(DMGetGlobalVector(dac, &vecc)); 1276 PetscCall(DMGetGlobalVector(daf, &vecf)); 1277 PetscCall(VecScatterCreate(vecf, isf, vecc, NULL, inject)); 1278 PetscCall(DMRestoreGlobalVector(dac, &vecc)); 1279 PetscCall(DMRestoreGlobalVector(daf, &vecf)); 1280 PetscCall(ISDestroy(&isf)); 1281 PetscFunctionReturn(0); 1282 } 1283 1284 PetscErrorCode DMCreateInjection_DA(DM dac, DM daf, Mat *mat) { 1285 PetscInt dimc, Mc, Nc, Pc, mc, nc, pc, dofc, sc, dimf, Mf, Nf, Pf, mf, nf, pf, doff, sf; 1286 DMBoundaryType bxc, byc, bzc, bxf, byf, bzf; 1287 DMDAStencilType stc, stf; 1288 VecScatter inject = NULL; 1289 1290 PetscFunctionBegin; 1291 PetscValidHeaderSpecific(dac, DM_CLASSID, 1); 1292 PetscValidHeaderSpecific(daf, DM_CLASSID, 2); 1293 PetscValidPointer(mat, 3); 1294 1295 PetscCall(DMDAGetInfo(dac, &dimc, &Mc, &Nc, &Pc, &mc, &nc, &pc, &dofc, &sc, &bxc, &byc, &bzc, &stc)); 1296 PetscCall(DMDAGetInfo(daf, &dimf, &Mf, &Nf, &Pf, &mf, &nf, &pf, &doff, &sf, &bxf, &byf, &bzf, &stf)); 1297 PetscCheck(dimc == dimf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Dimensions of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, dimc, dimf); 1298 PetscCheck(dofc == doff, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "DOF of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, dofc, doff); 1299 PetscCheck(sc == sf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Stencil width of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, sc, sf); 1300 PetscCheck(bxc == bxf && byc == byf && bzc == bzf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Boundary type different in two DMDAs"); 1301 PetscCheck(stc == stf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Stencil type different in two DMDAs"); 1302 PetscCheck(Mc >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coarse grid requires at least 2 points in x direction"); 1303 PetscCheck(dimc <= 1 || Nc >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coarse grid requires at least 2 points in y direction"); 1304 PetscCheck(dimc <= 2 || Pc >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coarse grid requires at least 2 points in z direction"); 1305 1306 if (dimc == 1) { 1307 PetscCall(DMCreateInjection_DA_1D(dac, daf, &inject)); 1308 } else if (dimc == 2) { 1309 PetscCall(DMCreateInjection_DA_2D(dac, daf, &inject)); 1310 } else if (dimc == 3) { 1311 PetscCall(DMCreateInjection_DA_3D(dac, daf, &inject)); 1312 } 1313 PetscCall(MatCreateScatter(PetscObjectComm((PetscObject)inject), inject, mat)); 1314 PetscCall(VecScatterDestroy(&inject)); 1315 PetscFunctionReturn(0); 1316 } 1317 1318 /*@ 1319 DMCreateAggregates - Deprecated, see DMDACreateAggregates. 1320 1321 Level: intermediate 1322 @*/ 1323 PetscErrorCode DMCreateAggregates(DM dac, DM daf, Mat *mat) { 1324 return DMDACreateAggregates(dac, daf, mat); 1325 } 1326 1327 /*@ 1328 DMDACreateAggregates - Gets the aggregates that map between 1329 grids associated with two DMDAs. 1330 1331 Collective on dmc 1332 1333 Input Parameters: 1334 + dmc - the coarse grid DMDA 1335 - dmf - the fine grid DMDA 1336 1337 Output Parameters: 1338 . rest - the restriction matrix (transpose of the projection matrix) 1339 1340 Level: intermediate 1341 1342 Note: This routine is not used by PETSc. 1343 It is not clear what its use case is and it may be removed in a future release. 1344 Users should contact petsc-maint@mcs.anl.gov if they plan to use it. 1345 1346 .seealso: `DMRefine()`, `DMCreateInjection()`, `DMCreateInterpolation()` 1347 @*/ 1348 PetscErrorCode DMDACreateAggregates(DM dac, DM daf, Mat *rest) { 1349 PetscInt dimc, Mc, Nc, Pc, mc, nc, pc, dofc, sc; 1350 PetscInt dimf, Mf, Nf, Pf, mf, nf, pf, doff, sf; 1351 DMBoundaryType bxc, byc, bzc, bxf, byf, bzf; 1352 DMDAStencilType stc, stf; 1353 PetscInt i, j, l; 1354 PetscInt i_start, j_start, l_start, m_f, n_f, p_f; 1355 PetscInt i_start_ghost, j_start_ghost, l_start_ghost, m_ghost, n_ghost, p_ghost; 1356 const PetscInt *idx_f; 1357 PetscInt i_c, j_c, l_c; 1358 PetscInt i_start_c, j_start_c, l_start_c, m_c, n_c, p_c; 1359 PetscInt i_start_ghost_c, j_start_ghost_c, l_start_ghost_c, m_ghost_c, n_ghost_c, p_ghost_c; 1360 const PetscInt *idx_c; 1361 PetscInt d; 1362 PetscInt a; 1363 PetscInt max_agg_size; 1364 PetscInt *fine_nodes; 1365 PetscScalar *one_vec; 1366 PetscInt fn_idx; 1367 ISLocalToGlobalMapping ltogmf, ltogmc; 1368 1369 PetscFunctionBegin; 1370 PetscValidHeaderSpecificType(dac, DM_CLASSID, 1, DMDA); 1371 PetscValidHeaderSpecificType(daf, DM_CLASSID, 2, DMDA); 1372 PetscValidPointer(rest, 3); 1373 1374 PetscCall(DMDAGetInfo(dac, &dimc, &Mc, &Nc, &Pc, &mc, &nc, &pc, &dofc, &sc, &bxc, &byc, &bzc, &stc)); 1375 PetscCall(DMDAGetInfo(daf, &dimf, &Mf, &Nf, &Pf, &mf, &nf, &pf, &doff, &sf, &bxf, &byf, &bzf, &stf)); 1376 PetscCheck(dimc == dimf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Dimensions of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, dimc, dimf); 1377 PetscCheck(dofc == doff, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "DOF of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, dofc, doff); 1378 PetscCheck(sc == sf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Stencil width of DMDA do not match %" PetscInt_FMT " %" PetscInt_FMT, sc, sf); 1379 PetscCheck(bxc == bxf && byc == byf && bzc == bzf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Boundary type different in two DMDAs"); 1380 PetscCheck(stc == stf, PetscObjectComm((PetscObject)daf), PETSC_ERR_ARG_INCOMP, "Stencil type different in two DMDAs"); 1381 1382 PetscCheck(Mf >= Mc, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coarse grid has more points than fine grid, Mc %" PetscInt_FMT ", Mf %" PetscInt_FMT, Mc, Mf); 1383 PetscCheck(Nf >= Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coarse grid has more points than fine grid, Nc %" PetscInt_FMT ", Nf %" PetscInt_FMT, Nc, Nf); 1384 PetscCheck(Pf >= Pc, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Coarse grid has more points than fine grid, Pc %" PetscInt_FMT ", Pf %" PetscInt_FMT, Pc, Pf); 1385 1386 if (Pc < 0) Pc = 1; 1387 if (Pf < 0) Pf = 1; 1388 if (Nc < 0) Nc = 1; 1389 if (Nf < 0) Nf = 1; 1390 1391 PetscCall(DMDAGetCorners(daf, &i_start, &j_start, &l_start, &m_f, &n_f, &p_f)); 1392 PetscCall(DMDAGetGhostCorners(daf, &i_start_ghost, &j_start_ghost, &l_start_ghost, &m_ghost, &n_ghost, &p_ghost)); 1393 1394 PetscCall(DMGetLocalToGlobalMapping(daf, <ogmf)); 1395 PetscCall(ISLocalToGlobalMappingGetIndices(ltogmf, &idx_f)); 1396 1397 PetscCall(DMDAGetCorners(dac, &i_start_c, &j_start_c, &l_start_c, &m_c, &n_c, &p_c)); 1398 PetscCall(DMDAGetGhostCorners(dac, &i_start_ghost_c, &j_start_ghost_c, &l_start_ghost_c, &m_ghost_c, &n_ghost_c, &p_ghost_c)); 1399 1400 PetscCall(DMGetLocalToGlobalMapping(dac, <ogmc)); 1401 PetscCall(ISLocalToGlobalMappingGetIndices(ltogmc, &idx_c)); 1402 1403 /* 1404 Basic idea is as follows. Here's a 2D example, suppose r_x, r_y are the ratios 1405 for dimension 1 and 2 respectively. 1406 Let (i,j) be a coarse grid node. All the fine grid nodes between r_x*i and r_x*(i+1) 1407 and r_y*j and r_y*(j+1) will be grouped into the same coarse grid agregate. 1408 Each specific dof on the fine grid is mapped to one dof on the coarse grid. 1409 */ 1410 1411 max_agg_size = (Mf / Mc + 1) * (Nf / Nc + 1) * (Pf / Pc + 1); 1412 1413 /* create the matrix that will contain the restriction operator */ 1414 PetscCall(MatCreateAIJ(PetscObjectComm((PetscObject)daf), m_c * n_c * p_c * dofc, m_f * n_f * p_f * doff, Mc * Nc * Pc * dofc, Mf * Nf * Pf * doff, max_agg_size, NULL, max_agg_size, NULL, rest)); 1415 1416 /* store nodes in the fine grid here */ 1417 PetscCall(PetscMalloc2(max_agg_size, &one_vec, max_agg_size, &fine_nodes)); 1418 for (i = 0; i < max_agg_size; i++) one_vec[i] = 1.0; 1419 1420 /* loop over all coarse nodes */ 1421 for (l_c = l_start_c; l_c < l_start_c + p_c; l_c++) { 1422 for (j_c = j_start_c; j_c < j_start_c + n_c; j_c++) { 1423 for (i_c = i_start_c; i_c < i_start_c + m_c; i_c++) { 1424 for (d = 0; d < dofc; d++) { 1425 /* convert to local "natural" numbering and then to PETSc global numbering */ 1426 a = idx_c[dofc * (m_ghost_c * n_ghost_c * (l_c - l_start_ghost_c) + m_ghost_c * (j_c - j_start_ghost_c) + (i_c - i_start_ghost_c))] + d; 1427 1428 fn_idx = 0; 1429 /* Corresponding fine points are all points (i_f, j_f, l_f) such that 1430 i_c*Mf/Mc <= i_f < (i_c+1)*Mf/Mc 1431 (same for other dimensions) 1432 */ 1433 for (l = l_c * Pf / Pc; l < PetscMin((l_c + 1) * Pf / Pc, Pf); l++) { 1434 for (j = j_c * Nf / Nc; j < PetscMin((j_c + 1) * Nf / Nc, Nf); j++) { 1435 for (i = i_c * Mf / Mc; i < PetscMin((i_c + 1) * Mf / Mc, Mf); i++) { 1436 fine_nodes[fn_idx] = idx_f[doff * (m_ghost * n_ghost * (l - l_start_ghost) + m_ghost * (j - j_start_ghost) + (i - i_start_ghost))] + d; 1437 fn_idx++; 1438 } 1439 } 1440 } 1441 /* add all these points to one aggregate */ 1442 PetscCall(MatSetValues(*rest, 1, &a, fn_idx, fine_nodes, one_vec, INSERT_VALUES)); 1443 } 1444 } 1445 } 1446 } 1447 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogmf, &idx_f)); 1448 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogmc, &idx_c)); 1449 PetscCall(PetscFree2(one_vec, fine_nodes)); 1450 PetscCall(MatAssemblyBegin(*rest, MAT_FINAL_ASSEMBLY)); 1451 PetscCall(MatAssemblyEnd(*rest, MAT_FINAL_ASSEMBLY)); 1452 PetscFunctionReturn(0); 1453 } 1454