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 0;
     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 0;
     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 0;
     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 0;
     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 0;
     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 0;
     45: }

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

     52:   for (f = 0; f < Nf; ++f) {
     53:     if (usePoly) {
     54:       switch (dir) {
     55:       case 0: funcs[f] = xfunc;break;
     56:       case 1: funcs[f] = yfunc;break;
     57:       case 2: funcs[f] = zfunc;break;
     58:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No function for direction %D", dir);
     59:       }
     60:     } else {
     61:       switch (dir) {
     62:       case 0: funcs[f] = xsin;break;
     63:       case 1: funcs[f] = ysin;break;
     64:       case 2: funcs[f] = zsin;break;
     65:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No function for direction %D", dir);
     66:       }
     67:     }
     68:   }
     69:   return(0);
     70: }

     72: static PetscErrorCode PCMGCreateCoarseSpaceDefault_Private(PC pc, PetscInt level, PCMGCoarseSpaceType cstype, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace)
     73: {
     74:   PetscBool         poly = cstype == PCMG_POLYNOMIAL ? PETSC_TRUE : PETSC_FALSE;
     75:   PetscErrorCode (**funcs)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar*,void*);
     76:   void            **ctxs;
     77:   PetscInt          dim, d, Nf, f, k;
     78:   PetscErrorCode    ierr;

     81:   DMGetCoordinateDim(dm, &dim);
     82:   DMGetNumFields(dm, &Nf);
     83:   if (Nc % dim) SETERRQ2(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "The number of coarse vectors %D must be divisible by the dimension %D", Nc, dim);
     84:   PetscMalloc2(Nf, &funcs, Nf, &ctxs);
     85:   if (!*coarseSpace) {PetscCalloc1(Nc, coarseSpace);}
     86:   for (k = 0; k < Nc/dim; ++k) {
     87:     for (f = 0; f < Nf; ++f) {ctxs[f] = &k;}
     88:     for (d = 0; d < dim; ++d) {
     89:       if (!(*coarseSpace)[k*dim+d]) {DMCreateGlobalVector(dm, &(*coarseSpace)[k*dim+d]);}
     90:       DMSetBasisFunction_Internal(Nf, poly, d, funcs);
     91:       DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, (*coarseSpace)[k*dim+d]);
     92:     }
     93:   }
     94:   PetscFree2(funcs, ctxs);
     95:   return(0);
     96: }

     98: static PetscErrorCode PCMGCreateCoarseSpace_Polynomial(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace)
     99: {

    103:   PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_POLYNOMIAL, dm, ksp, Nc, initialGuess, coarseSpace);
    104:   return(0);
    105: }

    107: PetscErrorCode PCMGCreateCoarseSpace_Harmonic(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace)
    108: {

    112:   PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_HARMONIC, dm, ksp, Nc, initialGuess, coarseSpace);
    113:   return(0);
    114: }

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

    119:   Input Parameters:
    120: + pc     - The PCMG
    121: . l      - The level
    122: . Nc     - The size of the space (number of vectors)
    123: - cspace - The space from level l-1, or NULL

    125:   Output Parameter:
    126: . space  - The space which must be accurately interpolated.

    128:   Level: developer

    130:   Note: This space is normally used to adapt the interpolator.

    132: .seealso: PCMGAdaptInterpolator_Private()
    133: */
    134: PetscErrorCode PCMGComputeCoarseSpace_Internal(PC pc, PetscInt l, PCMGCoarseSpaceType cstype, PetscInt Nc, const Vec cspace[], Vec *space[])
    135: {
    136:   PetscErrorCode (*coarseConstructor)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec*[]);
    137:   DM             dm;
    138:   KSP            smooth;

    142:   switch (cstype) {
    143:   case PCMG_POLYNOMIAL: coarseConstructor = &PCMGCreateCoarseSpace_Polynomial;break;
    144:   case PCMG_HARMONIC:   coarseConstructor = &PCMGCreateCoarseSpace_Harmonic;break;
    145:   case PCMG_EIGENVECTOR:
    146:     if (l > 0) {PCMGGetCoarseSpaceConstructor("BAMG_MEV", &coarseConstructor);}
    147:     else       {PCMGGetCoarseSpaceConstructor("BAMG_EV", &coarseConstructor);}
    148:     break;
    149:   case PCMG_GENERALIZED_EIGENVECTOR:
    150:     if (l > 0) {PCMGGetCoarseSpaceConstructor("BAMG_MGEV", &coarseConstructor);}
    151:     else       {PCMGGetCoarseSpaceConstructor("BAMG_GEV", &coarseConstructor);}
    152:     break;
    153:   default: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle coarse space type %D", cstype);
    154:   }
    155:   PCMGGetSmoother(pc, l, &smooth);
    156:   KSPGetDM(smooth, &dm);
    157:   (*coarseConstructor)(pc, l, dm, smooth, Nc, cspace, space);
    158:   return(0);
    159: }

    161: /*
    162:   PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level 1

    164:   Input Parameters:
    165: + pc      - The PCMG
    166: . l       - The level l
    167: . csmooth - The (coarse) smoother for level l-1
    168: . fsmooth - The (fine) smoother for level l
    169: . Nc      - The size of the subspace used for adaptation
    170: . cspace  - The (coarse) vectors in the subspace for level l-1
    171: - fspace  - The (fine) vectors in the subspace for level l

    173:   Level: developer

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

    177: .seealso: PCMGComputeCoarseSpace_Private()
    178: */
    179: PetscErrorCode PCMGAdaptInterpolator_Internal(PC pc, PetscInt l, KSP csmooth, KSP fsmooth, PetscInt Nc, Vec cspace[], Vec fspace[])
    180: {
    181:   PC_MG         *mg = (PC_MG *) pc->data;
    182:   DM             dm, cdm;
    183:   Mat            Interp, InterpAdapt;

    187:   /* There is no interpolator for the coarse level */
    188:   if (!l) return(0);
    189:   KSPGetDM(csmooth, &cdm);
    190:   KSPGetDM(fsmooth, &dm);
    191:   PCMGGetInterpolation(pc, l, &Interp);

    193:   DMAdaptInterpolator(cdm, dm, Interp, fsmooth, Nc, fspace, cspace, &InterpAdapt, pc);
    194:   if (mg->mespMonitor) {DMCheckInterpolator(dm, InterpAdapt, Nc, cspace, fspace, 0.5/* PETSC_SMALL */);}
    195:   PCMGSetInterpolation(pc, l, InterpAdapt);
    196:   PCMGSetRestriction(pc, l, InterpAdapt);
    197:   MatDestroy(&InterpAdapt);
    198:   return(0);
    199: }

    201: /*
    202:   PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted

    204:   Input Parameters:
    205: + pc - The PCMG
    206: - l  - The level l

    208:   Level: developer

    210:   Note: This routine recomputes the Galerkin triple product for the operator on level l.
    211: */
    212: PetscErrorCode PCMGRecomputeLevelOperators_Internal(PC pc, PetscInt l)
    213: {
    214:   Mat              fA, fB;                   /* The system and preconditioning operators on level l+1 */
    215:   Mat              A,  B;                    /* The system and preconditioning operators on level l */
    216:   Mat              Interp, Restrc;           /* The interpolation operator from level l to l+1, and restriction operator from level l+1 to l */
    217:   KSP              smooth, fsmooth;          /* The smoothers on levels l and l+1 */
    218:   PCMGGalerkinType galerkin;                 /* The Galerkin projection flag */
    219:   MatReuse         reuse = MAT_REUSE_MATRIX; /* The matrices are always assumed to be present already */
    220:   PetscBool        doA   = PETSC_FALSE;      /* Updates the system operator */
    221:   PetscBool        doB   = PETSC_FALSE;      /* Updates the preconditioning operator (A == B, then update B) */
    222:   PetscInt         n;                        /* The number of multigrid levels */
    223:   PetscErrorCode   ierr;

    226:   PCMGGetGalerkin(pc, &galerkin);
    227:   if (galerkin >= PC_MG_GALERKIN_NONE) return(0);
    228:   PCMGGetLevels(pc, &n);
    229:   /* Do not recompute operator for the finest grid */
    230:   if (l == n-1) return(0);
    231:   PCMGGetSmoother(pc, l,   &smooth);
    232:   KSPGetOperators(smooth, &A, &B);
    233:   PCMGGetSmoother(pc, l+1, &fsmooth);
    234:   KSPGetOperators(fsmooth, &fA, &fB);
    235:   PCMGGetInterpolation(pc, l+1, &Interp);
    236:   PCMGGetRestriction(pc, l+1, &Restrc);
    237:   if ((galerkin == PC_MG_GALERKIN_PMAT) ||  (galerkin == PC_MG_GALERKIN_BOTH))                doB = PETSC_TRUE;
    238:   if ((galerkin == PC_MG_GALERKIN_MAT)  || ((galerkin == PC_MG_GALERKIN_BOTH) && (fA != fB))) doA = PETSC_TRUE;
    239:   if (doA) {MatGalerkin(Restrc, fA, Interp, reuse, 1.0, &A);}
    240:   if (doB) {MatGalerkin(Restrc, fB, Interp, reuse, 1.0, &B);}
    241:   return(0);
    242: }