Actual source code: vi.c
1: #include <petsc/private/snesimpl.h>
2: #include <petscdm.h>
4: /*@C
5: SNESVISetComputeVariableBounds - Sets a function that is called to compute the bounds on variable for
6: (differential) variable inequalities.
8: Input Parameters:
9: + snes - the `SNES` context
10: - compute - function that computes the bounds
12: Calling sequence of `compute`:
13: + snes - the `SNES` context
14: . lower - vector to hold lower bounds
15: - higher - vector to hold upper bounds
17: Level: advanced
19: Notes:
20: Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS`, and semi-smooth `SNESVINEWTONSSLS` solvers.
22: For entries with no bounds you can set `PETSC_NINFINITY` or `PETSC_INFINITY`
24: You may use `SNESVISetVariableBounds()` to provide the bounds once if they will never change
26: If you have associated a `DM` with the `SNES` and provided a function to the `DM` via `DMSetVariableBounds()` that will be used automatically
27: to provide the bounds and you need not use this function.
29: See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets
31: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `DMSetVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`,
32: `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
33: @*/
34: PetscErrorCode SNESVISetComputeVariableBounds(SNES snes, PetscErrorCode (*compute)(SNES snes, Vec lower, Vec higher))
35: {
36: PetscErrorCode (*f)(SNES, PetscErrorCode (*)(SNES, Vec, Vec));
38: PetscFunctionBegin;
40: PetscCall(PetscObjectQueryFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", &f));
41: if (f) PetscUseMethod(snes, "SNESVISetComputeVariableBounds_C", (SNES, PetscErrorCode (*)(SNES, Vec, Vec)), (snes, compute));
42: else PetscCall(SNESVISetComputeVariableBounds_VI(snes, compute));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: PetscErrorCode SNESVISetComputeVariableBounds_VI(SNES snes, SNESVIComputeVariableBoundsFn *compute)
47: {
48: PetscFunctionBegin;
49: snes->ops->computevariablebounds = compute;
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode SNESVIMonitorResidual(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
54: {
55: Vec X, F, Finactive;
56: IS isactive;
58: PetscFunctionBegin;
60: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
61: PetscCall(SNESGetSolution(snes, &X));
62: PetscCall(SNESVIGetActiveSetIS(snes, X, F, &isactive));
63: PetscCall(VecDuplicate(F, &Finactive));
64: PetscCall(PetscObjectCompose((PetscObject)Finactive, "__Vec_bc_zero__", (PetscObject)snes));
65: PetscCall(VecCopy(F, Finactive));
66: PetscCall(VecISSet(Finactive, isactive, 0.0));
67: PetscCall(ISDestroy(&isactive));
68: PetscCall(PetscViewerPushFormat(vf->viewer, vf->format));
69: PetscCall(VecView(Finactive, vf->viewer));
70: PetscCall(PetscViewerPopFormat(vf->viewer));
71: PetscCall(PetscObjectCompose((PetscObject)Finactive, "__Vec_bc_zero__", NULL));
72: PetscCall(VecDestroy(&Finactive));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode SNESVIMonitorActive(SNES snes, PetscInt its, PetscReal fgnorm, PetscViewerAndFormat *vf)
77: {
78: Vec X, F, A;
79: IS isactive;
81: PetscFunctionBegin;
83: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
84: PetscCall(SNESGetSolution(snes, &X));
85: PetscCall(SNESVIGetActiveSetIS(snes, X, F, &isactive));
86: PetscCall(VecDuplicate(F, &A));
87: PetscCall(PetscObjectCompose((PetscObject)A, "__Vec_bc_zero__", (PetscObject)snes));
88: PetscCall(VecSet(A, 0.));
89: PetscCall(VecISSet(A, isactive, 1.));
90: PetscCall(ISDestroy(&isactive));
91: PetscCall(PetscViewerPushFormat(vf->viewer, vf->format));
92: PetscCall(VecView(A, vf->viewer));
93: PetscCall(PetscViewerPopFormat(vf->viewer));
94: PetscCall(PetscObjectCompose((PetscObject)A, "__Vec_bc_zero__", NULL));
95: PetscCall(VecDestroy(&A));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode SNESMonitorVI(SNES snes, PetscInt its, PetscReal fgnorm, void *dummy)
100: {
101: PetscViewer viewer = (PetscViewer)dummy;
102: const PetscScalar *x, *xl, *xu, *f;
103: PetscInt i, n, act[2] = {0, 0}, fact[2], N;
104: /* Number of components that actually hit the bounds (c.f. active variables) */
105: PetscInt act_bound[2] = {0, 0}, fact_bound[2];
106: PetscReal rnorm, fnorm, zerotolerance = snes->vizerotolerance;
107: double tmp;
109: PetscFunctionBegin;
111: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
112: PetscCall(VecGetSize(snes->vec_sol, &N));
113: PetscCall(VecGetArrayRead(snes->xl, &xl));
114: PetscCall(VecGetArrayRead(snes->xu, &xu));
115: PetscCall(VecGetArrayRead(snes->vec_sol, &x));
116: PetscCall(VecGetArrayRead(snes->vec_func, &f));
118: rnorm = 0.0;
119: for (i = 0; i < n; i++) {
120: if ((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0)) rnorm += PetscRealPart(PetscConj(f[i]) * f[i]);
121: else if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance && PetscRealPart(f[i]) > 0.0) act[0]++;
122: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance && PetscRealPart(f[i]) < 0.0) act[1]++;
123: else SETERRQ(PetscObjectComm((PetscObject)snes), PETSC_ERR_PLIB, "Can never get here");
124: }
126: for (i = 0; i < n; i++) {
127: if (PetscRealPart(x[i]) <= PetscRealPart(xl[i]) + zerotolerance) act_bound[0]++;
128: else if (PetscRealPart(x[i]) >= PetscRealPart(xu[i]) - zerotolerance) act_bound[1]++;
129: }
130: PetscCall(VecRestoreArrayRead(snes->vec_func, &f));
131: PetscCall(VecRestoreArrayRead(snes->xl, &xl));
132: PetscCall(VecRestoreArrayRead(snes->xu, &xu));
133: PetscCall(VecRestoreArrayRead(snes->vec_sol, &x));
134: PetscCallMPI(MPIU_Allreduce(&rnorm, &fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
135: PetscCallMPI(MPIU_Allreduce(act, fact, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
136: PetscCallMPI(MPIU_Allreduce(act_bound, fact_bound, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
137: fnorm = PetscSqrtReal(fnorm);
139: PetscCall(PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel));
140: if (snes->ntruebounds) tmp = ((double)(fact[0] + fact[1])) / ((double)snes->ntruebounds);
141: else tmp = 0.0;
142: PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " SNES VI Function norm %g Active lower constraints %" PetscInt_FMT "/%" PetscInt_FMT " upper constraints %" PetscInt_FMT "/%" PetscInt_FMT " Percent of total %g Percent of bounded %g\n", its, (double)fnorm, fact[0], fact_bound[0], fact[1], fact_bound[1], ((double)(fact[0] + fact[1])) / ((double)N), tmp));
144: PetscCall(PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*
149: Checks if J^T F = 0 which implies we've found a local minimum of the norm of the function,
150: || F(u) ||_2 but not a zero, F(u) = 0. In the case when one cannot compute J^T F we use the fact that
151: 0 = (J^T F)^T W = F^T J W iff W not in the null space of J. Thanks for Jorge More
152: for this trick. One assumes that the probability that W is in the null space of J is very, very small.
153: */
154: PetscErrorCode SNESVICheckLocalMin_Private(SNES snes, Mat A, Vec F, Vec W, PetscReal fnorm, PetscBool *ismin)
155: {
156: PetscReal a1;
157: PetscBool hastranspose;
159: PetscFunctionBegin;
160: *ismin = PETSC_FALSE;
161: PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, &hastranspose));
162: if (hastranspose) {
163: /* Compute || J^T F|| */
164: PetscCall(MatMultTranspose(A, F, W));
165: PetscCall(VecNorm(W, NORM_2, &a1));
166: PetscCall(PetscInfo(snes, "|| J^T F|| %g near zero implies found a local minimum\n", (double)(a1 / fnorm)));
167: if (a1 / fnorm < 1.e-4) *ismin = PETSC_TRUE;
168: } else {
169: Vec work;
170: PetscScalar result;
171: PetscReal wnorm;
173: PetscCall(VecSetRandom(W, NULL));
174: PetscCall(VecNorm(W, NORM_2, &wnorm));
175: PetscCall(VecDuplicate(W, &work));
176: PetscCall(MatMult(A, W, work));
177: PetscCall(VecDot(F, work, &result));
178: PetscCall(VecDestroy(&work));
179: a1 = PetscAbsScalar(result) / (fnorm * wnorm);
180: PetscCall(PetscInfo(snes, "(F^T J random)/(|| F ||*||random|| %g near zero implies found a local minimum\n", (double)a1));
181: if (a1 < 1.e-4) *ismin = PETSC_TRUE;
182: }
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*
187: SNESConvergedDefault_VI - Checks the convergence of the semismooth newton algorithm.
189: Notes:
190: The convergence criterion currently implemented is
191: merit < abstol
192: merit < rtol*merit_initial
193: */
194: PetscErrorCode SNESConvergedDefault_VI(SNES snes, PetscInt it, PetscReal xnorm, PetscReal gradnorm, PetscReal fnorm, SNESConvergedReason *reason, void *dummy)
195: {
196: PetscFunctionBegin;
198: PetscAssertPointer(reason, 6);
200: *reason = SNES_CONVERGED_ITERATING;
202: if (!it) {
203: /* set parameter for default relative tolerance convergence test */
204: snes->ttol = fnorm * snes->rtol;
205: }
206: if (fnorm != fnorm) {
207: PetscCall(PetscInfo(snes, "Failed to converged, function norm is NaN\n"));
208: *reason = SNES_DIVERGED_FUNCTION_NANORINF;
209: } else if (fnorm < snes->abstol && (it || !snes->forceiteration)) {
210: PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g\n", (double)fnorm, (double)snes->abstol));
211: *reason = SNES_CONVERGED_FNORM_ABS;
212: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
213: PetscCall(PetscInfo(snes, "Exceeded maximum number of function evaluations: %" PetscInt_FMT " > %" PetscInt_FMT "\n", snes->nfuncs, snes->max_funcs));
214: *reason = SNES_DIVERGED_FUNCTION_COUNT;
215: }
217: if (it && !*reason) {
218: if (fnorm < snes->ttol) {
219: PetscCall(PetscInfo(snes, "Converged due to function norm %g < %g (relative tolerance)\n", (double)fnorm, (double)snes->ttol));
220: *reason = SNES_CONVERGED_FNORM_RELATIVE;
221: }
222: }
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*
227: SNESVIProjectOntoBounds - Projects X onto the feasible region so that Xl[i] <= X[i] <= Xu[i] for i = 1...n.
229: Input Parameters:
230: . SNES - nonlinear solver context
232: Output Parameters:
233: . X - Bound projected X
235: */
237: PetscErrorCode SNESVIProjectOntoBounds(SNES snes, Vec X)
238: {
239: const PetscScalar *xl, *xu;
240: PetscScalar *x;
241: PetscInt i, n;
243: PetscFunctionBegin;
244: PetscCall(VecGetLocalSize(X, &n));
245: PetscCall(VecGetArray(X, &x));
246: PetscCall(VecGetArrayRead(snes->xl, &xl));
247: PetscCall(VecGetArrayRead(snes->xu, &xu));
249: for (i = 0; i < n; i++) {
250: if (PetscRealPart(x[i]) < PetscRealPart(xl[i])) x[i] = xl[i];
251: else if (PetscRealPart(x[i]) > PetscRealPart(xu[i])) x[i] = xu[i];
252: }
253: PetscCall(VecRestoreArray(X, &x));
254: PetscCall(VecRestoreArrayRead(snes->xl, &xl));
255: PetscCall(VecRestoreArrayRead(snes->xu, &xu));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*@
260: SNESVIGetActiveSetIS - Gets the global indices for the active set variables
262: Input Parameters:
263: + snes - the `SNES` context
264: . X - the `snes` solution vector
265: - F - the nonlinear function vector
267: Output Parameter:
268: . ISact - active set index set
270: Level: developer
272: Note:
273: See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets
275: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`
276: @*/
277: PetscErrorCode SNESVIGetActiveSetIS(SNES snes, Vec X, Vec F, IS *ISact)
278: {
279: Vec Xl = snes->xl, Xu = snes->xu;
280: const PetscScalar *x, *f, *xl, *xu;
281: PetscInt *idx_act, i, nlocal, nloc_isact = 0, ilow, ihigh, i1 = 0;
282: PetscReal zerotolerance = snes->vizerotolerance;
284: PetscFunctionBegin;
285: PetscCall(VecGetLocalSize(X, &nlocal));
286: PetscCall(VecGetOwnershipRange(X, &ilow, &ihigh));
287: PetscCall(VecGetArrayRead(X, &x));
288: PetscCall(VecGetArrayRead(Xl, &xl));
289: PetscCall(VecGetArrayRead(Xu, &xu));
290: PetscCall(VecGetArrayRead(F, &f));
291: /* Compute active set size */
292: for (i = 0; i < nlocal; i++) {
293: if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) nloc_isact++;
294: }
296: PetscCall(PetscMalloc1(nloc_isact, &idx_act));
298: /* Set active set indices */
299: for (i = 0; i < nlocal; i++) {
300: if (!((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0))) idx_act[i1++] = ilow + i;
301: }
303: /* Create active set IS */
304: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), nloc_isact, idx_act, PETSC_OWN_POINTER, ISact));
306: PetscCall(VecRestoreArrayRead(X, &x));
307: PetscCall(VecRestoreArrayRead(Xl, &xl));
308: PetscCall(VecRestoreArrayRead(Xu, &xu));
309: PetscCall(VecRestoreArrayRead(F, &f));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@
314: SNESVIComputeInactiveSetFnorm - Computes the function norm for variational inequalities on the inactive set
316: Input Parameters:
317: + snes - the `SNES` context
318: . F - the nonlinear function vector
319: - X - the `SNES` solution vector
321: Output Parameter:
322: . fnorm - the function norm
324: Level: developer
326: Note:
327: See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets
329: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESLineSearchSetVIFunctions()`
330: @*/
331: PetscErrorCode SNESVIComputeInactiveSetFnorm(SNES snes, Vec F, Vec X, PetscReal *fnorm)
332: {
333: const PetscScalar *x, *xl, *xu, *f;
334: PetscInt i, n;
335: PetscReal zerotolerance = snes->vizerotolerance;
337: PetscFunctionBegin;
339: PetscAssertPointer(fnorm, 4);
340: PetscCall(VecGetLocalSize(X, &n));
341: PetscCall(VecGetArrayRead(snes->xl, &xl));
342: PetscCall(VecGetArrayRead(snes->xu, &xu));
343: PetscCall(VecGetArrayRead(X, &x));
344: PetscCall(VecGetArrayRead(F, &f));
345: *fnorm = 0.0;
346: for (i = 0; i < n; i++) {
347: if ((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0)) *fnorm += PetscRealPart(PetscConj(f[i]) * f[i]);
348: }
349: PetscCall(VecRestoreArrayRead(F, &f));
350: PetscCall(VecRestoreArrayRead(snes->xl, &xl));
351: PetscCall(VecRestoreArrayRead(snes->xu, &xu));
352: PetscCall(VecRestoreArrayRead(X, &x));
353: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fnorm, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
354: *fnorm = PetscSqrtReal(*fnorm);
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: /*@
359: SNESVIComputeInactiveSetFtY - Computes the directional derivative for variational inequalities on the inactive set,
360: assuming that there exists some $G(x)$ for which the `SNESFunctionFn` $F(x) = grad G(x)$ (relevant for some line search algorithms)
362: Input Parameters:
363: + snes - the `SNES` context
364: . F - the nonlinear function vector
365: . X - the `SNES` solution vector
366: - Y - the direction vector
368: Output Parameter:
369: . fty - the directional derivative
371: Level: developer
373: Note:
374: See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets
376: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`
377: @*/
378: PetscErrorCode SNESVIComputeInactiveSetFtY(SNES snes, Vec F, Vec X, Vec Y, PetscScalar *fty)
379: {
380: const PetscScalar *x, *xl, *xu, *y, *f;
381: PetscInt i, n;
382: PetscReal zerotolerance = snes->vizerotolerance;
384: PetscFunctionBegin;
386: PetscAssertPointer(fty, 5);
387: PetscCall(VecGetLocalSize(X, &n));
388: PetscCall(VecGetArrayRead(F, &f));
389: PetscCall(VecGetArrayRead(X, &x));
390: PetscCall(VecGetArrayRead(snes->xl, &xl));
391: PetscCall(VecGetArrayRead(snes->xu, &xu));
392: PetscCall(VecGetArrayRead(Y, &y));
393: *fty = 0.0;
394: for (i = 0; i < n; i++) {
395: if ((PetscRealPart(x[i]) > PetscRealPart(xl[i]) + zerotolerance || (PetscRealPart(f[i]) <= 0.0)) && ((PetscRealPart(x[i]) < PetscRealPart(xu[i]) - zerotolerance) || PetscRealPart(f[i]) >= 0.0)) *fty += f[i] * PetscConj(y[i]);
396: }
397: PetscCall(VecRestoreArrayRead(F, &f));
398: PetscCall(VecRestoreArrayRead(X, &x));
399: PetscCall(VecRestoreArrayRead(snes->xl, &xl));
400: PetscCall(VecRestoreArrayRead(snes->xu, &xu));
401: PetscCall(VecRestoreArrayRead(Y, &y));
402: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fty, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)snes)));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: static PetscErrorCode SNESVIDMComputeVariableBounds(SNES snes, Vec xl, Vec xu)
407: {
408: PetscFunctionBegin;
409: PetscCall(DMComputeVariableBounds(snes->dm, xl, xu));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: // PetscClangLinter pragma disable: -fdoc-sowing-chars
414: /*
415: SNESSetUp_VI - Does setup common to all VI solvers -- basically makes sure bounds have been properly set up
416: of the SNESVI nonlinear solver.
418: Input Parameter:
419: . snes - the SNES context
421: Level: developer
423: Note:
424: For basic use of the SNES solvers, the user need not explicitly call
425: SNESSetUp(), since these actions will automatically occur during
426: the call to SNESSolve().
428: .seealso: `SNESSetUp()`
429: */
430: PetscErrorCode SNESSetUp_VI(SNES snes)
431: {
432: PetscInt i_start[3], i_end[3];
434: PetscFunctionBegin;
435: PetscCall(SNESSetWorkVecs(snes, 1));
436: PetscCall(SNESSetUpMatrices(snes));
438: if (!snes->ops->computevariablebounds && snes->dm) {
439: PetscBool flag;
440: PetscCall(DMHasVariableBounds(snes->dm, &flag));
441: if (flag) snes->ops->computevariablebounds = SNESVIDMComputeVariableBounds;
442: }
443: if (!snes->usersetbounds) {
444: if (snes->ops->computevariablebounds) {
445: if (!snes->xl) PetscCall(VecDuplicate(snes->work[0], &snes->xl));
446: if (!snes->xu) PetscCall(VecDuplicate(snes->work[0], &snes->xu));
447: PetscUseTypeMethod(snes, computevariablebounds, snes->xl, snes->xu);
448: } else if (!snes->xl && !snes->xu) {
449: /* If the lower and upper bound on variables are not set, set it to -Inf and Inf */
450: PetscCall(VecDuplicate(snes->work[0], &snes->xl));
451: PetscCall(VecSet(snes->xl, PETSC_NINFINITY));
452: PetscCall(VecDuplicate(snes->work[0], &snes->xu));
453: PetscCall(VecSet(snes->xu, PETSC_INFINITY));
454: } else {
455: /* Check if lower bound, upper bound and solution vector distribution across the processors is identical */
456: PetscCall(VecGetOwnershipRange(snes->work[0], i_start, i_end));
457: PetscCall(VecGetOwnershipRange(snes->xl, i_start + 1, i_end + 1));
458: PetscCall(VecGetOwnershipRange(snes->xu, i_start + 2, i_end + 2));
459: if ((i_start[0] != i_start[1]) || (i_start[0] != i_start[2]) || (i_end[0] != i_end[1]) || (i_end[0] != i_end[2]))
460: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Distribution of lower bound, upper bound and the solution vector should be identical across all the processors.");
461: }
462: }
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
465: PetscErrorCode SNESReset_VI(SNES snes)
466: {
467: PetscFunctionBegin;
468: PetscCall(VecDestroy(&snes->xl));
469: PetscCall(VecDestroy(&snes->xu));
470: snes->usersetbounds = PETSC_FALSE;
471: PetscFunctionReturn(PETSC_SUCCESS);
472: }
474: /*
475: SNESDestroy_VI - Destroys the private SNES_VI context that was created
476: with SNESCreate_VI().
478: Input Parameter:
479: . snes - the SNES context
481: Application Interface Routine: SNESDestroy()
482: */
483: PetscErrorCode SNESDestroy_VI(SNES snes)
484: {
485: PetscFunctionBegin;
486: PetscCall(PetscFree(snes->data));
488: /* clear composed functions */
489: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", NULL));
490: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", NULL));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*@
495: SNESVISetVariableBounds - Sets the lower and upper bounds for the solution vector. `xl` <= x <= `xu`. This allows solving
496: (differential) variable inequalities.
498: Input Parameters:
499: + snes - the `SNES` context.
500: . xl - lower bound.
501: - xu - upper bound.
503: Level: advanced
505: Notes:
506: If this routine is not called then the lower and upper bounds are set to
507: `PETSC_NINFINITY` and `PETSC_INFINITY` respectively during `SNESSetUp()`.
509: Problems with bound constraints can be solved with the reduced space, `SNESVINEWTONRSLS` or semi-smooth `SNESVINEWTONSSLS` solvers.
511: For particular components that have no bounds you can use `PETSC_NINFINITY` or `PETSC_INFINITY`
513: `SNESVISetComputeVariableBounds()` can be used to provide a function that computes the bounds. This should be used if you are using, for example, grid
514: sequencing and need bounds set for a variety of vectors
516: See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets
518: .seealso: [](sec_vi), `SNES`, `SNESVIGetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
519: @*/
520: PetscErrorCode SNESVISetVariableBounds(SNES snes, Vec xl, Vec xu)
521: {
522: PetscErrorCode (*f)(SNES, Vec, Vec);
524: PetscFunctionBegin;
528: PetscCall(PetscObjectQueryFunction((PetscObject)snes, "SNESVISetVariableBounds_C", &f));
529: if (f) PetscUseMethod(snes, "SNESVISetVariableBounds_C", (SNES, Vec, Vec), (snes, xl, xu));
530: else PetscCall(SNESVISetVariableBounds_VI(snes, xl, xu));
531: snes->usersetbounds = PETSC_TRUE;
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: PetscErrorCode SNESVISetVariableBounds_VI(SNES snes, Vec xl, Vec xu)
536: {
537: const PetscScalar *xxl, *xxu;
538: PetscInt i, n, cnt = 0;
540: PetscFunctionBegin;
541: PetscCall(SNESGetFunction(snes, &snes->vec_func, NULL, NULL));
542: PetscCheck(snes->vec_func, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() first");
543: {
544: PetscInt xlN, xuN, N;
545: PetscCall(VecGetSize(xl, &xlN));
546: PetscCall(VecGetSize(xu, &xuN));
547: PetscCall(VecGetSize(snes->vec_func, &N));
548: PetscCheck(xlN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths lower bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xlN, N);
549: PetscCheck(xuN == N, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Incompatible vector lengths: upper bound = %" PetscInt_FMT " solution vector = %" PetscInt_FMT, xuN, N);
550: }
551: PetscCall(PetscObjectReference((PetscObject)xl));
552: PetscCall(PetscObjectReference((PetscObject)xu));
553: PetscCall(VecDestroy(&snes->xl));
554: PetscCall(VecDestroy(&snes->xu));
555: snes->xl = xl;
556: snes->xu = xu;
557: PetscCall(VecGetLocalSize(xl, &n));
558: PetscCall(VecGetArrayRead(xl, &xxl));
559: PetscCall(VecGetArrayRead(xu, &xxu));
560: for (i = 0; i < n; i++) cnt += ((xxl[i] != PETSC_NINFINITY) || (xxu[i] != PETSC_INFINITY));
562: PetscCallMPI(MPIU_Allreduce(&cnt, &snes->ntruebounds, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)snes)));
563: PetscCall(VecRestoreArrayRead(xl, &xxl));
564: PetscCall(VecRestoreArrayRead(xu, &xxu));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*@
569: SNESVIGetVariableBounds - Gets the lower and upper bounds for the solution vector. `xl` <= x <= `xu`. These are used in solving
570: (differential) variable inequalities.
572: Input Parameters:
573: + snes - the `SNES` context.
574: . xl - lower bound (may be `NULL`)
575: - xu - upper bound (may be `NULL`)
577: Level: advanced
579: Note:
580: These vectors are owned by the `SNESVI` and should not be destroyed by the caller
582: .seealso: [](sec_vi), `SNES`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESSetFunctionDomainError()`, `SNESSetJacobianDomainError()`, `SNESVINEWTONRSLS`, `SNESVINEWTONSSLS`, `SNESSetType()`, `PETSC_NINFINITY`, `PETSC_INFINITY`
583: @*/
584: PetscErrorCode SNESVIGetVariableBounds(SNES snes, Vec *xl, Vec *xu)
585: {
586: PetscFunctionBegin;
587: PetscCheck(snes->usersetbounds, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Must set SNESVI bounds before calling SNESVIGetVariableBounds()");
588: if (xl) *xl = snes->xl;
589: if (xu) *xu = snes->xu;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: PetscErrorCode SNESSetFromOptions_VI(SNES snes, PetscOptionItems PetscOptionsObject)
594: {
595: PetscBool flg = PETSC_FALSE;
597: PetscFunctionBegin;
598: PetscOptionsHeadBegin(PetscOptionsObject, "SNES VI options");
599: PetscCall(PetscOptionsReal("-snes_vi_zero_tolerance", "Tolerance for considering x[] value to be on a bound", "None", snes->vizerotolerance, &snes->vizerotolerance, NULL));
600: PetscCall(PetscOptionsBool("-snes_vi_monitor", "Monitor all non-active variables", "SNESMonitorResidual", flg, &flg, NULL));
601: if (flg) PetscCall(SNESMonitorSet(snes, SNESMonitorVI, PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes)), NULL));
602: flg = PETSC_FALSE;
603: PetscCall(SNESMonitorSetFromOptions(snes, "-snes_vi_monitor_residual", "View residual at each iteration, using zero for active constraints", "SNESVIMonitorResidual", SNESVIMonitorResidual, NULL));
604: PetscCall(SNESMonitorSetFromOptions(snes, "-snes_vi_monitor_active", "View active set at each iteration, using zero for inactive dofs", "SNESVIMonitorActive", SNESVIMonitorActive, NULL));
605: PetscOptionsHeadEnd();
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }