Actual source code: nasm.c

      1: #include <petsc/private/snesimpl.h>
      2: #include <petscdm.h>

      4: typedef struct {
      5:   PetscInt   n;                   /* local subdomains */
      6:   SNES       *subsnes;            /* nonlinear solvers for each subdomain */
      7:   Vec        *x;                  /* solution vectors */
      8:   Vec        *xl;                 /* solution local vectors */
      9:   Vec        *y;                  /* step vectors */
     10:   Vec        *b;                  /* rhs vectors */
     11:   Vec        weight;              /* weighting for adding updates on overlaps, in global space */
     12:   VecScatter *oscatter;           /* scatter from global space to the subdomain global space */
     13:   VecScatter *oscatter_copy;      /* copy of the above */
     14:   VecScatter *iscatter;           /* scatter from global space to the nonoverlapping subdomain space */
     15:   VecScatter *gscatter;           /* scatter from global space to the subdomain local space */
     16:   PCASMType  type;                /* ASM type */
     17:   PetscBool  usesdm;              /* use the DM for setting up the subproblems */
     18:   PetscBool  finaljacobian;       /* compute the jacobian of the converged solution */
     19:   PetscReal  damping;             /* damping parameter for updates from the blocks */
     20:   PetscBool  weight_set;          /* use a weight in the overlap updates */

     22:   /* logging events */
     23:   PetscLogEvent eventrestrictinterp;
     24:   PetscLogEvent eventsubsolve;

     26:   PetscInt      fjtype;            /* type of computed jacobian */
     27:   Vec           xinit;             /* initial solution in case the final jacobian type is computed as first */
     28: } SNES_NASM;

     30: const char *const SNESNASMTypes[] = {"NONE","RESTRICT","INTERPOLATE","BASIC","PCASMType","PC_ASM_",NULL};
     31: const char *const SNESNASMFJTypes[] = {"FINALOUTER","FINALINNER","INITIAL"};

     33: static PetscErrorCode SNESReset_NASM(SNES snes)
     34: {
     35:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
     37:   PetscInt       i;

     40:   for (i=0; i<nasm->n; i++) {
     41:     if (nasm->xl) { VecDestroy(&nasm->xl[i]); }
     42:     if (nasm->x) { VecDestroy(&nasm->x[i]); }
     43:     if (nasm->y) { VecDestroy(&nasm->y[i]); }
     44:     if (nasm->b) { VecDestroy(&nasm->b[i]); }

     46:     if (nasm->subsnes) { SNESDestroy(&nasm->subsnes[i]); }
     47:     if (nasm->oscatter) { VecScatterDestroy(&nasm->oscatter[i]); }
     48:     if (nasm->oscatter_copy) { VecScatterDestroy(&nasm->oscatter_copy[i]); }
     49:     if (nasm->iscatter) { VecScatterDestroy(&nasm->iscatter[i]); }
     50:     if (nasm->gscatter) { VecScatterDestroy(&nasm->gscatter[i]); }
     51:   }

     53:   PetscFree(nasm->x);
     54:   PetscFree(nasm->xl);
     55:   PetscFree(nasm->y);
     56:   PetscFree(nasm->b);

     58:   if (nasm->xinit) {VecDestroy(&nasm->xinit);}

     60:   PetscFree(nasm->subsnes);
     61:   PetscFree(nasm->oscatter);
     62:   PetscFree(nasm->oscatter_copy);
     63:   PetscFree(nasm->iscatter);
     64:   PetscFree(nasm->gscatter);

     66:   if (nasm->weight_set) {
     67:     VecDestroy(&nasm->weight);
     68:   }

     70:   nasm->eventrestrictinterp = 0;
     71:   nasm->eventsubsolve = 0;
     72:   return(0);
     73: }

     75: static PetscErrorCode SNESDestroy_NASM(SNES snes)
     76: {

     80:   SNESReset_NASM(snes);
     81:   PetscFree(snes->data);
     82:   return(0);
     83: }

     85: static PetscErrorCode DMGlobalToLocalSubDomainDirichletHook_Private(DM dm,Vec g,InsertMode mode,Vec l,void *ctx)
     86: {
     88:   Vec            bcs = (Vec)ctx;

     91:   VecCopy(bcs,l);
     92:   return(0);
     93: }

     95: static PetscErrorCode SNESSetUp_NASM(SNES snes)
     96: {
     97:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
     99:   DM             dm,subdm;
    100:   DM             *subdms;
    101:   PetscInt       i;
    102:   const char     *optionsprefix;
    103:   Vec            F;
    104:   PetscMPIInt    size;
    105:   KSP            ksp;
    106:   PC             pc;

    109:   if (!nasm->subsnes) {
    110:     SNESGetDM(snes,&dm);
    111:     if (dm) {
    112:       nasm->usesdm = PETSC_TRUE;
    113:       DMCreateDomainDecomposition(dm,&nasm->n,NULL,NULL,NULL,&subdms);
    114:       if (!subdms) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM has no default decomposition defined.  Set subsolves manually with SNESNASMSetSubdomains().");
    115:       DMCreateDomainDecompositionScatters(dm,nasm->n,subdms,&nasm->iscatter,&nasm->oscatter,&nasm->gscatter);
    116:       PetscMalloc1(nasm->n, &nasm->oscatter_copy);
    117:       for (i=0; i<nasm->n; i++) {
    118:         VecScatterCopy(nasm->oscatter[i], &nasm->oscatter_copy[i]);
    119:       }

    121:       SNESGetOptionsPrefix(snes, &optionsprefix);
    122:       PetscMalloc1(nasm->n,&nasm->subsnes);
    123:       for (i=0; i<nasm->n; i++) {
    124:         SNESCreate(PETSC_COMM_SELF,&nasm->subsnes[i]);
    125:         PetscObjectIncrementTabLevel((PetscObject)nasm->subsnes[i], (PetscObject)snes, 1);
    126:         SNESAppendOptionsPrefix(nasm->subsnes[i],optionsprefix);
    127:         SNESAppendOptionsPrefix(nasm->subsnes[i],"sub_");
    128:         SNESSetDM(nasm->subsnes[i],subdms[i]);
    129:         MPI_Comm_size(PetscObjectComm((PetscObject)nasm->subsnes[i]),&size);
    130:         if (size == 1) {
    131:           SNESGetKSP(nasm->subsnes[i],&ksp);
    132:           KSPGetPC(ksp,&pc);
    133:           KSPSetType(ksp,KSPPREONLY);
    134:           PCSetType(pc,PCLU);
    135:         }
    136:         SNESSetFromOptions(nasm->subsnes[i]);
    137:         DMDestroy(&subdms[i]);
    138:       }
    139:       PetscFree(subdms);
    140:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Cannot construct local problems automatically without a DM!");
    141:   } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Must set subproblems manually if there is no DM!");
    142:   /* allocate the global vectors */
    143:   if (!nasm->x) {
    144:     PetscCalloc1(nasm->n,&nasm->x);
    145:   }
    146:   if (!nasm->xl) {
    147:     PetscCalloc1(nasm->n,&nasm->xl);
    148:   }
    149:   if (!nasm->y) {
    150:     PetscCalloc1(nasm->n,&nasm->y);
    151:   }
    152:   if (!nasm->b) {
    153:     PetscCalloc1(nasm->n,&nasm->b);
    154:   }

    156:   for (i=0; i<nasm->n; i++) {
    157:     SNESGetFunction(nasm->subsnes[i],&F,NULL,NULL);
    158:     if (!nasm->x[i]) {VecDuplicate(F,&nasm->x[i]);}
    159:     if (!nasm->y[i]) {VecDuplicate(F,&nasm->y[i]);}
    160:     if (!nasm->b[i]) {VecDuplicate(F,&nasm->b[i]);}
    161:     if (!nasm->xl[i]) {
    162:       SNESGetDM(nasm->subsnes[i],&subdm);
    163:       DMCreateLocalVector(subdm,&nasm->xl[i]);
    164:       DMGlobalToLocalHookAdd(subdm,DMGlobalToLocalSubDomainDirichletHook_Private,NULL,nasm->xl[i]);
    165:     }
    166:   }
    167:   if (nasm->finaljacobian) {
    168:     SNESSetUpMatrices(snes);
    169:     if (nasm->fjtype == 2) {
    170:       VecDuplicate(snes->vec_sol,&nasm->xinit);
    171:     }
    172:     for (i=0; i<nasm->n;i++) {
    173:       SNESSetUpMatrices(nasm->subsnes[i]);
    174:     }
    175:   }
    176:   return(0);
    177: }

    179: static PetscErrorCode SNESSetFromOptions_NASM(PetscOptionItems *PetscOptionsObject,SNES snes)
    180: {
    181:   PetscErrorCode    ierr;
    182:   PCASMType         asmtype;
    183:   PetscBool         flg,monflg;
    184:   SNES_NASM         *nasm = (SNES_NASM*)snes->data;

    187:   PetscOptionsHead(PetscOptionsObject,"Nonlinear Additive Schwarz options");
    188:   PetscOptionsEnum("-snes_nasm_type","Type of restriction/extension","",SNESNASMTypes,(PetscEnum)nasm->type,(PetscEnum*)&asmtype,&flg);
    189:   if (flg) {SNESNASMSetType(snes,asmtype);}
    190:   flg    = PETSC_FALSE;
    191:   monflg = PETSC_TRUE;
    192:   PetscOptionsReal("-snes_nasm_damping","The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)","SNESNASMSetDamping",nasm->damping,&nasm->damping,&flg);
    193:   if (flg) {SNESNASMSetDamping(snes,nasm->damping);}
    194:   PetscOptionsDeprecated("-snes_nasm_sub_view",NULL,"3.15","Use -snes_view ::ascii_info_detail");
    195:   PetscOptionsBool("-snes_nasm_finaljacobian","Compute the global jacobian of the final iterate (for ASPIN)","",nasm->finaljacobian,&nasm->finaljacobian,NULL);
    196:   PetscOptionsEList("-snes_nasm_finaljacobian_type","The type of the final jacobian computed.","",SNESNASMFJTypes,3,SNESNASMFJTypes[0],&nasm->fjtype,NULL);
    197:   PetscOptionsBool("-snes_nasm_log","Log times for subSNES solves and restriction","",monflg,&monflg,&flg);
    198:   if (flg) {
    199:     PetscLogEventRegister("SNESNASMSubSolve",((PetscObject)snes)->classid,&nasm->eventsubsolve);
    200:     PetscLogEventRegister("SNESNASMRestrict",((PetscObject)snes)->classid,&nasm->eventrestrictinterp);
    201:   }
    202:   PetscOptionsTail();
    203:   return(0);
    204: }

    206: static PetscErrorCode SNESView_NASM(SNES snes, PetscViewer viewer)
    207: {
    208:   SNES_NASM         *nasm = (SNES_NASM*)snes->data;
    209:   PetscErrorCode    ierr;
    210:   PetscMPIInt       rank,size;
    211:   PetscInt          i,N,bsz;
    212:   PetscBool         iascii,isstring;
    213:   PetscViewer       sviewer;
    214:   MPI_Comm          comm;
    215:   PetscViewerFormat format;
    216:   const char        *prefix;

    219:   PetscObjectGetComm((PetscObject)snes,&comm);
    220:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    221:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
    222:   MPI_Comm_rank(comm,&rank);
    223:   MPI_Comm_size(comm,&size);
    224:   MPIU_Allreduce(&nasm->n,&N,1,MPIU_INT,MPI_SUM,comm);
    225:   if (iascii) {
    226:     PetscViewerASCIIPrintf(viewer, "  total subdomain blocks = %D\n",N);
    227:     PetscViewerGetFormat(viewer,&format);
    228:     if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
    229:       if (nasm->subsnes) {
    230:         PetscViewerASCIIPrintf(viewer,"  Local solver information for first block on rank 0:\n");
    231:         SNESGetOptionsPrefix(snes,&prefix);
    232:         PetscViewerASCIIPrintf(viewer,"  Use -%ssnes_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
    233:         PetscViewerASCIIPushTab(viewer);
    234:         PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    235:         if (!rank) {
    236:           PetscViewerASCIIPushTab(viewer);
    237:           SNESView(nasm->subsnes[0],sviewer);
    238:           PetscViewerASCIIPopTab(viewer);
    239:         }
    240:         PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    241:         PetscViewerASCIIPopTab(viewer);
    242:       }
    243:     } else {
    244:       /* print the solver on each block */
    245:       PetscViewerASCIIPushSynchronized(viewer);
    246:       PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] number of local blocks = %D\n",(int)rank,nasm->n);
    247:       PetscViewerFlush(viewer);
    248:       PetscViewerASCIIPopSynchronized(viewer);
    249:       PetscViewerASCIIPrintf(viewer,"  Local solver information for each block is in the following SNES objects:\n");
    250:       PetscViewerASCIIPushTab(viewer);
    251:       PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
    252:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    253:       for (i=0; i<nasm->n; i++) {
    254:         VecGetLocalSize(nasm->x[i],&bsz);
    255:         PetscViewerASCIIPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
    256:         SNESView(nasm->subsnes[i],sviewer);
    257:         PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
    258:       }
    259:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    260:       PetscViewerFlush(viewer);
    261:       PetscViewerASCIIPopTab(viewer);
    262:     }
    263:   } else if (isstring) {
    264:     PetscViewerStringSPrintf(viewer," blocks=%D,type=%s",N,SNESNASMTypes[nasm->type]);
    265:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    266:     if (nasm->subsnes && !rank) {SNESView(nasm->subsnes[0],sviewer);}
    267:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
    268:   }
    269:   return(0);
    270: }

    272: /*@
    273:    SNESNASMSetType - Set the type of subdomain update used

    275:    Logically Collective on SNES

    277:    Input Parameters:
    278: +  SNES - the SNES context
    279: -  type - the type of update, PC_ASM_BASIC or PC_ASM_RESTRICT

    281:    Level: intermediate

    283: .seealso: SNESNASM, SNESNASMGetType(), PCASMSetType()
    284: @*/
    285: PetscErrorCode SNESNASMSetType(SNES snes,PCASMType type)
    286: {
    288:   PetscErrorCode (*f)(SNES,PCASMType);

    291:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetType_C",&f);
    292:   if (f) {(f)(snes,type);}
    293:   return(0);
    294: }

    296: static PetscErrorCode SNESNASMSetType_NASM(SNES snes,PCASMType type)
    297: {
    298:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    301:   if (type != PC_ASM_BASIC && type != PC_ASM_RESTRICT) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"SNESNASM only supports basic and restrict types");
    302:   nasm->type = type;
    303:   return(0);
    304: }

    306: /*@
    307:    SNESNASMGetType - Get the type of subdomain update used

    309:    Logically Collective on SNES

    311:    Input Parameters:
    312: .  SNES - the SNES context

    314:    Output Parameters:
    315: .  type - the type of update

    317:    Level: intermediate

    319: .seealso: SNESNASM, SNESNASMSetType(), PCASMGetType()
    320: @*/
    321: PetscErrorCode SNESNASMGetType(SNES snes,PCASMType *type)
    322: {

    326:   PetscUseMethod(snes,"SNESNASMGetType_C",(SNES,PCASMType*),(snes,type));
    327:   return(0);
    328: }

    330: static PetscErrorCode SNESNASMGetType_NASM(SNES snes,PCASMType *type)
    331: {
    332:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    335:   *type = nasm->type;
    336:   return(0);
    337: }

    339: /*@
    340:    SNESNASMSetSubdomains - Manually Set the context required to restrict and solve subdomain problems.

    342:    Not Collective

    344:    Input Parameters:
    345: +  SNES - the SNES context
    346: .  n - the number of local subdomains
    347: .  subsnes - solvers defined on the local subdomains
    348: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
    349: .  oscatter - scatters into the overlapping portions of the local subdomains
    350: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

    352:    Level: intermediate

    354: .seealso: SNESNASM, SNESNASMGetSubdomains()
    355: @*/
    356: PetscErrorCode SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
    357: {
    359:   PetscErrorCode (*f)(SNES,PetscInt,SNES*,VecScatter*,VecScatter*,VecScatter*);

    362:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",&f);
    363:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
    364:   return(0);
    365: }

    367: static PetscErrorCode SNESNASMSetSubdomains_NASM(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[])
    368: {
    369:   PetscInt       i;
    371:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    374:   if (snes->setupcalled) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESNASMSetSubdomains() should be called before calling SNESSetUp().");

    376:   /* tear down the previously set things */
    377:   SNESReset(snes);

    379:   nasm->n = n;
    380:   if (oscatter) {
    381:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)oscatter[i]);}
    382:   }
    383:   if (iscatter) {
    384:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)iscatter[i]);}
    385:   }
    386:   if (gscatter) {
    387:     for (i=0; i<n; i++) {PetscObjectReference((PetscObject)gscatter[i]);}
    388:   }
    389:   if (oscatter) {
    390:     PetscMalloc1(n,&nasm->oscatter);
    391:     PetscMalloc1(n,&nasm->oscatter_copy);
    392:     for (i=0; i<n; i++) {
    393:       nasm->oscatter[i] = oscatter[i];
    394:       VecScatterCopy(oscatter[i], &nasm->oscatter_copy[i]);
    395:     }
    396:   }
    397:   if (iscatter) {
    398:     PetscMalloc1(n,&nasm->iscatter);
    399:     for (i=0; i<n; i++) {
    400:       nasm->iscatter[i] = iscatter[i];
    401:     }
    402:   }
    403:   if (gscatter) {
    404:     PetscMalloc1(n,&nasm->gscatter);
    405:     for (i=0; i<n; i++) {
    406:       nasm->gscatter[i] = gscatter[i];
    407:     }
    408:   }

    410:   if (subsnes) {
    411:     PetscMalloc1(n,&nasm->subsnes);
    412:     for (i=0; i<n; i++) {
    413:       nasm->subsnes[i] = subsnes[i];
    414:     }
    415:   }
    416:   return(0);
    417: }

    419: /*@
    420:    SNESNASMGetSubdomains - Get the local subdomain context.

    422:    Not Collective

    424:    Input Parameters:
    425: .  SNES - the SNES context

    427:    Output Parameters:
    428: +  n - the number of local subdomains
    429: .  subsnes - solvers defined on the local subdomains
    430: .  iscatter - scatters into the nonoverlapping portions of the local subdomains
    431: .  oscatter - scatters into the overlapping portions of the local subdomains
    432: -  gscatter - scatters into the (ghosted) local vector of the local subdomain

    434:    Level: intermediate

    436: .seealso: SNESNASM, SNESNASMSetSubdomains()
    437: @*/
    438: PetscErrorCode SNESNASMGetSubdomains(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
    439: {
    441:   PetscErrorCode (*f)(SNES,PetscInt*,SNES**,VecScatter**,VecScatter**,VecScatter**);

    444:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",&f);
    445:   if (f) {(f)(snes,n,subsnes,iscatter,oscatter,gscatter);}
    446:   return(0);
    447: }

    449: static PetscErrorCode SNESNASMGetSubdomains_NASM(SNES snes,PetscInt *n,SNES *subsnes[],VecScatter *iscatter[],VecScatter *oscatter[],VecScatter *gscatter[])
    450: {
    451:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    454:   if (n) *n = nasm->n;
    455:   if (oscatter) *oscatter = nasm->oscatter;
    456:   if (iscatter) *iscatter = nasm->iscatter;
    457:   if (gscatter) *gscatter = nasm->gscatter;
    458:   if (subsnes)  *subsnes  = nasm->subsnes;
    459:   return(0);
    460: }

    462: /*@
    463:    SNESNASMGetSubdomainVecs - Get the processor-local subdomain vectors

    465:    Not Collective

    467:    Input Parameters:
    468: .  SNES - the SNES context

    470:    Output Parameters:
    471: +  n - the number of local subdomains
    472: .  x - The subdomain solution vector
    473: .  y - The subdomain step vector
    474: .  b - The subdomain RHS vector
    475: -  xl - The subdomain local vectors (ghosted)

    477:    Level: developer

    479: .seealso: SNESNASM, SNESNASMGetSubdomains()
    480: @*/
    481: PetscErrorCode SNESNASMGetSubdomainVecs(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b, Vec **xl)
    482: {
    484:   PetscErrorCode (*f)(SNES,PetscInt*,Vec**,Vec**,Vec**,Vec**);

    487:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",&f);
    488:   if (f) {(f)(snes,n,x,y,b,xl);}
    489:   return(0);
    490: }

    492: static PetscErrorCode SNESNASMGetSubdomainVecs_NASM(SNES snes,PetscInt *n,Vec **x,Vec **y,Vec **b,Vec **xl)
    493: {
    494:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    497:   if (n)  *n  = nasm->n;
    498:   if (x)  *x  = nasm->x;
    499:   if (y)  *y  = nasm->y;
    500:   if (b)  *b  = nasm->b;
    501:   if (xl) *xl = nasm->xl;
    502:   return(0);
    503: }

    505: /*@
    506:    SNESNASMSetComputeFinalJacobian - Schedules the computation of the global and subdomain Jacobians upon convergence

    508:    Collective on SNES

    510:    Input Parameters:
    511: +  SNES - the SNES context
    512: -  flg - indication of whether to compute the Jacobians or not

    514:    Level: developer

    516:    Notes:
    517:    This is used almost exclusively in the implementation of ASPIN, where the converged subdomain and global Jacobian
    518:    is needed at each linear iteration.

    520: .seealso: SNESNASM, SNESNASMGetSubdomains()
    521: @*/
    522: PetscErrorCode SNESNASMSetComputeFinalJacobian(SNES snes,PetscBool flg)
    523: {
    524:   PetscErrorCode (*f)(SNES,PetscBool);

    528:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",&f);
    529:   if (f) {(f)(snes,flg);}
    530:   return(0);
    531: }

    533: static PetscErrorCode SNESNASMSetComputeFinalJacobian_NASM(SNES snes,PetscBool flg)
    534: {
    535:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    538:   nasm->finaljacobian = flg;
    539:   return(0);
    540: }

    542: /*@
    543:    SNESNASMSetDamping - Sets the update damping for NASM

    545:    Logically collective on SNES

    547:    Input Parameters:
    548: +  SNES - the SNES context
    549: -  dmp - damping

    551:    Level: intermediate

    553:    Notes:
    554:     The new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)

    556: .seealso: SNESNASM, SNESNASMGetDamping()
    557: @*/
    558: PetscErrorCode SNESNASMSetDamping(SNES snes,PetscReal dmp)
    559: {
    560:   PetscErrorCode (*f)(SNES,PetscReal);

    564:   PetscObjectQueryFunction((PetscObject)snes,"SNESNASMSetDamping_C",(void (**)(void))&f);
    565:   if (f) {(f)(snes,dmp);}
    566:   return(0);
    567: }

    569: static PetscErrorCode SNESNASMSetDamping_NASM(SNES snes,PetscReal dmp)
    570: {
    571:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    574:   nasm->damping = dmp;
    575:   return(0);
    576: }

    578: /*@
    579:    SNESNASMGetDamping - Gets the update damping for NASM

    581:    Not Collective

    583:    Input Parameters:
    584: +  SNES - the SNES context
    585: -  dmp - damping

    587:    Level: intermediate

    589: .seealso: SNESNASM, SNESNASMSetDamping()
    590: @*/
    591: PetscErrorCode SNESNASMGetDamping(SNES snes,PetscReal *dmp)
    592: {

    596:   PetscUseMethod(snes,"SNESNASMGetDamping_C",(SNES,PetscReal*),(snes,dmp));
    597:   return(0);
    598: }

    600: static PetscErrorCode SNESNASMGetDamping_NASM(SNES snes,PetscReal *dmp)
    601: {
    602:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    605:   *dmp = nasm->damping;
    606:   return(0);
    607: }

    610: /*
    611:   Input Parameters:
    612: + snes - The solver
    613: . B - The RHS vector
    614: - X - The initial guess

    616:   Output Parameters:
    617: . Y - The solution update

    619:   TODO: All scatters should be packed into one
    620: */
    621: PetscErrorCode SNESNASMSolveLocal_Private(SNES snes,Vec B,Vec Y,Vec X)
    622: {
    623:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
    624:   SNES           subsnes;
    625:   PetscInt       i;
    626:   PetscReal      dmp;
    628:   Vec            Xl,Bl,Yl,Xlloc;
    629:   VecScatter     iscat,oscat,gscat,oscat_copy;
    630:   DM             dm,subdm;
    631:   PCASMType      type;

    634:   SNESNASMGetType(snes,&type);
    635:   SNESGetDM(snes,&dm);
    636:   VecSet(Y,0);
    637:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
    638:   for (i=0; i<nasm->n; i++) {
    639:     /* scatter the solution to the global solution and the local solution */
    640:     Xl      = nasm->x[i];
    641:     Xlloc   = nasm->xl[i];
    642:     oscat   = nasm->oscatter[i];
    643:     oscat_copy = nasm->oscatter_copy[i];
    644:     gscat   = nasm->gscatter[i];
    645:     VecScatterBegin(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
    646:     VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
    647:     if (B) {
    648:       /* scatter the RHS to the local RHS */
    649:       Bl   = nasm->b[i];
    650:       VecScatterBegin(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
    651:     }
    652:   }
    653:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}

    656:   if (nasm->eventsubsolve) {PetscLogEventBegin(nasm->eventsubsolve,snes,0,0,0);}
    657:   for (i=0; i<nasm->n; i++) {
    658:     Xl    = nasm->x[i];
    659:     Xlloc = nasm->xl[i];
    660:     Yl    = nasm->y[i];
    661:     subsnes = nasm->subsnes[i];
    662:     SNESGetDM(subsnes,&subdm);
    663:     iscat   = nasm->iscatter[i];
    664:     oscat   = nasm->oscatter[i];
    665:     oscat_copy = nasm->oscatter_copy[i];
    666:     gscat   = nasm->gscatter[i];
    667:     VecScatterEnd(oscat,X,Xl,INSERT_VALUES,SCATTER_FORWARD);
    668:     VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
    669:     if (B) {
    670:       Bl   = nasm->b[i];
    671:       VecScatterEnd(oscat_copy,B,Bl,INSERT_VALUES,SCATTER_FORWARD);
    672:     } else Bl = NULL;

    674:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
    675:     VecCopy(Xl,Yl);
    676:     SNESSolve(subsnes,Bl,Xl);
    677:     VecAYPX(Yl,-1.0,Xl);
    678:     VecScale(Yl, nasm->damping);
    679:     if (type == PC_ASM_BASIC) {
    680:       VecScatterBegin(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
    681:       VecScatterEnd(oscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
    682:     } else if (type == PC_ASM_RESTRICT) {
    683:       VecScatterBegin(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
    684:       VecScatterEnd(iscat,Yl,Y,ADD_VALUES,SCATTER_REVERSE);
    685:     } else SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"Only basic and restrict types are supported for SNESNASM");
    686:   }
    687:   if (nasm->eventsubsolve) {PetscLogEventEnd(nasm->eventsubsolve,snes,0,0,0);}
    688:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
    689:   if (nasm->weight_set) {
    690:     VecPointwiseMult(Y,Y,nasm->weight);
    691:   }
    692:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
    693:   SNESNASMGetDamping(snes,&dmp);
    694:   VecAXPY(X,dmp,Y);
    695:   return(0);
    696: }

    698: static PetscErrorCode SNESNASMComputeFinalJacobian_Private(SNES snes, Vec Xfinal)
    699: {
    700:   Vec            X = Xfinal;
    701:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;
    702:   SNES           subsnes;
    703:   PetscInt       i,lag = 1;
    705:   Vec            Xlloc,Xl,Fl,F;
    706:   VecScatter     oscat,gscat;
    707:   DM             dm,subdm;

    710:   if (nasm->fjtype == 2) X = nasm->xinit;
    711:   F = snes->vec_func;
    712:   if (snes->normschedule == SNES_NORM_NONE) {SNESComputeFunction(snes,X,F);}
    713:   SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
    714:   SNESGetDM(snes,&dm);
    715:   if (nasm->eventrestrictinterp) {PetscLogEventBegin(nasm->eventrestrictinterp,snes,0,0,0);}
    716:   if (nasm->fjtype != 1) {
    717:     for (i=0; i<nasm->n; i++) {
    718:       Xlloc = nasm->xl[i];
    719:       gscat = nasm->gscatter[i];
    720:       VecScatterBegin(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);
    721:     }
    722:   }
    723:   if (nasm->eventrestrictinterp) {PetscLogEventEnd(nasm->eventrestrictinterp,snes,0,0,0);}
    724:   for (i=0; i<nasm->n; i++) {
    725:     Fl      = nasm->subsnes[i]->vec_func;
    726:     Xl      = nasm->x[i];
    727:     Xlloc   = nasm->xl[i];
    728:     subsnes = nasm->subsnes[i];
    729:     oscat   = nasm->oscatter[i];
    730:     gscat   = nasm->gscatter[i];
    731:     if (nasm->fjtype != 1) {VecScatterEnd(gscat,X,Xlloc,INSERT_VALUES,SCATTER_FORWARD);}
    732:     SNESGetDM(subsnes,&subdm);
    733:     DMSubDomainRestrict(dm,oscat,gscat,subdm);
    734:     if (nasm->fjtype != 1) {
    735:       DMLocalToGlobalBegin(subdm,Xlloc,INSERT_VALUES,Xl);
    736:       DMLocalToGlobalEnd(subdm,Xlloc,INSERT_VALUES,Xl);
    737:     }
    738:     if (subsnes->lagjacobian == -1)    subsnes->lagjacobian = -2;
    739:     else if (subsnes->lagjacobian > 1) lag = subsnes->lagjacobian;
    740:     SNESComputeFunction(subsnes,Xl,Fl);
    741:     SNESComputeJacobian(subsnes,Xl,subsnes->jacobian,subsnes->jacobian_pre);
    742:     if (lag > 1) subsnes->lagjacobian = lag;
    743:   }
    744:   return(0);
    745: }

    747: static PetscErrorCode SNESSolve_NASM(SNES snes)
    748: {
    749:   Vec              F;
    750:   Vec              X;
    751:   Vec              B;
    752:   Vec              Y;
    753:   PetscInt         i;
    754:   PetscReal        fnorm = 0.0;
    755:   PetscErrorCode   ierr;
    756:   SNESNormSchedule normschedule;
    757:   SNES_NASM        *nasm = (SNES_NASM*)snes->data;

    761:   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

    763:   PetscCitationsRegister(SNESCitation,&SNEScite);
    764:   X = snes->vec_sol;
    765:   Y = snes->vec_sol_update;
    766:   F = snes->vec_func;
    767:   B = snes->vec_rhs;

    769:   PetscObjectSAWsTakeAccess((PetscObject)snes);
    770:   snes->iter   = 0;
    771:   snes->norm   = 0.;
    772:   PetscObjectSAWsGrantAccess((PetscObject)snes);
    773:   snes->reason = SNES_CONVERGED_ITERATING;
    774:   SNESGetNormSchedule(snes, &normschedule);
    775:   if (normschedule == SNES_NORM_ALWAYS || normschedule == SNES_NORM_INITIAL_ONLY || normschedule == SNES_NORM_INITIAL_FINAL_ONLY) {
    776:     /* compute the initial function and preconditioned update delX */
    777:     if (!snes->vec_func_init_set) {
    778:       SNESComputeFunction(snes,X,F);
    779:     } else snes->vec_func_init_set = PETSC_FALSE;

    781:     VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
    782:     SNESCheckFunctionNorm(snes,fnorm);
    783:     PetscObjectSAWsTakeAccess((PetscObject)snes);
    784:     snes->iter = 0;
    785:     snes->norm = fnorm;
    786:     PetscObjectSAWsGrantAccess((PetscObject)snes);
    787:     SNESLogConvergenceHistory(snes,snes->norm,0);
    788:     SNESMonitor(snes,0,snes->norm);

    790:     /* test convergence */
    791:     (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
    792:     if (snes->reason) return(0);
    793:   } else {
    794:     PetscObjectSAWsGrantAccess((PetscObject)snes);
    795:     SNESLogConvergenceHistory(snes,snes->norm,0);
    796:     SNESMonitor(snes,0,snes->norm);
    797:   }

    799:   /* Call general purpose update function */
    800:   if (snes->ops->update) {
    801:     (*snes->ops->update)(snes, snes->iter);
    802:   }
    803:   /* copy the initial solution over for later */
    804:   if (nasm->fjtype == 2) {VecCopy(X,nasm->xinit);}

    806:   for (i=0; i < snes->max_its; i++) {
    807:     SNESNASMSolveLocal_Private(snes,B,Y,X);
    808:     if (normschedule == SNES_NORM_ALWAYS || ((i == snes->max_its - 1) && (normschedule == SNES_NORM_INITIAL_FINAL_ONLY || normschedule == SNES_NORM_FINAL_ONLY))) {
    809:       SNESComputeFunction(snes,X,F);
    810:       VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F||  */
    811:       SNESCheckFunctionNorm(snes,fnorm);
    812:     }
    813:     /* Monitor convergence */
    814:     PetscObjectSAWsTakeAccess((PetscObject)snes);
    815:     snes->iter = i+1;
    816:     snes->norm = fnorm;
    817:     PetscObjectSAWsGrantAccess((PetscObject)snes);
    818:     SNESLogConvergenceHistory(snes,snes->norm,0);
    819:     SNESMonitor(snes,snes->iter,snes->norm);
    820:     /* Test for convergence */
    821:     if (normschedule == SNES_NORM_ALWAYS) {(*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);}
    822:     if (snes->reason) break;
    823:     /* Call general purpose update function */
    824:     if (snes->ops->update) {(*snes->ops->update)(snes, snes->iter);}
    825:   }
    826:   if (nasm->finaljacobian) {
    827:     SNESNASMComputeFinalJacobian_Private(snes,X);
    828:     SNESCheckJacobianDomainerror(snes);
    829:   }
    830:   if (normschedule == SNES_NORM_ALWAYS) {
    831:     if (i == snes->max_its) {
    832:       PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",snes->max_its);
    833:       if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
    834:     }
    835:   } else if (!snes->reason) snes->reason = SNES_CONVERGED_ITS; /* NASM is meant to be used as a preconditioner */
    836:   return(0);
    837: }

    839: /*MC
    840:   SNESNASM - Nonlinear Additive Schwarz

    842:    Options Database:
    843: +  -snes_nasm_log - enable logging events for the communication and solve stages
    844: .  -snes_nasm_type <basic,restrict> - type of subdomain update used
    845: .  -snes_asm_damping <dmp> - the new solution is obtained as old solution plus dmp times (sum of the solutions on the subdomains)
    846: .  -snes_nasm_finaljacobian - compute the local and global jacobians of the final iterate
    847: .  -snes_nasm_finaljacobian_type <finalinner,finalouter,initial> - pick state the jacobian is calculated at
    848: .  -sub_snes_ - options prefix of the subdomain nonlinear solves
    849: .  -sub_ksp_ - options prefix of the subdomain Krylov solver
    850: -  -sub_pc_ - options prefix of the subdomain preconditioner

    852:    Level: advanced

    854:    Developer Note: This is a non-Newton based nonlinear solver that does not directly require a Jacobian; hence the flag snes->usesksp is set to
    855:        false and SNESView() and -snes_view do not display a KSP object. However, if the flag nasm->finaljacobian is set (for example, if
    856:        NASM is used as a nonlinear preconditioner for  KSPASPIN) then SNESSetUpMatrices() is called to generate the Jacobian (needed by KSPASPIN)
    857:        and this utilizes the KSP for storing the matrices, but the KSP is never used for solving a linear system. Note that when SNESNASM is
    858:        used by SNESASPIN they share the same Jacobian matrices because SNESSetUp() (called on the outer SNES KSPASPIN) causes the inner SNES
    859:        object (in this case SNESNASM) to inherit the outer Jacobian matrices.

    861:    References:
    862: .  1. - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
    863:    SIAM Review, 57(4), 2015

    865: .seealso: SNESCreate(), SNES, SNESSetType(), SNESType (for list of available types), SNESNASMSetType(), SNESNASMGetType(), SNESNASMSetSubdomains(), SNESNASMGetSubdomains(), SNESNASMGetSubdomainVecs(), SNESNASMSetComputeFinalJacobian(), SNESNASMSetDamping(), SNESNASMGetDamping()
    866: M*/

    868: PETSC_EXTERN PetscErrorCode SNESCreate_NASM(SNES snes)
    869: {
    870:   SNES_NASM      *nasm;

    874:   PetscNewLog(snes,&nasm);
    875:   snes->data = (void*)nasm;

    877:   nasm->n        = PETSC_DECIDE;
    878:   nasm->subsnes  = NULL;
    879:   nasm->x        = NULL;
    880:   nasm->xl       = NULL;
    881:   nasm->y        = NULL;
    882:   nasm->b        = NULL;
    883:   nasm->oscatter = NULL;
    884:   nasm->oscatter_copy = NULL;
    885:   nasm->iscatter = NULL;
    886:   nasm->gscatter = NULL;
    887:   nasm->damping  = 1.;

    889:   nasm->type              = PC_ASM_BASIC;
    890:   nasm->finaljacobian     = PETSC_FALSE;
    891:   nasm->weight_set        = PETSC_FALSE;

    893:   snes->ops->destroy        = SNESDestroy_NASM;
    894:   snes->ops->setup          = SNESSetUp_NASM;
    895:   snes->ops->setfromoptions = SNESSetFromOptions_NASM;
    896:   snes->ops->view           = SNESView_NASM;
    897:   snes->ops->solve          = SNESSolve_NASM;
    898:   snes->ops->reset          = SNESReset_NASM;

    900:   snes->usesksp = PETSC_FALSE;
    901:   snes->usesnpc = PETSC_FALSE;

    903:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

    905:   nasm->fjtype              = 0;
    906:   nasm->xinit               = NULL;
    907:   nasm->eventrestrictinterp = 0;
    908:   nasm->eventsubsolve       = 0;

    910:   if (!snes->tolerancesset) {
    911:     snes->max_its   = 10000;
    912:     snes->max_funcs = 10000;
    913:   }

    915:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetType_C",SNESNASMSetType_NASM);
    916:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetType_C",SNESNASMGetType_NASM);
    917:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetSubdomains_C",SNESNASMSetSubdomains_NASM);
    918:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomains_C",SNESNASMGetSubdomains_NASM);
    919:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetDamping_C",SNESNASMSetDamping_NASM);
    920:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetDamping_C",SNESNASMGetDamping_NASM);
    921:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMGetSubdomainVecs_C",SNESNASMGetSubdomainVecs_NASM);
    922:   PetscObjectComposeFunction((PetscObject)snes,"SNESNASMSetComputeFinalJacobian_C",SNESNASMSetComputeFinalJacobian_NASM);
    923:   return(0);
    924: }

    926: /*@
    927:    SNESNASMGetSNES - Gets a subsolver

    929:    Not collective

    931:    Input Parameters:
    932: +  snes - the SNES context
    933: -  i - the number of the subsnes to get

    935:    Output Parameters:
    936: .  subsnes - the subsolver context

    938:    Level: intermediate

    940: .seealso: SNESNASM, SNESNASMGetNumber()
    941: @*/
    942: PetscErrorCode SNESNASMGetSNES(SNES snes,PetscInt i,SNES *subsnes)
    943: {
    944:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    947:   if (i < 0 || i >= nasm->n) SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_OUTOFRANGE,"No such subsolver");
    948:   *subsnes = nasm->subsnes[i];
    949:   return(0);
    950: }

    952: /*@
    953:    SNESNASMGetNumber - Gets number of subsolvers

    955:    Not collective

    957:    Input Parameters:
    958: .  snes - the SNES context

    960:    Output Parameters:
    961: .  n - the number of subsolvers

    963:    Level: intermediate

    965: .seealso: SNESNASM, SNESNASMGetSNES()
    966: @*/
    967: PetscErrorCode SNESNASMGetNumber(SNES snes,PetscInt *n)
    968: {
    969:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    972:   *n = nasm->n;
    973:   return(0);
    974: }

    976: /*@
    977:    SNESNASMSetWeight - Sets weight to use when adding overlapping updates

    979:    Collective

    981:    Input Parameters:
    982: +  snes - the SNES context
    983: -  weight - the weights to use (typically 1/N for each dof, where N is the number of patches it appears in)

    985:    Level: intermediate

    987: .seealso: SNESNASM
    988: @*/
    989: PetscErrorCode SNESNASMSetWeight(SNES snes,Vec weight)
    990: {
    991:   SNES_NASM      *nasm = (SNES_NASM*)snes->data;

    996:   VecDestroy(&nasm->weight);
    997:   nasm->weight_set = PETSC_TRUE;
    998:   nasm->weight     = weight;
    999:   PetscObjectReference((PetscObject)nasm->weight);

    1001:   return(0);
    1002: }