Actual source code: virs.c
petsc-3.11.4 2019-09-28
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: */
17: PetscErrorCode SNESVIGetInactiveSet(SNES snes,IS *inact)
18: {
19: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
22: *inact = vi->IS_inact;
23: return(0);
24: }
26: /*
27: 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
28: defined by the reduced space method.
30: Simple calls the regular DM interpolation and restricts it to operation on the variables not associated with active constraints.
32: <*/
33: typedef struct {
34: PetscInt n; /* size of vectors in the reduced DM space */
35: IS inactive;
37: PetscErrorCode (*createinterpolation)(DM,DM,Mat*,Vec*); /* DM's original routines */
38: PetscErrorCode (*coarsen)(DM, MPI_Comm, DM*);
39: PetscErrorCode (*createglobalvector)(DM,Vec*);
40: PetscErrorCode (*getinjection)(DM,DM,Mat*);
41: PetscErrorCode (*hascreateinjection)(DM,PetscBool*);
43: DM dm; /* when destroying this object we need to reset the above function into the base DM */
44: } DM_SNESVI;
46: /*
47: DMCreateGlobalVector_SNESVI - Creates global vector of the size of the reduced space
49: */
50: PetscErrorCode DMCreateGlobalVector_SNESVI(DM dm,Vec *vec)
51: {
53: PetscContainer isnes;
54: DM_SNESVI *dmsnesvi;
57: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
58: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Composed SNES is missing");
59: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
60: VecCreateMPI(PetscObjectComm((PetscObject)dm),dmsnesvi->n,PETSC_DETERMINE,vec);
61: return(0);
62: }
64: static PetscErrorCode DMHasCreateInjection_SNESVI(DM dm, PetscBool *flg)
65: {
69: *flg = PETSC_FALSE;
70: return(0);
71: }
73: /*
74: DMCreateInterpolation_SNESVI - Modifieds the interpolation obtained from the DM by removing all rows and columns associated with active constraints.
76: */
77: PetscErrorCode DMCreateInterpolation_SNESVI(DM dm1,DM dm2,Mat *mat,Vec *vec)
78: {
80: PetscContainer isnes;
81: DM_SNESVI *dmsnesvi1,*dmsnesvi2;
82: Mat interp;
85: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
86: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
87: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
88: PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
89: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing");
90: PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);
92: (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
93: MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
94: MatDestroy(&interp);
95: *vec = 0;
96: return(0);
97: }
99: static PetscErrorCode DMSetVI(DM,IS);
100: static PetscErrorCode DMDestroyVI(DM);
102: /*
103: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
105: */
106: PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
107: {
109: PetscContainer isnes;
110: DM_SNESVI *dmsnesvi1;
111: Vec finemarked,coarsemarked;
112: IS inactive;
113: Mat inject;
114: const PetscInt *index;
115: PetscInt n,k,cnt = 0,rstart,*coarseindex;
116: PetscScalar *marked;
119: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
120: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
121: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
123: /* get the original coarsen */
124: (*dmsnesvi1->coarsen)(dm1,comm,dm2);
126: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
127: /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
128: /* PetscObjectReference((PetscObject)*dm2);*/
130: /* need to set back global vectors in order to use the original injection */
131: DMClearGlobalVectors(dm1);
133: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
135: DMCreateGlobalVector(dm1,&finemarked);
136: DMCreateGlobalVector(*dm2,&coarsemarked);
138: /*
139: fill finemarked with locations of inactive points
140: */
141: ISGetIndices(dmsnesvi1->inactive,&index);
142: ISGetLocalSize(dmsnesvi1->inactive,&n);
143: VecSet(finemarked,0.0);
144: for (k=0; k<n; k++) {
145: VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
146: }
147: VecAssemblyBegin(finemarked);
148: VecAssemblyEnd(finemarked);
150: DMCreateInjection(*dm2,dm1,&inject);
151: MatRestrict(inject,finemarked,coarsemarked);
152: MatDestroy(&inject);
154: /*
155: create index set list of coarse inactive points from coarsemarked
156: */
157: VecGetLocalSize(coarsemarked,&n);
158: VecGetOwnershipRange(coarsemarked,&rstart,NULL);
159: VecGetArray(coarsemarked,&marked);
160: for (k=0; k<n; k++) {
161: if (marked[k] != 0.0) cnt++;
162: }
163: PetscMalloc1(cnt,&coarseindex);
164: cnt = 0;
165: for (k=0; k<n; k++) {
166: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
167: }
168: VecRestoreArray(coarsemarked,&marked);
169: ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);
171: DMClearGlobalVectors(dm1);
173: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
175: DMSetVI(*dm2,inactive);
177: VecDestroy(&finemarked);
178: VecDestroy(&coarsemarked);
179: ISDestroy(&inactive);
180: return(0);
181: }
183: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
184: {
188: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
189: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
190: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
191: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
192: dmsnesvi->dm->ops->getinjection = dmsnesvi->getinjection;
193: dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
194: /* need to clear out this vectors because some of them may not have a reference to the DM
195: but they are counted as having references to the DM in DMDestroy() */
196: DMClearGlobalVectors(dmsnesvi->dm);
198: ISDestroy(&dmsnesvi->inactive);
199: PetscFree(dmsnesvi);
200: return(0);
201: }
203: /*
204: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
205: be restricted to only those variables NOT associated with active constraints.
207: */
208: static PetscErrorCode DMSetVI(DM dm,IS inactive)
209: {
211: PetscContainer isnes;
212: DM_SNESVI *dmsnesvi;
215: if (!dm) return(0);
217: PetscObjectReference((PetscObject)inactive);
219: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
220: if (!isnes) {
221: PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
222: PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
223: PetscNew(&dmsnesvi);
224: PetscContainerSetPointer(isnes,(void*)dmsnesvi);
225: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
226: PetscContainerDestroy(&isnes);
228: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
229: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
230: dmsnesvi->coarsen = dm->ops->coarsen;
231: dm->ops->coarsen = DMCoarsen_SNESVI;
232: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
233: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
234: dmsnesvi->getinjection = dm->ops->getinjection;
235: dm->ops->getinjection = NULL;
236: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
237: dm->ops->hascreateinjection = DMHasCreateInjection_SNESVI;
238: } else {
239: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
240: ISDestroy(&dmsnesvi->inactive);
241: }
242: DMClearGlobalVectors(dm);
243: ISGetLocalSize(inactive,&dmsnesvi->n);
245: dmsnesvi->inactive = inactive;
246: dmsnesvi->dm = dm;
247: return(0);
248: }
250: /*
251: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
252: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
253: */
254: static PetscErrorCode DMDestroyVI(DM dm)
255: {
259: if (!dm) return(0);
260: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
261: return(0);
262: }
264: /* --------------------------------------------------------------------------------------------------------*/
267: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
268: {
272: SNESVIGetActiveSetIS(snes,X,F,ISact);
273: ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
274: return(0);
275: }
277: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
278: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv)
279: {
281: Vec v;
284: VecCreate(PetscObjectComm((PetscObject)snes),&v);
285: VecSetSizes(v,n,PETSC_DECIDE);
286: VecSetType(v,VECSTANDARD);
287: *newv = v;
288: return(0);
289: }
291: /* Resets the snes PC and KSP when the active set sizes change */
292: PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)
293: {
295: KSP snesksp;
298: SNESGetKSP(snes,&snesksp);
299: KSPReset(snesksp);
300: KSPResetFromOptions(snesksp);
302: /*
303: KSP kspnew;
304: PC pcnew;
305: MatSolverType stype;
308: KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
309: kspnew->pc_side = snesksp->pc_side;
310: kspnew->rtol = snesksp->rtol;
311: kspnew->abstol = snesksp->abstol;
312: kspnew->max_it = snesksp->max_it;
313: KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
314: KSPGetPC(kspnew,&pcnew);
315: PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
316: PCSetOperators(kspnew->pc,Amat,Pmat);
317: PCFactorGetMatSolverType(snesksp->pc,&stype);
318: PCFactorSetMatSolverType(kspnew->pc,stype);
319: KSPDestroy(&snesksp);
320: snes->ksp = kspnew;
321: PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);
322: KSPSetFromOptions(kspnew);*/
323: return(0);
324: }
326: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
327: implemented in this algorithm. It basically identifies the active constraints and does
328: a linear solve on the other variables (those not associated with the active constraints). */
329: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
330: {
331: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
332: PetscErrorCode ierr;
333: PetscInt maxits,i,lits;
334: SNESLineSearchReason lssucceed;
335: PetscReal fnorm,gnorm,xnorm=0,ynorm;
336: Vec Y,X,F;
337: KSPConvergedReason kspreason;
338: KSP ksp;
339: PC pc;
342: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
343: SNESGetKSP(snes,&ksp);
344: KSPGetPC(ksp,&pc);
345: PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH);
347: snes->numFailures = 0;
348: snes->numLinearSolveFailures = 0;
349: snes->reason = SNES_CONVERGED_ITERATING;
351: maxits = snes->max_its; /* maximum number of iterations */
352: X = snes->vec_sol; /* solution vector */
353: F = snes->vec_func; /* residual vector */
354: Y = snes->work[0]; /* work vectors */
356: SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
357: SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
358: SNESLineSearchSetUp(snes->linesearch);
360: PetscObjectSAWsTakeAccess((PetscObject)snes);
361: snes->iter = 0;
362: snes->norm = 0.0;
363: PetscObjectSAWsGrantAccess((PetscObject)snes);
365: SNESVIProjectOntoBounds(snes,X);
366: SNESComputeFunction(snes,X,F);
367: SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
368: VecNorm(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
369: SNESCheckFunctionNorm(snes,fnorm);
370: PetscObjectSAWsTakeAccess((PetscObject)snes);
371: snes->norm = fnorm;
372: PetscObjectSAWsGrantAccess((PetscObject)snes);
373: SNESLogConvergenceHistory(snes,fnorm,0);
374: SNESMonitor(snes,0,fnorm);
376: /* test convergence */
377: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
378: if (snes->reason) return(0);
381: for (i=0; i<maxits; i++) {
383: IS IS_act; /* _act -> active set _inact -> inactive set */
384: IS IS_redact; /* redundant active set */
385: VecScatter scat_act,scat_inact;
386: PetscInt nis_act,nis_inact;
387: Vec Y_act,Y_inact,F_inact;
388: Mat jac_inact_inact,prejac_inact_inact;
389: PetscBool isequal;
391: /* Call general purpose update function */
392: if (snes->ops->update) {
393: (*snes->ops->update)(snes, snes->iter);
394: }
395: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
396: SNESCheckJacobianDomainerror(snes);
398: /* Create active and inactive index sets */
400: /*original
401: SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
402: */
403: SNESVIGetActiveSetIS(snes,X,F,&IS_act);
405: if (vi->checkredundancy) {
406: (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
407: if (IS_redact) {
408: ISSort(IS_redact);
409: ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);
410: ISDestroy(&IS_redact);
411: } else {
412: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
413: }
414: } else {
415: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
416: }
419: /* Create inactive set submatrix */
420: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
422: if (0) { /* Dead code (temporary developer hack) */
423: IS keptrows;
424: MatFindNonzeroRows(jac_inact_inact,&keptrows);
425: if (keptrows) {
426: PetscInt cnt,*nrows,k;
427: const PetscInt *krows,*inact;
428: PetscInt rstart;
430: MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);
431: MatDestroy(&jac_inact_inact);
432: ISDestroy(&IS_act);
434: ISGetLocalSize(keptrows,&cnt);
435: ISGetIndices(keptrows,&krows);
436: ISGetIndices(vi->IS_inact,&inact);
437: PetscMalloc1(cnt,&nrows);
438: for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
439: ISRestoreIndices(keptrows,&krows);
440: ISRestoreIndices(vi->IS_inact,&inact);
441: ISDestroy(&keptrows);
442: ISDestroy(&vi->IS_inact);
444: ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);
445: ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);
446: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
447: }
448: }
449: DMSetVI(snes->dm,vi->IS_inact);
450: /* remove later */
452: /*
453: VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
454: VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
455: VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
456: VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
457: ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
458: */
460: /* Get sizes of active and inactive sets */
461: ISGetLocalSize(IS_act,&nis_act);
462: ISGetLocalSize(vi->IS_inact,&nis_inact);
464: /* Create active and inactive set vectors */
465: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
466: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
467: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);
469: /* Create scatter contexts */
470: VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
471: VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);
473: /* Do a vec scatter to active and inactive set vectors */
474: VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
475: VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
477: VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
478: VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
480: VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
481: VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
483: /* Active set direction = 0 */
484: VecSet(Y_act,0);
485: if (snes->jacobian != snes->jacobian_pre) {
486: MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
487: } else prejac_inact_inact = jac_inact_inact;
489: ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);
490: if (!isequal) {
491: SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
492: PCFieldSplitRestrictIS(pc,vi->IS_inact);
493: }
495: /* ISView(vi->IS_inact,0); */
496: /* ISView(IS_act,0);*/
497: /* MatView(snes->jacobian_pre,0); */
501: KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);
502: KSPSetUp(snes->ksp);
503: {
504: PC pc;
505: PetscBool flg;
506: KSPGetPC(snes->ksp,&pc);
507: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
508: if (flg) {
509: KSP *subksps;
510: PCFieldSplitGetSubKSP(pc,NULL,&subksps);
511: KSPGetPC(subksps[0],&pc);
512: PetscFree(subksps);
513: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
514: if (flg) {
515: PetscInt n,N = 101*101,j,cnts[3] = {0,0,0};
516: const PetscInt *ii;
518: ISGetSize(vi->IS_inact,&n);
519: ISGetIndices(vi->IS_inact,&ii);
520: for (j=0; j<n; j++) {
521: if (ii[j] < N) cnts[0]++;
522: else if (ii[j] < 2*N) cnts[1]++;
523: else if (ii[j] < 3*N) cnts[2]++;
524: }
525: ISRestoreIndices(vi->IS_inact,&ii);
527: PCBJacobiSetTotalBlocks(pc,3,cnts);
528: }
529: }
530: }
532: KSPSolve(snes->ksp,F_inact,Y_inact);
533: VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
534: VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
535: VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
536: VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
538: VecDestroy(&F_inact);
539: VecDestroy(&Y_act);
540: VecDestroy(&Y_inact);
541: VecScatterDestroy(&scat_act);
542: VecScatterDestroy(&scat_inact);
543: ISDestroy(&IS_act);
544: if (!isequal) {
545: ISDestroy(&vi->IS_inact_prev);
546: ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);
547: }
548: ISDestroy(&vi->IS_inact);
549: MatDestroy(&jac_inact_inact);
550: if (snes->jacobian != snes->jacobian_pre) {
551: MatDestroy(&prejac_inact_inact);
552: }
554: KSPGetConvergedReason(snes->ksp,&kspreason);
555: if (kspreason < 0) {
556: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
557: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
558: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
559: break;
560: }
561: }
563: KSPGetIterationNumber(snes->ksp,&lits);
564: snes->linear_its += lits;
565: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
566: /*
567: if (snes->ops->precheck) {
568: PetscBool changed_y = PETSC_FALSE;
569: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
570: }
572: if (PetscLogPrintInfo) {
573: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
574: }
575: */
576: /* Compute a (scaled) negative update in the line search routine:
577: Y <- X - lambda*Y
578: and evaluate G = function(Y) (depends on the line search).
579: */
580: VecCopy(Y,snes->vec_sol_update);
581: ynorm = 1; gnorm = fnorm;
582: SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
583: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
584: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
585: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
586: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
587: if (snes->domainerror) {
588: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
589: DMDestroyVI(snes->dm);
590: return(0);
591: }
592: if (lssucceed) {
593: if (++snes->numFailures >= snes->maxFailures) {
594: PetscBool ismin;
595: snes->reason = SNES_DIVERGED_LINE_SEARCH;
596: SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);
597: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
598: break;
599: }
600: }
601: DMDestroyVI(snes->dm);
602: /* Update function and solution vectors */
603: fnorm = gnorm;
604: /* Monitor convergence */
605: PetscObjectSAWsTakeAccess((PetscObject)snes);
606: snes->iter = i+1;
607: snes->norm = fnorm;
608: snes->xnorm = xnorm;
609: snes->ynorm = ynorm;
610: PetscObjectSAWsGrantAccess((PetscObject)snes);
611: SNESLogConvergenceHistory(snes,snes->norm,lits);
612: SNESMonitor(snes,snes->iter,snes->norm);
613: /* Test for convergence, xnorm = || X || */
614: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
615: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
616: if (snes->reason) break;
617: }
618: /* 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 */
619: DMDestroyVI(snes->dm);
620: if (i == maxits) {
621: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
622: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
623: }
624: return(0);
625: }
627: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
628: {
629: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
633: vi->checkredundancy = func;
634: vi->ctxP = ctx;
635: return(0);
636: }
638: #if defined(PETSC_HAVE_MATLAB_ENGINE)
639: #include <engine.h>
640: #include <mex.h>
641: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
643: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx)
644: {
645: PetscErrorCode ierr;
646: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
647: int nlhs = 1, nrhs = 5;
648: mxArray *plhs[1], *prhs[5];
649: long long int l1 = 0, l2 = 0, ls = 0;
650: PetscInt *indices=NULL;
658: /* Create IS for reduced active set of size 0, its size and indices will
659: bet set by the Matlab function */
660: ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);
661: /* call Matlab function in ctx */
662: PetscMemcpy(&ls,&snes,sizeof(snes));
663: PetscMemcpy(&l1,&is_act,sizeof(is_act));
664: PetscMemcpy(&l2,is_redact,sizeof(is_act));
665: prhs[0] = mxCreateDoubleScalar((double)ls);
666: prhs[1] = mxCreateDoubleScalar((double)l1);
667: prhs[2] = mxCreateDoubleScalar((double)l2);
668: prhs[3] = mxCreateString(sctx->funcname);
669: prhs[4] = sctx->ctx;
670: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");
671: mxGetScalar(plhs[0]);
672: mxDestroyArray(prhs[0]);
673: mxDestroyArray(prhs[1]);
674: mxDestroyArray(prhs[2]);
675: mxDestroyArray(prhs[3]);
676: mxDestroyArray(plhs[0]);
677: return(0);
678: }
680: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
681: {
682: PetscErrorCode ierr;
683: SNESMatlabContext *sctx;
686: /* currently sctx is memory bleed */
687: PetscNew(&sctx);
688: PetscStrallocpy(func,&sctx->funcname);
689: sctx->ctx = mxDuplicateArray(ctx);
690: SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
691: return(0);
692: }
694: #endif
696: /* -------------------------------------------------------------------------- */
697: /*
698: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
699: of the SNESVI nonlinear solver.
701: Input Parameter:
702: . snes - the SNES context
704: Application Interface Routine: SNESSetUp()
706: Notes:
707: For basic use of the SNES solvers, the user need not explicitly call
708: SNESSetUp(), since these actions will automatically occur during
709: the call to SNESSolve().
710: */
711: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
712: {
713: PetscErrorCode ierr;
714: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
715: PetscInt *indices;
716: PetscInt i,n,rstart,rend;
717: SNESLineSearch linesearch;
720: SNESSetUp_VI(snes);
722: /* Set up previous active index set for the first snes solve
723: vi->IS_inact_prev = 0,1,2,....N */
725: VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
726: VecGetLocalSize(snes->vec_sol,&n);
727: PetscMalloc1(n,&indices);
728: for (i=0; i < n; i++) indices[i] = rstart + i;
729: ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);
731: /* set the line search functions */
732: if (!snes->linesearch) {
733: SNESGetLineSearch(snes, &linesearch);
734: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
735: }
736: return(0);
737: }
738: /* -------------------------------------------------------------------------- */
739: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
740: {
741: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
742: PetscErrorCode ierr;
745: SNESReset_VI(snes);
746: ISDestroy(&vi->IS_inact_prev);
747: return(0);
748: }
750: /* -------------------------------------------------------------------------- */
751: /*MC
752: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
754: Options Database:
755: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
756: - -snes_vi_monitor - prints the number of active constraints at each iteration.
758: Level: beginner
760: References:
761: . 1. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
762: Applications, Optimization Methods and Software, 21 (2006).
764: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
766: M*/
767: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
768: {
769: PetscErrorCode ierr;
770: SNES_VINEWTONRSLS *vi;
773: snes->ops->reset = SNESReset_VINEWTONRSLS;
774: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
775: snes->ops->solve = SNESSolve_VINEWTONRSLS;
776: snes->ops->destroy = SNESDestroy_VI;
777: snes->ops->setfromoptions = SNESSetFromOptions_VI;
778: snes->ops->view = NULL;
779: snes->ops->converged = SNESConvergedDefault_VI;
781: snes->usesksp = PETSC_TRUE;
782: snes->usesnpc = PETSC_FALSE;
784: snes->alwayscomputesfinalresidual = PETSC_TRUE;
786: PetscNewLog(snes,&vi);
787: snes->data = (void*)vi;
788: vi->checkredundancy = NULL;
790: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
791: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
792: return(0);
793: }