diff --git a/dune/fem-dg/algorithm/evolution.hh b/dune/fem-dg/algorithm/evolution.hh
index 710889b5e5128039cea359cc1e4f51a2a81a59e0..c1f544e18302a7a8ab02d4eff8f75712c0ea611e 100644
--- a/dune/fem-dg/algorithm/evolution.hh
+++ b/dune/fem-dg/algorithm/evolution.hh
@@ -51,6 +51,7 @@ namespace Fem
{
protected:
const std::string keyPrefix_;
+ const Dune::Fem::ParameterReader parameter_;
public:
/**
@@ -58,8 +59,10 @@ namespace Fem
*
* \param keyPrefix the key prefix for the parameter file.
*/
- TimeSteppingParameters( const std::string keyPrefix = "femdg.stepper." )
- : keyPrefix_( keyPrefix )
+ TimeSteppingParameters( const std::string keyPrefix = "femdg.stepper.",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : keyPrefix_( keyPrefix ),
+ parameter_( parameter )
{}
/**
@@ -70,7 +73,7 @@ namespace Fem
*/
virtual double fixedTimeStep() const
{
- return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "fixedtimestep" , 0.0 );
+ return parameter_.getValue< double >( keyPrefix_ + "fixedtimestep" , 0.0 );
}
/**
@@ -81,7 +84,7 @@ namespace Fem
*/
virtual double fixedTimeStepEocLoopFactor() const
{
- return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "fixedtimestepeocloopfactor" , 1.0 );
+ return parameter_.getValue< double >( keyPrefix_ + "fixedtimestepeocloopfactor" , 1.0 );
}
/**
@@ -89,7 +92,7 @@ namespace Fem
*/
virtual double startTime() const
{
- return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "starttime" , 0.0 );
+ return parameter_.getValue< double >( keyPrefix_ + "starttime" , 0.0 );
}
/**
@@ -97,7 +100,7 @@ namespace Fem
*/
virtual double endTime() const
{
- return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "endtime"/*, 1.0 */);
+ return parameter_.getValue< double >( keyPrefix_ + "endtime"/*, 1.0 */);
}
/**
@@ -109,7 +112,7 @@ namespace Fem
*/
virtual int printCount() const
{
- return Dune::Fem::Parameter::getValue< int >( keyPrefix_ + "printcount" , -1 );
+ return parameter_.getValue< int >( keyPrefix_ + "printcount" , -1 );
}
/**
@@ -118,7 +121,7 @@ namespace Fem
*/
virtual double maxTimeStep() const
{
- return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "maxtimestep", std::numeric_limits<double>::max());
+ return parameter_.getValue< double >( keyPrefix_ + "maxtimestep", std::numeric_limits<double>::max());
}
/**
@@ -126,7 +129,7 @@ namespace Fem
*/
virtual int maximalTimeSteps () const
{
- return Dune::Fem::Parameter::getValue< int >( keyPrefix_ + "maximaltimesteps", std::numeric_limits<int>::max());
+ return parameter_.getValue< int >( keyPrefix_ + "maximaltimesteps", std::numeric_limits<int>::max());
}
/**
@@ -142,7 +145,7 @@ namespace Fem
*/
virtual bool stopAtEndTime() const
{
- return Dune::Fem::Parameter::getValue< bool >( keyPrefix_ + "stopatendtime", bool(false) );
+ return parameter_.getValue< bool >( keyPrefix_ + "stopatendtime", bool(false) );
}
};
diff --git a/dune/fem-dg/models/additional.hh b/dune/fem-dg/models/additional.hh
index 2ddba3a64041fe4a0b7a33cf152c5c88d03d38b7..44e8e744b66e4cb219e59d360a765d9025b57df7 100644
--- a/dune/fem-dg/models/additional.hh
+++ b/dune/fem-dg/models/additional.hh
@@ -145,17 +145,24 @@ struct Additional
{
{}
}
+ template< class Entity, class Point >
+ static void adjustAverageValue ( const Entity& entity, const Point &x, RangeType &u )
+ {
+ {}
+ }
static const int limitedDimRange = FunctionSpace :: dimRange;
static const bool hasAdvection = true;
static const bool hasDiffusion = false;
static const bool hasStiffSource = false;
static const bool hasNonStiffSource = false;
static const bool hasFlux = true;
+ static const bool threading = true;
static const Dune::Fem::Solver::Enum solverId = Dune::Fem::Solver::Enum::fem;
static const Dune::Fem::Formulation::Enum formId = Dune::Fem::Formulation::Enum::primal;
static const Dune::Fem::AdvectionLimiter::Enum limiterId = Dune::Fem::AdvectionLimiter::Enum::limited;
static const Dune::Fem::AdvectionFlux::Enum advFluxId = Dune::Fem::AdvectionFlux::Enum::llf;
static const Dune::Fem::DiffusionFlux::Enum diffFluxId = Dune::Fem::DiffusionFlux::Enum::none;
+ static const Dune::Fem::AdvectionLimiterFunction::Enum limiterFunctionId = Dune::Fem::AdvectionLimiterFunction::Enum::minmod;
};
// Model
diff --git a/dune/fem-dg/operator/adaptation/utility.hh b/dune/fem-dg/operator/adaptation/utility.hh
index cfc3a429e4036e93fcdf8754ed4e384505c09412..7d4fa9655fa3b6b14ece011500f693eeee57aea9 100644
--- a/dune/fem-dg/operator/adaptation/utility.hh
+++ b/dune/fem-dg/operator/adaptation/utility.hh
@@ -24,14 +24,20 @@ namespace Fem
protected:
const std::string keyPrefix_;
int markStrategy_;
-
-
+ const Dune::Fem::ParameterReader parameter_;
public:
//! constructor
- AdaptationParameters( const std::string keyPrefix = "fem.adaptation." )
+ AdaptationParameters( const std::string keyPrefix = "fem.adaptation.",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
: keyPrefix_( keyPrefix ),
- markStrategy_( getStrategy() )
+ markStrategy_( getStrategy() ),
+ parameter_(parameter)
+ {}
+ AdaptationParameters( const Dune::Fem::ParameterReader ¶meter )
+ : keyPrefix_( "fem.adaptation" ),
+ markStrategy_( getStrategy() ),
+ parameter_(parameter)
{}
////TODO:
diff --git a/dune/fem-dg/operator/dg/operatorbase.hh b/dune/fem-dg/operator/dg/operatorbase.hh
index 5cf4dcdc9a960a23882fe9343f3a4dfabe7fa79f..4033a0323b83feb37e30d0e2e98c305d71853056 100644
--- a/dune/fem-dg/operator/dg/operatorbase.hh
+++ b/dune/fem-dg/operator/dg/operatorbase.hh
@@ -142,9 +142,10 @@ namespace Fem
template< class ExtraParameterTupleImp >
DGAdvectionDiffusionOperatorBase( GridPartType& gridPart, const ModelType& model,
ExtraParameterTupleImp& tuple,
- const std::string name = "" )
+ const std::string name = "",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
: model_( model )
- , numflux_( model_ )
+ , numflux_( model_, parameter )
, gridPart_( gridPart )
, space_( gridPart_ )
, discreteModel_( model_, numflux_,
diff --git a/dune/fem-dg/operator/dg/primaloperator.hh b/dune/fem-dg/operator/dg/primaloperator.hh
index 5150e451405c67716a93cd6a26385bea07719a81..554f6769fda3412cd4be9ad4c5d3f791c9bba4a5 100644
--- a/dune/fem-dg/operator/dg/primaloperator.hh
+++ b/dune/fem-dg/operator/dg/primaloperator.hh
@@ -364,18 +364,19 @@ namespace Fem
DGLimitedAdvectionOperator( GridPartType& gridPart,
const ModelType& model,
ExtraParameterTupleImp tuple,
- const std::string name = "" )
+ const std::string name = "",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
: model_( model )
- , advflux_( model_ )
+ , advflux_( model_, parameter )
, gridPart_( gridPart )
, space_( gridPart_ )
, limiterSpace_( gridPart_ )
, uTmp_()
, fvSpc_()
, indicator_()
- , diffFlux_( gridPart_, model_, DGPrimalDiffusionFluxParameters( ParameterKey::generate( name, "dgdiffusionflux." ) ) )
+ , diffFlux_( gridPart_, model_, DGPrimalDiffusionFluxParameters( ParameterKey::generate( name, "dgdiffusionflux." ), parameter ) )
, discreteModel1_( model_, advflux_, diffFlux_ )
- , limitDiscreteModel_( model_ , space_.order() )
+ , limitDiscreteModel_( model_ , space_.order(), parameter )
, pass0_()
, pass1_( limitDiscreteModel_, pass0_, limiterSpace_ )
, pass2_( discreteModel1_, pass1_, space_ )
diff --git a/dune/fem-dg/operator/fluxes/advection/parameters.hh b/dune/fem-dg/operator/fluxes/advection/parameters.hh
index 011076986f80e2ec1442ca35d9959cc3fdd2063f..66622ffc4070961fdeeb9aaf112b550deaa12082 100644
--- a/dune/fem-dg/operator/fluxes/advection/parameters.hh
+++ b/dune/fem-dg/operator/fluxes/advection/parameters.hh
@@ -70,8 +70,14 @@ namespace Fem
*
* \param[in] keyPrefix key prefix for parameter file.
*/
- AdvectionFluxParameters( const std::string keyPrefix = "dgadvectionflux." )
- : keyPrefix_( keyPrefix )
+ AdvectionFluxParameters( const std::string keyPrefix = "dgadvectionflux.",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : keyPrefix_( keyPrefix ),
+ parameter_( parameter )
+ {}
+ AdvectionFluxParameters( const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : keyPrefix_( "dgadvectionflux." ),
+ parameter_( parameter )
{}
/**
@@ -94,12 +100,13 @@ namespace Fem
*/
virtual IdEnum getMethod() const
{
- const int i = Fem::Parameter::getEnum( keyPrefix_ + "method", AdvectionFlux::_strings );
+ const int i = parameter_.getEnum( keyPrefix_ + "method", AdvectionFlux::_strings );
return AdvectionFlux::_enums[i];
}
private:
const std::string keyPrefix_;
+ const Dune::Fem::ParameterReader parameter_;
};
diff --git a/dune/fem-dg/operator/fluxes/diffusion/parameters.hh b/dune/fem-dg/operator/fluxes/diffusion/parameters.hh
index 855a2e95ece64de8586e9713cc5146c6b27bfb56..06f8f24dbdf95fc81859da31595dfe7656bfebf0 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/parameters.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/parameters.hh
@@ -139,8 +139,10 @@ namespace Fem
*
* \param[in] keyPrefix key prefix for parameter file.
*/
- DGPrimalDiffusionFluxParameters( const std::string keyPrefix = "dgdiffusionflux." )
- : keyPrefix_( keyPrefix )
+ DGPrimalDiffusionFluxParameters( const std::string keyPrefix = "dgdiffusionflux.",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : keyPrefix_( keyPrefix ),
+ parameter_( parameter )
{}
/**
@@ -163,7 +165,7 @@ namespace Fem
*/
virtual IdEnum getMethod() const
{
- const int i = Fem::Parameter::getEnum( keyPrefix_ + "method", PrimalDiffusionFlux::_strings, 0 );
+ const int i = parameter_.getEnum( keyPrefix_ + "method", PrimalDiffusionFlux::_strings, 0 );
return PrimalDiffusionFlux::_enums[i];
}
@@ -187,20 +189,20 @@ namespace Fem
*/
virtual LiftingEnum getLifting() const
{
- const int i = Fem::Parameter::getEnum( keyPrefix_ + "lifting", PrimalDiffusionLifting::_strings, 0 );
+ const int i = parameter_.getEnum( keyPrefix_ + "lifting", PrimalDiffusionLifting::_strings, 0 );
return PrimalDiffusionLifting::_enums[i];
}
//! todo please doc me
virtual double penalty() const
{
- return Fem::Parameter::getValue<double>( keyPrefix_ + "penalty", 1.0 );
+ return parameter_.getValue<double>( keyPrefix_ + "penalty", 1.0 );
}
//! todo please doc me
virtual double liftfactor() const
{
- return Fem::Parameter::getValue<double>( keyPrefix_ + "liftfactor", 1.0 );
+ return parameter_.getValue<double>( keyPrefix_ + "liftfactor", 1.0 );
}
/**
@@ -208,19 +210,19 @@ namespace Fem
*/
virtual double theoryparameters() const
{
- return Fem::Parameter::getValue<double>( keyPrefix_ + "theoryparameters", 0. );
+ return parameter_.getValue<double>( keyPrefix_ + "theoryparameters", 0. );
}
//! todo please doc me
template <class DomainType>
void upwind( DomainType& upwd ) const
{
- Fem::Parameter::get(keyPrefix_ + "upwind", upwd, upwd);
+ parameter_.get(keyPrefix_ + "upwind", upwd, upwd);
}
private:
const std::string keyPrefix_;
-
+ const Dune::Fem::ParameterReader parameter_;
};
/**
diff --git a/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh b/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
index 9904cb89cf74caf7ed1214046cd121fc21df64e7..64511b6b97b4b1427293005c8a442ac2fce6ba26 100644
--- a/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
+++ b/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
@@ -72,6 +72,7 @@ namespace Fem
double crsTol_;
int finLevel_;
int crsLevel_;
+ const AdaptationParameters param;
const bool shockIndicatorAdaptivty_;
public:
@@ -81,13 +82,14 @@ namespace Fem
//! constructor
StandardLimiterDiscreteModel(const Model& mod,
const int polOrd,
- const AdaptationParameters& param = AdaptationParameters() )
- : BaseType(mod),
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : BaseType(mod, 1e-8, parameter),
indicator_( 0 ),
refTol_( -1 ),
crsTol_( -1 ),
finLevel_( 0 ),
crsLevel_( 0 ),
+ param(parameter),
shockIndicatorAdaptivty_( param.shockIndicator() )
{
if( shockIndicatorAdaptivty_ )
diff --git a/dune/fem-dg/operator/limiter/limiterutility.hh b/dune/fem-dg/operator/limiter/limiterutility.hh
index 1e692b9ff03bcccbf7c1823fbb69bcc8e7b6f64a..b313ec6531e481f6c053aa120108b80baa3f13f0 100644
--- a/dune/fem-dg/operator/limiter/limiterutility.hh
+++ b/dune/fem-dg/operator/limiter/limiterutility.hh
@@ -909,16 +909,16 @@ namespace Fem
struct LimiterFunctionBase
{
const double limitEps_;
- LimiterFunctionBase()
- : limitEps_( getEpsilon() )
+ LimiterFunctionBase( const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : limitEps_( getEpsilon(parameter) )
{
}
//! get tolerance for shock detector
- static double getEpsilon()
+ static double getEpsilon( const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
{
// default value is 1e-8
- return Parameter::getValue("femdg.limiter.limiteps", double(1e-8) );
+ return parameter.getValue("femdg.limiter.limiteps", double(1e-8) );
}
//! return epsilon for limiting
@@ -955,7 +955,8 @@ namespace Fem
using LimiterFunctionBase :: limitEps_ ;
typedef Field FieldType;
- MinModLimiter() : LimiterFunctionBase()
+ MinModLimiter( const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : LimiterFunctionBase(parameter)
{
printInfo( "minmod" );
}
@@ -983,7 +984,8 @@ namespace Fem
using LimiterFunctionBase :: limitEps_ ;
typedef Field FieldType;
- SuperBeeLimiter() : LimiterFunctionBase()
+ SuperBeeLimiter( const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : LimiterFunctionBase(parameter)
{
printInfo( "superbee" );
}
@@ -1012,7 +1014,8 @@ namespace Fem
using LimiterFunctionBase :: limitEps_ ;
typedef Field FieldType;
- VanLeerLimiter() : LimiterFunctionBase()
+ VanLeerLimiter( const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : LimiterFunctionBase(parameter)
{
printInfo( "vanLeer" );
}
diff --git a/dune/fem-dg/operator/limiter/limitpass.hh b/dune/fem-dg/operator/limiter/limitpass.hh
index abcd04a1baf2eb37f2f73e668ff3277a93be928c..7c2ea8618b326345f2bb7166898648d4519672c5 100644
--- a/dune/fem-dg/operator/limiter/limitpass.hh
+++ b/dune/fem-dg/operator/limiter/limitpass.hh
@@ -211,9 +211,10 @@ namespace Fem
/** \brief default limiter discrete model */
LimiterDefaultDiscreteModel(const Model& mod,
- const DomainFieldType veloEps = 1e-8)
+ const DomainFieldType veloEps = 1e-8,
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
: model_(mod)
- , limiterFunction_()
+ , limiterFunction_(parameter)
, velocity_(0)
, veloEps_( veloEps )
{
@@ -822,7 +823,7 @@ namespace Fem
ContextType cLeft( entity, intersection, faceQuadInner, volume_ );
// create quadrature of low order
- EvalType local( cLeft, entityValue, entityValue );
+ EvalType local( cLeft, entityValue, entityValue, 0 );
// check for boundary Value
if( discreteModel_.hasBoundaryValue( local ) )
diff --git a/dune/fem-dg/solver/dg.hh b/dune/fem-dg/solver/dg.hh
index be996224d99219490e8117a9e36642b247bb1164..637fbbf9fce0d447ad19a13949da8456d475d332 100644
--- a/dune/fem-dg/solver/dg.hh
+++ b/dune/fem-dg/solver/dg.hh
@@ -140,6 +140,41 @@ namespace Fem
std::cout << "cfl = " << double(tp_.factor()) << " " << tp_.time() << std::endl;
}
+ DGOperator( const DiscreteFunctionSpaceType& space,
+ const AdvectionModel &advectionModel,
+ const DiffusionModel &diffusionModel,
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : space_( space ),
+ extra_(),
+ tpPtr_( new TimeProviderType(space_.gridPart().comm(), parameter) ),
+ tp_( *tpPtr_ ),
+ model_( advectionModel, diffusionModel ),
+ fullOperator_( space.gridPart(), model_, extra_, name(), parameter ),
+ explOperator_( space.gridPart(), model_, extra_, name(), parameter ),
+ implOperator_( space.gridPart(), model_, extra_, name(), parameter ),
+ rkSolver_( tp_, fullOperator_, explOperator_, implOperator_, name(), parameter ),
+ initialized_( false )
+ {
+ //Dune::Fem::Parameter::append("fem.parallel.numberofthreads", std::to_string( Additional::nThreads ) );
+
+ const TimeSteppingParameters param("femdg.stepper.",parameter);
+ const double maxTimeStep = param.maxTimeStep();
+ fixedTimeStep_ = param.fixedTimeStep();
+
+ // start first time step with prescribed fixed time step
+ // if it is not 0 otherwise use the internal estimate
+ tp_.provideTimeStepEstimate(maxTimeStep);
+
+ // adjust fixed time step with timeprovider.factor()
+ fixedTimeStep_ /= tp_.factor() ;
+ if ( fixedTimeStep_ > 1e-20 )
+ tp_.init( fixedTimeStep_ );
+ else
+ tp_.init();
+
+ std::cout << "cfl = " << double(tp_.factor()) << " " << tp_.time() << std::endl;
+ }
+
DGOperator( TimeProviderType& tp,
const DiscreteFunctionSpaceType& space,
const AdvectionModel &advectionModel,
diff --git a/dune/fem-dg/solver/rungekuttasolver.hh b/dune/fem-dg/solver/rungekuttasolver.hh
index 878277c1a3eb9f2d5cf74e55fe33f982cfd6e9ac..bd9acd8b7a54f3f7c9fbe9845601eba6f566722f 100644
--- a/dune/fem-dg/solver/rungekuttasolver.hh
+++ b/dune/fem-dg/solver/rungekuttasolver.hh
@@ -29,11 +29,13 @@ namespace Fem
*
* \param[in] keyPrefix key prefix for parameter file.
*/
- RungeKuttaSolverParameters( const std::string keyPrefix = "fem.ode." )
- : DuneODE::ImplicitRungeKuttaSolverParameters( keyPrefix )
+ RungeKuttaSolverParameters( const std::string keyPrefix = "fem.ode.",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
+ : DuneODE::ImplicitRungeKuttaSolverParameters( keyPrefix, parameter )
{}
using DuneODE::ImplicitRungeKuttaSolverParameters::keyPrefix_;
+ using DuneODE::ImplicitRungeKuttaSolverParameters::parameter;
/**
* \brief returns a factor for the decision whether to use IMEX Runge-Kutta
@@ -50,7 +52,7 @@ namespace Fem
*/
virtual double explicitFactor() const
{
- return Fem::Parameter::getValue< double >( keyPrefix_ + "explicitfactor" , 1.0 );
+ return parameter().getValue< double >( keyPrefix_ + "explicitfactor" , 1.0 );
}
/**
@@ -82,7 +84,7 @@ namespace Fem
// methods
static const std::string odeSolver[] = { "EX", "IM" ,"IMEX", "IMEX+" };
std::string key( keyPrefix_ + "odesolver" );
- return Fem::Parameter::getEnum( key, odeSolver, 0 );
+ return parameter().getEnum( key, odeSolver, 0 );
}
/**
@@ -95,11 +97,10 @@ namespace Fem
{
std::string key( keyPrefix_ + "order");
if ( Fem :: Parameter :: exists( key ) )
- return Fem::Parameter::getValue< int > ( key );
+ return parameter().getValue< int > ( key );
else
return defaultRKOrder ;
}
-
};
@@ -162,7 +163,9 @@ namespace Fem
{
template < class OdeParameter >
static solverpair_t
- createExplicitSolver( Op& op, Fem::TimeProviderBase& tp, const int rkSteps, const OdeParameter& param, const std::string& name = "" )
+ createExplicitSolver( Op& op, Fem::TimeProviderBase& tp, const int rkSteps, const OdeParameter& param,
+ const Dune::Fem::ParameterReader ¶meter,
+ const std::string& name = "" )
{
typedef DuneODE :: ExplicitRungeKuttaSolver< DiscreteFunctionType > ExplicitOdeSolverType;
return solverpair_t( new ExplicitOdeSolverType( op, tp, rkSteps ), nullptr );
@@ -170,7 +173,9 @@ namespace Fem
template < class OdeParameter >
static solverpair_t
- createImplicitSolver( Op& op, Fem::TimeProviderBase& tp, const int rkSteps, const OdeParameter& param, const std::string& name = "" )
+ createImplicitSolver( Op& op, Fem::TimeProviderBase& tp, const int rkSteps, const OdeParameter& param,
+ const Dune::Fem::ParameterReader ¶meter,
+ const std::string& name = "" )
{
#ifdef COUNT_FLOPS
return solverpair_t();
@@ -187,14 +192,16 @@ namespace Fem
typedef typename NonlinearInverseOperatorType::ParametersType NonlinParametersType;
- return solverpair_t(new ImplicitOdeSolverType( *helmOp, tp, rkSteps, param, NonlinParametersType( ParameterKey::generate( name, "fem.solver.newton." ) ) ), helmOp );
+ return solverpair_t(new ImplicitOdeSolverType( *helmOp, tp, rkSteps, param, NonlinParametersType( ParameterKey::generate( name, "fem.solver.newton." ), parameter ) ), helmOp );
#endif
}
template < class ExplOp, class ImplOp, class OdeParameter >
static solverpair_t
createSemiImplicitSolver( ExplOp& explOp, ImplOp& implOp, Fem::TimeProviderBase& tp,
- const int rkSteps, const OdeParameter& param, const std::string& name = "" )
+ const int rkSteps, const OdeParameter& param,
+ const Dune::Fem::ParameterReader ¶meter,
+ const std::string& name = "" )
{
#ifdef COUNT_FLOPS
return solverpair_t();
@@ -213,7 +220,7 @@ namespace Fem
typedef typename NonlinearInverseOperatorType::ParametersType NonlinParametersType;
- return solverpair_t(new SemiImplicitOdeSolverType( explOp, *helmOp, tp, rkSteps, param, NonlinParametersType( ParameterKey::generate( name, "fem.solver.newton." ) ) ), helmOp );
+ return solverpair_t(new SemiImplicitOdeSolverType( explOp, *helmOp, tp, rkSteps, param, NonlinParametersType( ParameterKey::generate( name, "fem.solver.newton." ), parameter ) ), helmOp );
#endif
}
@@ -248,13 +255,14 @@ namespace Fem
OperatorType& op,
ExplicitOperatorType& advOp,
ImplicitOperatorType& diffOp,
- const std::string name = "" )
+ const std::string name = "",
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
: operator_( op ),
timeProvider_( tp ),
name_( name ),
explicitOperator_( advOp ),
implicitOperator_( diffOp ),
- param_( new RungeKuttaSolverParameters( ParameterKey::generate( name_, "fem.ode." ) ) ),
+ param_( new RungeKuttaSolverParameters( ParameterKey::generate( name_, "fem.ode." ), parameter ) ),
explicitSolver_( 0 ),
odeSolver_( 0 ),
helmholtzOperator_( 0 ),
@@ -271,11 +279,11 @@ namespace Fem
// create implicit or explicit ode solver
if( odeSolverType_ == 0 )
{
- solver = OdeSolversType :: createExplicitSolver( operator_, tp, rkSteps_, *param_, name_ );
+ solver = OdeSolversType :: createExplicitSolver( operator_, tp, rkSteps_, *param_, parameter, name_ );
}
else if (odeSolverType_ == 1)
{
- solver = OdeSolversType :: createImplicitSolver( operator_, tp, rkSteps_, *param_, name_ );
+ solver = OdeSolversType :: createImplicitSolver( operator_, tp, rkSteps_, *param_, parameter, name_ );
}
else if( odeSolverType_ > 1 )
{
@@ -285,11 +293,11 @@ namespace Fem
DUNE_THROW(Dune::InvalidStateException,"Advection and Diffusion operator are the same, therefore IMEX cannot work!");
}
- solver = OdeSolversType :: createSemiImplicitSolver( explicitOperator_, implicitOperator_, tp, rkSteps_, *param_, name_ );
+ solver = OdeSolversType :: createSemiImplicitSolver( explicitOperator_, implicitOperator_, tp, rkSteps_, *param_, parameter, name_ );
// IMEX+
if( odeSolverType_ == 3 )
- explicitSolver_ = OdeSolversType :: createExplicitSolver( operator_, tp, rkSteps_, *param_, name_ ).first;
+ explicitSolver_ = OdeSolversType :: createExplicitSolver( operator_, tp, rkSteps_, *param_, parameter, name_ ).first;
}
else
{
diff --git a/pydemo/euler/euler.py b/pydemo/euler/euler.py
index 792f420bc0fa2103791fde59f2b03baa391dfc1e..42f39309fb923b73c808bea3450dbcd38b5414ac 100644
--- a/pydemo/euler/euler.py
+++ b/pydemo/euler/euler.py
@@ -46,20 +46,20 @@ def CompressibleEuler(dim, gamma):
return (pL - pR)/(0.5*(pL + pR))
return Model
-def CompressibleEulerNeuman(dim, gamma):
+def CompressibleEulerNeuman(dim, gamma, bnd=range(1,5)):
class Model(CompressibleEuler(dim,gamma)):
- boundary = {range(1,5): lambda t,x,u,n: Model.F_c(t,x,u)*n}
+ boundary = {bnd: lambda t,x,u,n: Model.F_c(t,x,u)*n}
return Model
-def CompressibleEulerDirichlet(dim, gamma):
+def CompressibleEulerDirichlet(dim, gamma, bnd=range(1,5)):
class Model(CompressibleEuler(dim,gamma)):
- boundary = {range(1,5): lambda t,x,u: u}
+ boundary = {bnd: lambda t,x,u: u}
return Model
-def CompressibleEulerSlip(dim, gamma):
+def CompressibleEulerSlip(dim, gamma,bnd=range(1,5)):
class Model(CompressibleEuler(dim,gamma)):
def outflowFlux(t,x,u,n):
_,_, p = CompressibleEuler(dim,gamma).toPrim(u)
return as_vector([ 0, *(p*n), 0 ])
- boundary = {range(1,5): outflowFlux}
+ boundary = {bnd: outflowFlux}
return Model
def riemanProblem(Model,x,x0,UL,UR):
@@ -138,7 +138,8 @@ def vortex(dim=2,gamma=1.4):
u = S*x[1]*exp(f)/(2.*pi*R)
v = 1. - S*x[0]*exp(f)/(2.*pi*R)
p = rho / (gamma*M*M)
- Model = CompressibleEulerDirichlet(dim,gamma)
+ # Model = CompressibleEuler(dim,gamma)
+ Model = CompressibleEulerSlip(dim,gamma,bnd=(1,2))
return Model,\
Model.toCons( as_vector( [rho,u,v,p] )),\
[-10, -10], [10, 10], [20, 20], 100,\
diff --git a/pydemo/euler/euler.ufl b/pydemo/euler/euler.ufl
deleted file mode 100644
index defd0b142daa183763276901ac03fc5c27d69a76..0000000000000000000000000000000000000000
--- a/pydemo/euler/euler.ufl
+++ /dev/null
@@ -1,46 +0,0 @@
-gamma = 1.4
-class Compressible2DEuler:
- def pressure(U):
- return (U[3]-(U[1]**2+U[2]**2)/2)*(gamma-1)
- def toCons(U):
- return [U[0],U[0]*U[1],U[0]*U[2],U[3]/(gamma-1)+U[0]/2*(U[1]**2+U[2]**2)]
- def toPrim(U):
- return [U[0],U[1]/U[0],U[2]/U[0],Compressible2DEuler.pressure(U)]
- def F_c(U):
- rho, u,v, p = Compressible2DEuler.toPrim(U)
- rE = U[-1]
- res = as_matrix( [ [rho*u, rho*v],
- [rho*u**2 + p, rho*u*v],
- [rho*u*v, rho*v**2 + p],
- [(rE+p)*u, (rE+p)*v]
- ])
- return res
- def alpha(U,n):
- rho, u,v, p = Compressible2DEuler.toPrim(U)
- return abs(u*n[0]+v*n[1]) + sqrt(gamma*p/rho)
-class LLFFlux:
- def flux(Descr,U,n,dt):
- return jump(U)
- UL = U('-')
- UR = U('+')
- return UL-UR
- FL = Descr.F_c(UL)*n
- FR = Descr.F_c(UR)*n
- flux = (FL+FR)
- max_value = lambda a,b: (a+b+abs(a-b))/2.
- visc = max_value( Descr.alpha(UL,n), Descr.alpha(UR,n) )
- return flux + visc/dt*(UR-UL)
-
-Model = Compressible2DEuler
-NumFLux = LLFFlux
-
-space = Space(2,4)
-u = TrialFunction(space)
-v = TestFunction(space)
-dt = NamedConstant(triangle, "dt") # time step
-
-# just a starting point...
-def femdgModel(description):
- return inner(description.F_c(u),grad(v))*dx
-
-F = -dt*femdgModel(Model)
diff --git a/pydemo/euler/euler.ulf.old b/pydemo/euler/euler.ulf.old
deleted file mode 100644
index 128f055a0a9d3c8f88ee3d868ec82e6b9796e381..0000000000000000000000000000000000000000
--- a/pydemo/euler/euler.ulf.old
+++ /dev/null
@@ -1,53 +0,0 @@
-class Compressible2DEuler:
- def rhoeps(U):
- v = U[1]*U[1]
- for i in range(dim-1): v += U[i+2]*U[i+2]
- v *= 0.5/U[0]
- rE = U[dim]
- return rE - v
- def pressure(U):
- return (gamma-1)*Compressible2DEuler.rhoeps(U)
- def velo(U):
- return [U[i] for i in range(1,dim+1)]
- def toPrim(U):
- return [U[0], Compressible2DEuler.velo(U), U[dim]]
- def F_c(U):
- rho, v, p = Compressible2DEuler.toPrim(U)
- rE = U[dim]
- res = [ rho*v,
- [rho*v[0]*v[0] + p, rho*v[0]*v[1]],
- [rho*v[0]*v[1], rho*v[1]*v[1] + p],
- (rE+p)*v
- ]
- def maxLambda(U,n):
- rho, v, p = Compressible2DEuler.toPrim(U)
- return abs(dot(v,n)) + sqrt(gamma*p/rho)
- def velocity(U):
- return Compressible2DEuler.velo(U)
- def physical(U):
- return conditional( (U[0]>1e-8), conditional( Compressible2DEuler.rhoeps(U) > 1e-8 , 1, 0 ), 0 )
- def jump(U,V):
- pL = Compressible2DEuler.pressure(U)
- pR = Compressible2DEuler.pressure(V)
- return (pL - pR)/(0.5*(pL + pR))
-
- def outflowFlux(u,n):
- return Compressible2DEuler.F_c(u)*n
- def outflowValue(u):
- return u
-
- # boundaryFlux = {1: outflowFlux}
- boundaryValue = {1: outflowValue}
-
-Model = Compressible2DEuler
-
-from dune.ufl import Space, NamedConstant
-space = Space(2,4)
-u = TrialFunction(space)
-v = TestFunction(space)
-
-# just a starting point...
-def femdgModel(description):
- return inner(description.F_c(u),grad(v))*dx
-
-F = femdgModel(Model)
diff --git a/pydemo/euler/main.py b/pydemo/euler/main.py
index 2a9fe2c88c3b3d92d242b2cc17b1d985473b6606..a66ab8606054996aca79fbf6959a11289e6825cc 100644
--- a/pydemo/euler/main.py
+++ b/pydemo/euler/main.py
@@ -1,3 +1,5 @@
+import mpi4py.rc
+mpi4py.rc.threaded = True
from dune.fem import parameter
from dune.femdg.testing import run
@@ -10,12 +12,22 @@ dim = 2
gamma = 1.4
parameter.append({"fem.verboserank": -1})
-parameter.append({"fem.parallel.numberofthreads": 4})
-parameter.append("parameter")
-primitive=lambda Model,uh: {"pressure":Model.toPrim(uh)[2]}
+primitive=lambda Model,uh: {"pressure": Model.toPrim(uh)[2]}
+parameters = {"fem.ode.odesolver": "EX",
+ "fem.timeprovider.factor": 0.45,
+ "dgadvectionflux.method": "EULER-HLLC",
+ "femdg.limiter.limiteps": 1,
+ "femdg.limiter.admissiblefunctions": 1,
+ "femdg.limiter.tolerance": 1}
+#-----------------
+# femdg.limiter.admissiblefunctions:
+# 0 = only dg solution | 1 = only reconstruction | 2 = both
+#-----------------
+
run(*problem(),
- startLevel=0, polOrder=2, limiter="default",
- primitive=primitive, saveStep=0.1, subsamp=2,
- dt=None)
+ startLevel=0, polOrder=2, limiter="default",
+ primitive=primitive, saveStep=0.1, subsamp=2,
+ dt=None,threading=True,grid="alucube",
+ parameters=parameters)
diff --git a/pydemo/euler/paramBase b/pydemo/euler/paramBase
deleted file mode 100644
index 5b901d39786de0382089a6b317687431116b5152..0000000000000000000000000000000000000000
--- a/pydemo/euler/paramBase
+++ /dev/null
@@ -1,15 +0,0 @@
-# OMP THREADS
-#------------
-fem.verboserank: 0
-# number of threads used in OMP program
-fem.parallel.numberofthreads: 1
-# write diagnostics file (
-# 0 don't, 1 only speedup file, 2 write all runfiles
-# 3 only write 0, others at end, 4 all files at end for scaling)
-fem.parallel.diagnostics: 1
-# if true non-blocking communication is enabled
-femdg.nonblockingcomm: true
-fem.threads.verbose: true
-# if true thread 0 only communicates data but does not computation
-fem.threads.communicationthread: false
-fem.threads.partitioningmethod: sfc
diff --git a/pydemo/euler/paramSolver b/pydemo/euler/paramSolver
deleted file mode 100644
index fed5f13c3b14ca190d87aeac98a9e896fd9a0a83..0000000000000000000000000000000000000000
--- a/pydemo/euler/paramSolver
+++ /dev/null
@@ -1,27 +0,0 @@
-# SOLVER CONFIGURATION
-#---------------------
-
-fem.ode.odesolver: EX # ode solvers: EX, IM, IMEX
-# fem.ode.order: 3
-fem.ode.verbose: cfl # ode output: none, cfl, full
-fem.ode.cflincrease: 1.25
-fem.ode.miniterations: 95
-fem.ode.maxiterations: 105
-fem.ode.cflStart: 1.
-#fem.ode.cflMax: 5
-fem.timeprovider.factor: 0.45
-fem.timeprovider.updatestep: 1
-femdg.stepper.maxtimestep: 1e-3
-
-# parameter for the implicit solvers
-fem.solver.verbose: false
-fem.solver.gmres.restart: 15
-fem.solver.newton.verbose: false
-fem.solver.newton.linear.verbose: false
-fem.solver.newton.maxlineariterations: 1000
-fem.solver.newton.tolerance: 1e-10
-
-dgdiffusionflux.method: CDG2 # diffusion flux: CDG2, CDG, BR2, IP, NIPG, BO
-dgdiffusionflux.theoryparameters: 1 # scaling with theory parameters
-dgdiffusionflux.penalty: 0.
-dgdiffusionflux.liftfactor: 1.0
diff --git a/pydemo/euler/parameter b/pydemo/euler/parameter
deleted file mode 100644
index e781c1b166139f6859c839fdd7cb783b65cd037a..0000000000000000000000000000000000000000
--- a/pydemo/euler/parameter
+++ /dev/null
@@ -1,21 +0,0 @@
-# LIMITER SETTINGS
-#-----------------
-# 0 = only dg solution | 1 = only reconstruction | 2 = both
-femdg.limiter.admissiblefunctions: 1
-# tolerance for shock indicator
-femdg.limiter.tolerance: 1 # 1e-10
-# threshold for avoiding over-excessive limitation
-femdg.limiter.limiteps: 1e-8
-
-# GENERAL
-#--------
-paramfile: paramBase
-
-
-# choises are: LLF, HLL, HLLC, LLF2
-dgadvectionflux.method: LLF
-
-# SOLVER CONFIGURATION
-#---------------------
-fem.ode.odesolver: EX
-paramfile: paramSolver
diff --git a/pydemo/fulltest.py b/pydemo/fulltest.py
index 3d03e41ef30a9554f74c7a99603de996d5ea4317..4f50d7117b55701dc3b9a9c859c184f1f4e32fa1 100644
--- a/pydemo/fulltest.py
+++ b/pydemo/fulltest.py
@@ -9,3 +9,7 @@ parameter.append("parameter")
for p in shallowWater.problems + euler.problems:
run(*p(), startLevel=0, polOrder=2, limiter="default",
primitive=None, saveStep=None, subsamp=0)
+ dt=None,threading=True,grid="alusimplex")
+ run(*p(), startLevel=0, polOrder=2, limiter="default",
+ primitive=None, saveStep=None, subsamp=0)
+ dt=None,threading=True,grid="yasp")
diff --git a/pydemo/paramBase b/pydemo/paramThreading
similarity index 96%
rename from pydemo/paramBase
rename to pydemo/paramThreading
index 5b901d39786de0382089a6b317687431116b5152..c298d0adfa10d8ad2ceedaec44c6d9cd215ae61b 100644
--- a/pydemo/paramBase
+++ b/pydemo/paramThreading
@@ -1,6 +1,5 @@
# OMP THREADS
#------------
-fem.verboserank: 0
# number of threads used in OMP program
fem.parallel.numberofthreads: 1
# write diagnostics file (
diff --git a/pydemo/scalar/main.py b/pydemo/scalar/main.py
index f96632c76cbd174749d840e4d2b0a9658f87bc20..be0fc3fddb0887c72bb8f738ed436af886f5fdb6 100644
--- a/pydemo/scalar/main.py
+++ b/pydemo/scalar/main.py
@@ -3,13 +3,28 @@ from dune.femdg.testing import run
# from scalar import shockTransport as problem
# from scalar import sinProblem as problem
-from scalar import sinTransportProblem as problem
+# from scalar import sinTransportProblem as problem
# from scalar import pulse as problem
-# from scalar import diffusivePulse as problem
+from scalar import diffusivePulse as problem
parameter.append({"fem.verboserank": 0})
-parameter.append("parameter")
+
+parameters = {"fem.ode.odesolver": "IM", # EX, IM, IMEX
+ "fem.ode.order": 3,
+ "fem.ode.verbose": "none", # none, cfl, full
+ "fem.ode.cflincrease": 1.25,
+ "fem.ode.miniterations": 35,
+ "fem.ode.maxiterations": 100,
+ "fem.timeprovider.factor": 0.45,
+ "dgadvectionflux.method": "LLF",
+ "fem.solver.gmres.restart": 50,
+ "dgdiffusionflux.method": "CDG2", # CDG2, CDG, BR2, IP, NIPG, BO
+ "dgdiffusionflux.theoryparameters": 1, # scaling with theory parameters
+ "dgdiffusionflux.penalty": 0,
+ "dgdiffusionflux.liftfactor": 1}
run(*problem(),
startLevel=0, polOrder=2, limiter=None,
- primitive=None, saveStep=0.1, subsamp=0)
+ primitive=None, saveStep=0.01, subsamp=0,
+ dt=0.1,
+ parameters=parameters)
diff --git a/pydemo/scalar/paramBase b/pydemo/scalar/paramBase
deleted file mode 100644
index cba313dd29c863cffd04cb03c51fc44214167e79..0000000000000000000000000000000000000000
--- a/pydemo/scalar/paramBase
+++ /dev/null
@@ -1,15 +0,0 @@
-# OMP THREADS
-#------------
-fem.verboserank: 0
-# number of threads used in OMP program
-fem.parallel.numberofthreads: 4
-# write diagnostics file (
-# 0 don't, 1 only speedup file, 2 write all runfiles
-# 3 only write 0, others at end, 4 all files at end for scaling)
-fem.parallel.diagnostics: 1
-# if true non-blocking communication is enabled
-femdg.nonblockingcomm: true
-fem.threads.verbose: true
-# if true thread 0 only communicates data but does not computation
-fem.threads.communicationthread: false
-fem.threads.partitioningmethod: sfc
diff --git a/pydemo/scalar/paramSolver b/pydemo/scalar/paramSolver
deleted file mode 100644
index a87d3cb261e6ae9cbfb6e5982f5c0d1048b57306..0000000000000000000000000000000000000000
--- a/pydemo/scalar/paramSolver
+++ /dev/null
@@ -1,31 +0,0 @@
-# SOLVER CONFIGURATION
-#---------------------
-
-fem.ode.odesolver: IMEX # ode solvers: EX, IM, IMEX
-fem.ode.order: 3
-fem.ode.verbose: cfl # ode output: none, cfl, full
-fem.ode.cflincrease: 1.25
-fem.ode.miniterations: 35
-fem.ode.maxiterations: 100
-fem.ode.cflStart: 1.
-fem.ode.cflMax: 5
-fem.timeprovider.factor: 0.45
-fem.timeprovider.updatestep: 1
-
-# parameter for the implicit solvers
-# fem.differenceoperator.eps: 1e-12
-fem.solver.gmres.restart: 50
-fem.solver.newton.verbose: false
-fem.solver.newton.linear.verbose: false
-fem.solver.verbose: false
-fem.solver.newton.maxlineariterations: 100000
-fem.solver.newton.tolerance: 1e-07
-fem.solver.newton.linabstol: 2e-8
-fem.solver.newton.linreduction: 2e-8
-
-
-
-dgdiffusionflux.method: CDG2 # diffusion flux: CDG2, CDG, BR2, IP, NIPG, BO
-dgdiffusionflux.theoryparameters: 1 # scaling with theory parameters
-dgdiffusionflux.penalty: 0.
-dgdiffusionflux.liftfactor: 1.0
diff --git a/pydemo/scalar/parameter b/pydemo/scalar/parameter
deleted file mode 100644
index 3f04925934276d5319a83d2b99360d5c9b76db8d..0000000000000000000000000000000000000000
--- a/pydemo/scalar/parameter
+++ /dev/null
@@ -1,20 +0,0 @@
-# LIMITER SETTINGS
-#-----------------
-# 0 = only dg solution | 1 = only reconstruction | 2 = both
-femdg.limiter.admissiblefunctions: 1
-# tolerance for shock indicator
-femdg.limiter.tolerance: 1
-# threshold for avoiding over-excessive limitation
-femdg.limiter.limiteps: 1e-8
-
-# GENERAL
-#--------
-paramfile: paramBase
-
-
-# choises are: LLF, HLL, HLLC, LLF2
-dgadvectionflux.method: LLF
-
-# SOLVER CONFIGURATION
-#---------------------
-paramfile: paramSolver
diff --git a/pydemo/scalar/scalar.py b/pydemo/scalar/scalar.py
index 4f312c66113f9ff2c6ec7bb1c9f63b4f66594e5e..68c3ed80f106d71e3ea3e3cd372a7976e0c910af 100644
--- a/pydemo/scalar/scalar.py
+++ b/pydemo/scalar/scalar.py
@@ -11,7 +11,7 @@ class Transport1D:
def maxLambda(t,x,U,n):
return abs(n[0])
def velocity(t,x,U):
- return v
+ return as_vector([1,0])
def physical(U):
return 1
def jump(U,V):
diff --git a/pydemo/shallowWater/main.py b/pydemo/shallowWater/main.py
index 8483451bca820598141235ead1e0da03e8db1db2..d76b7c0c3e411cf42f5f0eb61cc5d4b9cd22b440 100644
--- a/pydemo/shallowWater/main.py
+++ b/pydemo/shallowWater/main.py
@@ -7,7 +7,14 @@ parameter.append({"fem.verboserank": -1})
parameter.append("parameter")
primitive=lambda Model,uh: {"freesurface":Model.toPrim(uh)[0]}
-
+parameters = {"fem.ode.odesolver": "EX",
+ "fem.timeprovider.factor": 0.45/2.,
+ "dgadvectionflux.method": "LLF",
+ "femdg.limiter.limiteps": 1,
+ "femdg.limiter.admissiblefunctions": 1,
+ "femdg.limiter.tolerance": 1}
run(*problem(1),
startLevel=0, polOrder=2, limiter="default",
- primitive=primitive, saveStep=0.1, subsamp=0)
+ primitive=primitive, saveStep=0.1, subsamp=0,
+ dt=None,threading=True, grid="alucube",
+ parameters=parameters)
diff --git a/python/dune/femdg/__init__.py b/python/dune/femdg/__init__.py
index e89deb8549b3e45529533633730a018b59dfbb7a..fe166b020858b8c497cecf3d191591af89fdc11d 100644
--- a/python/dune/femdg/__init__.py
+++ b/python/dune/femdg/__init__.py
@@ -1 +1,7 @@
from ._operators import *
+
+registry = {}
+
+registry["scheme"] = {
+ "femdg": femDGOperator
+ }
diff --git a/python/dune/femdg/_operators.py b/python/dune/femdg/_operators.py
index f17e0c5bb87f73db52379d16f63d6320b3696941..40d002d8a478962965901674db3beaf91912f33c 100644
--- a/python/dune/femdg/_operators.py
+++ b/python/dune/femdg/_operators.py
@@ -153,8 +153,9 @@ def generateMethod(struct,expr, cppType, name,
# create DG operator + solver (limiter = none,minmod,vanleer,superbee),
# (diffusionScheme = cdg2,br2,ip,nipg,...)
-def createFemDGSolver(Model, space,
- limiter="minmod", diffusionScheme = "cdg2", threading=False ):
+def femDGOperator(Model, space,
+ limiter="minmod", diffusionScheme = "cdg2", threading=False,
+ parameters={}):
import dune.create as create
if limiter is None or limiter is False:
@@ -480,12 +481,14 @@ def createFemDGSolver(Model, space,
constructor = Constructor(['const '+spaceType + ' &space',
'const '+advModelType + ' &advectionModel',
- 'const '+diffModelType + ' &diffusionModel'
+ 'const '+diffModelType + ' &diffusionModel',
+ 'const pybind11::dict ¶meters'
],
- ['return new DuneType(space, advectionModel, diffusionModel);'],
+ ['return new DuneType(space, advectionModel, diffusionModel, Dune::FemPy::pyParameter( parameters, std::make_shared< std::string >() ) );'],
['"space"_a',
'"advectionModel"_a',
'"diffusionModel"_a',
+ '"parameters"_a',
'pybind11::keep_alive< 1, 2 >()',
'pybind11::keep_alive< 1, 3 >()',
'pybind11::keep_alive< 1, 4 >()'])
@@ -499,8 +502,8 @@ def createFemDGSolver(Model, space,
# add method to set a fixed time step
setTimeStepSize = Method('setTimeStepSize', '&DuneType::setTimeStepSize')
# add method to solve (not requiring u_h_n)
- solve = Method('solve', '&DuneType::solve', extra=['"target"_a'])
+ solve = Method('step', '&DuneType::solve', extra=['"target"_a'])
return load(includes, typeName, constructor, setTimeStepSize, deltaT, applyLimiter, solve,
preamble=writer.writer.getvalue()).\
- Operator( space, advModel, diffModel )
+ Operator( space, advModel, diffModel, parameters=parameters )
diff --git a/python/dune/femdg/testing.py b/python/dune/femdg/testing.py
index 3e641f547df4bff2cfd9f605e5e1dd654678aae1..c6ada045786dc1751c869cac1349c2492f2c939e 100644
--- a/python/dune/femdg/testing.py
+++ b/python/dune/femdg/testing.py
@@ -2,16 +2,28 @@ import time, math
from dune.grid import structuredGrid, cartesianDomain, OutputType
import dune.create as create
from dune.fem.function import integrate
-from dune.femdg import createFemDGSolver
from dune.ufl import NamedConstant
from ufl import dot, SpatialCoordinate
def run(Model, initial, x0,x1,N, endTime, name, exact,
polOrder, limiter="default", startLevel=0,
primitive=None, saveStep=None, subsamp=0,
- dt=None):
- domain = cartesianDomain(x0,x1,N,periodic=[False,False])
- grid = create.grid("alucube",domain)
+ dt=None,grid="yasp",threading=True,
+ parameters={}):
+ periodic=[True,]*len(x0)
+ if hasattr(Model,"boundary"):
+ bnd=set()
+ for b in Model.boundary:
+ bnd.update(b)
+ for i in range(len(x0)):
+ if 2*i+1 in bnd:
+ assert(2*i+2 in bnd)
+ periodic[i] = False
+ if any(bnd):
+ print("Setting periodic boundaries",periodic,flush=True)
+
+ domain = cartesianDomain(x0,x1,N,periodic=periodic,overlap=0)
+ grid = create.grid(grid,domain)
grid.hierarchicalGrid.globalRefine(startLevel)
dimR = Model.dimension
t = 0
@@ -20,7 +32,7 @@ def run(Model, initial, x0,x1,N, endTime, name, exact,
space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
u_h = space.interpolate(initial, name='u_h')
- operator = createFemDGSolver( Model, space, limiter=limiter, threading=True )
+ operator = create.scheme("femDG",Model, space, limiter=limiter, threading=True, parameters=parameters )
operator.applyLimiter( u_h );
print("number of elements: ",grid.size(0),flush=True)
if saveStep is not None:
@@ -45,11 +57,11 @@ def run(Model, initial, x0,x1,N, endTime, name, exact,
while t < endTime:
if dt is not None:
operator.setTimeStepSize(dt)
- operator.solve(target=u_h)
+ operator.step(target=u_h)
dt = operator.deltaT()
if math.isnan( u_h.scalarProductDofs( u_h ) ):
grid.writeVTK(name, subsampling=subsamp, celldata=[u_h])
- print('ERROR: dofs invalid t =', t)
+ print('ERROR: dofs invalid t =', t,flush=True)
print('[',tcount,']','dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
exit(0)
t += dt
@@ -65,20 +77,20 @@ def run(Model, initial, x0,x1,N, endTime, name, exact,
vtk.write(name, count, outputType=OutputType.appendedraw)
saveTime += saveStep
- print("time loop:",time.time()-start)
- print("number of time steps ", tcount)
+ print("time loop:",time.time()-start,flush=True)
+ print("number of time steps ", tcount,flush=True)
if saveStep is not None:
try:
velo[0].setConstant("time",[t])
except:
pass
vtk.write(name, count, outputType=OutputType.appendedraw)
- else:
- if exact is not None:
- grid.writeVTK(name, subsampling=subsamp, celldata=[u_h, exact(t)])
- else:
- grid.writeVTK(name, subsampling=subsamp, celldata=[u_h])
+ # output the final result and compute error (if exact is available)
if exact is not None:
+ grid.writeVTK(name, subsampling=subsamp,
+ celldata={"solution":u_h, "exact":exact(t)})
error = integrate( grid, dot(u_h-exact(t),u_h-exact(t)), order=5 )
- print("error:", math.sqrt(error) )
+ print("error:", math.sqrt(error),flush=True )
+ else:
+ grid.writeVTK(name, subsampling=subsamp, celldata=[u_h])