Commit 0e52f00c authored by Timo Koch's avatar Timo Koch

Merge branch 'fix/angeli_postTimeStep' into 'master'

Fix calculation of L2 error in angeli test

See merge request !1411
parents 70c12576 9bd243b9
...@@ -27,12 +27,13 @@ ...@@ -27,12 +27,13 @@
#include <ctime> #include <ctime>
#include <iostream> #include <iostream>
#include <type_traits>
#include <tuple>
#include <dune/common/parallel/mpihelper.hh> #include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh> #include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh> #include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh> #include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
#include <dumux/assembly/staggeredfvassembler.hh> #include <dumux/assembly/staggeredfvassembler.hh>
#include <dumux/assembly/diffmethod.hh> #include <dumux/assembly/diffmethod.hh>
...@@ -47,6 +48,62 @@ ...@@ -47,6 +48,62 @@
#include "problem.hh" #include "problem.hh"
/*!
* \brief Creates analytical solution.
* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
* \param time the time at which to evaluate the analytical solution
* \param problem the problem for which to evaluate the analytical solution
*/
template<class Scalar, class Problem>
auto createAnalyticalSolution(const Scalar time, const Problem& problem)
{
const auto& fvGridGeometry = problem.fvGridGeometry();
using GridView = typename std::decay_t<decltype(fvGridGeometry)>::GridView;
static constexpr auto dim = GridView::dimension;
static constexpr auto dimWorld = GridView::dimensionworld;
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
std::vector<Scalar> analyticalPressure;
std::vector<VelocityVector> analyticalVelocity;
std::vector<VelocityVector> analyticalVelocityOnFace;
analyticalPressure.resize(fvGridGeometry.numCellCenterDofs());
analyticalVelocity.resize(fvGridGeometry.numCellCenterDofs());
analyticalVelocityOnFace.resize(fvGridGeometry.numFaceDofs());
using Indices = typename Problem::Indices;
for (const auto& element : elements(fvGridGeometry.gridView()))
{
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
auto ccDofIdx = scv.dofIndex();
auto ccDofPosition = scv.dofPosition();
auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition, time);
// velocities on faces
for (auto&& scvf : scvfs(fvGeometry))
{
const auto faceDofIdx = scvf.dofIndex();
const auto faceDofPosition = scvf.center();
const auto dirIdx = scvf.directionIndex();
const auto analyticalSolutionAtFace = problem.analyticalSolution(faceDofPosition, time);
analyticalVelocityOnFace[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
for(int dirIdx = 0; dirIdx < dim; ++dirIdx)
analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
}
}
return std::make_tuple(analyticalPressure, analyticalVelocity, analyticalVelocityOnFace);
}
int main(int argc, char** argv) try int main(int argc, char** argv) try
{ {
using namespace Dumux; using namespace Dumux;
...@@ -93,7 +150,7 @@ int main(int argc, char** argv) try ...@@ -93,7 +150,7 @@ int main(int argc, char** argv) try
// the problem (initial and boundary conditions) // the problem (initial and boundary conditions)
using Problem = GetPropType<TypeTag, Properties::Problem>; using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(fvGridGeometry); auto problem = std::make_shared<Problem>(fvGridGeometry);
problem->setTimeLoop(timeLoop); problem->updateTimeStepSize(timeLoop->timeStepSize());
// the solution vector // the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
...@@ -112,9 +169,11 @@ int main(int argc, char** argv) try ...@@ -112,9 +169,11 @@ int main(int argc, char** argv) try
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name()); StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using IOFields = GetPropType<TypeTag, Properties::IOFields>; using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter); // Add model specific output fields IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact");
vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact"); auto analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact"); vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
vtkWriter.addFaceField(std::get<2>(analyticalSolution), "faceVelocityExact");
vtkWriter.write(0.0); vtkWriter.write(0.0);
// the assembler with time loop for instationary problem // the assembler with time loop for instationary problem
...@@ -130,6 +189,7 @@ int main(int argc, char** argv) try ...@@ -130,6 +189,7 @@ int main(int argc, char** argv) try
NewtonSolver nonLinearSolver(assembler, linearSolver); NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop // time loop
const bool printL2Error = getParam<bool>("Problem.PrintL2Error");
timeLoop->start(); do timeLoop->start(); do
{ {
// set previous solution for storage evaluations // set previous solution for storage evaluations
...@@ -142,12 +202,29 @@ int main(int argc, char** argv) try ...@@ -142,12 +202,29 @@ int main(int argc, char** argv) try
xOld = x; xOld = x;
gridVariables->advanceTimeStep(); gridVariables->advanceTimeStep();
problem->createAnalyticalSolution(); if (printL2Error)
{
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
const auto l2error = L2Error::calculateL2Error(*problem, x);
const int numCellCenterDofs = fvGridGeometry->numCellCenterDofs();
const int numFaceDofs = fvGridGeometry->numFaceDofs();
std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
<< std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
<< std::scientific
<< "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
<< ", L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
<< ", L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
<< std::endl;
}
// advance to the time loop to the next step // advance to the time loop to the next step
timeLoop->advanceTimeStep(); timeLoop->advanceTimeStep();
problem->updateTime(timeLoop->time());
problem->postTimeStep(x); analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
// write vtk output // write vtk output
vtkWriter.write(timeLoop->time()); vtkWriter.write(timeLoop->time());
...@@ -157,6 +234,7 @@ int main(int argc, char** argv) try ...@@ -157,6 +234,7 @@ int main(int argc, char** argv) try
// set new dt as suggested by newton solver // set new dt as suggested by newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
problem->updateTimeStepSize(timeLoop->timeStepSize());
} while (!timeLoop->finished()); } while (!timeLoop->finished());
......
...@@ -89,7 +89,6 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag> ...@@ -89,7 +89,6 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag>
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>; using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>; using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>; using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>; using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
...@@ -103,41 +102,14 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag> ...@@ -103,41 +102,14 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag>
using VelocityVector = Dune::FieldVector<Scalar, dimWorld>; using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
public: public:
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
AngeliTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) AngeliTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
: ParentType(fvGridGeometry), eps_(1e-6) : ParentType(fvGridGeometry)
{ {
printL2Error_ = getParam<bool>("Problem.PrintL2Error");
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0); kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
} }
/*!
* \name Problem parameters
*/
// \{
bool shouldWriteRestartFile() const
{
return false;
}
void postTimeStep(const SolutionVector& curSol) const
{
if(printL2Error_)
{
using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
const auto l2error = L2Error::calculateL2Error(*this, curSol);
const int numCellCenterDofs = this->fvGridGeometry().numCellCenterDofs();
const int numFaceDofs = this->fvGridGeometry().numFaceDofs();
std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
<< std::setw(6) << numCellCenterDofs << " cc dofs and " << numFaceDofs << " face dofs (total: " << numCellCenterDofs + numFaceDofs << "): "
<< std::scientific
<< "L2(p) = " << l2error.first[Indices::pressureIdx] << " / " << l2error.second[Indices::pressureIdx]
<< ", L2(vx) = " << l2error.first[Indices::velocityXIdx] << " / " << l2error.second[Indices::velocityXIdx]
<< ", L2(vy) = " << l2error.first[Indices::velocityYIdx] << " / " << l2error.second[Indices::velocityYIdx]
<< std::endl;
}
}
/*! /*!
* \brief Returns the temperature within the domain in [K]. * \brief Returns the temperature within the domain in [K].
* *
...@@ -146,17 +118,6 @@ public: ...@@ -146,17 +118,6 @@ public:
Scalar temperature() const Scalar temperature() const
{ return 298.0; } { return 298.0; }
/*!
* \brief Returns the sources within the domain.
*
* \param globalPos The global position
*/
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{
return NumEqVector(0.0);
}
// \} // \}
/*! /*!
* \name Boundary conditions * \name Boundary conditions
...@@ -209,7 +170,7 @@ public: ...@@ -209,7 +170,7 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const PrimaryVariables dirichletAtPos(const GlobalPosition & globalPos) const
{ {
// use the values of the analytical solution // use the values of the analytical solution
return analyticalSolution(globalPos, time()); return analyticalSolution(globalPos, time_+timeStepSize_);
} }
/*! /*!
...@@ -223,7 +184,7 @@ public: ...@@ -223,7 +184,7 @@ public:
const Scalar x = globalPos[0]; const Scalar x = globalPos[0];
const Scalar y = globalPos[1]; const Scalar y = globalPos[1];
const Scalar t = time + timeLoop_->timeStepSize(); const Scalar t = time;
PrimaryVariables values; PrimaryVariables values;
...@@ -248,93 +209,31 @@ public: ...@@ -248,93 +209,31 @@ public:
*/ */
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{ {
return analyticalSolution(globalPos, -timeLoop_->timeStepSize()); return analyticalSolution(globalPos, 0);
} }
/*! /*!
* \brief Returns the analytical solution for the pressure. * \brief Updates the time
*/
auto& getAnalyticalPressureSolution() const
{
return analyticalPressure_;
}
/*!
* \brief Returns the analytical solution for the velocity.
*/
auto& getAnalyticalVelocitySolution() const
{
return analyticalVelocity_;
}
/*!
* \brief Returns the analytical solution for the velocity at the faces.
*/ */
auto& getAnalyticalVelocitySolutionOnFace() const void updateTime(const Scalar time)
{
return analyticalVelocityOnFace_;
}
void setTimeLoop(TimeLoopPtr timeLoop)
{ {
timeLoop_ = timeLoop; time_ = time;
createAnalyticalSolution();
} }
Scalar time() const /*!
{ * \brief Updates the time step size
return timeLoop_->time();
}
/*!
* \brief Adds additional VTK output data to the VTKWriter.
*
* Function is called by the output module on every write.
*/ */
void createAnalyticalSolution() void updateTimeStepSize(const Scalar timeStepSize)
{ {
analyticalPressure_.resize(this->fvGridGeometry().numCellCenterDofs()); timeStepSize_ = timeStepSize;
analyticalVelocity_.resize(this->fvGridGeometry().numCellCenterDofs());
analyticalVelocityOnFace_.resize(this->fvGridGeometry().numFaceDofs());
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
for (auto&& scv : scvs(fvGeometry))
{
auto ccDofIdx = scv.dofIndex();
auto ccDofPosition = scv.dofPosition();
auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition, time());
// velocities on faces
for (auto&& scvf : scvfs(fvGeometry))
{
const auto faceDofIdx = scvf.dofIndex();
const auto faceDofPosition = scvf.center();
const auto dirIdx = scvf.directionIndex();
const auto analyticalSolutionAtFace = analyticalSolution(faceDofPosition, time());
analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
}
analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
for(int dirIdx = 0; dirIdx < ModelTraits::dim(); ++dirIdx)
analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
}
}
} }
private: private:
static constexpr Scalar eps_ = 1e-6;
Scalar eps_;
Scalar kinematicViscosity_; Scalar kinematicViscosity_;
TimeLoopPtr timeLoop_; Scalar time_ = 0;
bool printL2Error_; Scalar timeStepSize_ = 0;
std::vector<Scalar> analyticalPressure_;
std::vector<VelocityVector> analyticalVelocity_;
std::vector<VelocityVector> analyticalVelocityOnFace_;
}; };
} // end namespace Dumux } // end namespace Dumux
......
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