<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
<style type="text/css" style="display:none;"><!-- P {margin-top:0;margin-bottom:0;} --></style>
</head>
<body dir="ltr">
<div id="divtagdefaultwrapper" dir="ltr" style="font-size: 12pt; color: rgb(0, 0, 0); font-family: Calibri, Helvetica, sans-serif, EmojiFont, "Apple Color Emoji", "Segoe UI Emoji", NotoColorEmoji, "Segoe UI Symbol", "Android Emoji", EmojiSymbols;">
<p></p>
<div>Dear Dune developers,</div>
<div>there seems to be an issue when using the Dune::GlobalIndexSet class (<dune/grid/utility/globalindexset.hh>) if we compute an SPgrid on a high number of threads.</div>
<div>Indeed, GlobalIndexSet seems to assume that each thread will have a specific "region" of the grid. Therefore, it gets the global index of the grid <span style="font-size:12pt">by using the local grid index + an offset.</span></div>
<div>However, SPgrids can be used with a high number of threads, which leads to a high number of overlaps and the cells owned by the different threads are intertwined. Consequently, Dune::GlobalIndexSet yields the wrong global index.</div>
<div><br>
</div>
<div>I adapted a dumux introductory exercise to showcase this problem (see zip folder attached and the summary given bellow). <span style="font-size:12pt">Even though it uses DuMux, I hope it is still informative.</span></div>
<div><br>
</div>
<div>I have the following question on this issue:</div>
<div>1) did I evaluate the issue causing the error in Dune::GlobalIndexSet correctly? </div>
<div>2) it seems that Dune::GlobalIndexSet is only used within parmetisgridpartitioner::repartition() which is itself not used by any other dune modules. Is this correct? </div>
<div>3) Consequently, we had planned on using a custom mapper in our implementation(/python binding) to go from '<span style="font-family:Calibri,Helvetica,sans-serif,EmojiFont,"Apple Color Emoji","Segoe UI Emoji",NotoColorEmoji,"Segoe UI Symbol","Android Emoji",EmojiSymbols; font-size:16px">globalIdSet'
to contiguous global indexes</span>. Once this is implemented are there other problems that might arise when using SPgrids with a high number of threads? </div>
<div><br>
</div>
<div>Many thanks for your help and with best regards</div>
<div><br>
</div>
<div>Mona</div>
<div><br>
</div>
<div>Summary of the example:</div>
<div><br>
</div>
<div>######### core code:</div>
<div><br>
</div>
<div><br>
</div>
<div>int main(int argc, char** argv)</div>
<div>{</div>
<div> using namespace Dumux;</div>
<div><br>
</div>
<div> </div>
<div> </div>
<div> // define the type tag for this problem</div>
<div> using TypeTag = Properties::TTag::Injection2pCC;</div>
<div><br>
</div>
<div> /// initialize MPI+X, finalize is done automatically on exit</div>
<div> Dumux::initialize(argc, argv);</div>
<div><br>
</div>
<div> // parse command line arguments and input file</div>
<div> Parameters::init(argc, argv);</div>
<div><br>
</div>
<div> // try to create a grid (from the given grid file or the input file)</div>
<div> GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;</div>
<div> gridManager.init();</div>
<div><br>
</div>
<div> ////////////////////////////////////////////////////////////</div>
<div> // run instationary non-linear problem on this grid</div>
<div> ////////////////////////////////////////////////////////////</div>
<div><br>
</div>
<div> // we compute on the leaf grid view</div>
<div> const auto& leafGridView = gridManager.grid().leafGridView();</div>
<div><br>
</div>
<div> // create the finite volume grid geometry</div>
<div> using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;</div>
<div> auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);</div>
<div> </div>
<div> // create the finite volume grid geometry</div>
<div> using SoilGridType = GetPropType<TypeTag, Properties::Grid>;</div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>using GridView = typename SoilGridType::Traits::LeafGridView;</span></div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>int dimWorld = GridView::dimensionworld;</span></div>
<div><span style="white-space:pre; white-space:normal"></span></div>
<div> </div>
<div> typedef typename GridView::Grid Grid;</div>
<div> typedef typename Grid::GlobalIdSet GlobalIdSet;</div>
<div><br>
</div>
<div> const GlobalIdSet & globalIdSet = gridManager.grid().leafGridView().grid().globalIdSet();/** retrieve globally unique Id set */</div>
<div><br>
</div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>/** local id**/</span></div>
<div> const typename GridView::IndexSet& indexSet = gridManager.grid().leafGridView().indexSet();</div>
<div><span style="white-space:normal"><span style="white-space:pre"></span> </span></div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>/** retrieve contiguous global unique Id set via GlobalIndexSet*/ </span></div>
<div> auto cellIdx = std::make_shared<Dune::GlobalIndexSet<GridView>>(gridManager.grid().leafGridView(), 0);</div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>std::cout<<"thread rank:"<<leafGridView.comm().rank()<<" dimWorld:"<<dimWorld <<"\n";</span></div>
<div> for (const auto& e : elements(gridGeometry->gridView())) {</div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>auto p = e.geometry().center(); // get cell center to compaire with indexes</span></div>
<div> auto gId2 = globalIdSet.id(e); // real global index (not contiguous)</div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>auto lId = indexSet.index(e); // local index</span></div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>int gIdx = cellIdx->index(e);//get the cell global index</span></div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>std::cout<<"rank "<<leafGridView.comm().rank()</span></div>
<div><span style="white-space:normal"><span style="white-space:pre"></span><<"<span style="white-space:pre">
</span>GlobalIndexSet_gIdx "<<gIdx<<" real_gIdx "<<gId2<<" lId "<<lId;</span></div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>std::cout<<" partitionType "<<e.partitionType() ;//check if cell is owned by the thread (i.e., e.partitionType() == interior)</span></div>
<div> for (int i=0; i< 3 ; i++) { // print cell center coordinates</div>
<div><span style="white-space:normal"><span style="white-space:pre"></span>std::cout<<", dim "<<i<<": "<<p[i];</span></div>
<div> }std::cout<<std::endl;</div>
<div> }<span style="white-space:pre"> </span></div>
<div> return 0;</div>
<div>} // end main</div>
<div><br>
</div>
<div><br>
</div>
<div>#########</div>
<div>inputs:</div>
<div><br>
</div>
<div>grid parameter:</div>
<div><br>
</div>
<div>[Grid]</div>
<div>LowerLeft = 0 0 0</div>
<div>UpperRight = 6 6 6</div>
<div>Cells = 2 2 2</div>
<div><br>
</div>
<div>#########</div>
<div>output:</div>
<div><br>
</div>
<div>for 1 thread:</div>
<div><br>
</div>
<div>rank 0 GlobalIndexSet_gIdx 0 real_gIdx 31 lId 0 partitionType interior, dim 0: 1.5, dim 1: 1.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 1 real_gIdx 33 lId 1 partitionType interior, dim 0: 4.5, dim 1: 1.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 2 real_gIdx 41 lId 2 partitionType interior, dim 0: 1.5, dim 1: 4.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 3 real_gIdx 43 lId 3 partitionType interior, dim 0: 4.5, dim 1: 4.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 4 real_gIdx 81 lId 4 partitionType interior, dim 0: 1.5, dim 1: 1.5, dim 2: 4.5</div>
<div>rank 0 GlobalIndexSet_gIdx 5 real_gIdx 83 lId 5 partitionType interior, dim 0: 4.5, dim 1: 1.5, dim 2: 4.5</div>
<div>rank 0 GlobalIndexSet_gIdx 6 real_gIdx 91 lId 6 partitionType interior, dim 0: 1.5, dim 1: 4.5, dim 2: 4.5</div>
<div>rank 0 GlobalIndexSet_gIdx 7 real_gIdx 93 lId 7 partitionType interior, dim 0: 4.5, dim 1: 4.5, dim 2: 4.5</div>
<div><br>
</div>
<div>for 2 threads (both threads give the same output except for the partitionType. Example for thread 0):</div>
<div>rank 0 GlobalIndexSet_gIdx 0 real_gIdx 31 lId 0 partitionType interior, dim 0: 1.5, dim 1: 1.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 4 real_gIdx 33 lId 1 partitionType overlap, dim 0: 4.5, dim 1: 1.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 1 real_gIdx 41 lId 2 partitionType interior, dim 0: 1.5, dim 1: 4.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 5 real_gIdx 43 lId 3 partitionType overlap, dim 0: 4.5, dim 1: 4.5, dim 2: 1.5</div>
<div>rank 0 GlobalIndexSet_gIdx 2 real_gIdx 81 lId 4 partitionType interior, dim 0: 1.5, dim 1: 1.5, dim 2: 4.5</div>
<div>rank 0 GlobalIndexSet_gIdx 6 real_gIdx 83 lId 5 partitionType overlap, dim 0: 4.5, dim 1: 1.5, dim 2: 4.5</div>
<div>rank 0 GlobalIndexSet_gIdx 3 real_gIdx 91 lId 6 partitionType interior, dim 0: 1.5, dim 1: 4.5, dim 2: 4.5</div>
<div>rank 0 GlobalIndexSet_gIdx 7 real_gIdx 93 lId 7 partitionType overlap, dim 0: 4.5, dim 1: 4.5, dim 2: 4.5</div>
<div><br>
</div>
<div><br>
</div>
<div>As can be seen, when we use several threads, the real_gIdx follows the expected order (from bottom left to top right), but the GlobalIndexSet_gIdx changes.</div>
<div>###########</div>
<br>
<p></p>
</div>
<br>
<font face="Arial" color="Black" size="1"><br>
------------------------------------------------------------------------------------------------<br>
------------------------------------------------------------------------------------------------<br>
Forschungszentrum Jülich GmbH<br>
52425 Jülich<br>
Sitz der Gesellschaft: Jülich<br>
Eingetragen im Handelsregister des Amtsgerichts Düren Nr. HR B 3498<br>
Vorsitzender des Aufsichtsrats: MinDir Stefan Müller<br>
Geschäftsführung: Prof. Dr. Astrid Lambrecht (Vorsitzende),<br>
Karsten Beneke (stellv. Vorsitzender), Dr. Ir. Pieter Jansens<br>
------------------------------------------------------------------------------------------------<br>
------------------------------------------------------------------------------------------------<br>
</font>
</body>
</html>