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