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