1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_GRID_TEST_CHECKCOMMUNICATE_HH
4 #define DUNE_GRID_TEST_CHECKCOMMUNICATE_HH
5
6 #include <iostream>
7 #include <fstream>
8 #include <vector>
9
10 #include <dune/common/fvector.hh>
11 #include <dune/geometry/referenceelements.hh>
12 #include <dune/grid/common/gridview.hh>
13 #include <dune/grid/common/capabilities.hh>
14 #include <dune/grid/common/datahandleif.hh>
15
16 /* Communication Test for Parallel Grids
17 * -------------------------------------
18 *
19 * For a fixed codimension c and a fixed upwind direction u, the test works
20 * as follows:
21 * 1) In the center of all upwind codim c subentities of the interioir codim 0
22 * leaf entites a function is stored. Also a flag is set to 1.
23 * The computation is also performed on the subentities of the physical
24 * boundary.
25 * -> For all leaf subentities of codim c the flag should be set to 1
26 * with the exception of the border subentities on the inflow
27 * processor boundary and in the ghost elements - on these
28 * subentities the flag is zero.
29 * 2) Exchange both the data and the flags.
30 * 3) Test if the flag for all leaf subentities of codim c is set to 1.
31 *
32 * Note: This test requires the normals on both sides of an intersection to
33 * sum up to zero, i.e., there is exactly one tangent plane to the grid
34 * in every point of the intersection (actually, the barycenter would be
35 * sufficient).
36 */
37
38
39 /*****
40 The exchange is done using the ExampleDataHandle given below.
41 Together with the function value and the flag the coordinates
42 of all corners of the subenties are transmitted, giving
43 the possibility for additional testing in the scatter/set
44 methods.
45 ******/
46 /*******************************************************************/
47 namespace
48 {
49
50 template< class Grid, int codim = Grid::dimension+1 >
51 struct NextCodim
52 {
53 static const bool canCommunicate = Dune::Capabilities::canCommunicate< Grid, codim-1 >::v;
54
55 static const int v = (canCommunicate ? codim-1 : NextCodim< Grid, codim-1 >::v);
56 };
57
58 template< class Grid >
59 struct NextCodim< Grid, 0 >
60 {
61 static const int v = -1;
62 };
63
64 }
65
66
67 template <class IndexSetImp,
68 class GlobalIdSetImp,
69 class DataVectorType >
70 class ExampleDataHandle
71 : public Dune::CommDataHandleIF< ExampleDataHandle< IndexSetImp, GlobalIdSetImp, DataVectorType >, typename DataVectorType::value_type >
72 {
73 const IndexSetImp & iset_;
74 const GlobalIdSetImp & ids_;
75 int cdim_;
76 DataVectorType & data1_;
77 DataVectorType & data2_;
78 public:
79 typedef typename DataVectorType :: value_type DataType;
ExampleDataHandle(const IndexSetImp & iset,const GlobalIdSetImp & ids,int cdim,DataVectorType & d1,DataVectorType & d2)80 ExampleDataHandle(const IndexSetImp & iset,
81 const GlobalIdSetImp & ids,
82 int cdim,
83 DataVectorType & d1, DataVectorType & d2) :
84 iset_(iset), ids_(ids) , cdim_(cdim), data1_(d1) , data2_(d2)
85 {}
86
87 //! returns true if data for this codim should be communicated
contains(int dim,int codim) const88 bool contains ([[maybe_unused]] int dim, int codim) const
89 {
90 return (codim==cdim_);
91 }
92
93 //! returns true if size per entity of given dim and codim is a constant
fixedSize(int dim,int codim) const94 bool fixedSize ([[maybe_unused]] int dim, [[maybe_unused]] int codim) const
95 {
96 // this problem is a fixed size problem,
97 // but to simulate also non-fixed size problems
98 // we set this to false, should work anyway
99 return false;
100 }
101
102 /*! how many objects of type DataType have to be sent for a given entity
103 Note: Only the sender side needs to know this size.
104 */
105 template<class EntityType>
size(EntityType & e) const106 size_t size (EntityType& e) const
107 {
108 // flag+data+coordinates
109 typedef typename EntityType::Geometry Geometry;
110 return 2+e.geometry().corners() * Geometry::coorddimension;
111 }
112
113 //! pack data from user to message buffer
114 template<class MessageBuffer, class EntityType>
gather(MessageBuffer & buff,const EntityType & e) const115 void gather (MessageBuffer& buff, const EntityType& e) const
116 {
117 int idx = iset_.index(e);
118
119 //typename GlobalIdSetImp :: IdType id = ids_.id( e );
120 //buff.write( id );
121 buff.write(data2_[idx]); // flag
122 buff.write(data1_[idx]); // data
123
124 // all corner coordinates
125 typedef typename EntityType::Geometry Geometry;
126 const Geometry &geometry = e.geometry();
127 for( int i = 0; i < geometry.corners(); ++i )
128 {
129 typedef Dune::FieldVector< typename Geometry::ctype, Geometry::coorddimension > Vector;
130 const Vector corner = geometry.corner( i );
131 for( int j = 0; j < Geometry::coorddimension; ++j )
132 buff.write( corner[ j ] );
133 }
134 }
135
136 /*! unpack data from message buffer to user
137 n is the number of objects sent by the sender
138 */
139 template<class MessageBuffer, class EntityType>
scatter(MessageBuffer & buff,const EntityType & e,size_t n)140 void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
141 {
142 using std::sqrt;
143 using Geometry = typename EntityType::Geometry;
144 using ctype = typename Geometry::ctype;
145
146 // define a tolerance for floating-point checks
147 const ctype tolerance = sqrt(std::numeric_limits< ctype >::epsilon());
148
149 // as this problem is a fixed size problem we can check the sizes
150 assert( n == size(e) );
151
152 // make sure that global id on all processors are the same
153 // here compare the id of the entity that sended the data with my id
154 //typename GlobalIdSetImp :: IdType id;
155 //buff.read( id );
156 //typename GlobalIdSetImp :: IdType myId = ids_.id( e );
157 //std::cout << id << " id | my Id = " << myId << "\n";
158 //assert( id == myId );
159
160 // else do normal scatter
161 int idx = iset_.index(e);
162 DataType x=0.0;
163 buff.read(x); // flag
164
165 // for ghost entities x > 0 must be true
166 assert( ( e.partitionType() == Dune::GhostEntity ) ? (x>=0.0) : 1);
167
168 if (x>=0)
169 { // only overwrite existing data if flag = 1, i.e.,
170 // the sending processor acctually computed the value
171 data2_[idx] = x;
172 x=0.;
173 buff.read(x); // correct function value
174 data1_[idx] = x;
175 }
176 else
177 {
178 x=0.;
179 buff.read(x);
180 }
181
182 // test if the sending/receiving entities are geometrically the same
183 const Geometry &geometry = e.geometry();
184 for( int i = 0; i < geometry.corners(); ++i )
185 {
186 typedef Dune::FieldVector< ctype, Geometry::coorddimension > Vector;
187 const Vector corner = geometry.corner( i );
188 for( int j = 0; j < Geometry::coorddimension; ++j )
189 {
190 buff.read(x);
191 if( std::abs( corner[ j ] - x ) > tolerance )
192 {
193 std::cerr << "ERROR in scatter: Vertex <" << i << "," << j << ">: "
194 << " this : (" << corner[ j ] << ")"
195 << " other : (" << x << ")"
196 << std::endl;
197 }
198 }
199 }
200 }
201
202 };
203
204
205 /*******************************************************************/
206 /*******************************************************************/
207 template< class GridView, int cdim, class OutputStream >
208 class CheckCommunication
209 {
210 typedef typename GridView :: Grid Grid;
211 typedef typename GridView :: IndexSet IndexSet;
212 typedef typename GridView :: IntersectionIterator IntersectionIterator;
213
214 typedef typename IntersectionIterator :: Intersection Intersection;
215
216 typedef typename GridView :: template Codim< 0 > :: Entity Entity;
217 typedef typename GridView :: template Codim< 0 > :: Iterator Iterator;
218
219 typedef typename GridView :: template Codim< cdim > :: Entity SubEntity;
220
221 enum { dimworld = Grid :: dimensionworld };
222 enum { dim = Grid :: dimension };
223
224 typedef typename Grid :: ctype ctype;
225
226 typedef Dune::FieldVector< ctype, dimworld > CoordinateVector;
227 typedef std :: vector< ctype > ArrayType;
228
229 CoordinateVector upwind_;
230 OutputStream &sout_;
231
232 const GridView &gridView_;
233 const IndexSet &indexSet_;
234 const int level_;
235
236 // the function
f(const CoordinateVector & x)237 ctype f ( const CoordinateVector &x )
238 {
239 CoordinateVector a( 1.0 );
240 a[0] = -0.5;
241 return a*x+1.5; //+cos(x*x);
242 }
243
244 // compute the data on the upwind entities
project(int dataSize,ArrayType & data,ArrayType & weight,int)245 void project ( int dataSize, ArrayType &data, ArrayType &weight, int /* rank */ )
246 {
247 // set initial data
248 for(int i=0 ; i<dataSize; ++i)
249 {
250 data[i] = 0.0;
251 weight[i] = -1.0;
252 }
253
254 Iterator end = gridView_.template end< 0 >();
255 for( Iterator it = gridView_.template begin< 0 >(); it != end ; ++it )
256 {
257 const Entity &entity = *it;
258
259 if( cdim == 0 )
260 {
261 CoordinateVector mid( 0.0 );
262 const int numVertices = entity.subEntities(dim);
263 for( int i = 0; i < numVertices; ++i )
264 mid += entity.geometry().corner( i );
265 mid /= ctype( numVertices );
266
267 int index = indexSet_.index( entity );
268 data[ index ] = f( mid );
269 weight[ index ] = 1.0;
270 }
271 else
272 {
273 const IntersectionIterator nend = gridView_.iend( entity );
274 for( IntersectionIterator nit = gridView_.ibegin( entity ); nit != nend; ++nit )
275 {
276 const Intersection &intersection = *nit;
277
278 const typename Intersection::LocalGeometry &geoInSelf = intersection.geometryInInside();
279
280 auto faceRefElement = referenceElement( geoInSelf );
281 const Dune::FieldVector< ctype, dim-1 > &bary = faceRefElement.position( 0, 0 );
282
283 const CoordinateVector normal = intersection.integrationOuterNormal( bary );
284 ctype calc = normal * upwind_;
285
286 // if testing level, then on non-conform grid also set values on
287 // intersections that are not boundary, but has no level
288 // neighbor
289 const bool proceedAnyway = (level_ < 0 ? false : !intersection.neighbor());
290 if( (calc > -1e-8) || intersection.boundary() || proceedAnyway )
291 {
292 auto insideRefElem = referenceElement( entity.geometry() );
293
294 const int indexInInside = intersection.indexInInside();
295 for( int i = 0; i < insideRefElem.size( indexInInside, 1, cdim ); ++i )
296 {
297 const int e = insideRefElem.subEntity( indexInInside, 1, i, cdim );
298 const int idx = indexSet_.subIndex( entity, e, cdim );
299 CoordinateVector cmid( 0.0 );
300 SubEntity subE = entity.template subEntity< cdim >( e );
301 const int c = subE.geometry().corners();
302 for( int j = 0; j < c; ++j )
303 cmid += subE.geometry().corner( j );
304 cmid /= ctype( c );
305
306 data[ idx ] = f( cmid );
307 weight[ idx ] = 1.0;
308 }
309
310 // on non-conforming grids the neighbor entities might not
311 // be the same as those on *it, therefore set data on neighbor
312 // as well
313 if( intersection.neighbor() )
314 {
315 Entity neigh = intersection.outside();
316
317 assert( (level_ < 0) ? (neigh.isLeaf()) : 1);
318 assert( (level_ < 0) ? 1 : (neigh.level() == level_) );
319
320 auto outsideRefElem = referenceElement( neigh.geometry() );
321
322 const int indexInOutside = intersection.indexInOutside();
323 for( int i = 0; i < outsideRefElem.size( indexInOutside, 1, cdim ); ++i )
324 {
325 const int e = outsideRefElem.subEntity( indexInOutside, 1, i, cdim );
326 const int idx = indexSet_.subIndex( neigh, e, cdim );
327 CoordinateVector cmid( 0.0 );
328 SubEntity subE = neigh.template subEntity< cdim >( e );
329 const int c = subE.geometry().corners();
330 for( int j = 0; j < c; ++j )
331 cmid += subE.geometry().corner( j );
332 cmid /= ctype( c );
333
334 data[ idx ] = f( cmid );
335 weight[ idx ] = 1.0;
336 }
337 }
338 }
339 }
340 }
341 }
342 }
343
344 // test if all flags are 1 and return the
345 // difference in the function values.
346 // if testweight is true an error is printed for each
347 // flag not equal to 1
test(int,ArrayType & data,ArrayType & weight,bool testweight)348 ctype test ( int /* dataSize */, ArrayType &data, ArrayType &weight, bool testweight )
349 {
350 const int rank = gridView_.comm().rank();
351 const int size = gridView_.comm().size();
352
353 // check whether there is an overlap or ghost region of
354 // cells for the current grid view.
355 if (size > 1 && cdim == 0 &&
356 gridView_.overlapSize(0) == 0 &&
357 gridView_.ghostSize(0) == 0)
358 {
359 sout_ << "<" << rank << "/test> Error in communication test.";
360 sout_ << " overlapSize+ghostSize:" << gridView_.overlapSize(0) + gridView_.ghostSize(0) << " (should not be 0)";
361 sout_ << " communicator size is:" << size;
362 sout_ << std :: endl;
363 }
364
365 //Variante MIT Geisterzellen
366 //typedef typename IndexSet :: template Codim<0> :: template Partition<All_Partition> :: Iterator IteratorType;
367
368 ctype maxerr = 0.0;
369 Iterator end = gridView_.template end< 0 >();
370 for( Iterator it = gridView_.template begin< 0 >(); it != end ; ++it )
371 {
372 const Entity &entity = *it;
373
374 CoordinateVector mid( 0.0 );
375 const int numVertices = entity.subEntities(dim);
376 for( int i = 0; i < numVertices; ++i )
377 mid += entity.geometry().corner( i );
378 mid /= ctype(numVertices);
379
380 if( cdim == 0 )
381 {
382 int index = indexSet_.index( entity );
383 ctype lerr = std::abs( f( mid ) - data[ index ] );
384 maxerr = std :: max( maxerr, lerr );
385 if( testweight && (weight[ index ] < 0) )
386 {
387 sout_ << "<" << rank << "/test> Error in communication test.";
388 sout_ << " weight:" << weight[ index ] << " (should be 0)";
389 sout_ << " value is : " << data[ index ];
390 sout_ << " index is: " << index;
391 sout_ << " level:" << entity.level();
392 sout_ << std :: endl;
393 }
394 }
395 else
396 {
397 const int numSubEntities = entity.subEntities(cdim);
398 for( int i=0; i < numSubEntities; ++i )
399 {
400 SubEntity subE = entity.template subEntity< cdim >( i );
401
402 const int index = indexSet_.index( subE );
403 CoordinateVector cmid( 0.0 );
404
405 const int numVertices = subE.geometry().corners();
406 for( int j = 0; j< numVertices; ++j )
407 cmid += subE.geometry().corner( j );
408 cmid /= ctype( numVertices );
409
410 ctype lerr = std::abs( f( cmid ) - data[ index ] );
411 maxerr = std::max( maxerr, lerr );
412 if( testweight && (weight[ index ] < 0) )
413 {
414 sout_ << "<" << rank << "/test> Error in communication test.";
415 sout_ << " weight:" << weight[ index ] << " should be zero!";
416 sout_ << " value is : " << data[ index ];
417 sout_ << " index is:" << index;
418 sout_ << " level: " << entity.level() ;
419 sout_ << std :: endl;
420
421 for( int j = 0; j < numVertices; )
422 {
423 auto refElem = referenceElement ( entity.geometry() );
424 const int vx = refElem.subEntity( i, cdim, j, dim );
425
426 sout_ << "index: " << indexSet_.subIndex( entity, vx, dim )
427 << " " << subE.geometry().corner( j );
428 (++j < numVertices ? sout_ << "/" : sout_ << std :: endl);
429 }
430 }
431 }
432 }
433 }
434 return maxerr;
435 }
436
437 // The main ''algorithm''
checkCommunication()438 bool checkCommunication ()
439 {
440 using std::sqrt;
441 // define a tolerance for floating-point checks
442 const ctype tolerance = sqrt(std::numeric_limits< ctype >::epsilon());
443
444 upwind_[ 0 ] = -0.1113;
445 int myrank = gridView_.comm().rank();
446
447 if( myrank == 0 )
448 {
449 std::cout << "TEST ";
450 (level_ < 0 ? std :: cout << "Leaf" : std :: cout << "Level<" << level_ << ">");
451 std :: cout << " communication for codim " << cdim << std :: endl;
452 }
453
454 const int dataSize = indexSet_.size( cdim );
455 ArrayType data(dataSize, 0.0);
456 ArrayType weight(dataSize, 0.0);
457 project( dataSize, data, weight, myrank );
458
459 ctype preresult = test( dataSize, data, weight, false );
460 sout_ << "Test before Communication on <" << myrank << "> " << preresult << std :: endl;
461 // Communicate
462 typedef typename Grid :: Traits :: GlobalIdSet GlobalIdSet;
463 ExampleDataHandle< IndexSet, GlobalIdSet, ArrayType >
464 dh( indexSet_, gridView_.grid().globalIdSet(), cdim, data, weight );
465
466 // call communication of grid
467 try
468 {
469 // call forward and backward communication
470 auto obj1 = gridView_.communicate( dh, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication );
471 if( ! obj1.ready() )
472 obj1.wait();
473
474 // make sure backward communication does the same, this should change nothing
475 auto obj2 = gridView_.communicate( dh, Dune::InteriorBorder_All_Interface, Dune::BackwardCommunication );
476 }
477 catch( const Dune::NotImplemented &exception )
478 {
479 if( myrank == 0 )
480 {
481 sout_ << "Error: Communication for codimension " << cdim
482 << " not implemented." << std :: endl;
483 sout_ << " (" << exception << ")" << std :: endl;
484 }
485 return false;
486 }
487
488 ctype result = test( dataSize, data, weight, true );
489 sout_ << "Test after Communication on <" << myrank << "> " << result << std :: endl;
490 return (std::abs(result) < tolerance);
491 }
492
493 public:
CheckCommunication(const GridView & gridView,OutputStream & sout,int level)494 CheckCommunication ( const GridView &gridView, OutputStream &sout, int level )
495 : upwind_( -1.0 ),
496 sout_( sout ),
497 gridView_( gridView ),
498 indexSet_( gridView_.indexSet() ),
499 level_( level )
500 {
501 // if no overlap and ghost is available we skip the check
502 const bool skipCheck = ( cdim == 0 ) ? (gridView_.overlapSize(0) == 0 && gridView_.ghostSize(0) == 0) : false ;
503
504 if( skipCheck )
505 {
506 std :: cerr << "Codim " << cdim << ": Test skiped because of empty set of overlap and ghosts !" << std :: endl;
507 }
508 else if ( ! checkCommunication() )
509 {
510 std :: cerr << "Error in communication test for codim "
511 << cdim << "!" << std :: endl;
512 }
513
514 // for automatic testing of all codims
515 CheckCommunication< GridView, NextCodim< Grid, cdim >::v, OutputStream >
516 test( gridView_, sout_, level_ );
517 }
518 };
519
520
521 template< class GridView, class OutputStream >
522 class CheckCommunication< GridView, -1, OutputStream >
523 {
524 public:
CheckCommunication(const GridView &,OutputStream &,int)525 CheckCommunication ( const GridView &, OutputStream &, int)
526 {}
527 };
528
529
530 template< class Grid, class OutputStream >
checkCommunication(const Grid & grid,int level,OutputStream & sout)531 void checkCommunication( const Grid &grid, int level, OutputStream &sout )
532 {
533 if( level < 0 )
534 {
535 typedef typename Grid::LeafGridView GridView;
536 GridView gridView = grid.leafGridView();
537 CheckCommunication< GridView, NextCodim< Grid >::v, OutputStream >
538 test( gridView, sout, level );
539 }
540 else
541 {
542 typedef typename Grid::LevelGridView GridView;
543 GridView gridView = grid.levelGridView( level );
544 CheckCommunication< GridView, NextCodim< Grid >::v, OutputStream >
545 test( gridView, sout, level );
546 }
547 }
548
549 #endif // DUNE_GRID_TEST_CHECKCOMMUNICATE_HH
550