Data Race Condition while using AddElementMatrixSymmetric

Main Forum Old-Forum - do not post here
1

I am still working on the implementation of the assembling in tensor product spaces. For that I need the local mass matrix of HCurl-functions. My resulting matrix seems mostly right, but there are a few entries that are wrong. It changes each time which entries are influenced. The error only occurs when I use a python script with python3. If I use the interface of jupyter notebook, my matrix is correct.
The local matices seem to be correct therefore I guess that it is a problem with AddElementMatrixSymmetric or rather a data race condition in IterateElements.
How can I slove it?

The code with this problem is as follows:

[code] template
void CalcElementMatrix_tpy(const FiniteElement & fely,
const ElementTransformation &eltrans_y, LocalHeap & lhy,
FlatMatrix elmaty)
{
const HCurlFiniteElement & fely_cast =
dynamic_cast<const HCurlFiniteElement&>(fely);
ELEMENT_TYPE eltypey = fely_cast.ElementType();

 int ndy = fely_cast.GetNDof();
 const IntegrationRule & ir_voly = 
                        SelectIntegrationRule (eltypey, 2*fely_cast.Order());

 MatrixFixWidth<Dy,double> shapey(ndy);
		
 elmaty = 0.0;
 for (int l = 0; l < ir_voly.Size(); l++)
 {
   HeapReset hr(lhy);
   const MappedIntegrationPoint<Dy,Dy> sip(ir_voly[l], eltrans_y);
   BaseMappedIntegrationPoint & mip = eltrans_y(ir_voly[l], lhy);

   fely_cast.CalcMappedShape (mip, shapey);

   elmaty  += ir_voly[l].Weight() *sip.GetMeasure()*  shapey * Trans (shapey);
 }

}

...

shared_ptr MyMaxwellAssemble(shared_ptr fes_in,
shared_ptr f_coef)
{
…

IterateElements(*fes_y, VOL, lhy, [&] (FESpace::Element ely, LocalHeap &lhy)
{
const ElementTransformation& eltrans_y = ma_y->GetTrafo(ely, lhy);
const FiniteElement& fely = fes_y->GetFE(ely,lhy);
FlatArray dofsy = ely.GetDofs();
FlatMatrix<> elmat_fesy(dofsy.Size(),lhy);

 CalcElementMatrix_tpy<2>(fely,eltrans_y,lhy,elmat_fesy);
 maty->AddElementMatrixSymmetric(dofsy,elmat_fesy);

});
…

}[/code]

2

Hi,
Do you use a standard HCurl space? IterateElements uses the coloring of the finite element space to decide which elements can be assembled at the same time.
Best
Christopher

3

Yes, I use the standard HCurl space:

meshx = Mesh(unit_square.GenerateMesh(maxh=0.5))
fesx = HCurl(meshx,order=1)
4

Hi,

Are you using the TaskManager ?
If not, IterateElements is not parallel anyway.
If yes, what happens without ?

Is the x-space or the y-space the HCurl ? The C++ part and the python part look inconsistent.

Your code snipped looks reasonable. It’s difficult to help with this information. Maybe you can upload a full example producing this random results.

Joachim

5

No, I do not parallelize it. That’s why I have no clue what’s the real problem.

Here a shortened version of my code, which still produces the error on my computer:

I am currently working in the ml-ngs file system. The two files which I’ve attached, I put into /1_Basic/4_utility_functions. The following lines I added in utility_function.hpp (in /1_Basic/4_utility_functions):

namespace mymaxassembleSMALL
{
  void MyMaxwellAssemble_small(shared_ptr<FESpace> fes_in,
            shared_ptr<CoefficientFunction> eps_coef,
            shared_ptr<CoefficientFunction> mu_coef,
            shared_ptr<CoefficientFunction> f_coef);
}

and added

m.def("MyMaxwellAssembleSMALL", &mymaxassembleSMALL::MyMaxwellAssemble_small);

in myngspy.cpp (in /1_Basic) to the list of exported functions in the end of the file. Of course I added “4_utility_functions/myMaxAssembling_small.cpp” in CMakeLists.txt to “add_ngsolve_python_module” as well.

Attachment: myMaxAssembling_small.cpp

Attachment: maxwell_testSMALL.py

6

Hi,
after creation of a sparse matrix, it is NOT zero-initialized.
You have to call

matrix.SetZero()
7

With this I get the error:
‘class std::shared_ptr<ngla::SparseMatrixSymmetric >’ has no member named ‘SetZero’

I used

for(int i=0; i<ndof_y;i++) { FlatVector<double> mat_rowvals = maty->GetRowValues(i); for (int k=0; k<mat_rowvals.Size(); k++) { mat_rowvals[k]=0.0; } }
as a work around, but I still get a similar wrong matrix. Different every time I run the program.
At least with the python script. (I use python 3 to execute it) If I use “MyMaxwellAssembleSMALL” in jupyter notebook, I get the correct matrix. (which is printed in my terminal, when I use this function)

8 
matrix -> SetZero()

if matrix is a shared pointer

9

Oh, you’re absolutely right. Thank you!
However the problem still remains. This doesn’t change the every changing entries of the matrix. It’s still different every time I run it.

10

With maty->SetZero() I get consistent results, can you repost the new version of your code?

This is because jupyter pipes only the python stdout and not the c++ one. You can use

pybind11::print(maty)

to get the output in jupyter (you need to include pybind11/pybind11.h to have that).

11

I cleaned out the cmake-files of ml-ngs and build/compiled everything anew.
Now it works. I have absolutly no idea why it didn’t work before since I haven’t changed the code.
Thank you very much for your time! I am sorry, that it took so long.

And thank you for the jupyter suggestion.