Data Race Condition while using AddElementMatrixSymmetric

jhauser
Topic Author
Offline
New Member

5 years 10 months ago #1363 by jhauser

Data Race Condition while using AddElementMatrixSymmetric was created by jhauser

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<int Dy>
  void CalcElementMatrix_tpy(const FiniteElement & fely, 
                            const ElementTransformation &eltrans_y, LocalHeap & lhy, 
                            FlatMatrix<double> elmaty)
  {
     const HCurlFiniteElement<Dy>  & fely_cast = 
                            dynamic_cast<const HCurlFiniteElement<Dy>&>(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<GridFunction> MyMaxwellAssemble(shared_ptr<FESpace> fes_in,
                                      shared_ptr<CoefficientFunction> 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<DofId> dofsy = ely.GetDofs();
     FlatMatrix<> elmat_fesy(dofsy.Size(),lhy); 
                 
     CalcElementMatrix_tpy<2>(fely,eltrans_y,lhy,elmat_fesy);
     maty->AddElementMatrixSymmetric(dofsy,elmat_fesy);
  });
....

}

christopher
Offline
Administrator

5 years 10 months ago #1364 by christopher

Replied by christopher on topic Data Race Condition while using AddElementMatrixSymmetric

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

jhauser
Topic Author
Offline
New Member

5 years 10 months ago - 5 years 10 months ago #1365 by jhauser

Replied by jhauser on topic Data Race Condition while using AddElementMatrixSymmetric

Yes, I use the standard HCurl space:

Code:

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

Last edit: 5 years 10 months ago by jhauser.

joachim
Offline
Administrator

5 years 10 months ago #1366 by joachim

Replied by joachim on topic Data Race Condition while using AddElementMatrixSymmetric

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

jhauser
Topic Author
Offline
New Member

5 years 10 months ago #1367 by jhauser

Replied by jhauser on topic Data Race Condition while using AddElementMatrixSymmetric

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):

Code:

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

Code:

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.

Attachments:

joachim
Offline
Administrator

5 years 10 months ago #1368 by joachim

Replied by joachim on topic Data Race Condition while using AddElementMatrixSymmetric

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

Code:

matrix.SetZero()

Time to create page: 0.106 seconds