Hello Sparsh,
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,
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 <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
--
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+unsubscr...@googlegroups.com
<mailto:dealii+unsubscr...@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>.
--
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+unsubscr...@googlegroups.com.
To view this discussion visit
https://groups.google.com/d/msgid/dealii/5bb10273-532f-4d91-b750-ce8e01c585c6%40gmail.com.
