Printing PETScWrappers::MPI::SparseMatrix to a single file
Sparsh Sinha
Hi everyone,
I am using dealii version 9.5.1. I have a basic question about printing ‘PETScWrappers::MPI::SparseMatrix system_matrix’ to a file.
I am trying to debug my program and found this (https://groups.google.com/g/dealii/c/BQjIt9XwwMI/m/xhUwz1yQAwAJ) solution to be useful. There are two solutions that have been provided in the link mentioned. First is to print by looping over all indices, and second is by using the 'print()' function. However, my program fails much later into the simulation and uses a large number of elements (DOF: 20K+). So I cannot use the first solution that leaves me with the second solution, which I have been able to implement into my program. But I run my program on multiple cores (16 or 32), which spreads system_matrix values to multiple files, since I would like to check system_matrix properties with other programs (e.g. MATLAB), so is it possible to print it into a single file. I also found ‘write_ascii()’ function to be useful but it prints on the screen, maybe it can be used in place of 'print()' if its output can be redirected to a file.
I have attached a minimial program below which implements both ‘print()’ and ‘write_ascii()’ function.
/////////////////////////////
using namespace dealii;
int main(int argc, char **argv)
{
try
{
using namespace dealii;
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
std::cout << "Running" << std::endl;
parallel::shared::Triangulation<2> triangulation(MPI_COMM_WORLD);
DoFHandler<2> dof_handler(triangulation);
FE_Q<2> fe_obj(1);
AffineConstraints<double> hanging_node_constraints;
const QGauss<2> quadrature_formula(fe_obj.degree + 1);
PETScWrappers::MPI::SparseMatrix system_matrix;
MPI_Comm mpi_communicator(MPI_COMM_WORLD);
const unsigned int n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator));
const unsigned int this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator));
ConditionalOStream pcout(std::cout, this_mpi_process == 0);
IndexSet locally_owned_dofs;
IndexSet locally_relevant_dofs;
std::ofstream matrixfile;
//Parallel processes can write in separate files using print function
const std::string matrixfilename = ("debug_matrix" +
Utilities::int_to_string(this_mpi_process, 3) +
".txt");
// This does not produce proper results, it overwrites
/*const std::string matrixfilename = ("debug_matrix.txt");*/
matrixfile.open(matrixfilename);
/////// Grid creation
Point<2> corner1, corner2;
std::vector<unsigned int> repetitions;
corner1 = Point<2>(0, 0);
corner2 = Point<2>(10, 10);
repetitions.push_back(10);
repetitions.push_back(10);
GridGenerator::subdivided_hyper_rectangle(triangulation,
repetitions ,
corner1 ,
corner2 ,
false );
///////
/////// Setup system
dof_handler.distribute_dofs(fe_obj);
locally_owned_dofs = dof_handler.locally_owned_dofs();
locally_relevant_dofs =
DoFTools::extract_locally_relevant_dofs(dof_handler);
hanging_node_constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler,
hanging_node_constraints);
DynamicSparsityPattern sparsity_pattern(locally_relevant_dofs);
DoFTools::make_sparsity_pattern(dof_handler,
sparsity_pattern,
hanging_node_constraints,
/*keep constrained dofs*/ false);
SparsityTools::distribute_sparsity_pattern(sparsity_pattern,
locally_owned_dofs,
mpi_communicator,
locally_relevant_dofs);
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
sparsity_pattern,
mpi_communicator);
///////
/////// Assemble system
system_matrix = 0;
FEValues<2> fe_values(fe_obj,
quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe_obj.n_dofs_per_cell();
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
if (cell->is_locally_owned())
{
cell_matrix = 0;
fe_values.reinit(cell);
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
cell_matrix(i, j) +=
fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point) *
fe_values.JxW(q_point); //
}
}
cell->get_dof_indices(local_dof_indices);
hanging_node_constraints.distribute_local_to_global(cell_matrix,
local_dof_indices,
system_matrix);
}
system_matrix.compress(VectorOperation::add);
/////// Write matrix
matrixfile << "\n System matrix values:" <<std::endl;
///// CORE QUERY
system_matrix.write_ascii(PETSC_VIEWER_ASCII_MATLAB);
system_matrix.print(matrixfile);
///////
matrixfile.close();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
//////////////////////
Is there any way to print the system_matrix to a single file either by using 'print()' or ‘write_ascii()’ function? Or do I need to use some other method for printing to a single file?
Best
Sparsh
Yann Jobic
I use the function "block_write" for sequential matrix/vectors, which
writes a binary file. It's way better than text files.
I used the function this way :
std::ofstream out5("MatrixSystem.bin");
system_matrix.block_write(out5);
std::ofstream out6("RHS.bin");
system_rhs.block_write(out6);
You may also want to look at the matlab functions to load the files.
https://github.com/dealii/dealii/wiki/Interfacing-With-Matlab
The matlab code is a bit buggy, but it helped me. I can share my scripts
if you want.
I don't know if it's working in parallel with the petsc types.
It may not be the best solution !
Regards,
Yann
Le 11/5/2024 à 5:43 AM, Sparsh Sinha a écrit :
> Hi everyone,
>
> printing ‘*PETScWrappers::MPI::SparseMatrix system_matrix*’ to a file.
> I am trying to debug my program and found this (https://
> groups.google.com/g/dealii/c/BQjIt9XwwMI/m/xhUwz1yQAwAJ>) solution to be
> mentioned. First is to print by looping over all indices, and second is
> I cannot use the first solution that leaves me with the second solution,
> which I have been able to implement into my program. But I run my
> program on multiple cores (16 or 32), which spreads system_matrix values
> to multiple files, since I would like to check system_matrix properties
> with other programs (e.g. MATLAB), so is it possible to print it into a
> other method for printing to a single file?
>
> Best
>
> Sparsh
>
> --
> The deal.II project is located at http://www.dealii.org/ <http://
> www.dealii.org/>
> For mailing list/forum options, see https://groups.google.com/d/forum/
> dealii?hl=en <https://groups.google.com/d/forum/dealii?hl=en>
> ---
> You received this message because you are subscribed to the Google
> Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion visit https://groups.google.com/d/msgid/
> dealii/23f65cfd-0d45-451f-8754-eb1e30e703fdn%40googlegroups.com
> <https://groups.google.com/d/msgid/dealii/23f65cfd-0d45-451f-8754-
> eb1e30e703fdn%40googlegroups.com?utm_medium=email&utm_source=footer>.
Wolfgang Bangerth
On 11/4/24 21:43, Sparsh Sinha wrote:
>
> *Is there any way to print the system_matrix to a single file either by
I would suggest to let every process print into its own file (say,
matrix-xxx.txt, where xxx is
std::to_string(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD))
). Then, after you run is finished, do
cat matrix-*.txt > matrix.txt
to collate all of these files into one.
Best
WB
Sparsh Sinha
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dealii/37c3fb95-48cc-4bc4-b584-548d00cfb9bf%40colostate.edu.