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: SETERRQ(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: SETERRQ(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;

 79:   DMGetCoordinateDim(dm, &dim);
 80:   DMGetNumFields(dm, &Nf);
 82:   PetscMalloc2(Nf, &funcs, Nf, &ctxs);
 83:   if (!*coarseSpace) PetscCalloc1(Nc, coarseSpace);
 84:   for (k = 0; k < Nc/dim; ++k) {
 85:     for (f = 0; f < Nf; ++f) {ctxs[f] = &k;}
 86:     for (d = 0; d < dim; ++d) {
 87:       if (!(*coarseSpace)[k*dim+d]) DMCreateGlobalVector(dm, &(*coarseSpace)[k*dim+d]);
 88:       DMSetBasisFunction_Internal(Nf, poly, d, funcs);
 89:       DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, (*coarseSpace)[k*dim+d]);
 90:     }
 91:   }
 92:   PetscFree2(funcs, ctxs);
 93:   return 0;
 94: }

 96: static PetscErrorCode PCMGCreateCoarseSpace_Polynomial(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace)
 97: {
 98:   PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_POLYNOMIAL, dm, ksp, Nc, initialGuess, coarseSpace);
 99:   return 0;
100: }

102: PetscErrorCode PCMGCreateCoarseSpace_Harmonic(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace)
103: {
104:   PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_HARMONIC, dm, ksp, Nc, initialGuess, coarseSpace);
105:   return 0;
106: }

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

111:   Input Parameters:
112: + pc     - The PCMG
113: . l      - The level
114: . Nc     - The size of the space (number of vectors)
115: - cspace - The space from level l-1, or NULL

117:   Output Parameter:
118: . space  - The space which must be accurately interpolated.

120:   Level: developer

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

124: .seealso: PCMGAdaptInterpolator_Private()
125: */
126: PetscErrorCode PCMGComputeCoarseSpace_Internal(PC pc, PetscInt l, PCMGCoarseSpaceType cstype, PetscInt Nc, const Vec cspace[], Vec *space[])
127: {
128:   PetscErrorCode (*coarseConstructor)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec*[]);
129:   DM             dm;
130:   KSP            smooth;

132:   switch (cstype) {
133:   case PCMG_POLYNOMIAL: coarseConstructor = &PCMGCreateCoarseSpace_Polynomial;break;
134:   case PCMG_HARMONIC:   coarseConstructor = &PCMGCreateCoarseSpace_Harmonic;break;
135:   case PCMG_EIGENVECTOR:
136:     if (l > 0) PCMGGetCoarseSpaceConstructor("BAMG_MEV", &coarseConstructor);
137:     else       PCMGGetCoarseSpaceConstructor("BAMG_EV", &coarseConstructor);
138:     break;
139:   case PCMG_GENERALIZED_EIGENVECTOR:
140:     if (l > 0) PCMGGetCoarseSpaceConstructor("BAMG_MGEV", &coarseConstructor);
141:     else       PCMGGetCoarseSpaceConstructor("BAMG_GEV", &coarseConstructor);
142:     break;
143:   default: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle coarse space type %D", cstype);
144:   }
145:   PCMGGetSmoother(pc, l, &smooth);
146:   KSPGetDM(smooth, &dm);
147:   (*coarseConstructor)(pc, l, dm, smooth, Nc, cspace, space);
148:   return 0;
149: }

151: /*
152:   PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level 1

154:   Input Parameters:
155: + pc      - The PCMG
156: . l       - The level l
157: . csmooth - The (coarse) smoother for level l-1
158: . fsmooth - The (fine) smoother for level l
159: . Nc      - The size of the subspace used for adaptation
160: . cspace  - The (coarse) vectors in the subspace for level l-1
161: - fspace  - The (fine) vectors in the subspace for level l

163:   Level: developer

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

167: .seealso: PCMGComputeCoarseSpace_Private()
168: */
169: PetscErrorCode PCMGAdaptInterpolator_Internal(PC pc, PetscInt l, KSP csmooth, KSP fsmooth, PetscInt Nc, Vec cspace[], Vec fspace[])
170: {
171:   PC_MG         *mg = (PC_MG *) pc->data;
172:   DM             dm, cdm;
173:   Mat            Interp, InterpAdapt;

175:   /* There is no interpolator for the coarse level */
176:   if (!l) return 0;
177:   KSPGetDM(csmooth, &cdm);
178:   KSPGetDM(fsmooth, &dm);
179:   PCMGGetInterpolation(pc, l, &Interp);

181:   DMAdaptInterpolator(cdm, dm, Interp, fsmooth, Nc, fspace, cspace, &InterpAdapt, pc);
182:   if (mg->mespMonitor) DMCheckInterpolator(dm, InterpAdapt, Nc, cspace, fspace, 0.5/* PETSC_SMALL */);
183:   PCMGSetInterpolation(pc, l, InterpAdapt);
184:   PCMGSetRestriction(pc, l, InterpAdapt);
185:   MatDestroy(&InterpAdapt);
186:   return 0;
187: }

189: /*
190:   PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted

192:   Input Parameters:
193: + pc - The PCMG
194: - l  - The level l

196:   Level: developer

198:   Note: This routine recomputes the Galerkin triple product for the operator on level l.
199: */
200: PetscErrorCode PCMGRecomputeLevelOperators_Internal(PC pc, PetscInt l)
201: {
202:   Mat              fA, fB;                   /* The system and preconditioning operators on level l+1 */
203:   Mat              A,  B;                    /* The system and preconditioning operators on level l */
204:   Mat              Interp, Restrc;           /* The interpolation operator from level l to l+1, and restriction operator from level l+1 to l */
205:   KSP              smooth, fsmooth;          /* The smoothers on levels l and l+1 */
206:   PCMGGalerkinType galerkin;                 /* The Galerkin projection flag */
207:   MatReuse         reuse = MAT_REUSE_MATRIX; /* The matrices are always assumed to be present already */
208:   PetscBool        doA   = PETSC_FALSE;      /* Updates the system operator */
209:   PetscBool        doB   = PETSC_FALSE;      /* Updates the preconditioning operator (A == B, then update B) */
210:   PetscInt         n;                        /* The number of multigrid levels */

212:   PCMGGetGalerkin(pc, &galerkin);
213:   if (galerkin >= PC_MG_GALERKIN_NONE) return 0;
214:   PCMGGetLevels(pc, &n);
215:   /* Do not recompute operator for the finest grid */
216:   if (l == n-1) return 0;
217:   PCMGGetSmoother(pc, l,   &smooth);
218:   KSPGetOperators(smooth, &A, &B);
219:   PCMGGetSmoother(pc, l+1, &fsmooth);
220:   KSPGetOperators(fsmooth, &fA, &fB);
221:   PCMGGetInterpolation(pc, l+1, &Interp);
222:   PCMGGetRestriction(pc, l+1, &Restrc);
223:   if ((galerkin == PC_MG_GALERKIN_PMAT) ||  (galerkin == PC_MG_GALERKIN_BOTH))                doB = PETSC_TRUE;
224:   if ((galerkin == PC_MG_GALERKIN_MAT)  || ((galerkin == PC_MG_GALERKIN_BOTH) && (fA != fB))) doA = PETSC_TRUE;
225:   if (doA) MatGalerkin(Restrc, fA, Interp, reuse, 1.0, &A);
226:   if (doB) MatGalerkin(Restrc, fB, Interp, reuse, 1.0, &B);
227:   return 0;
228: }