Actual source code: pcis.c

      1: #include <petsc/private/pcisimpl.h>

      3: static PetscErrorCode PCISSetUseStiffnessScaling_IS(PC pc, PetscBool use)
      4: {
      5:   PC_IS *pcis = (PC_IS *)pc->data;

      7:   PetscFunctionBegin;
      8:   pcis->use_stiffness_scaling = use;
      9:   PetscFunctionReturn(PETSC_SUCCESS);
     10: }

     12: /*@
     13:   PCISSetUseStiffnessScaling - Tells `PCIS` to construct partition of unity using
     14:   the local matrices' diagonal entries

     16:   Logically Collective

     18:   Input Parameters:
     19: + pc  - the preconditioning context
     20: - use - whether or not it should use matrix diagonal to build partition of unity.

     22:   Level: intermediate

     24: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
     25:           `PCISSetSubdomainScalingFactor()`,
     26:           `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
     27: @*/
     28: PetscErrorCode PCISSetUseStiffnessScaling(PC pc, PetscBool use)
     29: {
     30:   PetscFunctionBegin;
     33:   PetscTryMethod(pc, "PCISSetUseStiffnessScaling_C", (PC, PetscBool), (pc, use));
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: static PetscErrorCode PCISSetSubdomainDiagonalScaling_IS(PC pc, Vec scaling_factors)
     38: {
     39:   PC_IS *pcis = (PC_IS *)pc->data;

     41:   PetscFunctionBegin;
     42:   PetscCall(PetscObjectReference((PetscObject)scaling_factors));
     43:   PetscCall(VecDestroy(&pcis->D));
     44:   pcis->D = scaling_factors;
     45:   if (pc->setupcalled) {
     46:     PetscInt sn;

     48:     PetscCall(VecGetSize(pcis->D, &sn));
     49:     if (sn == pcis->n) {
     50:       PetscCall(VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
     51:       PetscCall(VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
     52:       PetscCall(VecDestroy(&pcis->D));
     53:       PetscCall(VecDuplicate(pcis->vec1_B, &pcis->D));
     54:       PetscCall(VecCopy(pcis->vec1_B, pcis->D));
     55:     } else PetscCheck(sn == pcis->n_B, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Invalid size for scaling vector. Expected %" PetscInt_FMT " (or full %" PetscInt_FMT "), found %" PetscInt_FMT, pcis->n_B, pcis->n, sn);
     56:   }
     57:   PetscFunctionReturn(PETSC_SUCCESS);
     58: }

     60: /*@
     61:   PCISSetSubdomainDiagonalScaling - Set diagonal scaling for `PCIS`.

     63:   Logically Collective

     65:   Input Parameters:
     66: + pc              - the preconditioning context
     67: - scaling_factors - scaling factors for the subdomain

     69:   Level: intermediate

     71:   Note:
     72:   Intended for use with jumping coefficients cases.

     74: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISScatterArrayNToVecB()`,
     75:           `PCISSetSubdomainScalingFactor()`, `PCISSetUseStiffnessScaling()`,
     76:           `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
     77: @*/
     78: PetscErrorCode PCISSetSubdomainDiagonalScaling(PC pc, Vec scaling_factors)
     79: {
     80:   PetscFunctionBegin;
     83:   PetscTryMethod(pc, "PCISSetSubdomainDiagonalScaling_C", (PC, Vec), (pc, scaling_factors));
     84:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

     87: static PetscErrorCode PCISSetSubdomainScalingFactor_IS(PC pc, PetscScalar scal)
     88: {
     89:   PC_IS *pcis = (PC_IS *)pc->data;

     91:   PetscFunctionBegin;
     92:   pcis->scaling_factor = scal;
     93:   if (pcis->D) PetscCall(VecSet(pcis->D, pcis->scaling_factor));
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: /*@
     98:   PCISSetSubdomainScalingFactor - Set scaling factor for `PCIS`.

    100:   Not Collective

    102:   Input Parameters:
    103: + pc   - the preconditioning context
    104: - scal - scaling factor for the subdomain

    106:   Level: intermediate

    108:   Note:
    109:   Intended for use with the jumping coefficients cases.

    111: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISScatterArrayNToVecB()`,
    112:           `PCISSetSubdomainDiagonalScaling()`, `PCISSetUseStiffnessScaling()`,
    113:           `PCISReset()`, `PCISInitialize()`, `PCISApplyInvSchur()`, `PCISApplySchur()`
    114: @*/
    115: PetscErrorCode PCISSetSubdomainScalingFactor(PC pc, PetscScalar scal)
    116: {
    117:   PetscFunctionBegin;
    119:   PetscTryMethod(pc, "PCISSetSubdomainScalingFactor_C", (PC, PetscScalar), (pc, scal));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: /*@
    124:   PCISSetUp - sets up the `PC_IS` portion of `PCNN` and `PCBDDC` preconditioner context as part of their setup process

    126:   Input Parameters:
    127: + pc              - the `PC` object, must be of type `PCNN` or `PCBDDC`
    128: . computematrices - Extract the blocks `A_II`, `A_BI`, `A_IB` and `A_BB` from the matrix
    129: - computesolvers  - Create the `KSP` for the local Dirichlet and Neumann problems

    131:   Level: advanced

    133: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
    134:           `PCISSetSubdomainScalingFactor()`,
    135:           `PCISReset()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
    136: @*/
    137: PetscErrorCode PCISSetUp(PC pc, PetscBool computematrices, PetscBool computesolvers)
    138: {
    139:   PC_IS    *pcis = (PC_IS *)pc->data;
    140:   Mat_IS   *matis;
    141:   MatReuse  reuse;
    142:   PetscBool flg, issbaij;

    144:   PetscFunctionBegin;
    145:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &flg));
    146:   PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Requires preconditioning matrix of type MATIS");
    147:   matis = (Mat_IS *)pc->pmat->data;
    148:   if (pc->useAmat) {
    149:     PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATIS, &flg));
    150:     PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Requires linear system matrix of type MATIS");
    151:   }

    153:   /* first time creation, get info on substructuring */
    154:   if (!pc->setupcalled) {
    155:     PetscInt  n_I;
    156:     PetscInt *idx_I_local, *idx_B_local, *idx_I_global, *idx_B_global, *count;
    157:     PetscInt  i;

    159:     /* get info on mapping */
    160:     PetscCall(PetscObjectReference((PetscObject)matis->rmapping));
    161:     PetscCall(ISLocalToGlobalMappingDestroy(&pcis->mapping));
    162:     pcis->mapping = matis->rmapping;
    163:     PetscCall(ISLocalToGlobalMappingGetSize(pcis->mapping, &pcis->n));
    164:     PetscCall(ISLocalToGlobalMappingGetInfo(pcis->mapping, &pcis->n_neigh, &pcis->neigh, &pcis->n_shared, &pcis->shared));

    166:     /* Identifying interior and interface nodes, in local numbering */
    167:     PetscCall(ISLocalToGlobalMappingGetNodeInfo(pcis->mapping, NULL, &count, NULL));
    168:     PetscCall(PetscMalloc1(pcis->n, &idx_I_local));
    169:     PetscCall(PetscMalloc1(pcis->n, &idx_B_local));
    170:     for (i = 0, pcis->n_B = 0, n_I = 0; i < pcis->n; i++) {
    171:       if (count[i] < 2) {
    172:         idx_I_local[n_I] = i;
    173:         n_I++;
    174:       } else {
    175:         idx_B_local[pcis->n_B] = i;
    176:         pcis->n_B++;
    177:       }
    178:     }
    179:     PetscCall(ISLocalToGlobalMappingRestoreNodeInfo(pcis->mapping, NULL, &count, NULL));

    181:     /* Getting the global numbering */
    182:     idx_B_global = PetscSafePointerPlusOffset(idx_I_local, n_I); /* Just avoiding allocating extra memory, since we have vacant space */
    183:     idx_I_global = PetscSafePointerPlusOffset(idx_B_local, pcis->n_B);
    184:     PetscCall(ISLocalToGlobalMappingApply(pcis->mapping, pcis->n_B, idx_B_local, idx_B_global));
    185:     PetscCall(ISLocalToGlobalMappingApply(pcis->mapping, n_I, idx_I_local, idx_I_global));

    187:     /* Creating the index sets */
    188:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pcis->n_B, idx_B_local, PETSC_COPY_VALUES, &pcis->is_B_local));
    189:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), pcis->n_B, idx_B_global, PETSC_COPY_VALUES, &pcis->is_B_global));
    190:     PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n_I, idx_I_local, PETSC_COPY_VALUES, &pcis->is_I_local));
    191:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), n_I, idx_I_global, PETSC_COPY_VALUES, &pcis->is_I_global));

    193:     /* Freeing memory */
    194:     PetscCall(PetscFree(idx_B_local));
    195:     PetscCall(PetscFree(idx_I_local));

    197:     /* Creating work vectors and arrays */
    198:     PetscCall(VecDuplicate(matis->x, &pcis->vec1_N));
    199:     PetscCall(VecDuplicate(pcis->vec1_N, &pcis->vec2_N));
    200:     PetscCall(VecCreate(PETSC_COMM_SELF, &pcis->vec1_D));
    201:     PetscCall(VecSetSizes(pcis->vec1_D, pcis->n - pcis->n_B, PETSC_DECIDE));
    202:     PetscCall(VecSetType(pcis->vec1_D, ((PetscObject)pcis->vec1_N)->type_name));
    203:     PetscCall(VecDuplicate(pcis->vec1_D, &pcis->vec2_D));
    204:     PetscCall(VecDuplicate(pcis->vec1_D, &pcis->vec3_D));
    205:     PetscCall(VecDuplicate(pcis->vec1_D, &pcis->vec4_D));
    206:     PetscCall(VecCreate(PETSC_COMM_SELF, &pcis->vec1_B));
    207:     PetscCall(VecSetSizes(pcis->vec1_B, pcis->n_B, PETSC_DECIDE));
    208:     PetscCall(VecSetType(pcis->vec1_B, ((PetscObject)pcis->vec1_N)->type_name));
    209:     PetscCall(VecDuplicate(pcis->vec1_B, &pcis->vec2_B));
    210:     PetscCall(VecDuplicate(pcis->vec1_B, &pcis->vec3_B));
    211:     PetscCall(MatCreateVecs(pc->pmat, &pcis->vec1_global, NULL));
    212:     PetscCall(PetscMalloc1(pcis->n, &pcis->work_N));
    213:     /* scaling vector */
    214:     if (!pcis->D) { /* it can happen that the user passed in a scaling vector via PCISSetSubdomainDiagonalScaling */
    215:       PetscCall(VecDuplicate(pcis->vec1_B, &pcis->D));
    216:       PetscCall(VecSet(pcis->D, pcis->scaling_factor));
    217:     }

    219:     /* Creating the scatter contexts */
    220:     PetscCall(VecScatterCreate(pcis->vec1_N, pcis->is_I_local, pcis->vec1_D, (IS)0, &pcis->N_to_D));
    221:     PetscCall(VecScatterCreate(pcis->vec1_global, pcis->is_I_global, pcis->vec1_D, (IS)0, &pcis->global_to_D));
    222:     PetscCall(VecScatterCreate(pcis->vec1_N, pcis->is_B_local, pcis->vec1_B, (IS)0, &pcis->N_to_B));
    223:     PetscCall(VecScatterCreate(pcis->vec1_global, pcis->is_B_global, pcis->vec1_B, (IS)0, &pcis->global_to_B));

    225:     /* map from boundary to local */
    226:     PetscCall(ISLocalToGlobalMappingCreateIS(pcis->is_B_local, &pcis->BtoNmap));
    227:   }

    229:   {
    230:     PetscInt sn;

    232:     PetscCall(VecGetSize(pcis->D, &sn));
    233:     if (sn == pcis->n) {
    234:       PetscCall(VecScatterBegin(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    235:       PetscCall(VecScatterEnd(pcis->N_to_B, pcis->D, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    236:       PetscCall(VecDestroy(&pcis->D));
    237:       PetscCall(VecDuplicate(pcis->vec1_B, &pcis->D));
    238:       PetscCall(VecCopy(pcis->vec1_B, pcis->D));
    239:     } else PetscCheck(sn == pcis->n_B, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Invalid size for scaling vector. Expected %" PetscInt_FMT " (or full %" PetscInt_FMT "), found %" PetscInt_FMT, pcis->n_B, pcis->n, sn);
    240:   }

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

    246:         [ A_II | A_IB ]
    247:     A = [------+------]
    248:         [ A_BI | A_BB ]
    249:   */
    250:   if (computematrices) {
    251:     PetscBool amat = (PetscBool)(pc->mat != pc->pmat && pc->useAmat);
    252:     PetscInt  bs, ibs;

    254:     reuse = MAT_INITIAL_MATRIX;
    255:     if (pcis->reusesubmatrices && pc->setupcalled) {
    256:       if (pc->flag == SAME_NONZERO_PATTERN) {
    257:         reuse = MAT_REUSE_MATRIX;
    258:       } else {
    259:         reuse = MAT_INITIAL_MATRIX;
    260:       }
    261:     }
    262:     if (reuse == MAT_INITIAL_MATRIX) {
    263:       PetscCall(MatDestroy(&pcis->A_II));
    264:       PetscCall(MatDestroy(&pcis->pA_II));
    265:       PetscCall(MatDestroy(&pcis->A_IB));
    266:       PetscCall(MatDestroy(&pcis->A_BI));
    267:       PetscCall(MatDestroy(&pcis->A_BB));
    268:     }

    270:     PetscCall(ISLocalToGlobalMappingGetBlockSize(pcis->mapping, &ibs));
    271:     PetscCall(MatGetBlockSize(matis->A, &bs));
    272:     PetscCall(MatCreateSubMatrix(matis->A, pcis->is_I_local, pcis->is_I_local, reuse, &pcis->pA_II));
    273:     if (amat) {
    274:       Mat_IS *amatis = (Mat_IS *)pc->mat->data;
    275:       PetscCall(MatCreateSubMatrix(amatis->A, pcis->is_I_local, pcis->is_I_local, reuse, &pcis->A_II));
    276:     } else {
    277:       PetscCall(PetscObjectReference((PetscObject)pcis->pA_II));
    278:       PetscCall(MatDestroy(&pcis->A_II));
    279:       pcis->A_II = pcis->pA_II;
    280:     }
    281:     PetscCall(MatSetBlockSize(pcis->A_II, bs == ibs ? bs : 1));
    282:     PetscCall(MatSetBlockSize(pcis->pA_II, bs == ibs ? bs : 1));
    283:     PetscCall(MatCreateSubMatrix(matis->A, pcis->is_B_local, pcis->is_B_local, reuse, &pcis->A_BB));
    284:     PetscCall(PetscObjectTypeCompare((PetscObject)matis->A, MATSEQSBAIJ, &issbaij));
    285:     if (!issbaij) {
    286:       PetscCall(MatCreateSubMatrix(matis->A, pcis->is_I_local, pcis->is_B_local, reuse, &pcis->A_IB));
    287:       PetscCall(MatCreateSubMatrix(matis->A, pcis->is_B_local, pcis->is_I_local, reuse, &pcis->A_BI));
    288:     } else {
    289:       Mat newmat;

    291:       PetscCall(MatConvert(matis->A, MATSEQBAIJ, MAT_INITIAL_MATRIX, &newmat));
    292:       PetscCall(MatCreateSubMatrix(newmat, pcis->is_I_local, pcis->is_B_local, reuse, &pcis->A_IB));
    293:       PetscCall(MatCreateSubMatrix(newmat, pcis->is_B_local, pcis->is_I_local, reuse, &pcis->A_BI));
    294:       PetscCall(MatDestroy(&newmat));
    295:     }
    296:     PetscCall(MatSetBlockSize(pcis->A_BB, bs == ibs ? bs : 1));
    297:   }

    299:   /* Creating scaling vector D */
    300:   PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_is_use_stiffness_scaling", &pcis->use_stiffness_scaling, NULL));
    301:   if (pcis->use_stiffness_scaling) {
    302:     PetscScalar *a;
    303:     PetscInt     i, n;

    305:     if (pcis->A_BB) {
    306:       PetscCall(MatGetDiagonal(pcis->A_BB, pcis->D));
    307:     } else {
    308:       PetscCall(MatGetDiagonal(matis->A, pcis->vec1_N));
    309:       PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD));
    310:       PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->D, INSERT_VALUES, SCATTER_FORWARD));
    311:     }
    312:     PetscCall(VecAbs(pcis->D));
    313:     PetscCall(VecGetLocalSize(pcis->D, &n));
    314:     PetscCall(VecGetArray(pcis->D, &a));
    315:     for (i = 0; i < n; i++)
    316:       if (PetscAbsScalar(a[i]) < PETSC_SMALL) a[i] = 1.0;
    317:     PetscCall(VecRestoreArray(pcis->D, &a));
    318:   }
    319:   PetscCall(VecSet(pcis->vec1_global, 0.0));
    320:   PetscCall(VecScatterBegin(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    321:   PetscCall(VecScatterEnd(pcis->global_to_B, pcis->D, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    322:   PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    323:   PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    324:   PetscCall(VecPointwiseDivide(pcis->D, pcis->D, pcis->vec1_B));
    325:   /* See historical note 01, at the bottom of this file. */

    327:   /* Creating the KSP contexts for the local Dirichlet and Neumann problems */
    328:   if (computesolvers) {
    329:     PC pc_ctx;

    331:     pcis->pure_neumann = matis->pure_neumann;
    332:     /* Dirichlet */
    333:     PetscCall(KSPCreate(PETSC_COMM_SELF, &pcis->ksp_D));
    334:     PetscCall(KSPSetNestLevel(pcis->ksp_D, pc->kspnestlevel));
    335:     PetscCall(KSPSetErrorIfNotConverged(pcis->ksp_D, pc->erroriffailure));
    336:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D, (PetscObject)pc, 1));
    337:     PetscCall(KSPSetOperators(pcis->ksp_D, pcis->A_II, pcis->A_II));
    338:     PetscCall(KSPSetOptionsPrefix(pcis->ksp_D, "is_localD_"));
    339:     PetscCall(KSPGetPC(pcis->ksp_D, &pc_ctx));
    340:     PetscCall(PCSetType(pc_ctx, PCLU));
    341:     PetscCall(KSPSetType(pcis->ksp_D, KSPPREONLY));
    342:     PetscCall(KSPSetFromOptions(pcis->ksp_D));
    343:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
    344:     PetscCall(KSPSetUp(pcis->ksp_D));
    345:     /* Neumann */
    346:     PetscCall(KSPCreate(PETSC_COMM_SELF, &pcis->ksp_N));
    347:     PetscCall(KSPSetNestLevel(pcis->ksp_N, pc->kspnestlevel));
    348:     PetscCall(KSPSetErrorIfNotConverged(pcis->ksp_N, pc->erroriffailure));
    349:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N, (PetscObject)pc, 1));
    350:     PetscCall(KSPSetOperators(pcis->ksp_N, matis->A, matis->A));
    351:     PetscCall(KSPSetOptionsPrefix(pcis->ksp_N, "is_localN_"));
    352:     PetscCall(KSPGetPC(pcis->ksp_N, &pc_ctx));
    353:     PetscCall(PCSetType(pc_ctx, PCLU));
    354:     PetscCall(KSPSetType(pcis->ksp_N, KSPPREONLY));
    355:     PetscCall(KSPSetFromOptions(pcis->ksp_N));
    356:     {
    357:       PetscBool damp_fixed = PETSC_FALSE, remove_nullspace_fixed = PETSC_FALSE, set_damping_factor_floating = PETSC_FALSE, not_damp_floating = PETSC_FALSE, not_remove_nullspace_floating = PETSC_FALSE;
    358:       PetscReal fixed_factor, floating_factor;

    360:       PetscCall(PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &fixed_factor, &damp_fixed));
    361:       if (!damp_fixed) fixed_factor = 0.0;
    362:       PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_damp_fixed", &damp_fixed, NULL));

    364:       PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_remove_nullspace_fixed", &remove_nullspace_fixed, NULL));

    366:       PetscCall(PetscOptionsGetReal(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &floating_factor, &set_damping_factor_floating));
    367:       if (!set_damping_factor_floating) floating_factor = 0.0;
    368:       PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_set_damping_factor_floating", &set_damping_factor_floating, NULL));
    369:       if (!set_damping_factor_floating) floating_factor = 1.e-12;

    371:       PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_damp_floating", &not_damp_floating, NULL));

    373:       PetscCall(PetscOptionsGetBool(((PetscObject)pc_ctx)->options, ((PetscObject)pc_ctx)->prefix, "-pc_is_not_remove_nullspace_floating", &not_remove_nullspace_floating, NULL));

    375:       if (pcis->pure_neumann) { /* floating subdomain */
    376:         if (!(not_damp_floating)) {
    377:           PetscCall(PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO));
    378:           PetscCall(PCFactorSetShiftAmount(pc_ctx, floating_factor));
    379:         }
    380:         if (!(not_remove_nullspace_floating)) {
    381:           MatNullSpace nullsp;
    382:           PetscCall(MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp));
    383:           PetscCall(MatSetNullSpace(matis->A, nullsp));
    384:           PetscCall(MatNullSpaceDestroy(&nullsp));
    385:         }
    386:       } else { /* fixed subdomain */
    387:         if (damp_fixed) {
    388:           PetscCall(PCFactorSetShiftType(pc_ctx, MAT_SHIFT_NONZERO));
    389:           PetscCall(PCFactorSetShiftAmount(pc_ctx, floating_factor));
    390:         }
    391:         if (remove_nullspace_fixed) {
    392:           MatNullSpace nullsp;
    393:           PetscCall(MatNullSpaceCreate(PETSC_COMM_SELF, PETSC_TRUE, 0, NULL, &nullsp));
    394:           PetscCall(MatSetNullSpace(matis->A, nullsp));
    395:           PetscCall(MatNullSpaceDestroy(&nullsp));
    396:         }
    397:       }
    398:     }
    399:     /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
    400:     PetscCall(KSPSetUp(pcis->ksp_N));
    401:   }
    402:   PetscFunctionReturn(PETSC_SUCCESS);
    403: }

    405: /*@
    406:   PCISReset - Removes all the `PC_IS` parts of the `PC` implementation data structure

    408:   Input Parameter:
    409: . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`

    411:   Level: advanced

    413: .seealso: [](ch_ksp), `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`, `PCISSetSubdomainScalingFactor()`,
    414:           `PCISInitialize()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
    415: @*/
    416: PetscErrorCode PCISReset(PC pc)
    417: {
    418:   PC_IS    *pcis = (PC_IS *)pc->data;
    419:   PetscBool correcttype;

    421:   PetscFunctionBegin;
    422:   if (!pc) PetscFunctionReturn(PETSC_SUCCESS);
    423:   PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &correcttype, PCBDDC, PCNN, ""));
    424:   PetscCheck(correcttype, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC must be of type PCNN or PCBDDC");
    425:   PetscCall(ISDestroy(&pcis->is_B_local));
    426:   PetscCall(ISDestroy(&pcis->is_I_local));
    427:   PetscCall(ISDestroy(&pcis->is_B_global));
    428:   PetscCall(ISDestroy(&pcis->is_I_global));
    429:   PetscCall(MatDestroy(&pcis->A_II));
    430:   PetscCall(MatDestroy(&pcis->pA_II));
    431:   PetscCall(MatDestroy(&pcis->A_IB));
    432:   PetscCall(MatDestroy(&pcis->A_BI));
    433:   PetscCall(MatDestroy(&pcis->A_BB));
    434:   PetscCall(VecDestroy(&pcis->D));
    435:   PetscCall(KSPDestroy(&pcis->ksp_N));
    436:   PetscCall(KSPDestroy(&pcis->ksp_D));
    437:   PetscCall(VecDestroy(&pcis->vec1_N));
    438:   PetscCall(VecDestroy(&pcis->vec2_N));
    439:   PetscCall(VecDestroy(&pcis->vec1_D));
    440:   PetscCall(VecDestroy(&pcis->vec2_D));
    441:   PetscCall(VecDestroy(&pcis->vec3_D));
    442:   PetscCall(VecDestroy(&pcis->vec4_D));
    443:   PetscCall(VecDestroy(&pcis->vec1_B));
    444:   PetscCall(VecDestroy(&pcis->vec2_B));
    445:   PetscCall(VecDestroy(&pcis->vec3_B));
    446:   PetscCall(VecDestroy(&pcis->vec1_global));
    447:   PetscCall(VecScatterDestroy(&pcis->global_to_D));
    448:   PetscCall(VecScatterDestroy(&pcis->N_to_B));
    449:   PetscCall(VecScatterDestroy(&pcis->N_to_D));
    450:   PetscCall(VecScatterDestroy(&pcis->global_to_B));
    451:   PetscCall(PetscFree(pcis->work_N));
    452:   if (pcis->n_neigh > -1) PetscCall(ISLocalToGlobalMappingRestoreInfo(pcis->mapping, &pcis->n_neigh, &pcis->neigh, &pcis->n_shared, &pcis->shared));
    453:   PetscCall(ISLocalToGlobalMappingDestroy(&pcis->mapping));
    454:   PetscCall(ISLocalToGlobalMappingDestroy(&pcis->BtoNmap));
    455:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", NULL));
    456:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", NULL));
    457:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", NULL));
    458:   PetscFunctionReturn(PETSC_SUCCESS);
    459: }

    461: /*@
    462:   PCISInitialize - initializes the `PC_IS` portion of `PCNN` and `PCBDDC` preconditioner context

    464:   Input Parameter:
    465: . pc - the `PC` object, must be of type `PCNN` or `PCBDDC`

    467:   Level: advanced

    469:   Note:
    470:   There is no preconditioner the `PCIS` prefixed routines provide functionality needed by `PCNN` or `PCBDDC`

    472: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
    473:           `PCISSetSubdomainScalingFactor()`,
    474:           `PCISReset()`, `PCISApplySchur()`, `PCISApplyInvSchur()`
    475: @*/
    476: PetscErrorCode PCISInitialize(PC pc)
    477: {
    478:   PC_IS    *pcis = (PC_IS *)pc->data;
    479:   PetscBool correcttype;

    481:   PetscFunctionBegin;
    482:   PetscCheck(pcis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC_IS context must be created by caller");
    483:   PetscCall(PetscObjectTypeCompareAny((PetscObject)pc, &correcttype, PCBDDC, PCNN, ""));
    484:   PetscCheck(correcttype, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PC must be of type PCNN or PCBDDC");
    485:   pcis->n_neigh          = -1;
    486:   pcis->scaling_factor   = 1.0;
    487:   pcis->reusesubmatrices = PETSC_TRUE;
    488:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetUseStiffnessScaling_C", PCISSetUseStiffnessScaling_IS));
    489:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainScalingFactor_C", PCISSetSubdomainScalingFactor_IS));
    490:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCISSetSubdomainDiagonalScaling_C", PCISSetSubdomainDiagonalScaling_IS));
    491:   PetscFunctionReturn(PETSC_SUCCESS);
    492: }

    494: /*@
    495:   PCISApplySchur - applies the Schur complement arising from the `MATIS` inside the `PCNN` preconditioner

    497:   Input Parameters:
    498: + pc     - preconditioner context
    499: . v      - vector to which the Schur complement is to be applied (it is NOT modified inside this function, UNLESS vec2_B is null)
    500: . vec1_B - location to store the result of Schur complement applied to chunk
    501: . vec2_B - workspace or `NULL`, `v` is used as workspace in that case
    502: . vec1_D - work space
    503: - vec2_D - work space

    505:   Level: advanced

    507: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
    508:           `PCISSetSubdomainScalingFactor()`, `PCISApplyInvSchur()`,
    509:           `PCISReset()`, `PCISInitialize()`
    510: @*/
    511: PetscErrorCode PCISApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D)
    512: {
    513:   PC_IS *pcis = (PC_IS *)pc->data;

    515:   PetscFunctionBegin;
    516:   if (!vec2_B) vec2_B = v;

    518:   PetscCall(MatMult(pcis->A_BB, v, vec1_B));
    519:   PetscCall(MatMult(pcis->A_IB, v, vec1_D));
    520:   PetscCall(KSPSolve(pcis->ksp_D, vec1_D, vec2_D));
    521:   PetscCall(KSPCheckSolve(pcis->ksp_D, pc, vec2_D));
    522:   PetscCall(MatMult(pcis->A_BI, vec2_D, vec2_B));
    523:   PetscCall(VecAXPY(vec1_B, -1.0, vec2_B));
    524:   PetscFunctionReturn(PETSC_SUCCESS);
    525: }

    527: /*@
    528:   PCISScatterArrayNToVecB - Scatters interface node values from a big array (of all local nodes, interior or interface,
    529:   including ghosts) into an interface vector, when in `SCATTER_FORWARD` mode, or vice-versa, when in `SCATTER_REVERSE`
    530:   mode.

    532:   Input Parameters:
    533: + pc      - preconditioner context
    534: . array_N - [when in `SCATTER_FORWARD` mode] Array to be scattered into the vector otherwise output array
    535: . imode   - insert mode, `ADD_VALUES` or `INSERT_VALUES`
    536: . smode   - scatter mode, `SCATTER_FORWARD` or `SCATTER_REVERSE` mode]
    537: - v_B     - [when in `SCATTER_REVERSE` mode] Vector to be scattered into the array, otherwise output vector

    539:   Level: advanced

    541:   Note:
    542:   The entries in the array that do not correspond to interface nodes remain unaltered.

    544: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`,
    545:           `PCISSetSubdomainScalingFactor()`, `PCISApplySchur()`, `PCISApplyInvSchur()`,
    546:           `PCISReset()`, `PCISInitialize()`, `InsertMode`
    547: @*/
    548: PetscErrorCode PCISScatterArrayNToVecB(PC pc, PetscScalar *array_N, Vec v_B, InsertMode imode, ScatterMode smode)
    549: {
    550:   PetscInt        i;
    551:   const PetscInt *idex;
    552:   PetscScalar    *array_B;
    553:   PC_IS          *pcis = (PC_IS *)pc->data;

    555:   PetscFunctionBegin;
    556:   PetscCall(VecGetArray(v_B, &array_B));
    557:   PetscCall(ISGetIndices(pcis->is_B_local, &idex));

    559:   if (smode == SCATTER_FORWARD) {
    560:     if (imode == INSERT_VALUES) {
    561:       for (i = 0; i < pcis->n_B; i++) array_B[i] = array_N[idex[i]];
    562:     } else { /* ADD_VALUES */
    563:       for (i = 0; i < pcis->n_B; i++) array_B[i] += array_N[idex[i]];
    564:     }
    565:   } else { /* SCATTER_REVERSE */
    566:     if (imode == INSERT_VALUES) {
    567:       for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] = array_B[i];
    568:     } else { /* ADD_VALUES */
    569:       for (i = 0; i < pcis->n_B; i++) array_N[idex[i]] += array_B[i];
    570:     }
    571:   }
    572:   PetscCall(ISRestoreIndices(pcis->is_B_local, &idex));
    573:   PetscCall(VecRestoreArray(v_B, &array_B));
    574:   PetscFunctionReturn(PETSC_SUCCESS);
    575: }

    577: /*@
    578:   PCISApplyInvSchur - Solves the Neumann problem related to applying the inverse of the Schur complement.

    580:   Input Parameters:
    581: + pc     - preconditioner context
    582: . b      - vector of local interface nodes (including ghosts)
    583: . x      - vector of local interface nodes (including ghosts); returns the application of the inverse of the Schur complement to `b`
    584: . vec1_N - vector of local nodes (interior and interface, including ghosts); used as work space
    585: - vec2_N - vector of local nodes (interior and interface, including ghosts); used as work space

    587:   Level: advanced

    589:   Note:
    590:   Solves the problem
    591: .vb
    592:   [ A_II  A_IB ] [ . ]   [ 0 ]
    593:   [            ] [   ] = [   ]
    594:   [ A_BI  A_BB ] [ x ]   [ b ]
    595: .ve

    597: .seealso: [](ch_ksp), `PCBDDC`, `PCNN`, `PCISSetUseStiffnessScaling()`, `PCISSetSubdomainDiagonalScaling()`, `PCISScatterArrayNToVecB()`,
    598:           `PCISSetSubdomainScalingFactor()`,
    599:           `PCISReset()`, `PCISInitialize()`
    600: @*/
    601: PetscErrorCode PCISApplyInvSchur(PC pc, Vec b, Vec x, Vec vec1_N, Vec vec2_N)
    602: {
    603:   PC_IS *pcis = (PC_IS *)pc->data;

    605:   PetscFunctionBegin;
    606:   /*
    607:     Neumann solvers.
    608:     Applying the inverse of the local Schur complement, i.e, solving a Neumann
    609:     Problem with zero at the interior nodes of the RHS and extracting the interface
    610:     part of the solution. inverse Schur complement is applied to b and the result
    611:     is stored in x.
    612:   */
    613:   /* Setting the RHS vec1_N */
    614:   PetscCall(VecSet(vec1_N, 0.0));
    615:   PetscCall(VecScatterBegin(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE));
    616:   PetscCall(VecScatterEnd(pcis->N_to_B, b, vec1_N, INSERT_VALUES, SCATTER_REVERSE));
    617:   /* Checking for consistency of the RHS */
    618:   {
    619:     PetscBool flg = PETSC_FALSE;
    620:     PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc_is_check_consistency", &flg, NULL));
    621:     if (flg) {
    622:       PetscScalar average;
    623:       PetscViewer viewer;
    624:       PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc), &viewer));

    626:       PetscCall(VecSum(vec1_N, &average));
    627:       average = average / ((PetscReal)pcis->n);
    628:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
    629:       if (pcis->pure_neumann) {
    630:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is floating. Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average)));
    631:       } else {
    632:         PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Subdomain %04d is fixed.    Average = % 1.14e\n", PetscGlobalRank, (double)PetscAbsScalar(average)));
    633:       }
    634:       PetscCall(PetscViewerFlush(viewer));
    635:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
    636:     }
    637:   }
    638:   /* Solving the system for vec2_N */
    639:   PetscCall(KSPSolve(pcis->ksp_N, vec1_N, vec2_N));
    640:   PetscCall(KSPCheckSolve(pcis->ksp_N, pc, vec2_N));
    641:   /* Extracting the local interface vector out of the solution */
    642:   PetscCall(VecScatterBegin(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD));
    643:   PetscCall(VecScatterEnd(pcis->N_to_B, vec2_N, x, INSERT_VALUES, SCATTER_FORWARD));
    644:   PetscFunctionReturn(PETSC_SUCCESS);
    645: }