Commit 9bd243b9 authored by Melanie Lipp's avatar Melanie Lipp Committed by Timo Koch

[test][freeflow] Move vtk output vectors from problem to main file

parent f5959801
...@@ -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;
...@@ -94,7 +151,6 @@ int main(int argc, char** argv) try ...@@ -94,7 +151,6 @@ int main(int argc, char** argv) try
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->updateTimeStepSize(timeLoop->timeStepSize()); problem->updateTimeStepSize(timeLoop->timeStepSize());
problem->createAnalyticalSolution();
// the solution vector // the solution vector
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>; using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
...@@ -113,9 +169,11 @@ int main(int argc, char** argv) try ...@@ -113,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
...@@ -144,8 +202,6 @@ int main(int argc, char** argv) try ...@@ -144,8 +202,6 @@ int main(int argc, char** argv) try
xOld = x; xOld = x;
gridVariables->advanceTimeStep(); gridVariables->advanceTimeStep();
problem->createAnalyticalSolution();
if (printL2Error) if (printL2Error)
{ {
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices; using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
...@@ -168,6 +224,7 @@ int main(int argc, char** argv) try ...@@ -168,6 +224,7 @@ int main(int argc, char** argv) try
// 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->updateTime(timeLoop->time());
analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
// write vtk output // write vtk output
vtkWriter.write(timeLoop->time()); vtkWriter.write(timeLoop->time());
......
...@@ -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,6 +102,8 @@ class AngeliTestProblem : public NavierStokesProblem<TypeTag> ...@@ -103,6 +102,8 @@ 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) : ParentType(fvGridGeometry)
{ {
...@@ -169,7 +170,7 @@ public: ...@@ -169,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_);
} }
/*! /*!
...@@ -183,7 +184,7 @@ public: ...@@ -183,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 + timeStepSize_; const Scalar t = time;
PrimaryVariables values; PrimaryVariables values;
...@@ -208,31 +209,7 @@ public: ...@@ -208,31 +209,7 @@ public:
*/ */
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{ {
return analyticalSolution(globalPos, -timeStepSize_); return analyticalSolution(globalPos, 0);
}
/*!
* \brief Returns the analytical solution for the pressure.
*/
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
{
return analyticalVelocityOnFace_;
} }
/*! /*!
...@@ -251,54 +228,12 @@ public: ...@@ -251,54 +228,12 @@ public:
timeStepSize_ = timeStepSize; timeStepSize_ = timeStepSize;
} }
/*!
* \brief Adds additional VTK output data to the VTKWriter.
*
* Function is called by the output module on every write.
*/
void createAnalyticalSolution()
{
analyticalPressure_.resize(this->fvGridGeometry().numCellCenterDofs());
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; static constexpr Scalar eps_ = 1e-6;
Scalar kinematicViscosity_; Scalar kinematicViscosity_;
Scalar time_ = 0; Scalar time_ = 0;
Scalar timeStepSize_ = 0; 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