Skip to content
Snippets Groups Projects
Commit a72f9b13 authored by Stefan Girke's avatar Stefan Girke
Browse files

cleaning up, operatortraits and upwindflux into single file

parent b0030b6a
No related branches found
No related tags found
No related merge requests found
......@@ -24,16 +24,16 @@ namespace Fem
typedef EvolutionAlgorithm< GridImp, ProblemTraits, polynomialOrder > BaseType ;
// type of Grid
typedef typename BaseType :: GridType GridType;
typedef typename BaseType :: GridType GridType;
// Choose a suitable GridView
typedef typename BaseType :: GridPartType GridPartType;
typedef typename BaseType :: GridPartType GridPartType;
// initial data type
typedef typename BaseType :: ProblemType ProblemType;
typedef typename BaseType :: ProblemType ProblemType;
// An analytical version of our model
typedef typename BaseType :: ModelType ModelType;
typedef typename BaseType :: ModelType ModelType;
// The DG space operator
// The first operator is sum of the other two
......@@ -51,17 +51,17 @@ namespace Fem
typedef typename BaseType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
// The ODE Solvers
typedef typename BaseType :: OdeSolverType OdeSolverType;
typedef typename BaseType :: OdeSolverType OdeSolverType;
typedef typename BaseType :: TimeProviderType TimeProviderType;
typedef typename BaseType :: TimeProviderType TimeProviderType;
typedef typename BaseType::OperatorTraits OperatorTraits;
typedef typename BaseType::OperatorTraits OperatorTraits;
typedef typename std::remove_pointer<typename std::tuple_element<0,typename BaseType::IndicatorTupleType>::type>::type DGIndicatorType;
typedef typename std::remove_pointer<typename std::tuple_element<1,typename BaseType::IndicatorTupleType>::type>::type GradientIndicatorType ;
// type of 64bit unsigned integer
typedef typename BaseType :: UInt64Type UInt64Type;
typedef typename BaseType :: UInt64Type UInt64Type;
typedef typename OperatorTraits :: ExtraParameterTupleType ExtraParameterTupleType;
......
......@@ -24,41 +24,41 @@ namespace Fem
typedef EvolutionAlgorithm< GridImp, ProblemTraits, polynomialOrder > BaseType ;
// type of Grid
typedef typename BaseType :: GridType GridType;
typedef typename BaseType :: GridType GridType;
// Choose a suitable GridView
typedef typename BaseType :: GridPartType GridPartType;
typedef typename BaseType :: GridPartType GridPartType;
// An analytical version of our model
typedef typename BaseType :: ModelType ModelType;
typedef typename BaseType :: ModelType ModelType;
// The DG space operator
// The first operator is sum of the other two
// The other two are needed for semi-implicit time discretization
typedef typename BaseType :: ExplicitOperatorType FullOperatorType;
typedef typename BaseType :: ExplicitOperatorType FullOperatorType;
typedef FullOperatorType ExplicitOperatorType;
typedef FullOperatorType ImplicitOperatorType;
typedef typename BaseType :: BasicLinearSolverType BasicLinearSolverType;
typedef typename BaseType :: BasicLinearSolverType BasicLinearSolverType;
// The discrete function for the unknown solution is defined in the DgOperator
typedef typename BaseType :: DiscreteFunctionType DiscreteFunctionType;
typedef typename BaseType :: DiscreteFunctionType DiscreteFunctionType;
// ... as well as the Space type
typedef typename BaseType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename BaseType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
// The ODE Solvers
typedef typename BaseType :: OdeSolverType OdeSolverType;
typedef typename BaseType :: OdeSolverType OdeSolverType;
typedef typename BaseType :: TimeProviderType TimeProviderType;
typedef typename BaseType :: TimeProviderType TimeProviderType;
typedef typename FullOperatorType :: ExtraParameterTupleType ExtraParameterTupleType;
typedef typename FullOperatorType :: ExtraParameterTupleType ExtraParameterTupleType;
typedef typename std::remove_pointer<typename std::tuple_element<0,typename BaseType::IndicatorTupleType>::type>::type DGIndicatorType;
typedef typename std::remove_pointer<typename std::tuple_element<1,typename BaseType::IndicatorTupleType>::type>::type GradientIndicatorType ;
// type of 64bit unsigned integer
typedef typename BaseType :: UInt64Type UInt64Type;
typedef typename BaseType :: UInt64Type UInt64Type;
using BaseType :: grid_;
using BaseType :: gridPart_;
......
......@@ -80,60 +80,6 @@ namespace Fem
};
class StepperParameters
: public Dune::Fem::LocalParameter< StepperParameters, StepperParameters >
{
protected:
const std::string keyPrefix_;
public:
StepperParameters( const std::string keyPrefix = "femdg.stepper." )
: keyPrefix_( keyPrefix )
{}
virtual double fixedTimeStep() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "fixedtimestep" , 0.0 );
}
virtual double fixedTimeStepEocLoopFactor() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "fixedtimestepeocloopfactor" , 1.0 );
}
virtual double startTime() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "starttime" , 0.0 );
}
virtual double endTime() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "endtime"/*, 1.0 */);
}
virtual int printCount() const
{
return Dune::Fem::Parameter::getValue< int >( keyPrefix_ + "printcount" , -1 );
}
virtual double maxTimeStep() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "maxtimestep", std::numeric_limits<double>::max());
}
virtual int maximalTimeSteps () const
{
return Dune::Fem::Parameter::getValue< int >( keyPrefix_ + "maximaltimesteps", std::numeric_limits<int>::max());
}
virtual bool stopAtEndTime() const
{
return Dune::Fem::Parameter::getValue< bool >( keyPrefix_ + "stopatendtime", bool(false) );
}
};
//! get memory in MB
inline double getMemoryUsage()
{
......@@ -153,18 +99,23 @@ namespace Fem
public:
// type of Grid
typedef typename Traits::GridType GridType;
typedef typename Traits::GridType GridType;
// type of IOTuple
typedef typename Traits::IOTupleType IOTupleType;
typedef typename Traits::IOTupleType IOTupleType;
typedef Fem::DataOutput<GridType, IOTupleType> DataOutputType;
// type of statistics monitor
typedef typename Traits::SolverMonitorType SolverMonitorType;
typedef typename Traits::SolverMonitorType SolverMonitorType;
typedef EocParameters EocParametersType;
//! constructor
AlgorithmBase ( GridType &grid, const std::string algorithmName = "" )
: grid_( grid ),
algorithmName_( algorithmName )
algorithmName_( algorithmName ),
eocParam_( Dune::ParameterKey::generate( "", "fem.eoc." ) )
{}
virtual ~AlgorithmBase () {}
......@@ -187,6 +138,12 @@ namespace Fem
virtual SolverMonitorType& monitor() = 0;
EocParametersType& eocParams()
{
return eocParam_;
}
//! finalize problem, i.e. calculated EOC ...
virtual void finalize ( const int eocloop ) = 0;
......@@ -198,6 +155,7 @@ namespace Fem
protected:
GridType &grid_;
const std::string algorithmName_;
EocParameters eocParam_;
};
......@@ -210,18 +168,17 @@ namespace Fem
template <class Algorithm>
void compute(Algorithm& algorithm)
{
const std::string name = "";
typedef typename Algorithm::GridType GridType;
typedef typename Algorithm::DataOutputType DataOutputType;
typedef typename Algorithm::IOTupleType IOTupleType;
typedef typename Algorithm::EocParametersType EocParametersType;
typedef typename Algorithm::GridType GridType;
GridType& grid = algorithm.grid();
typename Algorithm::IOTupleType dataTup = algorithm.dataTuple() ;
typedef Fem::DataOutput<GridType, typename Algorithm::IOTupleType> DataOutputType;
IOTupleType dataTup = algorithm.dataTuple() ;
DataOutputType dataOutput( grid, dataTup );
// initialize FemEoc if eocSteps > 1
EocParameters eocParam( Dune::ParameterKey::generate( "", "fem.eoc." ) );
EocParametersType& eocParam( algorithm.eocParams() );
FemEoc::clear();
FemEoc::initialize( eocParam.outputPath(), eocParam.fileName(), algorithm.description() );
......
......@@ -40,6 +40,62 @@ namespace Dune
namespace Fem
{
class StepperParameters
: public Dune::Fem::LocalParameter< StepperParameters, StepperParameters >
{
protected:
const std::string keyPrefix_;
public:
StepperParameters( const std::string keyPrefix = "femdg.stepper." )
: keyPrefix_( keyPrefix )
{}
virtual double fixedTimeStep() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "fixedtimestep" , 0.0 );
}
virtual double fixedTimeStepEocLoopFactor() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "fixedtimestepeocloopfactor" , 1.0 );
}
virtual double startTime() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "starttime" , 0.0 );
}
virtual double endTime() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "endtime"/*, 1.0 */);
}
virtual int printCount() const
{
return Dune::Fem::Parameter::getValue< int >( keyPrefix_ + "printcount" , -1 );
}
virtual double maxTimeStep() const
{
return Dune::Fem::Parameter::getValue< double >( keyPrefix_ + "maxtimestep", std::numeric_limits<double>::max());
}
virtual int maximalTimeSteps () const
{
return Dune::Fem::Parameter::getValue< int >( keyPrefix_ + "maximaltimesteps", std::numeric_limits<int>::max());
}
virtual bool stopAtEndTime() const
{
return Dune::Fem::Parameter::getValue< bool >( keyPrefix_ + "stopatendtime", bool(false) );
}
};
// EvolutionAlgorithmTraits
// -------------------------
......@@ -204,16 +260,14 @@ namespace Fem
typedef uint64_t UInt64Type ;
typedef StepperParameters StepperParametersType;
typedef EocParameters EocParametersType;
typedef AdaptationParameters AdaptationParametersType;
using BaseType::grid_;
using BaseType::eocParam_;
using BaseType::grid;
EvolutionAlgorithm ( GridType &grid, const std::string name = "" )
: BaseType( grid, name ),
param_( StepperParametersType( Dune::ParameterKey::generate( "", "femdg.stepper." ) ) ),
eocParam_( EocParametersType( Dune::ParameterKey::generate( "", "fem.eoc." ) ) ),
gridPart_( grid_ ),
space_( gridPart_ ),
solution_( "U_"+name, space() ),
......@@ -425,10 +479,10 @@ namespace Fem
void preInitializeStep ( TimeProviderType &tp, int loop )
{
// restoreData if checkpointing is enabled (default is disabled)
//bool newStart = ( eocParam_.steps() == 1) ? checkPointHandler_.restoreData( grid, tp ) : false;
bool newStart = false;
if( eocParam_.steps() == 1 )
checkPointHandler_.restoreData( tp );
bool newStart = ( eocParam_.steps() == 1) ? checkPointHandler_.restoreData( tp ) : false;
//bool newStart = false;
//if( eocParam_.steps() == 1 )
// checkPointHandler_.restoreData( tp );
initializeStep( tp, loop );
......@@ -550,7 +604,6 @@ namespace Fem
GridPartType gridPart_;
StepperParametersType param_;
EocParametersType eocParam_;
DiscreteFunctionSpaceType space_;
// the solution
......
#ifndef DUNE_FEM_DG_OPERATORTRAITS_HH
#define DUNE_FEM_DG_OPERATORTRAITS_HH
#include <dune/fem/quadrature/cachingquadrature.hh>
#include <dune/fem-dg/operator/fluxes/diffusionflux.hh>
#include <dune/fem-dg/operator/fluxes/upwindflux.hh>
#include <dune/fem-dg/operator/fluxes/eulerfluxes.hh>
namespace Dune {
// traits for the operator passes
template< class GridPart,
int polOrd,
class AnalyticalTraits,
class DiscreteFunctionImp,
class IndicatorFunctionImp,
class AdaptationHandlerImp,
class ExtraParameterTupleImp = std::tuple<>
>
struct OperatorTraits
{
typedef GridPart GridPartType;
typedef typename GridPartType::GridType GridType;
typedef AnalyticalTraits AnalyticalTraitsType;
typedef typename AnalyticalTraitsType::InitialDataType InitialDataType;
typedef typename AnalyticalTraitsType::ModelType ModelType ;
typedef UpwindFlux< ModelType > FluxType;
//typedef LLFFlux< ModelType > FluxType;
static const Dune::DGDiffusionFluxIdentifier PrimalDiffusionFluxId = Dune::method_general;
static const int polynomialOrder = polOrd == -1 ? 0 : polOrd;
typedef DiscreteFunctionImp DestinationType ;
typedef typename DestinationType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef Dune::Fem::CachingQuadrature< GridPartType, 0 > VolumeQuadratureType;
typedef Dune::Fem::CachingQuadrature< GridPartType, 1 > FaceQuadratureType;
typedef IndicatorFunctionImp IndicatorType;
typedef AdaptationHandlerImp AdaptationHandlerType ;
static const int limiterPolynomialOrder = polOrd == -1 ? 1 : polOrd;
typedef ExtraParameterTupleImp ExtraParameterTupleType;
};
} // end namespace Dune
#endif
#ifndef FEMDG_UPWIND_FLUX_HH
#define FEMDG_UPWIND_FLUX_HH
#include <dune/fem/version.hh>
#include <dune/fem/misc/fmatrixconverter.hh>
#include <dune/fem/io/parameter.hh>
#include <dune/fem-dg/operator/limiter/limitpass.hh>
// local includes
#include <dune/fem-dg/models/defaultmodel.hh>
#include <dune/fem-dg/pass/dgpass.hh>
/**
* @brief defines the advective flux
*/
template <class ModelType>
class UpwindFlux {
public:
typedef ModelType Model;
typedef typename Model::Traits Traits;
enum { dimRange = Model::dimRange };
typedef typename Model :: DomainType DomainType;
typedef typename Model :: RangeType RangeType;
typedef typename Model :: JacobianRangeType JacobianRangeType;
typedef typename Model :: FluxRangeType FluxRangeType;
typedef typename Model :: DiffusionRangeType DiffusionRangeType;
typedef typename Model :: FaceDomainType FaceDomainType;
typedef typename Model :: EntityType EntityType;
typedef typename Model :: IntersectionType IntersectionType;
public:
/**
* @brief Constructor
*/
UpwindFlux(const Model& mod) : model_(mod) {}
static std::string name () { return "UpwindFlux"; }
const Model& model() const {return model_;}
/**
* @brief evaluates the flux \f$g(u,v)\f$
*
* @return maximum wavespeed * normal
*/
template <class LocalEvaluation>
inline double numericalFlux( const LocalEvaluation& left,
const LocalEvaluation& right,
const RangeType& uLeft,
const RangeType& uRight,
const JacobianRangeType& jacLeft,
const JacobianRangeType& jacRight,
RangeType& gLeft,
RangeType& gRight ) const
{
const FaceDomainType& x = left.localPoint();
// get normal from intersection
const DomainType normal = left.intersection().integrationOuterNormal(x);
// get velocity
const DomainType v = model_.velocity( left );
const double upwind = normal * v;
if (upwind>0)
gLeft = uLeft;
else
gLeft = uRight;
gLeft *= upwind;
gRight = gLeft;
return std::abs(upwind);
}
protected:
const Model& model_;
};
#endif
......@@ -382,69 +382,4 @@ public:
const double tstepEps_;
};
/**
* @brief defines the advective flux
*/
template <class ModelType>
class UpwindFlux {
public:
typedef ModelType Model;
typedef typename Model::Traits Traits;
enum { dimRange = Model::dimRange };
typedef typename Model :: DomainType DomainType;
typedef typename Model :: RangeType RangeType;
typedef typename Model :: JacobianRangeType JacobianRangeType;
typedef typename Model :: FluxRangeType FluxRangeType;
typedef typename Model :: DiffusionRangeType DiffusionRangeType;
typedef typename Model :: FaceDomainType FaceDomainType;
typedef typename Model :: EntityType EntityType;
typedef typename Model :: IntersectionType IntersectionType;
public:
/**
* @brief Constructor
*/
UpwindFlux(const Model& mod) : model_(mod) {}
static std::string name () { return "UpwindFlux"; }
const Model& model() const {return model_;}
/**
* @brief evaluates the flux \f$g(u,v)\f$
*
* @return maximum wavespeed * normal
*/
template <class LocalEvaluation>
inline double numericalFlux( const LocalEvaluation& left,
const LocalEvaluation& right,
const RangeType& uLeft,
const RangeType& uRight,
const JacobianRangeType& jacLeft,
const JacobianRangeType& jacRight,
RangeType& gLeft,
RangeType& gRight ) const
{
const FaceDomainType& x = left.localPoint();
// get normal from intersection
const DomainType normal = left.intersection().integrationOuterNormal(x);
// get velocity
const DomainType v = model_.velocity( left );
const double upwind = normal * v;
if (upwind>0)
gLeft = uLeft;
else
gLeft = uRight;
gLeft *= upwind;
gRight = gLeft;
return std::abs(upwind);
}
protected:
const Model& model_;
};
#endif
......@@ -27,6 +27,7 @@
#include <dune/fem/operator/linear/spoperator.hh>
#include <dune/fem-dg/solver/linearsolvers.hh>
#include <dune/fem-dg/operator/dg/operatortraits.hh>
#include <dune/fem/misc/l2norm.hh>
......@@ -38,44 +39,6 @@
#include "models.hh"
// traits for the operator passes
template< class GridPart,
int polOrd,
class AnalyticalTraits,
class DiscreteFunctionImp,
class IndicatorFunctionImp,
class AdaptationHandlerImp,
class ExtraParameterTupleImp = std::tuple<>
>
struct OperatorTraits
{
typedef GridPart GridPartType;
typedef typename GridPartType::GridType GridType;
typedef AnalyticalTraits AnalyticalTraitsType;
typedef typename AnalyticalTraitsType::InitialDataType InitialDataType;
typedef typename AnalyticalTraitsType::ModelType ModelType ;
typedef UpwindFlux< ModelType > FluxType;
//typedef LLFFlux< ModelType > FluxType;
static const Dune::DGDiffusionFluxIdentifier PrimalDiffusionFluxId = Dune::method_general;
static const int polynomialOrder = polOrd == -1 ? 0 : polOrd;
typedef DiscreteFunctionImp DestinationType ;
typedef typename DestinationType :: DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef Dune::Fem::CachingQuadrature< GridPartType, 0 > VolumeQuadratureType;
typedef Dune::Fem::CachingQuadrature< GridPartType, 1 > FaceQuadratureType;
typedef IndicatorFunctionImp IndicatorType;
typedef AdaptationHandlerImp AdaptationHandlerType ;
static const int limiterPolynomialOrder = polOrd == -1 ? 1 : polOrd;
typedef ExtraParameterTupleImp ExtraParameterTupleType;
};
enum AdvectionDiffusionOperatorType
{
_unlimited = 0,
......@@ -337,9 +300,9 @@ public:
// --------- Operators using PASSES --------------------------
//============================================================
typedef OperatorTraits< GridPartType, polynomialOrder, AnalyticalTraitsType,
DiscreteFunctionType, IndicatorType,
AdaptationHandlerType, ExtraParameterTuple > OperatorTraitsType;
typedef Dune::OperatorTraits< GridPartType, polynomialOrder, AnalyticalTraitsType,
DiscreteFunctionType, IndicatorType,
AdaptationHandlerType, ExtraParameterTuple > OperatorTraitsType;
// TODO: advection/diffusion should not be precribed by model
static const int hasAdvection = AnalyticalTraitsType::ModelType::hasAdvection;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment