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", &section);
195:   SectionRealDuplicate(section, &sectionX);
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", &sectionX);
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,&gtol,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,&gtol,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: }