Actual source code: damgsnes.c

      1: #define PETSCSNES_DLL
      2:
     3:  #include petscda.h
     4:  #include petscmg.h
     5:  #include petscdmmg.h

      7: /*
      8:       period of -1 indicates update only on zeroth iteration of SNES
      9: */
     10: #define ShouldUpdate(l,it) (((dmmg[l-1]->updatejacobianperiod == -1) && (it == 0)) || \
     11:                             ((dmmg[l-1]->updatejacobianperiod >   0) && !(it % dmmg[l-1]->updatejacobianperiod)))
     12: /*
     13:    Evaluates the Jacobian on all of the grids. It is used by DMMG to provide the
     14:    ComputeJacobian() function that SNESSetJacobian() requires.
     15: */
     18: PetscErrorCode DMMGComputeJacobian_Multigrid(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
     19: {
     20:   DMMG           *dmmg = (DMMG*)ptr;
     22:   PetscInt       i,nlevels = dmmg[0]->nlevels,it;
     23:   KSP            ksp,lksp;
     24:   PC             pc;
     25:   PetscTruth     ismg;
     26:   Vec            W;
     27:   MatStructure   flg;

     30:   if (!dmmg) SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as user context which should contain DMMG");
     31:   SNESGetIterationNumber(snes,&it);

     33:   /* compute Jacobian on finest grid */
     34:   if (dmmg[nlevels-1]->updatejacobian && ShouldUpdate(nlevels,it)) {
     35:     (*DMMGGetFine(dmmg)->computejacobian)(snes,X,J,B,flag,DMMGGetFine(dmmg));
     36:   } else {
     37:     PetscLogInfo((0,"DMMGComputeJacobian_Multigrid:Skipping Jacobian, SNES iteration %D frequence %D level %D\n",it,dmmg[nlevels-1]->updatejacobianperiod,nlevels-1));
     38:     *flag = SAME_PRECONDITIONER;
     39:   }
     40:   MatSNESMFSetBase(DMMGGetFine(dmmg)->J,X);

     42:   /* create coarser grid Jacobians for preconditioner if multigrid is the preconditioner */
     43:   SNESGetKSP(snes,&ksp);
     44:   KSPGetPC(ksp,&pc);
     45:   PetscTypeCompare((PetscObject)pc,PCMG,&ismg);
     46:   if (ismg) {

     48:     PCMGGetSmoother(pc,nlevels-1,&lksp);
     49:     KSPSetOperators(lksp,DMMGGetFine(dmmg)->J,DMMGGetFine(dmmg)->B,*flag);

     51:     if (dmmg[0]->galerkin) {
     52:       for (i=nlevels-2; i>-1; i--) {
     53:         PetscTruth JeqB = (PetscTruth)( dmmg[i]->B == dmmg[i]->J);
     54:         MatDestroy(dmmg[i]->B);
     55:         MatPtAP(dmmg[i+1]->B,dmmg[i+1]->R,MAT_INITIAL_MATRIX,1.0,&dmmg[i]->B);
     56:         if (JeqB) dmmg[i]->J = dmmg[i]->B;
     57:         PCMGGetSmoother(pc,i,&lksp);
     58:         KSPSetOperators(lksp,dmmg[i]->J,dmmg[i]->B,*flag);
     59:       }
     60:     } else {
     61:       for (i=nlevels-1; i>0; i--) {
     62:         if (!dmmg[i-1]->w) {
     63:           VecDuplicate(dmmg[i-1]->x,&dmmg[i-1]->w);
     64:         }
     65:         W    = dmmg[i-1]->w;
     66:         /* restrict X to coarser grid */
     67:         MatRestrict(dmmg[i]->R,X,W);
     68:         X    = W;
     69:         /* scale to "natural" scaling for that grid */
     70:         VecPointwiseMult(X,X,dmmg[i]->Rscale);
     71:         /* tell the base vector for matrix free multiplies */
     72:         MatSNESMFSetBase(dmmg[i-1]->J,X);
     73:         /* compute Jacobian on coarse grid */
     74:         if (dmmg[i-1]->updatejacobian && ShouldUpdate(i,it)) {
     75:           (*dmmg[i-1]->computejacobian)(snes,X,&dmmg[i-1]->J,&dmmg[i-1]->B,&flg,dmmg[i-1]);
     76:           flg = SAME_NONZERO_PATTERN;
     77:         } else {
     78:           PetscLogInfo((0,"DMMGComputeJacobian_Multigrid:Skipping Jacobian, SNES iteration %D frequence %D level %D\n",it,dmmg[i-1]->updatejacobianperiod,i-1));
     79:           flg = SAME_PRECONDITIONER;
     80:         }
     81:         PCMGGetSmoother(pc,i-1,&lksp);
     82:         KSPSetOperators(lksp,dmmg[i-1]->J,dmmg[i-1]->B,flg);
     83:       }
     84:     }
     85:   }
     86:   return(0);
     87: }

     89: /* ---------------------------------------------------------------------------*/

     93: /*
     94:    DMMGFormFunction - This is a universal global FormFunction used by the DMMG code
     95:    when the user provides a local function.

     97:    Input Parameters:
     98: +  snes - the SNES context
     99: .  X - input vector
    100: -  ptr - optional user-defined context, as set by SNESSetFunction()

    102:    Output Parameter:
    103: .  F - function vector

    105:  */
    106: PetscErrorCode DMMGFormFunction(SNES snes,Vec X,Vec F,void *ptr)
    107: {
    108:   DMMG           dmmg = (DMMG)ptr;
    110:   Vec            localX;
    111:   DA             da = (DA)dmmg->dm;

    114:   DAGetLocalVector(da,&localX);
    115:   /*
    116:      Scatter ghost points to local vector, using the 2-step process
    117:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
    118:   */
    119:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    120:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    121:   DAFormFunction1(da,localX,F,dmmg->user);
    122:   DARestoreLocalVector(da,&localX);
    123:   return(0);
    124: }

    128: /*@C
    129:    SNESDAFormFunction - This is a universal function evaluation routine that
    130:    may be used with SNESSetFunction() as long as the user context has a DA
    131:    as its first record and the user has called DASetLocalFunction().

    133:    Collective on SNES

    135:    Input Parameters:
    136: +  snes - the SNES context
    137: .  X - input vector
    138: .  F - function vector
    139: -  ptr - pointer to a structure that must have a DA as its first entry. For example this
    140:          could be a DMMG

    142:    Level: intermediate

    144: .seealso: DASetLocalFunction(), DASetLocalJacobian(), DASetLocalAdicFunction(), DASetLocalAdicMFFunction(),
    145:           SNESSetFunction(), SNESSetJacobian()

    147: @*/
    148: PetscErrorCode PETSCSNES_DLLEXPORT SNESDAFormFunction(SNES snes,Vec X,Vec F,void *ptr)
    149: {
    151:   Vec            localX;
    152:   DA             da = *(DA*)ptr;

    155:   DAGetLocalVector(da,&localX);
    156:   /*
    157:      Scatter ghost points to local vector, using the 2-step process
    158:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
    159:   */
    160:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    161:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    162:   DAFormFunction1(da,localX,F,ptr);
    163:   if (PetscExceptionValue(ierr)) {
    164:     PetscErrorCode pDARestoreLocalVector(da,&localX);CHKERRQ(pierr);
    165:   }
    166:
    167:   DARestoreLocalVector(da,&localX);
    168:   return(0);
    169: }

    171: /* ---------------------------------------------------------------------------------------------------------------------------*/

    175: PetscErrorCode DMMGComputeJacobianWithFD(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
    176: {
    178:   DMMG           dmmg = (DMMG)ctx;
    179:
    181:   SNESDefaultComputeJacobianColor(snes,x1,J,B,flag,dmmg->fdcoloring);
    182:   return(0);
    183: }

    187: PetscErrorCode DMMGComputeJacobianWithMF(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
    188: {
    190:
    192:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    193:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    194:   return(0);
    195: }

    197: #if defined(PETSC_HAVE_ADIC)
    200: /*
    201:     DMMGComputeJacobianWithAdic - Evaluates the Jacobian via Adic when the user has provided
    202:     a local function evaluation routine.
    203: */
    204: PetscErrorCode DMMGComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    205: {
    206:   DMMG           dmmg = (DMMG) ptr;
    208:   Vec            localX;
    209:   DA             da = (DA) dmmg->dm;

    212:   DAGetLocalVector(da,&localX);
    213:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    214:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    215:   DAComputeJacobian1WithAdic(da,localX,*B,dmmg->user);
    216:   DARestoreLocalVector(da,&localX);
    217:   /* Assemble true Jacobian; if it is different */
    218:   if (*J != *B) {
    219:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    220:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    221:   }
    222:   MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
    223:   *flag = SAME_NONZERO_PATTERN;
    224:   return(0);
    225: }
    226: #endif

    230: /*
    231:     DMMGComputeJacobian - Evaluates the Jacobian when the user has provided
    232:     a local function evaluation routine.
    233: */
    234: PetscErrorCode DMMGComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    235: {
    236:   DMMG           dmmg = (DMMG) ptr;
    238:   Vec            localX;
    239:   DA             da = (DA) dmmg->dm;

    242:   DAGetLocalVector(da,&localX);
    243:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    244:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    245:   DAComputeJacobian1(da,localX,*B,dmmg->user);
    246:   DARestoreLocalVector(da,&localX);
    247:   /* Assemble true Jacobian; if it is different */
    248:   if (*J != *B) {
    249:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    250:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    251:   }
    252:   MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
    253:   *flag = SAME_NONZERO_PATTERN;
    254:   return(0);
    255: }

    257: #if defined(PETSC_HAVE_ADIC)
    260: /*@
    261:     SNESDAComputeJacobianWithAdic - This is a universal Jacobian evaluation routine
    262:     that may be used with SNESSetJacobian() as long as the user context has a DA as
    263:     its first record and DASetLocalAdicFunction() has been called.

    265:    Collective on SNES

    267:    Input Parameters:
    268: +  snes - the SNES context
    269: .  X - input vector
    270: .  J - Jacobian
    271: .  B - Jacobian used in preconditioner (usally same as J)
    272: .  flag - indicates if the matrix changed its structure
    273: -  ptr - optional user-defined context, as set by SNESSetFunction()

    275:    Level: intermediate

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

    279: @*/
    280: PetscErrorCode PETSCSNES_DLLEXPORT SNESDAComputeJacobianWithAdic(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    281: {
    282:   DA             da = *(DA*) ptr;
    284:   Vec            localX;

    287:   DAGetLocalVector(da,&localX);
    288:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    289:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    290:   DAComputeJacobian1WithAdic(da,localX,*B,ptr);
    291:   DARestoreLocalVector(da,&localX);
    292:   /* Assemble true Jacobian; if it is different */
    293:   if (*J != *B) {
    294:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    295:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    296:   }
    297:   MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
    298:   *flag = SAME_NONZERO_PATTERN;
    299:   return(0);
    300: }
    301: #endif

    305: /*
    306:     SNESDAComputeJacobianWithAdifor - This is a universal Jacobian evaluation routine
    307:     that may be used with SNESSetJacobian() from Fortran as long as the user context has
    308:     a DA as its first record and DASetLocalAdiforFunction() has been called.

    310:    Collective on SNES

    312:    Input Parameters:
    313: +  snes - the SNES context
    314: .  X - input vector
    315: .  J - Jacobian
    316: .  B - Jacobian used in preconditioner (usally same as J)
    317: .  flag - indicates if the matrix changed its structure
    318: -  ptr - optional user-defined context, as set by SNESSetFunction()

    320:    Level: intermediate

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

    324: */
    325: PetscErrorCode PETSCSNES_DLLEXPORT SNESDAComputeJacobianWithAdifor(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    326: {
    327:   DA             da = *(DA*) ptr;
    329:   Vec            localX;

    332:   DAGetLocalVector(da,&localX);
    333:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    334:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    335:   DAComputeJacobian1WithAdifor(da,localX,*B,ptr);
    336:   DARestoreLocalVector(da,&localX);
    337:   /* Assemble true Jacobian; if it is different */
    338:   if (*J != *B) {
    339:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    340:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    341:   }
    342:   MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
    343:   *flag = SAME_NONZERO_PATTERN;
    344:   return(0);
    345: }

    349: /*
    350:    SNESDAComputeJacobian - This is a universal Jacobian evaluation routine for a
    351:    locally provided Jacobian.

    353:    Collective on SNES

    355:    Input Parameters:
    356: +  snes - the SNES context
    357: .  X - input vector
    358: .  J - Jacobian
    359: .  B - Jacobian used in preconditioner (usally same as J)
    360: .  flag - indicates if the matrix changed its structure
    361: -  ptr - optional user-defined context, as set by SNESSetFunction()

    363:    Level: intermediate

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

    367: */
    368: PetscErrorCode PETSCSNES_DLLEXPORT SNESDAComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    369: {
    370:   DA             da = *(DA*) ptr;
    372:   Vec            localX;

    375:   DAGetLocalVector(da,&localX);
    376:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
    377:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
    378:   DAComputeJacobian1(da,localX,*B,ptr);
    379:   DARestoreLocalVector(da,&localX);
    380:   /* Assemble true Jacobian; if it is different */
    381:   if (*J != *B) {
    382:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    383:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
    384:   }
    385:   MatSetOption(*B,MAT_NEW_NONZERO_LOCATION_ERR);
    386:   *flag = SAME_NONZERO_PATTERN;
    387:   return(0);
    388: }

    392: PetscErrorCode DMMGSolveSNES(DMMG *dmmg,PetscInt level)
    393: {
    395:   PetscInt       nlevels = dmmg[0]->nlevels;

    398:   dmmg[0]->nlevels = level+1;
    399:   SNESSolve(dmmg[level]->snes,PETSC_NULL,dmmg[level]->x);
    400:   dmmg[0]->nlevels = nlevels;
    401:   return(0);
    402: }

    405: EXTERN PetscErrorCode PETSCSNES_DLLEXPORT NLFCreate_DAAD(NLF*);
    406: EXTERN PetscErrorCode PETSCSNES_DLLEXPORT NLFRelax_DAAD(NLF,MatSORType,PetscInt,Vec);
    407: EXTERN PetscErrorCode PETSCSNES_DLLEXPORT NLFDAADSetDA_DAAD(NLF,DA);
    408: EXTERN PetscErrorCode PETSCSNES_DLLEXPORT NLFDAADSetCtx_DAAD(NLF,void*);
    409: EXTERN PetscErrorCode PETSCSNES_DLLEXPORT NLFDAADSetResidual_DAAD(NLF,Vec);
    410: EXTERN PetscErrorCode PETSCSNES_DLLEXPORT NLFDAADSetNewtonIterations_DAAD(NLF,PetscInt);

    413: #if defined(PETSC_HAVE_ADIC)
    414:  #include src/ksp/pc/impls/mg/mgimpl.h
    415: /*
    416:           This is pre-beta FAS code. It's design should not be taken seriously!
    417: */
    420: PetscErrorCode DMMGSolveFAS(DMMG *dmmg,PetscInt level)
    421: {
    423:   PetscInt       i,j,k;
    424:   PetscReal      norm;
    425:   PetscScalar    zero = 0.0,mone = -1.0,one = 1.0;
    426:   MG             *mg;
    427:   PC             pc;

    430:   VecSet(dmmg[level]->r,zero);
    431:   for (j=1; j<=level; j++) {
    432:     if (!dmmg[j]->inject) {
    433:       DMGetInjection(dmmg[j-1]->dm,dmmg[j]->dm,&dmmg[j]->inject);
    434:     }
    435:   }

    437:   KSPGetPC(dmmg[level]->ksp,&pc);
    438:   mg   = ((MG*)pc->data);

    440:   for (i=0; i<100; i++) {

    442:     for (j=level; j>0; j--) {

    444:       /* Relax residual_fine - F(x_fine) = 0 */
    445:       for (k=0; k<dmmg[j]->presmooth; k++) {
    446:         NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
    447:       }

    449:       /* R*(residual_fine - F(x_fine)) */
    450:       DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
    451:       VecAYPX(dmmg[j]->w,mone,dmmg[j]->r);

    453:       if (j == level || dmmg[j]->monitorall) {
    454:         /* norm( residual_fine - f(x_fine) ) */
    455:         VecNorm(dmmg[j]->w,NORM_2,&norm);
    456:         if (j == level) {
    457:           if (norm < dmmg[level]->abstol) goto theend;
    458:           if (i == 0) {
    459:             dmmg[level]->rrtol = norm*dmmg[level]->rtol;
    460:           } else {
    461:             if (norm < dmmg[level]->rrtol) goto theend;
    462:           }
    463:         }
    464:       }

    466:       if (dmmg[j]->monitorall) {
    467:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
    468:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
    469:       }
    470:       MatRestrict(mg[j]->restrct,dmmg[j]->w,dmmg[j-1]->r);
    471:
    472:       /* F(R*x_fine) */
    473:       VecScatterBegin(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
    474:       VecScatterEnd(dmmg[j]->x,dmmg[j-1]->x,INSERT_VALUES,SCATTER_FORWARD,dmmg[j]->inject);
    475:       DMMGFormFunction(0,dmmg[j-1]->x,dmmg[j-1]->w,dmmg[j-1]);

    477:       /* residual_coarse = F(R*x_fine) + R*(residual_fine - F(x_fine)) */
    478:       VecAYPX(dmmg[j-1]->r,one,dmmg[j-1]->w);

    480:       /* save R*x_fine into b (needed when interpolating compute x back up */
    481:       VecCopy(dmmg[j-1]->x,dmmg[j-1]->b);
    482:     }

    484:     for (j=0; j<dmmg[0]->presmooth; j++) {
    485:       NLFRelax_DAAD(dmmg[0]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[0]->x);
    486:     }
    487:     if (dmmg[0]->monitorall){
    488:       DMMGFormFunction(0,dmmg[0]->x,dmmg[0]->w,dmmg[0]);
    489:       VecAXPY(dmmg[0]->w,mone,dmmg[0]->r);
    490:       VecNorm(dmmg[0]->w,NORM_2,&norm);
    491:       for (k=0; k<level+1; k++) {PetscPrintf(dmmg[0]->comm,"  ");}
    492:       PetscPrintf(dmmg[0]->comm,"FAS coarse grid function norm %g\n",norm);
    493:     }

    495:     for (j=1; j<=level; j++) {
    496:       /* x_fine = x_fine + R'*(x_coarse - R*x_fine) */
    497:       VecAXPY(dmmg[j-1]->x,mone,dmmg[j-1]->b);
    498:       MatInterpolateAdd(mg[j]->interpolate,dmmg[j-1]->x,dmmg[j]->x,dmmg[j]->x);

    500:       if (dmmg[j]->monitorall) {
    501:         /* norm( F(x_fine) - residual_fine ) */
    502:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
    503:         VecAXPY(dmmg[j]->w,mone,dmmg[j]->r);
    504:         VecNorm(dmmg[j]->w,NORM_2,&norm);
    505:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
    506:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
    507:       }

    509:       /* Relax residual_fine - F(x_fine)  = 0 */
    510:       for (k=0; k<dmmg[j]->postsmooth; k++) {
    511:         NLFRelax_DAAD(dmmg[j]->nlf,SOR_SYMMETRIC_SWEEP,1,dmmg[j]->x);
    512:       }

    514:       if (dmmg[j]->monitorall) {
    515:         /* norm( F(x_fine) - residual_fine ) */
    516:         DMMGFormFunction(0,dmmg[j]->x,dmmg[j]->w,dmmg[j]);
    517:         VecAXPY(dmmg[j]->w,mone,dmmg[j]->r);
    518:         VecNorm(dmmg[j]->w,NORM_2,&norm);
    519:         for (k=0; k<level-j+1; k++) {PetscPrintf(dmmg[j]->comm,"  ");}
    520:         PetscPrintf(dmmg[j]->comm,"FAS function norm %g\n",norm);
    521:       }
    522:     }

    524:     if (dmmg[level]->monitor){
    525:       DMMGFormFunction(0,dmmg[level]->x,dmmg[level]->w,dmmg[level]);
    526:       VecNorm(dmmg[level]->w,NORM_2,&norm);
    527:       PetscPrintf(dmmg[level]->comm,"%D FAS function norm %g\n",i,norm);
    528:     }
    529:   }
    530:   theend:
    531:   return(0);
    532: }
    533: #endif

    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_jacobian_period <p> - Indicates how often in the SNES solve the Jacobian is recomputed (on all levels)
    559:                                  as suggested by Florin Dobrian if p is -1 then Jacobian is computed only on first
    560:                                  SNES iteration (i.e. -1 is equivalent to infinity)

    562:     Level: advanced

    564: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNESLocal()

    566: @*/
    567: PetscErrorCode PETSCSNES_DLLEXPORT DMMGSetSNES(DMMG *dmmg,PetscErrorCode (*function)(SNES,Vec,Vec,void*),PetscErrorCode (*jacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*))
    568: {
    570:   PetscMPIInt    size;
    571:   PetscInt       i,nlevels = dmmg[0]->nlevels,period = 1;
    572:   PetscTruth     snesmonitor,mffdoperator,mffd,fdjacobian;
    573: #if defined(PETSC_HAVE_ADIC)
    574:   PetscTruth     mfadoperator,mfad,adjacobian;
    575: #endif
    576:   PetscViewer    ascii;
    577:   MPI_Comm       comm;

    580:   if (!dmmg)     SETERRQ(PETSC_ERR_ARG_NULL,"Passing null as DMMG");
    581:   if (!jacobian) jacobian = DMMGComputeJacobianWithFD;

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

    587:     PetscOptionsName("-dmmg_jacobian_fd","Compute sparse Jacobian explicitly with finite differencing","DMMGSetSNES",&fdjacobian);
    588:     if (fdjacobian) jacobian = DMMGComputeJacobianWithFD;
    589: #if defined(PETSC_HAVE_ADIC)
    590:     PetscOptionsName("-dmmg_jacobian_ad","Compute sparse Jacobian explicitly with ADIC (automatic differentiation)","DMMGSetSNES",&adjacobian);
    591:     if (adjacobian) jacobian = DMMGComputeJacobianWithAdic;
    592: #endif

    594:     PetscOptionsTruthGroupBegin("-dmmg_jacobian_mf_fd_operator","Apply Jacobian via matrix free finite differencing","DMMGSetSNES",&mffdoperator);
    595:     PetscOptionsTruthGroupEnd("-dmmg_jacobian_mf_fd","Apply Jacobian via matrix free finite differencing even in computing preconditioner","DMMGSetSNES",&mffd);
    596:     if (mffd) mffdoperator = PETSC_TRUE;
    597: #if defined(PETSC_HAVE_ADIC)
    598:     PetscOptionsTruthGroupBegin("-dmmg_jacobian_mf_ad_operator","Apply Jacobian via matrix free ADIC (automatic differentiation)","DMMGSetSNES",&mfadoperator);
    599:     PetscOptionsTruthGroupEnd("-dmmg_jacobian_mf_ad","Apply Jacobian via matrix free ADIC (automatic differentiation) even in computing preconditioner","DMMGSetSNES",&mfad);
    600:     if (mfad) mfadoperator = PETSC_TRUE;
    601: #endif
    602:   PetscOptionsEnd();

    604:   /* create solvers for each level */
    605:   for (i=0; i<nlevels; i++) {
    606:     SNESCreate(dmmg[i]->comm,&dmmg[i]->snes);
    607:     SNESGetKSP(dmmg[i]->snes,&dmmg[i]->ksp);
    608:     if (snesmonitor) {
    609:       PetscObjectGetComm((PetscObject)dmmg[i]->snes,&comm);
    610:       PetscViewerASCIIOpen(comm,"stdout",&ascii);
    611:       PetscViewerASCIISetTab(ascii,nlevels-i);
    612:       SNESSetMonitor(dmmg[i]->snes,SNESDefaultMonitor,ascii,(PetscErrorCode(*)(void*))PetscViewerDestroy);
    613:     }

    615:     if (mffdoperator) {
    616:       MatCreateSNESMF(dmmg[i]->snes,dmmg[i]->x,&dmmg[i]->J);
    617:       VecDuplicate(dmmg[i]->x,&dmmg[i]->work1);
    618:       VecDuplicate(dmmg[i]->x,&dmmg[i]->work2);
    619:       MatSNESMFSetFunction(dmmg[i]->J,dmmg[i]->work1,function,dmmg[i]);
    620:       if (mffd) {
    621:         dmmg[i]->B = dmmg[i]->J;
    622:         jacobian   = DMMGComputeJacobianWithMF;
    623:       }
    624: #if defined(PETSC_HAVE_ADIC)
    625:     } else if (mfadoperator) {
    626:       MatRegisterDAAD();
    627:       MatCreateDAAD((DA)dmmg[i]->dm,&dmmg[i]->J);
    628:       MatDAADSetCtx(dmmg[i]->J,dmmg[i]->user);
    629:       if (mfad) {
    630:         dmmg[i]->B = dmmg[i]->J;
    631:         jacobian   = DMMGComputeJacobianWithMF;
    632:       }
    633: #endif
    634:     }
    635:
    636:     if (!dmmg[i]->B) {
    637:       MPI_Comm_size(dmmg[i]->comm,&size);
    638:       DMGetMatrix(dmmg[i]->dm,MATAIJ,&dmmg[i]->B);
    639:     }
    640:     if (!dmmg[i]->J) {
    641:       dmmg[i]->J = dmmg[i]->B;
    642:     }

    644:     DMMGSetUpLevel(dmmg,dmmg[i]->ksp,i+1);
    645:
    646:     /*
    647:        if the number of levels is > 1 then we want the coarse solve in the grid sequencing to use LU
    648:        when possible
    649:     */
    650:     if (nlevels > 1 && i == 0) {
    651:       PC         pc;
    652:       KSP        cksp;
    653:       PetscTruth flg1,flg2,flg3;

    655:       KSPGetPC(dmmg[i]->ksp,&pc);
    656:       PCMGGetCoarseSolve(pc,&cksp);
    657:       KSPGetPC(cksp,&pc);
    658:       PetscTypeCompare((PetscObject)pc,PCILU,&flg1);
    659:       PetscTypeCompare((PetscObject)pc,PCSOR,&flg2);
    660:       PetscTypeCompare((PetscObject)pc,PETSC_NULL,&flg3);
    661:       if (flg1 || flg2 || flg3) {
    662:         PCSetType(pc,PCLU);
    663:       }
    664:     }

    666:     dmmg[i]->solve           = DMMGSolveSNES;
    667:     dmmg[i]->computejacobian = jacobian;
    668:     dmmg[i]->computefunction = function;
    669:   }

    671:   if (jacobian == DMMGComputeJacobianWithFD) {
    672:     ISColoring iscoloring;
    673:     for (i=0; i<nlevels; i++) {
    674:       DMGetColoring(dmmg[i]->dm,IS_COLORING_LOCAL,&iscoloring);
    675:       MatFDColoringCreate(dmmg[i]->B,iscoloring,&dmmg[i]->fdcoloring);
    676:       ISColoringDestroy(iscoloring);
    677:       MatFDColoringSetFunction(dmmg[i]->fdcoloring,(PetscErrorCode(*)(void))function,dmmg[i]);
    678:       MatFDColoringSetFromOptions(dmmg[i]->fdcoloring);
    679:     }
    680: #if defined(PETSC_HAVE_ADIC)
    681:   } else if (jacobian == DMMGComputeJacobianWithAdic) {
    682:     for (i=0; i<nlevels; i++) {
    683:       ISColoring iscoloring;
    684:       DMGetColoring(dmmg[i]->dm,IS_COLORING_GHOSTED,&iscoloring);
    685:       MatSetColoring(dmmg[i]->B,iscoloring);
    686:       ISColoringDestroy(iscoloring);
    687:     }
    688: #endif
    689:   }

    691:   for (i=0; i<nlevels; i++) {
    692:     SNESSetJacobian(dmmg[i]->snes,dmmg[i]->J,dmmg[i]->B,DMMGComputeJacobian_Multigrid,dmmg);
    693:     SNESSetFunction(dmmg[i]->snes,dmmg[i]->b,function,dmmg[i]);
    694:     SNESSetFromOptions(dmmg[i]->snes);
    695:   }

    697:   /* Create interpolation scaling */
    698:   for (i=1; i<nlevels; i++) {
    699:     DMGetInterpolationScale(dmmg[i-1]->dm,dmmg[i]->dm,dmmg[i]->R,&dmmg[i]->Rscale);
    700:   }

    702:   PetscOptionsGetInt(PETSC_NULL,"-dmmg_jacobian_period",&period,PETSC_NULL);
    703:   for (i=0; i<nlevels; i++) {
    704:     dmmg[i]->updatejacobian       = PETSC_TRUE;
    705:     dmmg[i]->updatejacobianperiod = period;
    706:   }

    708: #if defined(PETSC_HAVE_ADIC)
    709:   {
    710:     PetscTruth flg;
    711:     PetscOptionsHasName(PETSC_NULL,"-dmmg_fas",&flg);
    712:     if (flg) {
    713:       PetscInt newton_its;
    714:       PetscOptionsHasName(0,"-dmmg_fas_view",&flg);
    715:       for (i=0; i<nlevels; i++) {
    716:         NLFCreate_DAAD(&dmmg[i]->nlf);
    717:         NLFDAADSetDA_DAAD(dmmg[i]->nlf,(DA)dmmg[i]->dm);
    718:         NLFDAADSetCtx_DAAD(dmmg[i]->nlf,dmmg[i]->user);
    719:         NLFDAADSetResidual_DAAD(dmmg[i]->nlf,dmmg[i]->r);
    720:         VecDuplicate(dmmg[i]->b,&dmmg[i]->w);

    722:         dmmg[i]->monitor    = PETSC_FALSE;
    723:         PetscOptionsHasName(0,"-dmmg_fas_monitor",&dmmg[i]->monitor);
    724:         dmmg[i]->monitorall = PETSC_FALSE;
    725:         PetscOptionsHasName(0,"-dmmg_fas_monitor_all",&dmmg[i]->monitorall);
    726:         dmmg[i]->presmooth  = 2;
    727:         PetscOptionsGetInt(0,"-dmmg_fas_presmooth",&dmmg[i]->presmooth,0);
    728:         dmmg[i]->postsmooth = 2;
    729:         PetscOptionsGetInt(0,"-dmmg_fas_postsmooth",&dmmg[i]->postsmooth,0);
    730:         dmmg[i]->coarsesmooth = 2;
    731:         PetscOptionsGetInt(0,"-dmmg_fas_coarsesmooth",&dmmg[i]->coarsesmooth,0);

    733:         dmmg[i]->rtol = 1.e-8;
    734:         PetscOptionsGetReal(0,"-dmmg_fas_rtol",&dmmg[i]->rtol,0);
    735:         dmmg[i]->abstol = 1.e-50;
    736:         PetscOptionsGetReal(0,"-dmmg_fas_atol",&dmmg[i]->abstol,0);

    738:         newton_its = 2;
    739:         PetscOptionsGetInt(0,"-dmmg_fas_newton_its",&newton_its,0);
    740:         NLFDAADSetNewtonIterations_DAAD(dmmg[i]->nlf,newton_its);

    742:         if (flg) {
    743:           if (i == 0) {
    744:             PetscPrintf(dmmg[i]->comm,"FAS Solver Parameters\n");
    745:             PetscPrintf(dmmg[i]->comm,"  rtol %g atol %g\n",dmmg[i]->rtol,dmmg[i]->abstol);
    746:             PetscPrintf(dmmg[i]->comm,"             coarsesmooths %D\n",dmmg[i]->coarsesmooth);
    747:             PetscPrintf(dmmg[i]->comm,"             Newton iterations %D\n",newton_its);
    748:           } else {
    749:             PetscPrintf(dmmg[i]->comm,"  level %D   presmooths    %D\n",i,dmmg[i]->presmooth);
    750:             PetscPrintf(dmmg[i]->comm,"             postsmooths   %D\n",dmmg[i]->postsmooth);
    751:             PetscPrintf(dmmg[i]->comm,"             Newton iterations %D\n",newton_its);
    752:           }
    753:         }
    754:         dmmg[i]->solve = DMMGSolveFAS;
    755:       }
    756:     }
    757:   }
    758: #endif
    759:
    760:   return(0);
    761: }

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

    767:     Collective on DMMG

    769:    Synopsis:
    770:    PetscErrorCode DMMGSetSNESLocal(DMMG *dmmg,DALocalFunction1 function, DALocalFunction1 jacobian,
    771:                         DALocalFunction1 ad_function, DALocalFunction1 admf_function);

    773:     Input Parameter:
    774: +   dmmg - the context
    775: .   function - the function that defines the nonlinear system
    776: .   jacobian - function defines the local part of the Jacobian
    777: .   ad_function - the name of the function with an ad_ prefix. This is ignored if ADIC is
    778:                   not installed
    779: -   admf_function - the name of the function with an ad_ prefix. This is ignored if ADIC is
    780:                   not installed

    782:     Options Database Keys:
    783: +    -dmmg_snes_monitor
    784: .    -dmmg_jacobian_fd
    785: .    -dmmg_jacobian_ad
    786: .    -dmmg_jacobian_mf_fd_operator
    787: .    -dmmg_jacobian_mf_fd
    788: .    -dmmg_jacobian_mf_ad_operator
    789: .    -dmmg_jacobian_mf_ad
    790: -    -dmmg_jacobian_period <p> - Indicates how often in the SNES solve the Jacobian is recomputed (on all levels)
    791:                                  as suggested by Florin Dobrian if p is -1 then Jacobian is computed only on first
    792:                                  SNES iteration (i.e. -1 is equivalent to infinity)

    795:     Level: intermediate

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

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

    805: .seealso DMMGCreate(), DMMGDestroy, DMMGSetKSP(), DMMGSetSNES()

    807: M*/

    811: PetscErrorCode DMMGSetSNESLocal_Private(DMMG *dmmg,DALocalFunction1 function,DALocalFunction1 jacobian,DALocalFunction1 ad_function,DALocalFunction1 admf_function)
    812: {
    814:   PetscInt       i,nlevels = dmmg[0]->nlevels;
    815:   PetscErrorCode (*computejacobian)(SNES,Vec,Mat*,Mat*,MatStructure*,void*) = 0;

    819:   if (jacobian)         computejacobian = DMMGComputeJacobian;
    820: #if defined(PETSC_HAVE_ADIC)
    821:   else if (ad_function) computejacobian = DMMGComputeJacobianWithAdic;
    822: #endif

    824:   DMMGSetSNES(dmmg,DMMGFormFunction,computejacobian);
    825:   for (i=0; i<nlevels; i++) {
    826:     DASetLocalFunction((DA)dmmg[i]->dm,function);
    827:     DASetLocalJacobian((DA)dmmg[i]->dm,jacobian);
    828:     DASetLocalAdicFunction((DA)dmmg[i]->dm,ad_function);
    829:     DASetLocalAdicMFFunction((DA)dmmg[i]->dm,admf_function);
    830:   }
    831:   return(0);
    832: }

    836: static PetscErrorCode DMMGFunctioni(PetscInt i,Vec u,PetscScalar* r,void* ctx)
    837: {
    838:   DMMG           dmmg = (DMMG)ctx;
    839:   Vec            U = dmmg->lwork1;
    841:   VecScatter     gtol;

    844:   /* copy u into interior part of U */
    845:   DAGetScatter((DA)dmmg->dm,0,&gtol,0);
    846:   VecScatterBegin(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);
    847:   VecScatterEnd(u,U,INSERT_VALUES,SCATTER_FORWARD_LOCAL,gtol);
    848:   DAFormFunctioni1((DA)dmmg->dm,i,U,r,dmmg->user);
    849:   return(0);
    850: }

    854: static PetscErrorCode DMMGFunctioniBase(Vec u,void* ctx)
    855: {
    856:   DMMG           dmmg = (DMMG)ctx;
    857:   Vec            U = dmmg->lwork1;

    861:   DAGlobalToLocalBegin((DA)dmmg->dm,u,INSERT_VALUES,U);
    862:   DAGlobalToLocalEnd((DA)dmmg->dm,u,INSERT_VALUES,U);
    863:   return(0);
    864: }

    868: 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*))
    869: {
    871:   PetscInt       i,nlevels = dmmg[0]->nlevels;

    874:   for (i=0; i<nlevels; i++) {
    875:     DASetLocalFunctioni((DA)dmmg[i]->dm,functioni);
    876:     DASetLocalAdicFunctioni((DA)dmmg[i]->dm,adi);
    877:     DASetLocalAdicMFFunctioni((DA)dmmg[i]->dm,adimf);
    878:     MatSNESMFSetFunctioni(dmmg[i]->J,DMMGFunctioni);
    879:     MatSNESMFSetFunctioniBase(dmmg[i]->J,DMMGFunctioniBase);
    880:     DACreateLocalVector((DA)dmmg[i]->dm,&dmmg[i]->lwork1);
    881:   }
    882:   return(0);
    883: }

    886: #if defined(PETSC_HAVE_ADIC)
    888: #include "adic/ad_utils.h"

    893: PetscErrorCode PetscADView(PetscInt N,PetscInt nc,double *ptr,PetscViewer viewer)
    894: {
    895:   PetscInt       i,j,nlen  = PetscADGetDerivTypeSize();
    896:   char           *cptr = (char*)ptr;
    897:   double         *values;

    901:   for (i=0; i<N; i++) {
    902:     PetscPrintf(PETSC_COMM_SELF,"Element %D value %g derivatives: ",i,*(double*)cptr);
    903:     values = PetscADGetGradArray(cptr);
    904:     for (j=0; j<nc; j++) {
    905:       PetscPrintf(PETSC_COMM_SELF,"%g ",*values++);
    906:     }
    907:     PetscPrintf(PETSC_COMM_SELF,"\n");
    908:     cptr += nlen;
    909:   }

    911:   return(0);
    912: }

    914: #endif