Problem solving transient brick with MUMPS

Jack · January 19, 2025, 6:23pm

I want to solve a two dimensional unsteady Navier-Stokes equation about circular cylinder flow. Belowing is my code:

#define GMM_WARNING_LEVEL 1

#include "getfem/getfem_model_solvers.h"
#include "getfem/getfem_export.h"
#include "getfem/getfem_mesher.h"
#include "getfem/getfem_generic_assembly.h"
#include "getfem/getfem_import.h"
#include "getfem/getfem_derivatives.h"
#include "getfem/getfem_assembling.h"

#include "gmm/gmm.h"
#include "gmm/gmm_MUMPS_interface.h"
#include "mumps_seq/mpi.h"

#include <filesystem>
#include <complex.h>
#include <chrono>

using std::endl; using std::cout; using std::cerr;
using std::ends; using std::cin; using std::string;

using bgeot::size_type;
using bgeot::base_node;
using bgeot::base_small_vector;
using bgeot::scalar_type; /* = double */

typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
typedef gmm::rsvector<scalar_type> sparse_vector_type;
typedef gmm::row_matrix<sparse_vector_type> sparse_matrix_type;
typedef std::vector<scalar_type> plain_vector;

typedef std::complex<double> complex;
typedef getfem::modeling_standard_complex_plain_vector complex_plain_vector;
typedef gmm::dense_matrix<std::complex<double>> complex_dense_matrix;
typedef getfem::modeling_standard_complex_sparse_matrix complex_sparse_matrix;

typedef gmm::mumps_interf<double>::MUMPS_STRUC_C MUMPS_STRUC;

/* Global parameters----------------------------------------- */
double Re = 50.;    /* Reynolds number */
double nu = 1./Re;  /* viscosity */
double dt = 0.05;   /* time interval */
int nt = 1640;      /* number of time iteration steps */
int ntSave = 20;    /* saving result */
double residual = 1e-8;

/* Main program---------------------------------------------- */
int main(int argc, char *argv[]){
    MPI_Init(&argc, &argv);
try{
    /* Initialise------------------------------- */
    getfem::mesh msh;
    getfem::import_mesh("./mesh/cylFlow.msh2", "gmsh", msh);
    int N = msh.dim();
    GMM_ASSERT1(N==2, "mesh dim != 2.");

    int MESH_INLET = 101;
    int MESH_OUTLET = 102;
    int MESH_TOP = 103;
    int MESH_BOTTOM = 104;
    int MESH_CYL_WALL = 105;

    getfem::mesh_fem mf_p(msh);
    getfem::mesh_fem mf_u(msh);
    getfem::mesh_fem mf_v(msh);
    getfem::mesh_im mim(msh);

    mf_p.set_qdim(1);
    mf_u.set_qdim(1);
    mf_v.set_qdim(1);

    getfem::pintegration_method ppi = getfem::int_method_descriptor("IM_TRIANGLE(3)");
    mim.set_integration_method(msh.convex_index(), ppi);

    getfem::pfem pf_p = getfem::fem_descriptor("FEM_PK(2,1)");
    getfem::pfem pf_u = getfem::fem_descriptor("FEM_PK(2,2)");
    getfem::pfem pf_v = getfem::fem_descriptor("FEM_PK(2,2)");

    mf_p.set_finite_element(msh.convex_index(), pf_p);
    mf_u.set_finite_element(msh.convex_index(), pf_u);
    mf_v.set_finite_element(msh.convex_index(), pf_v);

    string outPath = "./out/";
    if (!std::filesystem::exists(outPath)){
        std::filesystem::create_directory(outPath);
        cout << "Create output path: "+outPath << endl;
    }
    std::ofstream of1(outPath+"info.txt");

    cout << "\ngetfem mesh regions: " << msh.regions_index() << endl;
    of1 << "\ngetfem mesh regions: " << msh.regions_index() << endl;

    /* Build model---------------------------- */
    getfem::model model;
    model.add_fem_variable("p", mf_p);
    model.add_fem_variable("u", mf_u, 2);   /* niter=2 */
    model.add_fem_variable("v", mf_v, 2);   /* niter=2 */
    dal::bit_vector transientBricks;

    int nbdof_p = mf_p.nb_dof();
    int nbdof_u = mf_u.nb_dof();
    int nbdof_v = mf_v.nb_dof();
    int nbdof_total = nbdof_p + nbdof_u + nbdof_v;

    gmm::sub_interval I_p = model.interval_of_variable("p");
    gmm::sub_interval I_u = model.interval_of_variable("u");
    gmm::sub_interval I_v = model.interval_of_variable("v");

    cout << "\nnbdof: \n"
        << "nbdof_p: " << nbdof_p << endl
        << "nbdof_u: " << nbdof_u << endl
        << "nbdof_v: " << nbdof_v << endl
        << "\nsub interval: \n"
        << "I_p: " << I_p << endl
        << "I_u: " << I_u << endl
        << "I_v: " << I_v << endl;
    of1 << "\nnbdof: \n"
        << "nbdof_p: " << nbdof_p << endl
        << "nbdof_u: " << nbdof_u << endl
        << "nbdof_v: " << nbdof_v << endl
        << "\nsub interval: \n"
        << "I_p: " << I_p << endl
        << "I_u: " << I_u << endl
        << "I_v: " << I_v << endl;

    /* Nonlinear term: (u.nabla)u */
    int ibNonlinear = getfem::add_nonlinear_term(model, mim, "Test_u*[u,v].Grad_u + Test_v*[u,v].Grad_v");
    transientBricks.add(ibNonlinear);

    /* Laplacian term: -1/Re Delta u */
    model.add_initialized_scalar_data("nu", nu);
    int ibLaplacian = getfem::add_linear_term(model, mim, "nu*(Grad_u.Grad_Test_u + Grad_v.Grad_Test_v)");
    transientBricks.add(ibLaplacian);

    /* Pressure term: nabla p */
    int ibPressure = getfem::add_source_term(model, mim, "p*(Grad_Test_u(1)+Grad_Test_v(2))");
    transientBricks.add(ibPressure);

    /* Continuity condition */
    int ibContinuity = getfem::add_linear_term(model, mim, "Test_p*(Grad_u(1)+Grad_v(2))");
    transientBricks.add(ibContinuity);

    /* Boundary conditions----------------------- */
    /* Inlet: u=1, v=0 */
    model.add_initialized_scalar_data("Uin", 1.);
    getfem::add_Dirichlet_condition_with_simplification(model, "u", MESH_INLET, "Uin");
    getfem::add_Dirichlet_condition_with_simplification(model, "v", MESH_INLET);

    /* Outlet: p=0 */
    getfem::add_Dirichlet_condition_with_simplification(model, "p", MESH_OUTLET);

    /* Top and Bottom: v=0 */
    getfem::add_Dirichlet_condition_with_simplification(model, "v", MESH_TOP);
    getfem::add_Dirichlet_condition_with_simplification(model, "v", MESH_BOTTOM);

    /* Cylinder wall: u=0, v=0 */
    getfem::add_Dirichlet_condition_with_simplification(model, "u", MESH_CYL_WALL);
    getfem::add_Dirichlet_condition_with_simplification(model, "v", MESH_CYL_WALL);

    /* Transient part-------------------------------- */
    model.add_initialized_scalar_data("dt", dt);
    model.add_initialized_scalar_data("alpha", 0.5);
    getfem::add_basic_d_on_dt_brick(model, mim, "u", "dt", "alpha");
    getfem::add_basic_d_on_dt_brick(model, mim, "v", "dt", "alpha");
    getfem::add_midpoint_dispatcher(model, transientBricks);

    /* Output info */
    of1 << "=========================================================\n";
    model.listbricks(of1);
    of1 << "=========================================================\n";
    model.listvar(of1);

    /* Solve----------------------------------------------- */
    // MUMPS_STRUC id;
    gmm::mumps_interf<double>::MUMPS_STRUC_C id;
    id.job = -1;    /* MUMPS initialize */
    id.par = 1;
    id.sym = 0;     /* symmetric matrix=false */
    gmm::mumps_interf<double>::mumps_c(id);

    /* set MUMPS parameters, see https://petsc.org/release/manualpages/Mat/MATSOLVERMUMPS/ */
    id.icntl[0] = -1;
    id.icntl[1] = -1;
    id.icntl[2] = -1;
    id.icntl[3] = 0;
    id.icntl[13] = 200;     /* increase ICNTL(14) */
    id.job = 6;             /* Request all phases: analysis + factorization + solve */
    gmm::mumps_interf<double>::mumps_c(id);
    
    gmm::iteration iter(residual, 0, 5e4);
    model.first_iter();
    /* Initial conditions */
    plain_vector U0, V0;
    gmm::resize(U0, mf_u.nb_dof());
    gmm::resize(V0, mf_v.nb_dof());
    gmm::fill(U0, 1.);
    gmm::clear(V0);

    gmm::copy(U0, model.set_real_variable("u", 1));
    gmm::copy(V0, model.set_real_variable("v", 1));
    
    /* Time iteration */
    double sol_t = 0.;
    plain_vector X(nbdof_total);    /* solution of each time step */
    plain_vector rhs(nbdof_total);
    sparse_matrix TGM(nbdof_total, nbdof_total);

    for (int i=0; i<nt; ++i){
        sol_t += dt;
        cout << "t = " << sol_t << endl;
        iter.init();
        model.assembly(getfem::model::BUILD_ALL);
        getfem::standard_solve(model, iter);
        // gmm::copy(model.real_tangent_matrix(), TGM);
        // gmm::copy(model.real_rhs(), rhs);
        // gmm::MUMPS_solve(TGM, X, rhs);
        // gmm::copy(gmm::sub_vector(X, I_p), model.set_real_variable("p"));
        // gmm::copy(gmm::sub_vector(X, I_u), model.set_real_variable("u"));
        // gmm::copy(gmm::sub_vector(X, I_v), model.set_real_variable("v"));

        /* Export */
        if (i%ntSave==0){
            char s[100];
            sprintf(s, "step%05d", i);
            getfem::vtk_export plot(outPath+string(s)+".vtk", true);
            plot.exporting(msh);
            plot.write_point_data(mf_p, model.real_variable("p"), "p");
            plot.write_point_data(mf_u, model.real_variable("u"), "u");
            plot.write_point_data(mf_v, model.real_variable("v"), "v");
            cout << "\nExporting done." << endl;
        }

        model.next_iter();
    }

    cout.precision(16);
    cout << "L2 norm of u: " << getfem::asm_L2_norm(mim, mf_u, model.real_variable("u")) << endl
        << "Linfty norm of u: " << gmm::vect_norm1(model.real_variable("u")) << endl
        << "L2 norm of v: " << getfem::asm_L2_norm(mim, mf_v, model.real_variable("v")) << endl
        << "Linfty norm of v: " << gmm::vect_norm1(model.real_variable("v")) << endl
        << "L2 norm of p: " << getfem::asm_L2_norm(mim, mf_p, model.real_variable("p")) << endl
        << "Linfty norm of p: " << gmm::vect_norm1(model.real_variable("p")) << endl;
}
GMM_STANDARD_CATCH_ERROR;
MPI_Finalize();
return 0;
}

After compiling the code, I ran the program, but got an error message:

============================================
|    A GMM error has been detected !!!     |
============================================
Error in ./gmm/gmm_MUMPS_interface.h, line 153 : 
Solve with MUMPS failed: error -9, increase ICNTL(14)

The debug call stack generated by lldb is:

__cxa_throw (@__cxa_throw:3)
___lldb_unnamed_symbol4801 (@___lldb_unnamed_symbol4801:79)
___lldb_unnamed_symbol11886 (@___lldb_unnamed_symbol11886:235)
___lldb_unnamed_symbol8322 (@___lldb_unnamed_symbol8322:14)
getfem::standard_solve(getfem::model&, gmm::iteration&, std::shared_ptr<getfem::abstract_linear_solver<gmm::col_matrix<gmm::rsvector<double>>, std::vector<double, std::allocator<double>>>>, getfem::abstract_newton_line_search&) (@getfem::standard_solve(getfem::model&, gmm::iteration&, std::shared_ptr<getfem::abstract_linear_solver<gmm::col_matrix<gmm::rsvector<double>>, std::vector<double, std::allocator<double>>>>, getfem::abstract_newton_line_search&):454)
getfem::standard_solve(getfem::model&, gmm::iteration&) (@getfem::standard_solve(getfem::model&, gmm::iteration&):54)
main (\mnt\d\codes\CFD\getfem\test\myProject\cylinderFlow\Re50\unsteadyNavierStokes\mycode.cc:213)
___lldb_unnamed_symbol3278 (@___lldb_unnamed_symbol3278:29)
__libc_start_main (@__libc_start_main:44)
_start (@_start:15)

It seems that I need to modify the MUMPS solver parameters utilized by getfem::standard_solve. I tried to set some parameters of gmm::mumps_interf<double>::MUMPS_STRUC_C id;, but it didn’t work.

I used to solve a steady NS equation by using getfem and didn’t encounter this kind of error. Could anyone please tell me the reason and how to revise my codes?

Thanks a lot!

Konstantinos.Poulios · January 20, 2025, 1:05pm

if you have increased ICNTL(14) inside the file “gmm_MUMPS_interface.h” and is still giving the same error, then your tangent matrix might not be sparse enough. Can it be that the system tangent matrix contains some garbage? Or that you have basis function with big overlaps?

I would suggest to solve a minimal problem where you can print and inspect the tangent matrix.

Jack · January 21, 2025, 8:35pm

Thanks for yourreply! I will try it and give a response.

Topic		Replies	Views
Parallelization Development	4	124	June 24, 2024
Probleme installation MUMPS sur Mac ox General	5	40	January 24, 2025
Fast solve of linear problems General	5	18	February 3, 2025
Export im_nodes scalars on boundary regions Questions and issues	4	22	October 13, 2024
Any examples of XFEM in action? Questions and issues	2	19	November 26, 2024

Problem solving transient brick with MUMPS

Related topics