Actual source code: virs.c


  2: #include <../src/snes/impls/vi/rs/virsimpl.h>
  3: #include <petsc/private/dmimpl.h>
  4: #include <petsc/private/vecimpl.h>

  6: /*@
  7:    SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear
  8:      system is solved on)

 10:    Input parameter:
 11: .  snes - the `SNES` context

 13:    Output parameter:
 14: .  inact - inactive set index set

 16:    Level: advanced

 18: .seealso: `SNESVINEWTONRSLS`
 19: @*/
 20: PetscErrorCode SNESVIGetInactiveSet(SNES snes, IS *inact)
 21: {
 22:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

 24:   *inact = vi->IS_inact;
 25:   return 0;
 26: }

 28: /*
 29:     Provides a wrapper to a DM to allow it to be used to generated the interpolation/restriction from the DM for the smaller matrices and vectors
 30:   defined by the reduced space method.

 32:     Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints.
 33: */
 34: typedef struct {
 35:   PetscInt n; /* size of vectors in the reduced DM space */
 36:   IS       inactive;

 38:   PetscErrorCode (*createinterpolation)(DM, DM, Mat *, Vec *); /* DM's original routines */
 39:   PetscErrorCode (*coarsen)(DM, MPI_Comm, DM *);
 40:   PetscErrorCode (*createglobalvector)(DM, Vec *);
 41:   PetscErrorCode (*createinjection)(DM, DM, Mat *);
 42:   PetscErrorCode (*hascreateinjection)(DM, PetscBool *);

 44:   DM dm; /* when destroying this object we need to reset the above function into the base DM */
 45: } DM_SNESVI;

 47: /*
 48:      DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
 49: */
 50: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm, Vec *vec)
 51: {
 52:   PetscContainer isnes;
 53:   DM_SNESVI     *dmsnesvi;

 55:   PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes);
 57:   PetscContainerGetPointer(isnes, (void **)&dmsnesvi);
 58:   VecCreateMPI(PetscObjectComm((PetscObject)dm), dmsnesvi->n, PETSC_DETERMINE, vec);
 59:   VecSetDM(*vec, dm);
 60:   return 0;
 61: }

 63: /*
 64:      DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
 65: */
 66: PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1, DM dm2, Mat *mat, Vec *vec)
 67: {
 68:   PetscContainer isnes;
 69:   DM_SNESVI     *dmsnesvi1, *dmsnesvi2;
 70:   Mat            interp;

 72:   PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes);
 74:   PetscContainerGetPointer(isnes, (void **)&dmsnesvi1);
 75:   PetscObjectQuery((PetscObject)dm2, "VI", (PetscObject *)&isnes);
 77:   PetscContainerGetPointer(isnes, (void **)&dmsnesvi2);

 79:   (*dmsnesvi1->createinterpolation)(dm1, dm2, &interp, NULL);
 80:   MatCreateSubMatrix(interp, dmsnesvi2->inactive, dmsnesvi1->inactive, MAT_INITIAL_MATRIX, mat);
 81:   MatDestroy(&interp);
 82:   *vec = NULL;
 83:   return 0;
 84: }

 86: /*
 87:      DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
 88: */
 89: PetscErrorCode DMCoarsen_SNESVI(DM dm1, MPI_Comm comm, DM *dm2)
 90: {
 91:   PetscContainer  isnes;
 92:   DM_SNESVI      *dmsnesvi1;
 93:   Vec             finemarked, coarsemarked;
 94:   IS              inactive;
 95:   Mat             inject;
 96:   const PetscInt *index;
 97:   PetscInt        n, k, cnt = 0, rstart, *coarseindex;
 98:   PetscScalar    *marked;

100:   PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes);
102:   PetscContainerGetPointer(isnes, (void **)&dmsnesvi1);

104:   /* get the original coarsen */
105:   (*dmsnesvi1->coarsen)(dm1, comm, dm2);

107:   /* not sure why this extra reference is needed, but without the dm2 disappears too early */
108:   /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
109:   /*  PetscObjectReference((PetscObject)*dm2);*/

111:   /* need to set back global vectors in order to use the original injection */
112:   DMClearGlobalVectors(dm1);

114:   dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;

116:   DMCreateGlobalVector(dm1, &finemarked);
117:   DMCreateGlobalVector(*dm2, &coarsemarked);

119:   /*
120:      fill finemarked with locations of inactive points
121:   */
122:   ISGetIndices(dmsnesvi1->inactive, &index);
123:   ISGetLocalSize(dmsnesvi1->inactive, &n);
124:   VecSet(finemarked, 0.0);
125:   for (k = 0; k < n; k++) VecSetValue(finemarked, index[k], 1.0, INSERT_VALUES);
126:   VecAssemblyBegin(finemarked);
127:   VecAssemblyEnd(finemarked);

129:   DMCreateInjection(*dm2, dm1, &inject);
130:   MatRestrict(inject, finemarked, coarsemarked);
131:   MatDestroy(&inject);

133:   /*
134:      create index set list of coarse inactive points from coarsemarked
135:   */
136:   VecGetLocalSize(coarsemarked, &n);
137:   VecGetOwnershipRange(coarsemarked, &rstart, NULL);
138:   VecGetArray(coarsemarked, &marked);
139:   for (k = 0; k < n; k++) {
140:     if (marked[k] != 0.0) cnt++;
141:   }
142:   PetscMalloc1(cnt, &coarseindex);
143:   cnt = 0;
144:   for (k = 0; k < n; k++) {
145:     if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
146:   }
147:   VecRestoreArray(coarsemarked, &marked);
148:   ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked), cnt, coarseindex, PETSC_OWN_POINTER, &inactive);

150:   DMClearGlobalVectors(dm1);

152:   dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;

154:   DMSetVI(*dm2, inactive);

156:   VecDestroy(&finemarked);
157:   VecDestroy(&coarsemarked);
158:   ISDestroy(&inactive);
159:   return 0;
160: }

162: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
163: {
164:   /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
165:   dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
166:   dmsnesvi->dm->ops->coarsen             = dmsnesvi->coarsen;
167:   dmsnesvi->dm->ops->createglobalvector  = dmsnesvi->createglobalvector;
168:   dmsnesvi->dm->ops->createinjection     = dmsnesvi->createinjection;
169:   dmsnesvi->dm->ops->hascreateinjection  = dmsnesvi->hascreateinjection;
170:   /* need to clear out this vectors because some of them may not have a reference to the DM
171:     but they are counted as having references to the DM in DMDestroy() */
172:   DMClearGlobalVectors(dmsnesvi->dm);

174:   ISDestroy(&dmsnesvi->inactive);
175:   PetscFree(dmsnesvi);
176:   return 0;
177: }

179: /*@
180:      DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
181:                be restricted to only those variables NOT associated with active constraints.

183:     Logically Collective

185:     Input Parameters:
186: +   dm - the `DM` object
187: -   inactive - an `IS` indicating which points are currently not active

189:     Level: intermediate

191: .seealso: `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`
192: @*/
193: PetscErrorCode DMSetVI(DM dm, IS inactive)
194: {
195:   PetscContainer isnes;
196:   DM_SNESVI     *dmsnesvi;

198:   if (!dm) return 0;

200:   PetscObjectReference((PetscObject)inactive);

202:   PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes);
203:   if (!isnes) {
204:     PetscContainerCreate(PetscObjectComm((PetscObject)dm), &isnes);
205:     PetscContainerSetUserDestroy(isnes, (PetscErrorCode(*)(void *))DMDestroy_SNESVI);
206:     PetscNew(&dmsnesvi);
207:     PetscContainerSetPointer(isnes, (void *)dmsnesvi);
208:     PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)isnes);
209:     PetscContainerDestroy(&isnes);

211:     dmsnesvi->createinterpolation = dm->ops->createinterpolation;
212:     dm->ops->createinterpolation  = DMCreateInterpolation_SNESVI;
213:     dmsnesvi->coarsen             = dm->ops->coarsen;
214:     dm->ops->coarsen              = DMCoarsen_SNESVI;
215:     dmsnesvi->createglobalvector  = dm->ops->createglobalvector;
216:     dm->ops->createglobalvector   = DMCreateGlobalVector_SNESVI;
217:     dmsnesvi->createinjection     = dm->ops->createinjection;
218:     dm->ops->createinjection      = NULL;
219:     dmsnesvi->hascreateinjection  = dm->ops->hascreateinjection;
220:     dm->ops->hascreateinjection   = NULL;
221:   } else {
222:     PetscContainerGetPointer(isnes, (void **)&dmsnesvi);
223:     ISDestroy(&dmsnesvi->inactive);
224:   }
225:   DMClearGlobalVectors(dm);
226:   ISGetLocalSize(inactive, &dmsnesvi->n);

228:   dmsnesvi->inactive = inactive;
229:   dmsnesvi->dm       = dm;
230:   return 0;
231: }

233: /*
234:      DMDestroyVI - Frees the DM_SNESVI object contained in the DM
235:          - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
236: */
237: PetscErrorCode DMDestroyVI(DM dm)
238: {
239:   if (!dm) return 0;
240:   PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)NULL);
241:   return 0;
242: }

244: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes, Vec X, Vec F, IS *ISact, IS *ISinact)
245: {
246:   SNESVIGetActiveSetIS(snes, X, F, ISact);
247:   ISComplement(*ISact, X->map->rstart, X->map->rend, ISinact);
248:   return 0;
249: }

251: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
252: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
253: {
254:   Vec v;

256:   VecCreate(PetscObjectComm((PetscObject)snes), &v);
257:   VecSetSizes(v, n, PETSC_DECIDE);
258:   VecSetType(v, VECSTANDARD);
259:   *newv = v;
260:   return 0;
261: }

263: /* Resets the snes PC and KSP when the active set sizes change */
264: PetscErrorCode SNESVIResetPCandKSP(SNES snes, Mat Amat, Mat Pmat)
265: {
266:   KSP snesksp;

268:   SNESGetKSP(snes, &snesksp);
269:   KSPReset(snesksp);
270:   KSPResetFromOptions(snesksp);

272:   /*
273:   KSP                    kspnew;
274:   PC                     pcnew;
275:   MatSolverType          stype;

277:   KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
278:   kspnew->pc_side = snesksp->pc_side;
279:   kspnew->rtol    = snesksp->rtol;
280:   kspnew->abstol    = snesksp->abstol;
281:   kspnew->max_it  = snesksp->max_it;
282:   KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
283:   KSPGetPC(kspnew,&pcnew);
284:   PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
285:   PCSetOperators(kspnew->pc,Amat,Pmat);
286:   PCFactorGetMatSolverType(snesksp->pc,&stype);
287:   PCFactorSetMatSolverType(kspnew->pc,stype);
288:   KSPDestroy(&snesksp);
289:   snes->ksp = kspnew;
290:    KSPSetFromOptions(kspnew);*/
291:   return 0;
292: }

294: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
295:    implemented in this algorithm. It basically identifies the active constraints and does
296:    a linear solve on the other variables (those not associated with the active constraints). */
297: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
298: {
299:   SNES_VINEWTONRSLS   *vi = (SNES_VINEWTONRSLS *)snes->data;
300:   PetscInt             maxits, i, lits;
301:   SNESLineSearchReason lssucceed;
302:   PetscReal            fnorm, gnorm, xnorm = 0, ynorm;
303:   Vec                  Y, X, F;
304:   KSPConvergedReason   kspreason;
305:   KSP                  ksp;
306:   PC                   pc;

308:   /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
309:   SNESGetKSP(snes, &ksp);
310:   KSPGetPC(ksp, &pc);
311:   PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH);

313:   snes->numFailures            = 0;
314:   snes->numLinearSolveFailures = 0;
315:   snes->reason                 = SNES_CONVERGED_ITERATING;

317:   maxits = snes->max_its;  /* maximum number of iterations */
318:   X      = snes->vec_sol;  /* solution vector */
319:   F      = snes->vec_func; /* residual vector */
320:   Y      = snes->work[0];  /* work vectors */

322:   SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
323:   SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
324:   SNESLineSearchSetUp(snes->linesearch);

326:   PetscObjectSAWsTakeAccess((PetscObject)snes);
327:   snes->iter = 0;
328:   snes->norm = 0.0;
329:   PetscObjectSAWsGrantAccess((PetscObject)snes);

331:   SNESVIProjectOntoBounds(snes, X);
332:   SNESComputeFunction(snes, X, F);
333:   SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
334:   VecNorm(X, NORM_2, &xnorm); /* xnorm <- ||x||  */
335:   SNESCheckFunctionNorm(snes, fnorm);
336:   PetscObjectSAWsTakeAccess((PetscObject)snes);
337:   snes->norm = fnorm;
338:   PetscObjectSAWsGrantAccess((PetscObject)snes);
339:   SNESLogConvergenceHistory(snes, fnorm, 0);
340:   SNESMonitor(snes, 0, fnorm);

342:   /* test convergence */
343:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
344:   if (snes->reason) return 0;

346:   for (i = 0; i < maxits; i++) {
347:     IS         IS_act;    /* _act -> active set _inact -> inactive set */
348:     IS         IS_redact; /* redundant active set */
349:     VecScatter scat_act, scat_inact;
350:     PetscInt   nis_act, nis_inact;
351:     Vec        Y_act, Y_inact, F_inact;
352:     Mat        jac_inact_inact, prejac_inact_inact;
353:     PetscBool  isequal;

355:     /* Call general purpose update function */
356:     PetscTryTypeMethod(snes, update, snes->iter);
357:     SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
358:     SNESCheckJacobianDomainerror(snes);

360:     /* Create active and inactive index sets */

362:     /*original
363:     SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
364:      */
365:     SNESVIGetActiveSetIS(snes, X, F, &IS_act);

367:     if (vi->checkredundancy) {
368:       (*vi->checkredundancy)(snes, IS_act, &IS_redact, vi->ctxP);
369:       if (IS_redact) {
370:         ISSort(IS_redact);
371:         ISComplement(IS_redact, X->map->rstart, X->map->rend, &vi->IS_inact);
372:         ISDestroy(&IS_redact);
373:       } else {
374:         ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact);
375:       }
376:     } else {
377:       ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact);
378:     }

380:     /* Create inactive set submatrix */
381:     MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact);

383:     if (0) { /* Dead code (temporary developer hack) */
384:       IS keptrows;
385:       MatFindNonzeroRows(jac_inact_inact, &keptrows);
386:       if (keptrows) {
387:         PetscInt        cnt, *nrows, k;
388:         const PetscInt *krows, *inact;
389:         PetscInt        rstart;

391:         MatGetOwnershipRange(jac_inact_inact, &rstart, NULL);
392:         MatDestroy(&jac_inact_inact);
393:         ISDestroy(&IS_act);

395:         ISGetLocalSize(keptrows, &cnt);
396:         ISGetIndices(keptrows, &krows);
397:         ISGetIndices(vi->IS_inact, &inact);
398:         PetscMalloc1(cnt, &nrows);
399:         for (k = 0; k < cnt; k++) nrows[k] = inact[krows[k] - rstart];
400:         ISRestoreIndices(keptrows, &krows);
401:         ISRestoreIndices(vi->IS_inact, &inact);
402:         ISDestroy(&keptrows);
403:         ISDestroy(&vi->IS_inact);

405:         ISCreateGeneral(PetscObjectComm((PetscObject)snes), cnt, nrows, PETSC_OWN_POINTER, &vi->IS_inact);
406:         ISComplement(vi->IS_inact, F->map->rstart, F->map->rend, &IS_act);
407:         MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact);
408:       }
409:     }
410:     DMSetVI(snes->dm, vi->IS_inact);
411:     /* remove later */

413:     /*
414:     VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
415:     VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
416:     VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
417:     VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
418:     ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
419:      */

421:     /* Get sizes of active and inactive sets */
422:     ISGetLocalSize(IS_act, &nis_act);
423:     ISGetLocalSize(vi->IS_inact, &nis_inact);

425:     /* Create active and inactive set vectors */
426:     SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact);
427:     SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act);
428:     SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact);

430:     /* Create scatter contexts */
431:     VecScatterCreate(Y, IS_act, Y_act, NULL, &scat_act);
432:     VecScatterCreate(Y, vi->IS_inact, Y_inact, NULL, &scat_inact);

434:     /* Do a vec scatter to active and inactive set vectors */
435:     VecScatterBegin(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD);
436:     VecScatterEnd(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD);

438:     VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD);
439:     VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD);

441:     VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD);
442:     VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD);

444:     /* Active set direction = 0 */
445:     VecSet(Y_act, 0);
446:     if (snes->jacobian != snes->jacobian_pre) {
447:       MatCreateSubMatrix(snes->jacobian_pre, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &prejac_inact_inact);
448:     } else prejac_inact_inact = jac_inact_inact;

450:     ISEqual(vi->IS_inact_prev, vi->IS_inact, &isequal);
451:     if (!isequal) {
452:       SNESVIResetPCandKSP(snes, jac_inact_inact, prejac_inact_inact);
453:       PCFieldSplitRestrictIS(pc, vi->IS_inact);
454:     }

456:     /*      ISView(vi->IS_inact,0); */
457:     /*      ISView(IS_act,0);*/
458:     /*      MatView(snes->jacobian_pre,0); */

460:     KSPSetOperators(snes->ksp, jac_inact_inact, prejac_inact_inact);
461:     KSPSetUp(snes->ksp);
462:     {
463:       PC        pc;
464:       PetscBool flg;
465:       KSPGetPC(snes->ksp, &pc);
466:       PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg);
467:       if (flg) {
468:         KSP *subksps;
469:         PCFieldSplitGetSubKSP(pc, NULL, &subksps);
470:         KSPGetPC(subksps[0], &pc);
471:         PetscFree(subksps);
472:         PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &flg);
473:         if (flg) {
474:           PetscInt        n, N = 101 * 101, j, cnts[3] = {0, 0, 0};
475:           const PetscInt *ii;

477:           ISGetSize(vi->IS_inact, &n);
478:           ISGetIndices(vi->IS_inact, &ii);
479:           for (j = 0; j < n; j++) {
480:             if (ii[j] < N) cnts[0]++;
481:             else if (ii[j] < 2 * N) cnts[1]++;
482:             else if (ii[j] < 3 * N) cnts[2]++;
483:           }
484:           ISRestoreIndices(vi->IS_inact, &ii);

486:           PCBJacobiSetTotalBlocks(pc, 3, cnts);
487:         }
488:       }
489:     }

491:     KSPSolve(snes->ksp, F_inact, Y_inact);
492:     VecScatterBegin(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE);
493:     VecScatterEnd(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE);
494:     VecScatterBegin(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE);
495:     VecScatterEnd(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE);

497:     VecDestroy(&F_inact);
498:     VecDestroy(&Y_act);
499:     VecDestroy(&Y_inact);
500:     VecScatterDestroy(&scat_act);
501:     VecScatterDestroy(&scat_inact);
502:     ISDestroy(&IS_act);
503:     if (!isequal) {
504:       ISDestroy(&vi->IS_inact_prev);
505:       ISDuplicate(vi->IS_inact, &vi->IS_inact_prev);
506:     }
507:     ISDestroy(&vi->IS_inact);
508:     MatDestroy(&jac_inact_inact);
509:     if (snes->jacobian != snes->jacobian_pre) MatDestroy(&prejac_inact_inact);

511:     KSPGetConvergedReason(snes->ksp, &kspreason);
512:     if (kspreason < 0) {
513:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
514:         PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures);
515:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
516:         break;
517:       }
518:     }

520:     KSPGetIterationNumber(snes->ksp, &lits);
521:     snes->linear_its += lits;
522:     PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits);
523:     /*
524:     if (snes->ops->precheck) {
525:       PetscBool changed_y = PETSC_FALSE;
526:       PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
527:     }

529:     if (PetscLogPrintInfo) SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
530:     */
531:     /* Compute a (scaled) negative update in the line search routine:
532:          Y <- X - lambda*Y
533:        and evaluate G = function(Y) (depends on the line search).
534:     */
535:     VecCopy(Y, snes->vec_sol_update);
536:     ynorm = 1;
537:     gnorm = fnorm;
538:     SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
539:     SNESLineSearchGetReason(snes->linesearch, &lssucceed);
540:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
541:     PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lssucceed);
542:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
543:     if (snes->domainerror) {
544:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
545:       DMDestroyVI(snes->dm);
546:       return 0;
547:     }
548:     if (lssucceed) {
549:       if (++snes->numFailures >= snes->maxFailures) {
550:         PetscBool ismin;
551:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
552:         SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin);
553:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
554:         break;
555:       }
556:     }
557:     DMDestroyVI(snes->dm);
558:     /* Update function and solution vectors */
559:     fnorm = gnorm;
560:     /* Monitor convergence */
561:     PetscObjectSAWsTakeAccess((PetscObject)snes);
562:     snes->iter  = i + 1;
563:     snes->norm  = fnorm;
564:     snes->xnorm = xnorm;
565:     snes->ynorm = ynorm;
566:     PetscObjectSAWsGrantAccess((PetscObject)snes);
567:     SNESLogConvergenceHistory(snes, snes->norm, lits);
568:     SNESMonitor(snes, snes->iter, snes->norm);
569:     /* Test for convergence, xnorm = || X || */
570:     if (snes->ops->converged != SNESConvergedSkip) VecNorm(X, NORM_2, &xnorm);
571:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
572:     if (snes->reason) break;
573:   }
574:   /* make sure that the VI information attached to the DM is removed if the for loop above was broken early due to some exceptional conditional */
575:   DMDestroyVI(snes->dm);
576:   if (i == maxits) {
577:     PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
578:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
579:   }
580:   return 0;
581: }

583: /*@C
584:    SNESVISetRedundancyCheck - Provide a function to check for any redundancy in the VI active set

586:    Logically Collective

588:    Input Parameters:
589: +  snes - the `SNESVINEWTONRSLS` context
590: .  func - the function to check of redundancies
591: -  ctx - optional context used by the function

593:    Level: advanced

595:    Note:
596:    Sometimes the inactive set will result in a non-singular sub-Jacobian problem that needs to be solved, this allows the user,
597:    when they know more about their specific problem to provide a function that removes the redundancy that results in the singular linear system

599: .seealso: `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
600:  @*/
601: PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES, IS, IS *, void *), void *ctx)
602: {
603:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

606:   vi->checkredundancy = func;
607:   vi->ctxP            = ctx;
608:   return 0;
609: }

611: #if defined(PETSC_HAVE_MATLAB)
612:   #include <engine.h>
613:   #include <mex.h>
614: typedef struct {
615:   char    *funcname;
616:   mxArray *ctx;
617: } SNESMatlabContext;

619: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, void *ctx)
620: {
621:   SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
622:   int                nlhs = 1, nrhs = 5;
623:   mxArray           *plhs[1], *prhs[5];
624:   long long int      l1 = 0, l2 = 0, ls = 0;
625:   PetscInt          *indices = NULL;


632:   /* Create IS for reduced active set of size 0, its size and indices will
633:    bet set by the Matlab function */
634:   ISCreateGeneral(PetscObjectComm((PetscObject)snes), 0, indices, PETSC_OWN_POINTER, is_redact);
635:   /* call Matlab function in ctx */
636:   PetscArraycpy(&ls, &snes, 1);
637:   PetscArraycpy(&l1, &is_act, 1);
638:   PetscArraycpy(&l2, is_redact, 1);
639:   prhs[0] = mxCreateDoubleScalar((double)ls);
640:   prhs[1] = mxCreateDoubleScalar((double)l1);
641:   prhs[2] = mxCreateDoubleScalar((double)l2);
642:   prhs[3] = mxCreateString(sctx->funcname);
643:   prhs[4] = sctx->ctx;
644:   mexCallMATLAB(nlhs, plhs, nrhs, prhs, "PetscSNESVIRedundancyCheckInternal");
645:   mxGetScalar(plhs[0]);
646:   mxDestroyArray(prhs[0]);
647:   mxDestroyArray(prhs[1]);
648:   mxDestroyArray(prhs[2]);
649:   mxDestroyArray(prhs[3]);
650:   mxDestroyArray(plhs[0]);
651:   return 0;
652: }

654: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
655: {
656:   SNESMatlabContext *sctx;

658:   /* currently sctx is memory bleed */
659:   PetscNew(&sctx);
660:   PetscStrallocpy(func, &sctx->funcname);
661:   sctx->ctx = mxDuplicateArray(ctx);
662:   SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx);
663:   return 0;
664: }

666: #endif

668: /*
669:    SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
670:    of the SNESVI nonlinear solver.

672:    Input Parameter:
673: .  snes - the SNES context

675:    Application Interface Routine: SNESSetUp()

677:    Note:
678:    For basic use of the SNES solvers, the user need not explicitly call
679:    SNESSetUp(), since these actions will automatically occur during
680:    the call to SNESSolve().
681:  */
682: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
683: {
684:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
685:   PetscInt          *indices;
686:   PetscInt           i, n, rstart, rend;
687:   SNESLineSearch     linesearch;

689:   SNESSetUp_VI(snes);

691:   /* Set up previous active index set for the first snes solve
692:    vi->IS_inact_prev = 0,1,2,....N */

694:   VecGetOwnershipRange(snes->vec_sol, &rstart, &rend);
695:   VecGetLocalSize(snes->vec_sol, &n);
696:   PetscMalloc1(n, &indices);
697:   for (i = 0; i < n; i++) indices[i] = rstart + i;
698:   ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev);

700:   /* set the line search functions */
701:   if (!snes->linesearch) {
702:     SNESGetLineSearch(snes, &linesearch);
703:     SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
704:   }
705:   return 0;
706: }

708: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
709: {
710:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

712:   SNESReset_VI(snes);
713:   ISDestroy(&vi->IS_inact_prev);
714:   return 0;
715: }

717: /*MC
718:       SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method

720:    Options Database Keys:
721: +   -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
722: -   -snes_vi_monitor - prints the number of active constraints at each iteration.

724:    Level: beginner

726:    References:
727: .  * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
728:      Applications, Optimization Methods and Software, 21 (2006).

730:    Note:
731:    At each set of this methods the algorithm produces an inactive set of variables that are constrained to their current values
732:    (because changing these values would result in those variables no longer satisfying the inequality constraints)
733:    and produces a step direction by solving the linear system arising from the Jacobian with the inactive variables removed. In other
734:    words on a reduced space of the solution space. Based on the Newton update it then adjusts the inactive sep for the next iteration.

736: .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTRDC`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
737: M*/
738: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
739: {
740:   SNES_VINEWTONRSLS *vi;
741:   SNESLineSearch     linesearch;

743:   snes->ops->reset          = SNESReset_VINEWTONRSLS;
744:   snes->ops->setup          = SNESSetUp_VINEWTONRSLS;
745:   snes->ops->solve          = SNESSolve_VINEWTONRSLS;
746:   snes->ops->destroy        = SNESDestroy_VI;
747:   snes->ops->setfromoptions = SNESSetFromOptions_VI;
748:   snes->ops->view           = NULL;
749:   snes->ops->converged      = SNESConvergedDefault_VI;

751:   snes->usesksp = PETSC_TRUE;
752:   snes->usesnpc = PETSC_FALSE;

754:   SNESGetLineSearch(snes, &linesearch);
755:   if (!((PetscObject)linesearch)->type_name) SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
756:   SNESLineSearchBTSetAlpha(linesearch, 0.0);

758:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

760:   PetscNew(&vi);
761:   snes->data          = (void *)vi;
762:   vi->checkredundancy = NULL;

764:   PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI);
765:   PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI);
766:   return 0;
767: }