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: for (k = 0; k < n; k++) PetscCall(VecSetValue(finemarked, index[k], 1.0, INSERT_VALUES));
131: PetscCall(VecAssemblyBegin(finemarked));
132: PetscCall(VecAssemblyEnd(finemarked));
134: PetscCall(DMCreateInjection(*dm2, dm1, &inject));
135: PetscCall(MatRestrict(inject, finemarked, coarsemarked));
136: PetscCall(MatDestroy(&inject));
138: /*
139: create index set list of coarse inactive points from coarsemarked
140: */
141: PetscCall(VecGetLocalSize(coarsemarked, &n));
142: PetscCall(VecGetOwnershipRange(coarsemarked, &rstart, NULL));
143: PetscCall(VecGetArray(coarsemarked, &marked));
144: for (k = 0; k < n; k++) {
145: if (marked[k] != 0.0) cnt++;
146: }
147: PetscCall(PetscMalloc1(cnt, &coarseindex));
148: cnt = 0;
149: for (k = 0; k < n; k++) {
150: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
151: }
152: PetscCall(VecRestoreArray(coarsemarked, &marked));
153: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked), cnt, coarseindex, PETSC_OWN_POINTER, &inactive));
155: PetscCall(DMClearGlobalVectors(dm1));
157: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
159: PetscCall(DMSetVI(*dm2, inactive));
161: PetscCall(VecDestroy(&finemarked));
162: PetscCall(VecDestroy(&coarsemarked));
163: PetscCall(ISDestroy(&inactive));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: static PetscErrorCode DMDestroy_SNESVI(PetscCtxRt ctx)
168: {
169: DM_SNESVI *dmsnesvi = *(DM_SNESVI **)ctx;
171: PetscFunctionBegin;
172: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
173: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
174: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
175: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
176: dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection;
177: dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
178: /* need to clear out this vectors because some of them may not have a reference to the DM
179: but they are counted as having references to the DM in DMDestroy() */
180: PetscCall(DMClearGlobalVectors(dmsnesvi->dm));
182: PetscCall(ISDestroy(&dmsnesvi->inactive));
183: PetscCall(PetscFree(dmsnesvi));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: /*@
188: DMSetVI - Marks a `DM` as associated with a VI problem. This causes the interpolation/restriction operators to
189: be restricted to only those variables NOT associated with active constraints.
191: Logically Collective
193: Input Parameters:
194: + dm - the `DM` object
195: - inactive - an `IS` indicating which points are currently not active
197: Level: intermediate
199: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`
200: @*/
201: PetscErrorCode DMSetVI(DM dm, IS inactive)
202: {
203: PetscContainer isnes;
204: DM_SNESVI *dmsnesvi;
206: PetscFunctionBegin;
207: if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
209: PetscCall(PetscObjectReference((PetscObject)inactive));
211: PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
212: if (!isnes) {
213: PetscCall(PetscNew(&dmsnesvi));
214: PetscCall(PetscObjectContainerCompose((PetscObject)dm, "VI", dmsnesvi, DMDestroy_SNESVI));
216: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
217: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
218: dmsnesvi->coarsen = dm->ops->coarsen;
219: dm->ops->coarsen = DMCoarsen_SNESVI;
220: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
221: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
222: dmsnesvi->createinjection = dm->ops->createinjection;
223: dm->ops->createinjection = NULL;
224: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
225: dm->ops->hascreateinjection = NULL;
226: } else {
227: PetscCall(PetscContainerGetPointer(isnes, &dmsnesvi));
228: PetscCall(ISDestroy(&dmsnesvi->inactive));
229: }
230: PetscCall(DMClearGlobalVectors(dm));
231: PetscCall(ISGetLocalSize(inactive, &dmsnesvi->n));
233: dmsnesvi->inactive = inactive;
234: dmsnesvi->dm = dm;
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: /*@
239: DMDestroyVI - Frees the `DM_SNESVI` object contained in the `DM` and resets any function pointers the reduced-space `SNESVI` code composed onto it
241: Not Collective
243: Input Parameter:
244: . dm - the `DM` from which the VI context should be removed (may be `NULL`)
246: Level: developer
248: .seealso: `DM`, `SNESVINEWTONRSLS`, `SNESVISetVariableBounds()`, `PetscObjectCompose()`
249: @*/
250: PetscErrorCode DMDestroyVI(DM dm)
251: {
252: PetscFunctionBegin;
253: if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
254: PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)NULL));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /* Create active and inactive set vectors. The local size of this vector is set and PETSc computes the global size */
259: static PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
260: {
261: Vec v;
263: PetscFunctionBegin;
264: PetscCall(VecCreate(PetscObjectComm((PetscObject)snes), &v));
265: PetscCall(VecSetSizes(v, n, PETSC_DECIDE));
266: PetscCall(VecSetType(v, VECSTANDARD));
267: *newv = v;
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /* Resets the snes PC and KSP when the active set sizes change */
272: static PetscErrorCode SNESVIResetPCandKSP(SNES snes, Mat Amat, Mat Pmat)
273: {
274: KSP snesksp;
276: PetscFunctionBegin;
277: PetscCall(SNESGetKSP(snes, &snesksp));
278: PetscCall(KSPReset(snesksp));
279: PetscCall(KSPResetFromOptions(snesksp));
281: /*
282: KSP kspnew;
283: PC pcnew;
284: MatSolverType stype;
286: PetscCall(KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew));
287: kspnew->pc_side = snesksp->pc_side;
288: kspnew->rtol = snesksp->rtol;
289: kspnew->abstol = snesksp->abstol;
290: kspnew->max_it = snesksp->max_it;
291: PetscCall(KSPSetType(kspnew,((PetscObject)snesksp)->type_name));
292: PetscCall(KSPGetPC(kspnew,&pcnew));
293: PetscCall(PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name));
294: PetscCall(PCSetOperators(kspnew->pc,Amat,Pmat));
295: PetscCall(PCFactorGetMatSolverType(snesksp->pc,&stype));
296: PetscCall(PCFactorSetMatSolverType(kspnew->pc,stype));
297: PetscCall(KSPDestroy(&snesksp));
298: snes->ksp = kspnew;
299: PetscCall(KSPSetFromOptions(kspnew));*/
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
304: implemented in this algorithm. It basically identifies the active constraints and does
305: a linear solve on the other variables (those not associated with the active constraints). */
306: static PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
307: {
308: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
309: PetscInt maxits, i, lits;
310: SNESLineSearchReason lsreason;
311: PetscReal fnorm, gnorm, xnorm = 0, ynorm;
312: Vec Y, X, F;
313: KSPConvergedReason kspreason;
314: KSP ksp;
315: PC pc;
316: PetscBool isnleqerr;
318: PetscFunctionBegin;
319: /* SNESLINESEARCHNLEQERR solves with the KSP on the full space, but the KSP below is set up for the reduced (inactive set) space */
320: PetscCall(PetscObjectTypeCompare((PetscObject)snes->linesearch, SNESLINESEARCHNLEQERR, &isnleqerr));
321: PetscCheck(!isnleqerr, PetscObjectComm((PetscObject)snes), PETSC_ERR_SUP, "SNESLINESEARCHNLEQERR cannot be used with SNESVINEWTONRSLS since the line search solves with the KSP on the full space while SNESVINEWTONRSLS uses the KSP on the reduced (inactive set) space; use another line search, for example SNESLINESEARCHBT, or use SNESVINEWTONSSLS");
323: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
324: PetscCall(SNESGetKSP(snes, &ksp));
325: PetscCall(KSPGetPC(ksp, &pc));
326: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
328: snes->numFailures = 0;
329: snes->numLinearSolveFailures = 0;
330: snes->reason = SNES_CONVERGED_ITERATING;
332: maxits = snes->max_its; /* maximum number of iterations */
333: X = snes->vec_sol; /* solution vector */
334: F = snes->vec_func; /* residual vector */
335: Y = snes->work[0]; /* work vectors */
337: PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm, SNESVIComputeInactiveSetFtY));
338: PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL));
339: PetscCall(SNESLineSearchSetUp(snes->linesearch));
341: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
342: snes->iter = 0;
343: snes->norm = 0.0;
344: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
346: PetscCall(SNESVIProjectOntoBounds(snes, X));
347: PetscCall(SNESComputeFunction(snes, X, F));
348: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
349: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- ||x|| */
350: SNESCheckFunctionDomainError(snes, fnorm);
351: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
352: snes->norm = fnorm;
353: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
354: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
356: /* test convergence */
357: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
358: PetscCall(SNESMonitor(snes, 0, fnorm));
359: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
361: for (i = 0; i < maxits; i++) {
362: IS IS_act; /* _act -> active set _inact -> inactive set */
363: IS IS_redact; /* redundant active set */
364: VecScatter scat_act, scat_inact;
365: PetscInt nis_act, nis_inact;
366: Vec Y_act, Y_inact, F_inact;
367: Mat jac_inact_inact, prejac_inact_inact;
368: PetscBool isequal;
370: /* Call general purpose update function */
371: PetscTryTypeMethod(snes, update, snes->iter);
372: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
373: SNESCheckJacobianDomainError(snes);
375: /* Create active and inactive index sets */
377: /*original
378: PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact));
379: */
380: PetscCall(SNESVIGetActiveSetIS(snes, X, F, &IS_act));
382: if (vi->checkredundancy) {
383: PetscCall((*vi->checkredundancy)(snes, IS_act, &IS_redact, vi->ctxP));
384: if (IS_redact) {
385: PetscCall(ISSort(IS_redact));
386: PetscCall(ISComplement(IS_redact, X->map->rstart, X->map->rend, &vi->IS_inact));
387: PetscCall(ISDestroy(&IS_redact));
388: } else {
389: PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
390: }
391: } else {
392: PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
393: }
395: /* Create inactive set submatrix */
396: PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
398: if (0) { /* Dead code (temporary developer hack) */
399: IS keptrows;
400: PetscCall(MatFindNonzeroRows(jac_inact_inact, &keptrows));
401: if (keptrows) {
402: PetscInt cnt, *nrows, k;
403: const PetscInt *krows, *inact;
404: PetscInt rstart;
406: PetscCall(MatGetOwnershipRange(jac_inact_inact, &rstart, NULL));
407: PetscCall(MatDestroy(&jac_inact_inact));
408: PetscCall(ISDestroy(&IS_act));
410: PetscCall(ISGetLocalSize(keptrows, &cnt));
411: PetscCall(ISGetIndices(keptrows, &krows));
412: PetscCall(ISGetIndices(vi->IS_inact, &inact));
413: PetscCall(PetscMalloc1(cnt, &nrows));
414: for (k = 0; k < cnt; k++) nrows[k] = inact[krows[k] - rstart];
415: PetscCall(ISRestoreIndices(keptrows, &krows));
416: PetscCall(ISRestoreIndices(vi->IS_inact, &inact));
417: PetscCall(ISDestroy(&keptrows));
418: PetscCall(ISDestroy(&vi->IS_inact));
420: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), cnt, nrows, PETSC_OWN_POINTER, &vi->IS_inact));
421: PetscCall(ISComplement(vi->IS_inact, F->map->rstart, F->map->rend, &IS_act));
422: PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
423: }
424: }
425: PetscCall(DMSetVI(snes->dm, vi->IS_inact));
426: /* remove later */
428: /*
429: PetscCall(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)vi->xu)->comm)));
430: PetscCall(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)vi->xl)->comm)));
431: PetscCall(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X))));
432: PetscCall(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F))));
433: PetscCall(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact))));
434: */
436: /* Get sizes of active and inactive sets */
437: PetscCall(ISGetLocalSize(IS_act, &nis_act));
438: PetscCall(ISGetLocalSize(vi->IS_inact, &nis_inact));
440: /* Create active and inactive set vectors */
441: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact));
442: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act));
443: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact));
445: /* Create scatter contexts */
446: PetscCall(VecScatterCreate(Y, IS_act, Y_act, NULL, &scat_act));
447: PetscCall(VecScatterCreate(Y, vi->IS_inact, Y_inact, NULL, &scat_inact));
449: /* Do a vec scatter to active and inactive set vectors */
450: PetscCall(VecScatterBegin(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
451: PetscCall(VecScatterEnd(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
453: PetscCall(VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
454: PetscCall(VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
456: PetscCall(VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
457: PetscCall(VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
459: /* Active set direction = 0 */
460: PetscCall(VecSet(Y_act, 0));
461: if (snes->jacobian != snes->jacobian_pre) PetscCall(MatCreateSubMatrix(snes->jacobian_pre, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &prejac_inact_inact));
462: else prejac_inact_inact = jac_inact_inact;
464: PetscCall(ISEqual(vi->IS_inact_prev, vi->IS_inact, &isequal));
465: if (!isequal) {
466: PetscCall(SNESVIResetPCandKSP(snes, jac_inact_inact, prejac_inact_inact));
467: PetscCall(PCFieldSplitRestrictIS(pc, vi->IS_inact));
468: }
470: /* PetscCall(ISView(vi->IS_inact,0)); */
471: /* PetscCall(ISView(IS_act,0));*/
472: /* ierr = MatView(snes->jacobian_pre,0); */
474: PetscCall(KSPSetOperators(snes->ksp, jac_inact_inact, prejac_inact_inact));
475: PetscCall(KSPSetUp(snes->ksp));
476: {
477: PC pc;
478: PetscBool flg;
479: PetscCall(KSPGetPC(snes->ksp, &pc));
480: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg));
481: if (flg) {
482: KSP *subksps;
483: PetscCall(PCFieldSplitGetSubKSP(pc, NULL, &subksps));
484: PetscCall(KSPGetPC(subksps[0], &pc));
485: PetscCall(PetscFree(subksps));
486: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &flg));
487: if (flg) {
488: PetscInt n, N = 101 * 101, j, cnts[3] = {0, 0, 0};
489: const PetscInt *ii;
491: PetscCall(ISGetSize(vi->IS_inact, &n));
492: PetscCall(ISGetIndices(vi->IS_inact, &ii));
493: for (j = 0; j < n; j++) {
494: if (ii[j] < N) cnts[0]++;
495: else if (ii[j] < 2 * N) cnts[1]++;
496: else if (ii[j] < 3 * N) cnts[2]++;
497: }
498: PetscCall(ISRestoreIndices(vi->IS_inact, &ii));
500: PetscCall(PCBJacobiSetTotalBlocks(pc, 3, cnts));
501: }
502: }
503: }
505: PetscCall(KSPSolve(snes->ksp, F_inact, Y_inact));
506: PetscCall(VecScatterBegin(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
507: PetscCall(VecScatterEnd(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
508: PetscCall(VecScatterBegin(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
509: PetscCall(VecScatterEnd(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
511: PetscCall(VecDestroy(&F_inact));
512: PetscCall(VecDestroy(&Y_act));
513: PetscCall(VecDestroy(&Y_inact));
514: PetscCall(VecScatterDestroy(&scat_act));
515: PetscCall(VecScatterDestroy(&scat_inact));
516: PetscCall(ISDestroy(&IS_act));
517: if (!isequal) {
518: PetscCall(ISDestroy(&vi->IS_inact_prev));
519: PetscCall(ISDuplicate(vi->IS_inact, &vi->IS_inact_prev));
520: }
521: PetscCall(ISDestroy(&vi->IS_inact));
522: PetscCall(MatDestroy(&jac_inact_inact));
523: if (snes->jacobian != snes->jacobian_pre) PetscCall(MatDestroy(&prejac_inact_inact));
525: PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
526: if (kspreason < 0) {
527: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
528: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures));
529: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
530: break;
531: }
532: }
534: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
535: snes->linear_its += lits;
536: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
537: /*
538: if (snes->ops->precheck) {
539: PetscBool changed_y = PETSC_FALSE;
540: PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
541: }
543: if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
544: */
545: /* Compute a (scaled) negative update in the line search routine:
546: Y <- X - lambda*Y
547: and evaluate G = function(Y) (depends on the line search).
548: */
549: PetscCall(VecCopy(Y, snes->vec_sol_update));
550: ynorm = 1;
551: gnorm = fnorm;
552: PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y));
553: PetscCall(DMDestroyVI(snes->dm));
554: if (snes->reason) break;
555: PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsreason));
556: if (lsreason) {
557: if (snes->stol * xnorm > ynorm) {
558: snes->reason = SNES_CONVERGED_SNORM_RELATIVE;
559: break;
560: } else if (lsreason == SNES_LINESEARCH_FAILED_FUNCTION_DOMAIN) {
561: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
562: break;
563: } else if (lsreason == SNES_LINESEARCH_FAILED_NANORINF) {
564: snes->reason = SNES_DIVERGED_FUNCTION_NANORINF;
565: break;
566: } else if (lsreason == SNES_LINESEARCH_FAILED_OBJECTIVE_DOMAIN) {
567: snes->reason = SNES_DIVERGED_OBJECTIVE_DOMAIN;
568: break;
569: } else if (lsreason == SNES_LINESEARCH_FAILED_JACOBIAN_DOMAIN) {
570: snes->reason = SNES_DIVERGED_JACOBIAN_DOMAIN;
571: break;
572: } else if (++snes->numFailures >= snes->maxFailures) {
573: PetscBool ismin;
575: snes->reason = SNES_DIVERGED_LINE_SEARCH;
576: PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin));
577: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
578: break;
579: }
580: }
581: PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
582: PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lsreason));
583: /* Update function and solution vectors */
584: fnorm = gnorm;
585: /* Monitor convergence */
586: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
587: snes->iter = i + 1;
588: snes->norm = fnorm;
589: snes->xnorm = xnorm;
590: snes->ynorm = ynorm;
591: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
592: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
593: /* Test for convergence, xnorm = || X || */
594: if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
595: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
596: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
597: if (snes->reason) break;
598: }
599: /* 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 */
600: PetscCall(DMDestroyVI(snes->dm));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*@
605: SNESVISetRedundancyCheck - Provide a function to check for any redundancy in the VI active set
607: Logically Collective
609: Input Parameters:
610: + snes - the `SNESVINEWTONRSLS` context
611: . func - the function to check of redundancies
612: - ctx - optional context used by the function
614: Calling sequence of func:
615: + snes - the `SNES` context
616: . is_act - the set of points in the active sets
617: . is_redact - output, the set of points in the non-redundant active set
618: - ctx - optional context
620: Level: advanced
622: Note:
623: Sometimes the inactive set will result in a singular sub-Jacobian problem that needs to be solved, this allows the user,
624: when they know more about their specific problem to provide a function that removes the redundancy that results in the singular linear system
626: See `SNESVINEWTONRSLS` for a concise description of the active and inactive sets
628: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
629: @*/
630: PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES snes, IS is_act, IS *is_redact, PetscCtx ctx), PetscCtx ctx)
631: {
632: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
634: PetscFunctionBegin;
636: vi->checkredundancy = func;
637: vi->ctxP = ctx;
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: #if PetscDefined(HAVE_MATLAB)
642: #include <engine.h>
643: #include <mex.h>
644: typedef struct {
645: char *funcname;
646: mxArray *ctx;
647: } SNESMatlabContext;
649: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, PetscCtx ctx)
650: {
651: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
652: int nlhs = 1, nrhs = 5;
653: mxArray *plhs[1], *prhs[5];
654: long long int l1 = 0, l2 = 0, ls = 0;
655: PetscInt *indices = NULL;
657: PetscFunctionBegin;
660: PetscAssertPointer(is_redact, 3);
661: PetscCheckSameComm(snes, 1, is_act, 2);
663: /* Create IS for reduced active set of size 0, its size and indices will
664: bet set by the MATLAB function */
665: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), 0, indices, PETSC_OWN_POINTER, is_redact));
666: /* call MATLAB function in ctx */
667: PetscCall(PetscArraycpy(&ls, &snes, 1));
668: PetscCall(PetscArraycpy(&l1, &is_act, 1));
669: PetscCall(PetscArraycpy(&l2, is_redact, 1));
670: prhs[0] = mxCreateDoubleScalar((double)ls);
671: prhs[1] = mxCreateDoubleScalar((double)l1);
672: prhs[2] = mxCreateDoubleScalar((double)l2);
673: prhs[3] = mxCreateString(sctx->funcname);
674: prhs[4] = sctx->ctx;
675: PetscCall(mexCallMATLAB(nlhs, plhs, nrhs, prhs, "PetscSNESVIRedundancyCheckInternal"));
676: PetscCall(mxGetScalar(plhs[0]));
677: mxDestroyArray(prhs[0]);
678: mxDestroyArray(prhs[1]);
679: mxDestroyArray(prhs[2]);
680: mxDestroyArray(prhs[3]);
681: mxDestroyArray(plhs[0]);
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
686: {
687: SNESMatlabContext *sctx;
689: PetscFunctionBegin;
690: /* currently sctx is memory bleed */
691: PetscCall(PetscNew(&sctx));
692: PetscCall(PetscStrallocpy(func, &sctx->funcname));
693: sctx->ctx = mxDuplicateArray(ctx);
694: PetscCall(SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx));
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: #endif
700: static PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
701: {
702: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
703: PetscInt *indices;
704: PetscInt i, n, rstart, rend;
705: SNESLineSearch linesearch;
707: PetscFunctionBegin;
708: PetscCall(SNESSetUp_VI(snes));
710: /* Set up previous active index set for the first snes solve
711: vi->IS_inact_prev = 0,1,2,....N */
713: PetscCall(VecGetOwnershipRange(snes->work[0], &rstart, &rend));
714: PetscCall(VecGetLocalSize(snes->work[0], &n));
715: PetscCall(PetscMalloc1(n, &indices));
716: for (i = 0; i < n; i++) indices[i] = rstart + i;
717: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev));
719: /* set the line search functions */
720: if (!snes->linesearch) {
721: PetscCall(SNESGetLineSearch(snes, &linesearch));
722: PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
723: }
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: static PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
728: {
729: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
731: PetscFunctionBegin;
732: PetscCall(SNESReset_VI(snes));
733: PetscCall(ISDestroy(&vi->IS_inact_prev));
734: PetscFunctionReturn(PETSC_SUCCESS);
735: }
737: /*MC
738: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
740: Options Database Keys:
741: + -snes_type (vinewtonssls|vinewtonrsls) - A semi-smooth solver or a reduced space active set method
742: . -snes_vi_zero_tolerance - Tolerance for considering $u_i$ value to be on a bound.
743: . -snes_vi_monitor - Prints the number of active constraints (inactive set points) at each iteration.
744: . -snes_vi_monitor_residual - View the residual vector at each iteration, using zero for active constraints (i.e. the inactive variables).
745: - -snes_vi_monitor_active - View the active set by outputting a one for vector components in the active set and zero for the inactive.
747: Level: beginner
749: Note:
750: Reduced-space (active set methods) work as follows at each iteration\:
751: - 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
752: 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
753: the upper bound $u_i = U_i$ with also $[F(u)]_i \le 0.$
754: - 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
755: from both the rows and columns.
756: - A line search is then used to update the active variables (the inactive set of variables are not changed).
758: 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
759: 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.
761: See {cite}`benson2006flexible`
763: .seealso: [](ch_snes), `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
764: M*/
765: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
766: {
767: SNES_VINEWTONRSLS *vi;
768: SNESLineSearch linesearch;
770: PetscFunctionBegin;
771: snes->ops->reset = SNESReset_VINEWTONRSLS;
772: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
773: snes->ops->solve = SNESSolve_VINEWTONRSLS;
774: snes->ops->destroy = SNESDestroy_VI;
775: snes->ops->setfromoptions = SNESSetFromOptions_VI;
776: snes->ops->view = NULL;
777: if (!snes->ops->converged || snes->ops->converged == SNESConvergedDefault) snes->ops->converged = SNESConvergedDefault_VI;
779: snes->usesksp = PETSC_TRUE;
780: snes->usesnpc = PETSC_FALSE;
782: PetscCall(SNESGetLineSearch(snes, &linesearch));
783: if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
784: PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));
786: snes->alwayscomputesfinalresidual = PETSC_TRUE;
788: PetscCall(PetscNew(&vi));
789: snes->data = (void *)vi;
791: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
792: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
793: PetscFunctionReturn(PETSC_SUCCESS);
794: }