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(PetscNew(&vi));
423:   snes->data = (void *)vi;

425:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
426:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
427:   PetscFunctionReturn(PETSC_SUCCESS);
428: }