Calculating mass matrix for vectorial base, memory leak?

[Karol Lewandowski]

Hello, 
I am trying to compute a mass matrix using a form integrator:

using OpMassPrec = FormsIntegrators<BoundaryEleOp>::Assembly<

USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 3>;

I noticed that for some runs of the program, the resulting matrix has NaN values. 
I suspect some kind of memory leak. The integrate function looks ok, any idea what might be the problem?

Link to the function.
https://bitbucket.org/likask/mofem-cephas/src/master/mofem/src/finite_elements/BiLinearFormsIntegrators.hpp#lines-531

[Lukasz Kaczmraczyk]

Hi Karol,

Can you give more context. This problem is 2d or 3d, such that we can reproduce the problem. Is resulting matrix has nan values, or element s of the base vector looks as a nan. 

We using that matrix for two other problems, with Hdiv space, and the error should pop out there.

Lukasz

[Karol Lewandowski]

It is 3D problem and it looks like the base vectors has nan values, or at least that is what debugger is showing me. 
Yes, I know that you are using this specialisation in few places, I am going to test that. 

[Lukasz Kaczmraczyk]

Hi,

If debugger is showing nan on line 563, or 565, where

++t_col_base

or 

++t_row_base;

then you can see nan if this is the last element in the array if you do increment you step outside, but this is the end of the loop at the same time. So it does not matter.

What base it is, for Hex, or Tet?

Lukasz

[Lukasz Kaczmraczyk]

Karol,

I will propose to run Valgrind on dev server, to check if this is a problem.

L.

[Karol Lewandowski]

Yes, I will test this. Knowing life, I expect the problem to not appear on the server xD 

[Lukasz Kaczmraczyk]

Karol, please give me code, and all that is needed to reproduce this problem.

[Karol Lewandowski]

karol/develop branch on multifield plasticity, line ElasticOperators.hpp:1185

in param file just change the path to the input file

and run the code by 

./multifield_plasticity 

[Lukasz Kaczmraczyk]

Karol, 

I just realized looking into your code, that this is the wrong specification in your context, that is why is not working, and you have a strange error. We should add error handling with manage misuse like that one.  

Look at OP defnition:

template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase> struct OpMassImpl;

You have vector base BASE_DIM = 3, and approximate vectorial field, so FIELD_DIM = 3, 

using OpMassPrec = FormsIntegrators<BoundaryEleOp>::Assembly<

USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 3>;

But this is not the case, in fact, you approximate tensorial field, SIGMA is tensor of rank 2, so your field in fact has dimension 9, since, is 3x3 components. We do not have this specialization implemented, but that would be easy to do.

Lukasz

On Thursday, September 30, 2021 at 1:03:16 PM UTC+1 inz.karol...@gmail.com wrote:

[Karol Lewandowski]

Ok, so it should be 

using OpMassPrec = FormsIntegrators<BoundaryEleOp>::Assembly<

USER_ASSEMBLE>::BiLinearForm<GAUSS>::OpMass<3, 9>;

[Karol Lewandowski]

No ok, we need a new template

[Lukasz Kaczmraczyk]

Correct. 

And we need to add an error message which checks the rank of the field, to throw an error in cases of misuse. In the OP doWork. Somthing like that. 

#ifndef NDEBUG

if(data.getFieldDofs()[0]->getNbOfCoeffs() != 1) 

  SETERRQ(

    PETSC_COMM_SELF, MOFEM_DATA_INONSISTENCY,  

    "Field should have vector base (e.g. Hdiv/Hcurl) and it has to have 1 coefficient for each base function");

#endif

So in this case, if you have a vector base i.e. BASE_DIM=3, if you have one coefficient for each base function, your field dimension is 3. However you have data.getFieldDofs()[0]->getNbOfCoeffs() ==3.

I expect, that this error will happen from time to time.

Lukasz

[Karol Lewandowski]

do we have NDEBUG flag defined in MoFEM? 

How will the assembly look like for OpMassImpl<3, 9, GAUSS, OpBase> specialisation? 

[Lukasz Kaczmraczyk]

Hello,

We do not have to use NDEBUG, but sometimes is to do check only for debugging. Here is not so important, since check is an outside loop over integration points. However, I will tend to make all checks like that in operators to make only when debugging.

For this, you have to do to thinks,  first take only diagonal from the element mass matrix, 

 FTensor::Index<'i', 3> i;

  auto get_t_vec = [&](const int rr) {

    std::array<double *, 3> ptrs;

    for (auto i = 0; i != 3; ++i)

      ptrs[i] = &OpBase::locMat(rr + i, i);

    FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_vec(&ptrs[0], &ptrs[1], &ptrs[2);;

    return t_vec;

  };

and then in the loop over gauss pts, and base functions on rows and columns, do

        t_vec(i) += alpha * (t_row_base(k) * t_col_base(l));

Very similar to what is here,

http://mofem.eng.gla.ac.uk/mofem/html/struct_mo_f_e_m_1_1_op_mass_impl_3_011_00_01_f_i_e_l_d___d_i_m_00_01_g_a_u_s_s_00_01_op_base_01_4.html

Regards,

L.

[Lukasz Kaczmraczyk]

Small fix, should be,

t_vec(i) += alpha * (t_row_base(k) * t_col_base(k));