Method MatVecMult_KMInv

src/IB/KrylovMobilitySolver.cpp:863–941

} // resetKSPPC

PetscErrorCode
KrylovMobilitySolver::MatVecMult_KMInv(Mat A, Vec x, Vec y)
{
    PetscFunctionBeginUser;
    void* p_ctx;
    MatShellGetContext(A, &p_ctx);
    auto solver = static_cast<KrylovMobilitySolver*>(p_ctx);
    Pointer<IBStrategy> ib_method_ops = solver->d_cib_strategy;

#if !defined(NDEBUG)
    TBOX_ASSERT(solver);
    TBOX_ASSERT(solver->d_petsc_mat);
    TBOX_ASSERT(ib_method_ops);
#endif

    // Some constants
    static const double gamma = solver->d_scale_spread;
    static const double beta = solver->d_scale_interp;
    static const double delta = solver->d_reg_mob_factor;
    const double half_time = 0.5 * (solver->d_new_time + solver->d_current_time);

    // Use homogeneous BCs with Stokes solver in the preconditioner.
    dynamic_cast<IBTK::LinearSolver*>(solver->d_LInv.getPointer())->setHomogeneousBc(true);

    // Set y:= [J L^-1 S + \delta]x
    // 1) Spread force.
    solver->d_samrai_temp[0]->setToScalar(0.0);
    solver->d_cib_strategy->setConstraintForce(x, half_time, gamma);
    ib_method_ops->spreadForce(solver->d_samrai_temp[0]->getComponentDescriptorIndex(0),
                               nullptr,
                               std::vector<Pointer<RefineSchedule<NDIM>>>(),
                               half_time);
    if (solver->d_normalize_spread_force)
    {
        solver->d_cib_strategy->subtractMeanConstraintForce(
            x, solver->d_samrai_temp[0]->getComponentDescriptorIndex(0), gamma);
    }
    // 2) Solve Stokes system.
    solver->d_LInv->solveSystem(*solver->d_samrai_temp[1], *solver->d_samrai_temp[0]);

    // 3a) Fill velocity ghost cells.
    int u_data_idx = solver->d_samrai_temp[1]->getComponentDescriptorIndex(0);
    using InterpolationTransactionComponent = IBTK::HierarchyGhostCellInterpolation::InterpolationTransactionComponent;
    std::vector<InterpolationTransactionComponent> transaction_comps;
    InterpolationTransactionComponent u_component(u_data_idx,
                                                  DATA_REFINE_TYPE,
                                                  USE_CF_INTERPOLATION,
                                                  DATA_COARSEN_TYPE,
                                                  BDRY_EXTRAP_TYPE,
                                                  CONSISTENT_TYPE_2_BDRY,
                                                  solver->d_u_bc_coefs,
                                                  solver->d_fill_pattern);
    transaction_comps.push_back(u_component);
    solver->d_hier_bdry_fill->resetTransactionComponents(transaction_comps);
    static const bool homogeneous_bc = true;
    solver->d_hier_bdry_fill->setHomogeneousBc(homogeneous_bc);
    solver->d_hier_bdry_fill->fillData(half_time);
    solver->d_hier_bdry_fill->resetTransactionComponents(solver->d_transaction_comps);