Actual source code: damgsnes.c

  1: /*$Id: damgsnes.c,v 1.51 2001/08/06 21:18:06 bsmith Exp $*/
  2: 
 3:  #include petscda.h
 4:  #include petscmg.h


  7: /*
  8:    Evaluates the Jacobian on all of the grids. It is used by DMMG to provide the 
  9:    ComputeJacobian() function that SNESSetJacobian() requires.
 10: */
 11: #undef __FUNCT__
 13: int DMMGComputeJacobian_Multigrid(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
 14: {
 15:   DMMG       *dmmg = (DMMG*)ptr;
 16:   int        ierr,i,nlevels = dmmg[0]->nlevels;
 17:   SLES       sles,lsles;
 18:   PC         pc;
 19:   PetscTruth ismg;
 20:   Vec        W;

 23:   if (!dmmg) SETERRQ(1,"Passing null as user context which should contain DMMG");

 25:   /* compute Jacobian on finest grid */
 26:   (*DMMGGetFine(dmmg)->computejacobian)(snes,X,J,B,flag,DMMGGetFine(dmmg));
 27:   MatSNESMFSetBase(DMMGGetFine(dmmg)->J,X);

 29:   /* create coarser grid Jacobians for preconditioner if multigrid is the preconditioner */
 30:   SNESGetSLES(snes,&sles);
 31:   SLESGetPC(sles,&pc);
 32:   PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
 33:   if (ismg) {

 35:     MGGetSmoother(pc,nlevels-1,&lsles);
 36:     SLESSetOperators(lsles,DMMGGetFine(dmmg)->J,DMMGGetFine(dmmg)->B,*flag);

 38:     for (i=nlevels-1; i>0; i--) {

 40:       if (!dmmg[i-1]->w) {
 41:         VecDuplicate(dmmg[i-1]->x,&dmmg[i-1]->w);
 42:       }

 44:       W    = dmmg[i-1]->w;
 45:       /* restrict X to coarser grid */
 46:       MatRestrict(dmmg[i]->R,X,W);
 47:       X    = W;

 49:       /* scale to "natural" scaling for that grid */
 50:       VecPointwiseMult(dmmg[i]->Rscale,X,X);

 52:       /* tell the base vector for matrix free multiplies */
 53:       MatSNESMFSetBase(dmmg[i-1]->J,X);

 55:       /* compute Jacobian on coarse grid */
 56:       (*dmmg[i-1]->computejacobian)(snes,X,&dmmg[i-1]->J,&dmmg[i-1]->B,flag,dmmg[i-1]);

 58:       MGGetSmoother(pc,i-1,&lsles);
 59:       SLESSetOperators(lsles,dmmg[i-1]->J,dmmg[i-1]->B,*flag);
 60:     }
 61:   }
 62:   return(0);
 63: }

 65: /* ---------------------------------------------------------------------------*/

 67: #undef __FUNCT__
 69: /* 
 70:    DMMGFormFunction - This is a universal global FormFunction used by the DMMG code
 71:    when the user provides a local function.

 73:    Input Parameters:
 74: +  snes - the SNES context
 75: .  X - input vector
 76: -  ptr - optional user-defined context, as set by SNESSetFunction()

 78:    Output Parameter:
 79: .  F - function vector

 81:  */
 82: int DMMGFormFunction(SNES snes,Vec X,Vec F,void *ptr)
 83: {
 84:   DMMG             dmmg = (DMMG)ptr;
 85:   int              ierr;
 86:   Vec              localX;
 87:   DA               da = (DA)dmmg->dm;

 90:   DAGetLocalVector(da,&localX);
 91:   /*
 92:      Scatter ghost points to local vector, using the 2-step process
 93:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
 94:   */
 95:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
 96:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
 97:   DAFormFunction1(da,localX,F,dmmg->user);
 98:   DARestoreLocalVector(da,&localX);
 99:   return(0);
100: }

102: #undef __FUNCT__
104: /*@C 
105:    SNESDAFormFunction - This is a universal function evaluation routine that
106:    may be used with SNESSetFunction() as long as the user context has a DA
107:    as its first record and the user has called DASetLocalFunction().

109:    Collective on SNES

111:    Input Parameters:
112: +  snes - the SNES context
113: .  X - input vector
114: .  F - function vector
115: -  ptr - pointer to a structure that must have a DA as its first entry. For example this 
116:          could be a DMMG

118:    Level: intermediate

120: .seealso: DASetLocalFunction(), DASetLocalJacobian(), DASetLocalAdicFunction(), DASetLocalAdicMFFunction(),
121:           SNESSetFunction(), SNESSetJacobian()

123: @*/
124: int SNESDAFormFunction(SNES snes,Vec X,Vec F,void *ptr)
125: {
126:   int              ierr;
127:   Vec              localX;
128:   DA               da = *(DA*)ptr;

131:   DAGetLocalVector(da,&localX);
132:   /*
133:      Scatter ghost points to local vector, using the 2-step process
134:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
135:   */
136:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
137:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
138:   DAFormFunction1(da,localX,F,ptr);
139:   DARestoreLocalVector(da,&localX);
140:   return(0);
141: }

143: /* ---------------------------------------------------------------------------------------------------------------------------*/

145: #undef __FUNCT__
147: int DMMGComputeJacobianWithFD(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
148: {
149:   int  ierr;
150:   DMMG dmmg = (DMMG)ctx;
151: 
153:   SNESDefaultComputeJacobianColor(snes,x1,J,B,flag,dmmg->fdcoloring);
154:   return(0);
155: }

157: #undef __FUNCT__
159: int DMMGComputeJacobianWithMF(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
160: {
161:   int  ierr;
162: 
164:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
165:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
166:   return(0);
167: }

169: #undef __FUNCT__
171: /*
172:     DMMGComputeJacobianWithAdic - Evaluates the Jacobian via Adic when the user has provided
173:     a local function evaluation routine.
174: */
175: int DMMGComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
176: {
177:   DMMG             dmmg = (DMMG) ptr;
178:   int              ierr;
179:   Vec              localX;
180:   DA               da = (DA) dmmg->dm;

183:   DAGetLocalVector(da,&localX);
184:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
185:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
186:   DAComputeJacobian1WithAdic(da,localX,*B,dmmg->user);
187:   DARestoreLocalVector(da,&localX);
188:   /* Assemble true Jacobian; if it is different */
189:   if (*J != *B) {
190:     ierr  = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
191:     ierr  = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
192:   }
193:   ierr  = MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
194:   *flag = SAME_NONZERO_PATTERN;
195:   return(0);
196: }

198: #undef __FUNCT__
200: /*@
201:     SNESDAComputeJacobianWithAdic - This is a universal Jacobian evaluation routine
202:     that may be used with SNESSetJacobian() as long as the user context has a DA as
203:     its first record and DASetLocalAdicFunction() has been called.  

205:    Collective on SNES

207:    Input Parameters:
208: +  snes - the SNES context
209: .  X - input vector
210: .  J - Jacobian
211: .  B - Jacobian used in preconditioner (usally same as J)
212: .  flag - indicates if the matrix changed its structure
213: -  ptr - optional user-defined context, as set by SNESSetFunction()

215:    Level: intermediate

217: .seealso: DASetLocalFunction(), DASetLocalAdicFunction(), SNESSetFunction(), SNESSetJacobian()

219: @*/
220: int SNESDAComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
221: {
222:   DA   da = *(DA*) ptr;
223:   int  ierr;
224:   Vec  localX;

227:   DAGetLocalVector(da,&localX);
228:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
229:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
230:   DAComputeJacobian1WithAdic(da,localX,*B,ptr);
231:   DARestoreLocalVector(da,&localX);
232:   /* Assemble true Jacobian; if it is different */
233:   if (*J != *B) {
234:     ierr  = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
235:     ierr  = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
236:   }
237:   ierr  = MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
238:   *flag = SAME_NONZERO_PATTERN;
239:   return(0);
240: }

242: #undef __FUNCT__
244: /*
245:     SNESDAComputeJacobianWithAdifor - This is a universal Jacobian evaluation routine
246:     that may be used with SNESSetJacobian() from Fortran as long as the user context has 
247:     a DA as its first record and DASetLocalAdiforFunction() has been called.  

249:    Collective on SNES

251:    Input Parameters:
252: +  snes - the SNES context
253: .  X - input vector
254: .  J - Jacobian
255: .  B - Jacobian used in preconditioner (usally same as J)
256: .  flag - indicates if the matrix changed its structure
257: -  ptr - optional user-defined context, as set by SNESSetFunction()

259:    Level: intermediate

261: .seealso: DASetLocalFunction(), DASetLocalAdicFunction(), SNESSetFunction(), SNESSetJacobian()

263: */
264: int SNESDAComputeJacobianWithAdifor(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
265: {
266:   DA   da = *(DA*) ptr;
267:   int  ierr;
268:   Vec  localX;

271:   DAGetLocalVector(da,&localX);
272:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
273:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
274:   DAComputeJacobian1WithAdifor(da,localX,*B,ptr);
275:   DARestoreLocalVector(da,&localX);
276:   /* Assemble true Jacobian; if it is different */
277:   if (*J != *B) {
278:     ierr  = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
279:     ierr  = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
280:   }
281:   ierr  = MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
282:   *flag = SAME_NONZERO_PATTERN;
283:   return(0);
284: }

286: #undef __FUNCT__
288: /*
289:    SNESDAComputeJacobian - This is a universal Jacobian evaluation routine for a
290:    locally provided Jacobian.

292:    Collective on SNES

294:    Input Parameters:
295: +  snes - the SNES context
296: .  X - input vector
297: .  J - Jacobian
298: .  B - Jacobian used in preconditioner (usally same as J)
299: .  flag - indicates if the matrix changed its structure
300: -  ptr - optional user-defined context, as set by SNESSetFunction()

302:    Level: intermediate

304: .seealso: DASetLocalFunction(), DASetLocalJacobian(), SNESSetFunction(), SNESSetJacobian()

306: */
307: int SNESDAComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
308: {
309:   DA   da = *(DA*) ptr;
310:   int  ierr;
311:   Vec  localX;

314:   DAGetLocalVector(da,&localX);
315:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
316:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
317:   DAComputeJacobian1(da,localX,*B,ptr);
318:   DARestoreLocalVector(da,&localX);
319:   /* Assemble true Jacobian; if it is different */
320:   if (*J != *B) {
321:     ierr  = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
322:     ierr  = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
323:   }
324:   ierr  = MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
325:   *flag = SAME_NONZERO_PATTERN;
326:   return(0);
327: }

329: #undef __FUNCT__  
331: int DMMGSolveSNES(DMMG *dmmg,int level)
332: {
333:   int  ierr,nlevels = dmmg[0]->nlevels,its;

336:   dmmg[0]->nlevels = level+1;
337:   ierr             = SNESSolve(dmmg[level]->snes,dmmg[level]->x,&its);
338:   dmmg[0]->nlevels = nlevels;
339:   return(0);
340: }

342: /* ===========================================================================================================*/

344: #undef __FUNCT__  
346: /*@C
347:     DMMGSetSNES - Sets the nonlinear function that defines the nonlinear set of equations
348:     to be solved using the grid hierarchy.

350:     Collective on DMMG

352:     Input Parameter:
353: +   dmmg - the context
354: .   function - the function that defines the nonlinear system
355: -   jacobian - optional function to compute Jacobian

357:     Options Database Keys:
358: +    -dmmg_snes_monitor
359: .    -dmmg_jacobian_fd
360: .    -dmmg_jacobian_ad
361: .    -dmmg_jacobian_mf_fd_operator
362: .    -dmmg_jacobian_mf_fd
363: .    -dmmg_jacobian_mf_ad_operator
364: -    -dmmg_jacobian_mf_ad

366:     Level: advanced

368: .seealso DMMGCreate(), DMMGDestroy, DMMGSetSLES(), DMMGSetSNESLocal()

370: @*/
371: int DMMGSetSNES(DMMG *dmmg,int (*function)(SNES,Vec,Vec,void*),int (*jacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*))
372: {
373:   int         ierr,i,nlevels = dmmg[0]->nlevels;
374:   PetscTruth  snesmonitor,mffdoperator,mffd,fdjacobian;
375: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
376:   PetscTruth  mfadoperator,mfad,adjacobian;
377: #endif
378:   SLES        sles;
379:   PetscViewer ascii;
380:   MPI_Comm    comm;

383:   if (!dmmg)     SETERRQ(1,"Passing null as DMMG");
384:   if (!jacobian) jacobian = DMMGComputeJacobianWithFD;

386:   PetscOptionsBegin(dmmg[0]->comm,PETSC_NULL,"DMMG Options","SNES");
387:     PetscOptionsName("-dmmg_snes_monitor","Monitor nonlinear convergence","SNESSetMonitor",&snesmonitor);


390:     PetscOptionsName("-dmmg_jacobian_fd","Compute sparse Jacobian explicitly with finite differencing","DMMGSetSNES",&fdjacobian);
391:     if (fdjacobian) jacobian = DMMGComputeJacobianWithFD;
392: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
393:     PetscOptionsName("-dmmg_jacobian_ad","Compute sparse Jacobian explicitly with ADIC (automatic differentiation)","DMMGSetSNES",&adjacobian);
394:     if (adjacobian) jacobian = DMMGComputeJacobianWithAdic;
395: #endif

397:     PetscOptionsLogicalGroupBegin("-dmmg_jacobian_mf_fd_operator","Apply Jacobian via matrix free finite differencing","DMMGSetSNES",&mffdoperator);
398:     PetscOptionsLogicalGroupEnd("-dmmg_jacobian_mf_fd","Apply Jacobian via matrix free finite differencing even in computing preconditioner","DMMGSetSNES",&mffd);
399:     if (mffd) mffdoperator = PETSC_TRUE;
400: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
401:     PetscOptionsLogicalGroupBegin("-dmmg_jacobian_mf_ad_operator","Apply Jacobian via matrix free ADIC (automatic differentiation)","DMMGSetSNES",&mfadoperator);
402:     PetscOptionsLogicalGroupEnd("-dmmg_jacobian_mf_ad","Apply Jacobian via matrix free ADIC (automatic differentiation) even in computing preconditioner","DMMGSetSNES",&mfad);
403:     if (mfad) mfadoperator = PETSC_TRUE;
404: #endif
405:   PetscOptionsEnd();

407:   /* create solvers for each level */
408:   for (i=0; i<nlevels; i++) {
409:     SNESCreate(dmmg[i]->comm,&dmmg[i]->snes);
410:     if (snesmonitor) {
411:       PetscObjectGetComm((PetscObject)dmmg[i]->snes,&comm);
412:       PetscViewerASCIIOpen(comm,"stdout",&ascii);
413:       PetscViewerASCIISetTab(ascii,nlevels-i);
414:       SNESSetMonitor(dmmg[i]->snes,SNESDefaultMonitor,ascii,(int(*)(void*))PetscViewerDestroy);
415:     }

417:     if (mffdoperator) {
418:       MatCreateSNESMF(dmmg[i]->snes,dmmg[i]->x,&dmmg[i]->J);
419:       VecDuplicate(dmmg[i]->x,&dmmg[i]->work1);
420:       VecDuplicate(dmmg[i]->x,&dmmg[i]->work2);
421:       MatSNESMFSetFunction(dmmg[i]->J,dmmg[i]->work1,function,dmmg[i]);
422:       if (mffd) {
423:         dmmg[i]->B = dmmg[i]->J;
424:         jacobian   = DMMGComputeJacobianWithMF;
425:       }
426: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
427:     } else if (mfadoperator) {
428:       MatRegisterDAAD();
429:       MatCreateDAAD((DA)dmmg[i]->dm,&dmmg[i]->J);
430:       MatDAADSetCtx(dmmg[i]->J,dmmg[i]->user);
431:       if (mfad) {
432:         dmmg[i]->B = dmmg[i]->J;
433:         jacobian   = DMMGComputeJacobianWithMF;
434:       }
435: #endif
436:     }
437: 
438:     if (!dmmg[i]->B) {
439:       DMGetMatrix(dmmg[i]->dm,MATMPIAIJ,&dmmg[i]->B);
440:     }
441:     if (!dmmg[i]->J) {
442:       dmmg[i]->J = dmmg[i]->B;
443:     }

445:     SNESGetSLES(dmmg[i]->snes,&sles);
446:     DMMGSetUpLevel(dmmg,sles,i+1);
447: 
448:     /*
449:        if the number of levels is > 1 then we want the coarse solve in the grid sequencing to use LU
450:        when possible 
451:     */
452:     if (nlevels > 1 && i == 0) {
453:       PC         pc;
454:       SLES       csles;
455:       PetscTruth flg1,flg2,flg3;

457:       SLESGetPC(sles,&pc);
458:       MGGetCoarseSolve(pc,&csles);
459:       SLESGetPC(csles,&pc);
460:       PetscTypeCompare((PetscObject)pc,PCILU,&flg1);
461:       PetscTypeCompare((PetscObject)pc,PCSOR,&flg2);
462:       PetscTypeCompare((PetscObject)pc,PETSC_NULL,&flg3);
463:       if (flg1 || flg2 || flg3) {
464:         PCSetType(pc,PCLU);
465:       }
466:     }

468:     SNESSetFromOptions(dmmg[i]->snes);
469:     dmmg[i]->solve           = DMMGSolveSNES;
470:     dmmg[i]->computejacobian = jacobian;
471:     dmmg[i]->computefunction = function;
472:   }


475:   if (jacobian == DMMGComputeJacobianWithFD) {
476:     ISColoring iscoloring;
477:     for (i=0; i<nlevels; i++) {
478:       DMGetColoring(dmmg[i]->dm,IS_COLORING_LOCAL,&iscoloring);
479:       MatFDColoringCreate(dmmg[i]->B,iscoloring,&dmmg[i]->fdcoloring);
480:       ISColoringDestroy(iscoloring);
481:       MatFDColoringSetFunction(dmmg[i]->fdcoloring,(int(*)(void))function,dmmg[i]);
482:       MatFDColoringSetFromOptions(dmmg[i]->fdcoloring);
483:     }
484: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
485:   } else if (jacobian == DMMGComputeJacobianWithAdic) {
486:     for (i=0; i<nlevels; i++) {
487:       ISColoring iscoloring;
488:       DMGetColoring(dmmg[i]->dm,IS_COLORING_GHOSTED,&iscoloring);
489:       MatSetColoring(dmmg[i]->B,iscoloring);
490:       ISColoringDestroy(iscoloring);
491:     }
492: #endif
493:   }

495:   for (i=0; i<nlevels; i++) {
496:     SNESSetJacobian(dmmg[i]->snes,dmmg[i]->J,dmmg[i]->B,DMMGComputeJacobian_Multigrid,dmmg);
497:     SNESSetFunction(dmmg[i]->snes,dmmg[i]->b,function,dmmg[i]);
498:   }

500:   /* Create interpolation scaling */
501:   for (i=1; i<nlevels; i++) {
502:     DMGetInterpolationScale(dmmg[i-1]->dm,dmmg[i]->dm,dmmg[i]->R,&dmmg[i]->Rscale);
503:   }

505:   return(0);
506: }

508: #undef __FUNCT__  
510: /*@C
511:     DMMGSetInitialGuess - Sets the function that computes an initial guess, if not given
512:     uses 0.

514:     Collective on DMMG and SNES

516:     Input Parameter:
517: +   dmmg - the context
518: -   guess - the function

520:     Level: advanced

522: .seealso DMMGCreate(), DMMGDestroy, DMMGSetSLES()

524: @*/
525: int DMMGSetInitialGuess(DMMG *dmmg,int (*guess)(SNES,Vec,void*))
526: {
527:   int i,nlevels = dmmg[0]->nlevels;

530:   for (i=0; i<nlevels; i++) {
531:     dmmg[i]->initialguess = guess;
532:   }
533:   return(0);
534: }


537: /*M
538:     DMMGSetSNESLocal - Sets the local user function that defines the nonlinear set of equations
539:     that will use the grid hierarchy and (optionally) its derivative.

541:     Collective on DMMG

543:    Synopsis:
544:    int DMMGSetSNESLocal(DMMG *dmmg,DALocalFunction1 function, DALocalFunction1 jacobian,
545:                         DALocalFunction1 ad_function, DALocalFunction1 admf_function);

547:     Input Parameter:
548: +   dmmg - the context
549: .   function - the function that defines the nonlinear system
550: .   jacobian - function defines the local part of the Jacobian (not currently supported)
551: .   ad_function - the name of the function with an ad_ prefix. This is ignored if ADIC is
552:                   not installed
553: -   admf_function - the name of the function with an ad_ prefix. This is ignored if ADIC is
554:                   not installed

556:     Options Database Keys:
557: +    -dmmg_snes_monitor
558: .    -dmmg_jacobian_fd
559: .    -dmmg_jacobian_ad
560: .    -dmmg_jacobian_mf_fd_operator
561: .    -dmmg_jacobian_mf_fd
562: .    -dmmg_jacobian_mf_ad_operator
563: -    -dmmg_jacobian_mf_ad


566:     Level: intermediate

568:     Notes: 
569:     If ADIC or ADIFOR have been installed, this routine can use ADIC or ADIFOR to compute
570:     the derivative; however, that function cannot call other functions except those in
571:     standard C math libraries.

573:     If ADIC/ADIFOR have not been installed and the Jacobian is not provided, this routine
574:     uses finite differencing to approximate the Jacobian.

576: .seealso DMMGCreate(), DMMGDestroy, DMMGSetSLES(), DMMGSetSNES()

578: M*/

580: #undef __FUNCT__  
582: int DMMGSetSNESLocal_Private(DMMG *dmmg,DALocalFunction1 function,DALocalFunction1 jacobian,DALocalFunction1 ad_function,DALocalFunction1 admf_function)
583: {
584:   int ierr,i,nlevels = dmmg[0]->nlevels;
585:   int (*computejacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*) = 0;


589:   if (jacobian)         computejacobian = SNESDAComputeJacobian;
590: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
591:   else if (ad_function) computejacobian = DMMGComputeJacobianWithAdic;
592: #endif

594:   DMMGSetSNES(dmmg,DMMGFormFunction,computejacobian);
595:   for (i=0; i<nlevels; i++) {
596:     DASetLocalFunction((DA)dmmg[i]->dm,function);
597:     DASetLocalJacobian((DA)dmmg[i]->dm,jacobian);
598:     DASetLocalAdicFunction((DA)dmmg[i]->dm,ad_function);
599:     DASetLocalAdicMFFunction((DA)dmmg[i]->dm,admf_function);
600:   }
601:   return(0);
602: }

604: #undef __FUNCT__  
606: static int DMMGFunctioni(int i,Vec u,PetscScalar* r,void* ctx)
607: {
608:   DMMG       dmmg = (DMMG)ctx;
609:   Vec        U = dmmg->lwork1;
610:   int        ierr;
611:   VecScatter gtol;

614:   /* copy u into interior part of U */
615:   DAGetScatter((DA)dmmg->dm,0,&gtol,0);
616:   VecScatterBegin(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);
617:   VecScatterEnd(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);

619:   DAFormFunctioni1((DA)dmmg->dm,i,U,r,dmmg->user);
620:   return(0);
621: }

623: #undef __FUNCT__  
625: static int DMMGFunctioniBase(Vec u,void* ctx)
626: {
627:   DMMG dmmg = (DMMG)ctx;
628:   Vec  U = dmmg->lwork1;
629:   int  ierr;

632:   DAGlobalToLocalBegin((DA)dmmg->dm,u,INSERT_VALUES,U);
633:   DAGlobalToLocalEnd((DA)dmmg->dm,u,INSERT_VALUES,U);
634:   return(0);
635: }

637: #undef __FUNCT__  
639: int DMMGSetSNESLocali_Private(DMMG *dmmg,int (*functioni)(DALocalInfo*,MatStencil*,void*,PetscScalar*,void*),int (*adi)(DALocalInfo*,MatStencil*,void*,void*,void*),int (*adimf)(DALocalInfo*,MatStencil*,void*,void*,void*))
640: {
641:   int ierr,i,nlevels = dmmg[0]->nlevels;

644: CHKMEMQ;
645:   for (i=0; i<nlevels; i++) {
646: CHKMEMQ;
647:     DASetLocalFunctioni((DA)dmmg[i]->dm,functioni);
648: CHKMEMQ;
649:     DASetLocalAdicFunctioni((DA)dmmg[i]->dm,adi);
650: CHKMEMQ;
651:     DASetLocalAdicMFFunctioni((DA)dmmg[i]->dm,adimf);
652: CHKMEMQ;
653:     MatSNESMFSetFunctioni(dmmg[i]->J,DMMGFunctioni);
654: CHKMEMQ;
655:     MatSNESMFSetFunctioniBase(dmmg[i]->J,DMMGFunctioniBase);
656: CHKMEMQ;
657:     DACreateLocalVector((DA)dmmg[i]->dm,&dmmg[i]->lwork1);
658:   }
659:   return(0);
660: }


663: #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
664: EXTERN_C_BEGIN
665: #include "adic/ad_utils.h"
666: EXTERN_C_END

668: #undef __FUNCT__  
670: int PetscADView(int N,int nc,double *ptr,PetscViewer viewer)
671: {
672:   int        i,j,nlen  = PetscADGetDerivTypeSize();
673:   char       *cptr = (char*)ptr;
674:   double     *values;

677:   for (i=0; i<N; i++) {
678:     printf("Element %d value %g derivatives: ",i,*(double*)cptr);
679:     values = PetscADGetGradArray(cptr);
680:     for (j=0; j<nc; j++) {
681:       printf("%g ",*values++);
682:     }
683:     printf("n");
684:     cptr += nlen;
685:   }

687:   return(0);
688: }

690: #endif