Actual source code: dmplexsnes.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/snesimpl.h>
3: #include <petscds.h>
4: #include <petsc/private/petscimpl.h>
5: #include <petsc/private/petscfeimpl.h>
7: #ifdef PETSC_HAVE_LIBCEED
8: #include <petscdmceed.h>
9: #include <petscdmplexceed.h>
10: #endif
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: }
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.
21: Collective
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 user context
30: Output Parameter:
31: . u - The solution with a continuum pressure integral of zero
33: Level: developer
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.
38: .seealso: [](ch_snes), `SNESConvergedCorrectPressure()`
39: */
40: static PetscErrorCode SNESCorrectDiscretePressure_Private(SNES snes, PetscInt pfield, MatNullSpace nullspace, Vec u, void *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;
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: }
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.
77: Logically Collective
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 user context
87: Output Parameter:
88: . reason - `SNES_CONVERGED_ITERATING`, `SNES_CONVERGED_ITS`, or `SNES_DIVERGED_FNORM_NAN`
90: Options Database Key:
91: . -snes_convergence_test correct_pressure - see `SNESSetFromOptions()`
93: Level: advanced
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.
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.
105: .seealso: [](ch_snes), `SNES`, `DM`, `SNESConvergedDefault()`, `SNESSetConvergenceTest()`, `DMSetNullSpaceConstructor()`
106: @*/
107: PetscErrorCode SNESConvergedCorrectPressure(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gnorm, PetscReal f, SNESConvergedReason *reason, void *ctx)
108: {
109: PetscBool monitorIntegral = PETSC_FALSE;
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;
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;
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: }
141: static PetscErrorCode DMSNESConvertPlex(DM dm, DM *plex, PetscBool copy)
142: {
143: PetscBool isPlex;
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: }
166: /*@C
167: SNESMonitorFields - Monitors the residual for each field separately
169: Collective
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`
177: Level: intermediate
179: Note:
180: This routine prints the residual norm at each iteration.
182: .seealso: [](ch_snes), `SNES`, `SNESMonitorSet()`, `SNESMonitorDefault()`
183: @*/
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;
194: PetscFunctionBegin;
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;
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: }
228: /********************* SNES callbacks **************************/
230: /*@
231: DMPlexSNESComputeObjectiveFEM - Sums the local objectives from the local input X using pointwise functions specified by the user
233: Input Parameters:
234: + dm - The mesh
235: . X - Local solution
236: - user - The user context
238: Output Parameter:
239: . obj - Local objective value
241: Level: developer
243: .seealso: `DM`, `DMPlexSNESComputeResidualFEM()`
244: @*/
245: PetscErrorCode DMPlexSNESComputeObjectiveFEM(DM dm, Vec X, PetscReal *obj, void *user)
246: {
247: PetscInt Nf, cellHeight, cStart, cEnd;
248: PetscScalar *cintegral;
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, user));
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: }
266: /*@
267: DMPlexSNESComputeResidualFEM - Sums the local residual into vector `F` from the local input `X` using pointwise functions specified by the user
269: Input Parameters:
270: + dm - The mesh
271: . X - Local solution
272: - user - The user context
274: Output Parameter:
275: . F - Local output vector
277: Level: developer
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.
282: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMSNESComputeJacobianAction()`
283: @*/
284: PetscErrorCode DMPlexSNESComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
285: {
286: DM plex;
287: IS allcellIS;
288: PetscInt Nds, s;
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;
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;
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(DMPlexComputeResidual_Internal(plex, key, cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
315: PetscCall(ISDestroy(&cellIS));
316: }
317: PetscCall(ISDestroy(&allcellIS));
318: PetscCall(DMDestroy(&plex));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
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`
325: Input Parameters:
326: + dm - The mesh
327: . X - Local solution
328: - user - The user context
330: Output Parameter:
331: . F - Local output vector
333: Level: developer
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.
338: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
339: @*/
340: PetscErrorCode DMPlexSNESComputeResidualDS(DM dm, Vec X, Vec F, void *user)
341: {
342: DM plex;
343: IS allcellIS;
344: PetscInt Nds, s;
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;
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;
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;
380: PetscCall(PetscDSGetWeakForm(ds, &wf));
381: for (k = 0; k < Nk; ++k) {
382: DMLabel label = reskeys[k].label;
383: PetscInt val = reskeys[k].value;
385: if (!label) {
386: PetscCall(PetscObjectReference((PetscObject)allcellIS));
387: cellIS = allcellIS;
388: } else {
389: IS pointIS;
391: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
392: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
393: PetscCall(ISDestroy(&pointIS));
394: }
395: PetscCall(DMPlexComputeResidual_Internal(plex, reskeys[k], cellIS, PETSC_MIN_REAL, X, NULL, 0.0, F, user));
396: PetscCall(ISDestroy(&cellIS));
397: }
398: PetscCall(PetscFree(reskeys));
399: }
400: }
401: PetscCall(ISDestroy(&allcellIS));
402: PetscCall(DMDestroy(&plex));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: /*@
407: DMPlexSNESComputeBoundaryFEM - Form the boundary values for the local input `X`
409: Input Parameters:
410: + dm - The mesh
411: - user - The user context
413: Output Parameter:
414: . X - Local solution
416: Level: developer
418: .seealso: [](ch_snes), `DM`, `DMPLEX`, `DMPlexComputeJacobianAction()`
419: @*/
420: PetscErrorCode DMPlexSNESComputeBoundaryFEM(DM dm, Vec X, void *user)
421: {
422: DM plex;
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: }
431: /*@
432: DMSNESComputeJacobianAction - Compute the action of the Jacobian J(`X`) on `Y`
434: Input Parameters:
435: + dm - The `DM`
436: . X - Local solution vector
437: . Y - Local input vector
438: - user - The user context
440: Output Parameter:
441: . F - local output vector
443: Level: developer
445: Note:
446: Users will typically use `DMSNESCreateJacobianMF()` followed by `MatMult()` instead of calling this routine directly.
448: This only works with `DMPLEX`
450: Developer Note:
451: This should be called `DMPlexSNESComputeJacobianAction()`
453: .seealso: [](ch_snes), `DM`, `DMSNESCreateJacobianMF()`, `DMPlexSNESComputeResidualFEM()`
454: @*/
455: PetscErrorCode DMSNESComputeJacobianAction(DM dm, Vec X, Vec Y, Vec F, void *user)
456: {
457: DM plex;
458: IS allcellIS;
459: PetscInt Nds, s;
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;
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;
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;
495: PetscCall(PetscDSGetWeakForm(ds, &wf));
496: for (k = 0; k < Nk; ++k) {
497: DMLabel label = jackeys[k].label;
498: PetscInt val = jackeys[k].value;
500: if (!label) {
501: PetscCall(PetscObjectReference((PetscObject)allcellIS));
502: cellIS = allcellIS;
503: } else {
504: IS pointIS;
506: PetscCall(DMLabelGetStratumIS(label, val, &pointIS));
507: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
508: PetscCall(ISDestroy(&pointIS));
509: }
510: PetscCall(DMPlexComputeJacobian_Action_Internal(plex, jackeys[k], cellIS, 0.0, 0.0, X, NULL, Y, F, user));
511: PetscCall(ISDestroy(&cellIS));
512: }
513: PetscCall(PetscFree(jackeys));
514: }
515: }
516: PetscCall(ISDestroy(&allcellIS));
517: PetscCall(DMDestroy(&plex));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*@
522: DMPlexSNESComputeJacobianFEM - Form the local portion of the Jacobian matrix `Jac` at the local solution `X` using pointwise functions specified by the user.
524: Input Parameters:
525: + dm - The `DM`
526: . X - Local input vector
527: - user - The user context
529: Output Parameters:
530: + Jac - Jacobian matrix
531: - JacP - approximate Jacobian from which the preconditioner will be built, often `Jac`
533: Level: developer
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.
539: .seealso: [](ch_snes), `DMPLEX`, `Mat`
540: @*/
541: PetscErrorCode DMPlexSNESComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, void *user)
542: {
543: DM plex;
544: IS allcellIS;
545: PetscBool hasJac, hasPrec;
546: PetscInt Nds, s;
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;
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;
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(DMPlexComputeJacobian_Internal(plex, key, cellIS, 0.0, 0.0, X, NULL, Jac, JacP, user));
579: PetscCall(ISDestroy(&cellIS));
580: }
581: PetscCall(ISDestroy(&allcellIS));
582: PetscCall(DMDestroy(&plex));
583: PetscFunctionReturn(PETSC_SUCCESS);
584: }
586: struct _DMSNESJacobianMFCtx {
587: DM dm;
588: Vec X;
589: void *ctx;
590: };
592: static PetscErrorCode DMSNESJacobianMF_Destroy_Private(Mat A)
593: {
594: struct _DMSNESJacobianMFCtx *ctx;
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: }
605: static PetscErrorCode DMSNESJacobianMF_Mult_Private(Mat A, Vec Y, Vec Z)
606: {
607: struct _DMSNESJacobianMFCtx *ctx;
609: PetscFunctionBegin;
610: PetscCall(MatShellGetContext(A, &ctx));
611: PetscCall(DMSNESComputeJacobianAction(ctx->dm, ctx->X, Y, Z, ctx->ctx));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@
616: DMSNESCreateJacobianMF - Create a `Mat` which computes the action of the Jacobian matrix-free
618: Collective
620: Input Parameters:
621: + dm - The `DM`
622: . X - The evaluation point for the Jacobian
623: - user - A user context, or `NULL`
625: Output Parameter:
626: . J - The `Mat`
628: Level: advanced
630: Notes:
631: Vec `X` is kept in `J`, so updating `X` then updates the evaluation point.
633: This only works for `DMPLEX`
635: .seealso: [](ch_snes), `DM`, `SNES`, `DMSNESComputeJacobianAction()`
636: @*/
637: PetscErrorCode DMSNESCreateJacobianMF(DM dm, Vec X, void *user, Mat *J)
638: {
639: struct _DMSNESJacobianMFCtx *ctx;
640: PetscInt n, N;
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, &ctx));
651: ctx->dm = dm;
652: ctx->X = X;
653: ctx->ctx = user;
654: PetscCall(MatShellSetContext(*J, ctx));
655: PetscCall(MatShellSetOperation(*J, MATOP_DESTROY, (void (*)(void))DMSNESJacobianMF_Destroy_Private));
656: PetscCall(MatShellSetOperation(*J, MATOP_MULT, (void (*)(void))DMSNESJacobianMF_Mult_Private));
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: static PetscErrorCode MatComputeNeumannOverlap_Plex(Mat J, PetscReal t, Vec X, Vec X_t, PetscReal s, IS ovl, void *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;
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: void *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;
700: PetscCall(MatISGetLocalMat(pJ, &locpJ));
701: PetscCall(MatCopy(locpJ, J, SAME_NONZERO_PATTERN));
702: }
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: /*@
707: DMPlexSetSNESLocalFEM - Use `DMPLEX`'s internal FEM routines to compute `SNES` boundary values, objective, residual, and Jacobian.
709: Input Parameters:
710: + dm - The `DM` object
711: . use_obj - Use the objective function callback
712: - ctx - The user context that will be passed to pointwise evaluation routines
714: Level: developer
716: .seealso: [](ch_snes),`DMPLEX`, `SNES`, `PetscDSAddBoundary()`, `PetscDSSetObjective()`, `PetscDSSetResidual()`, `PetscDSSetJacobian()`
717: @*/
718: PetscErrorCode DMPlexSetSNESLocalFEM(DM dm, PetscBool use_obj, void *ctx)
719: {
720: PetscBool useCeed;
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: }
738: /*@
739: DMSNESCheckDiscretization - Check the discretization error of the exact solution
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
748: Output Parameter:
749: . error - An array which holds the discretization error in each field, or `NULL`
751: Level: developer
753: Note:
754: The user must call `PetscDSSetExactSolution()` beforehand
756: Developer Note:
757: How is this related to `PetscConvEst`?
759: .seealso: [](ch_snes), `PetscDSSetExactSolution()`, `DNSNESCheckFromOptions()`, `DMSNESCheckResidual()`, `DMSNESCheckJacobian()`
760: @*/
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, void *ctx);
764: void **ectxs;
765: PetscReal *err;
766: MPI_Comm comm;
767: PetscInt Nf, f;
769: PetscFunctionBegin;
773: if (error) PetscAssertPointer(error, 6);
775: PetscCall(DMComputeExactSolution(dm, t, u, NULL));
776: PetscCall(VecViewFromOptions(u, NULL, "-vec_view"));
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;
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;
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: }
830: /*@
831: DMSNESCheckResidual - Check the residual of the exact solution
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
839: Output Parameter:
840: . residual - The residual norm of the exact solution, or `NULL`
842: Level: developer
844: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
845: @*/
846: PetscErrorCode DMSNESCheckResidual(SNES snes, DM dm, Vec u, PetscReal tol, PetscReal *residual)
847: {
848: MPI_Comm comm;
849: Vec r;
850: PetscReal res;
852: PetscFunctionBegin;
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: }
879: /*@
880: DMSNESCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
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
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`
892: Level: developer
894: .seealso: [](ch_snes), `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
895: @*/
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;
905: PetscFunctionBegin;
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, "Preconditioning Matrix"));
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;
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]));
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: }
1001: 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: }
1010: /*@
1011: DMSNESCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
1013: Input Parameters:
1014: + snes - the `SNES` object
1015: - u - representative `SNES` vector
1017: Level: developer
1019: Note:
1020: The user must call `PetscDSSetExactSolution()` before this call
1022: .seealso: [](ch_snes), `SNES`, `DM`
1023: @*/
1024: PetscErrorCode DMSNESCheckFromOptions(SNES snes, Vec u)
1025: {
1026: DM dm;
1027: Vec sol;
1028: PetscBool check;
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: }
1041: static PetscErrorCode lower_inf(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1042: {
1043: PetscFunctionBegin;
1044: for (PetscInt c = 0; c < Nc; ++c) u[c] = PETSC_NINFINITY;
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: static PetscErrorCode upper_inf(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
1049: {
1050: PetscFunctionBegin;
1051: for (PetscInt c = 0; c < Nc; ++c) u[c] = PETSC_INFINITY;
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }
1055: PetscErrorCode DMPlexSetSNESVariableBounds(DM dm, SNES snes)
1056: {
1057: PetscDS ds;
1058: Vec lb, ub;
1059: PetscSimplePointFunc *lfuncs, *ufuncs;
1060: void **lctxs, **uctxs;
1061: PetscBool hasBound = PETSC_FALSE;
1062: PetscInt Nf;
1064: PetscFunctionBegin;
1065: // TODO Generalize for multiple DSes
1066: PetscCall(DMGetDS(dm, &ds));
1067: PetscCall(PetscDSGetNumFields(ds, &Nf));
1068: PetscCall(PetscMalloc4(Nf, &lfuncs, Nf, &lctxs, Nf, &ufuncs, Nf, &uctxs));
1069: for (PetscInt f = 0; f < Nf; ++f) {
1070: PetscCall(PetscDSGetLowerBound(ds, f, &lfuncs[f], &lctxs[f]));
1071: PetscCall(PetscDSGetUpperBound(ds, f, &ufuncs[f], &uctxs[f]));
1072: if (lfuncs[f] || ufuncs[f]) hasBound = PETSC_TRUE;
1073: if (!lfuncs[f]) lfuncs[f] = lower_inf;
1074: if (!ufuncs[f]) ufuncs[f] = upper_inf;
1075: }
1076: if (!hasBound) goto end;
1077: PetscCall(DMCreateGlobalVector(dm, &lb));
1078: PetscCall(DMCreateGlobalVector(dm, &ub));
1079: PetscCall(DMProjectFunction(dm, 0., lfuncs, lctxs, INSERT_VALUES, lb));
1080: PetscCall(DMProjectFunction(dm, 0., ufuncs, uctxs, INSERT_VALUES, ub));
1081: PetscCall(VecViewFromOptions(lb, NULL, "-dm_plex_snes_lb_view"));
1082: PetscCall(VecViewFromOptions(ub, NULL, "-dm_plex_snes_ub_view"));
1083: PetscCall(SNESVISetVariableBounds(snes, lb, ub));
1084: PetscCall(VecDestroy(&lb));
1085: PetscCall(VecDestroy(&ub));
1086: end:
1087: PetscCall(PetscFree4(lfuncs, lctxs, ufuncs, uctxs));
1088: PetscFunctionReturn(PETSC_SUCCESS);
1089: }