Actual source code: virs.c
petsc-3.14.6 2021-03-30
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 (*createinjection)(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: /*
65: 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: {
71: PetscContainer isnes;
72: DM_SNESVI *dmsnesvi1,*dmsnesvi2;
73: Mat interp;
76: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
77: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
78: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
79: PetscObjectQuery((PetscObject)dm2,"VI",(PetscObject*)&isnes);
80: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm2),PETSC_ERR_PLIB,"Composed VI data structure is missing");
81: PetscContainerGetPointer(isnes,(void**)&dmsnesvi2);
83: (*dmsnesvi1->createinterpolation)(dm1,dm2,&interp,NULL);
84: MatCreateSubMatrix(interp,dmsnesvi2->inactive,dmsnesvi1->inactive,MAT_INITIAL_MATRIX,mat);
85: MatDestroy(&interp);
86: *vec = NULL;
87: return(0);
88: }
90: static PetscErrorCode DMSetVI(DM,IS);
91: static PetscErrorCode DMDestroyVI(DM);
93: /*
94: DMCoarsen_SNESVI - Computes the regular coarsened DM then computes additional information about its inactive set
96: */
97: PetscErrorCode DMCoarsen_SNESVI(DM dm1,MPI_Comm comm,DM *dm2)
98: {
100: PetscContainer isnes;
101: DM_SNESVI *dmsnesvi1;
102: Vec finemarked,coarsemarked;
103: IS inactive;
104: Mat inject;
105: const PetscInt *index;
106: PetscInt n,k,cnt = 0,rstart,*coarseindex;
107: PetscScalar *marked;
110: PetscObjectQuery((PetscObject)dm1,"VI",(PetscObject*)&isnes);
111: if (!isnes) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_PLIB,"Composed VI data structure is missing");
112: PetscContainerGetPointer(isnes,(void**)&dmsnesvi1);
114: /* get the original coarsen */
115: (*dmsnesvi1->coarsen)(dm1,comm,dm2);
117: /* not sure why this extra reference is needed, but without the dm2 disappears too early */
118: /* Updating the KSPCreateVecs() to avoid using DMGetGlobalVector() when matrix is available removes the need for this reference? */
119: /* PetscObjectReference((PetscObject)*dm2);*/
121: /* need to set back global vectors in order to use the original injection */
122: DMClearGlobalVectors(dm1);
124: dm1->ops->createglobalvector = dmsnesvi1->createglobalvector;
126: DMCreateGlobalVector(dm1,&finemarked);
127: DMCreateGlobalVector(*dm2,&coarsemarked);
129: /*
130: fill finemarked with locations of inactive points
131: */
132: ISGetIndices(dmsnesvi1->inactive,&index);
133: ISGetLocalSize(dmsnesvi1->inactive,&n);
134: VecSet(finemarked,0.0);
135: for (k=0; k<n; k++) {
136: VecSetValue(finemarked,index[k],1.0,INSERT_VALUES);
137: }
138: VecAssemblyBegin(finemarked);
139: VecAssemblyEnd(finemarked);
141: DMCreateInjection(*dm2,dm1,&inject);
142: MatRestrict(inject,finemarked,coarsemarked);
143: MatDestroy(&inject);
145: /*
146: create index set list of coarse inactive points from coarsemarked
147: */
148: VecGetLocalSize(coarsemarked,&n);
149: VecGetOwnershipRange(coarsemarked,&rstart,NULL);
150: VecGetArray(coarsemarked,&marked);
151: for (k=0; k<n; k++) {
152: if (marked[k] != 0.0) cnt++;
153: }
154: PetscMalloc1(cnt,&coarseindex);
155: cnt = 0;
156: for (k=0; k<n; k++) {
157: if (marked[k] != 0.0) coarseindex[cnt++] = k + rstart;
158: }
159: VecRestoreArray(coarsemarked,&marked);
160: ISCreateGeneral(PetscObjectComm((PetscObject)coarsemarked),cnt,coarseindex,PETSC_OWN_POINTER,&inactive);
162: DMClearGlobalVectors(dm1);
164: dm1->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
166: DMSetVI(*dm2,inactive);
168: VecDestroy(&finemarked);
169: VecDestroy(&coarsemarked);
170: ISDestroy(&inactive);
171: return(0);
172: }
174: PetscErrorCode DMDestroy_SNESVI(DM_SNESVI *dmsnesvi)
175: {
179: /* reset the base methods in the DM object that were changed when the DM_SNESVI was reset */
180: dmsnesvi->dm->ops->createinterpolation = dmsnesvi->createinterpolation;
181: dmsnesvi->dm->ops->coarsen = dmsnesvi->coarsen;
182: dmsnesvi->dm->ops->createglobalvector = dmsnesvi->createglobalvector;
183: dmsnesvi->dm->ops->createinjection = dmsnesvi->createinjection;
184: dmsnesvi->dm->ops->hascreateinjection = dmsnesvi->hascreateinjection;
185: /* need to clear out this vectors because some of them may not have a reference to the DM
186: but they are counted as having references to the DM in DMDestroy() */
187: DMClearGlobalVectors(dmsnesvi->dm);
189: ISDestroy(&dmsnesvi->inactive);
190: PetscFree(dmsnesvi);
191: return(0);
192: }
194: /*
195: DMSetVI - Marks a DM as associated with a VI problem. This causes the interpolation/restriction operators to
196: be restricted to only those variables NOT associated with active constraints.
198: */
199: static PetscErrorCode DMSetVI(DM dm,IS inactive)
200: {
202: PetscContainer isnes;
203: DM_SNESVI *dmsnesvi;
206: if (!dm) return(0);
208: PetscObjectReference((PetscObject)inactive);
210: PetscObjectQuery((PetscObject)dm,"VI",(PetscObject*)&isnes);
211: if (!isnes) {
212: PetscContainerCreate(PetscObjectComm((PetscObject)dm),&isnes);
213: PetscContainerSetUserDestroy(isnes,(PetscErrorCode (*)(void*))DMDestroy_SNESVI);
214: PetscNew(&dmsnesvi);
215: PetscContainerSetPointer(isnes,(void*)dmsnesvi);
216: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)isnes);
217: PetscContainerDestroy(&isnes);
219: dmsnesvi->createinterpolation = dm->ops->createinterpolation;
220: dm->ops->createinterpolation = DMCreateInterpolation_SNESVI;
221: dmsnesvi->coarsen = dm->ops->coarsen;
222: dm->ops->coarsen = DMCoarsen_SNESVI;
223: dmsnesvi->createglobalvector = dm->ops->createglobalvector;
224: dm->ops->createglobalvector = DMCreateGlobalVector_SNESVI;
225: dmsnesvi->createinjection = dm->ops->createinjection;
226: dm->ops->createinjection = NULL;
227: dmsnesvi->hascreateinjection = dm->ops->hascreateinjection;
228: dm->ops->hascreateinjection = NULL;
229: } else {
230: PetscContainerGetPointer(isnes,(void**)&dmsnesvi);
231: ISDestroy(&dmsnesvi->inactive);
232: }
233: DMClearGlobalVectors(dm);
234: ISGetLocalSize(inactive,&dmsnesvi->n);
236: dmsnesvi->inactive = inactive;
237: dmsnesvi->dm = dm;
238: return(0);
239: }
241: /*
242: DMDestroyVI - Frees the DM_SNESVI object contained in the DM
243: - also resets the function pointers in the DM for createinterpolation() etc to use the original DM
244: */
245: static PetscErrorCode DMDestroyVI(DM dm)
246: {
250: if (!dm) return(0);
251: PetscObjectCompose((PetscObject)dm,"VI",(PetscObject)NULL);
252: return(0);
253: }
255: /* --------------------------------------------------------------------------------------------------------*/
258: PetscErrorCode SNESCreateIndexSets_VINEWTONRSLS(SNES snes,Vec X,Vec F,IS *ISact,IS *ISinact)
259: {
263: SNESVIGetActiveSetIS(snes,X,F,ISact);
264: ISComplement(*ISact,X->map->rstart,X->map->rend,ISinact);
265: return(0);
266: }
268: /* Create active and inactive set vectors. The local size of this vector is set and petsc computes the global size */
269: PetscErrorCode SNESCreateSubVectors_VINEWTONRSLS(SNES snes,PetscInt n,Vec *newv)
270: {
272: Vec v;
275: VecCreate(PetscObjectComm((PetscObject)snes),&v);
276: VecSetSizes(v,n,PETSC_DECIDE);
277: VecSetType(v,VECSTANDARD);
278: *newv = v;
279: return(0);
280: }
282: /* Resets the snes PC and KSP when the active set sizes change */
283: PetscErrorCode SNESVIResetPCandKSP(SNES snes,Mat Amat,Mat Pmat)
284: {
286: KSP snesksp;
289: SNESGetKSP(snes,&snesksp);
290: KSPReset(snesksp);
291: KSPResetFromOptions(snesksp);
293: /*
294: KSP kspnew;
295: PC pcnew;
296: MatSolverType stype;
299: KSPCreate(PetscObjectComm((PetscObject)snes),&kspnew);
300: kspnew->pc_side = snesksp->pc_side;
301: kspnew->rtol = snesksp->rtol;
302: kspnew->abstol = snesksp->abstol;
303: kspnew->max_it = snesksp->max_it;
304: KSPSetType(kspnew,((PetscObject)snesksp)->type_name);
305: KSPGetPC(kspnew,&pcnew);
306: PCSetType(kspnew->pc,((PetscObject)snesksp->pc)->type_name);
307: PCSetOperators(kspnew->pc,Amat,Pmat);
308: PCFactorGetMatSolverType(snesksp->pc,&stype);
309: PCFactorSetMatSolverType(kspnew->pc,stype);
310: KSPDestroy(&snesksp);
311: snes->ksp = kspnew;
312: PetscLogObjectParent((PetscObject)snes,(PetscObject)kspnew);
313: KSPSetFromOptions(kspnew);*/
314: return(0);
315: }
317: /* Variational Inequality solver using reduce space method. No semismooth algorithm is
318: implemented in this algorithm. It basically identifies the active constraints and does
319: a linear solve on the other variables (those not associated with the active constraints). */
320: PetscErrorCode SNESSolve_VINEWTONRSLS(SNES snes)
321: {
322: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
323: PetscErrorCode ierr;
324: PetscInt maxits,i,lits;
325: SNESLineSearchReason lssucceed;
326: PetscReal fnorm,gnorm,xnorm=0,ynorm;
327: Vec Y,X,F;
328: KSPConvergedReason kspreason;
329: KSP ksp;
330: PC pc;
333: /* Multigrid must use Galerkin for coarse grids with active set/reduced space methods; cannot rediscretize on coarser grids*/
334: SNESGetKSP(snes,&ksp);
335: KSPGetPC(ksp,&pc);
336: PCMGSetGalerkin(pc,PC_MG_GALERKIN_BOTH);
338: snes->numFailures = 0;
339: snes->numLinearSolveFailures = 0;
340: snes->reason = SNES_CONVERGED_ITERATING;
342: maxits = snes->max_its; /* maximum number of iterations */
343: X = snes->vec_sol; /* solution vector */
344: F = snes->vec_func; /* residual vector */
345: Y = snes->work[0]; /* work vectors */
347: SNESLineSearchSetVIFunctions(snes->linesearch, SNESVIProjectOntoBounds, SNESVIComputeInactiveSetFnorm);
348: SNESLineSearchSetVecs(snes->linesearch, X, NULL, NULL, NULL, NULL);
349: SNESLineSearchSetUp(snes->linesearch);
351: PetscObjectSAWsTakeAccess((PetscObject)snes);
352: snes->iter = 0;
353: snes->norm = 0.0;
354: PetscObjectSAWsGrantAccess((PetscObject)snes);
356: SNESVIProjectOntoBounds(snes,X);
357: SNESComputeFunction(snes,X,F);
358: SNESVIComputeInactiveSetFnorm(snes,F,X,&fnorm);
359: VecNorm(X,NORM_2,&xnorm); /* xnorm <- ||x|| */
360: SNESCheckFunctionNorm(snes,fnorm);
361: PetscObjectSAWsTakeAccess((PetscObject)snes);
362: snes->norm = fnorm;
363: PetscObjectSAWsGrantAccess((PetscObject)snes);
364: SNESLogConvergenceHistory(snes,fnorm,0);
365: SNESMonitor(snes,0,fnorm);
367: /* test convergence */
368: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
369: if (snes->reason) return(0);
372: for (i=0; i<maxits; i++) {
374: IS IS_act; /* _act -> active set _inact -> inactive set */
375: IS IS_redact; /* redundant active set */
376: VecScatter scat_act,scat_inact;
377: PetscInt nis_act,nis_inact;
378: Vec Y_act,Y_inact,F_inact;
379: Mat jac_inact_inact,prejac_inact_inact;
380: PetscBool isequal;
382: /* Call general purpose update function */
383: if (snes->ops->update) {
384: (*snes->ops->update)(snes, snes->iter);
385: }
386: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
387: SNESCheckJacobianDomainerror(snes);
389: /* Create active and inactive index sets */
391: /*original
392: SNESVICreateIndexSets_RS(snes,X,F,&IS_act,&vi->IS_inact);
393: */
394: SNESVIGetActiveSetIS(snes,X,F,&IS_act);
396: if (vi->checkredundancy) {
397: (*vi->checkredundancy)(snes,IS_act,&IS_redact,vi->ctxP);
398: if (IS_redact) {
399: ISSort(IS_redact);
400: ISComplement(IS_redact,X->map->rstart,X->map->rend,&vi->IS_inact);
401: ISDestroy(&IS_redact);
402: } else {
403: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
404: }
405: } else {
406: ISComplement(IS_act,X->map->rstart,X->map->rend,&vi->IS_inact);
407: }
410: /* Create inactive set submatrix */
411: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
413: if (0) { /* Dead code (temporary developer hack) */
414: IS keptrows;
415: MatFindNonzeroRows(jac_inact_inact,&keptrows);
416: if (keptrows) {
417: PetscInt cnt,*nrows,k;
418: const PetscInt *krows,*inact;
419: PetscInt rstart;
421: MatGetOwnershipRange(jac_inact_inact,&rstart,NULL);
422: MatDestroy(&jac_inact_inact);
423: ISDestroy(&IS_act);
425: ISGetLocalSize(keptrows,&cnt);
426: ISGetIndices(keptrows,&krows);
427: ISGetIndices(vi->IS_inact,&inact);
428: PetscMalloc1(cnt,&nrows);
429: for (k=0; k<cnt; k++) nrows[k] = inact[krows[k]-rstart];
430: ISRestoreIndices(keptrows,&krows);
431: ISRestoreIndices(vi->IS_inact,&inact);
432: ISDestroy(&keptrows);
433: ISDestroy(&vi->IS_inact);
435: ISCreateGeneral(PetscObjectComm((PetscObject)snes),cnt,nrows,PETSC_OWN_POINTER,&vi->IS_inact);
436: ISComplement(vi->IS_inact,F->map->rstart,F->map->rend,&IS_act);
437: MatCreateSubMatrix(snes->jacobian,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);
438: }
439: }
440: DMSetVI(snes->dm,vi->IS_inact);
441: /* remove later */
443: /*
444: VecView(vi->xu,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xu))->comm));
445: VecView(vi->xl,PETSC_VIEWER_BINARY_(((PetscObject)(vi->xl))->comm));
446: VecView(X,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)X)));
447: VecView(F,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)F)));
448: ISView(vi->IS_inact,PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)vi->IS_inact)));
449: */
451: /* Get sizes of active and inactive sets */
452: ISGetLocalSize(IS_act,&nis_act);
453: ISGetLocalSize(vi->IS_inact,&nis_inact);
455: /* Create active and inactive set vectors */
456: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&F_inact);
457: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_act,&Y_act);
458: SNESCreateSubVectors_VINEWTONRSLS(snes,nis_inact,&Y_inact);
460: /* Create scatter contexts */
461: VecScatterCreate(Y,IS_act,Y_act,NULL,&scat_act);
462: VecScatterCreate(Y,vi->IS_inact,Y_inact,NULL,&scat_inact);
464: /* Do a vec scatter to active and inactive set vectors */
465: VecScatterBegin(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
466: VecScatterEnd(scat_inact,F,F_inact,INSERT_VALUES,SCATTER_FORWARD);
468: VecScatterBegin(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
469: VecScatterEnd(scat_act,Y,Y_act,INSERT_VALUES,SCATTER_FORWARD);
471: VecScatterBegin(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
472: VecScatterEnd(scat_inact,Y,Y_inact,INSERT_VALUES,SCATTER_FORWARD);
474: /* Active set direction = 0 */
475: VecSet(Y_act,0);
476: if (snes->jacobian != snes->jacobian_pre) {
477: MatCreateSubMatrix(snes->jacobian_pre,vi->IS_inact,vi->IS_inact,MAT_INITIAL_MATRIX,&prejac_inact_inact);
478: } else prejac_inact_inact = jac_inact_inact;
480: ISEqual(vi->IS_inact_prev,vi->IS_inact,&isequal);
481: if (!isequal) {
482: SNESVIResetPCandKSP(snes,jac_inact_inact,prejac_inact_inact);
483: PCFieldSplitRestrictIS(pc,vi->IS_inact);
484: }
486: /* ISView(vi->IS_inact,0); */
487: /* ISView(IS_act,0);*/
488: /* MatView(snes->jacobian_pre,0); */
492: KSPSetOperators(snes->ksp,jac_inact_inact,prejac_inact_inact);
493: KSPSetUp(snes->ksp);
494: {
495: PC pc;
496: PetscBool flg;
497: KSPGetPC(snes->ksp,&pc);
498: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&flg);
499: if (flg) {
500: KSP *subksps;
501: PCFieldSplitGetSubKSP(pc,NULL,&subksps);
502: KSPGetPC(subksps[0],&pc);
503: PetscFree(subksps);
504: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
505: if (flg) {
506: PetscInt n,N = 101*101,j,cnts[3] = {0,0,0};
507: const PetscInt *ii;
509: ISGetSize(vi->IS_inact,&n);
510: ISGetIndices(vi->IS_inact,&ii);
511: for (j=0; j<n; j++) {
512: if (ii[j] < N) cnts[0]++;
513: else if (ii[j] < 2*N) cnts[1]++;
514: else if (ii[j] < 3*N) cnts[2]++;
515: }
516: ISRestoreIndices(vi->IS_inact,&ii);
518: PCBJacobiSetTotalBlocks(pc,3,cnts);
519: }
520: }
521: }
523: KSPSolve(snes->ksp,F_inact,Y_inact);
524: VecScatterBegin(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
525: VecScatterEnd(scat_act,Y_act,Y,INSERT_VALUES,SCATTER_REVERSE);
526: VecScatterBegin(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
527: VecScatterEnd(scat_inact,Y_inact,Y,INSERT_VALUES,SCATTER_REVERSE);
529: VecDestroy(&F_inact);
530: VecDestroy(&Y_act);
531: VecDestroy(&Y_inact);
532: VecScatterDestroy(&scat_act);
533: VecScatterDestroy(&scat_inact);
534: ISDestroy(&IS_act);
535: if (!isequal) {
536: ISDestroy(&vi->IS_inact_prev);
537: ISDuplicate(vi->IS_inact,&vi->IS_inact_prev);
538: }
539: ISDestroy(&vi->IS_inact);
540: MatDestroy(&jac_inact_inact);
541: if (snes->jacobian != snes->jacobian_pre) {
542: MatDestroy(&prejac_inact_inact);
543: }
545: KSPGetConvergedReason(snes->ksp,&kspreason);
546: if (kspreason < 0) {
547: if (++snes->numLinearSolveFailures >= snes->maxLinearSolveFailures) {
548: PetscInfo2(snes,"iter=%D, number linear solve failures %D greater than current SNES allowed, stopping solve\n",snes->iter,snes->numLinearSolveFailures);
549: snes->reason = SNES_DIVERGED_LINEAR_SOLVE;
550: break;
551: }
552: }
554: KSPGetIterationNumber(snes->ksp,&lits);
555: snes->linear_its += lits;
556: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
557: /*
558: if (snes->ops->precheck) {
559: PetscBool changed_y = PETSC_FALSE;
560: (*snes->ops->precheck)(snes,X,Y,snes->precheck,&changed_y);
561: }
563: if (PetscLogPrintInfo) {
564: SNESVICheckResidual_Private(snes,snes->jacobian,F,Y,G,W);
565: }
566: */
567: /* Compute a (scaled) negative update in the line search routine:
568: Y <- X - lambda*Y
569: and evaluate G = function(Y) (depends on the line search).
570: */
571: VecCopy(Y,snes->vec_sol_update);
572: ynorm = 1; gnorm = fnorm;
573: SNESLineSearchApply(snes->linesearch, X, F, &gnorm, Y);
574: SNESLineSearchGetReason(snes->linesearch, &lssucceed);
575: SNESLineSearchGetNorms(snes->linesearch, &xnorm, &gnorm, &ynorm);
576: PetscInfo4(snes,"fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n",(double)fnorm,(double)gnorm,(double)ynorm,(int)lssucceed);
577: if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
578: if (snes->domainerror) {
579: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
580: DMDestroyVI(snes->dm);
581: return(0);
582: }
583: if (lssucceed) {
584: if (++snes->numFailures >= snes->maxFailures) {
585: PetscBool ismin;
586: snes->reason = SNES_DIVERGED_LINE_SEARCH;
587: SNESVICheckLocalMin_Private(snes,snes->jacobian,F,X,gnorm,&ismin);
588: if (ismin) snes->reason = SNES_DIVERGED_LOCAL_MIN;
589: break;
590: }
591: }
592: DMDestroyVI(snes->dm);
593: /* Update function and solution vectors */
594: fnorm = gnorm;
595: /* Monitor convergence */
596: PetscObjectSAWsTakeAccess((PetscObject)snes);
597: snes->iter = i+1;
598: snes->norm = fnorm;
599: snes->xnorm = xnorm;
600: snes->ynorm = ynorm;
601: PetscObjectSAWsGrantAccess((PetscObject)snes);
602: SNESLogConvergenceHistory(snes,snes->norm,lits);
603: SNESMonitor(snes,snes->iter,snes->norm);
604: /* Test for convergence, xnorm = || X || */
605: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
606: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
607: if (snes->reason) break;
608: }
609: /* 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 */
610: DMDestroyVI(snes->dm);
611: if (i == maxits) {
612: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
613: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
614: }
615: return(0);
616: }
618: PetscErrorCode SNESVISetRedundancyCheck(SNES snes,PetscErrorCode (*func)(SNES,IS,IS*,void*),void *ctx)
619: {
620: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*)snes->data;
624: vi->checkredundancy = func;
625: vi->ctxP = ctx;
626: return(0);
627: }
629: #if defined(PETSC_HAVE_MATLAB_ENGINE)
630: #include <engine.h>
631: #include <mex.h>
632: typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext;
634: PetscErrorCode SNESVIRedundancyCheck_Matlab(SNES snes,IS is_act,IS *is_redact,void *ctx)
635: {
636: PetscErrorCode ierr;
637: SNESMatlabContext *sctx = (SNESMatlabContext*)ctx;
638: int nlhs = 1, nrhs = 5;
639: mxArray *plhs[1], *prhs[5];
640: long long int l1 = 0, l2 = 0, ls = 0;
641: PetscInt *indices=NULL;
649: /* Create IS for reduced active set of size 0, its size and indices will
650: bet set by the Matlab function */
651: ISCreateGeneral(PetscObjectComm((PetscObject)snes),0,indices,PETSC_OWN_POINTER,is_redact);
652: /* call Matlab function in ctx */
653: PetscArraycpy(&ls,&snes,1);
654: PetscArraycpy(&l1,&is_act,1);
655: PetscArraycpy(&l2,is_redact,1);
656: prhs[0] = mxCreateDoubleScalar((double)ls);
657: prhs[1] = mxCreateDoubleScalar((double)l1);
658: prhs[2] = mxCreateDoubleScalar((double)l2);
659: prhs[3] = mxCreateString(sctx->funcname);
660: prhs[4] = sctx->ctx;
661: mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESVIRedundancyCheckInternal");
662: mxGetScalar(plhs[0]);
663: mxDestroyArray(prhs[0]);
664: mxDestroyArray(prhs[1]);
665: mxDestroyArray(prhs[2]);
666: mxDestroyArray(prhs[3]);
667: mxDestroyArray(plhs[0]);
668: return(0);
669: }
671: PetscErrorCode SNESVISetRedundancyCheckMatlab(SNES snes,const char *func,mxArray *ctx)
672: {
673: PetscErrorCode ierr;
674: SNESMatlabContext *sctx;
677: /* currently sctx is memory bleed */
678: PetscNew(&sctx);
679: PetscStrallocpy(func,&sctx->funcname);
680: sctx->ctx = mxDuplicateArray(ctx);
681: SNESVISetRedundancyCheck(snes,SNESVIRedundancyCheck_Matlab,sctx);
682: return(0);
683: }
685: #endif
687: /* -------------------------------------------------------------------------- */
688: /*
689: SNESSetUp_VINEWTONRSLS - Sets up the internal data structures for the later use
690: of the SNESVI nonlinear solver.
692: Input Parameter:
693: . snes - the SNES context
695: Application Interface Routine: SNESSetUp()
697: Notes:
698: For basic use of the SNES solvers, the user need not explicitly call
699: SNESSetUp(), since these actions will automatically occur during
700: the call to SNESSolve().
701: */
702: PetscErrorCode SNESSetUp_VINEWTONRSLS(SNES snes)
703: {
704: PetscErrorCode ierr;
705: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
706: PetscInt *indices;
707: PetscInt i,n,rstart,rend;
708: SNESLineSearch linesearch;
711: SNESSetUp_VI(snes);
713: /* Set up previous active index set for the first snes solve
714: vi->IS_inact_prev = 0,1,2,....N */
716: VecGetOwnershipRange(snes->vec_sol,&rstart,&rend);
717: VecGetLocalSize(snes->vec_sol,&n);
718: PetscMalloc1(n,&indices);
719: for (i=0; i < n; i++) indices[i] = rstart + i;
720: ISCreateGeneral(PetscObjectComm((PetscObject)snes),n,indices,PETSC_OWN_POINTER,&vi->IS_inact_prev);
722: /* set the line search functions */
723: if (!snes->linesearch) {
724: SNESGetLineSearch(snes, &linesearch);
725: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
726: }
727: return(0);
728: }
729: /* -------------------------------------------------------------------------- */
730: PetscErrorCode SNESReset_VINEWTONRSLS(SNES snes)
731: {
732: SNES_VINEWTONRSLS *vi = (SNES_VINEWTONRSLS*) snes->data;
733: PetscErrorCode ierr;
736: SNESReset_VI(snes);
737: ISDestroy(&vi->IS_inact_prev);
738: return(0);
739: }
741: /* -------------------------------------------------------------------------- */
742: /*MC
743: SNESVINEWTONRSLS - Reduced space active set solvers for variational inequalities based on Newton's method
745: Options Database:
746: + -snes_type <vinewtonssls,vinewtonrsls> - a semi-smooth solver, a reduced space active set method
747: - -snes_vi_monitor - prints the number of active constraints at each iteration.
749: Level: beginner
751: References:
752: . 1. - T. S. Munson, and S. Benson. Flexible Complementarity Solvers for Large Scale
753: Applications, Optimization Methods and Software, 21 (2006).
755: .seealso: SNESVISetVariableBounds(), SNESVISetComputeVariableBounds(), SNESCreate(), SNES, SNESSetType(), SNESVINEWTONSSLS, SNESNEWTONTR, SNESLineSearchSetType(),SNESLineSearchSetPostCheck(), SNESLineSearchSetPreCheck()
757: M*/
758: PETSC_EXTERN PetscErrorCode SNESCreate_VINEWTONRSLS(SNES snes)
759: {
760: PetscErrorCode ierr;
761: SNES_VINEWTONRSLS *vi;
762: SNESLineSearch linesearch;
765: snes->ops->reset = SNESReset_VINEWTONRSLS;
766: snes->ops->setup = SNESSetUp_VINEWTONRSLS;
767: snes->ops->solve = SNESSolve_VINEWTONRSLS;
768: snes->ops->destroy = SNESDestroy_VI;
769: snes->ops->setfromoptions = SNESSetFromOptions_VI;
770: snes->ops->view = NULL;
771: snes->ops->converged = SNESConvergedDefault_VI;
773: snes->usesksp = PETSC_TRUE;
774: snes->usesnpc = PETSC_FALSE;
776: SNESGetLineSearch(snes, &linesearch);
777: if (!((PetscObject)linesearch)->type_name) {
778: SNESLineSearchSetType(linesearch, SNESLINESEARCHBT);
779: }
780: SNESLineSearchBTSetAlpha(linesearch, 0.0);
782: snes->alwayscomputesfinalresidual = PETSC_TRUE;
784: PetscNewLog(snes,&vi);
785: snes->data = (void*)vi;
786: vi->checkredundancy = NULL;
788: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetVariableBounds_C",SNESVISetVariableBounds_VI);
789: PetscObjectComposeFunction((PetscObject)snes,"SNESVISetComputeVariableBounds_C",SNESVISetComputeVariableBounds_VI);
790: return(0);
791: }