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: *inact = vi->IS_inact;
25: return 0;
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: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm, Vec *vec)
51: {
52: PetscContainer isnes;
53: DM_SNESVI *dmsnesvi;
55: PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes);
57: PetscContainerGetPointer(isnes, (void **)&dmsnesvi);
58: VecCreateMPI(PetscObjectComm((PetscObject)dm), dmsnesvi->n, PETSC_DETERMINE, vec);
59: VecSetDM(*vec, dm);
60: return 0;
61: }
63: /*
64: DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
65: */
66: PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1, DM dm2, Mat *mat, Vec *vec)
67: {
68: PetscContainer isnes;
69: DM_SNESVI *dmsnesvi1, *dmsnesvi2;
70: Mat interp;
72: PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes);
74: PetscContainerGetPointer(isnes, (void **)&dmsnesvi1);
75: PetscObjectQuery((PetscObject)dm2, "VI", (PetscObject *)&isnes);
77: PetscContainerGetPointer(isnes, (void **)&dmsnesvi2);
79: (*dmsnesvi1->createinterpolation)(dm1, dm2, &interp, NULL);
80: MatCreateSubMatrix(interp, dmsnesvi2->inactive, dmsnesvi1->inactive, MAT_INITIAL_MATRIX, mat);
81: MatDestroy(&interp);
82: *vec = NULL;
83: return 0;
84: }
86: /*
87: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
88: */
89: PetscErrorCode DMCoarsen_SNESVI(DM dm1, MPI_Comm comm, DM *dm2)
90: {
91: PetscContainer isnes;
92: DM_SNESVI *dmsnesvi1;
93: Vec finemarked, coarsemarked;
94: IS inactive;
95: Mat inject;
96: const PetscInt *index;
97: PetscInt n, k, cnt = 0, rstart, *coarseindex;
98: PetscScalar *marked;
100: PetscObjectQuery((PetscObject)dm1, "VI", (PetscObject *)&isnes);
102: PetscContainerGetPointer(isnes, (void **)&dmsnesvi1);
104: /* get the original coarsen */
105: (*dmsnesvi1->coarsen)(dm1, comm, dm2);
107: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
108: /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
109: /* PetscObjectReference((PetscObject)*dm2);*/
111: /* need to set back global vectors in order to use the original injection */
112: DMClearGlobalVectors(dm1);
114: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
116: DMCreateGlobalVector(dm1, &finemarked);
117: DMCreateGlobalVector(*dm2, &coarsemarked);
119: /*
120: fill finemarked with locations of inactive points
121: */
122: ISGetIndices(dmsnesvi1->inactive, &index);
123: ISGetLocalSize(dmsnesvi1->inactive, &n);
124: VecSet(finemarked, 0.0);
125: for (k = 0; k < n; k++) VecSetValue(finemarked, index[k], 1.0, INSERT_VALUES);
126: VecAssemblyBegin(finemarked);
127: VecAssemblyEnd(finemarked);
129: DMCreateInjection(*dm2, dm1, &inject);
130: MatRestrict(inject, finemarked, coarsemarked);
131: MatDestroy(&inject);
133: /*
134: create index set list of coarse inactive points from coarsemarked
135: */
136: VecGetLocalSize(coarsemarked, &n);
137: VecGetOwnershipRange(coarsemarked, &rstart, NULL);
138: VecGetArray(coarsemarked, &marked);
139: for (k = 0; k < n; k++) {
140: if (marked[k] != 0.0) cnt++;
141: }
142: PetscMalloc1(cnt, &coarseindex);
143: cnt = 0;
144: for (k = 0; k < n; k++) {
145: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
146: }
147: VecRestoreArray(coarsemarked, &marked);
148: ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked), cnt, coarseindex, PETSC_OWN_POINTER, &inactive);
150: DMClearGlobalVectors(dm1);
152: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
154: DMSetVI(*dm2, inactive);
156: VecDestroy(&finemarked);
157: VecDestroy(&coarsemarked);
158: ISDestroy(&inactive);
159: return 0;
160: }
162: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
163: {
164: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
165: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
166: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
167: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
168: dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection;
169: dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
170: /* need to clear out this vectors because some of them may not have a reference to the DM
171: but they are counted as having references to the DM in DMDestroy() */
172: DMClearGlobalVectors(dmsnesvi->dm);
174: ISDestroy(&dmsnesvi->inactive);
175: PetscFree(dmsnesvi);
176: return 0;
177: }
179: /*@
180: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
181: be restricted to only those variables NOT associated with active constraints.
183: Logically Collective
185: Input Parameters:
186: + dm - the `DM` object
187: - inactive - an `IS` indicating which points are currently not active
189: Level: intermediate
191: .seealso: `SNESVINEWTONRSLS`, `SNESVIGetInactiveSet()`
192: @*/
193: PetscErrorCode DMSetVI(DM dm, IS inactive)
194: {
195: PetscContainer isnes;
196: DM_SNESVI *dmsnesvi;
198: if (!dm) return 0;
200: PetscObjectReference((PetscObject)inactive);
202: PetscObjectQuery((PetscObject)dm, "VI", (PetscObject *)&isnes);
203: if (!isnes) {
204: PetscContainerCreate(PetscObjectComm((PetscObject)dm), &isnes);
205: PetscContainerSetUserDestroy(isnes, (PetscErrorCode(*)(void *))DMDestroy_SNESVI);
206: PetscNew(&dmsnesvi);
207: PetscContainerSetPointer(isnes, (void *)dmsnesvi);
208: PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)isnes);
209: PetscContainerDestroy(&isnes);
211: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
212: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
213: dmsnesvi->coarsen = dm->ops->coarsen;
214: dm->ops->coarsen = DMCoarsen_SNESVI;
215: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
216: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
217: dmsnesvi->createinjection = dm->ops->createinjection;
218: dm->ops->createinjection = NULL;
219: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
220: dm->ops->hascreateinjection = NULL;
221: } else {
222: PetscContainerGetPointer(isnes, (void **)&dmsnesvi);
223: ISDestroy(&dmsnesvi->inactive);
224: }
225: DMClearGlobalVectors(dm);
226: ISGetLocalSize(inactive, &dmsnesvi->n);
228: dmsnesvi->inactive = inactive;
229: dmsnesvi->dm = dm;
230: return 0;
231: }
233: /*
234: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
235: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
236: */
237: PetscErrorCode DMDestroyVI(DM dm)
238: {
239: if (!dm) return 0;
240: PetscObjectCompose((PetscObject)dm, "VI", (PetscObject)NULL);
241: return 0;
242: }
244: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes, Vec X, Vec F, IS *ISact, IS *ISinact)
245: {
246: SNESVIGetActiveSetIS(snes, X, F, ISact);
247: ISComplement(*ISact, X->map->rstart, X->map->rend, ISinact);
248: return 0;
249: }
251: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
252: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes, PetscInt n, Vec *newv)
253: {
254: Vec v;
256: VecCreate(PetscObjectComm((PetscObject)snes), &v);
257: VecSetSizes(v, n, PETSC_DECIDE);
258: VecSetType(v, VECSTANDARD);
259: *newv = v;
260: return 0;
261: }
263: /* Resets the snes PC and KSP when the active set sizes change */
264: PetscErrorCode SNESVIResetPCandKSP(SNES snes, Mat Amat, Mat Pmat)
265: {
266: KSP snesksp;
268: SNESGetKSP(snes, &snesksp);
269: KSPReset(snesksp);
270: KSPResetFromOptions(snesksp);
272: /*
273: KSP kspnew;
274: PC pcnew;
275: MatSolverType stype;
277: KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
278: kspnew->pc_side = snesksp->pc_side;
279: kspnew->rtol = snesksp->rtol;
280: kspnew->abstol = snesksp->abstol;
281: kspnew->max_it = snesksp->max_it;
282: KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
283: KSPGetPC(kspnew,&pcnew);
284: PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
285: PCSetOperators(kspnew->pc,Amat,Pmat);
286: PCFactorGetMatSolverType(snesksp->pc,&stype);
287: PCFactorSetMatSolverType(kspnew->pc,stype);
288: KSPDestroy(&snesksp);
289: snes->ksp = kspnew;
290: KSPSetFromOptions(kspnew);*/
291: return 0;
292: }
294: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
295: implemented in this algorithm. It basically identifies the active constraints and does
296: a linear solve on the other variables (those not associated with the active constraints). */
297: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
298: {
299: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
300: PetscInt maxits, i, lits;
301: SNESLineSearchReason lssucceed;
302: PetscReal fnorm, gnorm, xnorm = 0, ynorm;
303: Vec Y, X, F;
304: KSPConvergedReason kspreason;
305: KSP ksp;
306: PC pc;
308: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
309: SNESGetKSP(snes, &ksp);
310: KSPGetPC(ksp, &pc);
311: PCMGSetGalerkin(pc, PC_MG_GALERKIN_BOTH);
313: snes->numFailures = 0;
314: snes->numLinearSolveFailures = 0;
315: snes->reason = SNES_CONVERGED_ITERATING;
317: maxits = snes->max_its; /* maximum number of iterations */
318: X = snes->vec_sol; /* solution vector */
319: F = snes->vec_func; /* residual vector */
320: Y = snes->work[0]; /* work vectors */
322: SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
323: SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
324: SNESLineSearchSetUp(snes->linesearch);
326: PetscObjectSAWsTakeAccess((PetscObject)snes);
327: snes->iter = 0;
328: snes->norm = 0.0;
329: PetscObjectSAWsGrantAccess((PetscObject)snes);
331: SNESVIProjectOntoBounds(snes, X);
332: SNESComputeFunction(snes, X, F);
333: SNESVIComputeInactiveSetFnorm(snes, F, X, &fnorm);
334: VecNorm(X, NORM_2, &xnorm); /* xnorm <- ||x|| */
335: SNESCheckFunctionNorm(snes, fnorm);
336: PetscObjectSAWsTakeAccess((PetscObject)snes);
337: snes->norm = fnorm;
338: PetscObjectSAWsGrantAccess((PetscObject)snes);
339: SNESLogConvergenceHistory(snes, fnorm, 0);
340: SNESMonitor(snes, 0, fnorm);
342: /* test convergence */
343: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
344: if (snes->reason) return 0;
346: for (i = 0; i < maxits; i++) {
347: IS IS_act; /* _act -> active set _inact -> inactive set */
348: IS IS_redact; /* redundant active set */
349: VecScatter scat_act, scat_inact;
350: PetscInt nis_act, nis_inact;
351: Vec Y_act, Y_inact, F_inact;
352: Mat jac_inact_inact, prejac_inact_inact;
353: PetscBool isequal;
355: /* Call general purpose update function */
356: PetscTryTypeMethod(snes, update, snes->iter);
357: SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
358: SNESCheckJacobianDomainerror(snes);
360: /* Create active and inactive index sets */
362: /*original
363: SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
364: */
365: SNESVIGetActiveSetIS(snes, X, F, &IS_act);
367: if (vi->checkredundancy) {
368: (*vi->checkredundancy)(snes, IS_act, &IS_redact, vi->ctxP);
369: if (IS_redact) {
370: ISSort(IS_redact);
371: ISComplement(IS_redact, X->map->rstart, X->map->rend, &vi->IS_inact);
372: ISDestroy(&IS_redact);
373: } else {
374: ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact);
375: }
376: } else {
377: ISComplement(IS_act, X->map->rstart, X->map->rend, &vi->IS_inact);
378: }
380: /* Create inactive set submatrix */
381: MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact);
383: if (0) { /* Dead code (temporary developer hack) */
384: IS keptrows;
385: MatFindNonzeroRows(jac_inact_inact, &keptrows);
386: if (keptrows) {
387: PetscInt cnt, *nrows, k;
388: const PetscInt *krows, *inact;
389: PetscInt rstart;
391: MatGetOwnershipRange(jac_inact_inact, &rstart, NULL);
392: MatDestroy(&jac_inact_inact);
393: ISDestroy(&IS_act);
395: ISGetLocalSize(keptrows, &cnt);
396: ISGetIndices(keptrows, &krows);
397: ISGetIndices(vi->IS_inact, &inact);
398: PetscMalloc1(cnt, &nrows);
399: for (k = 0; k < cnt; k++) nrows[k] = inact[krows[k] - rstart];
400: ISRestoreIndices(keptrows, &krows);
401: ISRestoreIndices(vi->IS_inact, &inact);
402: ISDestroy(&keptrows);
403: ISDestroy(&vi->IS_inact);
405: ISCreateGeneral(PetscObjectComm((PetscObject)snes), cnt, nrows, PETSC_OWN_POINTER, &vi->IS_inact);
406: ISComplement(vi->IS_inact, F->map->rstart, F->map->rend, &IS_act);
407: MatCreateSubMatrix(snes->jacobian, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &jac_inact_inact);
408: }
409: }
410: DMSetVI(snes->dm, vi->IS_inact);
411: /* remove later */
413: /*
414: VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
415: VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
416: VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
417: VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
418: ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
419: */
421: /* Get sizes of active and inactive sets */
422: ISGetLocalSize(IS_act, &nis_act);
423: ISGetLocalSize(vi->IS_inact, &nis_inact);
425: /* Create active and inactive set vectors */
426: SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &F_inact);
427: SNESCreateSubVectors_VINEWTONRSLS(snes, nis_act, &Y_act);
428: SNESCreateSubVectors_VINEWTONRSLS(snes, nis_inact, &Y_inact);
430: /* Create scatter contexts */
431: VecScatterCreate(Y, IS_act, Y_act, NULL, &scat_act);
432: VecScatterCreate(Y, vi->IS_inact, Y_inact, NULL, &scat_inact);
434: /* Do a vec scatter to active and inactive set vectors */
435: VecScatterBegin(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD);
436: VecScatterEnd(scat_inact, F, F_inact, INSERT_VALUES, SCATTER_FORWARD);
438: VecScatterBegin(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD);
439: VecScatterEnd(scat_act, Y, Y_act, INSERT_VALUES, SCATTER_FORWARD);
441: VecScatterBegin(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD);
442: VecScatterEnd(scat_inact, Y, Y_inact, INSERT_VALUES, SCATTER_FORWARD);
444: /* Active set direction = 0 */
445: VecSet(Y_act, 0);
446: if (snes->jacobian != snes->jacobian_pre) {
447: MatCreateSubMatrix(snes->jacobian_pre, vi->IS_inact, vi->IS_inact, MAT_INITIAL_MATRIX, &prejac_inact_inact);
448: } else prejac_inact_inact = jac_inact_inact;
450: ISEqual(vi->IS_inact_prev, vi->IS_inact, &isequal);
451: if (!isequal) {
452: SNESVIResetPCandKSP(snes, jac_inact_inact, prejac_inact_inact);
453: PCFieldSplitRestrictIS(pc, vi->IS_inact);
454: }
456: /* ISView(vi->IS_inact,0); */
457: /* ISView(IS_act,0);*/
458: /* MatView(snes->jacobian_pre,0); */
460: KSPSetOperators(snes->ksp, jac_inact_inact, prejac_inact_inact);
461: KSPSetUp(snes->ksp);
462: {
463: PC pc;
464: PetscBool flg;
465: KSPGetPC(snes->ksp, &pc);
466: PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &flg);
467: if (flg) {
468: KSP *subksps;
469: PCFieldSplitGetSubKSP(pc, NULL, &subksps);
470: KSPGetPC(subksps[0], &pc);
471: PetscFree(subksps);
472: PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &flg);
473: if (flg) {
474: PetscInt n, N = 101 * 101, j, cnts[3] = {0, 0, 0};
475: const PetscInt *ii;
477: ISGetSize(vi->IS_inact, &n);
478: ISGetIndices(vi->IS_inact, &ii);
479: for (j = 0; j < n; j++) {
480: if (ii[j] < N) cnts[0]++;
481: else if (ii[j] < 2 * N) cnts[1]++;
482: else if (ii[j] < 3 * N) cnts[2]++;
483: }
484: ISRestoreIndices(vi->IS_inact, &ii);
486: PCBJacobiSetTotalBlocks(pc, 3, cnts);
487: }
488: }
489: }
491: KSPSolve(snes->ksp, F_inact, Y_inact);
492: VecScatterBegin(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE);
493: VecScatterEnd(scat_act, Y_act, Y, INSERT_VALUES, SCATTER_REVERSE);
494: VecScatterBegin(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE);
495: VecScatterEnd(scat_inact, Y_inact, Y, INSERT_VALUES, SCATTER_REVERSE);
497: VecDestroy(&F_inact);
498: VecDestroy(&Y_act);
499: VecDestroy(&Y_inact);
500: VecScatterDestroy(&scat_act);
501: VecScatterDestroy(&scat_inact);
502: ISDestroy(&IS_act);
503: if (!isequal) {
504: ISDestroy(&vi->IS_inact_prev);
505: ISDuplicate(vi->IS_inact, &vi->IS_inact_prev);
506: }
507: ISDestroy(&vi->IS_inact);
508: MatDestroy(&jac_inact_inact);
509: if (snes->jacobian != snes->jacobian_pre) MatDestroy(&prejac_inact_inact);
511: KSPGetConvergedReason(snes->ksp, &kspreason);
512: if (kspreason < 0) {
513: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
514: PetscInfo(snes, "iter=%" PetscInt_FMT ", number linear solve failures %" PetscInt_FMT " greater than current SNES allowed, stopping solve\n", snes->iter, snes->numLinearSolveFailures);
515: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
516: break;
517: }
518: }
520: KSPGetIterationNumber(snes->ksp, &lits);
521: snes->linear_its += lits;
522: PetscInfo(snes, "iter=%" PetscInt_FMT ", linear solve iterations=%" PetscInt_FMT "\n", snes->iter, lits);
523: /*
524: if (snes->ops->precheck) {
525: PetscBool changed_y = PETSC_FALSE;
526: PetscUseTypeMethod(snes,precheck ,X,Y,snes->precheck,&changed_y);
527: }
529: if (PetscLogPrintInfo) SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
530: */
531: /* Compute a (scaled) negative update in the line search routine:
532: Y <- X - lambda*Y
533: and evaluate G = function(Y) (depends on the line search).
534: */
535: VecCopy(Y, snes->vec_sol_update);
536: ynorm = 1;
537: gnorm = fnorm;
538: SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
539: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
540: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
541: PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lssucceed);
542: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
543: if (snes->domainerror) {
544: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
545: DMDestroyVI(snes->dm);
546: return 0;
547: }
548: if (lssucceed) {
549: if (++snes->numFailures >= snes->maxFailures) {
550: PetscBool ismin;
551: snes->reason = SNES_DIVERGED_LINE_SEARCH;
552: SNESVICheckLocalMin_Private(snes, snes->jacobian, F, X, gnorm, &ismin);
553: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
554: break;
555: }
556: }
557: DMDestroyVI(snes->dm);
558: /* Update function and solution vectors */
559: fnorm = gnorm;
560: /* Monitor convergence */
561: PetscObjectSAWsTakeAccess((PetscObject)snes);
562: snes->iter = i + 1;
563: snes->norm = fnorm;
564: snes->xnorm = xnorm;
565: snes->ynorm = ynorm;
566: PetscObjectSAWsGrantAccess((PetscObject)snes);
567: SNESLogConvergenceHistory(snes, snes->norm, lits);
568: SNESMonitor(snes, snes->iter, snes->norm);
569: /* Test for convergence, xnorm = || X || */
570: if (snes->ops->converged != SNESConvergedSkip) VecNorm(X, NORM_2, &xnorm);
571: PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
572: if (snes->reason) break;
573: }
574: /* 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 */
575: DMDestroyVI(snes->dm);
576: if (i == maxits) {
577: PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
578: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
579: }
580: return 0;
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 non-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: `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;
606: vi->checkredundancy = func;
607: vi->ctxP = ctx;
608: return 0;
609: }
611: #if defined(PETSC_HAVE_MATLAB)
612: #include <engine.h>
613: #include <mex.h>
614: typedef struct {
615: char *funcname;
616: mxArray *ctx;
617: } SNESMatlabContext;
619: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes, IS is_act, IS *is_redact, void *ctx)
620: {
621: SNESMatlabContext *sctx = (SNESMatlabContext *)ctx;
622: int nlhs = 1, nrhs = 5;
623: mxArray *plhs[1], *prhs[5];
624: long long int l1 = 0, l2 = 0, ls = 0;
625: PetscInt *indices = NULL;
632: /* Create IS for reduced active set of size 0, its size and indices will
633: bet set by the Matlab function */
634: ISCreateGeneral(PetscObjectComm((PetscObject)snes), 0, indices, PETSC_OWN_POINTER, is_redact);
635: /* call Matlab function in ctx */
636: PetscArraycpy(&ls, &snes, 1);
637: PetscArraycpy(&l1, &is_act, 1);
638: PetscArraycpy(&l2, is_redact, 1);
639: prhs[0] = mxCreateDoubleScalar((double)ls);
640: prhs[1] = mxCreateDoubleScalar((double)l1);
641: prhs[2] = mxCreateDoubleScalar((double)l2);
642: prhs[3] = mxCreateString(sctx->funcname);
643: prhs[4] = sctx->ctx;
644: mexCallMATLAB(nlhs, plhs, nrhs, prhs, "PetscSNESVIRedundancyCheckInternal");
645: mxGetScalar(plhs[0]);
646: mxDestroyArray(prhs[0]);
647: mxDestroyArray(prhs[1]);
648: mxDestroyArray(prhs[2]);
649: mxDestroyArray(prhs[3]);
650: mxDestroyArray(plhs[0]);
651: return 0;
652: }
654: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes, const char *func, mxArray *ctx)
655: {
656: SNESMatlabContext *sctx;
658: /* currently sctx is memory bleed */
659: PetscNew(&sctx);
660: PetscStrallocpy(func, &sctx->funcname);
661: sctx->ctx = mxDuplicateArray(ctx);
662: SNESVISetRedundancyCheck(snes, SNESVIRedundancyCheck_Matlab, sctx);
663: return 0;
664: }
666: #endif
668: /*
669: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
670: of the SNESVI nonlinear solver.
672: Input Parameter:
673: . snes - the SNES context
675: Application Interface Routine: SNESSetUp()
677: Note:
678: For basic use of the SNES solvers, the user need not explicitly call
679: SNESSetUp(), since these actions will automatically occur during
680: the call to SNESSolve().
681: */
682: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
683: {
684: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
685: PetscInt *indices;
686: PetscInt i, n, rstart, rend;
687: SNESLineSearch linesearch;
689: SNESSetUp_VI(snes);
691: /* Set up previous active index set for the first snes solve
692: vi->IS_inact_prev = 0,1,2,....N */
694: VecGetOwnershipRange(snes->vec_sol, &rstart, &rend);
695: VecGetLocalSize(snes->vec_sol, &n);
696: PetscMalloc1(n, &indices);
697: for (i = 0; i < n; i++) indices[i] = rstart + i;
698: ISCreateGeneral(PetscObjectComm((PetscObject)snes), n, indices, PETSC_OWN_POINTER, &vi->IS_inact_prev);
700: /* set the line search functions */
701: if (!snes->linesearch) {
702: SNESGetLineSearch(snes, &linesearch);
703: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
704: }
705: return 0;
706: }
708: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
709: {
710: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS *)snes->data;
712: SNESReset_VI(snes);
713: ISDestroy(&vi->IS_inact_prev);
714: return 0;
715: }
717: /*MC
718: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
720: Options Database Keys:
721: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
722: - -snes_vi_monitor - prints the number of active constraints at each iteration.
724: Level: beginner
726: References:
727: . * - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
728: Applications, Optimization Methods and Software, 21 (2006).
730: Note:
731: At each set of this methods the algorithm produces an inactive set of variables that are constrained to their current values
732: (because changing these values would result in those variables no longer satisfying the inequality constraints)
733: and produces a step direction by solving the linear system arising from the Jacobian with the inactive variables removed. In other
734: words on a reduced space of the solution space. Based on the Newton update it then adjusts the inactive sep for the next iteration.
736: .seealso: `SNESVISetVariableBounds()`, `SNESVISetComputeVariableBounds()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESVINEWTONSSLS`, `SNESNEWTONTRDC`, `SNESLineSearchSetType()`, `SNESLineSearchSetPostCheck()`, `SNESLineSearchSetPreCheck()`, `SNESVIGetInactiveSet()`, `DMSetVI()`, `SNESVISetRedundancyCheck()`
737: M*/
738: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
739: {
740: SNES_VINEWTONRSLS *vi;
741: SNESLineSearch linesearch;
743: snes->ops->reset = SNESReset_VINEWTONRSLS;
744: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
745: snes->ops->solve = SNESSolve_VINEWTONRSLS;
746: snes->ops->destroy = SNESDestroy_VI;
747: snes->ops->setfromoptions = SNESSetFromOptions_VI;
748: snes->ops->view = NULL;
749: snes->ops->converged = SNESConvergedDefault_VI;
751: snes->usesksp = PETSC_TRUE;
752: snes->usesnpc = PETSC_FALSE;
754: SNESGetLineSearch(snes, &linesearch);
755: if (!((PetscObject)linesearch)->type_name) SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
756: SNESLineSearchBTSetAlpha(linesearch, 0.0);
758: snes->alwayscomputesfinalresidual = PETSC_TRUE;
760: PetscNew(&vi);
761: snes->data = (void *)vi;
762: vi->checkredundancy = NULL;
764: PetscObjectComposeFunction((PetscObject)snes, "SNESVISetVariableBounds_C", SNESVISetVariableBounds_VI);
765: PetscObjectComposeFunction((PetscObject)snes, "SNESVISetComputeVariableBounds_C", SNESVISetComputeVariableBounds_VI);
766: return 0;
767: }