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:   PetscFunctionBegin;
 25:   *inact = vi->IS_inact;
 26:   PetscFunctionReturn(PETSC_SUCCESS);
 27: }

 29: /*
 30:     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
 31:   defined by the reduced space method.

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

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

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

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

 56:   PetscFunctionBegin;
 57:   PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
 58:   PetscCheck(isnes, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Composed SNES is missing");
 59:   PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi));
 60:   PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)dm), dmsnesvi->n, PETSC_DETERMINE, vec));
 61:   PetscCall(VecSetDM(*vec, dm));
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

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

 74:   PetscFunctionBegin;
 75:   PetscCall(PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes));
 76:   PetscCheck(isnes, PetscObjectComm((PetscObject)dm1), PETSC_ERR_PLIB, "Composed VI data structure is missing");
 77:   PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi1));
 78:   PetscCall(PetscObjectQuery((PetscObject)dm2, "VI", (PetscObject *)&isnes));
 79:   PetscCheck(isnes, PetscObjectComm((PetscObject)dm2), PETSC_ERR_PLIB, "Composed VI data structure is missing");
 80:   PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi2));

 82:   PetscCall((*dmsnesvi1->createinterpolation)(dm1, dm2, &interp, NULL));
 83:   PetscCall(MatCreateSubMatrix(interp, dmsnesvi2->inactive, dmsnesvi1->inactive, MAT_INITIAL_MATRIX, mat));
 84:   PetscCall(MatDestroy(&interp));
 85:   *vec = NULL;
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

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

103:   PetscFunctionBegin;
104:   PetscCall(PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes));
105:   PetscCheck(isnes, PetscObjectComm((PetscObject)dm1), PETSC_ERR_PLIB, "Composed VI data structure is missing");
106:   PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi1));

108:   /* get the original coarsen */
109:   PetscCall((*dmsnesvi1->coarsen)(dm1, comm, dm2));

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

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

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

120:   PetscCall(DMCreateGlobalVector(dm1, &finemarked));
121:   PetscCall(DMCreateGlobalVector(*dm2, &coarsemarked));

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

133:   PetscCall(DMCreateInjection(*dm2, dm1, &inject));
134:   PetscCall(MatRestrict(inject, finemarked, coarsemarked));
135:   PetscCall(MatDestroy(&inject));

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

154:   PetscCall(DMClearGlobalVectors(dm1));

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

158:   PetscCall(DMSetVI(*dm2, inactive));

160:   PetscCall(VecDestroy(&finemarked));
161:   PetscCall(VecDestroy(&coarsemarked));
162:   PetscCall(ISDestroy(&inactive));
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

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

179:   PetscCall(ISDestroy(&dmsnesvi->inactive));
180:   PetscCall(PetscFree(dmsnesvi));
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

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

188:     Logically Collective

190:     Input Parameters:
191: +   dm - the `DM` object
192: -   inactive - an `IS` indicating which points are currently not active

194:     Level: intermediate

196: .seealso: `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`
197: @*/
198: PetscErrorCode DMSetVI(DM dm, IS inactive)
199: {
200:   PetscContainer isnes;
201:   DM_SNESVI     *dmsnesvi;

203:   PetscFunctionBegin;
204:   if (!dm) PetscFunctionReturn(PETSC_SUCCESS);

206:   PetscCall(PetscObjectReference((PetscObject)inactive));

208:   PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
209:   if (!isnes) {
210:     PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)dm), &isnes));
211:     PetscCall(PetscContainerSetUserDestroy(isnes, (PetscErrorCode(*)(void *))DMDestroy_SNESVI));
212:     PetscCall(PetscNew(&dmsnesvi));
213:     PetscCall(PetscContainerSetPointer(isnes, (void *)dmsnesvi));
214:     PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)isnes));
215:     PetscCall(PetscContainerDestroy(&isnes));

217:     dmsnesvi->createinterpolation = dm->ops->createinterpolation;
218:     dm->ops->createinterpolation  = DMCreateInterpolation_SNESVI;
219:     dmsnesvi->coarsen             = dm->ops->coarsen;
220:     dm->ops->coarsen              = DMCoarsen_SNESVI;
221:     dmsnesvi->createglobalvector  = dm->ops->createglobalvector;
222:     dm->ops->createglobalvector   = DMCreateGlobalVector_SNESVI;
223:     dmsnesvi->createinjection     = dm->ops->createinjection;
224:     dm->ops->createinjection      = NULL;
225:     dmsnesvi->hascreateinjection  = dm->ops->hascreateinjection;
226:     dm->ops->hascreateinjection   = NULL;
227:   } else {
228:     PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi));
229:     PetscCall(ISDestroy(&dmsnesvi->inactive));
230:   }
231:   PetscCall(DMClearGlobalVectors(dm));
232:   PetscCall(ISGetLocalSize(inactive, &dmsnesvi->n));

234:   dmsnesvi->inactive = inactive;
235:   dmsnesvi->dm       = dm;
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: /*
240:      DMDestroyVI - Frees the DM_SNESVI object contained in the DM
241:          - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
242: */
243: PetscErrorCode DMDestroyVI(DM dm)
244: {
245:   PetscFunctionBegin;
246:   if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
247:   PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)NULL));
248:   PetscFunctionReturn(PETSC_SUCCESS);
249: }

251: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes, Vec X, Vec F, IS *ISact, IS *ISinact)
252: {
253:   PetscFunctionBegin;
254:   PetscCall(SNESVIGetActiveSetIS(snes, X, F, ISact));
255:   PetscCall(ISComplement(*ISact, X->map->rstart, X->map->rend, ISinact));
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

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

264:   PetscFunctionBegin;
265:   PetscCall(VecCreate(PetscObjectComm((PetscObject)snes), &v));
266:   PetscCall(VecSetSizes(v, n, PETSC_DECIDE));
267:   PetscCall(VecSetType(v, VECSTANDARD));
268:   *newv = v;
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: /* Resets the snes PC and KSP when the active set sizes change */
273: PetscErrorCode SNESVIResetPCandKSP(SNES snes, Mat Amat, Mat Pmat)
274: {
275:   KSP snesksp;

277:   PetscFunctionBegin;
278:   PetscCall(SNESGetKSP(snes, &snesksp));
279:   PetscCall(KSPReset(snesksp));
280:   PetscCall(KSPResetFromOptions(snesksp));

282:   /*
283:   KSP                    kspnew;
284:   PC                     pcnew;
285:   MatSolverType          stype;

287:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew));
288:   kspnew->pc_side = snesksp->pc_side;
289:   kspnew->rtol    = snesksp->rtol;
290:   kspnew->abstol    = snesksp->abstol;
291:   kspnew->max_it  = snesksp->max_it;
292:   PetscCall(KSPSetType(kspnew,((PetscObject)snesksp)->type_name));
293:   PetscCall(KSPGetPC(kspnew,&pcnew));
294:   PetscCall(PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name));
295:   PetscCall(PCSetOperators(kspnew->pc,Amat,Pmat));
296:   PetscCall(PCFactorGetMatSolverType(snesksp->pc,&stype));
297:   PetscCall(PCFactorSetMatSolverType(kspnew->pc,stype));
298:   PetscCall(KSPDestroy(&snesksp));
299:   snes->ksp = kspnew;
300:    PetscCall(KSPSetFromOptions(kspnew));*/
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
305:    implemented in this algorithm. It basically identifies the active constraints and does
306:    a linear solve on the other variables (those not associated with the active constraints). */
307: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
308: {
309:   SNES_VINEWTONRSLS   *vi = (SNES_VINEWTONRSLS *)snes->data;
310:   PetscInt             maxits, i, lits;
311:   SNESLineSearchReason lssucceed;
312:   PetscReal            fnorm, gnorm, xnorm = 0, ynorm;
313:   Vec                  Y, X, F;
314:   KSPConvergedReason   kspreason;
315:   KSP                  ksp;
316:   PC                   pc;

318:   PetscFunctionBegin;
319:   /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
320:   PetscCall(SNESGetKSP(snes, &ksp));
321:   PetscCall(KSPGetPC(ksp, &pc));
322:   PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));

324:   snes->numFailures            = 0;
325:   snes->numLinearSolveFailures = 0;
326:   snes->reason                 = SNES_CONVERGED_ITERATING;

328:   maxits = snes->max_its;  /* maximum number of iterations */
329:   X      = snes->vec_sol;  /* solution vector */
330:   F      = snes->vec_func; /* residual vector */
331:   Y      = snes->work[0];  /* work vectors */

333:   PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm));
334:   PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL));
335:   PetscCall(SNESLineSearchSetUp(snes->linesearch));

337:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
338:   snes->iter = 0;
339:   snes->norm = 0.0;
340:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

342:   PetscCall(SNESVIProjectOntoBounds(snes, X));
343:   PetscCall(SNESComputeFunction(snes, X, F));
344:   PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
345:   PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- ||x||  */
346:   SNESCheckFunctionNorm(snes, fnorm);
347:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
348:   snes->norm = fnorm;
349:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
350:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
351:   PetscCall(SNESMonitor(snes, 0, fnorm));

353:   /* test convergence */
354:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
355:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

357:   for (i = 0; i < maxits; i++) {
358:     IS         IS_act;    /* _act -> active set _inact -> inactive set */
359:     IS         IS_redact; /* redundant active set */
360:     VecScatter scat_act, scat_inact;
361:     PetscInt   nis_act, nis_inact;
362:     Vec        Y_act, Y_inact, F_inact;
363:     Mat        jac_inact_inact, prejac_inact_inact;
364:     PetscBool  isequal;

366:     /* Call general purpose update function */
367:     PetscTryTypeMethod(snes, update, snes->iter);
368:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
369:     SNESCheckJacobianDomainerror(snes);

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

373:     /*original
374:     PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact));
375:      */
376:     PetscCall(SNESVIGetActiveSetIS(snes, X, F, &IS_act));

378:     if (vi->checkredundancy) {
379:       PetscCall((*vi->checkredundancy)(snes, IS_act, &IS_redact, vi->ctxP));
380:       if (IS_redact) {
381:         PetscCall(ISSort(IS_redact));
382:         PetscCall(ISComplement(IS_redact, X->map->rstart, X->map->rend, &vi->IS_inact));
383:         PetscCall(ISDestroy(&IS_redact));
384:       } else {
385:         PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
386:       }
387:     } else {
388:       PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
389:     }

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

394:     if (0) { /* Dead code (temporary developer hack) */
395:       IS keptrows;
396:       PetscCall(MatFindNonzeroRows(jac_inact_inact, &keptrows));
397:       if (keptrows) {
398:         PetscInt        cnt, *nrows, k;
399:         const PetscInt *krows, *inact;
400:         PetscInt        rstart;

402:         PetscCall(MatGetOwnershipRange(jac_inact_inact, &rstart, NULL));
403:         PetscCall(MatDestroy(&jac_inact_inact));
404:         PetscCall(ISDestroy(&IS_act));

406:         PetscCall(ISGetLocalSize(keptrows, &cnt));
407:         PetscCall(ISGetIndices(keptrows, &krows));
408:         PetscCall(ISGetIndices(vi->IS_inact, &inact));
409:         PetscCall(PetscMalloc1(cnt, &nrows));
410:         for (k = 0; k < cnt; k++) nrows[k] = inact[krows[k] - rstart];
411:         PetscCall(ISRestoreIndices(keptrows, &krows));
412:         PetscCall(ISRestoreIndices(vi->IS_inact, &inact));
413:         PetscCall(ISDestroy(&keptrows));
414:         PetscCall(ISDestroy(&vi->IS_inact));

416:         PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), cnt, nrows, PETSC_OWN_POINTER, &vi->IS_inact));
417:         PetscCall(ISComplement(vi->IS_inact, F->map->rstart, F->map->rend, &IS_act));
418:         PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
419:       }
420:     }
421:     PetscCall(DMSetVI(snes->dm, vi->IS_inact));
422:     /* remove later */

424:     /*
425:     PetscCall(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm)));
426:     PetscCall(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm)));
427:     PetscCall(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X))));
428:     PetscCall(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F))));
429:     PetscCall(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact))));
430:      */

432:     /* Get sizes of active and inactive sets */
433:     PetscCall(ISGetLocalSize(IS_act, &nis_act));
434:     PetscCall(ISGetLocalSize(vi->IS_inact, &nis_inact));

436:     /* Create active and inactive set vectors */
437:     PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact));
438:     PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act));
439:     PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact));

441:     /* Create scatter contexts */
442:     PetscCall(VecScatterCreate(Y, IS_act, Y_act, NULL, &scat_act));
443:     PetscCall(VecScatterCreate(Y, vi->IS_inact, Y_inact, NULL, &scat_inact));

445:     /* Do a vec scatter to active and inactive set vectors */
446:     PetscCall(VecScatterBegin(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
447:     PetscCall(VecScatterEnd(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));

449:     PetscCall(VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
450:     PetscCall(VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));

452:     PetscCall(VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
453:     PetscCall(VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));

455:     /* Active set direction = 0 */
456:     PetscCall(VecSet(Y_act, 0));
457:     if (snes->jacobian != snes->jacobian_pre) {
458:       PetscCall(MatCreateSubMatrix(snes->jacobian_pre, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &prejac_inact_inact));
459:     } else prejac_inact_inact = jac_inact_inact;

461:     PetscCall(ISEqual(vi->IS_inact_prev, vi->IS_inact, &isequal));
462:     if (!isequal) {
463:       PetscCall(SNESVIResetPCandKSP(snes, jac_inact_inact, prejac_inact_inact));
464:       PetscCall(PCFieldSplitRestrictIS(pc, vi->IS_inact));
465:     }

467:     /*      PetscCall(ISView(vi->IS_inact,0)); */
468:     /*      PetscCall(ISView(IS_act,0));*/
469:     /*      ierr = MatView(snes->jacobian_pre,0); */

471:     PetscCall(KSPSetOperators(snes->ksp, jac_inact_inact, prejac_inact_inact));
472:     PetscCall(KSPSetUp(snes->ksp));
473:     {
474:       PC        pc;
475:       PetscBool flg;
476:       PetscCall(KSPGetPC(snes->ksp, &pc));
477:       PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg));
478:       if (flg) {
479:         KSP *subksps;
480:         PetscCall(PCFieldSplitGetSubKSP(pc, NULL, &subksps));
481:         PetscCall(KSPGetPC(subksps[0], &pc));
482:         PetscCall(PetscFree(subksps));
483:         PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &flg));
484:         if (flg) {
485:           PetscInt        n, N = 101 * 101, j, cnts[3] = {0, 0, 0};
486:           const PetscInt *ii;

488:           PetscCall(ISGetSize(vi->IS_inact, &n));
489:           PetscCall(ISGetIndices(vi->IS_inact, &ii));
490:           for (j = 0; j < n; j++) {
491:             if (ii[j] < N) cnts[0]++;
492:             else if (ii[j] < 2 * N) cnts[1]++;
493:             else if (ii[j] < 3 * N) cnts[2]++;
494:           }
495:           PetscCall(ISRestoreIndices(vi->IS_inact, &ii));

497:           PetscCall(PCBJacobiSetTotalBlocks(pc, 3, cnts));
498:         }
499:       }
500:     }

502:     PetscCall(KSPSolve(snes->ksp, F_inact, Y_inact));
503:     PetscCall(VecScatterBegin(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
504:     PetscCall(VecScatterEnd(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
505:     PetscCall(VecScatterBegin(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
506:     PetscCall(VecScatterEnd(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));

508:     PetscCall(VecDestroy(&F_inact));
509:     PetscCall(VecDestroy(&Y_act));
510:     PetscCall(VecDestroy(&Y_inact));
511:     PetscCall(VecScatterDestroy(&scat_act));
512:     PetscCall(VecScatterDestroy(&scat_inact));
513:     PetscCall(ISDestroy(&IS_act));
514:     if (!isequal) {
515:       PetscCall(ISDestroy(&vi->IS_inact_prev));
516:       PetscCall(ISDuplicate(vi->IS_inact, &vi->IS_inact_prev));
517:     }
518:     PetscCall(ISDestroy(&vi->IS_inact));
519:     PetscCall(MatDestroy(&jac_inact_inact));
520:     if (snes->jacobian != snes->jacobian_pre) PetscCall(MatDestroy(&prejac_inact_inact));

522:     PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
523:     if (kspreason < 0) {
524:       if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
525:         PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures));
526:         snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
527:         break;
528:       }
529:     }

531:     PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
532:     snes->linear_its += lits;
533:     PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
534:     /*
535:     if (snes->ops->precheck) {
536:       PetscBool changed_y = PETSC_FALSE;
537:       PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
538:     }

540:     if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
541:     */
542:     /* Compute a (scaled) negative update in the line search routine:
543:          Y <- X - lambda*Y
544:        and evaluate G = function(Y) (depends on the line search).
545:     */
546:     PetscCall(VecCopy(Y, snes->vec_sol_update));
547:     ynorm = 1;
548:     gnorm = fnorm;
549:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y));
550:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
551:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
552:     PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lssucceed));
553:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
554:     if (snes->domainerror) {
555:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
556:       PetscCall(DMDestroyVI(snes->dm));
557:       PetscFunctionReturn(PETSC_SUCCESS);
558:     }
559:     if (lssucceed) {
560:       if (++snes->numFailures >= snes->maxFailures) {
561:         PetscBool ismin;
562:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
563:         PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin));
564:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
565:         break;
566:       }
567:     }
568:     PetscCall(DMDestroyVI(snes->dm));
569:     /* Update function and solution vectors */
570:     fnorm = gnorm;
571:     /* Monitor convergence */
572:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
573:     snes->iter  = i + 1;
574:     snes->norm  = fnorm;
575:     snes->xnorm = xnorm;
576:     snes->ynorm = ynorm;
577:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
578:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
579:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
580:     /* Test for convergence, xnorm = || X || */
581:     if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
582:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
583:     if (snes->reason) break;
584:   }
585:   /* 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 */
586:   PetscCall(DMDestroyVI(snes->dm));
587:   if (i == maxits) {
588:     PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
589:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
590:   }
591:   PetscFunctionReturn(PETSC_SUCCESS);
592: }

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

597:    Logically Collective

599:    Input Parameters:
600: +  snes - the `SNESVINEWTONRSLS` context
601: .  func - the function to check of redundancies
602: -  ctx - optional context used by the function

604:    Level: advanced

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

610: .seealso: `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
611:  @*/
612: PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES, IS, IS *, void *), void *ctx)
613: {
614:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

616:   PetscFunctionBegin;
618:   vi->checkredundancy = func;
619:   vi->ctxP            = ctx;
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: #if defined(PETSC_HAVE_MATLAB)
624:   #include <engine.h>
625:   #include <mex.h>
626: typedef struct {
627:   char    *funcname;
628:   mxArray *ctx;
629: } SNESMatlabContext;

631: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, void *ctx)
632: {
633:   SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
634:   int                nlhs = 1, nrhs = 5;
635:   mxArray           *plhs[1], *prhs[5];
636:   long long int      l1 = 0, l2 = 0, ls = 0;
637:   PetscInt          *indices = NULL;

639:   PetscFunctionBegin;
643:   PetscCheckSameComm(snes, 1, is_act, 2);

645:   /* Create IS for reduced active set of size 0, its size and indices will
646:    bet set by the Matlab function */
647:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), 0, indices, PETSC_OWN_POINTER, is_redact));
648:   /* call Matlab function in ctx */
649:   PetscCall(PetscArraycpy(&ls, &snes, 1));
650:   PetscCall(PetscArraycpy(&l1, &is_act, 1));
651:   PetscCall(PetscArraycpy(&l2, is_redact, 1));
652:   prhs[0] = mxCreateDoubleScalar((double)ls);
653:   prhs[1] = mxCreateDoubleScalar((double)l1);
654:   prhs[2] = mxCreateDoubleScalar((double)l2);
655:   prhs[3] = mxCreateString(sctx->funcname);
656:   prhs[4] = sctx->ctx;
657:   PetscCall(mexCallMATLAB(nlhs, plhs, nrhs, prhs, "PetscSNESVIRedundancyCheckInternal"));
658:   PetscCall(mxGetScalar(plhs[0]));
659:   mxDestroyArray(prhs[0]);
660:   mxDestroyArray(prhs[1]);
661:   mxDestroyArray(prhs[2]);
662:   mxDestroyArray(prhs[3]);
663:   mxDestroyArray(plhs[0]);
664:   PetscFunctionReturn(PETSC_SUCCESS);
665: }

667: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
668: {
669:   SNESMatlabContext *sctx;

671:   PetscFunctionBegin;
672:   /* currently sctx is memory bleed */
673:   PetscCall(PetscNew(&sctx));
674:   PetscCall(PetscStrallocpy(func, &sctx->funcname));
675:   sctx->ctx = mxDuplicateArray(ctx);
676:   PetscCall(SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx));
677:   PetscFunctionReturn(PETSC_SUCCESS);
678: }

680: #endif

682: /*
683:    SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
684:    of the SNESVI nonlinear solver.

686:    Input Parameter:
687: .  snes - the SNES context

689:    Application Interface Routine: SNESSetUp()

691:    Note:
692:    For basic use of the SNES solvers, the user need not explicitly call
693:    SNESSetUp(), since these actions will automatically occur during
694:    the call to SNESSolve().
695:  */
696: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
697: {
698:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
699:   PetscInt          *indices;
700:   PetscInt           i, n, rstart, rend;
701:   SNESLineSearch     linesearch;

703:   PetscFunctionBegin;
704:   PetscCall(SNESSetUp_VI(snes));

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

709:   PetscCall(VecGetOwnershipRange(snes->vec_sol, &rstart, &rend));
710:   PetscCall(VecGetLocalSize(snes->vec_sol, &n));
711:   PetscCall(PetscMalloc1(n, &indices));
712:   for (i = 0; i < n; i++) indices[i] = rstart + i;
713:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev));

715:   /* set the line search functions */
716:   if (!snes->linesearch) {
717:     PetscCall(SNESGetLineSearch(snes, &linesearch));
718:     PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
719:   }
720:   PetscFunctionReturn(PETSC_SUCCESS);
721: }

723: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
724: {
725:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

727:   PetscFunctionBegin;
728:   PetscCall(SNESReset_VI(snes));
729:   PetscCall(ISDestroy(&vi->IS_inact_prev));
730:   PetscFunctionReturn(PETSC_SUCCESS);
731: }

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

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

740:    Level: beginner

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

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

752: .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
753: M*/
754: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
755: {
756:   SNES_VINEWTONRSLS *vi;
757:   SNESLineSearch     linesearch;

759:   PetscFunctionBegin;
760:   snes->ops->reset          = SNESReset_VINEWTONRSLS;
761:   snes->ops->setup          = SNESSetUp_VINEWTONRSLS;
762:   snes->ops->solve          = SNESSolve_VINEWTONRSLS;
763:   snes->ops->destroy        = SNESDestroy_VI;
764:   snes->ops->setfromoptions = SNESSetFromOptions_VI;
765:   snes->ops->view           = NULL;
766:   snes->ops->converged      = SNESConvergedDefault_VI;

768:   snes->usesksp = PETSC_TRUE;
769:   snes->usesnpc = PETSC_FALSE;

771:   PetscCall(SNESGetLineSearch(snes, &linesearch));
772:   if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
773:   PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));

775:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

777:   PetscCall(PetscNew(&vi));
778:   snes->data          = (void *)vi;
779:   vi->checkredundancy = NULL;

781:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
782:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
783:   PetscFunctionReturn(PETSC_SUCCESS);
784: }