Actual source code: pcis.c

      1: #define PETSCKSP_DLL

     3:  #include src/ksp/pc/impls/is/pcis.h

      5: /* -------------------------------------------------------------------------- */
      6: /*
      7:    PCISSetUp -
      8: */
     11: PetscErrorCode PETSCKSP_DLLEXPORT PCISSetUp(PC pc)
     12: {
     13:   PC_IS           *pcis = (PC_IS*)(pc->data);
     14:   Mat_IS          *matis = (Mat_IS*)pc->mat->data;
     15:   PetscInt        i;
     16:   PetscErrorCode  ierr;
     17:   PetscTruth      flg;
     18:
     20:   PetscTypeCompare((PetscObject)pc->mat,MATIS,&flg);
     21:   if (!flg){
     22:     SETERRQ(PETSC_ERR_ARG_WRONG,"Preconditioner type of Neumann Neumman requires matrix of type MATIS");
     23:   }

     25:   pcis->pure_neumann = matis->pure_neumann;

     27:   /*
     28:     Creating the local vector vec1_N, containing the inverse of the number
     29:     of subdomains to which each local node (either owned or ghost)
     30:     pertains. To accomplish that, we scatter local vectors of 1's to
     31:     a global vector (adding the values); scatter the result back to
     32:     local vectors and finally invert the result.
     33:   */
     34:   {
     35:     Vec    counter;
     36:     PetscScalar one=1.0, zero=0.0;
     37:     VecDuplicate(matis->x,&pcis->vec1_N);
     38:     MatGetVecs(pc->pmat,&counter,0); /* temporary auxiliar vector */
     39:     VecSet(counter,zero);
     40:     VecSet(pcis->vec1_N,one);
     41:     VecScatterBegin(pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE,matis->ctx);
     42:     VecScatterEnd  (pcis->vec1_N,counter,ADD_VALUES,SCATTER_REVERSE,matis->ctx);
     43:     VecScatterBegin(counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD,matis->ctx);
     44:     VecScatterEnd  (counter,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD,matis->ctx);
     45:     VecDestroy(counter);
     46:   }
     47:   /*
     48:     Creating local and global index sets for interior and
     49:     inteface nodes. Notice that interior nodes have D[i]==1.0.
     50:   */
     51:   {
     52:     PetscInt     n_I;
     53:     PetscInt    *idx_I_local,*idx_B_local,*idx_I_global,*idx_B_global;
     54:     PetscScalar *array;
     55:     /* Identifying interior and interface nodes, in local numbering */
     56:     VecGetSize(pcis->vec1_N,&pcis->n);
     57:     VecGetArray(pcis->vec1_N,&array);
     58:     PetscMalloc(pcis->n*sizeof(PetscInt),&idx_I_local);
     59:     PetscMalloc(pcis->n*sizeof(PetscInt),&idx_B_local);
     60:     for (i=0, pcis->n_B=0, n_I=0; i<pcis->n; i++) {
     61:       if (array[i] == 1.0) { idx_I_local[n_I]       = i; n_I++;       }
     62:       else                 { idx_B_local[pcis->n_B] = i; pcis->n_B++; }
     63:     }
     64:     /* Getting the global numbering */
     65:     idx_B_global = idx_I_local + n_I; /* Just avoiding allocating extra memory, since we have vacant space */
     66:     idx_I_global = idx_B_local + pcis->n_B;
     67:     ISLocalToGlobalMappingApply(matis->mapping,pcis->n_B,idx_B_local,idx_B_global);
     68:     ISLocalToGlobalMappingApply(matis->mapping,n_I,      idx_I_local,idx_I_global);
     69:     /* Creating the index sets. */
     70:     ISCreateGeneral(MPI_COMM_SELF,pcis->n_B,idx_B_local, &pcis->is_B_local);
     71:     ISCreateGeneral(MPI_COMM_SELF,pcis->n_B,idx_B_global,&pcis->is_B_global);
     72:     ISCreateGeneral(MPI_COMM_SELF,n_I      ,idx_I_local, &pcis->is_I_local);
     73:     ISCreateGeneral(MPI_COMM_SELF,n_I      ,idx_I_global,&pcis->is_I_global);
     74:     /* Freeing memory and restoring arrays */
     75:     PetscFree(idx_B_local);
     76:     PetscFree(idx_I_local);
     77:     VecRestoreArray(pcis->vec1_N,&array);
     78:   }

     80:   /*
     81:     Extracting the blocks A_II, A_BI, A_IB and A_BB from A. If the numbering
     82:     is such that interior nodes come first than the interface ones, we have

     84:     [           |      ]
     85:     [    A_II   | A_IB ]
     86:     A = [           |      ]
     87:     [-----------+------]
     88:     [    A_BI   | A_BB ]
     89:   */

     91:   MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,PETSC_DECIDE,MAT_INITIAL_MATRIX,&pcis->A_II);
     92:   MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_B_local,PETSC_DECIDE,MAT_INITIAL_MATRIX,&pcis->A_IB);
     93:   MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_I_local,PETSC_DECIDE,MAT_INITIAL_MATRIX,&pcis->A_BI);
     94:   MatGetSubMatrix(matis->A,pcis->is_B_local,pcis->is_B_local,PETSC_DECIDE,MAT_INITIAL_MATRIX,&pcis->A_BB);

     96:   /*
     97:     Creating work vectors and arrays
     98:   */
     99:   /* pcis->vec1_N has already been created */
    100:   VecDuplicate(pcis->vec1_N,&pcis->vec2_N);
    101:   VecCreateSeq(PETSC_COMM_SELF,pcis->n-pcis->n_B,&pcis->vec1_D);
    102:   VecDuplicate(pcis->vec1_D,&pcis->vec2_D);
    103:   VecDuplicate(pcis->vec1_D,&pcis->vec3_D);
    104:   VecCreateSeq(PETSC_COMM_SELF,pcis->n_B,&pcis->vec1_B);
    105:   VecDuplicate(pcis->vec1_B,&pcis->vec2_B);
    106:   VecDuplicate(pcis->vec1_B,&pcis->vec3_B);
    107:   MatGetVecs(pc->pmat,&pcis->vec1_global,0);
    108:   PetscMalloc((pcis->n)*sizeof(PetscScalar),&pcis->work_N);

    110:   /* Creating the scatter contexts */
    111:   VecScatterCreate(pcis->vec1_global,pcis->is_I_global,pcis->vec1_D,(IS)0,&pcis->global_to_D);
    112:   VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);
    113:   VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);

    115:   /* Creating scaling "matrix" D, from information in vec1_N */
    116:   VecDuplicate(pcis->vec1_B,&pcis->D);
    117:   VecScatterBegin(pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);
    118:   VecScatterEnd  (pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);
    119:   VecReciprocal(pcis->D);

    121:   /* See historical note 01, at the bottom of this file. */

    123:   /*
    124:     Creating the KSP contexts for the local Dirichlet and Neumann problems.
    125:   */
    126:   {
    127:     PC  pc_ctx;
    128:     /* Dirichlet */
    129:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);
    130:     KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);
    131:     KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");
    132:     KSPGetPC(pcis->ksp_D,&pc_ctx);
    133:     PCSetType(pc_ctx,PCLU);
    134:     KSPSetType(pcis->ksp_D,KSPPREONLY);
    135:     KSPSetFromOptions(pcis->ksp_D);
    136:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
    137:     KSPSetUp(pcis->ksp_D);
    138:     /* Neumann */
    139:     KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);
    140:     KSPSetOperators(pcis->ksp_N,matis->A,matis->A,SAME_PRECONDITIONER);
    141:     KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");
    142:     KSPGetPC(pcis->ksp_N,&pc_ctx);
    143:     PCSetType(pc_ctx,PCLU);
    144:     KSPSetType(pcis->ksp_N,KSPPREONLY);
    145:     KSPSetFromOptions(pcis->ksp_N);
    146:     {
    147:       PetscTruth damp_fixed,
    148:                  remove_nullspace_fixed,
    149:                  set_damping_factor_floating,
    150:                  not_damp_floating,
    151:                  not_remove_nullspace_floating;
    152:       PetscReal  fixed_factor,
    153:                  floating_factor;

    155:       PetscOptionsGetReal(pc_ctx->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);
    156:       if (!damp_fixed) { fixed_factor = 0.0; }
    157:       PetscOptionsHasName(pc_ctx->prefix,"-pc_is_damp_fixed",&damp_fixed);

    159:       PetscOptionsHasName(pc_ctx->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed);

    161:       PetscOptionsGetReal(pc_ctx->prefix,"-pc_is_set_damping_factor_floating",
    162:                               &floating_factor,&set_damping_factor_floating);
    163:       if (!set_damping_factor_floating) { floating_factor = 0.0; }
    164:       PetscOptionsHasName(pc_ctx->prefix,"-pc_is_set_damping_factor_floating",&set_damping_factor_floating);
    165:       if (!set_damping_factor_floating) { floating_factor = 1.e-12; }

    167:       PetscOptionsHasName(pc_ctx->prefix,"-pc_is_not_damp_floating",&not_damp_floating);

    169:       PetscOptionsHasName(pc_ctx->prefix,"-pc_is_not_remove_nullspace_floating",&not_remove_nullspace_floating);

    171:       if (pcis->pure_neumann) {  /* floating subdomain */
    172:         if (!(not_damp_floating)) {
    173:           PCFactorSetShiftNonzero(pc_ctx,floating_factor);
    174:         }
    175:         if (!(not_remove_nullspace_floating)){
    176:           MatNullSpace nullsp;
    177:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,PETSC_NULL,&nullsp);
    178:           KSPSetNullSpace(pcis->ksp_N,nullsp);
    179:           MatNullSpaceDestroy(nullsp);
    180:         }
    181:       } else {  /* fixed subdomain */
    182:         if (damp_fixed) {
    183:           PCFactorSetShiftNonzero(pc_ctx,fixed_factor);
    184:         }
    185:         if (remove_nullspace_fixed) {
    186:           MatNullSpace nullsp;
    187:           MatNullSpaceCreate(PETSC_COMM_SELF,PETSC_TRUE,0,PETSC_NULL,&nullsp);
    188:           KSPSetNullSpace(pcis->ksp_N,nullsp);
    189:           MatNullSpaceDestroy(nullsp);
    190:         }
    191:       }
    192:     }
    193:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
    194:     KSPSetUp(pcis->ksp_N);
    195:   }

    197:   ISLocalToGlobalMappingGetInfo(((Mat_IS*)(pc->mat->data))->mapping,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
    198:   pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_TRUE;
    199:   return(0);
    200: }

    202: /* -------------------------------------------------------------------------- */
    203: /*
    204:    PCISDestroy -
    205: */
    208: PetscErrorCode PETSCKSP_DLLEXPORT PCISDestroy(PC pc)
    209: {
    210:   PC_IS          *pcis = (PC_IS*)(pc->data);

    214:   if (pcis->is_B_local)  {ISDestroy(pcis->is_B_local);}
    215:   if (pcis->is_I_local)  {ISDestroy(pcis->is_I_local);}
    216:   if (pcis->is_B_global) {ISDestroy(pcis->is_B_global);}
    217:   if (pcis->is_I_global) {ISDestroy(pcis->is_I_global);}
    218:   if (pcis->A_II)        {MatDestroy(pcis->A_II);}
    219:   if (pcis->A_IB)        {MatDestroy(pcis->A_IB);}
    220:   if (pcis->A_BI)        {MatDestroy(pcis->A_BI);}
    221:   if (pcis->A_BB)        {MatDestroy(pcis->A_BB);}
    222:   if (pcis->D)           {VecDestroy(pcis->D);}
    223:   if (pcis->ksp_N)      {KSPDestroy(pcis->ksp_N);}
    224:   if (pcis->ksp_D)      {KSPDestroy(pcis->ksp_D);}
    225:   if (pcis->vec1_N)      {VecDestroy(pcis->vec1_N);}
    226:   if (pcis->vec2_N)      {VecDestroy(pcis->vec2_N);}
    227:   if (pcis->vec1_D)      {VecDestroy(pcis->vec1_D);}
    228:   if (pcis->vec2_D)      {VecDestroy(pcis->vec2_D);}
    229:   if (pcis->vec3_D)      {VecDestroy(pcis->vec3_D);}
    230:   if (pcis->vec1_B)      {VecDestroy(pcis->vec1_B);}
    231:   if (pcis->vec2_B)      {VecDestroy(pcis->vec2_B);}
    232:   if (pcis->vec3_B)      {VecDestroy(pcis->vec3_B);}
    233:   if (pcis->vec1_global) {VecDestroy(pcis->vec1_global);}
    234:   if (pcis->work_N)      {PetscFree(pcis->work_N);}
    235:   if (pcis->global_to_D) {VecScatterDestroy(pcis->global_to_D);}
    236:   if (pcis->N_to_B)      {VecScatterDestroy(pcis->N_to_B);}
    237:   if (pcis->global_to_B) {VecScatterDestroy(pcis->global_to_B);}
    238:   if (pcis->ISLocalToGlobalMappingGetInfoWasCalled) {
    239:     ISLocalToGlobalMappingRestoreInfo((ISLocalToGlobalMapping)0,&(pcis->n_neigh),&(pcis->neigh),&(pcis->n_shared),&(pcis->shared));
    240:   }
    241:   return(0);
    242: }

    244: /* -------------------------------------------------------------------------- */
    245: /*
    246:    PCISCreate -
    247: */
    250: PetscErrorCode PETSCKSP_DLLEXPORT PCISCreate(PC pc)
    251: {
    252:   PC_IS *pcis = (PC_IS*)(pc->data);

    255:   pcis->is_B_local  = 0;
    256:   pcis->is_I_local  = 0;
    257:   pcis->is_B_global = 0;
    258:   pcis->is_I_global = 0;
    259:   pcis->A_II        = 0;
    260:   pcis->A_IB        = 0;
    261:   pcis->A_BI        = 0;
    262:   pcis->A_BB        = 0;
    263:   pcis->D           = 0;
    264:   pcis->ksp_N      = 0;
    265:   pcis->ksp_D      = 0;
    266:   pcis->vec1_N      = 0;
    267:   pcis->vec2_N      = 0;
    268:   pcis->vec1_D      = 0;
    269:   pcis->vec2_D      = 0;
    270:   pcis->vec3_D      = 0;
    271:   pcis->vec1_B      = 0;
    272:   pcis->vec2_B      = 0;
    273:   pcis->vec3_B      = 0;
    274:   pcis->vec1_global = 0;
    275:   pcis->work_N      = 0;
    276:   pcis->global_to_D = 0;
    277:   pcis->N_to_B      = 0;
    278:   pcis->global_to_B = 0;
    279:   pcis->ISLocalToGlobalMappingGetInfoWasCalled = PETSC_FALSE;
    280:   return(0);
    281: }

    283: /* -------------------------------------------------------------------------- */
    284: /*
    285:    PCISApplySchur -

    287:    Input parameters:
    288: .  pc - preconditioner context
    289: .  v - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)

    291:    Output parameters:
    292: .  vec1_B - result of Schur complement applied to chunk
    293: .  vec2_B - garbage (used as work space), or null (and v is used as workspace)
    294: .  vec1_D - garbage (used as work space)
    295: .  vec2_D - garbage (used as work space)

    297: */
    300: PetscErrorCode PETSCKSP_DLLEXPORT PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
    301: {
    303:   PetscScalar    m_one = -1.0;
    304:   PC_IS          *pcis = (PC_IS*)(pc->data);

    307:   if (!vec2_B) { vec2_B = v; }

    309:   MatMult(pcis->A_BB,v,vec1_B);
    310:   MatMult(pcis->A_IB,v,vec1_D);
    311:   KSPSolve(pcis->ksp_D,vec1_D,vec2_D);
    312:   MatMult(pcis->A_BI,vec2_D,vec2_B);
    313:   VecAXPY(vec1_B,m_one,vec2_B);
    314:   return(0);
    315: }

    317: /* -------------------------------------------------------------------------- */
    318: /*
    319:    PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
    320:    including ghosts) into an interface vector, when in SCATTER_FORWARD mode, or vice-versa, when in SCATTER_REVERSE
    321:    mode.

    323:    Input parameters:
    324: .  pc - preconditioner context
    325: .  array_N - [when in SCATTER_FORWARD mode] Array to be scattered into the vector
    326: .  v_B - [when in SCATTER_REVERSE mode] Vector to be scattered into the array

    328:    Output parameter:
    329: .  array_N - [when in SCATTER_REVERSE mode] Array to receive the scattered vector
    330: .  v_B - [when in SCATTER_FORWARD mode] Vector to receive the scattered array

    332:    Notes:
    333:    The entries in the array that do not correspond to interface nodes remain unaltered.
    334: */
    337: PetscErrorCode PETSCKSP_DLLEXPORT PCISScatterArrayNToVecB (PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode, PC pc)
    338: {
    339:   PetscInt       i, *idex;
    341:   PetscScalar    *array_B;
    342:   PC_IS          *pcis = (PC_IS*)(pc->data);

    345:   VecGetArray(v_B,&array_B);
    346:   ISGetIndices(pcis->is_B_local,&idex);

    348:   if (smode == SCATTER_FORWARD) {
    349:     if (imode == INSERT_VALUES) {
    350:       for (i=0; i<pcis->n_B; i++) { array_B[i]  = array_N[idex[i]]; }
    351:     } else {  /* ADD_VALUES */
    352:       for (i=0; i<pcis->n_B; i++) { array_B[i] += array_N[idex[i]]; }
    353:     }
    354:   } else {  /* SCATTER_REVERSE */
    355:     if (imode == INSERT_VALUES) {
    356:       for (i=0; i<pcis->n_B; i++) { array_N[idex[i]]  = array_B[i]; }
    357:     } else {  /* ADD_VALUES */
    358:       for (i=0; i<pcis->n_B; i++) { array_N[idex[i]] += array_B[i]; }
    359:     }
    360:   }
    361:   ISRestoreIndices(pcis->is_B_local,&idex);
    362:   VecRestoreArray(v_B,&array_B);
    363:   return(0);
    364: }

    366: /* -------------------------------------------------------------------------- */
    367: /*
    368:    PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.
    369:    More precisely, solves the problem:
    370:                                         [ A_II  A_IB ] [ . ]   [ 0 ]
    371:                                         [            ] [   ] = [   ]
    372:                                         [ A_BI  A_BB ] [ x ]   [ b ]

    374:    Input parameters:
    375: .  pc - preconditioner context
    376: .  b - vector of local interface nodes (including ghosts)

    378:    Output parameters:
    379: .  x - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur
    380:        complement to b
    381: .  vec1_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)
    382: .  vec2_N - vector of local nodes (interior and interface, including ghosts); returns garbage (used as work space)

    384: */
    387: PetscErrorCode PETSCKSP_DLLEXPORT PCISApplyInvSchur (PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
    388: {
    390:   PC_IS          *pcis = (PC_IS*)(pc->data);
    391:   PetscScalar    zero  = 0.0;

    394:   /*
    395:     Neumann solvers.
    396:     Applying the inverse of the local Schur complement, i.e, solving a Neumann
    397:     Problem with zero at the interior nodes of the RHS and extracting the interface
    398:     part of the solution. inverse Schur complement is applied to b and the result
    399:     is stored in x.
    400:   */
    401:   /* Setting the RHS vec1_N */
    402:   VecSet(vec1_N,zero);
    403:   VecScatterBegin(b,vec1_N,INSERT_VALUES,SCATTER_REVERSE,pcis->N_to_B);
    404:   VecScatterEnd  (b,vec1_N,INSERT_VALUES,SCATTER_REVERSE,pcis->N_to_B);
    405:   /* Checking for consistency of the RHS */
    406:   {
    407:     PetscTruth flg;
    408:     PetscOptionsHasName(PETSC_NULL,"-pc_is_check_consistency",&flg);
    409:     if (flg) {
    410:       PetscScalar average;
    411:       VecSum(vec1_N,&average);
    412:       average = average / ((PetscReal)pcis->n);
    413:       if (pcis->pure_neumann) {
    414:         PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_(pc->comm),"Subdomain %04d is floating. Average = % 1.14e\n",
    415:                                              PetscGlobalRank,PetscAbsScalar(average));
    416:       } else {
    417:         PetscViewerASCIISynchronizedPrintf(PETSC_VIEWER_STDOUT_(pc->comm),"Subdomain %04d is fixed.    Average = % 1.14e\n",
    418:                                              PetscGlobalRank,PetscAbsScalar(average));
    419:       }
    420:       PetscViewerFlush(PETSC_VIEWER_STDOUT_(pc->comm));
    421:     }
    422:   }
    423:   /* Solving the system for vec2_N */
    424:   KSPSolve(pcis->ksp_N,vec1_N,vec2_N);
    425:   /* Extracting the local interface vector out of the solution */
    426:   VecScatterBegin(vec2_N,x,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);
    427:   VecScatterEnd  (vec2_N,x,INSERT_VALUES,SCATTER_FORWARD,pcis->N_to_B);
    428:   return(0);
    429: }