Method initializeStokesSolver

src/IB/CIBSaddlePointSolver.cpp:633–777

} // initializeKSP

void
CIBSaddlePointSolver::initializeStokesSolver(const SAMRAIVectorReal<NDIM, double>& sol_vec,
                                             const SAMRAIVectorReal<NDIM, double>& rhs_vec)
{
    Pointer<PatchHierarchy<NDIM>> patch_hier = sol_vec.getPatchHierarchy();
    const int coarsest_ln = sol_vec.getCoarsestLevelNumber();
    const int finest_ln = sol_vec.getFinestLevelNumber();

    // Set the nullspace of the LInv and subdomain solvers
    const double rho = d_ins_integrator->getStokesSpecifications()->getRho();
    const bool has_velocity_nullspace = d_normalize_velocity && IBTK::abs_equal_eps(rho, 0.0);
    const bool has_pressure_nullspace = d_normalize_pressure;

    for (const auto& nul_vec : d_nul_vecs)
    {
        if (nul_vec) free_vector_components(*nul_vec);
    }
    const int n_nul_vecs = (has_pressure_nullspace ? 1 : 0) + (has_velocity_nullspace ? NDIM : 0);
    d_nul_vecs.resize(n_nul_vecs);

    for (const auto& U_nul_vec : d_U_nul_vecs)
    {
        if (U_nul_vec) free_vector_components(*U_nul_vec);
    }
    const int n_U_nul_vecs = (has_velocity_nullspace ? NDIM : 0);
    d_U_nul_vecs.resize(n_U_nul_vecs);

    if (has_velocity_nullspace)
    {
        for (unsigned int k = 0; k < NDIM; ++k)
        {
            d_nul_vecs[k] = sol_vec.cloneVector(d_object_name + "::nul_vec_U_" + std::to_string(k));
            d_nul_vecs[k]->allocateVectorData(d_current_time);
            d_nul_vecs[k]->setToScalar(0.0);

            SAMRAIVectorReal<NDIM, double> svr_u(
                d_object_name + "::U_nul_vec_U_" + std::to_string(k), patch_hier, coarsest_ln, finest_ln);
            svr_u.addComponent(sol_vec.getComponentVariable(0),
                               sol_vec.getComponentDescriptorIndex(0),
                               sol_vec.getControlVolumeIndex(0));

            d_U_nul_vecs[k] = svr_u.cloneVector(svr_u.getName());
            d_U_nul_vecs[k]->allocateVectorData(d_current_time);
            d_U_nul_vecs[k]->setToScalar(0.0);
            for (int ln = coarsest_ln; ln <= finest_ln; ++ln)
            {
                Pointer<PatchLevel<NDIM>> level = patch_hier->getPatchLevel(ln);
                for (PatchLevel<NDIM>::Iterator p(level); p; p++)
                {
                    Pointer<Patch<NDIM>> patch = level->getPatch(p());
                    Pointer<SideData<NDIM, double>> nul_data =
                        patch->getPatchData(d_nul_vecs[k]->getComponentDescriptorIndex(0));
                    nul_data->getArrayData(k).fillAll(1.0);
                    Pointer<SideData<NDIM, double>> U_nul_data =
                        patch->getPatchData(d_U_nul_vecs[k]->getComponentDescriptorIndex(0));
                    U_nul_data->getArrayData(k).fillAll(1.0);
                }
            }