Munmap_chunk(): invalid pointer

[Мария Бронзова]

Hello everyone, 

I am trying to run the following script, but getting the error: 'munmap_chunk(): invalid pointer'. The eclipse debugger doesn't help me realize, where exactly the error occurs. 

In case you have an idea, how to correct it, please, let me know. I would be very gratefull!

Kind regards,

Mariia 

#include <deal.II/base/quadrature_lib.h>

#include <deal.II/base/logstream.h>

#include <deal.II/base/function.h>

#include <deal.II/lac/block_vector.h>

#include <deal.II/lac/full_matrix.h>

#include <deal.II/lac/block_sparse_matrix.h>

#include <deal.II/lac/solver_cg.h>

#include <deal.II/lac/precondition.h>

#include <deal.II/lac/linear_operator.h>

#include <deal.II/lac/packaged_operation.h>

#include <deal.II/grid/tria.h>

#include <deal.II/grid/grid_generator.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_dgq.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/matrix_tools.h>

#include <deal.II/numerics/data_out.h>

#include <deal.II/fe/fe_q.h>

#include <deal.II/lac/sparse_direct.h>

//#include <deal.II/fe/fe_raviart_thomas.h>

#include <fstream>

#include <iostream>

#include <math.h>

#include <complex>

#include <deal.II/base/tensor_function.h>

namespace Poroelasticity

{

  using namespace dealii;

  template <int dim>

  SymmetricTensor<4, dim, std::complex<double>> get_stress_strain_tensor(const std::complex<double> lambda,

                       const std::complex<double> mu)

  {

  const std::complex<double> nothing(0.0,0.0);

  SymmetricTensor<4, dim, std::complex<double>> stress_strain_tensor;

  for (unsigned int i=0; i< dim; ++i)

  for (unsigned int j=0; j< dim; ++j)

  for (unsigned int k=0; k< dim; ++k)

  for (unsigned int l=0; l< dim; ++l)

  stress_strain_tensor[i][j][k][l]= (((i==k) && (j==l) ? mu : nothing) +

              ((i==l) && (j==k) ? mu : nothing) +

((i==j) && (k==l) ? lambda : nothing));

  return stress_strain_tensor;

  }

  template <int dim>

  inline SymmetricTensor<2, dim> get_strain(const FEValues<dim> &fe_values,

                  const unsigned int shape_func,

const unsigned int q_point)

  {

  SymmetricTensor<2, dim> strain_tensor;

  for (unsigned int i=0; i< dim; ++i)

  strain_tensor[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i];

  for (unsigned int i=0; i< dim; ++i)

  for (unsigned int j=i+1; j< dim; ++j)

  strain_tensor[i][j] =

  (fe_values.shape_grad_component(shape_func, q_point, i)[j] +

   fe_values.shape_grad_component(shape_func, q_point, j)[i]) /

   2;

  return strain_tensor;

  }

  template <int dim>

  class Poroelasticity

  {

  public:

    Poroelasticity(const unsigned int degree);

    void run();

  private:

    void make_grid_and_dofs();

    void assemble_system();

    void solve();

    void output_results() const;

    const unsigned int degree;

    Triangulation<dim> triangulation;

    FESystem<dim>      fe;

    DoFHandler<dim>    dof_handler;

    BlockSparsityPattern      sparsity_pattern;

    BlockSparseMatrix<std::complex<double>> system_matrix;

    BlockVector<std::complex<double>> solution;

    BlockVector<std::complex<double>> system_rhs;

    static const SymmetricTensor<4, dim> stress_strain_tensor;

static constexpr unsigned int start_freq=0;

static constexpr unsigned int end_freq=2; //[Hz]

  };

template<int dim>

  class RightHandSide : public Function<dim>

  {

  public:

  RightHandSide()

  : Function<dim>(1)

  {}

  virtual double value(const Point<dim> &p,

         const unsigned int component = 0) const override;

  

  };

  template <int dim>

  double RightHandSide<dim>::value(const Point<dim> & p,

               const unsigned int component) const

  {

  return 0;

  }

  template <int dim>

  Poroelasticity<dim>::Poroelasticity(const unsigned int degree)

    : degree(degree)

  , triangulation(Triangulation<dim>::maximum_smoothing)

    , fe(FE_Q<dim>(degree+1),dim, FE_Q<dim>(degree),1)

    , dof_handler(triangulation)

  {}

  template <int dim>

  void Poroelasticity<dim>::make_grid_and_dofs()

  {

    GridGenerator::hyper_cube(triangulation, -1, 1);

    triangulation.refine_global(1);

    dof_handler.distribute_dofs(fe);

    DoFRenumbering::Cuthill_McKee(dof_handler);

    std::vector<unsigned int> block_component(dim+1, 0);

    block_component[dim]=1;

    DoFRenumbering::component_wise(dof_handler, block_component);

    const std::vector<types::global_dof_index> dofs_per_block =

          DoFTools::count_dofs_per_fe_block(dof_handler, block_component);

    const unsigned int n_u = dofs_per_block[0],

                       n_p = dofs_per_block[1];

    std::cout << "Number of active cells: " << triangulation.n_active_cells()

              << std::endl

              << "Total number of cells: " << triangulation.n_cells()

              << std::endl

              << "Number of degrees of freedom: " << dof_handler.n_dofs()

              << " (" << n_u << '+' << n_p << ')' << std::endl;

    BlockDynamicSparsityPattern dsp(2, 2);

    dsp.block(0, 0).reinit(n_u, n_u);

    dsp.block(1, 0).reinit(n_p, n_u);

    dsp.block(0, 1).reinit(n_u, n_p);

    dsp.block(1, 1).reinit(n_p, n_p);

    dsp.collect_sizes();

    DoFTools::make_sparsity_pattern(dof_handler, dsp);

    sparsity_pattern.copy_from(dsp);

    system_matrix.reinit(sparsity_pattern);

    solution.reinit(2);

    solution.block(0).reinit(n_u);

    solution.block(1).reinit(n_p);

    solution.collect_sizes();

    system_rhs.reinit(2);

    system_rhs.block(0).reinit(n_u);

    system_rhs.block(1).reinit(n_p);

    system_rhs.collect_sizes();

  }

  template <int dim>

  void Poroelasticity<dim>::assemble_system()

  {

    QGauss<dim>     quadrature_formula(degree + 2);

    QGauss<dim - 1> face_quadrature_formula(degree + 2);

    FEValues<dim>     fe_values(fe,

                            quadrature_formula,

                            update_values | update_gradients |

                              update_quadrature_points | update_JxW_values);

    FEFaceValues<dim> fe_face_values(fe,

                                     face_quadrature_formula,

                                     update_values | update_normal_vectors |

                                       update_quadrature_points |

                                       update_JxW_values);

    const unsigned int dofs_per_cell   = fe.n_dofs_per_cell();

    const unsigned int n_q_points      = quadrature_formula.size();

    std::cout <<"N_q_points="<< n_q_points<< std::endl;

    const unsigned int n_face_q_points = face_quadrature_formula.size();

    std::cout <<"N_face_q_points="<< n_face_q_points<< std::endl;

    constexpr double porosity = 0.96;

    constexpr double resistivity  = 87000; //[N*s/m^4]

constexpr double tortuosity  = 2.52;

constexpr double visc_length  = 3.7e-05; //[m]

constexpr double therm_length  = 0.000116; //[m]

constexpr double loss_factor  = 0.01;

constexpr double E_modulus = 1/43e5; //[Pa]

constexpr double poisson = 0.3;

constexpr double frame_density = 31; //[kg/m^3]

constexpr double fluid_density = 1.2; //[kg/m^3]

constexpr double Prandtl_number = 0.71; //[-]

const std::complex<double> i = {0,1};

const std::complex<double> shear_modulus = (220.+i*22.)*1e-4; //[N*m-2]

//constexpr double thermal_conductivity = 2.6*1e-2; //[w*m-1*k-1]

constexpr double specific_heat_ratio = 1.4; //[cp/cv]

constexpr double ambient_pressure = 1.0132*1e5; //[Pa]

std::vector<int> freq(end_freq-1);

std::vector<double> omega(freq.size());

std::vector<std::complex<double>> visc_factor(freq.size());

std::vector<std::complex<double>> dyn_tortuosity(freq.size());

std::vector<std::complex<double>> density11(freq.size());

std::vector<std::complex<double>> density12(freq.size());

std::vector<std::complex<double>> density22(freq.size());

std::vector<std::complex<double>> density0(freq.size());

std::vector<std::complex<double>> bulk_fluid(freq.size());

std::vector<std::complex<double>> A(freq.size());

std::vector<std::complex<double>> Q(freq.size());

std::vector<std::complex<double>> R(freq.size());

std::vector<std::complex<double>> gamma(freq.size());

const std::complex<double> bulk_sk_vacuum = 2.*shear_modulus*(poisson+1.)/(3.*(1.-2.*poisson)); //[]

const std::complex<double> bulk_elast_solid = poisson*E_modulus/((poisson+1.)*(1.-2.*poisson))+2.*shear_modulus/3.; //[]

FullMatrix<std::complex<double>> local_matrix(dofs_per_cell, dofs_per_cell);

//Vector<std::complex<double>>     local_rhs(dofs_per_cell);

//std::cout << freq.size() << std::endl;

for (unsigned int j=0; j<end_freq-start_freq; ++j)

{

//freq.push_back(j);

freq[j]=j+1;

//std::cout<< freq[j]<<std::endl;

omega[j]=2*M_PI*freq[j];

//std::cout<< omega[j]<<std::endl;

visc_factor[j]=std::sqrt(1.+i*4.*std::pow(tortuosity, 2)*loss_factor*fluid_density*omega[j]/

(std::pow(resistivity*visc_length*porosity, 2)));

dyn_tortuosity[j]=tortuosity-i*porosity*resistivity*visc_factor[j]/(fluid_density*omega[j]);

//std::cout<< dyn_tortuosity[j]<<std::endl;

density11[j]=(1.-porosity)*frame_density+porosity*fluid_density*(dyn_tortuosity[j]-1.);

density12[j]=-porosity*fluid_density*(dyn_tortuosity[j]-1.);

density22[j]=porosity*fluid_density*dyn_tortuosity[j];

density0[j]=density11[j]-std::pow(density12[j],2)/density22[j];

//std::cout<< density0[j]<<std::endl;

bulk_fluid[j]=ambient_pressure/(1.-(specific_heat_ratio-1.)/(specific_heat_ratio*(8.*loss_factor*

std::sqrt(1.+std::pow(therm_length/4.,2)*i*omega[j]*fluid_density*std::pow(Prandtl_number,2)/loss_factor)/

(fluid_density*std::pow(Prandtl_number*therm_length,2)*i*omega[j])+1.)));

//std::cout<< bulk_fluid[j]<<std::endl;

A[j]=((1-porosity)*(1-porosity-bulk_sk_vacuum/bulk_elast_solid)*

bulk_elast_solid+porosity*(bulk_elast_solid/bulk_fluid[j])*bulk_sk_vacuum)/

(1.-porosity-bulk_sk_vacuum/bulk_elast_solid+porosity*bulk_elast_solid/bulk_fluid[j])-2.*shear_modulus/3.;

//std::cout<< A[j]<<std::endl;

Q[j]=((1-porosity-bulk_sk_vacuum/bulk_elast_solid)*porosity*bulk_elast_solid)/

(1-porosity-bulk_sk_vacuum/bulk_elast_solid+porosity*bulk_elast_solid/bulk_fluid[j]);

//std::cout<< Q[j]<<std::endl;

R[j]=std::pow(porosity,2)*bulk_elast_solid/

(1-porosity-bulk_sk_vacuum/bulk_elast_solid+porosity*bulk_elast_solid/bulk_fluid[j]);

//std::cout<< R[j]<<std::endl;

gamma[j]=porosity*(density12[j]/density22[j] - Q[j]/R[j]);

//std::cout<< gamma[j]<<std::endl;

}

 

    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

    const RightHandSide<dim> right_hand_side;

    std::vector<double>         rhs_values(n_q_points);

    const FEValuesExtractors::Vector velocities(0);

    const FEValuesExtractors::Scalar pressure(dim);

    for (unsigned int k=0; k< end_freq-start_freq; ++k)

    {

for (const auto &cell : dof_handler.active_cell_iterators())

  {

fe_values.reinit(cell);

local_matrix = 0;

local_rhs    = 0;

right_hand_side.value_list(fe_values.get_quadrature_points(),rhs_values);

for (unsigned int q = 0; q < n_q_points; ++q)

  for (unsigned int i = 0; i < dofs_per_cell; ++i)

{

  const Tensor<1, dim> phi_i_u = fe_values[velocities].value(i, q);

  const double phi_i_p     = fe_values[pressure].value(i, q);

  const Tensor<1,dim> div_phi_i_p = fe_values[pressure].gradient(i, q);

  for (unsigned int j = 0; j < dofs_per_cell; ++j)

{

  const Tensor<1, dim> phi_j_u =

fe_values[velocities].value(j, q);

  const double phi_j_p = fe_values[pressure].value(j, q);

  const Tensor<1,dim> div_phi_j_p = fe_values[pressure].gradient(j, q);

  const SymmetricTensor<2, dim>

    eps_phi_i = get_strain(fe_values, i, q),

    eps_phi_j = get_strain(fe_values, j, q);

  const SymmetricTensor<4, dim, std::complex<double>>

    material_matrix = get_stress_strain_tensor<dim>(/*lambda = */ A[k]-std::pow(Q[k],2)/R[k],

                        /*mu = */ 2.*shear_modulus);

  local_matrix(i, j) +=

(eps_phi_i * material_matrix * eps_phi_j // transpose 

- std::pow(omega[k],2)*density0[k]*phi_i_u*phi_j_u  //

- gamma[k]*div_phi_i_p*phi_j_u

+ std::pow(porosity,2)*div_phi_i_p*div_phi_j_p/density22[k]

- std::pow(porosity,2)*phi_i_p*phi_j_p/R[k]

- gamma[k]*div_phi_i_p*phi_j_u)                //

* fe_values.JxW(q);

}

local_rhs(i) += rhs_values[q] * fe_values.JxW(q);

}

//}

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],

  local_matrix(i, j));

std::cout << "finished"<< std::endl;

for (unsigned int i = 0; i < dofs_per_cell; ++i)

  system_rhs(local_dof_indices[i]) += local_rhs(i);

std::cout << "rhs"<< std::endl;

    }

  }

  }

  template <int dim>

  void Poroelasticity<dim>::solve()

  {

  std::cout << "Solving linear system..";

  for (unsigned int k=0; k<end_freq-start_freq; ++k)

  {

SparseDirectUMFPACK A_direct;

//A_direct.initialize(all_local_matrices[k]);

//A_direct.vmult(solution, system_rhs);

//A_direct.solve(all_local_matrices[k], system_rhs);

//all_solutions[k-start_freq]=solution;

A_direct.solve(system_matrix, system_rhs);

  }

  }

  template <int dim>

  void Poroelasticity<dim>::run()

  {

    make_grid_and_dofs();

    assemble_system();

    //solve();

    //output_results();

  }

} 

int main()

{

  try

    {

      using namespace Step20;

      const unsigned int     fe_degree = 1;

      Poroelasticity<2> poroelasticity(fe_degree);

      poroelasticity.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;

}

[simon...@gmail.com]

Hi,
"Invalid pointer" suggest that this is some bug related to memory. One way to find such errors is to turn on the adress sanitizer. If you use gcc or clang, this can be done by configuring cmake this:

cmake -DCMAKE_BUILD_TYPE="Debug" -DCMAKE_CXX_FLAGS="-fsanitize=address" -DCMAKE_EXE_LINKER_FLAGS="-fsanitize=address" ..

and then recompiling. If you now run your program, it will crash with something like this:

=================================================================
==29575==ERROR: AddressSanitizer: heap-buffer-overflow on address 0x602000007a34 at pc 0x55a3f14c9cb4 bp 0x7ffd346a5f40 sp 0x7ffd346a5f30
WRITE of size 4 at 0x602000007a34 thread T0
#0 0x55a3f14c9cb3 in Poro::Poroelasticity<2>::assemble_system() /home/simon/workspace/sandbox/example.cc:249
#1 0x55a3f14c110f in Poro::Poroelasticity<2>::run() /home/simon/workspace/sandbox/example.cc:394
#2 0x55a3f14b517a in main /home/simon/workspace/sandbox/example.cc:411
#3 0x7fea8d0a80b2 in __libc_start_main (/lib/x86_64-linux-gnu/libc.so.6+0x270b2)
#4 0x55a3f14b4f5d in _start (/home/simon/workspace/sandbox/build/example+0x76f5d)


At the top entry #0, you can see that the code crashed at line 249,  which corresponds to the line:

freq[j]=j+1;

Here, freq had length 1 and you by mistake tried to access the j=1 entry.

Best,
Simon

[simon...@gmail.com]

Hi,

Good that you solved you problem.

Address sanitizer is great for finding this type of bugs, but note that you don't want to use it when you run the actual simulation.

In general, memory usage increases (by a  factor ~2 to ~3) and the code runs slower (by a factor ~2 to ~5).

Best,

Simon