C++ IBTK_CHKERRQ函数代码示例

void
PenaltyIBMethod::preprocessIntegrateData(double current_time, double new_time, int num_cycles)
{
    IBMethod::preprocessIntegrateData(current_time, new_time, num_cycles);

    const int coarsest_ln = 0;
    const int finest_ln = d_hierarchy->getFinestLevelNumber();

    // Look-up or allocate Lagangian data.
    d_K_data.resize(finest_ln + 1);
    d_M_data.resize(finest_ln + 1);
    d_Y_current_data.resize(finest_ln + 1);
    d_Y_new_data.resize(finest_ln + 1);
    d_V_current_data.resize(finest_ln + 1);
    d_V_new_data.resize(finest_ln + 1);
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;
        d_K_data[ln] = d_l_data_manager->getLData("K", ln);
        d_M_data[ln] = d_l_data_manager->getLData("M", ln);
        d_Y_current_data[ln] = d_l_data_manager->getLData("Y", ln);
        d_Y_new_data[ln] = d_l_data_manager->createLData("Y_new", ln, NDIM);
        d_V_current_data[ln] = d_l_data_manager->getLData("V", ln);
        d_V_new_data[ln] = d_l_data_manager->createLData("V_new", ln, NDIM);

        // Initialize Y^{n+1} and V^{n+1} to equal Y^{n} and V^{n}.
        int ierr;
        ierr = VecCopy(d_Y_current_data[ln]->getVec(), d_Y_new_data[ln]->getVec());
        IBTK_CHKERRQ(ierr);
        ierr = VecCopy(d_V_current_data[ln]->getVec(), d_V_new_data[ln]->getVec());
        IBTK_CHKERRQ(ierr);
    }
    return;
} // preprocessIntegrateData

void
IBImplicitModHelmholtzPETScLevelSolver::deallocateSolverState()
{
    if (!d_is_initialized) return;

    IBAMR_TIMER_START(t_deallocate_solver_state);

    // Deallocate PETSc objects.
    int ierr;
    ierr = KSPDestroy(d_petsc_ksp); IBTK_CHKERRQ(ierr);
    ierr = MatDestroy(d_petsc_mat); IBTK_CHKERRQ(ierr);
    ierr = VecDestroy(d_petsc_x); IBTK_CHKERRQ(ierr);
    ierr = VecDestroy(d_petsc_b); IBTK_CHKERRQ(ierr);
    d_dof_index_fill.setNull();

    d_petsc_ksp = PETSC_NULL;
    d_petsc_mat = PETSC_NULL;
    d_petsc_x = PETSC_NULL;
    d_petsc_b = PETSC_NULL;

    // Deallocate DOF index data.
    Pointer<PatchLevel<NDIM> > level = d_hierarchy->getPatchLevel(d_level_num);
    if (level->checkAllocated(d_dof_index_idx)) level->deallocatePatchData(d_dof_index_idx);

    // Indicate that the solver is NOT initialized.
    d_is_initialized = false;

    IBAMR_TIMER_STOP(t_deallocate_solver_state);
    return;
}// deallocateSolverState

void
PenaltyIBMethod::postprocessIntegrateData(double current_time, double new_time, int num_cycles)
{
    IBMethod::postprocessIntegrateData(current_time, new_time, num_cycles);

    const int coarsest_ln = 0;
    const int finest_ln = d_hierarchy->getFinestLevelNumber();

    // Reset time-dependent Lagrangian data.
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;
        int ierr;
        ierr = VecSwap(d_Y_current_data[ln]->getVec(), d_Y_new_data[ln]->getVec());
        IBTK_CHKERRQ(ierr);
        ierr = VecSwap(d_V_current_data[ln]->getVec(), d_V_new_data[ln]->getVec());
        IBTK_CHKERRQ(ierr);
    }

    // Deallocate Lagrangian scratch data.
    d_K_data.clear();
    d_M_data.clear();
    d_Y_current_data.clear();
    d_Y_new_data.clear();
    d_V_current_data.clear();
    d_V_new_data.clear();
    return;
} // postprocessIntegrateData

void CCPoissonPETScLevelSolver::initializeSolverStateSpecialized(const SAMRAIVectorReal<NDIM, double>& x,
                                                                 const SAMRAIVectorReal<NDIM, double>& /*b*/)
{
    // Allocate DOF index data.
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
    const int x_idx = x.getComponentDescriptorIndex(0);
    Pointer<CellDataFactory<NDIM, double> > x_fac = var_db->getPatchDescriptor()->getPatchDataFactory(x_idx);
    const int depth = x_fac->getDefaultDepth();
    Pointer<CellDataFactory<NDIM, int> > dof_index_fac =
        var_db->getPatchDescriptor()->getPatchDataFactory(d_dof_index_idx);
    dof_index_fac->setDefaultDepth(depth);
    Pointer<PatchLevel<NDIM> > level = d_hierarchy->getPatchLevel(d_level_num);
    if (!level->checkAllocated(d_dof_index_idx)) level->allocatePatchData(d_dof_index_idx);

    // Setup PETSc objects.
    int ierr;
    PETScVecUtilities::constructPatchLevelDOFIndices(d_num_dofs_per_proc, d_dof_index_idx, level);
    const int mpi_rank = SAMRAI_MPI::getRank();
    ierr = VecCreateMPI(PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_x);
    IBTK_CHKERRQ(ierr);
    ierr = VecCreateMPI(PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_b);
    IBTK_CHKERRQ(ierr);
    PETScMatUtilities::constructPatchLevelCCLaplaceOp(
        d_petsc_mat, d_poisson_spec, d_bc_coefs, d_solution_time, d_num_dofs_per_proc, d_dof_index_idx, level);
    d_petsc_pc = d_petsc_mat;
    d_petsc_ksp_ops_flag = SAME_PRECONDITIONER;
    d_data_synch_sched = PETScVecUtilities::constructDataSynchSchedule(x_idx, level);
    d_ghost_fill_sched = PETScVecUtilities::constructGhostFillSchedule(x_idx, level);
    return;
} // initializeSolverStateSpecialized

PetscErrorCode IBImplicitStaggeredHierarchyIntegrator::compositeIBJacobianApply(Vec x, Vec f)
{
    PetscErrorCode ierr;
    const double half_time = d_integrator_time + 0.5 * d_current_dt;

    Vec* component_sol_vecs;
    Vec* component_rhs_vecs;
    ierr = VecMultiVecGetSubVecs(x, &component_sol_vecs);
    IBTK_CHKERRQ(ierr);
    ierr = VecMultiVecGetSubVecs(f, &component_rhs_vecs);
    IBTK_CHKERRQ(ierr);

    Pointer<SAMRAIVectorReal<NDIM, double> > u =
        PETScSAMRAIVectorReal::getSAMRAIVector(component_sol_vecs[0]);
    Pointer<SAMRAIVectorReal<NDIM, double> > f_u =
        PETScSAMRAIVectorReal::getSAMRAIVector(component_rhs_vecs[0]);

    Pointer<Variable<NDIM> > u_var = d_ins_hier_integrator->getVelocityVariable();
    const int u_idx = u->getComponentDescriptorIndex(0);
    const int f_u_idx = f_u->getComponentDescriptorIndex(0);

    Vec X = component_sol_vecs[1];
    Vec R = component_rhs_vecs[1];

    // Evaluate the Eulerian terms.
    d_stokes_op->setHomogeneousBc(true);
    d_stokes_op->apply(*u, *f_u);

    d_ib_implicit_ops->computeLinearizedLagrangianForce(X, half_time);
    if (d_enable_logging)
        plog << d_object_name
             << "::integrateHierarchy(): spreading Lagrangian force to the Eulerian grid\n";
    d_hier_velocity_data_ops->setToScalar(d_f_idx, 0.0);
    d_u_phys_bdry_op->setPatchDataIndex(d_f_idx);
    d_ib_implicit_ops->spreadLinearizedForce(d_f_idx,
                                             d_u_phys_bdry_op,
                                             getProlongRefineSchedules(d_object_name + "::f"),
                                             half_time);
    d_hier_velocity_data_ops->subtract(f_u_idx, f_u_idx, d_f_idx);
    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(component_rhs_vecs[0]));
    IBTK_CHKERRQ(ierr);

    // Evaluate the Lagrangian terms.
    d_hier_velocity_data_ops->scale(d_u_idx, 0.5, u_idx);
    d_u_phys_bdry_op->setPatchDataIndex(d_u_idx);
    d_ib_implicit_ops->interpolateLinearizedVelocity(
        d_u_idx,
        getCoarsenSchedules(d_object_name + "::u::CONSERVATIVE_COARSEN"),
        getGhostfillRefineSchedules(d_object_name + "::u"),
        half_time);
    d_ib_implicit_ops->computeLinearizedResidual(X, R);

    // Ensure that PETSc sees that the state of the RHS vector has changed.
    // This is a nasty hack.
    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(f));
    IBTK_CHKERRQ(ierr);
    return ierr;
} // compositeIBJacobianApply

bool
PETScKrylovLinearSolver::solveSystem(SAMRAIVectorReal<NDIM, double>& x, SAMRAIVectorReal<NDIM, double>& b)
{
    IBTK_TIMER_START(t_solve_system);

#if !defined(NDEBUG)
    TBOX_ASSERT(d_A);
#endif
    int ierr;

    // Initialize the solver, when necessary.
    const bool deallocate_after_solve = !d_is_initialized;
    if (deallocate_after_solve) initializeSolverState(x, b);
#if !defined(NDEBUG)
    TBOX_ASSERT(d_petsc_ksp);
#endif
    resetKSPOptions();

    // Allocate scratch data.
    d_b->allocateVectorData();

    // Solve the system using a PETSc KSP object.
    d_b->copyVector(Pointer<SAMRAIVectorReal<NDIM, double> >(&b, false));
    d_A->setHomogeneousBc(d_homogeneous_bc);
    d_A->modifyRhsForBcs(*d_b);
    d_A->setHomogeneousBc(true);
    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_x, Pointer<SAMRAIVectorReal<NDIM, double> >(&x, false));
    PETScSAMRAIVectorReal::replaceSAMRAIVector(d_petsc_b, d_b);
    ierr = KSPSolve(d_petsc_ksp, d_petsc_b, d_petsc_x);
    IBTK_CHKERRQ(ierr);
    d_A->setHomogeneousBc(d_homogeneous_bc);
    d_A->imposeSolBcs(x);

    // Get iterations count and residual norm.
    ierr = KSPGetIterationNumber(d_petsc_ksp, &d_current_iterations);
    IBTK_CHKERRQ(ierr);
    ierr = KSPGetResidualNorm(d_petsc_ksp, &d_current_residual_norm);
    IBTK_CHKERRQ(ierr);
    d_A->setHomogeneousBc(d_homogeneous_bc);

    // Determine the convergence reason.
    KSPConvergedReason reason;
    ierr = KSPGetConvergedReason(d_petsc_ksp, &reason);
    IBTK_CHKERRQ(ierr);
    const bool converged = (static_cast<int>(reason) > 0);
    if (d_enable_logging) reportKSPConvergedReason(reason, plog);

    // Dealocate scratch data.
    d_b->deallocateVectorData();

    // Deallocate the solver, when necessary.
    if (deallocate_after_solve) deallocateSolverState();

    IBTK_TIMER_STOP(t_solve_system);
    return converged;
} // solveSystem

PetscErrorCode
IBImplicitStaggeredHierarchyIntegrator::compositeIBPCApply_SAMRAI(PC pc, Vec x, Vec y)
{
    PetscErrorCode ierr;
    void* ctx;
    ierr = PCShellGetContext(pc, &ctx);
    IBTK_CHKERRQ(ierr);
    IBImplicitStaggeredHierarchyIntegrator* ib_integrator =
        static_cast<IBImplicitStaggeredHierarchyIntegrator*>(ctx);
    ierr = ib_integrator->compositeIBPCApply(x, y);
    IBTK_CHKERRQ(ierr);
    return ierr;
} // compositeIBPCApply_SAMRAI

void GeneralizedIBMethod::preprocessIntegrateData(double current_time,
                                                  double new_time,
                                                  int num_cycles)
{
    d_ib_force_and_torque_fcn_needs_init =
        d_ib_force_fcn_needs_init || d_ib_force_and_torque_fcn_needs_init;
    IBMethod::preprocessIntegrateData(current_time, new_time, num_cycles);

    const int coarsest_ln = 0;
    const int finest_ln = d_hierarchy->getFinestLevelNumber();
    const double start_time = d_ib_solver->getStartTime();

    if (d_ib_force_and_torque_fcn)
    {
        if (d_ib_force_and_torque_fcn_needs_init)
        {
            const bool initial_time =
                MathUtilities<double>::equalEps(current_time, start_time);
            resetLagrangianForceAndTorqueFunction(current_time, initial_time);
            d_ib_force_and_torque_fcn_needs_init = false;
        }
    }

    // Look-up or allocate Lagangian data.
    d_D_current_data.resize(finest_ln + 1);
    d_D_new_data.resize(finest_ln + 1);
    d_N_current_data.resize(finest_ln + 1);
    d_N_new_data.resize(finest_ln + 1);
    d_W_current_data.resize(finest_ln + 1);
    d_W_new_data.resize(finest_ln + 1);
    for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
    {
        if (!d_l_data_manager->levelContainsLagrangianData(ln)) continue;
        d_D_current_data[ln] = d_l_data_manager->getLData("D", ln);
        d_D_new_data[ln] = d_l_data_manager->createLData("D_new", ln, NDIM * NDIM);
        d_N_current_data[ln] = d_l_data_manager->createLData("N", ln, NDIM);
        d_N_new_data[ln] = d_l_data_manager->createLData("N_new", ln, NDIM);
        d_W_current_data[ln] = d_l_data_manager->getLData("W", ln);
        d_W_new_data[ln] = d_l_data_manager->createLData("W_new", ln, NDIM);

        // Initialize D^{n+1} to equal D^{n}, and initialize W^{n+1} to equal
        // W^{n}.
        int ierr;
        ierr = VecCopy(d_D_current_data[ln]->getVec(), d_D_new_data[ln]->getVec());
        IBTK_CHKERRQ(ierr);
        ierr = VecCopy(d_W_current_data[ln]->getVec(), d_W_new_data[ln]->getVec());
        IBTK_CHKERRQ(ierr);
    }
    return;
} // preprocessIntegrateData

void StaggeredStokesPETScLevelSolver::initializeSolverStateSpecialized(
    const SAMRAIVectorReal<NDIM, double>& x,
    const SAMRAIVectorReal<NDIM, double>& /*b*/)
{
    // Allocate DOF index data.
    Pointer<PatchLevel<NDIM> > level = d_hierarchy->getPatchLevel(d_level_num);
    if (!level->checkAllocated(d_u_dof_index_idx)) level->allocatePatchData(d_u_dof_index_idx);
    if (!level->checkAllocated(d_p_dof_index_idx)) level->allocatePatchData(d_p_dof_index_idx);

    // Setup PETSc objects.
    int ierr;
    StaggeredStokesPETScVecUtilities::constructPatchLevelDOFIndices(
        d_num_dofs_per_proc, d_u_dof_index_idx, d_p_dof_index_idx, level);
    const int mpi_rank = SAMRAI_MPI::getRank();
    ierr = VecCreateMPI(
               PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_x);
    IBTK_CHKERRQ(ierr);
    ierr = VecCreateMPI(
               PETSC_COMM_WORLD, d_num_dofs_per_proc[mpi_rank], PETSC_DETERMINE, &d_petsc_b);
    IBTK_CHKERRQ(ierr);
    StaggeredStokesPETScMatUtilities::constructPatchLevelMACStokesOp(d_petsc_mat,
            d_U_problem_coefs,
            d_U_bc_coefs,
            d_new_time,
            d_num_dofs_per_proc,
            d_u_dof_index_idx,
            d_p_dof_index_idx,
            level);
    ierr = MatDuplicate(d_petsc_mat, MAT_COPY_VALUES, &d_petsc_pc);
    IBTK_CHKERRQ(ierr);
    HierarchyDataOpsManager<NDIM>* hier_ops_manager =
        HierarchyDataOpsManager<NDIM>::getManager();
    Pointer<HierarchyDataOpsInteger<NDIM> > hier_p_dof_index_ops =
        hier_ops_manager->getOperationsInteger(d_p_dof_index_var, d_hierarchy, true);
    hier_p_dof_index_ops->resetLevels(d_level_num, d_level_num);
    const int min_p_idx = hier_p_dof_index_ops->min(
                              d_p_dof_index_idx); // NOTE: HierarchyDataOpsInteger::max() is broken
    ierr = MatZeroRowsColumns(d_petsc_pc, 1, &min_p_idx, 1.0, NULL, NULL);
    IBTK_CHKERRQ(ierr);
    d_petsc_ksp_ops_flag = SAME_PRECONDITIONER;
    const int u_idx = x.getComponentDescriptorIndex(0);
    const int p_idx = x.getComponentDescriptorIndex(1);
    d_data_synch_sched =
        StaggeredStokesPETScVecUtilities::constructDataSynchSchedule(u_idx, p_idx, level);
    d_ghost_fill_sched =
        StaggeredStokesPETScVecUtilities::constructGhostFillSchedule(u_idx, p_idx, level);
    return;
} // initializeSolverStateSpecialized

PetscErrorCode VecAYPX_SAMRAI(Vec y, const PetscScalar alpha, Vec x)
{
    IBTK_TIMER_START(t_vec_aypx);
#if !defined(NDEBUG)
    TBOX_ASSERT(x);
    TBOX_ASSERT(y);
#endif
    static const bool interior_only = false;
    if (MathUtilities<double>::equalEps(alpha, 1.0))
    {
        PSVR_CAST2(y)->add(PSVR_CAST2(x), PSVR_CAST2(y), interior_only);
    }
    else if (MathUtilities<double>::equalEps(alpha, -1.0))
    {
        PSVR_CAST2(y)->subtract(PSVR_CAST2(x), PSVR_CAST2(y), interior_only);
    }
    else
    {
        PSVR_CAST2(y)->axpy(alpha, PSVR_CAST2(y), PSVR_CAST2(x), interior_only);
    }
    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(y));
    IBTK_CHKERRQ(ierr);
    IBTK_TIMER_STOP(t_vec_aypx);
    PetscFunctionReturn(0);
} // VecAYPX

Pointer<SAMRAIVectorReal<NDIM, double> > PETScNewtonKrylovSolver::getFunctionVector() const
{
    Vec petsc_f;
    int ierr = SNESGetFunction(d_petsc_snes, &petsc_f, NULL, NULL);
    IBTK_CHKERRQ(ierr);
    return PETScSAMRAIVectorReal::getSAMRAIVector(petsc_f);
} // getFunctionVector

Pointer<SAMRAIVectorReal<NDIM, double> > PETScSNESJacobianJOWrapper::getBaseVector() const
{
    Vec petsc_x;
    int ierr = SNESGetSolution(d_petsc_snes, &petsc_x);
    IBTK_CHKERRQ(ierr);
    return PETScSAMRAIVectorReal::getSAMRAIVector(petsc_x);
} // getBaseVector

PetscErrorCode VecSetRandom_SAMRAI(Vec x, PetscRandom rctx)
{
    IBTK_TIMER_START(t_vec_set_random);
#if !defined(NDEBUG)
    TBOX_ASSERT(x);
#endif
    PetscScalar lo, hi;
    int ierr;
    ierr = PetscRandomGetInterval(rctx, &lo, &hi);
    IBTK_CHKERRQ(ierr);
    PSVR_CAST2(x)->setRandomValues(hi - lo, lo);
    ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(x));
    IBTK_CHKERRQ(ierr);
    IBTK_TIMER_STOP(t_vec_set_random);
    PetscFunctionReturn(0);
} // VecSetRandom

PetscErrorCode
VecMAXPY_SAMRAI(
    Vec y,
    PetscInt nv,
    const PetscScalar* alpha,
    Vec* x)
{
    IBTK_TIMER_START(t_vec_maxpy);
#ifdef DEBUG_CHECK_ASSERTIONS
    TBOX_ASSERT(y != PETSC_NULL);
    for (PetscInt i = 0; i < nv; ++i)
    {
        TBOX_ASSERT(x[i] != PETSC_NULL);
    }
#endif
    static const bool interior_only = false;
    for (PetscInt i = 0; i < nv; ++i)
    {
        if (MathUtilities<double>::equalEps(alpha[i],1.0))
        {
            PSVR_CAST2(y)->add(PSVR_CAST2(x[i]), PSVR_CAST2(y), interior_only);
        }
        else if (MathUtilities<double>::equalEps(alpha[i],-1.0))
        {
            PSVR_CAST2(y)->subtract(PSVR_CAST2(y), PSVR_CAST2(x[i]), interior_only);
        }
        else
        {
            PSVR_CAST2(y)->axpy(alpha[i], PSVR_CAST2(x[i]), PSVR_CAST2(y), interior_only);
        }
    }
    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(y)); IBTK_CHKERRQ(ierr);
    IBTK_TIMER_STOP(t_vec_maxpy);
    PetscFunctionReturn(0);
}// VecMAXPY

PetscErrorCode
VecWAXPY_SAMRAI(
    Vec w,
    PetscScalar alpha,
    Vec x,
    Vec y)
{
    IBTK_TIMER_START(t_vec_waxpy);
#ifdef DEBUG_CHECK_ASSERTIONS
    TBOX_ASSERT(x != PETSC_NULL);
    TBOX_ASSERT(y != PETSC_NULL);
    TBOX_ASSERT(w != PETSC_NULL);
#endif
    static const bool interior_only = false;
    if (MathUtilities<double>::equalEps(alpha,1.0))
    {
        PSVR_CAST2(w)->add(PSVR_CAST2(x), PSVR_CAST2(y), interior_only);
    }
    else if (MathUtilities<double>::equalEps(alpha,-1.0))
    {
        PSVR_CAST2(w)->subtract(PSVR_CAST2(y), PSVR_CAST2(x), interior_only);
    }
    else
    {
        PSVR_CAST2(w)->axpy(alpha, PSVR_CAST2(x), PSVR_CAST2(y), interior_only);
    }
    int ierr = PetscObjectStateIncrease(reinterpret_cast<PetscObject>(w)); IBTK_CHKERRQ(ierr);
    IBTK_TIMER_STOP(t_vec_waxpy);
    PetscFunctionReturn(0);
}// VecWAXPY