dune-fem 2.8.0
Loading...
Searching...
No Matches
molgalerkin.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_SCHEMES_MOLGALERKIN_HH
2#define DUNE_FEM_SCHEMES_MOLGALERKIN_HH
3
4// fem includes
8
9namespace Dune
10{
11
12 namespace Fem
13 {
14
15#if HAVE_PETSC
16 //forward declaration of PetscLinearOperator
17 template <typename DomainFunction, typename RangeFunction >
18 class PetscLinearOperator ;
19#endif
20
21 // GalerkinOperator
22 // ----------------
23
24 template< class Integrands, class DomainFunction, class RangeFunction = DomainFunction >
26 : public virtual Operator< DomainFunction, RangeFunction >
27 {
28 typedef DomainFunction DomainFunctionType;
29 typedef RangeFunction RangeFunctionType;
30
31
32 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value, "DomainFunction and RangeFunction must be defined on the same grid part." );
33
34 typedef typename RangeFunctionType::GridPartType GridPartType;
36
37 typedef Impl::GalerkinOperator< Integrands > GalerkinOperatorImplType;
38 typedef typename RangeFunctionType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
39
40 typedef typename GalerkinOperatorImplType::template QuadratureSelector<
42
44
45 template< class... Args >
46 explicit MOLGalerkinOperator ( const GridPartType &gridPart, Args &&... args )
48 impl_( gridPart, std::forward< Args >( args )... ),
50 communicate_( true )
51 {
52 }
53
54 void setCommunicate( const bool communicate )
55 {
56 communicate_ = communicate;
58 {
59 std::cout << "MOLGalerkinOperator::setCommunicate: communicate was disabled!" << std::endl;
60 }
61 }
62
63 void setQuadratureOrders(unsigned int interior, unsigned int surface)
64 {
65 size_t size = impl_.size();
66 for( size_t i=0; i<size; ++i )
67 impl_[ i ].setQuadratureOrders(interior,surface);
68 }
69
70 virtual void operator() ( const DomainFunctionType &u, RangeFunctionType &w ) const final override
71 {
72 evaluate( u, w );
73 }
74
75 template< class GridFunction >
76 void operator() ( const GridFunction &u, RangeFunctionType &w ) const
77 {
78 evaluate( u, w );
79 }
80
81 const GridPartType &gridPart () const { return impl().gridPart(); }
82
83 typedef Integrands ModelType;
84 typedef Integrands DirichletModelType;
85 ModelType &model() const { return impl().model(); }
86 const GalerkinOperatorImplType& impl() const { return (*impl_); }
87
88 std::size_t gridSizeInterior () const { return gridSizeInterior_; }
89
90 protected:
91 // update number of interior elements as sum over threads
92 std::size_t gatherGridSizeInterior () const
93 {
94 std::size_t gridSizeInterior = 0;
95 const size_t size = MPIManager::numThreads();
96 for( size_t i=0; i<size; ++i )
98 return gridSizeInterior;
99 }
100
101 template <class Iterators>
102 void applyInverseMass( const Iterators& iterators, RangeFunctionType& w ) const
103 {
104 // temporary local function
106 TemporaryLocalFunctionType wLocal( w.space() );
107 LocalMassMatrixType localMassMatrix( w.space(), impl().interiorQuadratureOrder( w.space().order() ) );
108
109 // iterate over all elements (in the grid or per thread)
110 // thread safety is guaranteed through discontinuous data (spaces)
111 for( const auto& entity : iterators )
112 {
113 // obtain local function
114 auto guard = bindGuard( wLocal, entity );
115 // obtain local dofs
116 w.getLocalDofs( entity, wLocal.localDofVector() );
117
118 // apply inverse mass matrix
119 // TODO: add mass term if needed (from ufl expression)
120 localMassMatrix.applyInverse( entity, wLocal );
121
122 // overwrite dofs in discrete function
123 w.setLocalDofs( entity, wLocal.localDofVector() );
124 }
125 }
126
127 template< class GridFunction >
128 void evaluate( const GridFunction &u, RangeFunctionType &w ) const
129 {
131 w.clear();
132
133 std::shared_mutex mutex;
134
135 auto doEval = [this, &u, &w, &mutex] ()
136 {
137 // version with locking
138 this->impl().evaluate( u, w, this->iterators_, mutex );
139
140 // version without locking
141 //RangeFunctionType wTmp( w );
142 //this->impl().evaluate( u, wTmp, this->iterators_ );
143 //std::lock_guard guard ( mutex );
144 //w += wTmp;
145 };
146
147 bool singleThreadModeError = false ;
148
149 try {
150 // execute in parallel
151 MPIManager :: run ( doEval );
152
153 // update number of interior elements as sum over threads
155 }
156 catch ( const SingleThreadModeError& e )
157 {
158 singleThreadModeError = true;
159 }
160
161 // method of lines
162 auto doInvMass = [this, &w] ()
163 {
164 this->applyInverseMass( this->iterators_, w );
165 };
166
167 if( ! singleThreadModeError )
168 {
169 try {
170 // execute in parallel
171 MPIManager :: run ( doInvMass );
172 }
173 catch ( const SingleThreadModeError& e )
174 {
175 singleThreadModeError = true;
176 }
177 }
178
179 // if error occurred, redo the whole evaluation
180 if( singleThreadModeError )
181 {
182 // reset w from previous entries
183 w.clear();
184 // re-run in single thread mode if previous attempt failed
185 impl().evaluate( u, w, iterators_ );
187
188 // update number of interior elements
189 gridSizeInterior_ = impl().gridSizeInterior();
190 }
191
192 // synchronize data
193 if( communicate_ )
194 w.communicate();
195 }
196
197 // GalerkinOperator implementation (see galerkin.hh)
200
201 mutable std::size_t gridSizeInterior_;
203 };
204
205
206
207 // DifferentiableGalerkinOperator
208 // ------------------------------
209
210 template< class Integrands, class JacobianOperator >
212 : public MOLGalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
213 public DifferentiableOperator< JacobianOperator >
214 {
216
217 public:
218 typedef JacobianOperator JacobianOperatorType;
219
222 typedef typename DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType;
223 typedef typename RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType;
224
226
227 template< class... Args >
229 const RangeDiscreteFunctionSpaceType &rSpace,
230 Args &&... args )
231 : BaseType( rSpace.gridPart(), std::forward< Args >( args )... ),
232 dSpace_(dSpace),
233 rSpace_(rSpace),
234 stencilDAN_(dSpace,rSpace), stencilD_(dSpace,rSpace),
235 domainSpaceSequence_(dSpace.sequence()),
236 rangeSpaceSequence_(rSpace.sequence())
237 {}
238
239 virtual void jacobian ( const DomainFunctionType &u, JacobianOperatorType &jOp ) const final override
240 {
241 // assemble Jacobian, same as GalerkinOperator
242 assemble( u, jOp );
243 }
244
245 template< class GridFunction >
246 void jacobian ( const GridFunction &u, JacobianOperatorType &jOp ) const
247 {
248 assemble( u, jOp );
249 }
250
252 {
253 return dSpace_;
254 }
256 {
257 return rSpace_;
258 }
259
261
262 protected:
267
268 void prepare( JacobianOperatorType& jOp ) const
269 {
270 if ( domainSpaceSequence_ != domainSpace().sequence()
271 || rangeSpaceSequence_ != rangeSpace().sequence() )
272 {
273 domainSpaceSequence_ = domainSpace().sequence();
274 rangeSpaceSequence_ = rangeSpace().sequence();
275 if( impl().model().hasSkeleton() )
277 else
279 }
280 if( impl().model().hasSkeleton() )
281 jOp.reserve( stencilDAN_ );
282 else
283 jOp.reserve( stencilD_ );
284 // set all entries to zero
285 jOp.clear();
286 }
287
288 template < class GridFunction >
289 void assemble( const GridFunction &u, JacobianOperatorType &jOp ) const
290 {
291 prepare( jOp );
293
294 bool singleThreadModeError = false;
295 std::shared_mutex mutex;
296
297 auto doAssemble = [this, &u, &jOp, &mutex] ()
298 {
299 // assemble Jacobian, same as GalerkinOperator
300 this->impl().assemble( u, jOp, this->iterators_, mutex );
301 };
302
303 try {
304 // execute in parallel
305 MPIManager :: run ( doAssemble );
306
307 // update number of interior elements as sum over threads
309 }
310 catch ( const SingleThreadModeError& e )
311 {
312 singleThreadModeError = true;
313 }
314
315 if (!singleThreadModeError)
316 {
318
319 // method of lines
320 auto doInvMass = [this, &getSet] ()
321 {
322 applyInverseMass( this->iterators_, getSet, this->impl().model().hasSkeleton() );
323 };
324
325 try {
326 // execute in parallel
327 MPIManager :: run ( doInvMass );
328 }
329 catch ( const SingleThreadModeError& e )
330 {
331 singleThreadModeError = true;
332 }
333 }
334
335 if (singleThreadModeError)
336 {
337 // redo matrix assembly since it failed
338 jOp.clear();
339 impl().assemble( u, jOp, iterators_ );
340
341 // update number of interior elements
342 gridSizeInterior_ = impl().gridSizeInterior();
343
344 GetSetLocalMatrixImpl getSet( jOp );
345
346 // apply inverse mass
347 applyInverseMass( iterators_, getSet, impl().model().hasSkeleton() );
348 }
349
350 // note: assembly done without local contributions so need
351 // to call flush assembly
352 jOp.flushAssembly();
353 }
354
355
357 {
360
361 template <class Entity, class LocalMatrix>
362 void getLocalMatrix( const Entity& domainEntity, const Entity& rangeEntity,
363 LocalMatrix& lop ) const
364 {
365 jOp_.getLocalMatrix( domainEntity, rangeEntity, lop );
366 }
367
368 template <class Entity, class LocalMatrix>
369 void setLocalMatrix( const Entity& domainEntity, const Entity& rangeEntity,
370 const LocalMatrix& lop )
371 {
372 jOp_.setLocalMatrix( domainEntity, rangeEntity, lop );
373 }
374 };
375
377 {
379 typedef std::shared_mutex mutex_t;
381
383
384 template <class Entity, class LocalMatrix>
385 void getLocalMatrix( const Entity& domainEntity, const Entity& rangeEntity,
386 LocalMatrix& lop ) const
387 {
388 std::lock_guard< mutex_t > lock( mutex_ );
389 BaseType::getLocalMatrix( domainEntity, rangeEntity, lop );
390 }
391
392 template <class Entity, class LocalMatrix>
393 void setLocalMatrix( const Entity& domainEntity, const Entity& rangeEntity,
394 const LocalMatrix& lop )
395 {
396 std::lock_guard< mutex_t > lock( mutex_ );
397 BaseType::setLocalMatrix( domainEntity, rangeEntity, lop );
398 }
399 };
400
401 template <class JacOp>
403 {
406 };
407
408#if HAVE_PETSC
409 //- PetscLinearOperator does not work with threaded applyInverseMass
410 //- because the mix of getLocalMatrix and setLocalMatrix seems to be
411 //- problematic. Therefore we used the locked version.
412 template <class DomainFunction, class RangeFunction>
413 struct GetSetLocalMatrix< PetscLinearOperator< DomainFunction, RangeFunction > > : public GetSetLocalMatrixImplLocked
414 {
415 typedef GetSetLocalMatrixImplLocked BaseType;
416 GetSetLocalMatrix( JacobianOperatorType &jOp ) : BaseType( jOp ) {}
417 };
418#endif
419
420 template <class Iterators, class GetSetLocalMatrixOp >
421 void applyInverseMass ( const Iterators& iterators, GetSetLocalMatrixOp& getSet,
422 const bool hasSkeleton ) const
423 {
425 JacobianOperatorType& jOp = getSet.jOp_;
426
427 LocalMassMatrixType localMassMatrix( jOp.rangeSpace(), impl().interiorQuadratureOrder( jOp.rangeSpace().order() ) );
428
430 lop(jOp.domainSpace(), jOp.rangeSpace());
431
432 // multiply with inverse mass matrix
433 for( const auto& inside : iterators )
434 {
435 // scale diagonal
436 {
437 auto guard = bindGuard( lop, inside,inside );
438 lop.bind(inside,inside);
439 getSet.getLocalMatrix( inside, inside, lop );
440 localMassMatrix.leftMultiplyInverse( lop );
441 getSet.setLocalMatrix( inside, inside, lop );
442 }
443
444 if( hasSkeleton )
445 {
446 for( const auto &intersection : intersections( gridPart(), inside ) )
447 {
448 // scale off-diagonal
449 if( intersection.neighbor() )
450 {
451 const auto& outside = intersection.outside();
452 auto guard = bindGuard( lop, outside,inside );
453 getSet.getLocalMatrix( outside, inside, lop );
454 localMassMatrix.leftMultiplyInverse( lop );
455 getSet.setLocalMatrix( outside, inside, lop );
456 }
457 }
458 }
459 }
460 }
461
464
468 };
469
470
471
472 // AutomaticDifferenceGalerkinOperator
473 // -----------------------------------
474
475 template< class Integrands, class DomainFunction, class RangeFunction >
477 : public MOLGalerkinOperator< Integrands, DomainFunction, RangeFunction >,
478 public AutomaticDifferenceOperator< DomainFunction, RangeFunction >
479 {
482
483 public:
485
486 template< class... Args >
488 : BaseType( gridPart, std::forward< Args >( args )... ), AutomaticDifferenceOperatorType()
489 {}
490 };
491
492
493
494 // ModelDifferentiableGalerkinOperator
495 // -----------------------------------
496
497 template < class LinearOperator, class ModelIntegrands >
499 : public MOLDifferentiableGalerkinOperator< ModelIntegrands, LinearOperator >
500 {
502
503 typedef typename ModelIntegrands::ModelType ModelType;
504
506 typedef typename LinearOperator::RangeSpaceType DiscreteFunctionSpaceType;
507
509 : BaseType( dfSpace.gridPart(), model )
510 {}
511
512 template< class GridFunction >
513 void apply ( const GridFunction &u, RangeFunctionType &w ) const
514 {
515 (*this)( u, w );
516 }
517
518 template< class GridFunction >
519 void apply ( const GridFunction &u, LinearOperator &jOp ) const
520 {
521 (*this).jacobian( u, jOp );
522 }
523 };
524
525
526 // MethodOfLinesScheme
527 // -------------------
528
529 template< class Integrands, class LinearOperator, class InverseOperator, bool addDirichletBC>
530 using MethodOfLinesScheme = Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC,
532
533 } // namespace Fem
534
535} // namespace Dune
536
537#endif // #ifndef DUNE_FEM_SCHEMES_MOLGALERKIN_HH
STL namespace.
Definition: bindguard.hh:11
Impl::GalerkinSchemeImpl< Integrands, LinearOperator, InverseOperator, addDirichletBC, MOLDifferentiableGalerkinOperator > MethodOfLinesScheme
Definition: molgalerkin.hh:531
static auto bindGuard(Object &object, Args &&... args) -> std::enable_if_t< isBindable< Object, Args... >::value, BindGuard< Object > >
Definition: bindguard.hh:67
A temporary function carrying values for one entity.
Definition: temporary.hh:208
static bool verbose()
obtain the cached value for fem.verbose
Definition: io/parameter.hh:445
Exception thrown when a code segment that is supposed to be only accessed in single thread mode is ac...
Definition: mpimanager.hh:43
static int numThreads()
return number of current threads
Definition: mpimanager.hh:421
void update()
update internal list of iterators
Definition: threaditerator.hh:88
size_t size() const
return number of threads
Definition: threadsafevalue.hh:48
void applyInverse(MassCaller &caller, const EntityType &entity, const BasisFunctionSet &basisFunctionSet, LocalFunction &lf) const
Definition: localmassmatrix.hh:285
void leftMultiplyInverse(LocalMatrix &localMatrix) const
compute M^-1 * localMatrix
Definition: localmassmatrix.hh:339
Local Mass Matrix for arbitrary spaces.
Definition: localmassmatrix.hh:909
operator providing a Jacobian through automatic differentiation
Definition: automaticdifferenceoperator.hh:86
abstract differentiable operator
Definition: differentiableoperator.hh:29
abstract operator
Definition: operator.hh:34
DomainFunction DomainFunctionType
type of discrete function in the operator's domain
Definition: operator.hh:36
abstract affine-linear operator
Definition: operator.hh:87
void update()
clear previously computed entries such that a re-compute happens when used again
Definition: stencil.hh:136
A local matrix with a small array as storage.
Definition: temporarylocalmatrix.hh:100
void bind(const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity)
initialize the local matrix to entities
Definition: temporarylocalmatrix.hh:156
Definition: molgalerkin.hh:27
ThreadSafeValue< GalerkinOperatorImplType > impl_
Definition: molgalerkin.hh:199
bool communicate_
Definition: molgalerkin.hh:202
LocalMassMatrix< DiscreteFunctionSpaceType, InteriorQuadratureType > LocalMassMatrixType
Definition: molgalerkin.hh:43
ThreadIteratorType iterators_
Definition: molgalerkin.hh:198
void applyInverseMass(const Iterators &iterators, RangeFunctionType &w) const
Definition: molgalerkin.hh:102
ModelType & model() const
Definition: molgalerkin.hh:85
const GalerkinOperatorImplType & impl() const
Definition: molgalerkin.hh:86
std::size_t gridSizeInterior_
Definition: molgalerkin.hh:201
const GridPartType & gridPart() const
Definition: molgalerkin.hh:81
void setCommunicate(const bool communicate)
Definition: molgalerkin.hh:54
RangeFunction RangeFunctionType
Definition: molgalerkin.hh:29
void evaluate(const GridFunction &u, RangeFunctionType &w) const
Definition: molgalerkin.hh:128
virtual void operator()(const DomainFunctionType &u, RangeFunctionType &w) const final override
application operator
Definition: molgalerkin.hh:70
ThreadIterator< GridPartType > ThreadIteratorType
Definition: molgalerkin.hh:35
std::size_t gridSizeInterior() const
Definition: molgalerkin.hh:88
RangeFunctionType::GridPartType GridPartType
Definition: molgalerkin.hh:34
DomainFunction DomainFunctionType
Definition: molgalerkin.hh:28
Impl::GalerkinOperator< Integrands > GalerkinOperatorImplType
Definition: molgalerkin.hh:37
RangeFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
Definition: molgalerkin.hh:38
MOLGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: molgalerkin.hh:46
void setQuadratureOrders(unsigned int interior, unsigned int surface)
Definition: molgalerkin.hh:63
Integrands DirichletModelType
Definition: molgalerkin.hh:84
Integrands ModelType
Definition: molgalerkin.hh:83
std::size_t gatherGridSizeInterior() const
Definition: molgalerkin.hh:92
GalerkinOperatorImplType::template QuadratureSelector< DiscreteFunctionSpaceType >::InteriorQuadratureType InteriorQuadratureType
Definition: molgalerkin.hh:41
Definition: molgalerkin.hh:214
DiagonalAndNeighborStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > stencilDAN_
Definition: molgalerkin.hh:465
int domainSpaceSequence_
Definition: molgalerkin.hh:467
ThreadIteratorType iterators_
Definition: molgalerkin.hh:198
MOLDifferentiableGalerkinOperator(const DomainDiscreteFunctionSpaceType &dSpace, const RangeDiscreteFunctionSpaceType &rSpace, Args &&... args)
Definition: molgalerkin.hh:228
DiagonalStencil< DomainDiscreteFunctionSpaceType, RangeDiscreteFunctionSpaceType > stencilD_
Definition: molgalerkin.hh:466
const RangeDiscreteFunctionSpaceType & rSpace_
Definition: molgalerkin.hh:463
const GalerkinOperatorImplType & impl() const
Definition: molgalerkin.hh:86
std::size_t gridSizeInterior_
Definition: molgalerkin.hh:201
virtual void jacobian(const DomainFunctionType &u, JacobianOperatorType &jOp) const final override
obtain linearization
Definition: molgalerkin.hh:239
const GridPartType & gridPart() const
Definition: molgalerkin.hh:81
const DomainDiscreteFunctionSpaceType & dSpace_
Definition: molgalerkin.hh:462
RangeFunctionType::DiscreteFunctionSpaceType RangeDiscreteFunctionSpaceType
Definition: molgalerkin.hh:223
void prepare(JacobianOperatorType &jOp) const
Definition: molgalerkin.hh:268
void assemble(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: molgalerkin.hh:289
void applyInverseMass(const Iterators &iterators, GetSetLocalMatrixOp &getSet, const bool hasSkeleton) const
Definition: molgalerkin.hh:421
const DomainDiscreteFunctionSpaceType & domainSpace() const
Definition: molgalerkin.hh:251
DomainFunctionType::DiscreteFunctionSpaceType DomainDiscreteFunctionSpaceType
Definition: molgalerkin.hh:222
int rangeSpaceSequence_
Definition: molgalerkin.hh:467
BaseType::GridPartType GridPartType
Definition: molgalerkin.hh:225
JacobianOperator JacobianOperatorType
Definition: molgalerkin.hh:218
BaseType::RangeFunctionType RangeFunctionType
Definition: molgalerkin.hh:221
BaseType::DomainFunctionType DomainFunctionType
Definition: molgalerkin.hh:220
void jacobian(const GridFunction &u, JacobianOperatorType &jOp) const
Definition: molgalerkin.hh:246
const RangeDiscreteFunctionSpaceType & rangeSpace() const
Definition: molgalerkin.hh:255
std::size_t gatherGridSizeInterior() const
Definition: molgalerkin.hh:92
GetSetLocalMatrixImpl(JacobianOperatorType &jOp)
Definition: molgalerkin.hh:359
void setLocalMatrix(const Entity &domainEntity, const Entity &rangeEntity, const LocalMatrix &lop)
Definition: molgalerkin.hh:369
void getLocalMatrix(const Entity &domainEntity, const Entity &rangeEntity, LocalMatrix &lop) const
Definition: molgalerkin.hh:362
JacobianOperatorType & jOp_
Definition: molgalerkin.hh:358
std::shared_mutex mutex_t
Definition: molgalerkin.hh:379
GetSetLocalMatrixImpl BaseType
Definition: molgalerkin.hh:378
GetSetLocalMatrixImplLocked(JacobianOperatorType &jOp)
Definition: molgalerkin.hh:382
void getLocalMatrix(const Entity &domainEntity, const Entity &rangeEntity, LocalMatrix &lop) const
Definition: molgalerkin.hh:385
void setLocalMatrix(const Entity &domainEntity, const Entity &rangeEntity, const LocalMatrix &lop)
Definition: molgalerkin.hh:393
GetSetLocalMatrixImpl BaseType
Definition: molgalerkin.hh:404
GetSetLocalMatrix(JacobianOperatorType &jOp)
Definition: molgalerkin.hh:405
BaseType::GridPartType GridPartType
Definition: molgalerkin.hh:484
MOLAutomaticDifferenceGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: molgalerkin.hh:487
MOLModelDifferentiableGalerkinOperator(ModelType &model, const DiscreteFunctionSpaceType &dfSpace)
Definition: molgalerkin.hh:508
LinearOperator::DomainFunctionType RangeFunctionType
Definition: molgalerkin.hh:505
void apply(const GridFunction &u, RangeFunctionType &w) const
Definition: molgalerkin.hh:513
MOLDifferentiableGalerkinOperator< ModelIntegrands, LinearOperator > BaseType
Definition: molgalerkin.hh:501
LinearOperator::RangeSpaceType DiscreteFunctionSpaceType
Definition: molgalerkin.hh:506
void apply(const GridFunction &u, LinearOperator &jOp) const
Definition: molgalerkin.hh:519
ModelIntegrands::ModelType ModelType
Definition: molgalerkin.hh:503