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
--
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/23f65cfd-0d45-451f-8754-eb1e30e703fdn%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h> // Quadrature
#include <deal.II/base/function.h> // Funtion class
#include <deal.II/base/logstream.h> // Logging dealii
#include <deal.II/base/tensor.h> // Tensor
#include <deal.II/base/symmetric_tensor.h> // Implements symmetric tensors of rank 2 and 4
// Vector and Matrix libraries
#include <deal.II/lac/vector.h> // Vector
#include <deal.II/lac/full_matrix.h> // Full Dense matrix, whether 0 or non 0
#include <deal.II/lac/sparse_matrix.h> // Sparse Matrix, only stored non 0
#include <deal.II/lac/dynamic_sparsity_pattern.h> // Dynamic matrix, grows, at end needs to be converted to Sparse Matrix
#include <deal.II/lac/solver_cg.h> // Conjugate Gradient solver
#include <deal.II/lac/precondition.h> // Preconditions, simplifies linear systems for solving
#include <deal.II/lac/affine_constraints.h> // Applies constraints on dofs boundary or refinement
// Parallel headers
#include <deal.II/base/multithread_info.h> // Useful in multithread programming
#include <deal.II/lac/petsc_vector.h> // Parallelize vector and operations
#include <deal.II/lac/petsc_sparse_matrix.h> // Partitions and parallelize matrices
#include <deal.II/lac/petsc_solver.h> // Parallelize solver
#include <deal.II/lac/petsc_precondition.h> // Again parallel precondition
#include <deal.II/lac/sparsity_tools.h> // Performs function on sparsity patterns, renumbering, etc
// Triangulation objects, generators, tools, iterators
#include <deal.II/grid/tria.h> // Triangulation object
#include <deal.II/grid/grid_generator.h> // Grid Generator
#include <deal.II/grid/grid_refinement.h> // Grid refinement and functions
#include <deal.II/grid/manifold_lib.h> // Manifold lib, Polar, Spherical, Flat, etc.
#include <deal.II/grid/grid_tools.h> // Grid manipulation tools
#include <deal.II/grid/tria_accessor.h> // Triagulation object accessor, e.g. vertices, child, etc
#include <deal.II/grid/tria_iterator.h> // Iterator to face, cells, vertices, etc.
#include <deal.II/grid/grid_out.h>
// And a header that implements filters for iterators looping over all
// cells. We will use this when selecting only those cells for output that are
// owned by the present process in a parallel program:
#include <deal.II/grid/filtered_iterator.h> // Filters interator using predicate, etc.
// DOF and tolls
#include <deal.II/dofs/dof_handler.h> // For a triangulation enumerates, and some functions
#include <deal.II/dofs/dof_accessor.h> // Accesses data for edges, faces and cells
#include <deal.II/dofs/dof_tools.h> // DOF manipulations tools, e.g. renumbering
#include <deal.II/dofs/dof_renumbering.h> // Renumbering algorithm
// FE and Q values and systems
#include <deal.II/fe/fe_values.h> // To assemble matrices, vectors. Link FE, Q, maps
#include <deal.II/fe/fe_system.h> // Vector values elements
#include <deal.II/fe/fe_q.h> // Implements scalar lagrange FE Q points
#include <deal.II/fe/fe_dgq.h> // Same thing but for discontinuous FE
#include <deal.II/fe/mapping_q.h> // Polynomial mappings Qp to get higher order mapping
// Tools, error and output
#include <deal.II/numerics/vector_tools.h> // Operations on vectors
#include <deal.II/numerics/matrix_tools.h> // Tools to apply on matrices, including BC
#include <deal.II/numerics/data_out.h> // Data out
#include <deal.II/numerics/error_estimator.h> // Error estimation per cell
// Solvers and others
#include <deal.II/lac/identity_matrix.h> // Identity matrix
#include <deal.II/lac/solver_minres.h> // Minimal residual method for symmetric matrices
#include <deal.II/lac/solver_gmres.h> // Generalized Minimal Residual Method, norm stops
#include <deal.II/lac/sparse_direct.h> // Sparse direct solver
// Parallel stuff
#include <deal.II/distributed/shared_tria.h> // Parallel triangulation every processor knows about every cell of the global mesh
#include <deal.II/base/conditional_ostream.h> // Conditional ostream pcout on mpi id
#include <deal.II/base/utilities.h> // Mainly MPI abstraction
// To compute rotation matrices of the local coordinate systems at specific
// points in the domain.
#include <deal.II/physics/transformations.h> // Assist in transformation of tensor quantities
// C++ headers
#include <fstream> // read/ write to file
#include <cmath> // cmath pow etc
#include <iostream> // cin cout cerr
#include <stdio.h> // printf scanf
#include <sstream> // string stream class
#include <unistd.h> // standard symbolic constants and types
#include <fcntl.h> // file control options
#include <cstdio> // C stdio
#include <iomanip> // manipulation i/o
#include "stdlib.h" // malloc free
#include <chrono> // time
#include <random> // random no.
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;
}
