Actual source code: viss.c
1: #include <../src/snes/impls/vi/ss/vissimpl.h>
3: /*@
4: SNESVIComputeMeritFunction - Evaluates the merit function for the mixed complementarity problem.
6: Input Parameter:
7: . phi - the `Vec` holding the evaluation of the semismooth function
9: Output Parameters:
10: + merit - the merit function 1/2 ||phi||^2
11: - phinorm - the two-norm of the vector, ||phi||
13: Level: developer
15: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONSSLS`, `SNESVIComputeFunction()`
16: @*/
17: PetscErrorCode SNESVIComputeMeritFunction(Vec phi, PetscReal *merit, PetscReal *phinorm)
18: {
19: PetscFunctionBegin;
20: PetscCall(VecNormBegin(phi, NORM_2, phinorm));
21: PetscCall(VecNormEnd(phi, NORM_2, phinorm));
22: *merit = 0.5 * (*phinorm) * (*phinorm);
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static inline PetscScalar Phi(PetscScalar a, PetscScalar b)
27: {
28: return a + b - PetscSqrtScalar(a * a + b * b);
29: }
31: static inline PetscScalar DPhi(PetscScalar a, PetscScalar b)
32: {
33: if ((PetscAbsScalar(a) >= 1.e-6) || (PetscAbsScalar(b) >= 1.e-6)) return 1.0 - a / PetscSqrtScalar(a * a + b * b);
34: else return .5;
35: }
37: /*@
38: SNESVIComputeFunction - Provides the function that reformulates a system of nonlinear equations in mixed complementarity form to a system of nonlinear
39: equations in semismooth form.
41: Input Parameters:
42: + snes - the `SNES` context
43: . X - current iterate
44: - functx - user defined function context
46: Output Parameter:
47: . phi - the evaluation of semismooth function at `X`
49: Level: developer
51: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONSSLS`, `SNESVIComputeMeritFunction()`
52: @*/
53: PetscErrorCode SNESVIComputeFunction(SNES snes, Vec X, Vec phi, void *functx)
54: {
55: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
56: Vec Xl = snes->xl, Xu = snes->xu, F = snes->vec_func;
57: PetscScalar *phi_arr, *f_arr, *l, *u;
58: const PetscScalar *x_arr;
59: PetscInt i, nlocal;
61: PetscFunctionBegin;
62: PetscCall((*vi->computeuserfunction)(snes, X, F, functx));
63: PetscCall(VecGetLocalSize(X, &nlocal));
64: PetscCall(VecGetArrayRead(X, &x_arr));
65: PetscCall(VecGetArray(F, &f_arr));
66: PetscCall(VecGetArray(Xl, &l));
67: PetscCall(VecGetArray(Xu, &u));
68: PetscCall(VecGetArray(phi, &phi_arr));
70: for (i = 0; i < nlocal; i++) {
71: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
72: phi_arr[i] = f_arr[i];
73: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
74: phi_arr[i] = -Phi(u[i] - x_arr[i], -f_arr[i]);
75: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
76: phi_arr[i] = Phi(x_arr[i] - l[i], f_arr[i]);
77: } else if (l[i] == u[i]) {
78: phi_arr[i] = l[i] - x_arr[i];
79: } else { /* both bounds on variable */
80: phi_arr[i] = Phi(x_arr[i] - l[i], -Phi(u[i] - x_arr[i], -f_arr[i]));
81: }
82: }
84: PetscCall(VecRestoreArrayRead(X, &x_arr));
85: PetscCall(VecRestoreArray(F, &f_arr));
86: PetscCall(VecRestoreArray(Xl, &l));
87: PetscCall(VecRestoreArray(Xu, &u));
88: PetscCall(VecRestoreArray(phi, &phi_arr));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: /*
93: SNESVIComputeBsubdifferentialVectors - Computes the diagonal shift (Da) and row scaling (Db) vectors needed for the
94: the semismooth jacobian.
95: */
96: static PetscErrorCode SNESVIComputeBsubdifferentialVectors(SNES snes, Vec X, Vec F, Mat jac, Vec Da, Vec Db)
97: {
98: PetscScalar *l, *u, *x, *f, *da, *db, da1, da2, db1, db2;
99: PetscInt i, nlocal;
101: PetscFunctionBegin;
102: PetscCall(VecGetArray(X, &x));
103: PetscCall(VecGetArray(F, &f));
104: PetscCall(VecGetArray(snes->xl, &l));
105: PetscCall(VecGetArray(snes->xu, &u));
106: PetscCall(VecGetArray(Da, &da));
107: PetscCall(VecGetArray(Db, &db));
108: PetscCall(VecGetLocalSize(X, &nlocal));
110: for (i = 0; i < nlocal; i++) {
111: if ((PetscRealPart(l[i]) <= PETSC_NINFINITY) && (PetscRealPart(u[i]) >= PETSC_INFINITY)) { /* no constraints on variable */
112: da[i] = 0;
113: db[i] = 1;
114: } else if (PetscRealPart(l[i]) <= PETSC_NINFINITY) { /* upper bound on variable only */
115: da[i] = DPhi(u[i] - x[i], -f[i]);
116: db[i] = DPhi(-f[i], u[i] - x[i]);
117: } else if (PetscRealPart(u[i]) >= PETSC_INFINITY) { /* lower bound on variable only */
118: da[i] = DPhi(x[i] - l[i], f[i]);
119: db[i] = DPhi(f[i], x[i] - l[i]);
120: } else if (l[i] == u[i]) { /* fixed variable */
121: da[i] = 1;
122: db[i] = 0;
123: } else { /* upper and lower bounds on variable */
124: da1 = DPhi(x[i] - l[i], -Phi(u[i] - x[i], -f[i]));
125: db1 = DPhi(-Phi(u[i] - x[i], -f[i]), x[i] - l[i]);
126: da2 = DPhi(u[i] - x[i], -f[i]);
127: db2 = DPhi(-f[i], u[i] - x[i]);
128: da[i] = da1 + db1 * da2;
129: db[i] = db1 * db2;
130: }
131: }
133: PetscCall(VecRestoreArray(X, &x));
134: PetscCall(VecRestoreArray(F, &f));
135: PetscCall(VecRestoreArray(snes->xl, &l));
136: PetscCall(VecRestoreArray(snes->xu, &u));
137: PetscCall(VecRestoreArray(Da, &da));
138: PetscCall(VecRestoreArray(Db, &db));
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*
143: SNESVIComputeJacobian - Computes the jacobian of the semismooth function.The Jacobian for the semismooth function is an element of the B-subdifferential of the Fischer-Burmeister function for complementarity problems.
145: Input Parameters:
146: . Da - Diagonal shift vector for the semismooth jacobian.
147: . Db - Row scaling vector for the semismooth jacobian.
149: Output Parameters:
150: . jac - semismooth jacobian
151: . jac_pre - optional preconditioning matrix
153: Note:
154: The semismooth jacobian matrix is given by
155: jac = Da + Db*jacfun
156: where Db is the row scaling matrix stored as a vector,
157: Da is the diagonal perturbation matrix stored as a vector
158: and jacfun is the jacobian of the original nonlinear function.
159: */
160: static PetscErrorCode SNESVIComputeJacobian(Mat jac, Mat jac_pre, Vec Da, Vec Db)
161: {
162: /* Do row scaling and add diagonal perturbation */
163: PetscFunctionBegin;
164: PetscCall(MatDiagonalScale(jac, Db, NULL));
165: PetscCall(MatDiagonalSet(jac, Da, ADD_VALUES));
166: if (jac != jac_pre) { /* If jac and jac_pre are different */
167: PetscCall(MatDiagonalScale(jac_pre, Db, NULL));
168: PetscCall(MatDiagonalSet(jac_pre, Da, ADD_VALUES));
169: }
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: // PetscClangLinter pragma disable: -fdoc-sowing-chars
174: /*
175: SNESVIComputeMeritFunctionGradient - Computes the gradient of the merit function psi.
177: Input Parameters:
178: phi - semismooth function.
179: H - semismooth jacobian
181: Output Parameter:
182: dpsi - merit function gradient
184: Note:
185: The merit function gradient is computed as follows
186: dpsi = H^T*phi
187: */
188: static PetscErrorCode SNESVIComputeMeritFunctionGradient(Mat H, Vec phi, Vec dpsi)
189: {
190: PetscFunctionBegin;
191: PetscCall(MatMultTranspose(H, phi, dpsi));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: static PetscErrorCode SNESSolve_VINEWTONSSLS(SNES snes)
196: {
197: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
198: PetscInt maxits, i, lits;
199: SNESLineSearchReason lssucceed;
200: PetscReal gnorm, xnorm = 0, ynorm;
201: Vec Y, X, F;
202: KSPConvergedReason kspreason;
203: DM dm;
204: DMSNES sdm;
206: PetscFunctionBegin;
207: PetscCall(SNESGetDM(snes, &dm));
208: PetscCall(DMGetDMSNES(dm, &sdm));
210: vi->computeuserfunction = sdm->ops->computefunction;
211: sdm->ops->computefunction = SNESVIComputeFunction;
213: snes->numFailures = 0;
214: snes->numLinearSolveFailures = 0;
215: snes->reason = SNES_CONVERGED_ITERATING;
217: maxits = snes->max_its; /* maximum number of iterations */
218: X = snes->vec_sol; /* solution vector */
219: F = snes->vec_func; /* residual vector */
220: Y = snes->work[0]; /* work vectors */
222: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
223: snes->iter = 0;
224: snes->norm = 0.0;
225: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
227: PetscCall(SNESVIProjectOntoBounds(snes, X));
228: PetscCall(SNESComputeFunction(snes, X, vi->phi));
229: if (snes->domainerror) {
230: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
231: sdm->ops->computefunction = vi->computeuserfunction;
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
234: /* Compute Merit function */
235: PetscCall(SNESVIComputeMeritFunction(vi->phi, &vi->merit, &vi->phinorm));
237: PetscCall(VecNormBegin(X, NORM_2, &xnorm)); /* xnorm <- ||x|| */
238: PetscCall(VecNormEnd(X, NORM_2, &xnorm));
239: SNESCheckFunctionNorm(snes, vi->merit);
241: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
242: snes->norm = vi->phinorm;
243: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
244: PetscCall(SNESLogConvergenceHistory(snes, vi->phinorm, 0));
246: /* test convergence */
247: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, vi->phinorm));
248: PetscCall(SNESMonitor(snes, 0, vi->phinorm));
249: if (snes->reason) {
250: sdm->ops->computefunction = vi->computeuserfunction;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: for (i = 0; i < maxits; i++) {
255: /* Call general purpose update function */
256: PetscTryTypeMethod(snes, update, snes->iter);
258: /* Solve J Y = Phi, where J is the semismooth jacobian */
260: /* Get the jacobian -- note that the function must be the original function for snes_fd and snes_fd_color to work for this*/
261: sdm->ops->computefunction = vi->computeuserfunction;
262: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
263: SNESCheckJacobianDomainerror(snes);
264: sdm->ops->computefunction = SNESVIComputeFunction;
266: /* Get the diagonal shift and row scaling vectors */
267: PetscCall(SNESVIComputeBsubdifferentialVectors(snes, X, F, snes->jacobian, vi->Da, vi->Db));
268: /* Compute the semismooth jacobian */
269: PetscCall(SNESVIComputeJacobian(snes->jacobian, snes->jacobian_pre, vi->Da, vi->Db));
270: /* Compute the merit function gradient */
271: PetscCall(SNESVIComputeMeritFunctionGradient(snes->jacobian, vi->phi, vi->dpsi));
272: PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
273: PetscCall(KSPSolve(snes->ksp, vi->phi, Y));
274: PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
276: if (kspreason < 0) {
277: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
278: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures));
279: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
280: break;
281: }
282: }
283: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
284: snes->linear_its += lits;
285: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
286: /*
287: if (snes->ops->precheck) {
288: PetscBool changed_y = PETSC_FALSE;
289: PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
290: }
292: if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
293: */
294: /* Compute a (scaled) negative update in the line search routine:
295: Y <- X - lambda*Y
296: and evaluate G = function(Y) (depends on the line search).
297: */
298: PetscCall(VecCopy(Y, snes->vec_sol_update));
299: ynorm = 1;
300: gnorm = vi->phinorm;
301: PetscCall(SNESLineSearchApply(snes->linesearch, X, vi->phi, &gnorm, Y));
302: PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
303: PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
304: PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)vi->phinorm, (double)gnorm, (double)ynorm, (int)lssucceed));
305: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
306: if (snes->domainerror) {
307: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
308: sdm->ops->computefunction = vi->computeuserfunction;
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
311: if (lssucceed) {
312: if (++snes->numFailures >= snes->maxFailures) {
313: PetscBool ismin;
314: snes->reason = SNES_DIVERGED_LINE_SEARCH;
315: PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, vi->phi, X, gnorm, &ismin));
316: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
317: break;
318: }
319: }
320: /* Update function and solution vectors */
321: vi->phinorm = gnorm;
322: vi->merit = 0.5 * vi->phinorm * vi->phinorm;
323: /* Monitor convergence */
324: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
325: snes->iter = i + 1;
326: snes->norm = vi->phinorm;
327: snes->xnorm = xnorm;
328: snes->ynorm = ynorm;
329: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
330: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
331: /* Test for convergence, xnorm = || X || */
332: if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
333: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, vi->phinorm));
334: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
335: if (snes->reason) break;
336: }
337: sdm->ops->computefunction = vi->computeuserfunction;
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: static PetscErrorCode SNESSetUp_VINEWTONSSLS(SNES snes)
342: {
343: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
345: PetscFunctionBegin;
346: PetscCall(SNESSetUp_VI(snes));
347: PetscCall(VecDuplicate(snes->vec_sol, &vi->dpsi));
348: PetscCall(VecDuplicate(snes->vec_sol, &vi->phi));
349: PetscCall(VecDuplicate(snes->vec_sol, &vi->Da));
350: PetscCall(VecDuplicate(snes->vec_sol, &vi->Db));
351: PetscCall(VecDuplicate(snes->vec_sol, &vi->z));
352: PetscCall(VecDuplicate(snes->vec_sol, &vi->t));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode SNESReset_VINEWTONSSLS(SNES snes)
357: {
358: SNES_VINEWTONSSLS *vi = (SNES_VINEWTONSSLS *)snes->data;
360: PetscFunctionBegin;
361: PetscCall(SNESReset_VI(snes));
362: PetscCall(VecDestroy(&vi->dpsi));
363: PetscCall(VecDestroy(&vi->phi));
364: PetscCall(VecDestroy(&vi->Da));
365: PetscCall(VecDestroy(&vi->Db));
366: PetscCall(VecDestroy(&vi->z));
367: PetscCall(VecDestroy(&vi->t));
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: static PetscErrorCode SNESSetFromOptions_VINEWTONSSLS(SNES snes, PetscOptionItems *PetscOptionsObject)
372: {
373: PetscFunctionBegin;
374: PetscCall(SNESSetFromOptions_VI(snes, PetscOptionsObject));
375: PetscOptionsHeadBegin(PetscOptionsObject, "SNES semismooth method options");
376: PetscOptionsHeadEnd();
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*MC
381: SNESVINEWTONSSLS - Semi-smooth solver for variational inequalities based on Newton's method
383: Options Database Keys:
384: + -snes_type <vinewtonssls,vinewtonrsls> a semi-smooth solver, a reduced space active set method
385: - -snes_vi_monitor - prints the number of active constraints at each iteration.
387: Level: beginner
389: Notes:
390: This family of algorithms is much like an interior point method.
392: The reduced space active set solvers `SNESVINEWTONRSLS` provide an alternative approach that does not result in extremely ill-conditioned linear systems
394: See {cite}`munson.facchinei.ea:semismooth` and {cite}`benson2006flexible`
396: .seealso: [](ch_snes), `SNESVINEWTONRSLS`, `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONRSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`
397: M*/
398: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONSSLS(SNES snes)
399: {
400: SNES_VINEWTONSSLS *vi;
401: SNESLineSearch linesearch;
403: PetscFunctionBegin;
404: snes->ops->reset = SNESReset_VINEWTONSSLS;
405: snes->ops->setup = SNESSetUp_VINEWTONSSLS;
406: snes->ops->solve = SNESSolve_VINEWTONSSLS;
407: snes->ops->destroy = SNESDestroy_VI;
408: snes->ops->setfromoptions = SNESSetFromOptions_VINEWTONSSLS;
409: snes->ops->view = NULL;
411: snes->usesksp = PETSC_TRUE;
412: snes->usesnpc = PETSC_FALSE;
414: PetscCall(SNESGetLineSearch(snes, &linesearch));
415: if (!((PetscObject)linesearch)->type_name) {
416: PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
417: PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));
418: }
420: snes->alwayscomputesfinalresidual = PETSC_FALSE;
422: PetscCall(SNESParametersInitialize(snes));
424: PetscCall(PetscNew(&vi));
425: snes->data = (void *)vi;
427: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
428: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
429: PetscFunctionReturn(PETSC_SUCCESS);
430: }