1 #include <petsc/private/pcmgimpl.h> /*I "petscksp.h" I*/ 2 #include <petscdm.h> 3 4 static PetscErrorCode xfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 5 { 6 PetscInt k = *((PetscInt *)ctx), c; 7 8 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[0], k); 9 return PETSC_SUCCESS; 10 } 11 static PetscErrorCode yfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 12 { 13 PetscInt k = *((PetscInt *)ctx), c; 14 15 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[1], k); 16 return PETSC_SUCCESS; 17 } 18 static PetscErrorCode zfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 19 { 20 PetscInt k = *((PetscInt *)ctx), c; 21 22 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[2], k); 23 return PETSC_SUCCESS; 24 } 25 static PetscErrorCode xsin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 26 { 27 PetscInt k = *((PetscInt *)ctx), c; 28 29 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI * (k + 1) * coords[0]); 30 return PETSC_SUCCESS; 31 } 32 static PetscErrorCode ysin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 33 { 34 PetscInt k = *((PetscInt *)ctx), c; 35 36 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI * (k + 1) * coords[1]); 37 return PETSC_SUCCESS; 38 } 39 static PetscErrorCode zsin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx) 40 { 41 PetscInt k = *((PetscInt *)ctx), c; 42 43 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI * (k + 1) * coords[2]); 44 return PETSC_SUCCESS; 45 } 46 47 PetscErrorCode DMSetBasisFunction_Internal(PetscInt Nf, PetscBool usePoly, PetscInt dir, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) 48 { 49 PetscInt f; 50 51 PetscFunctionBeginUser; 52 for (f = 0; f < Nf; ++f) { 53 if (usePoly) { 54 switch (dir) { 55 case 0: 56 funcs[f] = xfunc; 57 break; 58 case 1: 59 funcs[f] = yfunc; 60 break; 61 case 2: 62 funcs[f] = zfunc; 63 break; 64 default: 65 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No function for direction %" PetscInt_FMT, dir); 66 } 67 } else { 68 switch (dir) { 69 case 0: 70 funcs[f] = xsin; 71 break; 72 case 1: 73 funcs[f] = ysin; 74 break; 75 case 2: 76 funcs[f] = zsin; 77 break; 78 default: 79 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No function for direction %" PetscInt_FMT, dir); 80 } 81 } 82 } 83 PetscFunctionReturn(PETSC_SUCCESS); 84 } 85 86 static PetscErrorCode PCMGCreateCoarseSpaceDefault_Private(PC pc, PetscInt level, PCMGCoarseSpaceType cstype, DM dm, KSP ksp, PetscInt Nc, Mat initialGuess, Mat *coarseSpace) 87 { 88 PetscBool poly = cstype == PCMG_ADAPT_POLYNOMIAL ? PETSC_TRUE : PETSC_FALSE; 89 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *); 90 void **ctxs; 91 PetscInt dim, d, Nf, f, k, m, M; 92 Vec tmp; 93 94 PetscFunctionBegin; 95 Nc = Nc < 0 ? 6 : Nc; 96 PetscCall(DMGetCoordinateDim(dm, &dim)); 97 PetscCall(DMGetNumFields(dm, &Nf)); 98 PetscCheck(Nc % dim == 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "The number of coarse vectors %" PetscInt_FMT " must be divisible by the dimension %" PetscInt_FMT, Nc, dim); 99 PetscCall(PetscMalloc2(Nf, &funcs, Nf, &ctxs)); 100 PetscCall(DMGetGlobalVector(dm, &tmp)); 101 PetscCall(VecGetSize(tmp, &M)); 102 PetscCall(VecGetLocalSize(tmp, &m)); 103 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)pc), m, PETSC_DECIDE, M, Nc, NULL, coarseSpace)); 104 PetscCall(DMRestoreGlobalVector(dm, &tmp)); 105 for (k = 0; k < Nc / dim; ++k) { 106 for (f = 0; f < Nf; ++f) ctxs[f] = &k; 107 for (d = 0; d < dim; ++d) { 108 PetscCall(MatDenseGetColumnVecWrite(*coarseSpace, k * dim + d, &tmp)); 109 PetscCall(DMSetBasisFunction_Internal(Nf, poly, d, funcs)); 110 PetscCall(DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, tmp)); 111 PetscCall(MatDenseRestoreColumnVecWrite(*coarseSpace, k * dim + d, &tmp)); 112 } 113 } 114 PetscCall(PetscFree2(funcs, ctxs)); 115 PetscFunctionReturn(PETSC_SUCCESS); 116 } 117 118 static PetscErrorCode PCMGCreateCoarseSpace_Polynomial(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, Mat initialGuess, Mat *coarseSpace) 119 { 120 PetscFunctionBegin; 121 PetscCall(PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_ADAPT_POLYNOMIAL, dm, ksp, Nc, initialGuess, coarseSpace)); 122 PetscFunctionReturn(PETSC_SUCCESS); 123 } 124 125 static PetscErrorCode PCMGCreateCoarseSpace_Harmonic(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, Mat initialGuess, Mat *coarseSpace) 126 { 127 PetscFunctionBegin; 128 PetscCall(PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_ADAPT_HARMONIC, dm, ksp, Nc, initialGuess, coarseSpace)); 129 PetscFunctionReturn(PETSC_SUCCESS); 130 } 131 132 /* 133 PCMGComputeCoarseSpace_Internal - Compute vectors on level l that must be accurately interpolated. 134 135 Input Parameters: 136 + pc - The `PCMG` 137 . l - The level 138 . Nc - The number of vectors requested 139 - cspace - The initial guess for the space, or `NULL` 140 141 Output Parameter: 142 . space - The space which must be accurately interpolated. 143 144 Level: developer 145 146 Note: This space is normally used to adapt the interpolator. If Nc is negative, an adaptive choice can be made. 147 148 .seealso: [](ch_ksp), `PCMGAdaptInterpolator_Private()` 149 */ 150 PetscErrorCode PCMGComputeCoarseSpace_Internal(PC pc, PetscInt l, PCMGCoarseSpaceType cstype, PetscInt Nc, Mat cspace, Mat *space) 151 { 152 PetscErrorCode (*coarseConstructor)(PC, PetscInt, DM, KSP, PetscInt, Mat, Mat *) = NULL; 153 DM dm; 154 KSP smooth; 155 156 PetscFunctionBegin; 157 *space = NULL; 158 switch (cstype) { 159 case PCMG_ADAPT_POLYNOMIAL: 160 coarseConstructor = &PCMGCreateCoarseSpace_Polynomial; 161 break; 162 case PCMG_ADAPT_HARMONIC: 163 coarseConstructor = &PCMGCreateCoarseSpace_Harmonic; 164 break; 165 case PCMG_ADAPT_EIGENVECTOR: 166 Nc = Nc < 0 ? 6 : Nc; 167 if (l > 0) PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_MEV", &coarseConstructor)); 168 else PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_EV", &coarseConstructor)); 169 break; 170 case PCMG_ADAPT_GENERALIZED_EIGENVECTOR: 171 Nc = Nc < 0 ? 6 : Nc; 172 if (l > 0) PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_MGEV", &coarseConstructor)); 173 else PetscCall(PCMGGetCoarseSpaceConstructor("BAMG_GEV", &coarseConstructor)); 174 break; 175 case PCMG_ADAPT_GDSW: 176 coarseConstructor = &PCMGGDSWCreateCoarseSpace_Private; 177 break; 178 case PCMG_ADAPT_NONE: 179 break; 180 default: 181 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot handle coarse space type %d", cstype); 182 } 183 if (coarseConstructor) { 184 PetscCall(PCMGGetSmoother(pc, l, &smooth)); 185 PetscCall(KSPGetDM(smooth, &dm)); 186 PetscCall((*coarseConstructor)(pc, l, dm, smooth, Nc, cspace, space)); 187 } 188 PetscFunctionReturn(PETSC_SUCCESS); 189 } 190 191 /* 192 PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level l 193 194 Input Parameters: 195 + pc - The PCMG 196 . l - The level l 197 . csmooth - The (coarse) smoother for level l-1 198 . fsmooth - The (fine) smoother for level l 199 . cspace - The (coarse) vectors in the subspace for level l-1 200 - fspace - The (fine) vectors in the subspace for level l 201 202 Level: developer 203 204 Note: This routine resets the interpolation and restriction for level l. 205 206 .seealso: [](ch_ksp), `PCMGComputeCoarseSpace_Private()` 207 */ 208 PetscErrorCode PCMGAdaptInterpolator_Internal(PC pc, PetscInt l, KSP csmooth, KSP fsmooth, Mat cspace, Mat fspace) 209 { 210 PC_MG *mg = (PC_MG *)pc->data; 211 DM dm, cdm; 212 Mat Interp, InterpAdapt; 213 214 PetscFunctionBegin; 215 /* There is no interpolator for the coarse level */ 216 if (!l) PetscFunctionReturn(PETSC_SUCCESS); 217 PetscCall(KSPGetDM(csmooth, &cdm)); 218 PetscCall(KSPGetDM(fsmooth, &dm)); 219 PetscCall(PCMGGetInterpolation(pc, l, &Interp)); 220 if (Interp == fspace && !cspace) PetscFunctionReturn(PETSC_SUCCESS); 221 PetscCall(DMAdaptInterpolator(cdm, dm, Interp, fsmooth, fspace, cspace, &InterpAdapt, pc)); 222 if (mg->mespMonitor) PetscCall(DMCheckInterpolator(dm, InterpAdapt, cspace, fspace, 0.5 /* PETSC_SMALL */)); 223 PetscCall(PCMGSetInterpolation(pc, l, InterpAdapt)); 224 PetscCall(PCMGSetRestriction(pc, l, InterpAdapt)); /* MATT: Remove????? */ 225 PetscCall(MatDestroy(&InterpAdapt)); 226 PetscFunctionReturn(PETSC_SUCCESS); 227 } 228 229 /* 230 PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted 231 232 Note: This routine recomputes the Galerkin triple product for the operator on level l. 233 */ 234 PetscErrorCode PCMGRecomputeLevelOperators_Internal(PC pc, PetscInt l) 235 { 236 Mat fA, fB; /* The system and preconditioning operators on level l+1 */ 237 Mat A, B; /* The system and preconditioning operators on level l */ 238 Mat Interp, Restrc; /* The interpolation operator from level l to l+1, and restriction operator from level l+1 to l */ 239 KSP smooth, fsmooth; /* The smoothers on levels l and l+1 */ 240 PCMGGalerkinType galerkin; /* The Galerkin projection flag */ 241 MatReuse reuse = MAT_REUSE_MATRIX; /* The matrices are always assumed to be present already */ 242 PetscBool doA = PETSC_FALSE; /* Updates the system operator */ 243 PetscBool doB = PETSC_FALSE; /* Updates the preconditioning operator (A == B, then update B) */ 244 PetscInt n; /* The number of multigrid levels */ 245 246 PetscFunctionBegin; 247 PetscCall(PCMGGetGalerkin(pc, &galerkin)); 248 if (galerkin >= PC_MG_GALERKIN_NONE) PetscFunctionReturn(PETSC_SUCCESS); 249 PetscCall(PCMGGetLevels(pc, &n)); 250 /* Do not recompute operator for the finest grid */ 251 if (l == n - 1) PetscFunctionReturn(PETSC_SUCCESS); 252 PetscCall(PCMGGetSmoother(pc, l, &smooth)); 253 PetscCall(KSPGetOperators(smooth, &A, &B)); 254 PetscCall(PCMGGetSmoother(pc, l + 1, &fsmooth)); 255 PetscCall(KSPGetOperators(fsmooth, &fA, &fB)); 256 PetscCall(PCMGGetInterpolation(pc, l + 1, &Interp)); 257 PetscCall(PCMGGetRestriction(pc, l + 1, &Restrc)); 258 if ((galerkin == PC_MG_GALERKIN_PMAT) || (galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE; 259 if ((galerkin == PC_MG_GALERKIN_MAT) || ((galerkin == PC_MG_GALERKIN_BOTH) && (fA != fB))) doA = PETSC_TRUE; 260 if (doA) PetscCall(MatGalerkin(Restrc, fA, Interp, reuse, 1.0, &A)); 261 if (doB) PetscCall(MatGalerkin(Restrc, fB, Interp, reuse, 1.0, &B)); 262 PetscFunctionReturn(PETSC_SUCCESS); 263 } 264