C++ SAMRAIVectorReal类代码示例

void
INSStaggeredVCStokesOperator::modifyRhsForInhomogeneousBc(
    SAMRAIVectorReal<NDIM,double>& y)
{
    // Set y := y - A*0, i.e., shift the right-hand-side vector to account for
    // inhomogeneous boundary conditions.
    if (!d_homogeneous_bc)
    {
        d_correcting_rhs = true;

        Pointer<SAMRAIVectorReal<NDIM,double> > x = y.cloneVector("");
        x->allocateVectorData();
        x->setToScalar(0.0);

        Pointer<SAMRAIVectorReal<NDIM,double> > b = y.cloneVector("");
        b->allocateVectorData();
        b->setToScalar(0.0);

        apply(*x,*b);
        y.subtract(Pointer<SAMRAIVectorReal<NDIM,double> >(&y, false), b);

        x->freeVectorComponents();
        x.setNull();

        b->freeVectorComponents();
        b.setNull();

        d_correcting_rhs = false;
    }
    return;
}// modifyRhsForInhomogeneousBc

void CCPoissonPETScLevelSolver::setupKSPVecs(Vec& petsc_x,
                                             Vec& petsc_b,
                                             SAMRAIVectorReal<NDIM, double>& x,
                                             SAMRAIVectorReal<NDIM, double>& b,
                                             Pointer<PatchLevel<NDIM> > patch_level)
{
    if (!d_initial_guess_nonzero) copyToPETScVec(petsc_x, x, patch_level);
    const int b_idx = b.getComponentDescriptorIndex(0);
    Pointer<CellVariable<NDIM, double> > b_var = b.getComponentVariable(0);
    VariableDatabase<NDIM>* var_db = VariableDatabase<NDIM>::getDatabase();
    int b_adj_idx = var_db->registerClonedPatchDataIndex(b_var, b_idx);
    patch_level->allocatePatchData(b_adj_idx);
    for (PatchLevel<NDIM>::Iterator p(patch_level); p; p++)
    {
        Pointer<Patch<NDIM> > patch = patch_level->getPatch(p());
        Pointer<CellData<NDIM, double> > b_data = patch->getPatchData(b_idx);
        Pointer<CellData<NDIM, double> > b_adj_data = patch->getPatchData(b_adj_idx);
        b_adj_data->copy(*b_data);
        if (!patch->getPatchGeometry()->intersectsPhysicalBoundary()) continue;
        PoissonUtilities::adjustCCBoundaryRhsEntries(
            patch, *b_adj_data, d_poisson_spec, d_bc_coefs, d_solution_time, d_homogeneous_bc);
    }
    PETScVecUtilities::copyToPatchLevelVec(petsc_b, b_adj_idx, d_dof_index_idx, patch_level);
    patch_level->deallocatePatchData(b_adj_idx);
    var_db->removePatchDataIndex(b_adj_idx);
    return;
} // setupKSPVecs

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

    int ierr;

    if (d_enable_logging) plog << d_object_name << "::solveSystem():" << std::endl;

    // Initialize the solver, when necessary.
    const bool deallocate_after_solve = !d_is_initialized;
    if (deallocate_after_solve) initializeSolverState(x,b);
#if 0 // XXXX
    // Configure solver.
    ierr = KSPSetTolerances(d_petsc_ksp, d_rel_residual_tol, d_abs_residual_tol, PETSC_DEFAULT, d_max_iterations); IBTK_CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(d_petsc_ksp, d_initial_guess_nonzero ? PETSC_TRUE : PETSC_FALSE); IBTK_CHKERRQ(ierr);
#endif

    // Solve the system.
    Pointer<PatchLevel<NDIM> > patch_level = d_hierarchy->getPatchLevel(d_level_num);
    const int x_idx = x.getComponentDescriptorIndex(0);
    Pointer<SideVariable<NDIM,double> > x_var = x.getComponentVariable(0);
    const int b_idx = b.getComponentDescriptorIndex(0);
    Pointer<SideVariable<NDIM,double> > b_var = b.getComponentVariable(0);
    if (d_initial_guess_nonzero) PETScVecUtilities::copyToPatchLevelVec(d_petsc_x, x_idx, x_var, patch_level);
    PETScVecUtilities::copyToPatchLevelVec(d_petsc_b, b_idx, b_var, patch_level);
    PETScVecUtilities::constrainPatchLevelVec(d_petsc_b, d_dof_index_idx, d_dof_index_var, patch_level, d_dof_index_fill);
    ierr = KSPSolve(d_petsc_ksp, d_petsc_b, d_petsc_x); IBTK_CHKERRQ(ierr);
    PETScVecUtilities::copyFromPatchLevelVec(d_petsc_x, x_idx, x_var, patch_level);

    typedef SideDataSynchronization::SynchronizationTransactionComponent SynchronizationTransactionComponent; // XXXX
    SynchronizationTransactionComponent x_synch_transaction = SynchronizationTransactionComponent(x_idx, "CONSERVATIVE_COARSEN");
    Pointer<SideDataSynchronization> side_synch_op = new SideDataSynchronization();
    side_synch_op->initializeOperatorState(x_synch_transaction, x.getPatchHierarchy());
    side_synch_op->synchronizeData(0.0);

    // Log solver info.
    KSPConvergedReason reason;
    ierr = KSPGetConvergedReason(d_petsc_ksp, &reason); IBTK_CHKERRQ(ierr);
    const bool converged = reason > 0;
    if (d_enable_logging)
    {
        plog << d_object_name << "::solveSystem(): solver " << (converged ? "converged" : "diverged") << "\n"
             << "iterations = " << d_current_its << "\n"
             << "residual norm = " << d_current_residual_norm << std::endl;
    }

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

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

void ConvectiveOperator::apply(SAMRAIVectorReal<NDIM, double>& x,
                               SAMRAIVectorReal<NDIM, double>& y)
{
    // Get the vector components.
    const int Q_idx = x.getComponentDescriptorIndex(0);
    const int N_idx = y.getComponentDescriptorIndex(0);

    // Compute the action of the operator.
    applyConvectiveOperator(Q_idx, N_idx);
    return;
} // apply

void PETScSNESFunctionGOWrapper::initializeOperatorState(const SAMRAIVectorReal<NDIM, double>& in,
                                                         const SAMRAIVectorReal<NDIM, double>& out)
{
    if (d_is_initialized) deallocateOperatorState();
    d_x = in.cloneVector("");
    d_y = out.cloneVector("");
    MPI_Comm comm;
    int ierr = PetscObjectGetComm(reinterpret_cast<PetscObject>(d_petsc_snes), &comm);
    IBTK_CHKERRQ(ierr);
    d_petsc_x = PETScSAMRAIVectorReal::createPETScVector(d_x, comm);
    d_petsc_y = PETScSAMRAIVectorReal::createPETScVector(d_y, comm);
    d_is_initialized = true;
    return;
} // initializeOperatorState

Pointer<SAMRAIVectorReal<NDIM,double> >
FACPreconditionerStrategy::getLevelSAMRAIVectorReal(
    const SAMRAIVectorReal<NDIM,double>& vec,
    int level_num) const
{
    std::ostringstream name_str;
    name_str << vec.getName() << "::level_" << level_num;
    Pointer<SAMRAIVectorReal<NDIM,double> > level_vec = new SAMRAIVectorReal<NDIM,double>(name_str.str(), vec.getPatchHierarchy(), level_num, level_num);
    for (int comp = 0; comp < vec.getNumberOfComponents(); ++comp)
    {
        level_vec->addComponent(vec.getComponentVariable(comp), vec.getComponentDescriptorIndex(comp), vec.getControlVolumeIndex(comp));
    }
    return level_vec;
}// getLevelSAMRAIVectorReal

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

    if (d_enable_logging) plog << d_object_name << "::solveSystem():" << std::endl;

    // Initialize the solver, when necessary.
    const bool deallocate_after_solve = !d_is_initialized;
    if (deallocate_after_solve) initializeSolverState(x,b);

    // Solve the system using the hypre solver.
    static const int comp = 0;
    const int x_idx = x.getComponentDescriptorIndex(comp);
    const int b_idx = b.getComponentDescriptorIndex(comp);

    bool converged = true;
    IntVector<NDIM> chkbrd_mode_id;
#if (NDIM > 2)
    for (chkbrd_mode_id(2) = 0; chkbrd_mode_id(2) < 2; ++chkbrd_mode_id(2))
    {
#endif
        for (chkbrd_mode_id(1) = 0; chkbrd_mode_id(1) < 2; ++chkbrd_mode_id(1))
        {
            for (chkbrd_mode_id(0) = 0; chkbrd_mode_id(0) < 2; ++chkbrd_mode_id(0))
            {
                bool converged_mode = solveSystem(x_idx, b_idx, chkbrd_mode_id);
                if (d_enable_logging)
                {
                    plog << d_object_name << "::solveSystem(): solver " << (converged_mode ? "converged" : "diverged") << "\n"
                         << "chkbrd_mode_id = " << chkbrd_mode_id << "\n"
                         << "iterations = " << d_current_its << "\n"
                         << "residual norm = " << d_current_residual_norm << std::endl;
                }
                converged = converged && converged_mode;
            }
        }
#if (NDIM > 2)
    }
#endif

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

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

void GeneralOperator::applyAdd(SAMRAIVectorReal<NDIM, double>& x,
                               SAMRAIVectorReal<NDIM, double>& y,
                               SAMRAIVectorReal<NDIM, double>& z)
{
    // Guard against the case that y == z.
    Pointer<SAMRAIVectorReal<NDIM, double> > zz = z.cloneVector(z.getName());
    zz->allocateVectorData();
    zz->copyVector(Pointer<SAMRAIVectorReal<NDIM, double> >(&z, false));
    apply(x, *zz);
    z.add(Pointer<SAMRAIVectorReal<NDIM, double> >(&y, false), zz);
    zz->deallocateVectorData();
    zz->freeVectorComponents();
    zz.setNull();
    return;
} // applyAdd

void
PETScPCLSWrapper::initializeSolverState(const SAMRAIVectorReal<NDIM, double>& x,
                                        const SAMRAIVectorReal<NDIM, double>& b)
{
    if (d_is_initialized) deallocateSolverState();
    d_x = x.cloneVector("");
    d_b = b.cloneVector("");
    MPI_Comm comm;
    int ierr = PetscObjectGetComm(reinterpret_cast<PetscObject>(d_petsc_pc), &comm);
    IBTK_CHKERRQ(ierr);
    d_petsc_x = PETScSAMRAIVectorReal::createPETScVector(d_x, comm);
    d_petsc_b = PETScSAMRAIVectorReal::createPETScVector(d_b, comm);
    d_is_initialized = true;
    return;
} // initializeSolverState

void
SCPoissonPETScLevelSolver::copyToPETScVec(Vec& petsc_x, SAMRAIVectorReal<NDIM, double>& x)
{
    const int x_idx = x.getComponentDescriptorIndex(0);
    PETScVecUtilities::copyToPatchLevelVec(petsc_x, x_idx, d_dof_index_idx, d_level);
    return;
} // copyToPETScVec

void
INSStaggeredVCStokesOperator::initializeOperatorState(
    const SAMRAIVectorReal<NDIM,double>& in,
    const SAMRAIVectorReal<NDIM,double>& /*out*/)
{
    IBAMR_TIMER_START(t_initialize_operator_state);

    if (d_is_initialized) deallocateOperatorState();

    d_x_scratch = in.cloneVector("INSStaggeredVCStokesOperator::x_scratch");
    d_x_scratch->allocateVectorData();

    typedef HierarchyGhostCellInterpolation::InterpolationTransactionComponent InterpolationTransactionComponent;
    InterpolationTransactionComponent U_scratch_component(d_x_scratch->getComponentDescriptorIndex(0), U_DATA_COARSEN_TYPE, BDRY_EXTRAP_TYPE, CONSISTENT_TYPE_2_BDRY);
    InterpolationTransactionComponent P_scratch_component(d_x_scratch->getComponentDescriptorIndex(1), P_DATA_COARSEN_TYPE, BDRY_EXTRAP_TYPE, CONSISTENT_TYPE_2_BDRY);
    InterpolationTransactionComponent mu_component(d_mu_data_idx, MU_DATA_COARSEN_TYPE, BDRY_EXTRAP_TYPE, CONSISTENT_TYPE_2_BDRY);

    std::vector<InterpolationTransactionComponent> U_P_MU_components(3);
    U_P_MU_components[0] = U_scratch_component;
    U_P_MU_components[1] = P_scratch_component;
    U_P_MU_components[2] = mu_component;

    d_U_P_MU_bdry_fill_op = new HierarchyGhostCellInterpolation();
    d_U_P_MU_bdry_fill_op->initializeOperatorState(U_P_MU_components, d_x_scratch->getPatchHierarchy());

    d_is_initialized = true;

    IBAMR_TIMER_STOP(t_initialize_operator_state);
    return;
}// initializeOperatorState

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

void
LinearOperator::modifyRhsForBcs(SAMRAIVectorReal<NDIM, double>& y)
{
    if (d_homogeneous_bc) return;

    // Set y := y - A*0, i.e., shift the right-hand-side vector to account for
    // inhomogeneous boundary conditions.
    Pointer<SAMRAIVectorReal<NDIM, double> > x = y.cloneVector("");
    Pointer<SAMRAIVectorReal<NDIM, double> > b = y.cloneVector("");
    x->allocateVectorData();
    b->allocateVectorData();
    x->setToScalar(0.0);
    apply(*x, *b);
    y.subtract(Pointer<SAMRAIVectorReal<NDIM, double> >(&y, false), b);
    x->freeVectorComponents();
    b->freeVectorComponents();
    return;
} // modifyRhsForBcs

void
VCSCViscousPETScLevelSolver::setupKSPVecs(Vec& petsc_x,
                                          Vec& petsc_b,
                                          SAMRAIVectorReal<NDIM, double>& x,
                                          SAMRAIVectorReal<NDIM, double>& b)
{
    if (d_initial_guess_nonzero) copyToPETScVec(petsc_x, x);
    const bool level_zero = (d_level_num == 0);
    const int x_idx = x.getComponentDescriptorIndex(0);
    const int b_idx = b.getComponentDescriptorIndex(0);
    const auto b_adj_idx = d_cached_eulerian_data.getCachedPatchDataIndex(b_idx);
    for (PatchLevel<NDIM>::Iterator p(d_level); p; p++)
    {
        Pointer<Patch<NDIM> > patch = d_level->getPatch(p());
        Pointer<PatchGeometry<NDIM> > pgeom = patch->getPatchGeometry();
        Pointer<SideData<NDIM, double> > x_data = patch->getPatchData(x_idx);
        Pointer<SideData<NDIM, double> > b_data = patch->getPatchData(b_idx);
        Pointer<SideData<NDIM, double> > b_adj_data = patch->getPatchData(b_adj_idx);
        b_adj_data->copy(*b_data);
        const bool at_physical_bdry = pgeom->intersectsPhysicalBoundary();
        if (at_physical_bdry)
        {
            PoissonUtilities::adjustVCSCViscousOpRHSAtPhysicalBoundary(*b_adj_data,
                                                                       patch,
                                                                       d_poisson_spec,
                                                                       1.0,
                                                                       d_bc_coefs,
                                                                       d_solution_time,
                                                                       d_homogeneous_bc,
                                                                       d_mu_interp_type);
        }
        const Array<BoundaryBox<NDIM> >& type_1_cf_bdry =
            level_zero ? Array<BoundaryBox<NDIM> >() :
                         d_cf_boundary->getBoundaries(patch->getPatchNumber(), /* boundary type */ 1, d_mu_interp_type);
        const bool at_cf_bdry = type_1_cf_bdry.size() > 0;
        if (at_cf_bdry)
        {
            PoissonUtilities::adjustVCSCViscousOpRHSAtCoarseFineBoundary(
                *b_adj_data, *x_data, patch, d_poisson_spec, 1.0, type_1_cf_bdry);
        }
    }
    PETScVecUtilities::copyToPatchLevelVec(petsc_b, b_adj_idx, d_dof_index_idx, d_level);
    return;
} // setupKSPVecs

void StaggeredStokesPETScLevelSolver::copyToPETScVec(Vec& petsc_x,
        SAMRAIVectorReal<NDIM, double>& x,
        Pointer<PatchLevel<NDIM> > patch_level)
{
    const int u_idx = x.getComponentDescriptorIndex(0);
    const int p_idx = x.getComponentDescriptorIndex(1);
    StaggeredStokesPETScVecUtilities::copyToPatchLevelVec(
        petsc_x, u_idx, d_u_dof_index_idx, p_idx, d_p_dof_index_idx, patch_level);
    return;
} // copyToPETScVec