Actual source code: damgsnes.c
1: #define PETSCSNES_DLL
2:
3: #include petscda.h
4: #include src/dm/da/daimpl.h
5: /* It appears that preprocessor directives are not respected by bfort */
6: #ifdef PETSC_HAVE_SIEVE
7: #include petscmesh.h
8: #endif
9: #include petscmg.h
10: #include petscdmmg.h
12: #if defined(PETSC_HAVE_ADIC)
19: #endif
22: EXTERN PetscErrorCode NLFCreate_DAAD(NLF*);
23: EXTERN PetscErrorCode NLFDAADSetDA_DAAD(NLF,DA);
24: EXTERN PetscErrorCode NLFDAADSetCtx_DAAD(NLF,void*);
25: EXTERN PetscErrorCode NLFDAADSetResidual_DAAD(NLF,Vec);
26: EXTERN PetscErrorCode NLFDAADSetNewtonIterations_DAAD(NLF,PetscInt);
29: /*
30: period of -1 indicates update only on zeroth iteration of SNES
31: */
32: #define ShouldUpdate(l,it) (((dmmg[l-1]->updatejacobianperiod == -1) && (it == 0)) || \
33: ((dmmg[l-1]->updatejacobianperiod > 0) && !(it % dmmg[l-1]->updatejacobianperiod)))
34: /*
35: Evaluates the Jacobian on all of the grids. It is used by DMMG to provide the
36: ComputeJacobian() function that SNESSetJacobian() requires.
37: */
40: PetscErrorCode DMMGComputeJacobian_Multigrid(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
41: {
42: DMMG *dmmg = (DMMG*)ptr;
44: PetscInt i,nlevels = dmmg[0]->nlevels,it;
45: KSP ksp,lksp;
46: PC pc;
47: PetscTruth ismg,galerkin = PETSC_FALSE;
48: Vec W;
49: MatStructure flg;
52: if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as user context which should contain DMMG");
53: SNESGetIterationNumber(snes,&it);
55: /* compute Jacobian on finest grid */
56: if (dmmg[nlevels-1]->updatejacobian && ShouldUpdate(nlevels,it)) {
57: (*DMMGGetFine(dmmg)->computejacobian)(snes,X,J,B,flag,DMMGGetFine(dmmg));
58: } else {
59: PetscInfo3(0,"Skipping Jacobian, SNES iteration %D frequence %D level %D\n",it,dmmg[nlevels-1]->updatejacobianperiod,nlevels-1);
60: *flag = SAME_PRECONDITIONER;
61: }
62: MatSNESMFSetBase(DMMGGetFine(dmmg)->J,X);
64: /* create coarser grid Jacobians for preconditioner if multigrid is the preconditioner */
65: SNESGetKSP(snes,&ksp);
66: KSPGetPC(ksp,&pc);
67: PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
68: if (ismg) {
69: PCMGGetGalerkin(pc,&galerkin);
70: }
72: if (!galerkin) {
73: for (i=nlevels-1; i>0; i--) {
74: if (!dmmg[i-1]->w) {
75: VecDuplicate(dmmg[i-1]->x,&dmmg[i-1]->w);
76: }
77: W = dmmg[i-1]->w;
78: /* restrict X to coarser grid */
79: MatRestrict(dmmg[i]->R,X,W);
80: X = W;
81: /* scale to "natural" scaling for that grid */
82: VecPointwiseMult(X,X,dmmg[i]->Rscale);
83: /* tell the base vector for matrix free multiplies */
84: MatSNESMFSetBase(dmmg[i-1]->J,X);
85: /* compute Jacobian on coarse grid */
86: if (dmmg[i-1]->updatejacobian && ShouldUpdate(i,it)) {
87: (*dmmg[i-1]->computejacobian)(snes,X,&dmmg[i-1]->J,&dmmg[i-1]->B,&flg,dmmg[i-1]);
88: flg = SAME_NONZERO_PATTERN;
89: } else {
90: PetscInfo3(0,"Skipping Jacobian, SNES iteration %D frequence %D level %D\n",it,dmmg[i-1]->updatejacobianperiod,i-1);
91: flg = SAME_PRECONDITIONER;
92: }
93: if (ismg) {
94: PCMGGetSmoother(pc,i-1,&lksp);
95: KSPSetOperators(lksp,dmmg[i-1]->J,dmmg[i-1]->B,flg);
96: }
97: }
98: }
99: return(0);
100: }
102: /* ---------------------------------------------------------------------------*/
107: /*
108: DMMGFormFunction - This is a universal global FormFunction used by the DMMG code
109: when the user provides a local function.
111: Input Parameters:
112: + snes - the SNES context
113: . X - input vector
114: - ptr - optional user-defined context, as set by SNESSetFunction()
116: Output Parameter:
117: . F - function vector
119: */
120: PetscErrorCode DMMGFormFunction(SNES snes,Vec X,Vec F,void *ptr)
121: {
122: DMMG dmmg = (DMMG)ptr;
124: Vec localX;
125: DA da = (DA)dmmg->dm;
128: DAGetLocalVector(da,&localX);
129: /*
130: Scatter ghost points to local vector, using the 2-step process
131: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
132: */
133: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
134: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
135: DAFormFunction1(da,localX,F,dmmg->user);
136: DARestoreLocalVector(da,&localX);
137: return(0);
138: }
142: PetscErrorCode DMMGFormFunctionGhost(SNES snes,Vec X,Vec F,void *ptr)
143: {
144: DMMG dmmg = (DMMG)ptr;
146: Vec localX, localF;
147: DA da = (DA)dmmg->dm;
150: DAGetLocalVector(da,&localX);
151: DAGetLocalVector(da,&localF);
152: /*
153: Scatter ghost points to local vector, using the 2-step process
154: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
155: */
156: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
157: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
158: VecSet(F, 0.0);
159: VecSet(localF, 0.0);
160: DAFormFunction1(da,localX,localF,dmmg->user);
161: DALocalToGlobalBegin(da,localF,F);
162: DALocalToGlobalEnd(da,localF,F);
163: DARestoreLocalVector(da,&localX);
164: DARestoreLocalVector(da,&localF);
165: return(0);
166: }
168: #ifdef PETSC_HAVE_SIEVE
171: /*
172: DMMGFormFunctionMesh - This is a universal global FormFunction used by the DMMG code
173: when the user provides a local function.
175: Input Parameters:
176: + snes - the SNES context
177: . X - input vector
178: - ptr - This is the DMMG object
180: Output Parameter:
181: . F - function vector
183: */
184: PetscErrorCode DMMGFormFunctionMesh(SNES snes, Vec X, Vec F, void *ptr)
185: {
186: DMMG dmmg = (DMMG) ptr;
187: Mesh mesh = (Mesh) dmmg->dm;
188: SectionReal sectionX, section;
189: Vec localVec;
190: VecScatter scatter;
194: MeshGetSectionReal(mesh, "default", §ion);
195: SectionRealDuplicate(section, §ionX);
196: /*
197: Scatter ghost points to local vector, using the 2-step process
198: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
199: */
200: MeshGetGlobalScatter(mesh, &scatter);
201: MeshCreateLocalVector(mesh, sectionX, &localVec);
202: VecScatterBegin(X, localVec, INSERT_VALUES, SCATTER_REVERSE, scatter);
203: VecScatterEnd(X, localVec, INSERT_VALUES, SCATTER_REVERSE, scatter);
204: VecDestroy(localVec);
205: MeshFormFunction(mesh, sectionX, section, dmmg->user);
206: MeshCreateLocalVector(mesh, section, &localVec);
207: VecScatterBegin(localVec, F, INSERT_VALUES, SCATTER_FORWARD, scatter);
208: VecScatterEnd(localVec, F, INSERT_VALUES, SCATTER_FORWARD, scatter);
209: VecDestroy(localVec);
210: SectionRealDestroy(sectionX);
211: SectionRealDestroy(section);
212: return(0);
213: }
217: /*
218: DMMGComputeJacobianMesh - This is a universal global FormJacobian used by the DMMG code
219: when the user provides a local function.
221: Input Parameters:
222: + snes - the SNES context
223: . X - input vector
224: - ptr - This is the DMMG object
226: Output Parameter:
227: . F - function vector
229: */
230: PetscErrorCode DMMGComputeJacobianMesh(SNES snes, Vec X, Mat *J, Mat *B, MatStructure *flag, void *ptr)
231: {
232: DMMG dmmg = (DMMG) ptr;
233: Mesh mesh = (Mesh) dmmg->dm;
234: SectionReal sectionX;
235: Vec localVec;
236: VecScatter scatter;
240: MeshGetSectionReal(mesh, "default", §ionX);
241: /*
242: Scatter ghost points to local vector, using the 2-step process
243: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
244: */
245: MeshGetGlobalScatter(mesh, &scatter);
246: MeshCreateLocalVector(mesh, sectionX, &localVec);
247: VecScatterBegin(X, localVec, INSERT_VALUES, SCATTER_REVERSE, scatter);
248: VecScatterEnd(X, localVec, INSERT_VALUES, SCATTER_REVERSE, scatter);
249: VecDestroy(localVec);
250: MeshFormJacobian(mesh, sectionX, *B, dmmg->user);
251: /* Assemble true Jacobian; if it is different */
252: if (*J != *B) {
253: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
254: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
255: }
256: MatSetOption(*B, MAT_NEW_NONZERO_LOCATION_ERR);
257: *flag = SAME_NONZERO_PATTERN;
258: SectionRealDestroy(sectionX);
259: return(0);
260: }
261: #endif
265: /*
266: DMMGFormFunctionFD - This is a universal global FormFunction used by the DMMG code
267: when the user provides a local function used to compute the Jacobian via FD.
269: Input Parameters:
270: + snes - the SNES context
271: . X - input vector
272: - ptr - optional user-defined context, as set by SNESSetFunction()
274: Output Parameter:
275: . F - function vector
277: */
278: PetscErrorCode DMMGFormFunctionFD(SNES snes,Vec X,Vec F,void *ptr)
279: {
280: DMMG dmmg = (DMMG)ptr;
282: Vec localX;
283: DA da = (DA)dmmg->dm;
284: PetscInt N,n;
285:
287: /* determine whether X=localX */
288: DAGetLocalVector(da,&localX);
289: VecGetSize(X,&N);
290: VecGetSize(localX,&n);
292: if (n != N){ /* X != localX */
293: /* Scatter ghost points to local vector, using the 2-step process
294: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
295: */
296: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
297: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
298: } else {
299: DARestoreLocalVector(da,&localX);
300: localX = X;
301: }
302: DAFormFunction(da,dmmg->lfj,localX,F,dmmg->user);
303: if (n != N){
304: DARestoreLocalVector(da,&localX);
305: }
306: return(0);
307: }
311: /*@C
312: SNESDAFormFunction - This is a universal function evaluation routine that
313: may be used with SNESSetFunction() as long as the user context has a DA
314: as its first record and the user has called DASetLocalFunction().
316: Collective on SNES
318: Input Parameters:
319: + snes - the SNES context
320: . X - input vector
321: . F - function vector
322: - ptr - pointer to a structure that must have a DA as its first entry. For example this
323: could be a DMMG, this ptr must have been passed into SNESDAFormFunction() as the context
325: Level: intermediate
327: .seealso: DASetLocalFunction(), DASetLocalJacobian(), DASetLocalAdicFunction(), DASetLocalAdicMFFunction(),
328: SNESSetFunction(), SNESSetJacobian()
330: @*/
331: PetscErrorCode SNESDAFormFunction(SNES snes,Vec X,Vec F,void *ptr)
332: {
334: Vec localX;
335: DA da = *(DA*)ptr;
336: PetscInt N,n;
337:
339: if (!da) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Looks like you called SNESSetFromFuntion(snes,SNESDAFormFunction,) without the DA context");
341: /* determine whether X=localX */
342: DAGetLocalVector(da,&localX);
343: VecGetSize(X,&N);
344: VecGetSize(localX,&n);
345:
346:
347: if (n != N){ /* X != localX */
348: /* Scatter ghost points to local vector, using the 2-step process
349: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
350: */
351: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
352: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
353: } else {
354: DARestoreLocalVector(da,&localX);
355: localX = X;
356: }
357: DAFormFunction1(da,localX,F,ptr);
358: if (n != N){
359: if (PetscExceptionValue(ierr)) {
360: PetscErrorCode pDARestoreLocalVector(da,&localX);CHKERRQ(pierr);
361: }
362:
363: DARestoreLocalVector(da,&localX);
364: }
365: return(0);
366: }
368: /* ------------------------------------------------------------------------------*/
369: #include include/private/matimpl.h
372: PetscErrorCode DMMGComputeJacobianWithFD(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
373: {
375: DMMG dmmg = (DMMG)ctx;
376: MatFDColoring color = (MatFDColoring)dmmg->fdcoloring;
377:
379: if (color->ctype == IS_COLORING_GHOSTED){
380: DA da=(DA)dmmg->dm;
381: Vec x1_loc;
382: DAGetLocalVector(da,&x1_loc);
383: DAGlobalToLocalBegin(da,x1,INSERT_VALUES,x1_loc);
384: DAGlobalToLocalEnd(da,x1,INSERT_VALUES,x1_loc);
385: SNESDefaultComputeJacobianColor(snes,x1_loc,J,B,flag,dmmg->fdcoloring);
386: DARestoreLocalVector(da,&x1_loc);
387: } else {
388: SNESDefaultComputeJacobianColor(snes,x1,J,B,flag,dmmg->fdcoloring);
389: }
390: return(0);
391: }
395: PetscErrorCode DMMGComputeJacobianWithMF(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
396: {
398:
400: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
401: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
402: return(0);
403: }
407: /*
408: DMMGComputeJacobian - Evaluates the Jacobian when the user has provided
409: a local function evaluation routine.
410: */
411: PetscErrorCode DMMGComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
412: {
413: DMMG dmmg = (DMMG) ptr;
415: Vec localX;
416: DA da = (DA) dmmg->dm;
419: DAGetLocalVector(da,&localX);
420: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
421: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
422: DAComputeJacobian1(da,localX,*B,dmmg->user);
423: DARestoreLocalVector(da,&localX);
424: /* Assemble true Jacobian; if it is different */
425: if (*J != *B) {
426: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
427: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
428: }
429: MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
430: *flag = SAME_NONZERO_PATTERN;
431: return(0);
432: }
436: /*
437: SNESDAComputeJacobianWithAdifor - This is a universal Jacobian evaluation routine
438: that may be used with SNESSetJacobian() from Fortran as long as the user context has
439: a DA as its first record and DASetLocalAdiforFunction() has been called.
441: Collective on SNES
443: Input Parameters:
444: + snes - the SNES context
445: . X - input vector
446: . J - Jacobian
447: . B - Jacobian used in preconditioner (usally same as J)
448: . flag - indicates if the matrix changed its structure
449: - ptr - optional user-defined context, as set by SNESSetFunction()
451: Level: intermediate
453: .seealso: DASetLocalFunction(), DASetLocalAdicFunction(), SNESSetFunction(), SNESSetJacobian()
455: */
456: PetscErrorCode SNESDAComputeJacobianWithAdifor(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
457: {
458: DA da = *(DA*) ptr;
460: Vec localX;
463: DAGetLocalVector(da,&localX);
464: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
465: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
466: DAComputeJacobian1WithAdifor(da,localX,*B,ptr);
467: DARestoreLocalVector(da,&localX);
468: /* Assemble true Jacobian; if it is different */
469: if (*J != *B) {
470: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
471: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
472: }
473: MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
474: *flag = SAME_NONZERO_PATTERN;
475: return(0);
476: }
480: /*
481: SNESDAComputeJacobian - This is a universal Jacobian evaluation routine for a
482: locally provided Jacobian.
484: Collective on SNES
486: Input Parameters:
487: + snes - the SNES context
488: . X - input vector
489: . J - Jacobian
490: . B - Jacobian used in preconditioner (usally same as J)
491: . flag - indicates if the matrix changed its structure
492: - ptr - optional user-defined context, as set by SNESSetFunction()
494: Level: intermediate
496: .seealso: DASetLocalFunction(), DASetLocalJacobian(), SNESSetFunction(), SNESSetJacobian()
498: */
499: PetscErrorCode SNESDAComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
500: {
501: DA da = *(DA*) ptr;
503: Vec localX;
506: DAGetLocalVector(da,&localX);
507: DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
508: DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
509: DAComputeJacobian1(da,localX,*B,ptr);
510: DARestoreLocalVector(da,&localX);
511: /* Assemble true Jacobian; if it is different */
512: if (*J != *B) {
513: MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
514: MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
515: }
516: MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
517: *flag = SAME_NONZERO_PATTERN;
518: return(0);
519: }
523: PetscErrorCode DMMGSolveSNES(DMMG *dmmg,PetscInt level)
524: {
526: PetscInt nlevels = dmmg[0]->nlevels;
529: dmmg[0]->nlevels = level+1;
530: SNESSolve(dmmg[level]->snes,PETSC_NULL,dmmg[level]->x);
531: dmmg[0]->nlevels = nlevels;
532: return(0);
533: }
535: /* ===========================================================================================================*/
539: /*@C
540: DMMGSetSNES - Sets the nonlinear function that defines the nonlinear set of equations
541: to be solved using the grid hierarchy.
543: Collective on DMMG
545: Input Parameter:
546: + dmmg - the context
547: . function - the function that defines the nonlinear system
548: - jacobian - optional function to compute Jacobian
550: Options Database Keys:
551: + -dmmg_snes_monitor
552: . -dmmg_jacobian_fd
553: . -dmmg_jacobian_ad
554: . -dmmg_jacobian_mf_fd_operator
555: . -dmmg_jacobian_mf_fd
556: . -dmmg_jacobian_mf_ad_operator
557: . -dmmg_jacobian_mf_ad
558: . -dmmg_iscoloring_type
559: - -dmmg_jacobian_period <p> - Indicates how often in the SNES solve the Jacobian is recomputed (on all levels)
560: as suggested by Florin Dobrian if p is -1 then Jacobian is computed only on first
561: SNES iteration (i.e. -1 is equivalent to infinity)
563: Level: advanced
565: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNESLocal()
567: @*/
568: PetscErrorCode DMMGSetSNES(DMMG *dmmg,PetscErrorCode (*function)(SNES,Vec,Vec,void*),PetscErrorCode (*jacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*))
569: {
570: PetscErrorCode ierr;
571: PetscInt i,nlevels = dmmg[0]->nlevels,period = 1,isctype=0;
572: PetscTruth snesmonitor,mffdoperator,mffd,fdjacobian;
573: #if defined(PETSC_HAVE_ADIC)
574: PetscTruth mfadoperator,mfad,adjacobian;
575: #endif
576: PetscViewerASCIIMonitor ascii;
577: MPI_Comm comm;
578: PetscMPIInt size;
579: ISColoringType ctype;
582: if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");
583: if (!jacobian) jacobian = DMMGComputeJacobianWithFD;
585: PetscOptionsBegin(dmmg[0]->comm,PETSC_NULL,"DMMG Options","SNES");
586: PetscOptionsName("-dmmg_snes_monitor","Monitor nonlinear convergence","SNESMonitorSet",&snesmonitor);
589: PetscOptionsName("-dmmg_jacobian_fd","Compute sparse Jacobian explicitly with finite differencing","DMMGSetSNES",&fdjacobian);
590: if (fdjacobian) jacobian = DMMGComputeJacobianWithFD;
591: #if defined(PETSC_HAVE_ADIC)
592: PetscOptionsName("-dmmg_jacobian_ad","Compute sparse Jacobian explicitly with ADIC (automatic differentiation)","DMMGSetSNES",&adjacobian);
593: if (adjacobian) jacobian = DMMGComputeJacobianWithAdic;
594: #endif
596: PetscOptionsTruthGroupBegin("-dmmg_jacobian_mf_fd_operator","Apply Jacobian via matrix free finite differencing","DMMGSetSNES",&mffdoperator);
597: PetscOptionsTruthGroupEnd("-dmmg_jacobian_mf_fd","Apply Jacobian via matrix free finite differencing even in computing preconditioner","DMMGSetSNES",&mffd);
598: if (mffd) mffdoperator = PETSC_TRUE;
599: #if defined(PETSC_HAVE_ADIC)
600: PetscOptionsTruthGroupBegin("-dmmg_jacobian_mf_ad_operator","Apply Jacobian via matrix free ADIC (automatic differentiation)","DMMGSetSNES",&mfadoperator);
601: PetscOptionsTruthGroupEnd("-dmmg_jacobian_mf_ad","Apply Jacobian via matrix free ADIC (automatic differentiation) even in computing preconditioner","DMMGSetSNES",&mfad);
602: if (mfad) mfadoperator = PETSC_TRUE;
603: #endif
604: MPI_Comm_size(dmmg[0]->comm,&size);
605: if (size == 1){
606: isctype = 0;
607: } else {
608: isctype = 1;
609: }
610: PetscOptionsEList("-dmmg_iscoloring_type","Type of ISColoring","None",ISColoringTypes,2,ISColoringTypes[isctype],&isctype,PETSC_NULL);
611:
612: PetscOptionsEnd();
614: /* create solvers for each level */
615: for (i=0; i<nlevels; i++) {
616: SNESCreate(dmmg[i]->comm,&dmmg[i]->snes);
617: SNESGetKSP(dmmg[i]->snes,&dmmg[i]->ksp);
618: if (snesmonitor) {
619: PetscObjectGetComm((PetscObject)dmmg[i]->snes,&comm);
620: PetscViewerASCIIMonitorCreate(comm,"stdout",nlevels-i,&ascii);
621: SNESMonitorSet(dmmg[i]->snes,SNESMonitorDefault,ascii,(PetscErrorCode(*)(void*))PetscViewerASCIIMonitorDestroy);
622: }
624: if (mffdoperator) {
625: MatCreateSNESMF(dmmg[i]->snes,dmmg[i]->x,&dmmg[i]->J);
626: VecDuplicate(dmmg[i]->x,&dmmg[i]->work1);
627: VecDuplicate(dmmg[i]->x,&dmmg[i]->work2);
628: MatSNESMFSetFunction(dmmg[i]->J,dmmg[i]->work1,function,dmmg[i]);
629: if (mffd) {
630: dmmg[i]->B = dmmg[i]->J;
631: jacobian = DMMGComputeJacobianWithMF;
632: }
633: #if defined(PETSC_HAVE_ADIC)
634: } else if (mfadoperator) {
635: MatRegisterDAAD();
636: MatCreateDAAD((DA)dmmg[i]->dm,&dmmg[i]->J);
637: MatDAADSetCtx(dmmg[i]->J,dmmg[i]->user);
638: if (mfad) {
639: dmmg[i]->B = dmmg[i]->J;
640: jacobian = DMMGComputeJacobianWithMF;
641: }
642: #endif
643: }
644:
645: if (!dmmg[i]->B) {
646: DMGetMatrix(dmmg[i]->dm,dmmg[i]->mtype,&dmmg[i]->B);
647: }
648: if (!dmmg[i]->J) {
649: dmmg[i]->J = dmmg[i]->B;
650: }
652: DMMGSetUpLevel(dmmg,dmmg[i]->ksp,i+1);
653:
654: /*
655: if the number of levels is > 1 then we want the coarse solve in the grid sequencing to use LU
656: when possible
657: */
658: if (nlevels > 1 && i == 0) {
659: PC pc;
660: KSP cksp;
661: PetscTruth flg1,flg2,flg3;
663: KSPGetPC(dmmg[i]->ksp,&pc);
664: PCMGGetCoarseSolve(pc,&cksp);
665: KSPGetPC(cksp,&pc);
666: PetscTypeCompare((PetscObject)pc,PCILU,&flg1);
667: PetscTypeCompare((PetscObject)pc,PCSOR,&flg2);
668: PetscTypeCompare((PetscObject)pc,PETSC_NULL,&flg3);
669: if (flg1 || flg2 || flg3) {
670: PCSetType(pc,PCLU);
671: }
672: }
674: dmmg[i]->solve = DMMGSolveSNES;
675: dmmg[i]->computejacobian = jacobian;
676: dmmg[i]->computefunction = function;
677: }
679: if (jacobian == DMMGComputeJacobianWithFD) {
680: ISColoring iscoloring;
681: ctype = (ISColoringType)isctype;
683: for (i=0; i<nlevels; i++) {
684: DMGetColoring(dmmg[i]->dm,ctype,&iscoloring);
685: MatFDColoringCreate(dmmg[i]->B,iscoloring,&dmmg[i]->fdcoloring);
686: ISColoringDestroy(iscoloring);
687: if (function == DMMGFormFunction) function = DMMGFormFunctionFD;
688: MatFDColoringSetFunction(dmmg[i]->fdcoloring,(PetscErrorCode(*)(void))function,dmmg[i]);
689: MatFDColoringSetFromOptions(dmmg[i]->fdcoloring);
690: }
691: #if defined(PETSC_HAVE_ADIC)
692: } else if (jacobian == DMMGComputeJacobianWithAdic) {
693: for (i=0; i<nlevels; i++) {
694: ISColoring iscoloring;
695: DMGetColoring(dmmg[i]->dm,IS_COLORING_GHOSTED,&iscoloring);
696: MatSetColoring(dmmg[i]->B,iscoloring);
697: ISColoringDestroy(iscoloring);
698: }
699: #endif
700: }
701: for (i=0; i<nlevels; i++) {
702: SNESSetJacobian(dmmg[i]->snes,dmmg[i]->J,dmmg[i]->B,DMMGComputeJacobian_Multigrid,dmmg);
703: SNESSetFunction(dmmg[i]->snes,dmmg[i]->b,function,dmmg[i]);
704: SNESSetFromOptions(dmmg[i]->snes);
705: }
707: /* Create interpolation scaling */
708: for (i=1; i<nlevels; i++) {
709: DMGetInterpolationScale(dmmg[i-1]->dm,dmmg[i]->dm,dmmg[i]->R,&dmmg[i]->Rscale);
710: }
712: PetscOptionsGetInt(PETSC_NULL,"-dmmg_jacobian_period",&period,PETSC_NULL);
713: for (i=0; i<nlevels; i++) {
714: dmmg[i]->updatejacobian = PETSC_TRUE;
715: dmmg[i]->updatejacobianperiod = period;
716: }
718: #if defined(PETSC_HAVE_ADIC)
719: {
720: PetscTruth flg;
721: PetscOptionsHasName(PETSC_NULL,"-dmmg_fas",&flg);
722: if (flg) {
723: PetscTruth block = PETSC_FALSE;
724: PetscTruth ngmres = PETSC_FALSE;
725: PetscInt newton_its;
726: PetscOptionsHasName(0,"-dmmg_fas_view",&flg);
727: for (i=0; i<nlevels; i++) {
728: NLFCreate_DAAD(&dmmg[i]->nlf);
729: NLFDAADSetDA_DAAD(dmmg[i]->nlf,(DA)dmmg[i]->dm);
730: NLFDAADSetCtx_DAAD(dmmg[i]->nlf,dmmg[i]->user);
731: NLFDAADSetResidual_DAAD(dmmg[i]->nlf,dmmg[i]->r);
732: VecDuplicate(dmmg[i]->b,&dmmg[i]->w);
734: dmmg[i]->monitor = PETSC_FALSE;
735: PetscOptionsHasName(0,"-dmmg_fas_monitor",&dmmg[i]->monitor);
736: dmmg[i]->monitorall = PETSC_FALSE;
737: PetscOptionsHasName(0,"-dmmg_fas_monitor_all",&dmmg[i]->monitorall);
738: dmmg[i]->presmooth = 2;
739: PetscOptionsGetInt(0,"-dmmg_fas_presmooth",&dmmg[i]->presmooth,0);
740: dmmg[i]->postsmooth = 2;
741: PetscOptionsGetInt(0,"-dmmg_fas_postsmooth",&dmmg[i]->postsmooth,0);
742: dmmg[i]->coarsesmooth = 2;
743: PetscOptionsGetInt(0,"-dmmg_fas_coarsesmooth",&dmmg[i]->coarsesmooth,0);
745: dmmg[i]->rtol = 1.e-8;
746: PetscOptionsGetReal(0,"-dmmg_fas_rtol",&dmmg[i]->rtol,0);
747: dmmg[i]->abstol = 1.e-50;
748: PetscOptionsGetReal(0,"-dmmg_fas_atol",&dmmg[i]->abstol,0);
750: newton_its = 2;
751: PetscOptionsGetInt(0,"-dmmg_fas_newton_its",&newton_its,0);
752: NLFDAADSetNewtonIterations_DAAD(dmmg[i]->nlf,newton_its);
754: if (flg) {
755: if (i == 0) {
756: PetscPrintf(dmmg[i]->comm,"FAS Solver Parameters\n");
757: PetscPrintf(dmmg[i]->comm," rtol %G atol %G\n",dmmg[i]->rtol,dmmg[i]->abstol);
758: PetscPrintf(dmmg[i]->comm," coarsesmooths %D\n",dmmg[i]->coarsesmooth);
759: PetscPrintf(dmmg[i]->comm," Newton iterations %D\n",newton_its);
760: } else {
761: PetscPrintf(dmmg[i]->comm," level %D presmooths %D\n",i,dmmg[i]->presmooth);
762: PetscPrintf(dmmg[i]->comm," postsmooths %D\n",dmmg[i]->postsmooth);
763: PetscPrintf(dmmg[i]->comm," Newton iterations %D\n",newton_its);
764: }
765: }
766: PetscOptionsHasName(0,"-dmmg_fas_block",&block);
767: PetscOptionsHasName(0,"-dmmg_fas_ngmres",&ngmres);
768: if (block) {
769: dmmg[i]->solve = DMMGSolveFASb;
770: if (flg) {
771: PetscPrintf(dmmg[i]->comm," using point-block smoothing\n");
772: }
773: } else if(ngmres) {
774: dmmg[i]->solve = DMMGSolveFAS_NGMRES;
775: if (flg) {
776: PetscPrintf(dmmg[i]->comm," using non-linear gmres\n");
777: }
778: } else {
779: dmmg[i]->solve = DMMGSolveFAS4;
780: }
781: }
782: }
783: }
784: #endif
785:
786: return(0);
787: }
791: /*@
792: DMMGSetSNESLocalFD - Sets the local user function that is used to approximately compute the Jacobian
793: via finite differences.
795: Collective on DMMG
797: Input Parameter:
798: + dmmg - the context
799: - function - the function that defines the nonlinear system
801: Level: intermediate
803: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNES(), DMMGSetSNESLocal()
805: @*/
806: PetscErrorCode DMMGSetSNESLocalFD(DMMG *dmmg,DALocalFunction1 function)
807: {
808: PetscInt i,nlevels = dmmg[0]->nlevels;
811: for (i=0; i<nlevels; i++) {
812: dmmg[i]->lfj = (PetscErrorCode (*)(void))function;
813: }
814: return(0);
815: }
818: /*M
819: DMMGSetSNESLocal - Sets the local user function that defines the nonlinear set of equations
820: that will use the grid hierarchy and (optionally) its derivative.
822: Collective on DMMG
824: Synopsis:
825: PetscErrorCode DMMGSetSNESLocal(DMMG *dmmg,DALocalFunction1 function, DALocalFunction1 jacobian,
826: DALocalFunction1 ad_function, DALocalFunction1 admf_function);
828: Input Parameter:
829: + dmmg - the context
830: . function - the function that defines the nonlinear system
831: . jacobian - function defines the local part of the Jacobian
832: . ad_function - the name of the function with an ad_ prefix. This is ignored if ADIC is
833: not installed
834: - admf_function - the name of the function with an ad_ prefix. This is ignored if ADIC is
835: not installed
837: Options Database Keys:
838: + -dmmg_snes_monitor
839: . -dmmg_jacobian_fd
840: . -dmmg_jacobian_ad
841: . -dmmg_jacobian_mf_fd_operator
842: . -dmmg_jacobian_mf_fd
843: . -dmmg_jacobian_mf_ad_operator
844: . -dmmg_jacobian_mf_ad
845: - -dmmg_jacobian_period <p> - Indicates how often in the SNES solve the Jacobian is recomputed (on all levels)
846: as suggested by Florin Dobrian if p is -1 then Jacobian is computed only on first
847: SNES iteration (i.e. -1 is equivalent to infinity)
850: Level: intermediate
852: Notes:
853: If ADIC or ADIFOR have been installed, this routine can use ADIC or ADIFOR to compute
854: the derivative; however, that function cannot call other functions except those in
855: standard C math libraries.
857: If ADIC/ADIFOR have not been installed and the Jacobian is not provided, this routine
858: uses finite differencing to approximate the Jacobian.
860: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNES()
862: M*/
866: PetscErrorCode DMMGSetSNESLocal_Private(DMMG *dmmg,DALocalFunction1 function,DALocalFunction1 jacobian,DALocalFunction1 ad_function,DALocalFunction1 admf_function)
867: {
869: PetscInt i,nlevels = dmmg[0]->nlevels,cookie;
870: PetscErrorCode (*computejacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*) = 0;
874: if (jacobian) computejacobian = DMMGComputeJacobian;
875: #if defined(PETSC_HAVE_ADIC)
876: else if (ad_function) computejacobian = DMMGComputeJacobianWithAdic;
877: #endif
878: CHKMEMQ;
879: PetscObjectGetCookie((PetscObject) dmmg[0]->dm,&cookie);
880: if (cookie == DA_COOKIE) {
881: PetscTruth flag;
882: PetscOptionsHasName(PETSC_NULL, "-dmmg_form_function_ghost", &flag);
883: if (flag) {
884: DMMGSetSNES(dmmg,DMMGFormFunctionGhost,computejacobian);
885: } else {
886: DMMGSetSNES(dmmg,DMMGFormFunction,computejacobian);
887: }
888: for (i=0; i<nlevels; i++) {
889: DASetLocalFunction((DA)dmmg[i]->dm,function);
890: dmmg[i]->lfj = (PetscErrorCode (*)(void))function;
891: DASetLocalJacobian((DA)dmmg[i]->dm,jacobian);
892: DASetLocalAdicFunction((DA)dmmg[i]->dm,ad_function);
893: DASetLocalAdicMFFunction((DA)dmmg[i]->dm,admf_function);
894: }
895: CHKMEMQ;
896: } else {
897: #ifdef PETSC_HAVE_SIEVE
898: DMMGSetSNES(dmmg, DMMGFormFunctionMesh, DMMGComputeJacobianMesh);
899: for (i=0; i<nlevels; i++) {
900: MeshSetLocalFunction((Mesh) dmmg[i]->dm, (PetscErrorCode (*)(Mesh,SectionReal,SectionReal,void*)) function);
901: dmmg[i]->lfj = (PetscErrorCode (*)(void)) function;
902: MeshSetLocalJacobian((Mesh) dmmg[i]->dm, (PetscErrorCode (*)(Mesh,SectionReal,Mat,void*)) jacobian);
903: }
904: CHKMEMQ;
905: #endif
906: }
907: CHKMEMQ;
908: return(0);
909: }
913: PetscErrorCode DMMGFunctioni(PetscInt i,Vec u,PetscScalar* r,void* ctx)
914: {
915: DMMG dmmg = (DMMG)ctx;
916: Vec U = dmmg->lwork1;
918: VecScatter gtol;
921: /* copy u into interior part of U */
922: DAGetScatter((DA)dmmg->dm,0,>ol,0);
923: VecScatterBegin(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);
924: VecScatterEnd(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);
925: DAFormFunctioni1((DA)dmmg->dm,i,U,r,dmmg->user);
926: return(0);
927: }
931: PetscErrorCode DMMGFunctionib(PetscInt i,Vec u,PetscScalar* r,void* ctx)
932: {
933: DMMG dmmg = (DMMG)ctx;
934: Vec U = dmmg->lwork1;
936: VecScatter gtol;
939: /* copy u into interior part of U */
940: DAGetScatter((DA)dmmg->dm,0,>ol,0);
941: VecScatterBegin(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);
942: VecScatterEnd(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);
943: DAFormFunctionib1((DA)dmmg->dm,i,U,r,dmmg->user);
944: return(0);
945: }
949: PetscErrorCode DMMGFunctioniBase(Vec u,void* ctx)
950: {
951: DMMG dmmg = (DMMG)ctx;
952: Vec U = dmmg->lwork1;
956: DAGlobalToLocalBegin((DA)dmmg->dm,u,INSERT_VALUES,U);
957: DAGlobalToLocalEnd((DA)dmmg->dm,u,INSERT_VALUES,U);
958: return(0);
959: }
963: PetscErrorCode DMMGSetSNESLocali_Private(DMMG *dmmg,PetscErrorCode (*functioni)(DALocalInfo*,MatStencil*,void*,PetscScalar*,void*),PetscErrorCode (*adi)(DALocalInfo*,MatStencil*,void*,void*,void*),PetscErrorCode (*adimf)(DALocalInfo*,MatStencil*,void*,void*,void*))
964: {
966: PetscInt i,nlevels = dmmg[0]->nlevels;
969: for (i=0; i<nlevels; i++) {
970: DASetLocalFunctioni((DA)dmmg[i]->dm,functioni);
971: DASetLocalAdicFunctioni((DA)dmmg[i]->dm,adi);
972: DASetLocalAdicMFFunctioni((DA)dmmg[i]->dm,adimf);
973: MatSNESMFSetFunctioni(dmmg[i]->J,DMMGFunctioni);
974: MatSNESMFSetFunctioniBase(dmmg[i]->J,DMMGFunctioniBase);
975: if (!dmmg[i]->lwork1) {
976: DACreateLocalVector((DA)dmmg[i]->dm,&dmmg[i]->lwork1);
977: }
978: }
979: return(0);
980: }
984: PetscErrorCode DMMGSetSNESLocalib_Private(DMMG *dmmg,PetscErrorCode (*functioni)(DALocalInfo*,MatStencil*,void*,PetscScalar*,void*),PetscErrorCode (*adi)(DALocalInfo*,MatStencil*,void*,void*,void*),PetscErrorCode (*adimf)(DALocalInfo*,MatStencil*,void*,void*,void*))
985: {
987: PetscInt i,nlevels = dmmg[0]->nlevels;
990: for (i=0; i<nlevels; i++) {
991: DASetLocalFunctionib((DA)dmmg[i]->dm,functioni);
992: DASetLocalAdicFunctionib((DA)dmmg[i]->dm,adi);
993: DASetLocalAdicMFFunctionib((DA)dmmg[i]->dm,adimf);
994: if (!dmmg[i]->lwork1) {
995: DACreateLocalVector((DA)dmmg[i]->dm,&dmmg[i]->lwork1);
996: }
997: }
998: return(0);
999: }
1001: static PetscErrorCode (*localfunc)(void) = 0;
1005: /*
1006: Uses the DM object to call the user provided function with the correct calling
1007: sequence.
1008: */
1009: PetscErrorCode DMMGInitialGuess_Local(DMMG dmmg,Vec x)
1010: {
1014: (*dmmg->dm->ops->forminitialguess)(dmmg->dm,localfunc,x,0);
1015: return(0);
1016: }
1020: /*@C
1021: DMMGSetInitialGuessLocal - sets code to compute the initial guess for each level
1023: Collective on DMMG
1025: Input Parameter:
1026: + dmmg - the context
1027: - localguess - the function that computes the initial guess
1029: Level: intermediate
1031: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNES(), DMMGSetInitialGuess(), DMMGSetSNESLocal()
1033: @*/
1034: PetscErrorCode DMMGSetInitialGuessLocal(DMMG *dmmg,PetscErrorCode (*localguess)(void))
1035: {
1039: localfunc = localguess; /* stash into ugly static for now */
1041: DMMGSetInitialGuess(dmmg,DMMGInitialGuess_Local);
1042: return(0);
1043: }