1 #include <petsc/private/pcmgimpl.h> /*I "petscksp.h" I*/
2 #include <petscdm.h>
3
xfunc(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)4 static PetscErrorCode xfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
5 {
6 PetscInt k = *((PetscInt *)ctx), c;
7
8 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[0], k);
9 return PETSC_SUCCESS;
10 }
yfunc(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)11 static PetscErrorCode yfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
12 {
13 PetscInt k = *((PetscInt *)ctx), c;
14
15 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[1], k);
16 return PETSC_SUCCESS;
17 }
zfunc(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)18 static PetscErrorCode zfunc(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
19 {
20 PetscInt k = *((PetscInt *)ctx), c;
21
22 for (c = 0; c < Nc; ++c) u[c] = PetscPowRealInt(coords[2], k);
23 return PETSC_SUCCESS;
24 }
xsin(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)25 static PetscErrorCode xsin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
26 {
27 PetscInt k = *((PetscInt *)ctx), c;
28
29 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI * (k + 1) * coords[0]);
30 return PETSC_SUCCESS;
31 }
ysin(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)32 static PetscErrorCode ysin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
33 {
34 PetscInt k = *((PetscInt *)ctx), c;
35
36 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI * (k + 1) * coords[1]);
37 return PETSC_SUCCESS;
38 }
zsin(PetscInt dim,PetscReal time,const PetscReal coords[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)39 static PetscErrorCode zsin(PetscInt dim, PetscReal time, const PetscReal coords[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
40 {
41 PetscInt k = *((PetscInt *)ctx), c;
42
43 for (c = 0; c < Nc; ++c) u[c] = PetscSinReal(PETSC_PI * (k + 1) * coords[2]);
44 return PETSC_SUCCESS;
45 }
46
DMSetBasisFunction_Internal(PetscInt Nf,PetscBool usePoly,PetscInt dir,PetscErrorCode (** funcs)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar *,void *))47 PetscErrorCode DMSetBasisFunction_Internal(PetscInt Nf, PetscBool usePoly, PetscInt dir, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *))
48 {
49 PetscInt f;
50
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 }
85
PCMGCreateCoarseSpaceDefault_Private(PC pc,PetscInt level,PCMGCoarseSpaceType cstype,DM dm,KSP ksp,PetscInt Nc,Mat initialGuess,Mat * coarseSpace)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;
93
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 }
117
PCMGCreateCoarseSpace_Polynomial(PC pc,PetscInt level,DM dm,KSP ksp,PetscInt Nc,Mat initialGuess,Mat * coarseSpace)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 }
124
PCMGCreateCoarseSpace_Harmonic(PC pc,PetscInt level,DM dm,KSP ksp,PetscInt Nc,Mat initialGuess,Mat * coarseSpace)125 static 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 }
131
132 /*
133 PCMGComputeCoarseSpace_Internal - Compute vectors on level l that must be accurately interpolated.
134
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`
140
141 Output Parameter:
142 . space - The space which must be accurately interpolated.
143
144 Level: developer
145
146 Note: This space is normally used to adapt the interpolator. If Nc is negative, an adaptive choice can be made.
147
148 .seealso: [](ch_ksp), `PCMGAdaptInterpolator_Private()`
149 */
PCMGComputeCoarseSpace_Internal(PC pc,PetscInt l,PCMGCoarseSpaceType cstype,PetscInt Nc,Mat cspace,Mat * space)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;
155
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 }
190
191 /*
192 PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level l
193
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
201
202 Level: developer
203
204 Note: This routine resets the interpolation and restriction for level l.
205
206 .seealso: [](ch_ksp), `PCMGComputeCoarseSpace_Private()`
207 */
PCMGAdaptInterpolator_Internal(PC pc,PetscInt l,KSP csmooth,KSP fsmooth,Mat cspace,Mat fspace)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;
213
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 }
228
229 /*
230 PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted
231
232 Note: This routine recomputes the Galerkin triple product for the operator on level l.
233 */
PCMGRecomputeLevelOperators_Internal(PC pc,PetscInt l)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 */
245
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 }
264