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(PetscContainerCreate(PetscObjectComm((PetscObject)dm), &isnes));
212: PetscCall(PetscContainerSetUserDestroy(isnes, (PetscErrorCode(*)(void *))DMDestroy_SNESVI));
213: PetscCall(PetscNew(&dmsnesvi));
214: PetscCall(PetscContainerSetPointer(isnes, (void *)dmsnesvi));
215: PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)isnes));
216: PetscCall(PetscContainerDestroy(&isnes));
218: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
219: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
220: dmsnesvi->coarsen = dm->ops->coarsen;
221: dm->ops->coarsen = DMCoarsen_SNESVI;
222: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
223: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
224: dmsnesvi->createinjection = dm->ops->createinjection;
225: dm->ops->createinjection = NULL;
226: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
227: dm->ops->hascreateinjection = NULL;
228: } else {
229: PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi));
230: PetscCall(ISDestroy(&dmsnesvi->inactive));
231: }
232: PetscCall(DMClearGlobalVectors(dm));
233: PetscCall(ISGetLocalSize(inactive, &dmsnesvi->n));
235: dmsnesvi->inactive = inactive;
236: dmsnesvi->dm = dm;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*
241: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
242: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
243: */
244: PetscErrorCode DMDestroyVI(DM dm)
245: {
246: PetscFunctionBegin;
247: if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
248: PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)NULL));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
253: static PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
254: {
255: Vec v;
257: PetscFunctionBegin;
258: PetscCall(VecCreate(PetscObjectComm((PetscObject)snes), &v));
259: PetscCall(VecSetSizes(v, n, PETSC_DECIDE));
260: PetscCall(VecSetType(v, VECSTANDARD));
261: *newv = v;
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /* Resets the snes PC and KSP when the active set sizes change */
266: static PetscErrorCode SNESVIResetPCandKSP(SNES snes, Mat Amat, Mat Pmat)
267: {
268: KSP snesksp;
270: PetscFunctionBegin;
271: PetscCall(SNESGetKSP(snes, &snesksp));
272: PetscCall(KSPReset(snesksp));
273: PetscCall(KSPResetFromOptions(snesksp));
275: /*
276: KSP kspnew;
277: PC pcnew;
278: MatSolverType stype;
280: PetscCall(KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew));
281: kspnew->pc_side = snesksp->pc_side;
282: kspnew->rtol = snesksp->rtol;
283: kspnew->abstol = snesksp->abstol;
284: kspnew->max_it = snesksp->max_it;
285: PetscCall(KSPSetType(kspnew,((PetscObject)snesksp)->type_name));
286: PetscCall(KSPGetPC(kspnew,&pcnew));
287: PetscCall(PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name));
288: PetscCall(PCSetOperators(kspnew->pc,Amat,Pmat));
289: PetscCall(PCFactorGetMatSolverType(snesksp->pc,&stype));
290: PetscCall(PCFactorSetMatSolverType(kspnew->pc,stype));
291: PetscCall(KSPDestroy(&snesksp));
292: snes->ksp = kspnew;
293: PetscCall(KSPSetFromOptions(kspnew));*/
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
298: implemented in this algorithm. It basically identifies the active constraints and does
299: a linear solve on the other variables (those not associated with the active constraints). */
300: static PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
301: {
302: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
303: PetscInt maxits, i, lits;
304: SNESLineSearchReason lssucceed;
305: PetscReal fnorm, gnorm, xnorm = 0, ynorm;
306: Vec Y, X, F;
307: KSPConvergedReason kspreason;
308: KSP ksp;
309: PC pc;
311: PetscFunctionBegin;
312: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
313: PetscCall(SNESGetKSP(snes, &ksp));
314: PetscCall(KSPGetPC(ksp, &pc));
315: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
317: snes->numFailures = 0;
318: snes->numLinearSolveFailures = 0;
319: snes->reason = SNES_CONVERGED_ITERATING;
321: maxits = snes->max_its; /* maximum number of iterations */
322: X = snes->vec_sol; /* solution vector */
323: F = snes->vec_func; /* residual vector */
324: Y = snes->work[0]; /* work vectors */
326: PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm));
327: PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL));
328: PetscCall(SNESLineSearchSetUp(snes->linesearch));
330: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
331: snes->iter = 0;
332: snes->norm = 0.0;
333: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
335: PetscCall(SNESVIProjectOntoBounds(snes, X));
336: PetscCall(SNESComputeFunction(snes, X, F));
337: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
338: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- ||x|| */
339: SNESCheckFunctionNorm(snes, fnorm);
340: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
341: snes->norm = fnorm;
342: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
343: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
345: /* test convergence */
346: PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
347: PetscCall(SNESMonitor(snes, 0, fnorm));
348: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
350: for (i = 0; i < maxits; i++) {
351: IS IS_act; /* _act -> active set _inact -> inactive set */
352: IS IS_redact; /* redundant active set */
353: VecScatter scat_act, scat_inact;
354: PetscInt nis_act, nis_inact;
355: Vec Y_act, Y_inact, F_inact;
356: Mat jac_inact_inact, prejac_inact_inact;
357: PetscBool isequal;
359: /* Call general purpose update function */
360: PetscTryTypeMethod(snes, update, snes->iter);
361: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
362: SNESCheckJacobianDomainerror(snes);
364: /* Create active and inactive index sets */
366: /*original
367: PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact));
368: */
369: PetscCall(SNESVIGetActiveSetIS(snes, X, F, &IS_act));
371: if (vi->checkredundancy) {
372: PetscCall((*vi->checkredundancy)(snes, IS_act, &IS_redact, vi->ctxP));
373: if (IS_redact) {
374: PetscCall(ISSort(IS_redact));
375: PetscCall(ISComplement(IS_redact, X->map->rstart, X->map->rend, &vi->IS_inact));
376: PetscCall(ISDestroy(&IS_redact));
377: } else {
378: PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
379: }
380: } else {
381: PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
382: }
384: /* Create inactive set submatrix */
385: PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
387: if (0) { /* Dead code (temporary developer hack) */
388: IS keptrows;
389: PetscCall(MatFindNonzeroRows(jac_inact_inact, &keptrows));
390: if (keptrows) {
391: PetscInt cnt, *nrows, k;
392: const PetscInt *krows, *inact;
393: PetscInt rstart;
395: PetscCall(MatGetOwnershipRange(jac_inact_inact, &rstart, NULL));
396: PetscCall(MatDestroy(&jac_inact_inact));
397: PetscCall(ISDestroy(&IS_act));
399: PetscCall(ISGetLocalSize(keptrows, &cnt));
400: PetscCall(ISGetIndices(keptrows, &krows));
401: PetscCall(ISGetIndices(vi->IS_inact, &inact));
402: PetscCall(PetscMalloc1(cnt, &nrows));
403: for (k = 0; k < cnt; k++) nrows[k] = inact[krows[k] - rstart];
404: PetscCall(ISRestoreIndices(keptrows, &krows));
405: PetscCall(ISRestoreIndices(vi->IS_inact, &inact));
406: PetscCall(ISDestroy(&keptrows));
407: PetscCall(ISDestroy(&vi->IS_inact));
409: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), cnt, nrows, PETSC_OWN_POINTER, &vi->IS_inact));
410: PetscCall(ISComplement(vi->IS_inact, F->map->rstart, F->map->rend, &IS_act));
411: PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
412: }
413: }
414: PetscCall(DMSetVI(snes->dm, vi->IS_inact));
415: /* remove later */
417: /*
418: PetscCall(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)vi->xu)->comm)));
419: PetscCall(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)vi->xl)->comm)));
420: PetscCall(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X))));
421: PetscCall(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F))));
422: PetscCall(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact))));
423: */
425: /* Get sizes of active and inactive sets */
426: PetscCall(ISGetLocalSize(IS_act, &nis_act));
427: PetscCall(ISGetLocalSize(vi->IS_inact, &nis_inact));
429: /* Create active and inactive set vectors */
430: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact));
431: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act));
432: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact));
434: /* Create scatter contexts */
435: PetscCall(VecScatterCreate(Y, IS_act, Y_act, NULL, &scat_act));
436: PetscCall(VecScatterCreate(Y, vi->IS_inact, Y_inact, NULL, &scat_inact));
438: /* Do a vec scatter to active and inactive set vectors */
439: PetscCall(VecScatterBegin(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
440: PetscCall(VecScatterEnd(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
442: PetscCall(VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
443: PetscCall(VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
445: PetscCall(VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
446: PetscCall(VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
448: /* Active set direction = 0 */
449: PetscCall(VecSet(Y_act, 0));
450: if (snes->jacobian != snes->jacobian_pre) {
451: PetscCall(MatCreateSubMatrix(snes->jacobian_pre, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &prejac_inact_inact));
452: } else prejac_inact_inact = jac_inact_inact;
454: PetscCall(ISEqual(vi->IS_inact_prev, vi->IS_inact, &isequal));
455: if (!isequal) {
456: PetscCall(SNESVIResetPCandKSP(snes, jac_inact_inact, prejac_inact_inact));
457: PetscCall(PCFieldSplitRestrictIS(pc, vi->IS_inact));
458: }
460: /* PetscCall(ISView(vi->IS_inact,0)); */
461: /* PetscCall(ISView(IS_act,0));*/
462: /* ierr = MatView(snes->jacobian_pre,0); */
464: PetscCall(KSPSetOperators(snes->ksp, jac_inact_inact, prejac_inact_inact));
465: PetscCall(KSPSetUp(snes->ksp));
466: {
467: PC pc;
468: PetscBool flg;
469: PetscCall(KSPGetPC(snes->ksp, &pc));
470: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg));
471: if (flg) {
472: KSP *subksps;
473: PetscCall(PCFieldSplitGetSubKSP(pc, NULL, &subksps));
474: PetscCall(KSPGetPC(subksps[0], &pc));
475: PetscCall(PetscFree(subksps));
476: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &flg));
477: if (flg) {
478: PetscInt n, N = 101 * 101, j, cnts[3] = {0, 0, 0};
479: const PetscInt *ii;
481: PetscCall(ISGetSize(vi->IS_inact, &n));
482: PetscCall(ISGetIndices(vi->IS_inact, &ii));
483: for (j = 0; j < n; j++) {
484: if (ii[j] < N) cnts[0]++;
485: else if (ii[j] < 2 * N) cnts[1]++;
486: else if (ii[j] < 3 * N) cnts[2]++;
487: }
488: PetscCall(ISRestoreIndices(vi->IS_inact, &ii));
490: PetscCall(PCBJacobiSetTotalBlocks(pc, 3, cnts));
491: }
492: }
493: }
495: PetscCall(KSPSolve(snes->ksp, F_inact, Y_inact));
496: PetscCall(VecScatterBegin(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
497: PetscCall(VecScatterEnd(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
498: PetscCall(VecScatterBegin(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
499: PetscCall(VecScatterEnd(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
501: PetscCall(VecDestroy(&F_inact));
502: PetscCall(VecDestroy(&Y_act));
503: PetscCall(VecDestroy(&Y_inact));
504: PetscCall(VecScatterDestroy(&scat_act));
505: PetscCall(VecScatterDestroy(&scat_inact));
506: PetscCall(ISDestroy(&IS_act));
507: if (!isequal) {
508: PetscCall(ISDestroy(&vi->IS_inact_prev));
509: PetscCall(ISDuplicate(vi->IS_inact, &vi->IS_inact_prev));
510: }
511: PetscCall(ISDestroy(&vi->IS_inact));
512: PetscCall(MatDestroy(&jac_inact_inact));
513: if (snes->jacobian != snes->jacobian_pre) PetscCall(MatDestroy(&prejac_inact_inact));
515: PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
516: if (kspreason < 0) {
517: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
518: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures));
519: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
520: break;
521: }
522: }
524: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
525: snes->linear_its += lits;
526: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
527: /*
528: if (snes->ops->precheck) {
529: PetscBool changed_y = PETSC_FALSE;
530: PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
531: }
533: if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
534: */
535: /* Compute a (scaled) negative update in the line search routine:
536: Y <- X - lambda*Y
537: and evaluate G = function(Y) (depends on the line search).
538: */
539: PetscCall(VecCopy(Y, snes->vec_sol_update));
540: ynorm = 1;
541: gnorm = fnorm;
542: PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y));
543: PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
544: PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
545: PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lssucceed));
546: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
547: if (snes->domainerror) {
548: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
549: PetscCall(DMDestroyVI(snes->dm));
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
552: if (lssucceed) {
553: if (++snes->numFailures >= snes->maxFailures) {
554: PetscBool ismin;
555: snes->reason = SNES_DIVERGED_LINE_SEARCH;
556: PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin));
557: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
558: break;
559: }
560: }
561: PetscCall(DMDestroyVI(snes->dm));
562: /* Update function and solution vectors */
563: fnorm = gnorm;
564: /* Monitor convergence */
565: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
566: snes->iter = i + 1;
567: snes->norm = fnorm;
568: snes->xnorm = xnorm;
569: snes->ynorm = ynorm;
570: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
571: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
572: /* Test for convergence, xnorm = || X || */
573: if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
574: PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
575: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
576: if (snes->reason) break;
577: }
578: /* 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 */
579: PetscCall(DMDestroyVI(snes->dm));
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: /*@C
584: SNESVISetRedundancyCheck - Provide a function to check for any redundancy in the VI active set
586: Logically Collective
588: Input Parameters:
589: + snes - the `SNESVINEWTONRSLS` context
590: . func - the function to check of redundancies
591: - ctx - optional context used by the function
593: Level: advanced
595: Note:
596: Sometimes the inactive set will result in a singular sub-Jacobian problem that needs to be solved, this allows the user,
597: when they know more about their specific problem to provide a function that removes the redundancy that results in the singular linear system
599: .seealso: [](ch_snes), `SNES`, `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
600: @*/
601: PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES, IS, IS *, void *), void *ctx)
602: {
603: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
605: PetscFunctionBegin;
607: vi->checkredundancy = func;
608: vi->ctxP = ctx;
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: #if defined(PETSC_HAVE_MATLAB)
613: #include <engine.h>
614: #include <mex.h>
615: typedef struct {
616: char *funcname;
617: mxArray *ctx;
618: } SNESMatlabContext;
620: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, void *ctx)
621: {
622: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
623: int nlhs = 1, nrhs = 5;
624: mxArray *plhs[1], *prhs[5];
625: long long int l1 = 0, l2 = 0, ls = 0;
626: PetscInt *indices = NULL;
628: PetscFunctionBegin;
631: PetscAssertPointer(is_redact, 3);
632: PetscCheckSameComm(snes, 1, is_act, 2);
634: /* Create IS for reduced active set of size 0, its size and indices will
635: bet set by the MATLAB function */
636: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), 0, indices, PETSC_OWN_POINTER, is_redact));
637: /* call MATLAB function in ctx */
638: PetscCall(PetscArraycpy(&ls, &snes, 1));
639: PetscCall(PetscArraycpy(&l1, &is_act, 1));
640: PetscCall(PetscArraycpy(&l2, is_redact, 1));
641: prhs[0] = mxCreateDoubleScalar((double)ls);
642: prhs[1] = mxCreateDoubleScalar((double)l1);
643: prhs[2] = mxCreateDoubleScalar((double)l2);
644: prhs[3] = mxCreateString(sctx->funcname);
645: prhs[4] = sctx->ctx;
646: PetscCall(mexCallMATLAB(nlhs, plhs, nrhs, prhs, "PetscSNESVIRedundancyCheckInternal"));
647: PetscCall(mxGetScalar(plhs[0]));
648: mxDestroyArray(prhs[0]);
649: mxDestroyArray(prhs[1]);
650: mxDestroyArray(prhs[2]);
651: mxDestroyArray(prhs[3]);
652: mxDestroyArray(plhs[0]);
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
657: {
658: SNESMatlabContext *sctx;
660: PetscFunctionBegin;
661: /* currently sctx is memory bleed */
662: PetscCall(PetscNew(&sctx));
663: PetscCall(PetscStrallocpy(func, &sctx->funcname));
664: sctx->ctx = mxDuplicateArray(ctx);
665: PetscCall(SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: #endif
671: static PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
672: {
673: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
674: PetscInt *indices;
675: PetscInt i, n, rstart, rend;
676: SNESLineSearch linesearch;
678: PetscFunctionBegin;
679: PetscCall(SNESSetUp_VI(snes));
681: /* Set up previous active index set for the first snes solve
682: vi->IS_inact_prev = 0,1,2,....N */
684: PetscCall(VecGetOwnershipRange(snes->vec_sol, &rstart, &rend));
685: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
686: PetscCall(PetscMalloc1(n, &indices));
687: for (i = 0; i < n; i++) indices[i] = rstart + i;
688: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev));
690: /* set the line search functions */
691: if (!snes->linesearch) {
692: PetscCall(SNESGetLineSearch(snes, &linesearch));
693: PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
694: }
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: static PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
699: {
700: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
702: PetscFunctionBegin;
703: PetscCall(SNESReset_VI(snes));
704: PetscCall(ISDestroy(&vi->IS_inact_prev));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: /*MC
709: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
711: Options Database Keys:
712: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver or a reduced space active set method
713: - -snes_vi_monitor - prints the number of active constraints at each iteration.
715: Level: beginner
717: Note:
718: At each set of this methods the algorithm produces an inactive set of variables that are constrained to their current values
719: (because changing these values would result in those variables no longer satisfying the inequality constraints)
720: and produces a step direction by solving the linear system arising from the Jacobian with the inactive variables removed. In other
721: words on a reduced space of the solution space. Based on the Newton update it then adjusts the inactive set for the next iteration.
723: See {cite}`benson2006flexible`
725: .seealso: [](ch_snes), `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
726: M*/
727: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
728: {
729: SNES_VINEWTONRSLS *vi;
730: SNESLineSearch linesearch;
732: PetscFunctionBegin;
733: snes->ops->reset = SNESReset_VINEWTONRSLS;
734: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
735: snes->ops->solve = SNESSolve_VINEWTONRSLS;
736: snes->ops->destroy = SNESDestroy_VI;
737: snes->ops->setfromoptions = SNESSetFromOptions_VI;
738: snes->ops->view = NULL;
739: snes->ops->converged = SNESConvergedDefault_VI;
741: snes->usesksp = PETSC_TRUE;
742: snes->usesnpc = PETSC_FALSE;
744: PetscCall(SNESGetLineSearch(snes, &linesearch));
745: if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
746: PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));
748: snes->alwayscomputesfinalresidual = PETSC_TRUE;
750: PetscCall(PetscNew(&vi));
751: snes->data = (void *)vi;
752: vi->checkredundancy = NULL;
754: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
755: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }