Commit 51135416 authored by Melanie Lipp's avatar Melanie Lipp
Browse files

[freeflow][analyticalsolution] Separate printError from building analytical...

[freeflow][analyticalsolution] Separate printError from building analytical solution vectors for vtk. Get rid of direct TypeTag-dependency.
parent 7c21b79d
Pipeline #7694 waiting for manual action with stages
......@@ -39,38 +39,15 @@
#include "errors.hh"
namespace Dumux {
/*!
* \ingroup NavierStokesModel
* \brief Helps to deal with analytical solutions
*/
template<class TypeTag>
class NavierStokesAnalyticalSolutionFacility
template<class Scalar, class Problem, class GridGeometry, class PrimaryVariables, class Indices>
class NavierStokesErrorTmp
{
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using GridView = typename GridGeometry::GridView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>;
using CellCenterSolutionVector = GetPropType<TypeTag, Properties::CellCenterSolutionVector>;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
public:
NavierStokesAnalyticalSolutionFacility(std::shared_ptr<const Problem> problem)
NavierStokesErrorTmp(std::shared_ptr<const Problem> problem)
: problem_(problem)
{
name_ = getParam<std::string>("Problem.Name");
......@@ -81,50 +58,9 @@ public:
logFile << ",,,,,abs L2,,,rel L2,,,Linf,,,"<< std::endl;
logFile << "time/refinement step,corresponding paraview files,cc dofs, face dofs, all dofs, p,u,v,p,u,v,p,u,v" << std::endl;
logFile.close();
createAnalyticalSolution();
}
/*!
* \brief Return the analytical solution of the problem at a scvf
*/
Scalar instationaryAnalyticalVelocitySolutionAtPos(const SubControlVolumeFace& scvf, const Scalar time) const
{
unsigned int dirIdx = scvf.directionIndex();
return problem_->instationaryAnalyticalSolution(scvf.center(), time)[Indices::velocity(dirIdx)];
}
/*!
* \brief Return the analytical solution of the problem at a scvf
*/
Scalar instationaryAnalyticalNonDirVelocitySolutionAtPos(const SubControlVolumeFace& scvf, const Scalar time) const
{
unsigned int nonDirIdx = (scvf.directionIndex() == 0)?1:0;
return problem_->instationaryAnalyticalSolution(scvf.center(), time)[Indices::velocity(nonDirIdx)];
}
/*!
* \brief Return the analytical solution of the problem at a given position
*
* \param globalPos The global position
*/
Scalar analyticalVelocitySolutionAtPos(const SubControlVolumeFace& scvf) const
{
unsigned int dirIdx = scvf.directionIndex();
return problem_->analyticalSolution(scvf.center())[Indices::velocity(dirIdx)]; ;
}
/*!
* \brief Return the analytical solution of the problem at a given position
*
* \param globalPos The global position
*/
Scalar analyticalNonDirVelocitySolutionAtPos(const SubControlVolumeFace& scvf) const
{
unsigned int nonDirIdx = (scvf.directionIndex() == 0)?1:0;
return problem_->analyticalSolution(scvf.center())[Indices::velocity(nonDirIdx)];
}
template<class SolutionVector>
void printErrors(const SolutionVector& curSol) const
{
PrimaryVariables l2NormAbs(0.0);
......@@ -136,18 +72,19 @@ public:
logFile << ", 00001.vtu ,";
logFile.close();
NavierStokesTestErrors<TypeTag>::calculateErrors(
NavierStokesTestErrors<Scalar, GridGeometry, PrimaryVariables, Indices>::calculateErrors(
l2NormAbs,
l2NormRel,
lInfinityNorm,
(problem_->gridGeometry()),
curSol,
[&](const SubControlVolume& scv) { return problem_->analyticalSolution(scv.dofPosition())[Indices::pressureIdx]; },
[&](const SubControlVolumeFace& scvf) { return analyticalVelocitySolutionAtPos(scvf);});
[&](const SubControlVolumeFace& scvf) { return problem_->analyticalSolution(scvf.center())[Indices::velocity(scvf.directionIndex())];});
NavierStokesTestErrors<TypeTag>::printErrors(l2NormAbs, l2NormRel, lInfinityNorm, (problem_->gridGeometry()), name_);
NavierStokesTestErrors<Scalar, GridGeometry, PrimaryVariables, Indices>::printErrors(l2NormAbs, l2NormRel, lInfinityNorm, (problem_->gridGeometry()), name_);
}
template<class SolutionVector>
void printErrors(const SolutionVector& curSol, Scalar time, unsigned int timeIndex) const
{
PrimaryVariables l2NormAbs(0.0);
......@@ -159,21 +96,48 @@ public:
logFile << time << " , 0000" << timeIndex+1 << ".vtu, ";
logFile.close();
NavierStokesTestErrors<TypeTag>::calculateErrors(
NavierStokesTestErrors<Scalar, GridGeometry, PrimaryVariables, Indices>::calculateErrors(
l2NormAbs,
l2NormRel,
lInfinityNorm,
(problem_->gridGeometry()),
curSol,
[&](const SubControlVolume& scv) { return problem_->instationaryAnalyticalSolution(scv.dofPosition(), time)[Indices::pressureIdx]; },
[&](const SubControlVolumeFace& scvf) { return instationaryAnalyticalVelocitySolutionAtPos(scvf, time);});
[&](const SubControlVolumeFace& scvf) { return problem_->instationaryAnalyticalSolution(scvf.center(), time)[Indices::velocity(scvf.directionIndex())];});
NavierStokesTestErrors<Scalar, GridGeometry, PrimaryVariables, Indices>::printErrors(l2NormAbs, l2NormRel, lInfinityNorm, (problem_->gridGeometry()), name_);
}
private:
std::shared_ptr<const Problem> problem_;
std::string name_;
};
/*!
* \ingroup NavierStokesModel
* \brief Helps to deal with analytical solutions
*/
template<class Scalar, class FaceSolutionVector, class CellCenterSolutionVector, class GridGeometry, class Problem, class Indices>
class NavierStokesAnalyticalSolutionFacility
{
using GridView = typename GridGeometry::GridView;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
NavierStokesTestErrors<TypeTag>::printErrors(l2NormAbs, l2NormRel, lInfinityNorm, (problem_->gridGeometry()), name_);
public:
NavierStokesAnalyticalSolutionFacility(std::shared_ptr<const Problem> problem)
: problem_(problem)
{
createAnalyticalSolution();//TODO figure out for instationary
}
void createAnalyticalSolution(Scalar time)
{
auto instationaryAnalyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return instationaryAnalyticalVelocitySolutionAtPos(scvf, time); };
auto instationaryAnalyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return problem_->instationaryAnalyticalSolution(scvf.center(), time)[Indices::velocity(scvf.directionIndex())]; };
createAnalyticalSolution_(
[&](const SubControlVolume& scv) { return problem_->instationaryAnalyticalSolution(scv.dofPosition(), time)[Indices::pressureIdx]; },
instationaryAnalyticalVelocitySolutionAtPosLambdaFunction,
......@@ -182,7 +146,7 @@ public:
void createAnalyticalSolution()
{
auto analyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return analyticalVelocitySolutionAtPos(scvf); };
auto analyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return problem_->analyticalSolution(scvf.center())[Indices::velocity(scvf.directionIndex())]; };
createAnalyticalSolution_(
[&](const SubControlVolume& scv) { return problem_->analyticalSolution(scv.dofPosition())[Indices::pressureIdx]; },
analyticalVelocitySolutionAtPosLambdaFunction,
......@@ -191,7 +155,7 @@ public:
auto getAnalyticalSolution(Scalar time) const
{
auto instationaryAnalyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return instationaryAnalyticalVelocitySolutionAtPos(scvf, time); };
auto instationaryAnalyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return problem_->instationaryAnalyticalSolution(scvf.center(), time)[Indices::velocity(scvf.directionIndex())]; };
return getAnalyticalSolution_(
[&](const SubControlVolume& scv) { return problem_->instationaryAnalyticalSolution(scv.dofPosition(), time)[Indices::pressureIdx]; },
instationaryAnalyticalVelocitySolutionAtPosLambdaFunction,
......@@ -200,7 +164,7 @@ public:
auto getAnalyticalSolution() const
{
auto analyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return analyticalVelocitySolutionAtPos(scvf); };
auto analyticalVelocitySolutionAtPosLambdaFunction = [&](const SubControlVolumeFace& scvf) { return problem_->analyticalSolution(scvf.center())[Indices::velocity(scvf.directionIndex())]; };
return getAnalyticalSolution_(
[&](const SubControlVolume& scv) { return problem_->analyticalSolution(scv.dofPosition())[Indices::pressureIdx]; },
analyticalVelocitySolutionAtPosLambdaFunction,
......@@ -337,8 +301,6 @@ private:
std::vector<VelocityVector> analyticalVelocity_;
std::vector<VelocityVector> analyticalVelocityOnFace_;
FaceSolutionVector analyticalFaceSolutionVector_;
std::string name_;
};
} // end namespace Dumux
......
......@@ -84,7 +84,12 @@ int main(int argc, char** argv)
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
const Dumux::NavierStokesAnalyticalSolutionFacility<TypeTag> analyticalSolFacility(problem);
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FaceSolutionVector = GetPropType<TypeTag, Properties::FaceSolutionVector>;
using CellCenterSolutionVector = GetPropType<TypeTag, Properties::CellCenterSolutionVector>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
const Dumux::NavierStokesAnalyticalSolutionFacility<Scalar, FaceSolutionVector, CellCenterSolutionVector, GridGeometry, Problem, Indices> analyticalSolFacility(problem);
// the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
......@@ -128,7 +133,12 @@ int main(int argc, char** argv)
// write vtk output
if (getParam<bool>("Problem.PrintErrors", false))
analyticalSolFacility.printErrors(x);
{
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
const Dumux::NavierStokesErrorTmp<Scalar, Problem, GridGeometry, PrimaryVariables, Indices> error(problem);
error.printErrors(x);
}
vtkWriter.write(1.0);
......
......@@ -57,10 +57,8 @@ class DoneaTestProblem : public NavierStokesProblem<TypeTag>
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
public:
DoneaTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
......
......@@ -33,14 +33,9 @@ namespace Dumux {
* \ingroup NavierStokesTests
* \brief Routines to calculate the discrete L2 error
*/
template<class TypeTag>
template<class Scalar, class GridGeometry, class PrimaryVariables, class Indices>
class NavierStokesTestErrors
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using Indices = typename ModelTraits::Indices;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using Extrusion = Extrusion_t<GridGeometry>;
......@@ -52,7 +47,7 @@ public:
* \param problem The object specifying the problem which ought to be simulated
* \param curSol Vector containing the current solution
*/
template<class GridGeometry, class SolutionVector, class LambdaA, class LambdaB>
template<class SolutionVector, class LambdaA, class LambdaB>
static void calculateErrors(PrimaryVariables& l2NormAbs,
PrimaryVariables& l2NormRel,
PrimaryVariables& lInfinityNorm,
......@@ -115,7 +110,6 @@ public:
lInfinityNorm[Indices::velocityYIdx] = maxVError;
}
template<class GridGeometry>
static void printErrors(const PrimaryVariables& l2NormAbs,
const PrimaryVariables& l2NormRel,
const PrimaryVariables& lInfinityNorm,
......@@ -146,8 +140,10 @@ private:
template<class SolutionVector, class LambdaA>
static void processPressure_(const SolutionVector& curSol, const LambdaA& instationaryAnalyticalPressureSolutionAtPosLambdaFunction, const SubControlVolume& scv, Scalar& sumPReference, Scalar& sumPError, Scalar& maxPError)
{
static constexpr auto dim = GridGeometry::GridView::dimension;
const auto analyticalSolutionCellCenter = instationaryAnalyticalPressureSolutionAtPosLambdaFunction(scv);
const auto numericalSolutionCellCenter = curSol[GridGeometry::cellCenterIdx()][scv.dofIndex()][Indices::pressureIdx - ModelTraits::dim()];
const auto numericalSolutionCellCenter = curSol[GridGeometry::cellCenterIdx()][scv.dofIndex()][Indices::pressureIdx - dim];
const Scalar pError = absDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter);
const Scalar pReference = absDiff_(analyticalSolutionCellCenter, 0.0);
......
Markdown is supported
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