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:   Note:
     18:   See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets

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

     26:   PetscFunctionBegin;
     27:   *inact = vi->IS_inact;
     28:   PetscFunctionReturn(PETSC_SUCCESS);
     29: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    122:   PetscCall(DMCreateGlobalVector(dm1, &finemarked));
    123:   PetscCall(DMCreateGlobalVector(*dm2, &coarsemarked));

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

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

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

    156:   PetscCall(DMClearGlobalVectors(dm1));

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

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

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

    168: static PetscErrorCode DMDestroy_SNESVI(PetscCtxRt ctx)
    169: {
    170:   DM_SNESVI *dmsnesvi = *(DM_SNESVI **)ctx;

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

    183:   PetscCall(ISDestroy(&dmsnesvi->inactive));
    184:   PetscCall(PetscFree(dmsnesvi));
    185:   PetscFunctionReturn(PETSC_SUCCESS);
    186: }

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

    192:   Logically Collective

    194:   Input Parameters:
    195: + dm       - the `DM` object
    196: - inactive - an `IS` indicating which points are currently not active

    198:   Level: intermediate

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

    207:   PetscFunctionBegin;
    208:   if (!dm) PetscFunctionReturn(PETSC_SUCCESS);

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

    212:   PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
    213:   if (!isnes) {
    214:     PetscCall(PetscNew(&dmsnesvi));
    215:     PetscCall(PetscObjectContainerCompose((PetscObject)dm, "VI", dmsnesvi, DMDestroy_SNESVI));

    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, &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: /* Create active and inactive set vectors. The local size of this vector is set and PETSc computes the global size */
    252: static PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
    253: {
    254:   Vec v;

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

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

    269:   PetscFunctionBegin;
    270:   PetscCall(SNESGetKSP(snes, &snesksp));
    271:   PetscCall(KSPReset(snesksp));
    272:   PetscCall(KSPResetFromOptions(snesksp));

    274:   /*
    275:   KSP                    kspnew;
    276:   PC                     pcnew;
    277:   MatSolverType          stype;

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

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

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

    316:   snes->numFailures            = 0;
    317:   snes->numLinearSolveFailures = 0;
    318:   snes->reason                 = SNES_CONVERGED_ITERATING;

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

    325:   PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm, SNESVIComputeInactiveSetFtY));
    326:   PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL));
    327:   PetscCall(SNESLineSearchSetUp(snes->linesearch));

    329:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    330:   snes->iter = 0;
    331:   snes->norm = 0.0;
    332:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

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

    344:   /* test convergence */
    345:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
    346:   PetscCall(SNESMonitor(snes, 0, fnorm));
    347:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

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

    358:     /* Call general purpose update function */
    359:     PetscTryTypeMethod(snes, update, snes->iter);
    360:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
    361:     SNESCheckJacobianDomainError(snes);

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

    365:     /*original
    366:     PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact));
    367:      */
    368:     PetscCall(SNESVIGetActiveSetIS(snes, X, F, &IS_act));

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

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

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

    394:         PetscCall(MatGetOwnershipRange(jac_inact_inact, &rstart, NULL));
    395:         PetscCall(MatDestroy(&jac_inact_inact));
    396:         PetscCall(ISDestroy(&IS_act));

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

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

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

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

    428:     /* Create active and inactive set vectors */
    429:     PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact));
    430:     PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act));
    431:     PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact));

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

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

    441:     PetscCall(VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
    442:     PetscCall(VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));

    444:     PetscCall(VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
    445:     PetscCall(VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));

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

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

    458:     /*      PetscCall(ISView(vi->IS_inact,0)); */
    459:     /*      PetscCall(ISView(IS_act,0));*/
    460:     /*      ierr = MatView(snes->jacobian_pre,0); */

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

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

    488:           PetscCall(PCBJacobiSetTotalBlocks(pc, 3, cnts));
    489:         }
    490:       }
    491:     }

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

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

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

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

    531:     if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
    532:     */
    533:     /* Compute a (scaled) negative update in the line search routine:
    534:          Y <- X - lambda*Y
    535:        and evaluate G = function(Y) (depends on the line search).
    536:     */
    537:     PetscCall(VecCopy(Y, snes->vec_sol_update));
    538:     ynorm = 1;
    539:     gnorm = fnorm;
    540:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y));
    541:     PetscCall(DMDestroyVI(snes->dm));
    542:     if (snes->reason) break;
    543:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsreason));
    544:     if (lsreason) {
    545:       if (snes->stol * xnorm > ynorm) {
    546:         snes->reason = SNES_CONVERGED_SNORM_RELATIVE;
    547:         break;
    548:       } else if (lsreason == SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN) {
    549:         snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
    550:         break;
    551:       } else if (lsreason == SNES_LINESEARCH_FAILED_NANORINF) {
    552:         snes->reason = SNES_DIVERGED_FUNCTION_NANORINF;
    553:         break;
    554:       } else if (lsreason == SNES_LINESEARCH_FAILED_OBJECTIVE_DOMAIN) {
    555:         snes->reason = SNES_DIVERGED_OBJECTIVE_DOMAIN;
    556:         break;
    557:       } else if (lsreason == SNES_LINESEARCH_FAILED_JACOBIAN_DOMAIN) {
    558:         snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN;
    559:         break;
    560:       } else if (++snes->numFailures >= snes->maxFailures) {
    561:         PetscBool ismin;

    563:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
    564:         PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin));
    565:         if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
    566:         break;
    567:       }
    568:     }
    569:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
    570:     PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lsreason));
    571:     /* Update function and solution vectors */
    572:     fnorm = gnorm;
    573:     /* Monitor convergence */
    574:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    575:     snes->iter  = i + 1;
    576:     snes->norm  = fnorm;
    577:     snes->xnorm = xnorm;
    578:     snes->ynorm = ynorm;
    579:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    580:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
    581:     /* Test for convergence, xnorm = || X || */
    582:     if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
    583:     PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
    584:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
    585:     if (snes->reason) break;
    586:   }
    587:   /* 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 */
    588:   PetscCall(DMDestroyVI(snes->dm));
    589:   PetscFunctionReturn(PETSC_SUCCESS);
    590: }

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

    595:   Logically Collective

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

    602:   Calling sequence of func:
    603: + snes      - the `SNES` context
    604: . is_act    - the set of points in the active sets
    605: . is_redact - output, the set of points in the non-redundant active set
    606: - ctx       - optional context

    608:   Level: advanced

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

    614:   See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets

    616: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
    617:  @*/
    618: PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES snes, IS is_act, IS *is_redact, PetscCtx ctx), PetscCtx ctx)
    619: {
    620:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

    622:   PetscFunctionBegin;
    624:   vi->checkredundancy = func;
    625:   vi->ctxP            = ctx;
    626:   PetscFunctionReturn(PETSC_SUCCESS);
    627: }

    629: #if defined(PETSC_HAVE_MATLAB)
    630:   #include <engine.h>
    631:   #include <mex.h>
    632: typedef struct {
    633:   char    *funcname;
    634:   mxArray *ctx;
    635: } SNESMatlabContext;

    637: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, PetscCtx ctx)
    638: {
    639:   SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
    640:   int                nlhs = 1, nrhs = 5;
    641:   mxArray           *plhs[1], *prhs[5];
    642:   long long int      l1 = 0, l2 = 0, ls = 0;
    643:   PetscInt          *indices = NULL;

    645:   PetscFunctionBegin;
    648:   PetscAssertPointer(is_redact, 3);
    649:   PetscCheckSameComm(snes, 1, is_act, 2);

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

    673: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
    674: {
    675:   SNESMatlabContext *sctx;

    677:   PetscFunctionBegin;
    678:   /* currently sctx is memory bleed */
    679:   PetscCall(PetscNew(&sctx));
    680:   PetscCall(PetscStrallocpy(func, &sctx->funcname));
    681:   sctx->ctx = mxDuplicateArray(ctx);
    682:   PetscCall(SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx));
    683:   PetscFunctionReturn(PETSC_SUCCESS);
    684: }

    686: #endif

    688: static PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
    689: {
    690:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
    691:   PetscInt          *indices;
    692:   PetscInt           i, n, rstart, rend;
    693:   SNESLineSearch     linesearch;

    695:   PetscFunctionBegin;
    696:   PetscCall(SNESSetUp_VI(snes));

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

    701:   PetscCall(VecGetOwnershipRange(snes->work[0], &rstart, &rend));
    702:   PetscCall(VecGetLocalSize(snes->work[0], &n));
    703:   PetscCall(PetscMalloc1(n, &indices));
    704:   for (i = 0; i < n; i++) indices[i] = rstart + i;
    705:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev));

    707:   /* set the line search functions */
    708:   if (!snes->linesearch) {
    709:     PetscCall(SNESGetLineSearch(snes, &linesearch));
    710:     PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
    711:   }
    712:   PetscFunctionReturn(PETSC_SUCCESS);
    713: }

    715: static PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
    716: {
    717:   SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;

    719:   PetscFunctionBegin;
    720:   PetscCall(SNESReset_VI(snes));
    721:   PetscCall(ISDestroy(&vi->IS_inact_prev));
    722:   PetscFunctionReturn(PETSC_SUCCESS);
    723: }

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

    728:    Options Database Keys:
    729: +  -snes_type (vinewtonssls|vinewtonrsls) - A semi-smooth solver or a reduced space active set method
    730: .  -snes_vi_zero_tolerance                - Tolerance for considering $u_i$ value to be on a bound.
    731: .  -snes_vi_monitor                       - Prints the number of active constraints (inactive set points) at each iteration.
    732: .  -snes_vi_monitor_residual              - View the residual vector at each iteration, using zero for active constraints (i.e. the inactive variables).
    733: -  -snes_vi_monitor_active                - View the active set by outputting a one for vector components in the active set and zero for the inactive.

    735:    Level: beginner

    737:    Note:
    738:    Reduced-space (active set methods) work as follows at each iteration\:
    739:    - The algorithm produces an inactive set of variables, that is a list of variables whose values will not be changed in the current iteration, i.e. they
    740:      are to be constrained to their current values. These are all the variables that are on the lower bound, that is $u_i = L_i$ with also $[F(u)]_i \ge 0$ or
    741:      the upper bound $u_i = U_i$ with also $[F(u)]_i \le 0.$
    742:    - A step direction is obtained by solving the linear system arising from the Jacobian used in Newton's method but with the inactive variables removed
    743:      from both the rows and columns.
    744:    - A line search is then used to update the active variables (the inactive set of variables are not changed).

    746:    The inactive set is chosen based on the sign of $[F(u)]_i$ because this gives exactly the set of points that would be moved outside of the domain given
    747:    an infinitesimal Newton (or even Richardson) step and our goal is to remain within the bounds, that is, to continue to satisfy the inequality constraints.

    749:    See {cite}`benson2006flexible`

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

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

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

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

    774:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    776:   PetscCall(SNESParametersInitialize(snes));

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

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