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