[Dune-devel] SOCIS 2014: FEM assembly with threads

[Agnese, Marco]

Hi DUNErs,
as suggested by Markus I have written a very simple code which uses threads to assemble the stiffness matrix and the RHS of a FEM scheme for a 1D Poisson problem.

You can find the code here: https://github.com/magnese/dune-common/blob/threads/dune/common/test/threadstest.cc

The code constructs a vector of grid points. The mesh is "partitioned" among threads. Each thread assembles a local stiffness matrix and than copy the values in the global stiffness matrix A. The only race condition which could happen is for the grid points which are the boundaries of the partitions. Therefore if we are using N threads the number "critical" points is N-1. To avoid a race condition I have created a mutex for each "critical" point (which is marked as shared). If an entry of type shared need to be inserted, assemble() takes the corresponding lock. This choice works in this case because the number of critical point is very small compared to the size of the vector. In general this approach is very inefficient. 

Now I am writing a more general class to obtain something similar of what we have for ParallelIndexSet. 
Let's suppose that we have a vector v and we want to call several threads which will read and write on this vector. The class should work like this:
1) the user specifies all the positions for each thread which can be accessed by the thread
2) the class will find the which part of the vector are critical and will create several mutex to address the problem (computed only once)
3) the user give to the class the function to be called by the threads. This function can only read from the global vector and return a buffer of values which need to be inserted in the global vector
4) when the buffers are created, the values are inserted in the vectors directly if they aren't in a critical position, otherwise using the mutexs.
Using the class, the user avoids to directly deal with threads and mutex. 

Obviously the difficult part it to implement an efficient and general lock mechanism. As soon as I have some code I will let you know. 

Therefore my objective right now is to implement this class and test it with the assemble(). What do you think? I would appreciate very much some feedbacks to understand if I am moving in the right direction or not. 

Cheers 
Marco.