Hi again, I took a look at the modified mwe you sent -- immediately noticing that you read though the whole thing, which is amazing, considering that this is a google group! So, thanks again for your time and effort, it is really appreciated!
I have been busy doing the following with our little toy problem: (i) Upgraded to deal.ii v9.2.0 with candi, together with PETSc 3.13.1 (and Trilinos 12.18.1). (ii) Further stripped down the mwe and deleted unnecessary parts and moved all but the mmult() out of the 'main loop'. So, I still observe the following: using Trilinos, all is fine. and using PETSc, using the call leftMatrix.mmult(result, rightMatrix) seems fine, whereas calling leftMatrix.mmult(result, rightMatrix, vector) still results in a memory leak for me. Could you verify that again for me please? If this is doing fine for you I will just multiply the rows of rightMatrix 'by hand' before using the call without the vector. Thanks in advance & kind regards, Richard -- 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 [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/030ba907-3f33-4831-9e7d-9ff018db5c5eo%40googlegroups.com.
/* This code is licensed under the "GNU GPL version 2 or later". See license.txt or https://www.gnu.org/licenses/gpl-2.0.html */ // Include files #include <deal.II/base/conditional_ostream.h> #include <deal.II/base/function.h> #include <deal.II/base/index_set.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/dofs/dof_renumbering.h> #include <deal.II/dofs/dof_tools.h> #include <deal.II/fe/fe_q.h> #include <deal.II/fe/fe_system.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/lac/block_sparse_matrix.h> #include <deal.II/lac/generic_linear_algebra.h> #include <deal.II/lac/sparsity_tools.h> #include <deal.II/numerics/vector_tools.h> #include <deal.II/distributed/tria.h> namespace LA { using namespace dealii::LinearAlgebraPETSc; // using namespace dealii::LinearAlgebraTrilinos; } // deal.II using namespace dealii; // Define flow problem. template <int dim> class flow_problem { public: flow_problem(const FESystem<dim> &fe); ~flow_problem(); void run(); private: // Create system matrix, rhs and distribute degrees of freedom. void setup_system(); const FESystem<dim> &fe; MPI_Comm mpi_communicator; parallel::distributed::Triangulation<dim> triangulation; DoFHandler<dim> dof_handler; ConditionalOStream pcout; // Distributed vectors. LA::MPI::BlockSparseMatrix system_matrix; LA::MPI::SparseMatrix B_inv_diag_D_C; AffineConstraints<double> constraints; std::vector<IndexSet> own_partitioning, rel_partitioning; }; // We are going to use the finite element discretization: // Q_(k+1)^c for solid and fluid displacements and velocities // and Q_k^dc for the pressure. template <int dim> flow_problem<dim>::flow_problem(const FESystem<dim> &fe) : fe(fe) , mpi_communicator(MPI_COMM_WORLD) , triangulation(mpi_communicator, typename Triangulation<dim>::MeshSmoothing( Triangulation<dim>::smoothing_on_refinement | Triangulation<dim>::smoothing_on_coarsening)) , dof_handler(triangulation) , pcout(std::cout, (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) {} // Destructor. template <int dim> flow_problem<dim>::~flow_problem() {} // System setup template <int dim> void flow_problem<dim>::setup_system() { // Get some mesh. GridGenerator::hyper_cube(triangulation, -0.5, 0.5, false); triangulation.refine_global(4); // Clear existing matrices. system_matrix.clear(); B_inv_diag_D_C.clear(); // Distribute DoFs. dof_handler.distribute_dofs(fe); // Cuthill_McKee renumbering. DoFRenumbering::Cuthill_McKee(dof_handler); // Component numbers. std::vector<unsigned int> block_component(dim + 1, 0); // velocity block_component[dim] = 1; // pressure // Count DoFs per block. std::vector<unsigned int> dofs_per_block(2); DoFTools::count_dofs_per_block(dof_handler, dofs_per_block, block_component); const unsigned int n_v = dofs_per_block[0], n_p = dofs_per_block[1]; // Renumber by block components. DoFRenumbering::component_wise(dof_handler, block_component); // Indexsets. IndexSet own_dofs, rel_dofs; own_dofs = dof_handler.locally_owned_dofs(); DoFTools::extract_locally_relevant_dofs(dof_handler, rel_dofs); // Block system; velocities and pressure in seperate blocks. own_partitioning.resize(2); own_partitioning[0] = own_dofs.get_view(0, n_v); own_partitioning[1] = own_dofs.get_view(n_v, n_v + n_p); rel_partitioning.resize(2); rel_partitioning[0] = rel_dofs.get_view(0, n_v); rel_partitioning[1] = rel_dofs.get_view(n_v, n_v + n_p); // Get zero Dirichlet constraints on vector space. { constraints.clear(); constraints.reinit(rel_dofs); DoFTools::make_hanging_node_constraints(dof_handler, constraints); std::vector<bool> component_mask(dim + 1, true); component_mask[dim] = false; unsigned int dirichlet_boundary_id = 0; VectorTools::interpolate_boundary_values(dof_handler, dirichlet_boundary_id, Functions::ZeroFunction<dim>(dim + 1), constraints, component_mask); constraints.close(); } // Reinit system_matrix. const unsigned int n_blocks = 2; { BlockDynamicSparsityPattern bdsp(n_blocks, n_blocks); std::vector<unsigned int> n_dof_per_block(n_blocks); n_dof_per_block[0] = n_v; n_dof_per_block[1] = n_p; for (unsigned int i = 0; i < n_blocks; i++) { for (unsigned int j = 0; j < n_blocks; j++) { bdsp.block(i, j).reinit(n_dof_per_block[i], n_dof_per_block[j]); } } bdsp.compress(); bdsp.collect_sizes(); DoFTools::make_sparsity_pattern(dof_handler, bdsp, constraints, true /*keep constraints for mmult*/, Utilities::MPI::this_mpi_process( mpi_communicator)); SparsityTools::distribute_sparsity_pattern( bdsp, dof_handler.locally_owned_dofs(), mpi_communicator, rel_dofs); system_matrix.reinit(own_partitioning, bdsp, mpi_communicator); } } // Run. template <int dim> void flow_problem<dim>::run() { // Setup matrices. setup_system(); // Fill the blocks with some dummy values. system_matrix = 0; const unsigned int dofs_per_cell = fe.dofs_per_cell; std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell); typename DoFHandler<dim>::active_cell_iterator cell = dof_handler .begin_active(), endc = dof_handler.end(); for (; cell != endc; ++cell) { if (cell->is_locally_owned()) { cell->get_dof_indices(local_dof_indices); for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { system_matrix.add(local_dof_indices[i], local_dof_indices[j], 1.0); } } } } system_matrix.compress(VectorOperation::add); // Get the 'diagonal' of the D block. LA::MPI::Vector inv_diag_D; inv_diag_D.reinit(own_partitioning[1],mpi_communicator); inv_diag_D = 0; inv_diag_D.add(1.0); // Outer loop to repeatedly perform the mmult(). const unsigned int max_k = 1000; for (unsigned int k = 0; k < max_k; k++) { pcout << "k = " << k << " / " << max_k << std::endl; // This causes a memory leak: system_matrix.block(0, 1).mmult(B_inv_diag_D_C, system_matrix.block(1, 0), inv_diag_D); // This does not cause a memory leak: system_matrix.block(0, 1).mmult(B_inv_diag_D_C, system_matrix.block(1, 0)); } } // Main function. int main(int argc, char *argv[]) { try { Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); deallog.depth_console(0); const unsigned int dim = 3, vel_degree = 1, press_degree = 1; FESystem<dim> fe(FE_Q<dim>(vel_degree), dim, FE_Q<dim>(press_degree), 1); flow_problem<dim> flow_problem(fe); flow_problem.run(); } 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; }
