Actual source code: mgadapt.c

  1: #include <petsc/private/pcmgimpl.h>
  2: #include <petscdm.h>

  4: static PetscErrorCode xfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, void *ctx)
  5: {
  6:   PetscInt k = *((PetscInt *)ctx), c;

  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;

 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;

 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;

 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;

 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;

 43:   for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI * (k + 1) * coords[2]);
 44:   return PETSC_SUCCESS;
 45: }

 47: PetscErrorCode DMSetBasisFunction_Internal(PetscInt Nf, PetscBool usePoly, PetscInt dir, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))
 48: {
 49:   PetscInt f;

 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: }

 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;

 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: }

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: }

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: }

132: /*
133:   PCMGComputeCoarseSpace_Internal - Compute vectors on level l that must be accurately interpolated.

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`

141:   Output Parameter:
142: . space  - The space which must be accurately interpolated.

144:   Level: developer

146:   Note: This space is normally used to adapt the interpolator. If Nc is negative, an adaptive choice can be made.

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;

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: }

191: /*
192:   PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level l

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

202:   Level: developer

204:   Note: This routine resets the interpolation and restriction for level l.

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;

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: }

229: /*
230:   PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted

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

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: }