Actual source code: virs.c

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

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

      9:   Input Parameter:
     10: . snes - the `SNES` context

     12:   Output Parameter:
     13: . inact - inactive set index set

     15:   Level: advanced

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

     23:   PetscFunctionBegin;
     24:   *inact = vi->IS_inact;
     25:   PetscFunctionReturn(PETSC_SUCCESS);
     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: static PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm, Vec *vec)
     51: {
     52:   PetscContainer isnes;
     53:   DM_SNESVI     *dmsnesvi;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    153:   PetscCall(DMClearGlobalVectors(dm1));

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

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

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

    165: static PetscErrorCode DMDestroy_SNESVI(void **ctx)
    166: {
    167:   DM_SNESVI *dmsnesvi = (DM_SNESVI *)*ctx;

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

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

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

    189:   Logically Collective

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

    195:   Level: intermediate

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

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

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

    209:   PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
    210:   if (!isnes) {
    211:     PetscCall(PetscNew(&dmsnesvi));
    212:     PetscCall(PetscObjectContainerCompose((PetscObject)dm, "VI", dmsnesvi, DMDestroy_SNESVI));

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

    231:   dmsnesvi->inactive = inactive;
    232:   dmsnesvi->dm       = dm;
    233:   PetscFunctionReturn(PETSC_SUCCESS);
    234: }

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

    248: /* Create active and inactive set vectors. The local size of this vector is set and PETSc computes the global size */
    249: static PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
    250: {
    251:   Vec v;

    253:   PetscFunctionBegin;
    254:   PetscCall(VecCreate(PetscObjectComm((PetscObject)snes), &v));
    255:   PetscCall(VecSetSizes(v, n, PETSC_DECIDE));
    256:   PetscCall(VecSetType(v, VECSTANDARD));
    257:   *newv = v;
    258:   PetscFunctionReturn(PETSC_SUCCESS);
    259: }

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

    266:   PetscFunctionBegin;
    267:   PetscCall(SNESGetKSP(snes, &snesksp));
    268:   PetscCall(KSPReset(snesksp));
    269:   PetscCall(KSPResetFromOptions(snesksp));

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

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

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

    307:   PetscFunctionBegin;
    308:   /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
    309:   PetscCall(SNESGetKSP(snes, &ksp));
    310:   PetscCall(KSPGetPC(ksp, &pc));
    311:   PetscCall(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:   PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm, SNESVIComputeInactiveSetFtY));
    323:   PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL));
    324:   PetscCall(SNESLineSearchSetUp(snes->linesearch));

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

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

    341:   /* test convergence */
    342:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
    343:   PetscCall(SNESMonitor(snes, 0, fnorm));
    344:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

    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:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
    358:     SNESCheckJacobianDomainerror(snes);

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

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

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

    380:     /* Create inactive set submatrix */
    381:     PetscCall(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:       PetscCall(MatFindNonzeroRows(jac_inact_inact, &keptrows));
    386:       if (keptrows) {
    387:         PetscInt        cnt, *nrows, k;
    388:         const PetscInt *krows, *inact;
    389:         PetscInt        rstart;

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

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

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

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

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

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

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

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

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

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

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

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

    455:     /*      PetscCall(ISView(vi->IS_inact,0)); */
    456:     /*      PetscCall(ISView(IS_act,0));*/
    457:     /*      ierr = MatView(snes->jacobian_pre,0); */

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

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

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

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

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

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

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

    528:     if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
    529:     */
    530:     /* Compute a (scaled) negative update in the line search routine:
    531:          Y <- X - lambda*Y
    532:        and evaluate G = function(Y) (depends on the line search).
    533:     */
    534:     PetscCall(VecCopy(Y, snes->vec_sol_update));
    535:     ynorm = 1;
    536:     gnorm = fnorm;
    537:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y));
    538:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
    539:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
    540:     PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lssucceed));
    541:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
    542:     if (snes->domainerror) {
    543:       snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
    544:       PetscCall(DMDestroyVI(snes->dm));
    545:       PetscFunctionReturn(PETSC_SUCCESS);
    546:     }
    547:     if (lssucceed) {
    548:       if (++snes->numFailures >= snes->maxFailures) {
    549:         PetscBool ismin;
    550:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
    551:         PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin));
    552:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
    553:         break;
    554:       }
    555:     }
    556:     PetscCall(DMDestroyVI(snes->dm));
    557:     /* Update function and solution vectors */
    558:     fnorm = gnorm;
    559:     /* Monitor convergence */
    560:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    561:     snes->iter  = i + 1;
    562:     snes->norm  = fnorm;
    563:     snes->xnorm = xnorm;
    564:     snes->ynorm = ynorm;
    565:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    566:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
    567:     /* Test for convergence, xnorm = || X || */
    568:     if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
    569:     PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
    570:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
    571:     if (snes->reason) break;
    572:   }
    573:   /* 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 */
    574:   PetscCall(DMDestroyVI(snes->dm));
    575:   PetscFunctionReturn(PETSC_SUCCESS);
    576: }

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

    581:   Logically Collective

    583:   Input Parameters:
    584: + snes - the `SNESVINEWTONRSLS` context
    585: . func - the function to check of redundancies
    586: - ctx  - optional context used by the function

    588:   Level: advanced

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

    594: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
    595:  @*/
    596: PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES, IS, IS *, void *), void *ctx)
    597: {
    598:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

    600:   PetscFunctionBegin;
    602:   vi->checkredundancy = func;
    603:   vi->ctxP            = ctx;
    604:   PetscFunctionReturn(PETSC_SUCCESS);
    605: }

    607: #if defined(PETSC_HAVE_MATLAB)
    608:   #include <engine.h>
    609:   #include <mex.h>
    610: typedef struct {
    611:   char    *funcname;
    612:   mxArray *ctx;
    613: } SNESMatlabContext;

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

    623:   PetscFunctionBegin;
    626:   PetscAssertPointer(is_redact, 3);
    627:   PetscCheckSameComm(snes, 1, is_act, 2);

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

    651: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
    652: {
    653:   SNESMatlabContext *sctx;

    655:   PetscFunctionBegin;
    656:   /* currently sctx is memory bleed */
    657:   PetscCall(PetscNew(&sctx));
    658:   PetscCall(PetscStrallocpy(func, &sctx->funcname));
    659:   sctx->ctx = mxDuplicateArray(ctx);
    660:   PetscCall(SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx));
    661:   PetscFunctionReturn(PETSC_SUCCESS);
    662: }

    664: #endif

    666: static PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
    667: {
    668:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
    669:   PetscInt          *indices;
    670:   PetscInt           i, n, rstart, rend;
    671:   SNESLineSearch     linesearch;

    673:   PetscFunctionBegin;
    674:   PetscCall(SNESSetUp_VI(snes));

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

    679:   PetscCall(VecGetOwnershipRange(snes->work[0], &rstart, &rend));
    680:   PetscCall(VecGetLocalSize(snes->work[0], &n));
    681:   PetscCall(PetscMalloc1(n, &indices));
    682:   for (i = 0; i < n; i++) indices[i] = rstart + i;
    683:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev));

    685:   /* set the line search functions */
    686:   if (!snes->linesearch) {
    687:     PetscCall(SNESGetLineSearch(snes, &linesearch));
    688:     PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
    689:   }
    690:   PetscFunctionReturn(PETSC_SUCCESS);
    691: }

    693: static PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
    694: {
    695:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

    697:   PetscFunctionBegin;
    698:   PetscCall(SNESReset_VI(snes));
    699:   PetscCall(ISDestroy(&vi->IS_inact_prev));
    700:   PetscFunctionReturn(PETSC_SUCCESS);
    701: }

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

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

    710:    Level: beginner

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

    718:    See {cite}`benson2006flexible`

    720: .seealso: [](ch_snes), `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
    721: M*/
    722: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
    723: {
    724:   SNES_VINEWTONRSLS *vi;
    725:   SNESLineSearch     linesearch;

    727:   PetscFunctionBegin;
    728:   snes->ops->reset          = SNESReset_VINEWTONRSLS;
    729:   snes->ops->setup          = SNESSetUp_VINEWTONRSLS;
    730:   snes->ops->solve          = SNESSolve_VINEWTONRSLS;
    731:   snes->ops->destroy        = SNESDestroy_VI;
    732:   snes->ops->setfromoptions = SNESSetFromOptions_VI;
    733:   snes->ops->view           = NULL;
    734:   snes->ops->converged      = SNESConvergedDefault_VI;

    736:   snes->usesksp = PETSC_TRUE;
    737:   snes->usesnpc = PETSC_FALSE;

    739:   PetscCall(SNESGetLineSearch(snes, &linesearch));
    740:   if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
    741:   PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));

    743:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    745:   PetscCall(SNESParametersInitialize(snes));

    747:   PetscCall(PetscNew(&vi));
    748:   snes->data          = (void *)vi;
    749:   vi->checkredundancy = NULL;

    751:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
    752:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
    753:   PetscFunctionReturn(PETSC_SUCCESS);
    754: }