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