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