Poisson equation solver not converging in C++, but solvable in MATLAB

21 views

Skip to first unread message

皮卡丘卡皮卡丘

unread,

Jun 18, 2025, 11:51:40 PMJun 18

to IBAMR Users

Hello everyone,

I am implementing a fiber orientation modeling strategy inspired by the paper "Image-based immersed boundary model of the aortic root", where the valve leaflet fiber direction is derived from the gradient of a scalar potential field u obtained by solving the Poisson equation.
I wrote the C++ code using libMesh (main parts shown below):

void apply_poisson_bc(libMesh::LinearImplicitSystem& system, const std::string& sys_name)
{
    using namespace libMesh;

    static ConstFunction<Number> const_func_one(60.0);
    static ConstFunction<Number> const_func_zero(0.0);

    if (sys_name == "Poisson_u")
    {
        
        // commissure 边界 ID=2，值=1
        std::set<boundary_id_type> commissure_bdry{3};
        system.get_dof_map().add_dirichlet_boundary(
            DirichletBoundary(commissure_bdry, {0}, &const_func_one)
        );
        
        // 自由边 边界 ID=1，值=0
        std::set<boundary_id_type> free_edge_bdry{4};
        system.get_dof_map().add_dirichlet_boundary(
            DirichletBoundary(free_edge_bdry, {0}, &const_func_zero)
        );
        
    }
    else if (sys_name == "Poisson_v")
    {
        // 上表面 边界 ID=3，值=0
        std::set<boundary_id_type> top_bdry{3};
        system.get_dof_map().add_dirichlet_boundary(
            DirichletBoundary(top_bdry, {0}, &const_func_zero)
        );

        // 下表面 边界 ID=4，值=1
        std::set<boundary_id_type> bottom_bdry{4};
        system.get_dof_map().add_dirichlet_boundary(
            DirichletBoundary(bottom_bdry, {0}, &const_func_one)
        );
    }
}


// 求解泊松方程
void solve_poisson(libMesh::EquationSystems& es, const std::string& sys_name)
{
    auto& system = es.get_system<libMesh::LinearImplicitSystem>(sys_name);
    const libMesh::MeshBase& mesh = es.get_mesh();
    const libMesh::DofMap& dof_map = system.get_dof_map();

    // 总自由度数目
    unsigned int total_dofs = dof_map.n_dofs();
    libMesh::out << "Total DOFs: " << total_dofs << std::endl;

    unsigned int constrained_dofs = dof_map.n_constrained_dofs();
    libMesh::out << "Total constrained DOFs: " << constrained_dofs << std::endl;

    system.matrix->zero();
    system.rhs->zero();

    int constrained_count = 0;  // 计数器

    const unsigned int dim = mesh.mesh_dimension();
    libMesh::FEType fe_type = system.variable_type(0);

    // 构建 FE 对象与积分规则
    auto fe = libMesh::FEBase::build(dim, fe_type);
    libMesh::QGauss qrule(dim, fe_type.order);
    fe->attach_quadrature_rule(&qrule);

    const std::vector<std::vector<libMesh::Real> >& phi = fe->get_phi(); // 形函数
    const std::vector<std::vector<libMesh::VectorValue<libMesh::Real> > >& dphi = fe->get_dphi(); // 形函数梯度

    const std::vector<libMesh::Real>& JxW = fe->get_JxW();

    // 装配矩阵
    for (auto el = mesh.active_local_elements_begin();
         el != mesh.active_local_elements_end(); ++el)
    {
        const libMesh::Elem* elem = *el;
        fe->reinit(elem);

        std::vector<libMesh::dof_id_type> dof_indices;
        dof_map.dof_indices(elem, dof_indices);

        const unsigned int n_dofs = dof_indices.size();

        libMesh::DenseMatrix<libMesh::Number> Ke(n_dofs, n_dofs);
        libMesh::DenseVector<libMesh::Number> Fe(n_dofs);

        for (unsigned int qp = 0; qp < qrule.n_points(); ++qp)
        {
            for (unsigned int i = 0; i < n_dofs; ++i)
            {
                for (unsigned int j = 0; j < n_dofs; ++j)
                {
                    Ke(i, j) += dphi[i][qp] * dphi[j][qp] * JxW[qp];
                }
                // 右端项若有 source term，可以加 source[i][qp]*JxW[qp]
                Fe(i) += 1e5 * phi[i][qp] * JxW[qp];
            }
        }
        */
        // 应用边界约束条件
        dof_map.heterogenously_constrain_element_matrix_and_vector(Ke, Fe, dof_indices);

        system.matrix->add_matrix(Ke, dof_indices);
        system.rhs->add_vector(Fe, dof_indices);    
    }
    
    system.matrix->close();
    system.rhs->close();
    system.matrix->print_matlab("K.mat");
    system.rhs->print_matlab("Fe.mat");
    // 打印全局 RHS 向量
    libMesh::out << "Global RHS: min=" << system.rhs->min()
              << ", max=" << system.rhs->max() << "\n";

    // 求解线性系统
    
    system.solve();
    system.update();


    // 打印求解器收敛情况和残差信息
    auto solver = system.get_linear_solver();

    if (solver)
    {
        int reason = solver->get_converged_reason();
        libMesh::out << "Converged reason = " << reason << "\n";
        libMesh::out << "Iterations = " << system.n_linear_iterations()
                    << ", Final residual = " << system.final_linear_residual() << "\n";
    }
}

and in the main code is :

        // 创建独立的 EquationSystems 对象，管理泊松方程
        std::unique_ptr<libMesh::EquationSystems> poisson_eq_sys = std::make_unique<libMesh::EquationSystems>(mesh1);

        // 定义有限元变量的顺序和族
        const libMesh::Order order = libMesh::FIRST;
        const libMesh::FEFamily family = libMesh::LAGRANGE;

        // 添加两个线性隐式系统：Poisson_u 和 Poisson_v
        auto& sys_u = poisson_eq_sys->add_system<libMesh::LinearImplicitSystem>("Poisson_u");
        auto& sys_v = poisson_eq_sys->add_system<libMesh::LinearImplicitSystem>("Poisson_v");

        // 添加变量
        sys_u.add_variable("u", order, family);
        sys_v.add_variable("v", order, family);    

        // 应用边界条件
        apply_poisson_bc(poisson_eq_sys->get_system<libMesh::LinearImplicitSystem>("Poisson_u"), "Poisson_u");
        apply_poisson_bc(poisson_eq_sys->get_system<libMesh::LinearImplicitSystem>("Poisson_v"), "Poisson_v");
        // 初始化方程系统（变量添加后）
        poisson_eq_sys->init();
        poisson_eq_sys->print_info();
        // 求解方程
        solve_poisson(*poisson_eq_sys, "Poisson_u");
        solve_poisson(*poisson_eq_sys, "Poisson_v");

        // 访问解向量
        auto& system_u = poisson_eq_sys->get_system<libMesh::LinearImplicitSystem>("Poisson_u");
        const auto& sol_u_ptr = system_u.current_local_solution;
    // 例如打印前10个解
        unsigned int vec_size = sol_u_ptr->size();
        for (unsigned int i = 0; i < std::min(vec_size, 10u); ++i)
        {
            std::cout << "sol_u[" << i << "] = " << (*sol_u_ptr)(i) << std::endl;
        }

Problem:

When calling system.solve(), the program hangs or stalls silently — no error, but also no residual or convergence output, as :
EquationSystems
n_systems()=2
System #0, "Poisson_u"
Type "LinearImplicit"
Variables="u"
Finite Element Types="LAGRANGE"
Approximation Orders="FIRST"
n_dofs()=2949
n_local_dofs()=2949
n_constrained_dofs()=244
n_local_constrained_dofs()=244
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 11.9447
Average Off-Processor Bandwidth <= 0
Maximum On-Processor Bandwidth <= 19
Maximum Off-Processor Bandwidth <= 0
DofMap Constraints
Number of DoF Constraints = 244
Number of Heterogenous Constraints= 122
Average DoF Constraint Length= 0
System #1, "Poisson_v"
Type "LinearImplicit"
Variables="v"
Finite Element Types="LAGRANGE"
Approximation Orders="FIRST"
n_dofs()=2949
n_local_dofs()=2949
n_constrained_dofs()=244
n_local_constrained_dofs()=244
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 11.9447
Average Off-Processor Bandwidth <= 0
Maximum On-Processor Bandwidth <= 19
Maximum Off-Processor Bandwidth <= 0
DofMap Constraints
Number of DoF Constraints = 244
Number of Heterogenous Constraints= 122
Average DoF Constraint Length= 0

Total DOFs: 2949
Total constrained DOFs: 244
Global RHS: min=0, max=133028
Converged reason = -11
Iterations = 0, Final residual = inf
Total DOFs: 2949
Total constrained DOFs: 244
Global RHS: min=0, max=133062
Converged reason = -11
Iterations = 0, Final residual = inf
sol_u[0] = inf
sol_u[1] = inf
sol_u[2] = inf
sol_u[3] = inf
sol_u[4] = inf
sol_u[5] = inf
sol_u[6] = inf
sol_u[7] = inf
sol_u[8] = inf
sol_u[9] = inf
I’ve tried using various solvers (CG, GMRES) and preconditioners (AMG, ILU), but nothing helps.

I exported the system matrix and RHS to MATLAB, and solving with phi = A\b gives the correct solution immediately. What could be the potential reasons why LinearImplicitSystem::solve() would hang or fail silently, especially when the same matrix and RHS solve fine in MATLAB? Is there anything I might be missing in the C++ setup or solver configuration?

I would really appreciate any suggestions or insights into what might be going wrong here. Thanks so much in advance for your help!

Best,
peng

Chang Wei

unread,

Jun 20, 2025, 2:45:33 AMJun 20

to IBAMR Users

Hi,

I did a quick test of your code and I think I have a rough idea of what’s going wrong. You shouldn't encapsulate the finite element assembly and solve steps entirely inside the solve_poisson function—this disrupts the expected workflow.

Generally, the workflow is: declare the system, assemble it, apply boundary conditions, initialize the system, and then solve. I suggest checking out the Poisson equation example provided by libMesh, which follows this standard computational workflow and should align well with your needs. You can find it at path_to_libMesh_src/examples/introduction/introduction_ex5

As for why MATLAB is able to compute the result, I suspect it's because the equation system in your code is not properly initialized, even if the assembly is done correctly. This is just a hypothesis, but it's likely that libMesh doesn’t solve the system using the raw assembled matrix directly—there may be some internal transformations or setup steps involved.

皮卡丘卡皮卡丘

unread,

Jun 20, 2025, 2:53:05 AMJun 20

to IBAMR Users

Thank you so much for your prompt and insightful response!

You're absolutely right, encapsulating the entire assembly and solve process inside the solve_poisson function did break the expected workflow. I refactored the code to follow the standard steps: declaring the system, assembling it, applying boundary conditions, initializing, and then solving, and it now works perfectly.

I also checked the example you mentioned (introduction_ex5), and it really helped clarify the intended structure. I appreciate the pointer.