hi, Jean-Paul,
here i attach a minimal test solving a Possion equations with trilinos's
amesos's direct solver, the memory
increases with the time step. Could you please check if there is anything
wired in that? Thanks in advance,
best,
lailai
On Tuesday, 31 January 2017 01:42:06 UTC-5, Jean-Paul Pelteret wrote:
>
> Dear Lailai,
>
> Internally we store the Amesos solver in a smart pointer
> <https://github.com/dealii/dealii/blob/master/include/deal.II/lac/trilinos_solver.h#L737>,
>
> so when the solver goes out of scope the solver is destroyed. At this time,
> the Amesos solver should clean up after itself - in its class
> documentation
> <https://trilinos.org/docs/dev/packages/amesos/doc/html/classAmesos__BaseSolver.html>
>
> I can find no mention of a "delete" or "destroy" function. So I don't think
> its likely that this is causing a memory leak, but if you're able to
> produce a minimal test case then we could investigate further. I suppose I
> should ask, which version of Trilinos are you using?
>
> Regards,
> Jean-Paul
>
>
> On Tuesday, January 31, 2017 at 3:41:47 AM UTC+1, femFluid wrote:
>>
>> hi, dear all,
>>
>> i am using trilino's direct solver interface to solve a time-dependent
>> problem. I realized that
>> when the number of iterations increases, the memory usage also increases.
>> This increase
>> does not arise when turning the solver off, so just assembling the
>> matrices.
>>
>> At the beginning, i used the following subroutine at every time step,
>>
>> {
>> deallog.push("DirectKLU");
>> TrilinosWrappers::SolverDirect::AdditionalData data;
>> data.solver_type = "Amesos_Klu";
>> SolverControl solver_control (1000, 1e-10);
>> TrilinosWrappers::SolverDirect solver(solver_control, data);
>> solver.solve (matrix, solution, rhs);
>> thick_constraints.distribute (thick_solution);
>> }
>>
>> then i realized that the memory might be allocated every time when
>> defining TrilinosWrappers::SolverDirect solver(.....), which is not
>> released.
>> Then I predefine the solver in a namespace, which is the used by calling
>> solver.solve (matrix, solution, rhs) at each step. Still, the memory
>> leakage
>> occurs.
>>
>> I also tried to use "Amesos_Mumps", the problem still exists. I tried to
>> see the documentation of Amesos, i realized that when amesos solver
>> finished solving, there is a procedure 'delete solver' after. May I ask
>> whether
>> such a member function exists in dealII wrappers? Thanks in advance,
>>
>> best,
>> lailai
>>
>
--
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].
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.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/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/numerics/solution_transfer.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_block_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_block_vector.h>
#include <deal.II/lac/trilinos_precondition.h>
#include <deal.II/lac/trilinos_solver.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <limits>
#include <string>
namespace Solver_
{
using namespace dealii;
TrilinosWrappers::SolverDirect::AdditionalData data_klu ("Amesos_Klu");
SolverControl solver_control_klu (1000, 1e-10);
TrilinosWrappers::SolverDirect solver_klu(solver_control_klu, data_klu);
}
namespace UserGeometry // currently only for 2d simulations.
{
const double x_inlet = -1;
const double x_outlet = 1;
const double y_bottom = -1;
const double y_top = 1;
}
namespace Name_Space
{
using namespace dealii;
template <int dim>
class Problem
{
public:
Problem ();
void run ();
private:
void setup_dofs ();
void assemble_p_system ();
void output_results () const;
void solve_p ();
Triangulation<dim> triangulation;
double global_Omega_diameter;
const unsigned int p_degree;
FE_Q<dim> p_fe;
DoFHandler<dim> p_dof_handler;
ConstraintMatrix p_constraints;
TrilinosWrappers::SparseMatrix p_matrix;
TrilinosWrappers::MPI::Vector p_rhs;
TrilinosWrappers::MPI::Vector p_solution,
old_p_solution;
double time_step;
double old_time_step;
unsigned int timestep_number;
bool rebuild_p_matrix;
};
template <int dim>
Problem<dim>::Problem ()
:
triangulation (Triangulation<dim>::maximum_smoothing),
p_degree (1),
p_fe (p_degree),
p_dof_handler (triangulation),
time_step (0),
old_time_step (0),
timestep_number (0),
rebuild_p_matrix (true)
{}
template <int dim>
void Problem<dim>::setup_dofs ()
{
{
p_dof_handler.distribute_dofs (p_fe);
p_constraints.clear ();
DoFTools::make_hanging_node_constraints (p_dof_handler,
p_constraints);
VectorTools::interpolate_boundary_values (p_dof_handler,
1,
ConstantFunction<2>(1.),
p_constraints);
p_constraints.close ();
}
const unsigned int
n_p = p_dof_handler.n_dofs();
{
p_matrix.clear ();
DynamicSparsityPattern dsp (n_p, n_p);
DoFTools::make_sparsity_pattern (p_dof_handler, dsp,
p_constraints, false);
p_matrix.reinit (dsp);
}
IndexSet p_partitioning = complete_index_set (n_p);
p_solution.reinit (p_partitioning, MPI_COMM_WORLD);
old_p_solution.reinit (p_partitioning, MPI_COMM_WORLD);
p_rhs.reinit (p_partitioning, MPI_COMM_WORLD);
p_solution = 0;
old_p_solution = 0;
}
template <int dim>
void Problem<dim>::assemble_p_system ()
{
if (rebuild_p_matrix == true)
{
p_matrix = 0;
}
p_rhs = 0;
QGauss<2> quadrature_formula(2);
FEValues<2> fe_values (p_fe, quadrature_formula,
update_values | update_gradients | update_JxW_values);
const unsigned int dofs_per_cell = p_fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
DoFHandler<2>::active_cell_iterator
cell = p_dof_handler.begin_active(),
endc = p_dof_handler.end();
for (; cell!=endc; ++cell)
{
fe_values.reinit (cell);
cell_matrix = 0;
cell_rhs = 0;
for (unsigned int q_index=0; q_index<n_q_points; ++q_index)
{
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_grad (i, q_index) *
fe_values.shape_grad (j, q_index) *
fe_values.JxW (q_index));
for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_rhs(i) += (fe_values.shape_value (i, q_index) *
1 *
fe_values.JxW (q_index));
}
cell->get_dof_indices (local_dof_indices);
if (rebuild_p_matrix == true)
p_constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
p_matrix,
p_rhs);
else
p_constraints.distribute_local_to_global (cell_rhs,
local_dof_indices,
p_rhs);
}
rebuild_p_matrix = false;
}
template <int dim>
void Problem<dim>::output_results () const
{
if (timestep_number % 1000 != 0)
return;
DataOut<dim> data_out;
data_out.add_data_vector (p_dof_handler, p_solution, "p");
data_out.build_patches (p_degree);
std::ostringstream filename;
filename << "solution-" << Utilities::int_to_string(timestep_number, 4) << ".vtk";
std::ofstream output (filename.str().c_str());
data_out.write_vtk (output);
}
// template <int dim>
// void SeaIceRheologyProblem<dim>::solve_stokes ()
// {
// std::cout << " Solving Stokes U and V..." << std::endl;
// {
// deallog.push("DirectKLU");
// Solver_::solver_klu.solve (stokes_matrix, stokes_solution, stokes_rhs);
// stokes_constraints.distribute (stokes_solution);
// }
//
// }
template <int dim>
void Problem<dim>::solve_p ()
{
std::cout << " Solving Pressure p..." << std::endl;
std::cout << " Assembling For p..." << std::endl << std::flush;
assemble_p_system ();
//std::cout << " Solving Stokes U and V..." << std::endl;
{
deallog.push("DirectKLU");
Solver_::solver_klu.solve (p_matrix, p_solution, p_rhs);
// p_constraints.distribute (p_solution);
}
}
template <int dim>
void Problem<dim>::run ()
{
{
std::vector<unsigned int> subdivisions (dim, 1);
subdivisions[0] = 1;
const Point<dim> bottom_left = (dim == 2 ?
Point<dim>(UserGeometry::x_inlet,UserGeometry::y_bottom) :
Point<dim>(-2,0,-1));
const Point<dim> top_right = (dim == 2 ?
Point<dim>(UserGeometry::x_outlet,UserGeometry::y_top) :
Point<dim>(-1,1,0));
GridGenerator::subdivided_hyper_rectangle (triangulation,
subdivisions,
bottom_left,
top_right);
}
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end(); ++cell)
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
if (cell->face(f)->center()[0] == UserGeometry::x_inlet)
cell->face(f)->set_all_boundary_ids(1);
if (cell->face(f)->center()[dim-1] == UserGeometry::y_top ||
cell->face(f)->center()[dim-1] == UserGeometry::y_bottom)
cell->face(f)->set_all_boundary_ids(2);
if (cell->face(f)->center()[0] == UserGeometry::x_outlet)
cell->face(f)->set_all_boundary_ids(3);
}
triangulation.refine_global (4);
setup_dofs ();
for (unsigned int timestep = 0; timestep<100000;
++timestep)
{
solve_p ();
timestep_number++;
output_results ();
std::cout<<timestep_number<<" "<<std::endl;
}
}
}
int main (int argc, char *argv[])
{
try
{
using namespace dealii;
using namespace Name_Space;
Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv,
numbers::invalid_unsigned_int);
AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)==1,
ExcMessage("This program can only be run in serial, use ./step-31"));
Problem<2> flow_problem;
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;
}
