From 131e84ecf392ca1f4c75c2af266e6b176e0b7cff Mon Sep 17 00:00:00 2001
From: Melanie Lipp <melanie.lipp@iws.uni-stuttgart.de>
Date: Wed, 25 Aug 2021 12:53:10 +0200
Subject: [PATCH] [test][freeflow] Move createAnalyticalSolution to file and
reuse for several tests. Change signature of analyticalSolution to include
time for angeli test.
---
.../navierstokes/analyticalsolutionvectors.hh | 158 ++++++++++++++++++
test/freeflow/navierstokes/angeli/main.cc | 74 ++------
test/freeflow/navierstokes/angeli/problem.hh | 29 ++--
test/freeflow/navierstokes/channel/1d/main.cc | 17 +-
.../navierstokes/channel/1d/problem.hh | 67 --------
test/freeflow/navierstokes/channel/2d/main.cc | 17 +-
.../navierstokes/channel/2d/problem.hh | 69 --------
test/freeflow/navierstokes/donea/main.cc | 17 +-
test/freeflow/navierstokes/donea/problem.hh | 67 --------
test/freeflow/navierstokes/kovasznay/main.cc | 17 +-
.../navierstokes/kovasznay/problem.hh | 71 --------
test/freeflow/navierstokes/sincos/main.cc | 79 +++------
test/freeflow/navierstokes/sincos/problem.hh | 2 -
13 files changed, 267 insertions(+), 417 deletions(-)
create mode 100644 test/freeflow/navierstokes/analyticalsolutionvectors.hh
diff --git a/test/freeflow/navierstokes/analyticalsolutionvectors.hh b/test/freeflow/navierstokes/analyticalsolutionvectors.hh
new file mode 100644
index 0000000000..fa41298e19
--- /dev/null
+++ b/test/freeflow/navierstokes/analyticalsolutionvectors.hh
@@ -0,0 +1,158 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup NavierStokesTests
+ * \copydoc Dumux::NavierStokesAnalyticalSolutionVectors
+ */
+#ifndef DUMUX_TEST_ANALYTICALSOLVECTORS_HH
+#define DUMUX_TEST_ANALYTICALSOLVECTORS_HH
+
+namespace Dumux {
+/*!
+ * \ingroup NavierStokesTests
+ * \brief Creates and provides the analytical solution in a form that can be written into the vtk output files
+ */
+template<class Scalar, class FaceSolutionVector, class CellCenterSolutionVector, class GridGeometry, class Problem, class Indices>
+class NavierStokesAnalyticalSolutionVectors
+{
+ 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>;
+
+public:
+ NavierStokesAnalyticalSolutionVectors(std::shared_ptr<const Problem> problem)
+ : problem_(problem)
+ { }
+
+ /*!
+ * \brief Creates the instationary analytical solution in a form that can be written into the vtk output files
+ */
+ void update(Scalar time)
+ {
+ auto instationaryAnalyticalPressureSolution = [&](const SubControlVolume& scv)
+ {
+ return problem_->analyticalSolution(scv.dofPosition(), time)[Indices::pressureIdx];
+ };
+
+ auto instationaryAnalyticalVelocitySolution = [&](const SubControlVolumeFace& scvf)
+ {
+ return problem_->analyticalSolution(scvf.center(), time)[Indices::velocity(scvf.directionIndex())];
+ };
+
+ auto instationaryAnalyticalVelocitySolutionAtScvCenter = [&](const SubControlVolume& scv)
+ {
+ return problem_->analyticalSolution(scv.center(), time);
+ };
+
+ createAnalyticalSolution_(instationaryAnalyticalPressureSolution, instationaryAnalyticalVelocitySolution, instationaryAnalyticalVelocitySolutionAtScvCenter);
+ }
+
+ /*!
+ * \brief Creates the stationary analytical solution in a form that can be written into the vtk output files
+ */
+ void update()
+ {
+ auto analyticalPressureSolution = [&](const SubControlVolume& scv)
+ {
+ return problem_->analyticalSolution(scv.dofPosition())[Indices::pressureIdx];
+ };
+
+ auto analyticalVelocitySolution = [&](const SubControlVolumeFace& scvf)
+ {
+ return problem_->analyticalSolution(scvf.center())[Indices::velocity(scvf.directionIndex())];
+ };
+
+ auto analyticalVelocitySolutionAtScvCenter = [&](const SubControlVolume& scv)
+ {
+ return problem_->analyticalSolution(scv.center());
+ };
+
+ createAnalyticalSolution_(analyticalPressureSolution, analyticalVelocitySolution, analyticalVelocitySolutionAtScvCenter);
+ }
+
+ /*!
+ * \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_;
+ }
+
+private:
+ /*!
+ * \brief Creates the analytical solution in a form that can be written into the vtk output files
+ */
+ template <class LambdaA, class LambdaB, class LambdaC>
+ void createAnalyticalSolution_(const LambdaA& analyticalPressureSolution,
+ const LambdaB& analyticalVelocitySolution,
+ const LambdaC& analyticalVelocitySolutionAtScvCenter)
+ {
+ analyticalPressure_.resize((problem_->gridGeometry()).numCellCenterDofs());
+ analyticalVelocity_.resize((problem_->gridGeometry()).numCellCenterDofs());
+ analyticalVelocityOnFace_.resize((problem_->gridGeometry()).numFaceDofs());
+
+ for (const auto& element : elements((problem_->gridGeometry()).gridView()))
+ {
+ auto fvGeometry = localView((problem_->gridGeometry()));
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ // velocities on faces
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ analyticalVelocityOnFace_[scvf.dofIndex()][scvf.directionIndex()] = analyticalVelocitySolution(scvf);
+ }
+
+ analyticalPressure_[scv.dofIndex()] = analyticalPressureSolution(scv);
+
+ for(int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
+ analyticalVelocity_[scv.dofIndex()][dirIdx] = analyticalVelocitySolutionAtScvCenter(scv)[Indices::velocity(dirIdx)];
+ }
+ }
+ }
+
+ std::shared_ptr<const Problem> problem_;
+
+ CellCenterSolutionVector analyticalPressure_;
+ std::vector<VelocityVector> analyticalVelocity_;
+ std::vector<VelocityVector> analyticalVelocityOnFace_;
+};
+} // end namespace Dumux
+
+#endif
diff --git a/test/freeflow/navierstokes/angeli/main.cc b/test/freeflow/navierstokes/angeli/main.cc
index 6233bb7531..d1769df4af 100644
--- a/test/freeflow/navierstokes/angeli/main.cc
+++ b/test/freeflow/navierstokes/angeli/main.cc
@@ -46,60 +46,7 @@
#include "properties.hh"
-/*!
-* \brief Creates analytical solution.
-* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
-* \param problem the problem for which to evaluate the analytical solution
-*/
-template<class Problem, class Scalar = double>
-auto createAnalyticalSolution(const Problem& problem)
-{
- const auto& gridGeometry = problem.gridGeometry();
- using GridView = typename std::decay_t<decltype(gridGeometry)>::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(gridGeometry.numCellCenterDofs());
- analyticalVelocity.resize(gridGeometry.numCellCenterDofs());
- analyticalVelocityOnFace.resize(gridGeometry.numFaceDofs());
-
- using Indices = typename Problem::Indices;
- auto fvGeometry = localView(gridGeometry);
- for (const auto& element : elements(gridGeometry.gridView()))
- {
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- auto ccDofIdx = scv.dofIndex();
- auto ccDofPosition = scv.dofPosition();
- auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
-
- // 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);
- 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);
-}
+#include "../analyticalsolutionvectors.hh"
int main(int argc, char** argv)
{
@@ -162,15 +109,22 @@ int main(int argc, char** argv)
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
+ // create analytical solution vectors
+ 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;
+
+ Dumux::NavierStokesAnalyticalSolutionVectors<Scalar, FaceSolutionVector, CellCenterSolutionVector, GridGeometry, Problem, Indices> analyticalSolVectors(problem);
+ analyticalSolVectors.update(tStart);
+
// initialize the vtk output module
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
-
- auto analyticalSolution = createAnalyticalSolution(*problem);
- vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
- vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
- vtkWriter.addFaceField(std::get<2>(analyticalSolution), "faceVelocityExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
+ vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
// the assembler with time loop for instationary problem
@@ -218,7 +172,7 @@ int main(int argc, char** argv)
}
// update the analytical solution
- analyticalSolution = createAnalyticalSolution(*problem);
+ analyticalSolVectors.update(timeLoop->time() + timeLoop->timeStepSize());
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
diff --git a/test/freeflow/navierstokes/angeli/problem.hh b/test/freeflow/navierstokes/angeli/problem.hh
index 2ae0f471cf..b8117d40c4 100644
--- a/test/freeflow/navierstokes/angeli/problem.hh
+++ b/test/freeflow/navierstokes/angeli/problem.hh
@@ -62,12 +62,10 @@ class AngeliTestProblem : 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 SubControlVolume = typename GridGeometry::SubControlVolume;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
public:
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
@@ -79,6 +77,7 @@ public:
rho_ = getParam<Scalar>("Component.LiquidDensity", 1.0);
useVelocityAveragingForDirichlet_ = getParam<bool>("Problem.UseVelocityAveragingForDirichlet", false);
useVelocityAveragingForInitial_ = getParam<bool>("Problem.UseVelocityAveragingForInitial", false);
+ tStart_ = getParam<Scalar>("TimeLoop.TStart", 0.0);
}
/*!
@@ -143,9 +142,9 @@ public:
PrimaryVariables dirichlet(const Element& element, const SubControlVolumeFace& scvf) const
{
if (useVelocityAveragingForDirichlet_)
- return averagedVelocity_(scvf);
+ return averagedVelocity_(scvf, time_);
else
- return analyticalSolution(scvf.center());
+ return analyticalSolution(scvf.center(), time_);
}
/*!
@@ -158,7 +157,7 @@ public:
PrimaryVariables dirichlet(const Element& element, const SubControlVolume& scv) const
{
PrimaryVariables priVars(0.0);
- priVars[Indices::pressureIdx] = analyticalSolution(scv.center())[Indices::pressureIdx];
+ priVars[Indices::pressureIdx] = analyticalSolution(scv.center(), time_)[Indices::pressureIdx];
return priVars;
}
@@ -166,11 +165,10 @@ public:
* \brief Returns the analytical solution of the problem at a given time and position.
* \param globalPos The global position
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar t) const
{
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
- const Scalar t = time_;
PrimaryVariables values;
@@ -181,6 +179,12 @@ public:
return values;
}
+ //TODO remove as soon as error calculation is able to deal with time-dependent analytical solution
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ {
+ return analyticalSolution(globalPos, time_);
+ }
+
// \}
/*!
@@ -194,7 +198,7 @@ public:
PrimaryVariables initial(const SubControlVolume& scv) const
{
PrimaryVariables priVars(0.0);
- priVars[Indices::pressureIdx] = analyticalSolution(scv.center())[Indices::pressureIdx];
+ priVars[Indices::pressureIdx] = analyticalSolution(scv.center(), tStart_)[Indices::pressureIdx];
return priVars;
}
@@ -211,9 +215,9 @@ public:
PrimaryVariables initial(const SubControlVolumeFace& scvf) const
{
if (useVelocityAveragingForInitial_)
- return averagedVelocity_(scvf);
+ return averagedVelocity_(scvf, tStart_);
else
- return analyticalSolution(scvf.center());
+ return analyticalSolution(scvf.center(), tStart_);
}
// \}
@@ -227,7 +231,7 @@ public:
}
private:
- PrimaryVariables averagedVelocity_(const SubControlVolumeFace& scvf) const
+ PrimaryVariables averagedVelocity_(const SubControlVolumeFace& scvf, Scalar t) const
{
PrimaryVariables priVars(0.0);
const auto geo = scvf.geometry();
@@ -236,7 +240,7 @@ private:
{
const auto w = qp.weight()*geo.integrationElement(qp.position());
const auto globalPos = geo.global(qp.position());
- const auto sol = analyticalSolution(globalPos);
+ const auto sol = analyticalSolution(globalPos, t);
priVars[Indices::velocityXIdx] += sol[Indices::velocityXIdx]*w;
priVars[Indices::velocityYIdx] += sol[Indices::velocityYIdx]*w;
}
@@ -249,6 +253,7 @@ private:
Scalar time_ = 0;
bool useVelocityAveragingForDirichlet_;
bool useVelocityAveragingForInitial_;
+ Scalar tStart_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokes/channel/1d/main.cc b/test/freeflow/navierstokes/channel/1d/main.cc
index f1aa15c446..7e8521fac5 100644
--- a/test/freeflow/navierstokes/channel/1d/main.cc
+++ b/test/freeflow/navierstokes/channel/1d/main.cc
@@ -44,6 +44,8 @@
#include "properties.hh"
+#include "../../analyticalsolutionvectors.hh"
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -91,13 +93,22 @@ int main(int argc, char** argv)
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
+ // create analytical solution vectors
+ 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;
+
+ Dumux::NavierStokesAnalyticalSolutionVectors<Scalar, FaceSolutionVector, CellCenterSolutionVector, GridGeometry, Problem, Indices> analyticalSolVectors(problem);
+ analyticalSolVectors.update();
+
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
- vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact");
- vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact");
- vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
+ vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
// the assembler with time loop for instationary problem
diff --git a/test/freeflow/navierstokes/channel/1d/problem.hh b/test/freeflow/navierstokes/channel/1d/problem.hh
index 19dc7edbf1..66b9823b1e 100644
--- a/test/freeflow/navierstokes/channel/1d/problem.hh
+++ b/test/freeflow/navierstokes/channel/1d/problem.hh
@@ -73,7 +73,6 @@ public:
printL2Error_ = getParam<bool>("Problem.PrintL2Error");
density_ = getParam<Scalar>("Component.LiquidDensity");
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity");
- createAnalyticalSolution_();
}
/*!
@@ -285,76 +284,10 @@ public:
return analyticalSolution(globalPos);
}
- /*!
- * \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_;
- }
-
private:
-
- /*!
- * \brief Adds additional VTK output data to the VTKWriter.
- *
- * Function is called by the output module on every write.
- */
- void createAnalyticalSolution_()
- {
- analyticalPressure_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocity_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocityOnFace_.resize(this->gridGeometry().numFaceDofs());
- auto fvGeometry = localView(this->gridGeometry());
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto ccDofIdx = scv.dofIndex();
- const auto ccDofPosition = scv.dofPosition();
- const auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition);
-
- // 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);
- 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)];
- }
- }
- }
-
bool printL2Error_;
Scalar density_;
Scalar kinematicViscosity_;
- std::vector<Scalar> analyticalPressure_;
- std::vector<GlobalPosition> analyticalVelocity_;
- std::vector<GlobalPosition> analyticalVelocityOnFace_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokes/channel/2d/main.cc b/test/freeflow/navierstokes/channel/2d/main.cc
index 2c180790aa..f0e9352d90 100644
--- a/test/freeflow/navierstokes/channel/2d/main.cc
+++ b/test/freeflow/navierstokes/channel/2d/main.cc
@@ -45,6 +45,8 @@
#include "properties.hh"
+#include "../../analyticalsolutionvectors.hh"
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -108,6 +110,12 @@ int main(int argc, char** argv)
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
+ // create analytical solution vectors
+ 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;
+
// initialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
@@ -117,9 +125,12 @@ int main(int argc, char** argv)
if (problem->hasAnalyticalSolution())
{
- vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact");
- vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact");
- vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+ Dumux::NavierStokesAnalyticalSolutionVectors<Scalar, FaceSolutionVector, CellCenterSolutionVector, GridGeometry, Problem, Indices> analyticalSolVectors(problem);
+ analyticalSolVectors.update();
+
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
+ vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
}
vtkWriter.write(restartTime);
diff --git a/test/freeflow/navierstokes/channel/2d/problem.hh b/test/freeflow/navierstokes/channel/2d/problem.hh
index 96122cf361..a39a7b46eb 100644
--- a/test/freeflow/navierstokes/channel/2d/problem.hh
+++ b/test/freeflow/navierstokes/channel/2d/problem.hh
@@ -66,8 +66,6 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- static constexpr auto dimWorld = GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld;
- using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
@@ -102,8 +100,6 @@ public:
const bool isStationary = getParam<bool>("Problem.IsStationary", false);
hasAnalyticalSolution_ = isStationary && useVelocityProfile_ && (outletCondition_ != OutletCondition::doNothing);
- if (hasAnalyticalSolution_)
- createAnalyticalSolution_();
}
/*!
@@ -338,73 +334,12 @@ public:
return timeLoop_->time();
}
-/*!
- * \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_;
- }
-
bool hasAnalyticalSolution()
{
return hasAnalyticalSolution_;
}
private:
-
- /*!
- * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
- */
- void createAnalyticalSolution_()
- {
- analyticalPressure_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocity_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocityOnFace_.resize(this->gridGeometry().numFaceDofs());
- auto fvGeometry = localView(this->gridGeometry());
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- auto ccDofIdx = scv.dofIndex();
- auto ccDofPosition = scv.dofPosition();
- auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition);
-
- // 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);
- 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)];
- }
- }
- }
-
bool isInlet_(const GlobalPosition& globalPos) const
{
return globalPos[0] < eps_;
@@ -424,10 +359,6 @@ private:
bool useVelocityProfile_;
bool hasAnalyticalSolution_;
TimeLoopPtr timeLoop_;
-
- std::vector<Scalar> analyticalPressure_;
- std::vector<VelocityVector> analyticalVelocity_;
- std::vector<VelocityVector> analyticalVelocityOnFace_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokes/donea/main.cc b/test/freeflow/navierstokes/donea/main.cc
index 177bb4f434..73bddb0732 100644
--- a/test/freeflow/navierstokes/donea/main.cc
+++ b/test/freeflow/navierstokes/donea/main.cc
@@ -45,6 +45,8 @@
#include "properties.hh"
+#include "../analyticalsolutionvectors.hh"
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -93,13 +95,22 @@ int main(int argc, char** argv)
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
+ // create analytical solution vectors
+ 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;
+
+ Dumux::NavierStokesAnalyticalSolutionVectors<Scalar, FaceSolutionVector, CellCenterSolutionVector, GridGeometry, Problem, Indices> analyticalSolVectors(problem);
+ analyticalSolVectors.update();
+
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
- vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact");
- vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact");
- vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
+ vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
// use the staggered FV assembler
diff --git a/test/freeflow/navierstokes/donea/problem.hh b/test/freeflow/navierstokes/donea/problem.hh
index e61ac2b385..7d8fb0ce0e 100644
--- a/test/freeflow/navierstokes/donea/problem.hh
+++ b/test/freeflow/navierstokes/donea/problem.hh
@@ -57,17 +57,14 @@ 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)
: ParentType(gridGeometry)
{
printL2Error_ = getParam<bool>("Problem.PrintL2Error");
- createAnalyticalSolution_();
mu_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
}
@@ -210,30 +207,6 @@ public:
return values;
}
- /*!
- * \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_;
- }
-
private:
Scalar f1_(Scalar x) const
{ return x*(1.0-x); /*x - x^2*/ }
@@ -283,47 +256,7 @@ private:
Scalar dxxV_ (Scalar x, Scalar y) const
{ return -f2_(y)*dddf2_(x); }
- /*!
- * \brief Adds additional VTK output data to the VTKWriter. Function is called by the output module on every write.
- */
- void createAnalyticalSolution_()
- {
- analyticalPressure_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocity_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocityOnFace_.resize(this->gridGeometry().numFaceDofs());
-
- auto fvGeometry = localView(this->gridGeometry());
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto ccDofIdx = scv.dofIndex();
- const auto ccDofPosition = scv.dofPosition();
- const auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition);
-
- // 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);
- 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)];
- }
- }
- }
-
bool printL2Error_;
- std::vector<Scalar> analyticalPressure_;
- std::vector<VelocityVector> analyticalVelocity_;
- std::vector<VelocityVector> analyticalVelocityOnFace_;
Scalar mu_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokes/kovasznay/main.cc b/test/freeflow/navierstokes/kovasznay/main.cc
index c9af40ea99..0b3f316cb1 100644
--- a/test/freeflow/navierstokes/kovasznay/main.cc
+++ b/test/freeflow/navierstokes/kovasznay/main.cc
@@ -46,6 +46,8 @@
#include "properties.hh"
+#include "../analyticalsolutionvectors.hh"
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -121,13 +123,22 @@ int main(int argc, char** argv)
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
+ // create analytical solution vectors
+ 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;
+
+ Dumux::NavierStokesAnalyticalSolutionVectors<Scalar, FaceSolutionVector, CellCenterSolutionVector, GridGeometry, Problem, Indices> analyticalSolVectors(problem);
+ analyticalSolVectors.update();
+
// intialize the vtk output module
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
- vtkWriter.addField(problem->getAnalyticalPressureSolution(), "pressureExact");
- vtkWriter.addField(problem->getAnalyticalVelocitySolution(), "velocityExact");
- vtkWriter.addFaceField(problem->getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
+ vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
vtkWriter.write(0.0);
// the assembler with time loop for instationary problem
diff --git a/test/freeflow/navierstokes/kovasznay/problem.hh b/test/freeflow/navierstokes/kovasznay/problem.hh
index 748d6b4654..a8bc9b170c 100644
--- a/test/freeflow/navierstokes/kovasznay/problem.hh
+++ b/test/freeflow/navierstokes/kovasznay/problem.hh
@@ -60,10 +60,8 @@ class KovasznayTestProblem : 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>;
static constexpr auto upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>();
@@ -78,8 +76,6 @@ public:
Scalar reynoldsNumber = 1.0 / kinematicViscosity_;
lambda_ = 0.5 * reynoldsNumber
- std::sqrt(reynoldsNumber * reynoldsNumber * 0.25 + 4.0 * M_PI * M_PI);
-
- createAnalyticalSolution_();
}
/*!
@@ -214,80 +210,13 @@ public:
return values;
}
- /*!
- * \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_;
- }
-
private:
-
- /*!
- * \brief Adds additional VTK output data to the VTKWriter.
- *
- * Function is called by the output module on every write.
- */
- void createAnalyticalSolution_()
- {
- analyticalPressure_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocity_.resize(this->gridGeometry().numCellCenterDofs());
- analyticalVelocityOnFace_.resize(this->gridGeometry().numFaceDofs());
-
- auto fvGeometry = localView(this->gridGeometry());
- for (const auto& element : elements(this->gridGeometry().gridView()))
- {
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto ccDofIdx = scv.dofIndex();
- const auto ccDofPosition = scv.dofPosition();
- const auto analyticalSolutionAtCc = analyticalSolution(ccDofPosition);
-
- // 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);
- 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)];
- }
- }
- }
-
static constexpr Scalar eps_=1e-6;
Scalar rho_;
Scalar kinematicViscosity_;
Scalar lambda_;
bool printL2Error_;
- std::vector<Scalar> analyticalPressure_;
- std::vector<VelocityVector> analyticalVelocity_;
- std::vector<VelocityVector> analyticalVelocityOnFace_;
};
} // end namespace Dumux
diff --git a/test/freeflow/navierstokes/sincos/main.cc b/test/freeflow/navierstokes/sincos/main.cc
index 90b37841fc..c965dc6c25 100644
--- a/test/freeflow/navierstokes/sincos/main.cc
+++ b/test/freeflow/navierstokes/sincos/main.cc
@@ -48,62 +48,14 @@
#include "properties.hh"
+#include "../analyticalsolutionvectors.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& gridGeometry = problem.gridGeometry();
- using GridView = typename std::decay_t<decltype(gridGeometry)>::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(gridGeometry.numCellCenterDofs());
- analyticalVelocity.resize(gridGeometry.numCellCenterDofs());
- analyticalVelocityOnFace.resize(gridGeometry.numFaceDofs());
-
- using Indices = typename Problem::Indices;
- auto fvGeometry = localView(gridGeometry);
- for (const auto& element : elements(gridGeometry.gridView()))
- {
- fvGeometry.bindElement(element);
- for (auto&& scv : scvs(fvGeometry))
- {
- const auto ccDofIdx = scv.dofIndex();
- const auto ccDofPosition = scv.dofPosition();
- const 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);
-}
-
template<class Problem>
auto createSource(const Problem& problem)
{
@@ -237,13 +189,26 @@ int main(int argc, char** argv)
vtkWriter.addField(sourceX, "sourceX");
vtkWriter.addField(sourceY, "sourceY");
- auto analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
- vtkWriter.addField(std::get<0>(analyticalSolution), "pressureExact");
- vtkWriter.addField(std::get<1>(analyticalSolution), "velocityExact");
- vtkWriter.addFaceField(std::get<2>(analyticalSolution), "faceVelocityExact");
+ const bool isStationary = getParam<bool>("Problem.IsStationary");
+
+ // create analytical solution vectors
+ 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;
+
+ Dumux::NavierStokesAnalyticalSolutionVectors<Scalar, FaceSolutionVector, CellCenterSolutionVector, GridGeometry, Problem, Indices> analyticalSolVectors(problem);
+ if (isStationary)
+ analyticalSolVectors.update();
+ else
+ analyticalSolVectors.update(0.0);
+
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
+ vtkWriter.addField(analyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
+ vtkWriter.addFaceField(analyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+
vtkWriter.write(0.0);
- const bool isStationary = getParam<bool>("Problem.IsStationary");
// the assembler with time loop for instationary problem
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
@@ -272,7 +237,7 @@ int main(int argc, char** argv)
}
// write vtk output
- analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
+ analyticalSolVectors.update();
vtkWriter.write(1.0);
timer.stop();
@@ -302,7 +267,7 @@ int main(int argc, char** argv)
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
problem->updateTime(timeLoop->time());
- analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
+ analyticalSolVectors.update(timeLoop->time());
// write vtk output
vtkWriter.write(timeLoop->time());
diff --git a/test/freeflow/navierstokes/sincos/problem.hh b/test/freeflow/navierstokes/sincos/problem.hh
index d197933e7d..3c43fbea8f 100644
--- a/test/freeflow/navierstokes/sincos/problem.hh
+++ b/test/freeflow/navierstokes/sincos/problem.hh
@@ -58,10 +58,8 @@ class SincosTestProblem : public NavierStokesProblem<TypeTag>
using SubControlVolume = typename GridGeometry::SubControlVolume;
using FVElementGeometry = typename GridGeometry::LocalView;
- 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:
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
--
GitLab