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