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

poisson uses new algorithm classes, now

parent 1c8f6f5b
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 959 additions and 691 deletions
dune/fem-dg/algorithm/base.hh
+ 1
− 1
View file @ 422f730a
......@@ -133,7 +133,7 @@ namespace Fem
//! return default data tuple for data output
virtual IOTupleType dataTuple() = 0;
// solve the problem for eoc loop 'loop' and return statistics
// solve the problem for eoc loop 'loop'
virtual void solve ( const int loop ) = 0;
virtual SolverMonitorType& monitor() = 0;
......
dune/fem-dg/algorithm/ellipticalgorithm.hh
+ 412
− 495
View file @ 422f730a
This diff is collapsed.
dune/fem-dg/algorithm/evolution.hh
+ 2
− 4
View file @ 422f730a
......@@ -127,7 +127,6 @@ namespace Fem
typedef typename AnalyticalTraits::InitialDataType InitialDataType;
// type of discrete function space and discrete function
typedef typename DiscreteTraits::InitialProjectorType InitialProjectorType;
// type of dg operator
......@@ -146,8 +145,6 @@ namespace Fem
typedef typename DiscreteTraits::IOTupleType IOTupleType;
// wrap operator
//typedef typename DgHelmHoltzOperatorType::JacobianOperatorType JacobianOperatorType;
//typedef typename ProblemTraits::template Solver< DgHelmHoltzOperatorType >::Type OdeSolverType;
typedef GridTimeProvider< GridType > TimeProviderType;
typedef typename DiscreteTraits::OdeSolverType OdeSolverType;
......@@ -286,8 +283,9 @@ namespace Fem
timeStepTimer_( Dune::FemTimer::addTo("max time/timestep") ),
fixedTimeStep_( param_.fixedTimeStep() )
{}
// return grid width of grid (overload in derived classes)
virtual double gridWidth () const { return 0.0; }
virtual double gridWidth () const { return GridWidth::calcGridWidth( gridPart_ ); }
// return size of grid
virtual UInt64Type gridSize () const
......
dune/fem-dg/algorithm/handler/solvermonitor.hh
+ 46
− 0
View file @ 422f730a
......@@ -54,6 +54,52 @@ namespace Fem
};
class DefaultSteadyStateSolverMonitorHandler
{
public:
typedef SolverMonitor<1> SolverMonitorType;
DefaultSteadyStateSolverMonitorHandler( const std::string keyPrefix = "" )
: solverMonitor_()
{}
void stepPrint()
{
std::cout << ", Newton: " << *solverMonitor_.newton_iterations
<< " ILS: " << *solverMonitor_.ils_iterations << std::endl;
}
template< class OdeSolverMonitorImp, class TimeProviderImp >
void step( OdeSolverMonitorImp& odeSolverMonitor, TimeProviderImp& tp )
{
}
template< class LinearSolverMonitorImp >
void finalize( const double gridWidth, const double gridSize, const LinearSolverMonitorImp& solver )
{
*solverMonitor_.gridWidth = gridWidth;
*solverMonitor_.elements = gridSize;
//*solverMonitor_.ils_iterations = invDgOperator_->iterations();
//*solverMonitor_.totalNewtonIterations_ = solver.iterations();
//*solverMonitor_.totalLinearSolverIterations_ = solver.linearIterations();
*solverMonitor_.newton_iterations = 0;
*solverMonitor_.ils_iterations = solver.iterations();
*solverMonitor_.max_newton_iterations = 0;
*solverMonitor_.max_ils_iterations = 0;
*solverMonitor_.operator_calls = 0;
}
SolverMonitorType& monitor()
{
return solverMonitor_;
}
private:
SolverMonitorType solverMonitor_;
};
class NoSolverMonitorHandler
{
......
dune/fem-dg/algorithm/steadystate.hh
+ 154
− 122
View file @ 422f730a
......@@ -13,198 +13,230 @@
#include <dune/fem-dg/algorithm/base.hh>
#include <dune/fem-dg/algorithm/handler/solvermonitor.hh>
namespace Dune
{
namespace Fem
{
template< class Grid,
template< class GridImp,
class ProblemTraits,
int polOrd,
ExtraParameterTuple = std::tuple<> >
class SolverMonitorHandlerImp >
struct SteadyStateTraits
{
enum { polynomialOrder = polOrd };
// type of Grid
typedef Grid GridType;
typedef GridImp GridType;
// Choose a suitable GridView
typedef AdaptiveLeafGridPart< GridType > HostGridPartType;
typedef AdaptiveLeafGridPart< GridType > HostGridPartType;
typedef typename ProblemTraits::GridPartType GridPartType;
typedef typename ProblemTraits::GridPartType GridPartType;
typedef typename ProblemTraits::AnalyticalTraits ModelTraits;
typedef typename ProblemTraits::template AnalyticalTraits< HostGridPartType > AnalyticalTraits;
typedef typename ProblemTraits::template DiscreteTraits< HostGridPartType, polynomialOrder > DiscreteTraits;
// obtain the problem dependent types, analytical context
typedef typename ModelTraits::ModelType ModelType;
typedef typename ModelTraits::ProblemType ProblemType;
typedef typename AnalyticalTraits::ModelType ModelType;
typedef typename AnalyticalTraits::ProblemType ProblemType;
struct OperatorTraits :
public Dune::PassTraits< ModelTraits, polynomialOrder, ModelType::dimRange >
{
static const int limiterPolynomialOrder = polynomialOrder;
typedef ExtraParameterTuple ExtraParameterTupleType;
};
// type of discrete function space and discrete function
typedef typename ProblemTraits::FunctionSpaceType FunctionSpaceType;
typedef typename DiscreteTraits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteTraits::DiscreteFunctionType DiscreteFunctionType;
//typedef typename DiscreteTraits::OperatorTraitsType OperatorTraits;
// type of discrete function space and discrete function
typedef typename OperatorTraits::DestinationType DiscreteFunctionType;
typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteTraits::ExtraParameterTuple ExtraParameterTuple;
typedef typename DiscreteTraits::FluxType FluxType;
typedef typename OperatorTraits::LinearOperatorType LinearOperatorType;
typedef typename OperatorTraits::LinearInverseOperatorType LinearInverseOperatorType;
//typedef typename DiscreteTraits::OdeSolverType OdeSolverType;
typedef typename DiscreteTraits::BasicLinearSolverType BasicLinearSolverType;
typedef Dune::DGAdvectionDiffusionOperator< OperatorTraits > FullOperatorType;
typedef DGPrimalMatrixAssembly< FullOperatorType > AssembledOperatorType;
typedef SolverMonitor<1> SolverMonitorType;
// tpye of jacobian operator used in the nested newton loops
typedef typename FullOperatorType::JacobianOperatorType JacobianOperatorType;
typedef typename DiscreteTraits::FullOperatorType FullOperatorType;
typedef typename DiscreteTraits::AssemblyOperatorType AssemblyOperatorType;
typedef typename DiscreteTraits::AssembledOperatorType AssembledOperatorType;
// type of IOTuple
typedef Dune::tuple< DiscreteFunctionType * > IOTupleType;
typedef typename DiscreteTraits::IOTupleType IOTupleType;
// type of restriction/prolongation projection for adaptive simulations
typedef Dune::Fem::RestrictProlongDefault< DiscreteFunctionType > RestrictionProlongationType;
typedef typename DiscreteTraits::RestrictionProlongationType RestrictionProlongationType;
typedef typename AnalyticalTraits::EOCErrorIDs EOCErrorIDs;
typedef SolverMonitorHandlerImp SolverMonitorHandlerType;
};
};
template< class Grid, class ProblemTraits, int polynomialOrder,
class SolverMonitorHandlerImp = DefaultSteadyStateSolverMonitorHandler >
class SteadyStateAlgorithm
: public AlgorithmBase< SteadyStateTraits< Grid, ProblemTraits, polynomialOrder,
SolverMonitorHandlerImp > >
{
typedef SteadyStateTraits< Grid, ProblemTraits, polynomialOrder, SolverMonitorHandlerImp > Traits;
typedef AlgorithmBase< Traits > BaseType;
public:
typedef typename BaseType::GridType GridType;
typedef typename BaseType::IOTupleType IOTupleType;
typedef typename BaseType::SolverMonitorType SolverMonitorType;
typedef typename Traits::HostGridPartType HostGridPartType;
template< class Grid, class ProblemTraits, int polynomialOrder, class ExtraParameterTuple = std::tuple< > >
class SteadyStateAlgorithm
: public AlgorithmBase< SteadyStateTraits< Grid, ProblemTraits, polynomialOrder, ExtraParameterTuple > >
{
// my traits class
typedef SteadyStateTraits< Grid, ProblemTraits, polynomialOrder, ExtraParameterTuple > Traits;
// initial data type
typedef typename Traits::ProblemType ProblemType;
typedef AlgorithmBase< Traits > BaseType;
typedef SteadyStateAlgorithm< Grid, ProblemTraits, polynomialOrder > This;
// An analytical version of our model
typedef typename Traits::ModelType ModelType;
public:
typedef typename BaseType::GridType GridType;
typedef typename BaseType::IOTupleType IOTupleType;
typedef typename BaseType::SolverMonitorType SolverMonitorType;
typedef typename Traits::GridPartType GridPartType;
typedef typename Traits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename Traits::DiscreteFunctionType DiscreteFunctionType;
typedef typename Traits::HostGridPartType HostGridPartType;
// The DG space operator
typedef typename Traits::FullOperatorType FullOperatorType;
// initial data type
typedef typename Traits::ProblemType ProblemType;
// type of steady state solver
typedef typename Traits::BasicLinearSolverType BasicLinearSolverType;
// An analytical version of our model
typedef typename Traits::ModelType ModelType;
// type of analytical traits
typedef typename Traits::AnalyticalTraits AnalyticalTraits;
typedef typename Traits::GridPartType GridPartType;
typedef typename Traits::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename Traits::DiscreteFunctionType DiscreteFunctionType;
// type of discrete traits
typedef typename Traits::DiscreteTraits DiscreteTraits;
// The DG space operator
typedef typename Traits::FullOperatorType FullOperatorType;
typedef typename Traits::RestrictionProlongationType RestrictionProlongationType;
// type of steady state solver
typedef typename Traits::InverseOperatorType InverseOperatorType;
typedef typename Traits::AssemblyOperatorType AssemblyOperatorType;
typedef typename Traits::AssembledOperatorType AssembledOperatorType;
// type of analytical traits
typedef typename Traits::AnalyticalTraits AnalyticalTraits;
// Error Handling
typedef typename Traits::EOCErrorIDs EOCErrorIDs;
// type of discrete traits
typedef typename Traits::DiscreteTraits DiscreteTraits;
typedef SolverMonitorHandlerImp SolverMonitorHandlerType;
// Error Handling
typedef typename AnalyticalTraits::EOCErrorIDs EOCErrorIDs;
typedef uint64_t UInt64Type ;
SteadyStateAlgorithm ( GridType &grid )
: BaseType( grid ),
problem_( AnalyticalTraits::problem() ),
model_( problem() ),
gridPart_( grid ),
space_( gridPart_ ),
solution_( "solution", space_ ),
dgOperator_( space_, model(), DiscreteTraits::quadOrder ),
eocIds_( AnalyticalTraits::initEoc() ),
description_( "Steady state scheme,\nOperator: " + dgOperator_.description() +"\nProblem: " + problem().description() )
{}
using BaseType::grid;
//! return reference to discrete space
DiscreteFunctionSpaceType &space () { return space_; }
SteadyStateAlgorithm ( GridType &grid, const std::string name = "" )
: BaseType( grid, name ),
problem_( ProblemTraits::problem() ),
model_( problem() ),
solverMonitorHandler_( "" ),
gridPart_( grid ),
space_( gridPart_ ),
solution_( "solution-" + name, space_ ),
rhs_( "rhs-" + name, space_ ),
eocIds_( AnalyticalTraits::initEoc() )
{
solution().clear();
}
//! returns data prefix for EOC loops ( default is loop )
virtual std::string dataPrefix () const
{
return problem().dataPrefix();
}
//! return reference to discrete space
DiscreteFunctionSpaceType &space () { return space_; }
IOTupleType dataTuple ()
{
return IOTupleType( &solution_ );
}
DiscreteFunctionType& solution ()
{
return solution_;
}
SolverMonitorType solve ( int step )
{
SolverMonitorType monitor;
DiscreteFunctionType rhs( "rhs", space_ );
std::string description () const { return problem().description(); }
rhs.clear();
solution_.clear();
//! returns data prefix for EOC loops ( default is loop )
virtual std::string dataPrefix () const
{
return problem().dataPrefix();
}
InverseOperatorType invOp( dgOperator_, NestedNewtonParameter() );
IOTupleType dataTuple ()
{
IOTupleType test = std::make_tuple( &solution_ );
return test;
}
virtual SolverMonitorType& monitor()
{
return solverMonitorHandler_.monitor();
}
// return grid width of grid (overload in derived classes)
virtual double gridWidth () const { return GridWidth::calcGridWidth( gridPart_ ); }
// return size of grid
virtual UInt64Type gridSize () const
{
UInt64Type grSize = grid().size( 0 );
return grid().comm().sum( grSize );
}
auto callBack = [ this ] ( const double f ) { dgOperator_.setRelaxFactor( f ); };
invOp( rhs, solution_, callBack );
virtual BasicLinearSolverType* createSolver() = 0;
monitor.gridwidth_ = GridWidth::calcGridWidth( gridPart_ );
monitor.totalNewtonIterations_ = invOp.iterations();
monitor.totalLinearSolverIterations_ = invOp.linearIterations();
virtual void solve ( int step )
{
rhs_.clear();
solution_.clear();
solver_.reset( this->createSolver() );
assert( solver_ );
monitor.finalize();
(*solver_)( rhs_, solution_ );
return monitor;
}
solverMonitorHandler_.finalize( gridWidth(), gridSize(), *solver_ );
}
//! finalize computation by calculating errors and EOCs
void finalize ( const int eocloop )
{
// submit error to the FEM EOC calculator
AnalyticalTraits::addEOCErrors( eocIds_, 0., solution_, model(), problem() );
}
//! finalize computation by calculating errors and EOCs
void finalize ( const int eocloop )
{
// submit error to the FEM EOC calculator
//AnalyticalTraits::addEOCErrors( eocIds_, solution_, model(), problem() );
}
std::string description () const { return description_; }
const ProblemType &problem () const
{
assert( problem_ );
return *problem_;
}
const ProblemType &problem () const
{
assert( problem_ );
return *problem_;
}
ProblemType &problem ()
{
assert( problem_ );
return *problem_;
}
ProblemType &problem ()
{
assert( problem_ );
return *problem_;
}
ModelType &model () { return model_; }
const ModelType &model () const { return model_; }
ModelType &model () { return model_; }
const ModelType &model () const { return model_; }
protected:
ProblemType *problem_;
ModelType model_;
protected:
ProblemType *problem_;
ModelType model_;
SolverMonitorHandlerType solverMonitorHandler_;
GridPartType gridPart_; // reference to grid part, i.e. the leaf grid
DiscreteFunctionSpaceType space_; // the discrete function space
DiscreteFunctionType solution_;
GridPartType gridPart_; // reference to grid part, i.e. the leaf grid
DiscreteFunctionSpaceType space_; // the discrete function space
DiscreteFunctionType solution_;
DiscreteFunctionType rhs_;
FullOperatorType dgOperator_;
EOCErrorIDs eocIds_;
std::unique_ptr< BasicLinearSolverType > solver_;
std::string description_;
};
EOCErrorIDs eocIds_;
};
} // namespace Fem
} // namespace Fem
} // namespace Dune
......
dune/fem-dg/assemble/primalmatrix.hh
+ 3
− 3
View file @ 422f730a
......@@ -29,9 +29,9 @@ class DGPrimalMatrixAssembly
typedef OperatorImp OperatorType;
typedef typename OperatorType::Traits::ModelType ModelType;
typedef typename OperatorType::Traits::DestinationType DestinationType;
typedef typename OperatorType::Traits::DiscreteModelType DiscreteModelType;
typedef typename DiscreteModelType::Selector SelectorType;
static const Dune::DGDiffusionFluxIdentifier DGDiffusionFluxIdentifier = OperatorType::Traits::PrimalDiffusionFluxId;
//typedef typename OperatorType::Traits::DiscreteModelType DiscreteModelType;
//typedef typename DiscreteModelType::Selector SelectorType;
//static const Dune::DGDiffusionFluxIdentifier DGDiffusionFluxIdentifier = OperatorType::Traits::PrimalDiffusionFluxId;
typedef typename DestinationType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionSpaceType::IteratorType IteratorType;
typedef typename IteratorType::Entity EntityType;
......
dune/fem-dg/test/poisson/main.cc
+ 1
− 1
View file @ 422f730a
......@@ -24,7 +24,7 @@ int main(int argc, char ** argv)
PoissonProblemCreator< GridType > problem;
// run simulation
Simulator::run( problem );
Dune::Simulator::run( problem );
}
catch (const Dune::Exception &e)
{
......
dune/fem-dg/test/poisson/parameter.in
+ 1
− 3
View file @ 422f730a
......@@ -15,8 +15,6 @@ fem.eoc.steps: 30
# polynomial order = 0,1,2,3,4
femdg.polynomialOrder: 2
# PROBLEM
# ----------
# problem number
......@@ -65,6 +63,6 @@ fem.io.grapedisplay: 0
# SOLVER
# ------
paramfile: ../parameter/paramSolver
istl.preconditioning.method: ilu-0
istl.preconditioning.method: amg-ilu-0
istl.preconditioning.iterations: 1
istl.preconditioning.relaxation: 1
dune/fem-dg/test/poisson/poissonproblem.hh
+ 0
− 2
View file @ 422f730a
......@@ -253,8 +253,6 @@ public:
std::ostringstream ofs;
ofs << "Problem: " << myName ;
ofs << ", End time: " << Dune:: Fem :: Parameter::getValue<double>("femdg.stepper.endtime");
return ofs.str();
}
......
dune/fem-dg/test/poisson/problemcreator.hh
+ 339
− 60
View file @ 422f730a
#ifndef FEMHOWTO_POISSONSTEPPER_HH
#define FEMHOWTO_POISSONSTEPPER_HH
#ifndef FEMDG_POISSONSTEPPER_HH
#define FEMDG_POISSONSTEPPER_HH
#include <config.h>
//#define WANT_ISTL 1
#ifndef NDEBUG
// enable fvector and fmatrix checking
#define DUNE_ISTL_WITH_CHECKING
#ifndef DIMRANGE
#define DIMRANGE 1
#endif
#ifndef POLORDER
#define POLORDER 1
#endif
#include "passtraits.hh"
// dune-grid includes
#include <dune/fem/io/parameter.hh>
#include <dune/grid/io/file/dgfparser/dgfparser.hh>
#include <dune/fem/misc/l2norm.hh>
// dune-fem includes
#include <dune/fem/io/parameter.hh>
//--------- GRID HELPER ---------------------
#include <dune/fem-dg/algorithm/gridinitializer.hh>
//--------- OPERATOR/SOLVER -----------------
#include <dune/fem-dg/assemble/primalmatrix.hh>
#include <dune/fem-dg/solver/linearsolvers.hh>
#include <dune/fem-dg/operator/dg/operatortraits.hh>
//--------- STEPPER -------------------------
#include <dune/fem-dg/algorithm/ellipticalgorithm.hh>
//--------- PROBLEMS ------------------------
#include "poissonproblem.hh"
//--------- MODELS --------------------------
#include "models.hh"
// dune-fem-dg includes
#include <dune/fem-dg/operator/fluxes/diffusionflux.hh>
enum AdvectionDiffusionOperatorType
{
_unlimited = 0,
_limited = 1,
};
#include <dune/fem-dg/operator/dg/dgoperatorchoice.hh>
#include <dune/fem-dg/operator/fluxes/noflux.hh>
#include <dune/fem-dg/assemble/primalmatrix.hh>
#include <dune/fem-dg/solver/linearsolvers.hh>
#include <dune/fem-dg/stepper/ellipticalgorithm.hh>
template< class Op, class DiffusionOp, class AdvectionOp, bool advection, bool diffusion >
struct OperatorChooser
{
typedef Op FullOperatorType;
typedef DiffusionOp ImplicitOperatorType;
typedef AdvectionOp ExplicitOperatorType;
};
template< class Op, class DiffusionOp, class AdvectionOp, bool advection >
struct OperatorChooser< Op, DiffusionOp, AdvectionOp, advection, false >
{
typedef AdvectionOp FullOperatorType;
typedef FullOperatorType ImplicitOperatorType;
typedef FullOperatorType ExplicitOperatorType;
};
template<class Op, class DiffusionOp, class AdvectionOp, bool diffusion >
struct OperatorChooser< Op, DiffusionOp, AdvectionOp, false, diffusion >
{
typedef DiffusionOp FullOperatorType;
typedef FullOperatorType ImplicitOperatorType;
typedef FullOperatorType ExplicitOperatorType;
};
template< class OperatorTraits, bool advection, bool diffusion, AdvectionDiffusionOperatorType op >
class AdvectionDiffusionOperators;
template< class OperatorTraits, bool advection, bool diffusion >
class AdvectionDiffusionOperators< OperatorTraits, advection, diffusion, _unlimited >
{
typedef Dune::DGAdvectionDiffusionOperator< OperatorTraits > DgType;
typedef Dune::DGAdvectionOperator< OperatorTraits > DgAdvectionType;
typedef Dune::DGDiffusionOperator< OperatorTraits > DgDiffusionType;
typedef OperatorChooser< DgType, DgAdvectionType, DgDiffusionType, advection, diffusion >
OperatorChooserType;
public:
typedef typename OperatorChooserType::FullOperatorType FullOperatorType;
typedef typename OperatorChooserType::ImplicitOperatorType ImplicitOperatorType;
typedef typename OperatorChooserType::ExplicitOperatorType ExplicitOperatorType;
};
template< class OperatorTraits, bool advection, bool diffusion >
class AdvectionDiffusionOperators< OperatorTraits, advection, diffusion, _limited >
{
typedef Dune::DGLimitedAdvectionDiffusionOperator< OperatorTraits > DgType;
typedef Dune::DGLimitedAdvectionOperator< OperatorTraits > DgAdvectionType;
typedef Dune::DGDiffusionOperator< OperatorTraits > DgDiffusionType;
typedef OperatorChooser< DgType, DgAdvectionType, DgDiffusionType, advection, diffusion >
OperatorChooserType;
public:
typedef typename OperatorChooserType::FullOperatorType FullOperatorType;
typedef typename OperatorChooserType::ImplicitOperatorType ImplicitOperatorType;
typedef typename OperatorChooserType::ExplicitOperatorType ExplicitOperatorType;
};
enum DiscreteFunctionSpacesType
{
_lagrange = 0,
_legendre = 1,
};
template< class FunctionSpaceImp, class GridPartImp, int polOrder, DiscreteFunctionSpacesType dfType, GalerkinType opType >
struct DiscreteFunctionSpaces;
template< class FunctionSpaceImp, class GridPartImp, int polOrder >
struct DiscreteFunctionSpaces< FunctionSpaceImp, GridPartImp, polOrder, _lagrange, cg >
{
typedef Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpaceImp, GridPartImp, polOrder, Dune::Fem::CachingStorage > type;
};
template< class FunctionSpaceImp, class GridPartImp, int polOrder >
struct DiscreteFunctionSpaces< FunctionSpaceImp, GridPartImp, polOrder, _legendre, dg >
{
typedef Dune::Fem::DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, Dune::Fem::CachingStorage > type;
};
//template< class TimeDependentFunctionImp, class DiscreteFunctionImp, GalerkinType gt >
//struct InitialProjectors;
//
//template< class TimeDependentFunctionImp, class DiscreteFunctionImp >
//struct InitialProjectors< TimeDependentFunctionImp, DiscreteFunctionImp, cg >
//{
// typedef Dune::Fem::LagrangeInterpolation< TimeDependentFunctionImp, DiscreteFunctionImp > type;
//};
//
//template< class TimeDependentFunctionImp, class DiscreteFunctionImp >
//struct InitialProjectors< TimeDependentFunctionImp, DiscreteFunctionImp, dg >
//{
// typedef Dune::Fem::L2Projection< TimeDependentFunctionImp, DiscreteFunctionImp > type;
//};
template< class DiscreteFunctionSpaceImp, SolverType solverType >
struct DiscreteFunctions;
template< class DiscreteFunctionSpaceImp >
struct DiscreteFunctions< DiscreteFunctionSpaceImp, fem >
{
typedef Dune::Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceImp > type;
typedef Dune::Fem::SparseRowLinearOperator< type, type > jacobian;
};
#if HAVE_DUNE_ISTL
template< class DiscreteFunctionSpaceImp >
struct DiscreteFunctions< DiscreteFunctionSpaceImp, istl >
{
typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpaceImp > type;
typedef Dune::Fem::ISTLLinearOperator< type, type > jacobian;
};
#endif
#if HAVE_PETSC
template< class DiscreteFunctionSpaceImp >
struct DiscreteFunctions< DiscreteFunctionSpaceImp, petsc >
{
typedef Dune::Fem::PetscDiscreteFunction< DiscreteFunctionSpaceImp > type;
typedef Dune::Fem::PetscLinearOperator< type, type > jacobian;
};
#endif
// local includes
#include "poissonproblem.hh"
#include "models.hh"
template <class GridType>
template< class GridImp >
struct PoissonProblemCreator
{
static const int dimRange = 1 ;
typedef Dune :: ProblemInterface<
Dune::Fem::FunctionSpace< double, double, GridType::dimension, dimRange> > ProblemType;
typedef GridImp GridType;
typedef Dune::Fem::DGAdaptiveLeafGridPart< GridType > HostGridPartType;
typedef HostGridPartType GridPartType;
template <class GridPart>
struct Traits
{
typedef ProblemType InitialDataType;
typedef PoissonModel< GridPart, InitialDataType > ModelType;
//typedef Dune::LLFAdvFlux< ModelType > FluxType;
typedef Dune::UpwindFlux< ModelType > FluxType;
//typedef Dune::NoFlux< ModelType > FluxType;
// choice of diffusion flux (see diffusionflux.hh for methods)
static const Dune :: DGDiffusionFluxIdentifier PrimalDiffusionFluxId
= Dune :: method_general ;
};
// define problem type here if interface should be avoided
typedef Dune::Fem::FunctionSpace< typename GridType::ctype, double, GridType::dimension, 1 >
FunctionSpaceType;
typedef Dune::ProblemInterface< FunctionSpaceType > ProblemInterfaceType;
static inline std::string diffusionFluxName()
template< class GridPart > // TODO: is this template parameter needed?
struct AnalyticalTraits
{
#if (defined PRIMALDG)
return Dune::Fem::Parameter::getValue< std::string >("dgdiffusionflux.method");
#else
return "LDG";
#endif
}
typedef ProblemInterfaceType ProblemType;
typedef ProblemInterfaceType InitialDataType;
typedef PoissonModel< GridPart, InitialDataType > ModelType;
//typedef typename InitialDataType::TimeDependentFunctionType TimeDependentFunctionType;
typedef std::vector< int > EOCErrorIDs;
static EOCErrorIDs initEoc ()
{
EOCErrorIDs ids;
ids.resize( 3 );
ids[ 0 ] = Dune::Fem::FemEoc::addEntry( std::string( "$L^2$-error" ) );
ids[ 1 ] = Dune::Fem::FemEoc::addEntry( std::string( "DG-error" ) );
ids[ 2 ] = Dune::Fem::FemEoc::addEntry( std::string( "sigma-norm" ) );
return ids;
}
template< class SolutionImp, class Model, class ExactFunction, class SigmaEstimator >
static void addEOCErrors ( const EOCErrorIDs &ids, SolutionImp &u, const Model &model, ExactFunction &f, const SigmaEstimator& sigma )
{
// calculate L2 - Norm
Dune::Fem::L2Norm< GridPartType > l2norm( u.space().gridPart() );
const double l2error = l2norm.distance( f, u );
Dune::Fem::DGNorm< GridPartType > dgnorm( u.space().gridPart() );
const double dgerror = dgnorm.distance( f, u );
Dune::Fem::H1Norm< GridPartType > sigmanorm( u.space().gridPart() );
const double sigmaerror = 0.0;// sigmanorm.distance( f, sigma );
Dune::Fem::FemEoc::setErrors( ids[ 0 ], l2error );
Dune::Fem::FemEoc::setErrors( ids[ 1 ], dgerror );
Dune::Fem::FemEoc::setErrors( ids[ 2 ], sigmaerror );
}
};
static inline std::string moduleName()
{
return "";
}
static inline Dune::GridPtr<GridType> initializeGrid()
static inline Dune::GridPtr<GridType>
initializeGrid()
{
// use default implementation
std::string filename = ParameterType::commonInputPath() + "/" +
// use default implementation
std::string filename = Dune::Fem::Parameter::commonInputPath() + "/" +
Dune::Fem::Parameter::getValue< std::string >(Dune::Fem::IOInterface::defaultGridKey(GridType :: dimension, false));
std::string description ("poisson-"+diffusionFluxName());
// initialize grid
Dune::GridPtr< GridType > gridptr = initialize< GridType >( description );
// use default implementation
Dune::GridPtr<GridType> gridptr = Dune::Fem::DefaultGridInitializer< GridType >::initialize();
std::string::size_type position = std::string::npos;
if( GridType::dimension == 2)
......@@ -193,24 +342,154 @@ struct PoissonProblemCreator
grid.adapt();
grid.postAdapt();
//Dune :: GrapeGridDisplay< GridType > grape( grid );
//grape.display ();
return gridptr ;
}
static ProblemType* problem( )
static ProblemInterfaceType* problem()
{
// choice of benchmark problem
int probNr = Dune::Fem::Parameter::getValue< int > ( "problem" );
return new Dune :: PoissonProblem< GridType,dimRange > ( probNr );
return new Dune :: PoissonProblem< GridType, 1 > ( probNr );
}
template <int polynomialOrder>
//Stepper Traits
template< class GridPart, int polOrd > // TODO: is GridPart as a template parameter needed?
struct DiscreteTraits
{
public:
typedef AnalyticalTraits< GridPartType > AnalyticalTraitsType;
static const int polynomialOrder = polOrd;
static inline std::string advectionFluxName()
{
return "LLF";
}
static inline std::string diffusionFluxName()
{
return Dune::Fem::Parameter::getValue< std::string >("dgdiffusionflux.method");
}
static const int quadOrder = polynomialOrder * 3 + 1;
static const SolverType solverType = istl ;
typedef typename DiscreteFunctionSpaces< FunctionSpaceType, GridPartType, polynomialOrder, _legendre, dg >::type DiscreteFunctionSpaceType;
typedef typename DiscreteFunctions< DiscreteFunctionSpaceType, solverType >::type DiscreteFunctionType;
typedef typename DiscreteFunctions< DiscreteFunctionSpaceType, solverType >::jacobian JacobianOperatorType;
//typedef typename InitialProjectors< typename AnalyticalTraitsType::TimeDependentFunctionType, DiscreteFunctionType, dg >::type InitialProjectorType;
typedef std::tuple<> ExtraParameterTuple;
typedef std::tuple< DiscreteFunctionType* > IOTupleType;
private:
typedef Dune::Fem::FunctionSpace< typename GridType::ctype, double, AnalyticalTraitsType::ModelType::dimDomain, 3> FVFunctionSpaceType;
typedef Dune::Fem::FiniteVolumeSpace<FVFunctionSpaceType,GridPartType, 0, Dune::Fem::SimpleStorage> IndicatorSpaceType;
public:
typedef Dune::Fem::AdaptiveDiscreteFunction<IndicatorSpaceType> IndicatorType;
//typedef DuneODE::OdeSolverInterface< DiscreteFunctionType > OdeSolverType;
// type of restriction/prolongation projection for adaptive simulations
typedef Dune::Fem::RestrictProlongDefault< DiscreteFunctionType > RestrictionProlongationType;
// type of linear solver for implicit ode
// typedef Dune::Fem::ParDGGeneralizedMinResInverseOperator< DiscreteFunctionType > BasicLinearSolverType;
typedef Dune::AdaptationHandler< GridType, FunctionSpaceType > AdaptationHandlerType;
//############################ poisson issues ############################
static const bool symmetricSolver = true ;
typedef Solvers<DiscreteFunctionSpaceType, solverType, symmetricSolver> SolversType;
typedef Dune::UpwindFlux< typename AnalyticalTraitsType::ModelType > FluxType;
typedef typename SolversType::LinearOperatorType FullOperatorType;
typedef typename SolversType::LinearInverseOperatorType BasicLinearSolverType;
//----------- passes! ------------------------
typedef Dune::OperatorTraits< GridPartType, polynomialOrder, AnalyticalTraitsType,
DiscreteFunctionType, IndicatorType,
AdaptationHandlerType, ExtraParameterTuple > OperatorTraitsType;
typedef Dune::DGAdvectionDiffusionOperator< OperatorTraitsType > AssemblyOperatorType;
typedef Dune::DGPrimalMatrixAssembly< AssemblyOperatorType > AssembledOperatorType;
//#########################################################################
// --------- Operators using PASSES --------------------------
//============================================================
//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;
//static const int hasDiffusion = AnalyticalTraitsType::ModelType::hasDiffusion;
//typedef AdvectionDiffusionOperators< OperatorTraitsType, hasAdvection, hasDiffusion, _unlimited > AdvectionDiffusionOperatorType;
//typedef typename AdvectionDiffusionOperatorType::FullOperatorType FullOperatorType;
//typedef typename AdvectionDiffusionOperatorType::ImplicitOperatorType ImplicitOperatorType;
//typedef typename AdvectionDiffusionOperatorType::ExplicitOperatorType ExplicitOperatorType;
// --------- Operators using PASSES --------------------------
//============================================================
//---------Adaptivity ----------------------------------------------
// advection = true , diffusion = true
//typedef Dune :: DGAdaptationIndicatorOperator< OperatorTraitsType, hasAdvection, hasDiffusion > DGIndicatorType;
// gradient estimator
//typedef Estimator< DiscreteFunctionType, typename AnalyticalTraitsType::ProblemType > GradientIndicatorType ;
//typedef std::tuple< DGIndicatorType*, GradientIndicatorType* > IndicatorTupleType;
// --------Adaptivity ----------------------------------------------
//------------- Limiter ---------------------------------------------
typedef FullOperatorType LimiterOperatorType;
//------------ Limiter ---------------------------------------------
};
template <int polOrd>
struct Stepper
{
// this should be ok but could lead to a henn-egg problem
typedef EllipticAlgorithm< GridType, PoissonProblemCreator<GridType>, polynomialOrder > Type;
typedef Dune::Fem::EllipticAlgorithm< GridType, PoissonProblemCreator<GridType>, polOrd > Type;
// advection diffusion stepper using passes
//typedef Dune::Fem::AdvectionDiffusionStepper< GridTypeImp,
// AdvectionDiffusionProblemCreator<GridTypeImp>,
// polOrd,
// NoDataWriter,
// NoDiagnostics,
// NoAdaptivity,
// NoEOCWriter,
// NoCheckPointing,
// NoLimiting > StepperType1;
//algorithm from Tobias without passes
//typedef Dune::Fem::Tobias::EvolutionAlgorithm< GridTypeImp,
// Tobias::AdvectionDiffusionProblemCreator<GridTypeImp>,
// polOrd > StepperType2
//
// combined stepper
// typedef Dune::Fem::CombinedStepper< StepperType1, StepperType2,
// CombinedDataWriter,
// CombinedDiagnostics,
// CombinedAdaptivity,
// CombinedEOCWriter,
// CombinedChecPointing,
// CombinedLimiting >
};
};
#ifndef COMBINED_PROBLEM_CREATOR
......@@ -218,4 +497,4 @@ struct PoissonProblemCreator
#endif
#define NEW_STEPPER_SELECTOR_USED
#endif // FEMHOWTO_POISSONSTEPPER_HH
#endif // FEMHOWTO_HEATSTEPPER_HH
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment