Skip to content
Snippets Groups Projects
Commit 49f43128 authored by Robert K's avatar Robert K
Browse files

[cleanup] Unified LDG and BR1 fluxes.

parent ad5b7464
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
dune/fem-dg/operator/fluxes/diffusion/averageflux.hh deleted 100644 → 0
+ 0
− 355
View file @ ad5b7464
#ifndef DUNE_FEM_DG_LDGAVERAGEFLUX_HH
#define DUNE_FEM_DG_LDGAVERAGEFLUX_HH
// local includes
#include "fluxbase.hh"
namespace Dune
{
namespace Fem
{
/**
* \brief A diffusion flux for the LDG scheme.
*
* \ingroup DiffusionFluxes
*/
template <class DiscreteFunctionSpaceImp,
class ModelImp,
class FluxParameterImp >
class LDGAverageDiffusionFlux :
public LDGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp >
{
typedef LDGDiffusionFluxBase< DiscreteFunctionSpaceImp, ModelImp, FluxParameterImp > BaseType;
public:
typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
enum { dimDomain = DiscreteFunctionSpaceType :: dimDomain };
enum { dimRange = DiscreteFunctionSpaceType :: dimRange };
typedef typename DiscreteFunctionSpaceType :: DomainType DomainType;
typedef typename DiscreteFunctionSpaceType :: RangeFieldType RangeFieldType;
typedef typename DiscreteFunctionSpaceType :: DomainFieldType DomainFieldType;
typedef typename DiscreteFunctionSpaceType :: RangeType RangeType;
typedef typename DiscreteFunctionSpaceType :: JacobianRangeType JacobianRangeType;
typedef FieldVector< DomainFieldType, dimDomain-1 > FaceDomainType;
typedef typename DiscreteFunctionSpaceType :: GridPartType GridPartType;
typedef typename GridPartType :: IntersectionIteratorType IntersectionIterator;
typedef typename IntersectionIterator :: Intersection Intersection;
typedef typename GridPartType :: GridType GridType;
typedef typename DiscreteFunctionSpaceType :: EntityType EntityType;
enum { dimGradRange = dimDomain * dimRange };
enum { polOrd = DiscreteFunctionSpaceType :: polynomialOrder };
typedef typename BaseType :: ParameterType ParameterType;
// type of gradient space
typedef typename DiscreteFunctionSpaceType ::
template ToNewDimRange< dimGradRange > :: Type DiscreteGradientSpaceType;
typedef typename DiscreteGradientSpaceType :: RangeType GradientRangeType;
typedef typename DiscreteGradientSpaceType :: JacobianRangeType GradientJacobianType;
// jacobians of the functions do not have to be evaluated for this flux
enum { evaluateJacobian = false };
protected:
using BaseType::determineDirection;
using BaseType::model_;
using BaseType::cflDiffinv_;
using BaseType::numericalFlux ;
using BaseType::parameter ;
using BaseType::nonconformingFactor_;
using BaseType::dimensionFactor_;
public:
/**
* \brief constructor
*/
LDGAverageDiffusionFlux(GridPartType& gridPart,
const ModelImp& mod,
const ParameterType& param ) :
BaseType( gridPart, mod, param ),
penalty_( parameter().penalty() ),
// Set CFL number for penalty term (compare diffusion in first pass)
penaltyTerm_( std::abs( penalty_ ) > 0 )
{
if( Fem::Parameter::verbose () )
{
std::cout << "LDGAverageDiffusionFlux: penalty = " << penalty_ << std::endl;
}
}
LDGAverageDiffusionFlux(const LDGAverageDiffusionFlux& other)
: BaseType( other ),
penalty_( other.penalty_ ),
penaltyTerm_( other.penaltyTerm_ )
{}
//! returns true if lifting has to be calculated
const bool hasLifting () const { return false; }
protected:
enum { realLDG = true };
double theta( const DomainType& normal ) const
{
if( realLDG )
{
// LDG thete is 1 or 0
if( determineDirection( normal ) )
return 1.0;
else
return 0.0;
}
else
// Average fluxes
return 0.5;
}
public:
/**
* \brief flux function on interfaces between cells
*
* \param left local evaluation
* \param right local evaluation
* \param uLeft DOF evaluation on this side of \c intersection
* \param uRight DOF evaluation on the other side of \c intersection
* \param gLeft result for this side of \c intersection
* \param gRight result for the other side of \c intersection
* \return wave speed estimate (multiplied with the integration element of the intersection).
* To estimate the time step |T|/wave is used
*/
template <class LocalEvaluation>
double gradientNumericalFlux(const LocalEvaluation& left,
const LocalEvaluation& right,
const RangeType& uLeft,
const RangeType& uRight,
GradientRangeType& gLeft,
GradientRangeType& gRight,
GradientJacobianType& gDiffLeft,
GradientJacobianType& gDiffRight) const
{
const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
// get factor for each side
const double thetaLeft = theta( normal );
const double thetaRight = 1.0 - thetaLeft;
GradientJacobianType diffmatrix;
double diffTimeStep = 0.0;
// select left or right or average
if( thetaLeft > 0 )
{
//TODO use diffusionTimeStep
diffTimeStep = 0;
/* central differences (might be suboptimal) */
model_.jacobian(left, uLeft, diffmatrix );
diffmatrix.mv(normal, gLeft );
gLeft *= thetaLeft ;
}
else
gLeft = 0;
if( thetaRight > 0 )
{
//const double diffStepRight = 0;
// right jacobian
model_.jacobian( right, uRight, diffmatrix );
diffmatrix.mv(normal, gRight);
// add to flux
gLeft.axpy( thetaRight, gRight );
// diffTimeStep = std::max( diffTimeStep, diffStepRight );
}
// copy flux
gRight = gLeft;
#ifndef NDEBUG
gDiffLeft = 0;
gDiffRight = 0;
#endif
// upper bound for the next time step length
return diffTimeStep; // * cflDiffinv_;
}
template <class LocalEvaluation>
double gradientBoundaryFlux(const LocalEvaluation& left,
const RangeType& uLeft,
const RangeType& uBnd,
GradientRangeType& gLeft,
GradientJacobianType& gDiffLeft) const
{
const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
// get factor for each side
const double thetaLeft = theta( normal );
const double thetaRight = 1.0 - thetaLeft;
GradientJacobianType diffmatrix;
// calculate uVal
RangeType uVal( 0 );
// select u from uLeft and uRight
if( thetaLeft > 0 )
uVal.axpy( thetaLeft , uLeft );
if( thetaRight > 0 )
uVal.axpy( thetaRight, uBnd );
// compute jacobian of u
model_.jacobian(left, uVal, diffmatrix );
diffmatrix.mv(normal, gLeft);
#ifndef NDEBUG
gDiffLeft = 0;
#endif
//TODO use diffusionTimeStep
const double diffTimeStep = 0; // * cflDiffinv_;
return diffTimeStep;
}
/**
* \brief flux function on interfaces between cells
*
* \param left local evaluation
* \param right local evaluation
* \param uLeft DOF evaluation on this side of \c intersection
* \param uRight DOF evaluation on the other side of \c intersection
* \param gLeft result for this side of \c intersection
* \param gRight result for the other side of \c intersection
* \return wave speed estimate (multiplied with the integration element of the intersection).
* To estimate the time step |T|/wave is used
*/
template <class LocalEvaluation>
double numericalFlux(const LocalEvaluation& left,
const LocalEvaluation& right,
const RangeType& uLeft,
const RangeType& uRight,
const JacobianRangeType& sigmaLeft,
const JacobianRangeType& sigmaRight,
RangeType& gLeft,
RangeType& gRight,
JacobianRangeType& gDiffLeft, // not used here (only for primal passes)
JacobianRangeType& gDiffRight )
{
const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
const double faceLengthSqr = normal.two_norm2();
/**********************************
* Diffusion sigma Flux (Pass 2) *
**********************************/
JacobianRangeType diffmatrix;
// diffusion for uLeft
model_.diffusion( left, uLeft, sigmaLeft, diffmatrix);
RangeType diffflux;
diffmatrix.mv(normal, diffflux);
// diffusion for uRight
model_.diffusion( right, uRight, sigmaRight, diffmatrix);
// add to diffflux
diffmatrix.umv(normal, diffflux);
diffflux *= 0.5;
const double faceVolumeEstimate = dimensionFactor_ *
(left.intersection().conforming() ? faceLengthSqr
: (nonconformingFactor_ * faceLengthSqr));
const double diffTimeLeft =
model_.diffusionTimeStep( left, faceVolumeEstimate, uLeft );
const double diffTimeRight =
model_.diffusionTimeStep( right, faceVolumeEstimate, uRight );
double diffTimeStep = std::max( diffTimeLeft, diffTimeRight );
// add penalty factor
const double factor = penalty_ * diffTimeStep ;
RangeType jump( uLeft );
jump -= uRight;
diffflux.axpy(factor, jump);
gLeft = diffflux;
gRight = diffflux;
#ifndef NDEBUG
gDiffLeft = 0;
gDiffRight = 0;
#endif
// timestep restict to diffusion timestep
// WARNING: reconsider this
diffTimeStep *= cflDiffinv_;
return diffTimeStep;
}
/**
* \brief same as numericalFlux() but for fluxes over boundary interfaces
*/
template <class LocalEvaluation>
double boundaryFlux(const LocalEvaluation& left,
const RangeType& uLeft,
const RangeType& uRight,
const JacobianRangeType& sigmaLeft,
RangeType& gLeft,
JacobianRangeType& gDiffLeft )
{
// get local point
const DomainType normal = left.intersection().integrationOuterNormal( left.localPosition() );
const double faceLengthSqr = normal.two_norm2();
/****************************/
/* Diffusion (Pass 2) */
/****************************/
JacobianRangeType diffmatrix;
model_.diffusion(left, uLeft, sigmaLeft, diffmatrix);
diffmatrix.mv(normal, gLeft);
const double faceVolumeEstimate = dimensionFactor_ * faceLengthSqr;
double diffTimeStep =
model_.diffusionTimeStep( left, faceVolumeEstimate, uLeft );
// add penalty term
const double factor = penalty_ * diffTimeStep;
RangeType jump( uLeft );
jump -= uRight;
gLeft.axpy(factor,jump);
#ifndef NDEBUG
gDiffLeft = 0;
#endif
diffTimeStep *= cflDiffinv_;
return diffTimeStep;
}
protected:
const double penalty_;
const bool penaltyTerm_;
}; // end LDGAverageDiffusionFlux
} // end namespace
} // end namespace
#endif
dune/fem-dg/operator/fluxes/diffusion/fluxes.hh
+ 27
− 11
View file @ 49f43128
......@@ -4,7 +4,6 @@
#include "fluxbase.hh"
#include "dgprimalfluxes.hh"
#include "ldgflux.hh"
#include "averageflux.hh"
namespace Dune
{
......@@ -227,9 +226,9 @@ namespace Fem
template <class DiscreteFunctionSpaceImp,
class Model>
class DGDualDiffusionFlux< DiscreteFunctionSpaceImp, Model, DualDiffusionFlux::Enum::average >
: public LDGAverageDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
: public LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
{
typedef LDGAverageDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
typedef LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
BaseType;
public:
......@@ -243,7 +242,7 @@ namespace Fem
DGDualDiffusionFlux( GridPartType& gridPart,
const Model& model,
const ParameterType& parameters = ParameterType() )
: BaseType( gridPart, model, parameters )
: BaseType( gridPart, model, parameters, DualDiffusionFlux::Enum::average )
{
}
};
......@@ -251,9 +250,9 @@ namespace Fem
template <class DiscreteFunctionSpaceImp,
class Model>
class DGDualDiffusionFlux< DiscreteFunctionSpaceImp, Model, DualDiffusionFlux::Enum::ldg >
: public LDGDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
: public LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
{
typedef LDGDiffusionFlux< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
typedef LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
BaseType;
public:
......@@ -267,18 +266,35 @@ namespace Fem
DGDualDiffusionFlux( GridPartType& gridPart,
const Model& model,
const ParameterType& parameters = ParameterType() )
: BaseType( gridPart, model, parameters )
: BaseType( gridPart, model, parameters, DualDiffusionFlux::Enum::ldg )
{
}
};
//TODO implement general case...
template <class DiscreteFunctionSpaceImp,
class Model>
class DGDualDiffusionFlux< DiscreteFunctionSpaceImp, Model, DualDiffusionFlux::Enum::general >
{};
: public LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
{
typedef LDGDiffusionFluxImpl< DiscreteFunctionSpaceImp, Model, DGDualDiffusionFluxParameters >
BaseType;
public:
typedef DiscreteFunctionSpaceImp DiscreteFunctionSpaceType;
typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
typedef typename BaseType::ParameterType ParameterType;
/**
* \brief constructor reading parameters
*/
DGDualDiffusionFlux( GridPartType& gridPart,
const Model& model,
const ParameterType& parameters = ParameterType() )
: BaseType( gridPart, model, parameters, DualDiffusionFlux::Enum::general )
{
}
};
} // end namespace
} // end namespace
} // end namespace Fem
} // end namespace Dune
#endif
dune/fem-dg/operator/fluxes/diffusion/ldgflux.hh
+ 187
− 271
View file @ 49f43128
This diff is collapsed.
dune/fem-dg/operator/fluxes/diffusion/parameters.hh
+ 4
− 4
View file @ 49f43128
......@@ -40,7 +40,7 @@ namespace Fem
//! Contains all known enums for primal diffusion fluxes which can be chosen via parameter file.
const Enum _enums[] = { Enum::cdg2, Enum::cdg, Enum::br2, Enum::ip, Enum::nipg, Enum::bo };
//! Contains all known names of primal diffusion fluxes which can be chosen via parameter file.
//! Contains all known names of primal DG diffusion fluxes which can be chosen via parameter file.
const std::string _strings[] = { "CDG2", "CDG" , "BR2", "IP" , "NIPG", "BO" };
//! Number of known primal diffusion fluxes which can be chosen via parameter file.
static const int _size = 6;
......@@ -50,15 +50,15 @@ namespace Fem
namespace DualDiffusionFlux
{
/**
* \brief Enum of all known primal diffusion flux implementations.
* \brief Enum of all known local DG diffusion flux implementations.
*
* \ingroup FemDGParameter
*/
enum class Enum
{
//! average theta flux.
//! Bassi-Rebay flux.
average,
//! ldg flux.
//! Local DG flux.
ldg,
//! general flux: Parameter selection is done via parameter file!
general,
......
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