Actual source code: virs.c
2: #include <../src/snes/impls/vi/rs/virsimpl.h>
3: #include <petsc/private/dmimpl.h>
4: #include <petsc/private/vecimpl.h>
6: /*@
7: SNESVIGetInactiveSet - Gets the global indices for the inactive set variables (these correspond to the degrees of freedom the linear
8: system is solved on)
10: Input parameter:
11: . snes - the `SNES` context
13: Output parameter:
14: . inact - inactive set index set
16: Level: advanced
18: .seealso: `SNESVINEWTONRSLS`
19: @*/
20: PetscErrorCode SNESVIGetInactiveSet(SNES snes, IS *inact)
21: {
22: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
24: PetscFunctionBegin;
25: *inact = vi->IS_inact;
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: /*
30: Provides a wrapper to a DM to allow it to be used to generated the interpolation/restriction from the DM for the smaller matrices and vectors
31: defined by the reduced space method.
33: Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints.
34: */
35: typedef struct {
36: PetscInt n; /* size of vectors in the reduced DM space */
37: IS inactive;
39: PetscErrorCode (*createinterpolation)(DM, DM, Mat *, Vec *); /* DM's original routines */
40: PetscErrorCode (*coarsen)(DM, MPI_Comm, DM *);
41: PetscErrorCode (*createglobalvector)(DM, Vec *);
42: PetscErrorCode (*createinjection)(DM, DM, Mat *);
43: PetscErrorCode (*hascreateinjection)(DM, PetscBool *);
45: DM dm; /* when destroying this object we need to reset the above function into the base DM */
46: } DM_SNESVI;
48: /*
49: DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
50: */
51: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm, Vec *vec)
52: {
53: PetscContainer isnes;
54: DM_SNESVI *dmsnesvi;
56: PetscFunctionBegin;
57: PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
58: PetscCheck(isnes, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Composed SNES is missing");
59: PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi));
60: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)dm), dmsnesvi->n, PETSC_DETERMINE, vec));
61: PetscCall(VecSetDM(*vec, dm));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /*
66: DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
67: */
68: PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1, DM dm2, Mat *mat, Vec *vec)
69: {
70: PetscContainer isnes;
71: DM_SNESVI *dmsnesvi1, *dmsnesvi2;
72: Mat interp;
74: PetscFunctionBegin;
75: PetscCall(PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes));
76: PetscCheck(isnes, PetscObjectComm((PetscObject)dm1), PETSC_ERR_PLIB, "Composed VI data structure is missing");
77: PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi1));
78: PetscCall(PetscObjectQuery((PetscObject)dm2, "VI", (PetscObject *)&isnes));
79: PetscCheck(isnes, PetscObjectComm((PetscObject)dm2), PETSC_ERR_PLIB, "Composed VI data structure is missing");
80: PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi2));
82: PetscCall((*dmsnesvi1->createinterpolation)(dm1, dm2, &interp, NULL));
83: PetscCall(MatCreateSubMatrix(interp, dmsnesvi2->inactive, dmsnesvi1->inactive, MAT_INITIAL_MATRIX, mat));
84: PetscCall(MatDestroy(&interp));
85: *vec = NULL;
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /*
90: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
91: */
92: PetscErrorCode DMCoarsen_SNESVI(DM dm1, MPI_Comm comm, DM *dm2)
93: {
94: PetscContainer isnes;
95: DM_SNESVI *dmsnesvi1;
96: Vec finemarked, coarsemarked;
97: IS inactive;
98: Mat inject;
99: const PetscInt *index;
100: PetscInt n, k, cnt = 0, rstart, *coarseindex;
101: PetscScalar *marked;
103: PetscFunctionBegin;
104: PetscCall(PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes));
105: PetscCheck(isnes, PetscObjectComm((PetscObject)dm1), PETSC_ERR_PLIB, "Composed VI data structure is missing");
106: PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi1));
108: /* get the original coarsen */
109: PetscCall((*dmsnesvi1->coarsen)(dm1, comm, dm2));
111: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
112: /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
113: /* PetscCall(PetscObjectReference((PetscObject)*dm2));*/
115: /* need to set back global vectors in order to use the original injection */
116: PetscCall(DMClearGlobalVectors(dm1));
118: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
120: PetscCall(DMCreateGlobalVector(dm1, &finemarked));
121: PetscCall(DMCreateGlobalVector(*dm2, &coarsemarked));
123: /*
124: fill finemarked with locations of inactive points
125: */
126: PetscCall(ISGetIndices(dmsnesvi1->inactive, &index));
127: PetscCall(ISGetLocalSize(dmsnesvi1->inactive, &n));
128: PetscCall(VecSet(finemarked, 0.0));
129: for (k = 0; k < n; k++) PetscCall(VecSetValue(finemarked, index[k], 1.0, INSERT_VALUES));
130: PetscCall(VecAssemblyBegin(finemarked));
131: PetscCall(VecAssemblyEnd(finemarked));
133: PetscCall(DMCreateInjection(*dm2, dm1, &inject));
134: PetscCall(MatRestrict(inject, finemarked, coarsemarked));
135: PetscCall(MatDestroy(&inject));
137: /*
138: create index set list of coarse inactive points from coarsemarked
139: */
140: PetscCall(VecGetLocalSize(coarsemarked, &n));
141: PetscCall(VecGetOwnershipRange(coarsemarked, &rstart, NULL));
142: PetscCall(VecGetArray(coarsemarked, &marked));
143: for (k = 0; k < n; k++) {
144: if (marked[k] != 0.0) cnt++;
145: }
146: PetscCall(PetscMalloc1(cnt, &coarseindex));
147: cnt = 0;
148: for (k = 0; k < n; k++) {
149: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
150: }
151: PetscCall(VecRestoreArray(coarsemarked, &marked));
152: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked), cnt, coarseindex, PETSC_OWN_POINTER, &inactive));
154: PetscCall(DMClearGlobalVectors(dm1));
156: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
158: PetscCall(DMSetVI(*dm2, inactive));
160: PetscCall(VecDestroy(&finemarked));
161: PetscCall(VecDestroy(&coarsemarked));
162: PetscCall(ISDestroy(&inactive));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
167: {
168: PetscFunctionBegin;
169: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
170: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
171: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
172: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
173: dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection;
174: dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
175: /* need to clear out this vectors because some of them may not have a reference to the DM
176: but they are counted as having references to the DM in DMDestroy() */
177: PetscCall(DMClearGlobalVectors(dmsnesvi->dm));
179: PetscCall(ISDestroy(&dmsnesvi->inactive));
180: PetscCall(PetscFree(dmsnesvi));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: /*@
185: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
186: be restricted to only those variables NOT associated with active constraints.
188: Logically Collective
190: Input Parameters:
191: + dm - the `DM` object
192: - inactive - an `IS` indicating which points are currently not active
194: Level: intermediate
196: .seealso: `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`
197: @*/
198: PetscErrorCode DMSetVI(DM dm, IS inactive)
199: {
200: PetscContainer isnes;
201: DM_SNESVI *dmsnesvi;
203: PetscFunctionBegin;
204: if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
206: PetscCall(PetscObjectReference((PetscObject)inactive));
208: PetscCall(PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes));
209: if (!isnes) {
210: PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)dm), &isnes));
211: PetscCall(PetscContainerSetUserDestroy(isnes, (PetscErrorCode(*)(void *))DMDestroy_SNESVI));
212: PetscCall(PetscNew(&dmsnesvi));
213: PetscCall(PetscContainerSetPointer(isnes, (void *)dmsnesvi));
214: PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)isnes));
215: PetscCall(PetscContainerDestroy(&isnes));
217: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
218: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
219: dmsnesvi->coarsen = dm->ops->coarsen;
220: dm->ops->coarsen = DMCoarsen_SNESVI;
221: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
222: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
223: dmsnesvi->createinjection = dm->ops->createinjection;
224: dm->ops->createinjection = NULL;
225: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
226: dm->ops->hascreateinjection = NULL;
227: } else {
228: PetscCall(PetscContainerGetPointer(isnes, (void **)&dmsnesvi));
229: PetscCall(ISDestroy(&dmsnesvi->inactive));
230: }
231: PetscCall(DMClearGlobalVectors(dm));
232: PetscCall(ISGetLocalSize(inactive, &dmsnesvi->n));
234: dmsnesvi->inactive = inactive;
235: dmsnesvi->dm = dm;
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*
240: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
241: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
242: */
243: PetscErrorCode DMDestroyVI(DM dm)
244: {
245: PetscFunctionBegin;
246: if (!dm) PetscFunctionReturn(PETSC_SUCCESS);
247: PetscCall(PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)NULL));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes, Vec X, Vec F, IS *ISact, IS *ISinact)
252: {
253: PetscFunctionBegin;
254: PetscCall(SNESVIGetActiveSetIS(snes, X, F, ISact));
255: PetscCall(ISComplement(*ISact, X->map->rstart, X->map->rend, ISinact));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
260: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
261: {
262: Vec v;
264: PetscFunctionBegin;
265: PetscCall(VecCreate(PetscObjectComm((PetscObject)snes), &v));
266: PetscCall(VecSetSizes(v, n, PETSC_DECIDE));
267: PetscCall(VecSetType(v, VECSTANDARD));
268: *newv = v;
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: /* Resets the snes PC and KSP when the active set sizes change */
273: PetscErrorCode SNESVIResetPCandKSP(SNES snes, Mat Amat, Mat Pmat)
274: {
275: KSP snesksp;
277: PetscFunctionBegin;
278: PetscCall(SNESGetKSP(snes, &snesksp));
279: PetscCall(KSPReset(snesksp));
280: PetscCall(KSPResetFromOptions(snesksp));
282: /*
283: KSP kspnew;
284: PC pcnew;
285: MatSolverType stype;
287: PetscCall(KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew));
288: kspnew->pc_side = snesksp->pc_side;
289: kspnew->rtol = snesksp->rtol;
290: kspnew->abstol = snesksp->abstol;
291: kspnew->max_it = snesksp->max_it;
292: PetscCall(KSPSetType(kspnew,((PetscObject)snesksp)->type_name));
293: PetscCall(KSPGetPC(kspnew,&pcnew));
294: PetscCall(PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name));
295: PetscCall(PCSetOperators(kspnew->pc,Amat,Pmat));
296: PetscCall(PCFactorGetMatSolverType(snesksp->pc,&stype));
297: PetscCall(PCFactorSetMatSolverType(kspnew->pc,stype));
298: PetscCall(KSPDestroy(&snesksp));
299: snes->ksp = kspnew;
300: PetscCall(KSPSetFromOptions(kspnew));*/
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
305: implemented in this algorithm. It basically identifies the active constraints and does
306: a linear solve on the other variables (those not associated with the active constraints). */
307: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
308: {
309: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
310: PetscInt maxits, i, lits;
311: SNESLineSearchReason lssucceed;
312: PetscReal fnorm, gnorm, xnorm = 0, ynorm;
313: Vec Y, X, F;
314: KSPConvergedReason kspreason;
315: KSP ksp;
316: PC pc;
318: PetscFunctionBegin;
319: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
320: PetscCall(SNESGetKSP(snes, &ksp));
321: PetscCall(KSPGetPC(ksp, &pc));
322: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH));
324: snes->numFailures = 0;
325: snes->numLinearSolveFailures = 0;
326: snes->reason = SNES_CONVERGED_ITERATING;
328: maxits = snes->max_its; /* maximum number of iterations */
329: X = snes->vec_sol; /* solution vector */
330: F = snes->vec_func; /* residual vector */
331: Y = snes->work[0]; /* work vectors */
333: PetscCall(SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm));
334: PetscCall(SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL));
335: PetscCall(SNESLineSearchSetUp(snes->linesearch));
337: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
338: snes->iter = 0;
339: snes->norm = 0.0;
340: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
342: PetscCall(SNESVIProjectOntoBounds(snes, X));
343: PetscCall(SNESComputeFunction(snes, X, F));
344: PetscCall(SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm));
345: PetscCall(VecNorm(X, NORM_2, &xnorm)); /* xnorm <- ||x|| */
346: SNESCheckFunctionNorm(snes, fnorm);
347: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
348: snes->norm = fnorm;
349: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
350: PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
351: PetscCall(SNESMonitor(snes, 0, fnorm));
353: /* test convergence */
354: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
355: if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
357: for (i = 0; i < maxits; i++) {
358: IS IS_act; /* _act -> active set _inact -> inactive set */
359: IS IS_redact; /* redundant active set */
360: VecScatter scat_act, scat_inact;
361: PetscInt nis_act, nis_inact;
362: Vec Y_act, Y_inact, F_inact;
363: Mat jac_inact_inact, prejac_inact_inact;
364: PetscBool isequal;
366: /* Call general purpose update function */
367: PetscTryTypeMethod(snes, update, snes->iter);
368: PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
369: SNESCheckJacobianDomainerror(snes);
371: /* Create active and inactive index sets */
373: /*original
374: PetscCall(SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact));
375: */
376: PetscCall(SNESVIGetActiveSetIS(snes, X, F, &IS_act));
378: if (vi->checkredundancy) {
379: PetscCall((*vi->checkredundancy)(snes, IS_act, &IS_redact, vi->ctxP));
380: if (IS_redact) {
381: PetscCall(ISSort(IS_redact));
382: PetscCall(ISComplement(IS_redact, X->map->rstart, X->map->rend, &vi->IS_inact));
383: PetscCall(ISDestroy(&IS_redact));
384: } else {
385: PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
386: }
387: } else {
388: PetscCall(ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact));
389: }
391: /* Create inactive set submatrix */
392: PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
394: if (0) { /* Dead code (temporary developer hack) */
395: IS keptrows;
396: PetscCall(MatFindNonzeroRows(jac_inact_inact, &keptrows));
397: if (keptrows) {
398: PetscInt cnt, *nrows, k;
399: const PetscInt *krows, *inact;
400: PetscInt rstart;
402: PetscCall(MatGetOwnershipRange(jac_inact_inact, &rstart, NULL));
403: PetscCall(MatDestroy(&jac_inact_inact));
404: PetscCall(ISDestroy(&IS_act));
406: PetscCall(ISGetLocalSize(keptrows, &cnt));
407: PetscCall(ISGetIndices(keptrows, &krows));
408: PetscCall(ISGetIndices(vi->IS_inact, &inact));
409: PetscCall(PetscMalloc1(cnt, &nrows));
410: for (k = 0; k < cnt; k++) nrows[k] = inact[krows[k] - rstart];
411: PetscCall(ISRestoreIndices(keptrows, &krows));
412: PetscCall(ISRestoreIndices(vi->IS_inact, &inact));
413: PetscCall(ISDestroy(&keptrows));
414: PetscCall(ISDestroy(&vi->IS_inact));
416: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), cnt, nrows, PETSC_OWN_POINTER, &vi->IS_inact));
417: PetscCall(ISComplement(vi->IS_inact, F->map->rstart, F->map->rend, &IS_act));
418: PetscCall(MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact));
419: }
420: }
421: PetscCall(DMSetVI(snes->dm, vi->IS_inact));
422: /* remove later */
424: /*
425: PetscCall(VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm)));
426: PetscCall(VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm)));
427: PetscCall(VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X))));
428: PetscCall(VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F))));
429: PetscCall(ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact))));
430: */
432: /* Get sizes of active and inactive sets */
433: PetscCall(ISGetLocalSize(IS_act, &nis_act));
434: PetscCall(ISGetLocalSize(vi->IS_inact, &nis_inact));
436: /* Create active and inactive set vectors */
437: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact));
438: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act));
439: PetscCall(SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact));
441: /* Create scatter contexts */
442: PetscCall(VecScatterCreate(Y, IS_act, Y_act, NULL, &scat_act));
443: PetscCall(VecScatterCreate(Y, vi->IS_inact, Y_inact, NULL, &scat_inact));
445: /* Do a vec scatter to active and inactive set vectors */
446: PetscCall(VecScatterBegin(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
447: PetscCall(VecScatterEnd(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD));
449: PetscCall(VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
450: PetscCall(VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD));
452: PetscCall(VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
453: PetscCall(VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD));
455: /* Active set direction = 0 */
456: PetscCall(VecSet(Y_act, 0));
457: if (snes->jacobian != snes->jacobian_pre) {
458: PetscCall(MatCreateSubMatrix(snes->jacobian_pre, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &prejac_inact_inact));
459: } else prejac_inact_inact = jac_inact_inact;
461: PetscCall(ISEqual(vi->IS_inact_prev, vi->IS_inact, &isequal));
462: if (!isequal) {
463: PetscCall(SNESVIResetPCandKSP(snes, jac_inact_inact, prejac_inact_inact));
464: PetscCall(PCFieldSplitRestrictIS(pc, vi->IS_inact));
465: }
467: /* PetscCall(ISView(vi->IS_inact,0)); */
468: /* PetscCall(ISView(IS_act,0));*/
469: /* ierr = MatView(snes->jacobian_pre,0); */
471: PetscCall(KSPSetOperators(snes->ksp, jac_inact_inact, prejac_inact_inact));
472: PetscCall(KSPSetUp(snes->ksp));
473: {
474: PC pc;
475: PetscBool flg;
476: PetscCall(KSPGetPC(snes->ksp, &pc));
477: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg));
478: if (flg) {
479: KSP *subksps;
480: PetscCall(PCFieldSplitGetSubKSP(pc, NULL, &subksps));
481: PetscCall(KSPGetPC(subksps[0], &pc));
482: PetscCall(PetscFree(subksps));
483: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &flg));
484: if (flg) {
485: PetscInt n, N = 101 * 101, j, cnts[3] = {0, 0, 0};
486: const PetscInt *ii;
488: PetscCall(ISGetSize(vi->IS_inact, &n));
489: PetscCall(ISGetIndices(vi->IS_inact, &ii));
490: for (j = 0; j < n; j++) {
491: if (ii[j] < N) cnts[0]++;
492: else if (ii[j] < 2 * N) cnts[1]++;
493: else if (ii[j] < 3 * N) cnts[2]++;
494: }
495: PetscCall(ISRestoreIndices(vi->IS_inact, &ii));
497: PetscCall(PCBJacobiSetTotalBlocks(pc, 3, cnts));
498: }
499: }
500: }
502: PetscCall(KSPSolve(snes->ksp, F_inact, Y_inact));
503: PetscCall(VecScatterBegin(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
504: PetscCall(VecScatterEnd(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE));
505: PetscCall(VecScatterBegin(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
506: PetscCall(VecScatterEnd(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE));
508: PetscCall(VecDestroy(&F_inact));
509: PetscCall(VecDestroy(&Y_act));
510: PetscCall(VecDestroy(&Y_inact));
511: PetscCall(VecScatterDestroy(&scat_act));
512: PetscCall(VecScatterDestroy(&scat_inact));
513: PetscCall(ISDestroy(&IS_act));
514: if (!isequal) {
515: PetscCall(ISDestroy(&vi->IS_inact_prev));
516: PetscCall(ISDuplicate(vi->IS_inact, &vi->IS_inact_prev));
517: }
518: PetscCall(ISDestroy(&vi->IS_inact));
519: PetscCall(MatDestroy(&jac_inact_inact));
520: if (snes->jacobian != snes->jacobian_pre) PetscCall(MatDestroy(&prejac_inact_inact));
522: PetscCall(KSPGetConvergedReason(snes->ksp, &kspreason));
523: if (kspreason < 0) {
524: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
525: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures));
526: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
527: break;
528: }
529: }
531: PetscCall(KSPGetIterationNumber(snes->ksp, &lits));
532: snes->linear_its += lits;
533: PetscCall(PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits));
534: /*
535: if (snes->ops->precheck) {
536: PetscBool changed_y = PETSC_FALSE;
537: PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
538: }
540: if (PetscLogPrintInfo) PetscCall(SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W));
541: */
542: /* Compute a (scaled) negative update in the line search routine:
543: Y <- X - lambda*Y
544: and evaluate G = function(Y) (depends on the line search).
545: */
546: PetscCall(VecCopy(Y, snes->vec_sol_update));
547: ynorm = 1;
548: gnorm = fnorm;
549: PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y));
550: PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
551: PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm));
552: PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lssucceed));
553: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
554: if (snes->domainerror) {
555: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
556: PetscCall(DMDestroyVI(snes->dm));
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
559: if (lssucceed) {
560: if (++snes->numFailures >= snes->maxFailures) {
561: PetscBool ismin;
562: snes->reason = SNES_DIVERGED_LINE_SEARCH;
563: PetscCall(SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin));
564: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
565: break;
566: }
567: }
568: PetscCall(DMDestroyVI(snes->dm));
569: /* Update function and solution vectors */
570: fnorm = gnorm;
571: /* Monitor convergence */
572: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
573: snes->iter = i + 1;
574: snes->norm = fnorm;
575: snes->xnorm = xnorm;
576: snes->ynorm = ynorm;
577: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
578: PetscCall(SNESLogConvergenceHistory(snes, snes->norm, lits));
579: PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
580: /* Test for convergence, xnorm = || X || */
581: if (snes->ops->converged != SNESConvergedSkip) PetscCall(VecNorm(X, NORM_2, &xnorm));
582: PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
583: if (snes->reason) break;
584: }
585: /* make sure that the VI information attached to the DM is removed if the for loop above was broken early due to some exceptional conditional */
586: PetscCall(DMDestroyVI(snes->dm));
587: if (i == maxits) {
588: PetscCall(PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits));
589: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
590: }
591: PetscFunctionReturn(PETSC_SUCCESS);
592: }
594: /*@C
595: SNESVISetRedundancyCheck - Provide a function to check for any redundancy in the VI active set
597: Logically Collective
599: Input Parameters:
600: + snes - the `SNESVINEWTONRSLS` context
601: . func - the function to check of redundancies
602: - ctx - optional context used by the function
604: Level: advanced
606: Note:
607: Sometimes the inactive set will result in a non-singular sub-Jacobian problem that needs to be solved, this allows the user,
608: when they know more about their specific problem to provide a function that removes the redundancy that results in the singular linear system
610: .seealso: `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`, `DMSetVI()`
611: @*/
612: PetscErrorCode SNESVISetRedundancyCheck(SNES snes, PetscErrorCode (*func)(SNES, IS, IS *, void *), void *ctx)
613: {
614: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
616: PetscFunctionBegin;
618: vi->checkredundancy = func;
619: vi->ctxP = ctx;
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: #if defined(PETSC_HAVE_MATLAB)
624: #include <engine.h>
625: #include <mex.h>
626: typedef struct {
627: char *funcname;
628: mxArray *ctx;
629: } SNESMatlabContext;
631: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, void *ctx)
632: {
633: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
634: int nlhs = 1, nrhs = 5;
635: mxArray *plhs[1], *prhs[5];
636: long long int l1 = 0, l2 = 0, ls = 0;
637: PetscInt *indices = NULL;
639: PetscFunctionBegin;
643: PetscCheckSameComm(snes, 1, is_act, 2);
645: /* Create IS for reduced active set of size 0, its size and indices will
646: bet set by the Matlab function */
647: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), 0, indices, PETSC_OWN_POINTER, is_redact));
648: /* call Matlab function in ctx */
649: PetscCall(PetscArraycpy(&ls, &snes, 1));
650: PetscCall(PetscArraycpy(&l1, &is_act, 1));
651: PetscCall(PetscArraycpy(&l2, is_redact, 1));
652: prhs[0] = mxCreateDoubleScalar((double)ls);
653: prhs[1] = mxCreateDoubleScalar((double)l1);
654: prhs[2] = mxCreateDoubleScalar((double)l2);
655: prhs[3] = mxCreateString(sctx->funcname);
656: prhs[4] = sctx->ctx;
657: PetscCall(mexCallMATLAB(nlhs, plhs, nrhs, prhs, "PetscSNESVIRedundancyCheckInternal"));
658: PetscCall(mxGetScalar(plhs[0]));
659: mxDestroyArray(prhs[0]);
660: mxDestroyArray(prhs[1]);
661: mxDestroyArray(prhs[2]);
662: mxDestroyArray(prhs[3]);
663: mxDestroyArray(plhs[0]);
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
668: {
669: SNESMatlabContext *sctx;
671: PetscFunctionBegin;
672: /* currently sctx is memory bleed */
673: PetscCall(PetscNew(&sctx));
674: PetscCall(PetscStrallocpy(func, &sctx->funcname));
675: sctx->ctx = mxDuplicateArray(ctx);
676: PetscCall(SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx));
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: #endif
682: /*
683: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
684: of the SNESVI nonlinear solver.
686: Input Parameter:
687: . snes - the SNES context
689: Application Interface Routine: SNESSetUp()
691: Note:
692: For basic use of the SNES solvers, the user need not explicitly call
693: SNESSetUp(), since these actions will automatically occur during
694: the call to SNESSolve().
695: */
696: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
697: {
698: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
699: PetscInt *indices;
700: PetscInt i, n, rstart, rend;
701: SNESLineSearch linesearch;
703: PetscFunctionBegin;
704: PetscCall(SNESSetUp_VI(snes));
706: /* Set up previous active index set for the first snes solve
707: vi->IS_inact_prev = 0,1,2,....N */
709: PetscCall(VecGetOwnershipRange(snes->vec_sol, &rstart, &rend));
710: PetscCall(VecGetLocalSize(snes->vec_sol, &n));
711: PetscCall(PetscMalloc1(n, &indices));
712: for (i = 0; i < n; i++) indices[i] = rstart + i;
713: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev));
715: /* set the line search functions */
716: if (!snes->linesearch) {
717: PetscCall(SNESGetLineSearch(snes, &linesearch));
718: PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
719: }
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
724: {
725: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
727: PetscFunctionBegin;
728: PetscCall(SNESReset_VI(snes));
729: PetscCall(ISDestroy(&vi->IS_inact_prev));
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: /*MC
734: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
736: Options Database Keys:
737: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
738: - -snes_vi_monitor - prints the number of active constraints at each iteration.
740: Level: beginner
742: References:
743: . * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
744: Applications, Optimization Methods and Software, 21 (2006).
746: Note:
747: At each set of this methods the algorithm produces an inactive set of variables that are constrained to their current values
748: (because changing these values would result in those variables no longer satisfying the inequality constraints)
749: and produces a step direction by solving the linear system arising from the Jacobian with the inactive variables removed. In other
750: words on a reduced space of the solution space. Based on the Newton update it then adjusts the inactive sep for the next iteration.
752: .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTR`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
753: M*/
754: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
755: {
756: SNES_VINEWTONRSLS *vi;
757: SNESLineSearch linesearch;
759: PetscFunctionBegin;
760: snes->ops->reset = SNESReset_VINEWTONRSLS;
761: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
762: snes->ops->solve = SNESSolve_VINEWTONRSLS;
763: snes->ops->destroy = SNESDestroy_VI;
764: snes->ops->setfromoptions = SNESSetFromOptions_VI;
765: snes->ops->view = NULL;
766: snes->ops->converged = SNESConvergedDefault_VI;
768: snes->usesksp = PETSC_TRUE;
769: snes->usesnpc = PETSC_FALSE;
771: PetscCall(SNESGetLineSearch(snes, &linesearch));
772: if (!((PetscObject)linesearch)->type_name) PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBT));
773: PetscCall(SNESLineSearchBTSetAlpha(linesearch, 0.0));
775: snes->alwayscomputesfinalresidual = PETSC_TRUE;
777: PetscCall(PetscNew(&vi));
778: snes->data = (void *)vi;
779: vi->checkredundancy = NULL;
781: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI));
782: PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI));
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }