<div dir="ltr">Hi again<div><br></div><div>I managed to create an index set and rebuild the RemoteIndices such that the code compiles and run without errors.</div><div><br></div><div>I am not sure exactly how to use the library though. The system has size 2*Sim.nc x 2*Sim.nc and is created with opm autodiff and hence has the <a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">SparseMatrix</a> format from Eigen. On each processor I create a Communication object and create an index set similar to Is from Figure 1 in <a href="https://www.dune-project.org/doc/comm/communication.pdf">communicator.pdf</a> with the following code.</div><div><br></div><div><div><span class="" style="white-space:pre"> </span>int rank;</div><div><span class="" style="white-space:pre"> </span>int procs;</div><div><span class="" style="white-space:pre"> </span>MPI_Comm_rank(MPI_COMM_WORLD, &rank);</div><div><span class="" style="white-space:pre"> </span>MPI_Comm_size(MPI_COMM_WORLD, &procs);</div></div><div><br></div><div><div><span style="white-space:pre"> </span>/* Create communication structure */</div><div><span class="" style="white-space:pre"> </span>Comm comm(MPI_COMM_WORLD);</div><div><span class="" style="white-space:pre"> </span></div><div><span class="" style="white-space:pre"> </span>typedef Dune::OwnerOverlapCopyAttributeSet::AttributeSet AS;</div><div><span class="" style="white-space:pre"> </span>typedef Dune::ParallelLocalIndex<AS> LocalIndex;</div><div><span class="" style="white-space:pre"> </span>typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet PIndexSet;</div><div><span class="" style="white-space:pre"> </span>typedef Dune::RemoteIndices<PIndexSet> RIndices;</div><div><br></div><div><span class="" style="white-space:pre"> </span>PIndexSet &IS = comm.indexSet();</div><div><span class="" style="white-space:pre"> </span>int elem_per_proc = 2*Sim.nc/procs; /* Truncated */</div><div><span class="" style="white-space:pre"> </span>int offset = 0;</div><div><br></div><div><span class="" style="white-space:pre"> </span>IS.beginResize();</div><div><br></div><div><span class="" style="white-space:pre"> </span>int start = elem_per_proc*rank;</div><div><br></div><div><span class="" style="white-space:pre"> </span>/* Last processor takes the rest */</div><div><span class="" style="white-space:pre"> </span>int end = elem_per_proc*(rank+1);</div><div><span class="" style="white-space:pre"> </span>if(rank == procs-1)</div><div><span class="" style="white-space:pre"> </span>end = 2*Sim.nc;</div><div><br></div><div><span class="" style="white-space:pre"> </span>/* Add the previous element as a ghost point if not on first processor */</div><div><span class="" style="white-space:pre"> </span>if(rank > 0){</div><div><span class="" style="white-space:pre"> </span>IS.add(start-1, LocalIndex(0, AS::overlap));</div><div><span class="" style="white-space:pre"> </span>offset = 1;</div><div><span class="" style="white-space:pre"> </span>}</div><div><br></div><div><span class="" style="white-space:pre"> </span>/* Indices owned by this process */</div><div><span class="" style="white-space:pre"> </span>for(int i = start; i < end; ++i){</div><div><span class="" style="white-space:pre"> </span>IS.add(i, LocalIndex(i-start+offset, AS::owner));</div><div><span class="" style="white-space:pre"> </span>}</div><div><br></div><div><span class="" style="white-space:pre"> </span>/* Add the following element as a ghost point if not on last processor */</div><div><span class="" style="white-space:pre"> </span>if(rank < procs-1)</div><div><span class="" style="white-space:pre"> </span>IS.add(end, LocalIndex(end-start+offset, AS::overlap));</div><div><br></div><div><span class="" style="white-space:pre"> </span>IS.endResize();</div><div><br></div><div><span class="" style="white-space:pre"> </span>RIndices &RI = comm.remoteIndices();</div><div><span class="" style="white-space:pre"> </span>RI.rebuild<true>();</div></div><div><br></div><div>When I need to solve the linear system I first create an OPM <a href="http://www.opm-project.org/documentation/master/opm-autodiff/html/class_opm_1_1_dune_matrix.php">DuneMatrix</a> which I then use to construct an object of the type on each processor</div><div><br></div><div>Dune::OverlappingSchwarzOperator<Opm::DuneMatrix, Vector, Vector, Comm><br></div><div><br></div><div>Then I construct the sequential and parallel preconditioners and the scalarproduct with the types</div><div><br></div><div><div>typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;</div><div>typedef Dune::BlockVector<Dune::FieldVector<double,1>> Vector;</div><div>typedef Dune::OverlappingSchwarzOperator<Opm::DuneMatrix, Vector, Vector, Comm> Matrix;</div><div><br></div><div>typedef Dune::OverlappingSchwarzScalarProduct<Vector, Comm> ScalarProduct;</div><div>typedef Dune::SeqILU0<Opm::DuneMatrix, Vector, Vector> SeqPrec;</div><div>typedef Dune::BlockPreconditioner<Vector, Vector, Comm, SeqPrec> ParPrec;</div></div><div><br></div><div>Lastly I call apply. I also tried to run all code but the part creating the communicator and the index set on the processor with rank 0 but then the program hangs when apply is called on the solver object.</div><div><br></div><div>I guess I am missing some synchronisation or communication between the processors? For instance, when apply is called, I guess it only updates the part that is owned by that processor. By the way, I realize that speedup attained from this approach may be limited since the system matrix is created sequentially. It is more to show the possibility of using Dune ISTL in parallel.</div></div><div class="gmail_extra"><br><div class="gmail_quote">On 11 June 2015 at 11:54, Tobias Ritschel <span dir="ltr"><<a href="mailto:tobiasritschel@gmail.com" target="_blank">tobiasritschel@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Thanks, I will have a look at the documentation. I solve for both pressure and saturation in each cell, so I I'll create the index set.</div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On 11 June 2015 at 11:36, Markus Blatt <span dir="ltr"><<a href="mailto:markus@dr-blatt.de" target="_blank">markus@dr-blatt.de</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span>On Thu, Jun 11, 2015 at 10:38:10AM +0200, Tobias Ritschel wrote:<br>
> Ah, this is clearly my error. I am using the constructor<br>
><br>
> OwnerOverlapCopyCommunication<br>
</span>> <<a href="http://www.dune-project.org/doc/doxygen/dune-istl-html/class_dune_1_1_owner_overlap_copy_communication.html#a93f446939bd4d6623bdf2e544f183f8c" target="_blank">http://www.dune-project.org/doc/doxygen/dune-istl-html/class_dune_1_1_owner_overlap_copy_communication.html#a93f446939bd4d6623bdf2e544f183f8c</a>><br>
> (MPI_Comm<br>
> comm_, SolverCategory::Category<br>
> <<a href="http://www.dune-project.org/doc/doxygen/dune-istl-html/struct_dune_1_1_solver_category.html#ae061380ac961490be6ee353cf0dc1733" target="_blank">http://www.dune-project.org/doc/doxygen/dune-istl-html/struct_dune_1_1_solver_category.html#ae061380ac961490be6ee353cf0dc1733</a>><br>
<span>> cat_=SolverCategory::overlapping,<br>
> bool freecomm_=false)<br>
><br>
> and used the created object without an index set. Is the information<br>
> automatically set up when constructing an IndexInfoFromGrid object, or<br>
> should they be set with AddLocalIndex and AddRemoteIndex? And is this the<br>
> most convenient way to set up the indices?<br>
<br>
</span>Don't use IndexInfoGrid. I assume that is a left over from ancient<br>
dune-disc times. You can access the underlying index set by<br>
OwnerOverlapCopyCommunication::indexSet() and use that to add the<br>
indices with the correct tag. Later on get the remote indices with<br>
OwnerOverlapCopyCommunication::remoteIndices() and call rebuild on it,<br>
<br>
You will find more information on the index sets and how they are used<br>
in dune-istl in the "Description of the parallel communication<br>
interface" on <a href="https://dune-project.org/doc/index.html" target="_blank">https://dune-project.org/doc/index.html</a>. You need one<br>
index per degree of freedom and use the attribute owner to mark a<br>
partitioning. In the easiest case you need to mark connected degrees<br>
of freedom that are not owner on the process as copy.<br>
<br>
If you are using dune-cornerpoint with one unknown per cell, then you<br>
can use ParallelISTLInformation::copyValuesTo from opm-core to do the job.<br>
<br>
Otherwise there is some manual work to do. E.g. compute the owner<br>
region using grid communication.<br>
<span><font color="#888888"><br>
Markus<br>
</font></span><span><br>
><br>
> On 11 June 2015 at 09:59, Markus Blatt <<a href="mailto:markus@dr-blatt.de" target="_blank">markus@dr-blatt.de</a>> wrote:<br>
><br>
> > Hi,<br>
> ><br>
> > On Wed, Jun 10, 2015 at 09:00:36PM +0200, Tobias Ritschel wrote:<br>
> > ><br>
> > > terminate called after throwing an instance of<br>
> > > 'Dune::InterfaceBuilder::RemotexIndicesStateError'<br>
> > ><br>
> > > As far as I can see, it is called when I call the member function apply<br>
> > on<br>
> > > the RestartedGMRES solver object.<br>
> > ><br>
> ><br>
> > This only happens if the information in the RemoteIndices is not in<br>
> > sync with the underlying IndexSets. This can be case if these were<br>
> > changed, but there was not call to RemoteIndices::rebuild().<br>
> > Thi is also reflected in the error message: "RemoteIndices is not in<br>
> > sync with the index set. Call RemoteIndices::rebuild first!"<br>
> ><br>
> > So the question is: How do you set up the index sets and the remote<br>
> > information that you are using in OwnerOverlapCopyCommunication?<br>
> ><br>
> > Cheers,<br>
> ><br>
> > Markus<br>
> ><br>
> ><br>
</span><div><div>> > Dr. Markus Blatt - HPC-Simulation-Software & Services<br>
> > <a href="http://www.dr-blatt.de" target="_blank">http://www.dr-blatt.de</a><br>
> > Hans-Bunte-Str. 8-10, 69123 Heidelberg, Germany<br>
> > Tel.: <a href="tel:%2B49%20%280%29%20160%2097590858" value="+4916097590858" target="_blank">+49 (0) 160 97590858</a><br>
> ><br>
> > _______________________________________________<br>
> > Dune mailing list<br>
> > <a href="mailto:Dune@dune-project.org" target="_blank">Dune@dune-project.org</a><br>
> > <a href="http://lists.dune-project.org/mailman/listinfo/dune" target="_blank">http://lists.dune-project.org/mailman/listinfo/dune</a><br>
> ><br>
> ><br>
<br>
--<br>
Do you need more support with DUNE or HPC in general?<br>
<br>
Dr. Markus Blatt - HPC-Simulation-Software & Services <a href="http://www.dr-blatt.de" target="_blank">http://www.dr-blatt.de</a><br>
Hans-Bunte-Str. 8-10, 69123 Heidelberg, Germany<br>
Tel.: <a href="tel:%2B49%20%280%29%20160%2097590858" value="+4916097590858" target="_blank">+49 (0) 160 97590858</a><br>
</div></div><br>_______________________________________________<br>
Dune mailing list<br>
<a href="mailto:Dune@dune-project.org" target="_blank">Dune@dune-project.org</a><br>
<a href="http://lists.dune-project.org/mailman/listinfo/dune" target="_blank">http://lists.dune-project.org/mailman/listinfo/dune</a><br>
<br></blockquote></div><br></div>
</div></div></blockquote></div><br></div>