1#ifndef DUNE_FEM_SCHEMES_MOLGALERKIN_HH
2#define DUNE_FEM_SCHEMES_MOLGALERKIN_HH
17 template <
typename DomainFunction,
typename RangeFunction >
18 class PetscLinearOperator ;
24 template<
class Integrands,
class DomainFunction,
class RangeFunction = DomainFunction >
26 :
public virtual Operator< DomainFunction, RangeFunction >
32 static_assert( std::is_same< typename DomainFunctionType::GridPartType, typename RangeFunctionType::GridPartType >::value,
"DomainFunction and RangeFunction must be defined on the same grid part." );
40 typedef typename GalerkinOperatorImplType::template QuadratureSelector<
45 template<
class... Args >
59 std::cout <<
"MOLGalerkinOperator::setCommunicate: communicate was disabled!" << std::endl;
66 for(
size_t i=0; i<size; ++i )
75 template<
class Gr
idFunction >
96 for(
size_t i=0; i<size; ++i )
101 template <
class Iterators>
106 TemporaryLocalFunctionType wLocal( w.space() );
111 for(
const auto& entity : iterators )
114 auto guard =
bindGuard( wLocal, entity );
116 w.getLocalDofs( entity, wLocal.localDofVector() );
123 w.setLocalDofs( entity, wLocal.localDofVector() );
127 template<
class Gr
idFunction >
133 std::shared_mutex mutex;
135 auto doEval = [
this, &u, &w, &mutex] ()
147 bool singleThreadModeError = false ;
151 MPIManager :: run ( doEval );
158 singleThreadModeError =
true;
162 auto doInvMass = [
this, &w] ()
167 if( ! singleThreadModeError )
171 MPIManager :: run ( doInvMass );
175 singleThreadModeError =
true;
180 if( singleThreadModeError )
210 template<
class Integrands,
class JacobianOperator >
212 :
public MOLGalerkinOperator< Integrands, typename JacobianOperator::DomainFunctionType, typename JacobianOperator::RangeFunctionType >,
227 template<
class... Args >
245 template<
class Gr
idFunction >
288 template <
class Gr
idFunction >
294 bool singleThreadModeError =
false;
295 std::shared_mutex mutex;
297 auto doAssemble = [
this, &u, &jOp, &mutex] ()
305 MPIManager :: run ( doAssemble );
312 singleThreadModeError =
true;
315 if (!singleThreadModeError)
320 auto doInvMass = [
this, &getSet] ()
327 MPIManager :: run ( doInvMass );
331 singleThreadModeError =
true;
335 if (singleThreadModeError)
361 template <
class Entity,
class LocalMatrix>
363 LocalMatrix& lop )
const
365 jOp_.getLocalMatrix( domainEntity, rangeEntity, lop );
368 template <
class Entity,
class LocalMatrix>
370 const LocalMatrix& lop )
372 jOp_.setLocalMatrix( domainEntity, rangeEntity, lop );
384 template <
class Entity,
class LocalMatrix>
386 LocalMatrix& lop )
const
388 std::lock_guard< mutex_t > lock(
mutex_ );
392 template <
class Entity,
class LocalMatrix>
394 const LocalMatrix& lop )
396 std::lock_guard< mutex_t > lock(
mutex_ );
401 template <
class JacOp>
412 template <
class DomainFunction,
class RangeFunction>
413 struct GetSetLocalMatrix< PetscLinearOperator< DomainFunction, RangeFunction > > :
public GetSetLocalMatrixImplLocked
415 typedef GetSetLocalMatrixImplLocked
BaseType;
420 template <
class Iterators,
class GetSetLocalMatrixOp >
422 const bool hasSkeleton )
const
427 LocalMassMatrixType localMassMatrix( jOp.rangeSpace(),
impl().interiorQuadratureOrder( jOp.rangeSpace().order() ) );
430 lop(jOp.domainSpace(), jOp.rangeSpace());
433 for(
const auto& inside : iterators )
437 auto guard =
bindGuard( lop, inside,inside );
438 lop.
bind(inside,inside);
439 getSet.getLocalMatrix( inside, inside, lop );
441 getSet.setLocalMatrix( inside, inside, lop );
446 for(
const auto &intersection : intersections(
gridPart(), inside ) )
449 if( intersection.neighbor() )
451 const auto& outside = intersection.outside();
452 auto guard =
bindGuard( lop, outside,inside );
453 getSet.getLocalMatrix( outside, inside, lop );
455 getSet.setLocalMatrix( outside, inside, lop );
475 template<
class Integrands,
class DomainFunction,
class RangeFunction >
486 template<
class... Args >
497 template <
class LinearOperator,
class ModelIntegrands >
512 template<
class Gr
idFunction >
518 template<
class Gr
idFunction >
521 (*this).jacobian( u, jOp );
529 template<
class Integrands,
class LinearOperator,
class InverseOperator,
bool addDirichletBC>
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
Definition: molgalerkin.hh:357
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
Definition: molgalerkin.hh:377
std::shared_mutex mutex_t
Definition: molgalerkin.hh:379
mutex_t mutex_
Definition: molgalerkin.hh:380
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
Definition: molgalerkin.hh:403
GetSetLocalMatrixImpl BaseType
Definition: molgalerkin.hh:404
GetSetLocalMatrix(JacobianOperatorType &jOp)
Definition: molgalerkin.hh:405
Definition: molgalerkin.hh:479
BaseType::GridPartType GridPartType
Definition: molgalerkin.hh:484
MOLAutomaticDifferenceGalerkinOperator(const GridPartType &gridPart, Args &&... args)
Definition: molgalerkin.hh:487
Definition: molgalerkin.hh:500
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