1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
3 #include <petscds.h>
4 #include <petsc/private/petscimpl.h>
5 #include <petsc/private/petscfeimpl.h>
6
7 #ifdef PETSC_HAVE_LIBCEED
8 #include <petscdmceed.h>
9 #include <petscdmplexceed.h>
10 #endif
11
pressure_Private(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar p[])12 static void pressure_Private(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar p[])
13 {
14 p[0] = u[uOff[1]];
15 }
16
17 /*
18 SNESCorrectDiscretePressure_Private - Add a vector in the nullspace to make the continuum integral of the pressure field equal to zero.
19 This is normally used only to evaluate convergence rates for the pressure accurately.
20
21 Collective
22
23 Input Parameters:
24 + snes - The `SNES`
25 . pfield - The field number for pressure
26 . nullspace - The pressure nullspace
27 . u - The solution vector
28 - ctx - An optional application context
29
30 Output Parameter:
31 . u - The solution with a continuum pressure integral of zero
32
33 Level: developer
34
35 Note:
36 If int(u) = a and int(n) = b, then int(u - a/b n) = a - a/b b = 0. We assume that the nullspace is a single vector given explicitly.
37
38 .seealso: [](ch_snes), `SNESConvergedCorrectPressure()`
39 */
SNESCorrectDiscretePressure_Private(SNES snes,PetscInt pfield,MatNullSpace nullspace,Vec u,PetscCtx ctx)40 static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, PetscCtx ctx)
41 {
42 DM dm;
43 PetscDS ds;
44 const Vec *nullvecs;
45 PetscScalar pintd, *intc, *intn;
46 MPI_Comm comm;
47 PetscInt Nf, Nv;
48
49 PetscFunctionBegin;
50 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
51 PetscCall(SNESGetDM(snes, &dm));
52 PetscCheck(dm, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a SNES DM");
53 PetscCheck(nullspace, comm, PETSC_ERR_ARG_WRONG, "Cannot compute test without a Jacobian nullspace");
54 PetscCall(DMGetDS(dm, &ds));
55 PetscCall(PetscDSSetObjective(ds, pfield, pressure_Private));
56 PetscCall(MatNullSpaceGetVecs(nullspace, NULL, &Nv, &nullvecs));
57 PetscCheck(Nv == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Can only handle a single null vector for pressure, not %" PetscInt_FMT, Nv);
58 PetscCall(VecDot(nullvecs[0], u, &pintd));
59 PetscCheck(PetscAbsScalar(pintd) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Discrete integral of pressure: %g", (double)PetscRealPart(pintd));
60 PetscCall(PetscDSGetNumFields(ds, &Nf));
61 PetscCall(PetscMalloc2(Nf, &intc, Nf, &intn));
62 PetscCall(DMPlexComputeIntegralFEM(dm, nullvecs[0], intn, ctx));
63 PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
64 PetscCall(VecAXPY(u, -intc[pfield] / intn[pfield], nullvecs[0]));
65 #if defined(PETSC_USE_DEBUG)
66 PetscCall(DMPlexComputeIntegralFEM(dm, u, intc, ctx));
67 PetscCheck(PetscAbsScalar(intc[pfield]) <= PETSC_SMALL, comm, PETSC_ERR_ARG_WRONG, "Continuum integral of pressure after correction: %g", (double)PetscRealPart(intc[pfield]));
68 #endif
69 PetscCall(PetscFree2(intc, intn));
70 PetscFunctionReturn(PETSC_SUCCESS);
71 }
72
73 /*@C
74 SNESConvergedCorrectPressure - The regular `SNES` convergence test that, up on convergence, adds a vector in the nullspace
75 to make the continuum integral of the pressure field equal to zero.
76
77 Logically Collective
78
79 Input Parameters:
80 + snes - the `SNES` context
81 . it - the iteration (0 indicates before any Newton steps)
82 . xnorm - 2-norm of current iterate
83 . gnorm - 2-norm of current step
84 . f - 2-norm of function at current iterate
85 - ctx - Optional application context
86
87 Output Parameter:
88 . reason - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FUNCTION_NANORINF`
89
90 Options Database Key:
91 . -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`
92
93 Level: advanced
94
95 Notes:
96 In order to use this convergence test, you must set up several PETSc structures. First fields must be added to the `DM`, and a `PetscDS`
97 must be created with discretizations of those fields. We currently assume that the pressure field has index 1.
98 The pressure field must have a nullspace, likely created using the `DMSetNullSpaceConstructor()` interface.
99 Last we must be able to integrate the pressure over the domain, so the `DM` attached to the SNES `must` be a `DMPLEX` at this time.
100
101 Developer Note:
102 This is a total misuse of the `SNES` convergence test handling system. It should be removed. Perhaps a `SNESSetPostSolve()` could
103 be constructed to handle this process.
104
105 .seealso: [](ch_snes), `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`
106 @*/
SNESConvergedCorrectPressure(SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason * reason,PetscCtx ctx)107 PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, PetscCtx ctx)
108 {
109 PetscBool monitorIntegral = PETSC_FALSE;
110
111 PetscFunctionBegin;
112 PetscCall(SNESConvergedDefault(snes, it, xnorm, gnorm, f, reason, ctx));
113 if (monitorIntegral) {
114 Mat J;
115 Vec u;
116 MatNullSpace nullspace;
117 const Vec *nullvecs;
118 PetscScalar pintd;
119
120 PetscCall(SNESGetSolution(snes, &u));
121 PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
122 PetscCall(MatGetNullSpace(J, &nullspace));
123 PetscCall(MatNullSpaceGetVecs(nullspace, NULL, NULL, &nullvecs));
124 PetscCall(VecDot(nullvecs[0], u, &pintd));
125 PetscCall(PetscInfo(snes, "SNES: Discrete integral of pressure: %g\n", (double)PetscRealPart(pintd)));
126 }
127 if (*reason > 0) {
128 Mat J;
129 Vec u;
130 MatNullSpace nullspace;
131 PetscInt pfield = 1;
132
133 PetscCall(SNESGetSolution(snes, &u));
134 PetscCall(SNESGetJacobian(snes, &J, NULL, NULL, NULL));
135 PetscCall(MatGetNullSpace(J, &nullspace));
136 PetscCall(SNESCorrectDiscretePressure_Private(snes, pfield, nullspace, u, ctx));
137 }
138 PetscFunctionReturn(PETSC_SUCCESS);
139 }
140
DMSNESConvertPlex(DM dm,DM * plex,PetscBool copy)141 static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
142 {
143 PetscBool isPlex;
144
145 PetscFunctionBegin;
146 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
147 if (isPlex) {
148 *plex = dm;
149 PetscCall(PetscObjectReference((PetscObject)dm));
150 } else {
151 PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex));
152 if (!*plex) {
153 PetscCall(DMConvert(dm, DMPLEX, plex));
154 PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex));
155 } else {
156 PetscCall(PetscObjectReference((PetscObject)*plex));
157 }
158 if (copy) {
159 PetscCall(DMCopyDMSNES(dm, *plex));
160 PetscCall(DMCopyAuxiliaryVec(dm, *plex));
161 }
162 }
163 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165
166 /*@C
167 SNESMonitorFields - Monitors the residual for each field separately
168
169 Collective
170
171 Input Parameters:
172 + snes - the `SNES` context, must have an attached `DM`
173 . its - iteration number
174 . fgnorm - 2-norm of residual
175 - vf - `PetscViewerAndFormat` of `PetscViewerType` `PETSCVIEWERASCII`
176
177 Level: intermediate
178
179 Note:
180 This routine prints the residual norm at each iteration.
181
182 .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
183 @*/
SNESMonitorFields(SNES snes,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat * vf)184 PetscErrorCode SNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
185 {
186 PetscViewer viewer = vf->viewer;
187 Vec res;
188 DM dm;
189 PetscSection s;
190 const PetscScalar *r;
191 PetscReal *lnorms, *norms;
192 PetscInt numFields, f, pStart, pEnd, p;
193
194 PetscFunctionBegin;
195 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 4);
196 PetscCall(SNESGetFunction(snes, &res, NULL, NULL));
197 PetscCall(SNESGetDM(snes, &dm));
198 PetscCall(DMGetLocalSection(dm, &s));
199 PetscCall(PetscSectionGetNumFields(s, &numFields));
200 PetscCall(PetscSectionGetChart(s, &pStart, &pEnd));
201 PetscCall(PetscCalloc2(numFields, &lnorms, numFields, &norms));
202 PetscCall(VecGetArrayRead(res, &r));
203 for (p = pStart; p < pEnd; ++p) {
204 for (f = 0; f < numFields; ++f) {
205 PetscInt fdof, foff, d;
206
207 PetscCall(PetscSectionGetFieldDof(s, p, f, &fdof));
208 PetscCall(PetscSectionGetFieldOffset(s, p, f, &foff));
209 for (d = 0; d < fdof; ++d) lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
210 }
211 }
212 PetscCall(VecRestoreArrayRead(res, &r));
213 PetscCallMPI(MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm)));
214 PetscCall(PetscViewerPushFormat(viewer, vf->format));
215 PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
216 PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES Function norm %14.12e [", its, (double)fgnorm));
217 for (f = 0; f < numFields; ++f) {
218 if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
219 PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)PetscSqrtReal(norms[f])));
220 }
221 PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
222 PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
223 PetscCall(PetscViewerPopFormat(viewer));
224 PetscCall(PetscFree2(lnorms, norms));
225 PetscFunctionReturn(PETSC_SUCCESS);
226 }
227
228 /********************* SNES callbacks **************************/
229
230 /*@
231 DMPlexSNESComputeObjectiveFEM - Sums the local objectives from the local input X using pointwise functions specified by the user
232
233 Input Parameters:
234 + dm - The mesh
235 . X - Local solution
236 - ctx - The application context
237
238 Output Parameter:
239 . obj - Local objective value
240
241 Level: developer
242
243 .seealso: `DM`, `DMPlexSNESComputeResidualFEM()`
244 @*/
DMPlexSNESComputeObjectiveFEM(DM dm,Vec X,PetscReal * obj,PetscCtx ctx)245 PetscErrorCode DMPlexSNESComputeObjectiveFEM(DM dm, Vec X, PetscReal *obj, PetscCtx ctx)
246 {
247 PetscInt Nf, cellHeight, cStart, cEnd;
248 PetscScalar *cintegral;
249
250 PetscFunctionBegin;
251 PetscCall(DMGetNumFields(dm, &Nf));
252 PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
253 PetscCall(DMPlexGetSimplexOrBoxCells(dm, cellHeight, &cStart, &cEnd));
254 PetscCall(PetscCalloc1((cEnd - cStart) * Nf, &cintegral));
255 PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
256 PetscCall(DMPlexComputeIntegral_Internal(dm, X, cStart, cEnd, cintegral, ctx));
257 /* Sum up values */
258 *obj = 0;
259 for (PetscInt c = cStart; c < cEnd; ++c)
260 for (PetscInt f = 0; f < Nf; ++f) *obj += PetscRealPart(cintegral[(c - cStart) * Nf + f]);
261 PetscCall(PetscLogEventBegin(DMPLEX_IntegralFEM, dm, 0, 0, 0));
262 PetscCall(PetscFree(cintegral));
263 PetscFunctionReturn(PETSC_SUCCESS);
264 }
265
266 /*@
267 DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user
268
269 Input Parameters:
270 + dm - The mesh
271 . X - Local solution
272 - ctx - The application context
273
274 Output Parameter:
275 . F - Local output vector
276
277 Level: developer
278
279 Note:
280 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
281
282 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMSNESComputeJacobianAction()`
283 @*/
DMPlexSNESComputeResidualFEM(DM dm,Vec X,Vec F,PetscCtx ctx)284 PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, PetscCtx ctx)
285 {
286 DM plex;
287 IS allcellIS;
288 PetscInt Nds, s;
289
290 PetscFunctionBegin;
291 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
292 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
293 PetscCall(DMGetNumDS(dm, &Nds));
294 for (s = 0; s < Nds; ++s) {
295 PetscDS ds;
296 IS cellIS;
297 PetscFormKey key;
298
299 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
300 key.value = 0;
301 key.field = 0;
302 key.part = 0;
303 if (!key.label) {
304 PetscCall(PetscObjectReference((PetscObject)allcellIS));
305 cellIS = allcellIS;
306 } else {
307 IS pointIS;
308
309 key.value = 1;
310 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
311 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
312 PetscCall(ISDestroy(&pointIS));
313 }
314 PetscCall(DMPlexComputeResidualByKey(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, ctx));
315 PetscCall(ISDestroy(&cellIS));
316 }
317 PetscCall(ISDestroy(&allcellIS));
318 PetscCall(DMDestroy(&plex));
319 PetscFunctionReturn(PETSC_SUCCESS);
320 }
321
322 /*@
323 DMPlexSNESComputeResidualDS - Sums the local residual into vector `F` from the local input `X` using all pointwise functions with unique keys in the `PetscDS`
324
325 Input Parameters:
326 + dm - The mesh
327 . X - Local solution
328 - ctx - The application context
329
330 Output Parameter:
331 . F - Local output vector
332
333 Level: developer
334
335 Note:
336 The residual is summed into `F`; the caller is responsible for using `VecZeroEntries()` or otherwise ensuring that any data in `F` is intentional.
337
338 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
339 @*/
DMPlexSNESComputeResidualDS(DM dm,Vec X,Vec F,PetscCtx ctx)340 PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, PetscCtx ctx)
341 {
342 DM plex;
343 IS allcellIS;
344 PetscInt Nds, s;
345
346 PetscFunctionBegin;
347 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
348 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
349 PetscCall(DMGetNumDS(dm, &Nds));
350 for (s = 0; s < Nds; ++s) {
351 PetscDS ds;
352 DMLabel label;
353 IS cellIS;
354
355 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
356 {
357 PetscWeakFormKind resmap[2] = {PETSC_WF_F0, PETSC_WF_F1};
358 PetscWeakForm wf;
359 PetscInt Nm = 2, m, Nk = 0, k, kp, off = 0;
360 PetscFormKey *reskeys;
361
362 /* Get unique residual keys */
363 for (m = 0; m < Nm; ++m) {
364 PetscInt Nkm;
365 PetscCall(PetscHMapFormGetSize(ds->wf->form[resmap[m]], &Nkm));
366 Nk += Nkm;
367 }
368 PetscCall(PetscMalloc1(Nk, &reskeys));
369 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[resmap[m]], &off, reskeys));
370 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
371 PetscCall(PetscFormKeySort(Nk, reskeys));
372 for (k = 0, kp = 1; kp < Nk; ++kp) {
373 if ((reskeys[k].label != reskeys[kp].label) || (reskeys[k].value != reskeys[kp].value)) {
374 ++k;
375 if (kp != k) reskeys[k] = reskeys[kp];
376 }
377 }
378 Nk = k;
379
380 PetscCall(PetscDSGetWeakForm(ds, &wf));
381 for (k = 0; k < Nk; ++k) {
382 DMLabel label = reskeys[k].label;
383 PetscInt val = reskeys[k].value;
384
385 if (!label) {
386 PetscCall(PetscObjectReference((PetscObject)allcellIS));
387 cellIS = allcellIS;
388 } else {
389 IS pointIS;
390
391 PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
392 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
393 PetscCall(ISDestroy(&pointIS));
394 }
395 PetscCall(DMPlexComputeResidualByKey(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, ctx));
396 PetscCall(ISDestroy(&cellIS));
397 }
398 PetscCall(PetscFree(reskeys));
399 }
400 }
401 PetscCall(ISDestroy(&allcellIS));
402 PetscCall(DMDestroy(&plex));
403 PetscFunctionReturn(PETSC_SUCCESS);
404 }
405
406 /*@
407 DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`
408
409 Input Parameters:
410 + dm - The mesh
411 - ctx - The application context
412
413 Output Parameter:
414 . X - Local solution
415
416 Level: developer
417
418 .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
419 @*/
DMPlexSNESComputeBoundaryFEM(DM dm,Vec X,PetscCtx ctx)420 PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, PetscCtx ctx)
421 {
422 DM plex;
423
424 PetscFunctionBegin;
425 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
426 PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, X, PETSC_MIN_REAL, NULL, NULL, NULL));
427 PetscCall(DMDestroy(&plex));
428 PetscFunctionReturn(PETSC_SUCCESS);
429 }
430
431 /*@
432 DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`
433
434 Input Parameters:
435 + dm - The `DM`
436 . X - Local solution vector
437 . Y - Local input vector
438 - ctx - The application context
439
440 Output Parameter:
441 . F - local output vector
442
443 Level: developer
444
445 Note:
446 Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
447
448 This only works with `DMPLEX`
449
450 Developer Note:
451 This should be called `DMPlexSNESComputeJacobianAction()`
452
453 .seealso: [](ch_snes), `DM`, `DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
454 @*/
DMSNESComputeJacobianAction(DM dm,Vec X,Vec Y,Vec F,PetscCtx ctx)455 PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, PetscCtx ctx)
456 {
457 DM plex;
458 IS allcellIS;
459 PetscInt Nds, s;
460
461 PetscFunctionBegin;
462 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
463 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
464 PetscCall(DMGetNumDS(dm, &Nds));
465 for (s = 0; s < Nds; ++s) {
466 PetscDS ds;
467 DMLabel label;
468 IS cellIS;
469
470 PetscCall(DMGetRegionNumDS(dm, s, &label, NULL, &ds, NULL));
471 {
472 PetscWeakFormKind jacmap[4] = {PETSC_WF_G0, PETSC_WF_G1, PETSC_WF_G2, PETSC_WF_G3};
473 PetscWeakForm wf;
474 PetscInt Nm = 4, m, Nk = 0, k, kp, off = 0;
475 PetscFormKey *jackeys;
476
477 /* Get unique Jacobian keys */
478 for (m = 0; m < Nm; ++m) {
479 PetscInt Nkm;
480 PetscCall(PetscHMapFormGetSize(ds->wf->form[jacmap[m]], &Nkm));
481 Nk += Nkm;
482 }
483 PetscCall(PetscMalloc1(Nk, &jackeys));
484 for (m = 0; m < Nm; ++m) PetscCall(PetscHMapFormGetKeys(ds->wf->form[jacmap[m]], &off, jackeys));
485 PetscCheck(off == Nk, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of keys %" PetscInt_FMT " should be %" PetscInt_FMT, off, Nk);
486 PetscCall(PetscFormKeySort(Nk, jackeys));
487 for (k = 0, kp = 1; kp < Nk; ++kp) {
488 if ((jackeys[k].label != jackeys[kp].label) || (jackeys[k].value != jackeys[kp].value)) {
489 ++k;
490 if (kp != k) jackeys[k] = jackeys[kp];
491 }
492 }
493 Nk = k;
494
495 PetscCall(PetscDSGetWeakForm(ds, &wf));
496 for (k = 0; k < Nk; ++k) {
497 DMLabel label = jackeys[k].label;
498 PetscInt val = jackeys[k].value;
499
500 if (!label) {
501 PetscCall(PetscObjectReference((PetscObject)allcellIS));
502 cellIS = allcellIS;
503 } else {
504 IS pointIS;
505
506 PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
507 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
508 PetscCall(ISDestroy(&pointIS));
509 }
510 PetscCall(DMPlexComputeJacobianActionByKey(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, ctx));
511 PetscCall(ISDestroy(&cellIS));
512 }
513 PetscCall(PetscFree(jackeys));
514 }
515 }
516 PetscCall(ISDestroy(&allcellIS));
517 PetscCall(DMDestroy(&plex));
518 PetscFunctionReturn(PETSC_SUCCESS);
519 }
520
521 /*@
522 DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.
523
524 Input Parameters:
525 + dm - The `DM`
526 . X - Local input vector
527 - ctx - The application context
528
529 Output Parameters:
530 + Jac - Jacobian matrix
531 - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`
532
533 Level: developer
534
535 Note:
536 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
537 like a GPU, or vectorize on a multicore machine.
538
539 .seealso: [](ch_snes), `DMPLEX`, `Mat`
540 @*/
DMPlexSNESComputeJacobianFEM(DM dm,Vec X,Mat Jac,Mat JacP,PetscCtx ctx)541 PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, PetscCtx ctx)
542 {
543 DM plex;
544 IS allcellIS;
545 PetscBool hasJac, hasPrec;
546 PetscInt Nds, s;
547
548 PetscFunctionBegin;
549 PetscCall(DMSNESConvertPlex(dm, &plex, PETSC_TRUE));
550 PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
551 PetscCall(DMGetNumDS(dm, &Nds));
552 for (s = 0; s < Nds; ++s) {
553 PetscDS ds;
554 IS cellIS;
555 PetscFormKey key;
556
557 PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
558 key.value = 0;
559 key.field = 0;
560 key.part = 0;
561 if (!key.label) {
562 PetscCall(PetscObjectReference((PetscObject)allcellIS));
563 cellIS = allcellIS;
564 } else {
565 IS pointIS;
566
567 key.value = 1;
568 PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
569 PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
570 PetscCall(ISDestroy(&pointIS));
571 }
572 if (!s) {
573 PetscCall(PetscDSHasJacobian(ds, &hasJac));
574 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
575 if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
576 PetscCall(MatZeroEntries(JacP));
577 }
578 PetscCall(DMPlexComputeJacobianByKey(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, ctx));
579 PetscCall(ISDestroy(&cellIS));
580 }
581 PetscCall(ISDestroy(&allcellIS));
582 PetscCall(DMDestroy(&plex));
583 PetscFunctionReturn(PETSC_SUCCESS);
584 }
585
586 struct _DMSNESJacobianMFCtx {
587 DM dm;
588 Vec X;
589 PetscCtx ctx;
590 };
591
DMSNESJacobianMF_Destroy_Private(Mat A)592 static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
593 {
594 struct _DMSNESJacobianMFCtx *ctx;
595
596 PetscFunctionBegin;
597 PetscCall(MatShellGetContext(A, &ctx));
598 PetscCall(MatShellSetContext(A, NULL));
599 PetscCall(DMDestroy(&ctx->dm));
600 PetscCall(VecDestroy(&ctx->X));
601 PetscCall(PetscFree(ctx));
602 PetscFunctionReturn(PETSC_SUCCESS);
603 }
604
DMSNESJacobianMF_Mult_Private(Mat A,Vec Y,Vec Z)605 static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
606 {
607 struct _DMSNESJacobianMFCtx *ctx;
608
609 PetscFunctionBegin;
610 PetscCall(MatShellGetContext(A, &ctx));
611 PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
612 PetscFunctionReturn(PETSC_SUCCESS);
613 }
614
615 /*@
616 DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
617
618 Collective
619
620 Input Parameters:
621 + dm - The `DM`
622 . X - The evaluation point for the Jacobian
623 - ctx - An application context, or `NULL`
624
625 Output Parameter:
626 . J - The `Mat`
627
628 Level: advanced
629
630 Notes:
631 Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.
632
633 This only works for `DMPLEX`
634
635 .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
636 @*/
DMSNESCreateJacobianMF(DM dm,Vec X,PetscCtx ctx,Mat * J)637 PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, PetscCtx ctx, Mat *J)
638 {
639 struct _DMSNESJacobianMFCtx *ictx;
640 PetscInt n, N;
641
642 PetscFunctionBegin;
643 PetscCall(MatCreate(PetscObjectComm((PetscObject)dm), J));
644 PetscCall(MatSetType(*J, MATSHELL));
645 PetscCall(VecGetLocalSize(X, &n));
646 PetscCall(VecGetSize(X, &N));
647 PetscCall(MatSetSizes(*J, n, n, N, N));
648 PetscCall(PetscObjectReference((PetscObject)dm));
649 PetscCall(PetscObjectReference((PetscObject)X));
650 PetscCall(PetscMalloc1(1, &ictx));
651 ictx->dm = dm;
652 ictx->X = X;
653 ictx->ctx = ctx;
654 PetscCall(MatShellSetContext(*J, ictx));
655 PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (PetscErrorCodeFn *)DMSNESJacobianMF_Destroy_Private));
656 PetscCall(MatShellSetOperation(*J, MATOP_MULT, (PetscErrorCodeFn *)DMSNESJacobianMF_Mult_Private));
657 PetscFunctionReturn(PETSC_SUCCESS);
658 }
659
MatComputeNeumannOverlap_Plex(Mat J,PetscReal t,Vec X,Vec X_t,PetscReal s,IS ovl,PetscCtx ctx)660 static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, PetscCtx ctx)
661 {
662 SNES snes;
663 Mat pJ;
664 DM ovldm, origdm;
665 DMSNES sdm;
666 PetscErrorCode (*bfun)(DM, Vec, void *);
667 PetscErrorCode (*jfun)(DM, Vec, Mat, Mat, void *);
668 void *bctx, *jctx;
669
670 PetscFunctionBegin;
671 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_MATIS", (PetscObject *)&pJ));
672 PetscCheck(pJ, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing overlapping Mat");
673 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Original_HPDDM", (PetscObject *)&origdm));
674 PetscCheck(origdm, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing original DM");
675 PetscCall(MatGetDM(pJ, &ovldm));
676 PetscCall(DMSNESGetBoundaryLocal(origdm, &bfun, &bctx));
677 PetscCall(DMSNESSetBoundaryLocal(ovldm, bfun, bctx));
678 PetscCall(DMSNESGetJacobianLocal(origdm, &jfun, &jctx));
679 PetscCall(DMSNESSetJacobianLocal(ovldm, jfun, jctx));
680 PetscCall(PetscObjectQuery((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject *)&snes));
681 if (!snes) {
682 PetscCall(SNESCreate(PetscObjectComm((PetscObject)ovl), &snes));
683 PetscCall(SNESSetDM(snes, ovldm));
684 PetscCall(PetscObjectCompose((PetscObject)ovl, "_DM_Overlap_HPDDM_SNES", (PetscObject)snes));
685 PetscCall(PetscObjectDereference((PetscObject)snes));
686 }
687 PetscCall(DMGetDMSNES(ovldm, &sdm));
688 PetscCall(VecLockReadPush(X));
689 {
690 PetscCtx ctx;
691 PetscErrorCode (*J)(SNES, Vec, Mat, Mat, void *);
692 PetscCall(DMSNESGetJacobian(ovldm, &J, &ctx));
693 PetscCallBack("SNES callback Jacobian", (*J)(snes, X, pJ, pJ, ctx));
694 }
695 PetscCall(VecLockReadPop(X));
696 /* this is a no-hop, just in case we decide to change the placeholder for the local Neumann matrix */
697 {
698 Mat locpJ;
699
700 PetscCall(MatISGetLocalMat(pJ, &locpJ));
701 PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
702 }
703 PetscFunctionReturn(PETSC_SUCCESS);
704 }
705
706 /*@
707 DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, objective, residual, and Jacobian.
708
709 Input Parameters:
710 + dm - The `DM` object
711 . use_obj - Use the objective function callback
712 - ctx - The application context that will be passed to pointwise evaluation routines
713
714 Level: developer
715
716 .seealso: [](ch_snes),`DMPLEX`, `SNES`, `PetscDSAddBoundary()`, `PetscDSSetObjective()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`
717 @*/
DMPlexSetSNESLocalFEM(DM dm,PetscBool use_obj,PetscCtx ctx)718 PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, PetscBool use_obj, PetscCtx ctx)
719 {
720 PetscBool useCeed;
721
722 PetscFunctionBegin;
723 PetscCall(DMPlexGetUseCeed(dm, &useCeed));
724 PetscCall(DMSNESSetBoundaryLocal(dm, DMPlexSNESComputeBoundaryFEM, ctx));
725 if (use_obj) PetscCall(DMSNESSetObjectiveLocal(dm, DMPlexSNESComputeObjectiveFEM, ctx));
726 if (useCeed) {
727 #ifdef PETSC_HAVE_LIBCEED
728 PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualCEED, ctx));
729 #else
730 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Cannot use CEED traversals without LibCEED. Rerun configure with --download-ceed");
731 #endif
732 } else PetscCall(DMSNESSetFunctionLocal(dm, DMPlexSNESComputeResidualFEM, ctx));
733 PetscCall(DMSNESSetJacobianLocal(dm, DMPlexSNESComputeJacobianFEM, ctx));
734 PetscCall(PetscObjectComposeFunction((PetscObject)dm, "MatComputeNeumannOverlap_C", MatComputeNeumannOverlap_Plex));
735 PetscFunctionReturn(PETSC_SUCCESS);
736 }
737
738 /*@
739 DMSNESCheckDiscretization - Check the discretization error of the exact solution
740
741 Input Parameters:
742 + snes - the `SNES` object
743 . dm - the `DM`
744 . t - the time
745 . u - a `DM` vector
746 - tol - A tolerance for the check, or -1 to print the results instead
747
748 Output Parameter:
749 . error - An array which holds the discretization error in each field, or `NULL`
750
751 Level: developer
752
753 Note:
754 The user must call `PetscDSSetExactSolution()` beforehand
755
756 Developer Note:
757 How is this related to `PetscConvEst`?
758
759 .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
760 @*/
DMSNESCheckDiscretization(SNES snes,DM dm,PetscReal t,Vec u,PetscReal tol,PetscReal error[])761 PetscErrorCode DMSNESCheckDiscretization(SNES snes, DM dm, PetscReal t, Vec u, PetscReal tol, PetscReal error[])
762 {
763 PetscErrorCode (**exacts)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, PetscCtx ctx);
764 void **ectxs;
765 PetscReal *err;
766 MPI_Comm comm;
767 PetscInt Nf, f;
768
769 PetscFunctionBegin;
770 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
771 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
772 PetscValidHeaderSpecific(u, VEC_CLASSID, 4);
773 if (error) PetscAssertPointer(error, 6);
774
775 PetscCall(DMComputeExactSolution(dm, t, u, NULL));
776 PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));
777
778 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
779 PetscCall(DMGetNumFields(dm, &Nf));
780 PetscCall(PetscCalloc3(Nf, &exacts, Nf, &ectxs, PetscMax(1, Nf), &err));
781 {
782 PetscInt Nds, s;
783
784 PetscCall(DMGetNumDS(dm, &Nds));
785 for (s = 0; s < Nds; ++s) {
786 PetscDS ds;
787 DMLabel label;
788 IS fieldIS;
789 const PetscInt *fields;
790 PetscInt dsNf, f;
791
792 PetscCall(DMGetRegionNumDS(dm, s, &label, &fieldIS, &ds, NULL));
793 PetscCall(PetscDSGetNumFields(ds, &dsNf));
794 PetscCall(ISGetIndices(fieldIS, &fields));
795 for (f = 0; f < dsNf; ++f) {
796 const PetscInt field = fields[f];
797 PetscCall(PetscDSGetExactSolution(ds, field, &exacts[field], &ectxs[field]));
798 }
799 PetscCall(ISRestoreIndices(fieldIS, &fields));
800 }
801 }
802 if (Nf > 1) {
803 PetscCall(DMComputeL2FieldDiff(dm, t, exacts, ectxs, u, err));
804 if (tol >= 0.0) {
805 for (f = 0; f < Nf; ++f) PetscCheck(err[f] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g for field %" PetscInt_FMT " exceeds tolerance %g", (double)err[f], f, (double)tol);
806 } else if (error) {
807 for (f = 0; f < Nf; ++f) error[f] = err[f];
808 } else {
809 PetscCall(PetscPrintf(comm, "L_2 Error: ["));
810 for (f = 0; f < Nf; ++f) {
811 if (f) PetscCall(PetscPrintf(comm, ", "));
812 PetscCall(PetscPrintf(comm, "%g", (double)err[f]));
813 }
814 PetscCall(PetscPrintf(comm, "]\n"));
815 }
816 } else {
817 PetscCall(DMComputeL2Diff(dm, t, exacts, ectxs, u, &err[0]));
818 if (tol >= 0.0) {
819 PetscCheck(err[0] <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Error %g exceeds tolerance %g", (double)err[0], (double)tol);
820 } else if (error) {
821 error[0] = err[0];
822 } else {
823 PetscCall(PetscPrintf(comm, "L_2 Error: %g\n", (double)err[0]));
824 }
825 }
826 PetscCall(PetscFree3(exacts, ectxs, err));
827 PetscFunctionReturn(PETSC_SUCCESS);
828 }
829
830 /*@
831 DMSNESCheckResidual - Check the residual of the exact solution
832
833 Input Parameters:
834 + snes - the `SNES` object
835 . dm - the `DM`
836 . u - a `DM` vector
837 - tol - A tolerance for the check, or -1 to print the results instead
838
839 Output Parameter:
840 . residual - The residual norm of the exact solution, or `NULL`
841
842 Level: developer
843
844 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
845 @*/
DMSNESCheckResidual(SNES snes,DM dm,Vec u,PetscReal tol,PetscReal * residual)846 PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
847 {
848 MPI_Comm comm;
849 Vec r;
850 PetscReal res;
851
852 PetscFunctionBegin;
853 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
854 PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
855 PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
856 if (residual) PetscAssertPointer(residual, 5);
857 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
858 PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
859 PetscCall(VecDuplicate(u, &r));
860 PetscCall(SNESComputeFunction(snes, u, r));
861 PetscCall(VecNorm(r, NORM_2, &res));
862 if (tol >= 0.0) {
863 PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
864 } else if (residual) {
865 *residual = res;
866 } else {
867 PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
868 PetscCall(VecFilter(r, 1.0e-10));
869 PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
870 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
871 PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", (PetscObject)snes));
872 PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
873 PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", NULL));
874 }
875 PetscCall(VecDestroy(&r));
876 PetscFunctionReturn(PETSC_SUCCESS);
877 }
878
879 /*@
880 DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
881
882 Input Parameters:
883 + snes - the `SNES` object
884 . dm - the `DM`
885 . u - a `DM` vector
886 - tol - A tolerance for the check, or -1 to print the results instead
887
888 Output Parameters:
889 + isLinear - Flag indicaing that the function looks linear, or `NULL`
890 - convRate - The rate of convergence of the linear model, or `NULL`
891
892 Level: developer
893
894 .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
895 @*/
DMSNESCheckJacobian(SNES snes,DM dm,Vec u,PetscReal tol,PetscBool * isLinear,PetscReal * convRate)896 PetscErrorCode DMSNESCheckJacobian(SNES snes, DM dm, Vec u, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
897 {
898 MPI_Comm comm;
899 PetscDS ds;
900 Mat J, M;
901 MatNullSpace nullspace;
902 PetscReal slope, intercept;
903 PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
904
905 PetscFunctionBegin;
906 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
907 if (dm) PetscValidHeaderSpecific(dm, DM_CLASSID, 2);
908 if (u) PetscValidHeaderSpecific(u, VEC_CLASSID, 3);
909 if (isLinear) PetscAssertPointer(isLinear, 5);
910 if (convRate) PetscAssertPointer(convRate, 6);
911 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
912 if (!dm) PetscCall(SNESGetDM(snes, &dm));
913 if (u) PetscCall(DMComputeExactSolution(dm, 0.0, u, NULL));
914 else PetscCall(SNESGetSolution(snes, &u));
915 /* Create and view matrices */
916 PetscCall(DMCreateMatrix(dm, &J));
917 PetscCall(DMGetDS(dm, &ds));
918 PetscCall(PetscDSHasJacobian(ds, &hasJac));
919 PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
920 if (hasJac && hasPrec) {
921 PetscCall(DMCreateMatrix(dm, &M));
922 PetscCall(SNESComputeJacobian(snes, u, J, M));
923 PetscCall(PetscObjectSetName((PetscObject)M, "Matrix used to construct preconditioner"));
924 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
925 PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
926 PetscCall(MatDestroy(&M));
927 } else {
928 PetscCall(SNESComputeJacobian(snes, u, J, J));
929 }
930 PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
931 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
932 PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
933 /* Check nullspace */
934 PetscCall(MatGetNullSpace(J, &nullspace));
935 if (nullspace) {
936 PetscBool isNull;
937 PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
938 PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
939 }
940 /* Taylor test */
941 {
942 PetscRandom rand;
943 Vec du, uhat, r, rhat, df;
944 PetscReal h;
945 PetscReal *es, *hs, *errors;
946 PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
947 PetscInt Nv, v;
948
949 /* Choose a perturbation direction */
950 PetscCall(PetscRandomCreate(comm, &rand));
951 PetscCall(VecDuplicate(u, &du));
952 PetscCall(VecSetRandom(du, rand));
953 PetscCall(PetscRandomDestroy(&rand));
954 PetscCall(VecDuplicate(u, &df));
955 PetscCall(MatMult(J, du, df));
956 /* Evaluate residual at u, F(u), save in vector r */
957 PetscCall(VecDuplicate(u, &r));
958 PetscCall(SNESComputeFunction(snes, u, r));
959 /* Look at the convergence of our Taylor approximation as we approach u */
960 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
961 PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
962 PetscCall(VecDuplicate(u, &uhat));
963 PetscCall(VecDuplicate(u, &rhat));
964 for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
965 PetscCall(VecWAXPY(uhat, h, du, u));
966 /* F(\hat u) \approx F(u) + J(u) (uhat - u) = F(u) + h * J(u) du */
967 PetscCall(SNESComputeFunction(snes, uhat, rhat));
968 PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
969 PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));
970
971 es[Nv] = errors[Nv] == 0 ? -16.0 : PetscLog10Real(errors[Nv]);
972 hs[Nv] = PetscLog10Real(h);
973 }
974 PetscCall(VecDestroy(&uhat));
975 PetscCall(VecDestroy(&rhat));
976 PetscCall(VecDestroy(&df));
977 PetscCall(VecDestroy(&r));
978 PetscCall(VecDestroy(&du));
979 for (v = 0; v < Nv; ++v) {
980 if ((tol >= 0) && (errors[v] > tol)) break;
981 else if (errors[v] > PETSC_SMALL) break;
982 }
983 if (v == Nv) isLin = PETSC_TRUE;
984 PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
985 PetscCall(PetscFree3(es, hs, errors));
986 /* Slope should be about 2 */
987 if (tol >= 0) {
988 PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
989 } else if (isLinear || convRate) {
990 if (isLinear) *isLinear = isLin;
991 if (convRate) *convRate = slope;
992 } else {
993 if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
994 else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
995 }
996 }
997 PetscCall(MatDestroy(&J));
998 PetscFunctionReturn(PETSC_SUCCESS);
999 }
1000
DMSNESCheck_Internal(SNES snes,DM dm,Vec u)1001 static PetscErrorCode DMSNESCheck_Internal(SNES snes, DM dm, Vec u)
1002 {
1003 PetscFunctionBegin;
1004 PetscCall(DMSNESCheckDiscretization(snes, dm, 0.0, u, -1.0, NULL));
1005 PetscCall(DMSNESCheckResidual(snes, dm, u, -1.0, NULL));
1006 PetscCall(DMSNESCheckJacobian(snes, dm, u, -1.0, NULL, NULL));
1007 PetscFunctionReturn(PETSC_SUCCESS);
1008 }
1009
1010 /*@
1011 DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1012
1013 Input Parameters:
1014 + snes - the `SNES` object
1015 - u - representative `SNES` vector
1016
1017 Level: developer
1018
1019 Note:
1020 The user must call `PetscDSSetExactSolution()` before this call
1021
1022 .seealso: [](ch_snes), `SNES`, `DM`
1023 @*/
DMSNESCheckFromOptions(SNES snes,Vec u)1024 PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1025 {
1026 DM dm;
1027 Vec sol;
1028 PetscBool check;
1029
1030 PetscFunctionBegin;
1031 PetscCall(PetscOptionsHasName(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-dmsnes_check", &check));
1032 if (!check) PetscFunctionReturn(PETSC_SUCCESS);
1033 PetscCall(SNESGetDM(snes, &dm));
1034 PetscCall(VecDuplicate(u, &sol));
1035 PetscCall(SNESSetSolution(snes, sol));
1036 PetscCall(DMSNESCheck_Internal(snes, dm, sol));
1037 PetscCall(VecDestroy(&sol));
1038 PetscFunctionReturn(PETSC_SUCCESS);
1039 }
1040
1041 /*@
1042 DMPlexSetSNESVariableBounds - Compute upper and lower bounds for the solution using pointsie functions from the `PetscDS`
1043
1044 Collective
1045
1046 Input Parameters:
1047 + dm - The `DM` object
1048 - snes - the `SNES` object
1049
1050 Level: intermediate
1051
1052 Notes:
1053 This calls `SNESVISetVariableBounds()` after generating the bounds vectors, so it only applied to `SNESVI` solves.
1054
1055 We project the actual bounds into the current finite element space so that they become more accurate with refinement.
1056
1057 .seealso: `SNESVISetVariableBounds()`, `SNESVI`, [](ch_snes), `DM`
1058 @*/
DMPlexSetSNESVariableBounds(DM dm,SNES snes)1059 PetscErrorCode DMPlexSetSNESVariableBounds(DM dm, SNES snes)
1060 {
1061 PetscDS ds;
1062 Vec lb, ub;
1063 PetscSimplePointFn **lfuncs, **ufuncs;
1064 void **lctxs, **uctxs;
1065 PetscBool hasBound, hasLower = PETSC_FALSE, hasUpper = PETSC_FALSE;
1066 PetscInt Nf;
1067
1068 PetscFunctionBegin;
1069 PetscCall(DMHasBound(dm, &hasBound));
1070 if (!hasBound) PetscFunctionReturn(PETSC_SUCCESS);
1071 // TODO Generalize for multiple DSes
1072 PetscCall(DMGetDS(dm, &ds));
1073 PetscCall(PetscDSGetNumFields(ds, &Nf));
1074 PetscCall(PetscMalloc4(Nf, &lfuncs, Nf, &lctxs, Nf, &ufuncs, Nf, &uctxs));
1075 for (PetscInt f = 0; f < Nf; ++f) {
1076 PetscCall(PetscDSGetLowerBound(ds, f, &lfuncs[f], &lctxs[f]));
1077 PetscCall(PetscDSGetUpperBound(ds, f, &ufuncs[f], &uctxs[f]));
1078 if (lfuncs[f]) hasLower = PETSC_TRUE;
1079 if (ufuncs[f]) hasUpper = PETSC_TRUE;
1080 }
1081 PetscCall(DMCreateGlobalVector(dm, &lb));
1082 PetscCall(DMCreateGlobalVector(dm, &ub));
1083 PetscCall(PetscObjectSetName((PetscObject)lb, "Lower Bound"));
1084 PetscCall(PetscObjectSetName((PetscObject)ub, "Upper Bound"));
1085 PetscCall(VecSet(lb, PETSC_NINFINITY));
1086 PetscCall(VecSet(ub, PETSC_INFINITY));
1087 if (hasLower) PetscCall(DMProjectFunction(dm, 0., lfuncs, lctxs, INSERT_VALUES, lb));
1088 if (hasUpper) PetscCall(DMProjectFunction(dm, 0., ufuncs, uctxs, INSERT_VALUES, ub));
1089 PetscCall(DMPlexInsertBounds(dm, PETSC_TRUE, 0., lb));
1090 PetscCall(DMPlexInsertBounds(dm, PETSC_FALSE, 0., ub));
1091 PetscCall(VecViewFromOptions(lb, NULL, "-dm_plex_snes_lb_view"));
1092 PetscCall(VecViewFromOptions(ub, NULL, "-dm_plex_snes_ub_view"));
1093 PetscCall(SNESVISetVariableBounds(snes, lb, ub));
1094 PetscCall(VecDestroy(&lb));
1095 PetscCall(VecDestroy(&ub));
1096 PetscCall(PetscFree4(lfuncs, lctxs, ufuncs, uctxs));
1097 PetscFunctionReturn(PETSC_SUCCESS);
1098 }
1099