diff --git a/test/freeflow/navierstokes/analyticalsolutionvectors.hh b/test/freeflow/navierstokes/analyticalsolutionvectors.hh
new file mode 100644
index 0000000000000000000000000000000000000000..eb3e2bb9850a4dc26cf3953560664933538b5fa5
--- /dev/null
+++ b/test/freeflow/navierstokes/analyticalsolutionvectors.hh
@@ -0,0 +1,116 @@
+// -*- 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 Problem, class Scalar = double>
+class NavierStokesAnalyticalSolutionVectors
+{
+ using GridGeometry = std::decay_t<decltype(std::declval<Problem>().gridGeometry())>;
+ static constexpr int dimWorld = GridGeometry::GridView::dimensionworld;
+ using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+ using Indices = typename Problem::Indices;
+
+public:
+ NavierStokesAnalyticalSolutionVectors(std::shared_ptr<const Problem> problem, Scalar tInitial = 0.0)
+ : problem_(problem)
+ { update(tInitial); }
+
+ /*!
+ * \brief Creates the analytical solution in a form that can be written into the vtk output files
+ */
+ void update(Scalar time = 0.0)
+ {
+ analyticalPressure_.resize(problem_->gridGeometry().numCellCenterDofs());
+ analyticalVelocity_.resize(problem_->gridGeometry().numCellCenterDofs());
+ analyticalVelocityOnFace_.resize(problem_->gridGeometry().numFaceDofs());
+
+ auto fvGeometry = localView(problem_->gridGeometry());
+
+ for (const auto& element : elements(problem_->gridGeometry().gridView()))
+ {
+ fvGeometry.bindElement(element);
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ // velocities on faces
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ analyticalVelocityOnFace_[scvf.dofIndex()][scvf.directionIndex()] = problem_->analyticalSolution(scvf.center(), time)[Indices::velocity(scvf.directionIndex())];
+ }
+
+ analyticalPressure_[scv.dofIndex()] = problem_->analyticalSolution(scv.dofPosition(), time)[Indices::pressureIdx];
+
+ for (int dirIdx = 0; dirIdx < dimWorld; ++dirIdx)
+ analyticalVelocity_[scv.dofIndex()][dirIdx] = problem_->analyticalSolution(scv.center(), time)[Indices::velocity(dirIdx)];
+ }
+ }
+ }
+
+ /*!
+ * \brief Returns the analytical solution for the pressure
+ */
+ const std::vector<Scalar>& getAnalyticalPressureSolution() const
+ {
+ return analyticalPressure_;
+ }
+
+ /*!
+ * \brief Returns the analytical solution for the velocity
+ */
+ const std::vector<VelocityVector>& getAnalyticalVelocitySolution() const
+ {
+ return analyticalVelocity_;
+ }
+
+ /*!
+ * \brief Returns the analytical solution for the velocity at the faces
+ */
+ const std::vector<VelocityVector>& getAnalyticalVelocitySolutionOnFace() const
+ {
+ return analyticalVelocityOnFace_;
+ }
+
+private:
+ std::shared_ptr<const Problem> problem_;
+
+ std::vector<Scalar> analyticalPressure_;
+ std::vector<VelocityVector> analyticalVelocity_;
+ std::vector<VelocityVector> analyticalVelocityOnFace_;
+};
+
+template<class Problem>
+NavierStokesAnalyticalSolutionVectors(std::shared_ptr<Problem> p)
+-> NavierStokesAnalyticalSolutionVectors<Problem>;
+
+template<class Problem, class Scalar>
+NavierStokesAnalyticalSolutionVectors(std::shared_ptr<Problem> p, Scalar t)
+-> NavierStokesAnalyticalSolutionVectors<Problem, Scalar>;
+} // end namespace Dumux
+
+#endif
diff --git a/test/freeflow/navierstokes/angeli/main.cc b/test/freeflow/navierstokes/angeli/main.cc
index 6233bb7531e88fddfce82583353457aa27428856..1ba0fb16dac369e679578bca348af6412625da5b 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 <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
int main(int argc, char** argv)
{
@@ -166,11 +113,10 @@ int main(int argc, char** argv)
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");
+ NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem, tStart);
+ 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
@@ -205,7 +151,7 @@ int main(int argc, char** argv)
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using L2Error = NavierStokesTestL2Error<Scalar, ModelTraits, PrimaryVariables>;
- const auto [l2errorAbs, l2errorRel] = L2Error::calculateL2Error(*problem, x);
+ const auto [l2errorAbs, l2errorRel] = L2Error::calculateL2Error(*problem, x, timeLoop->time() + timeLoop->timeStepSize());
const int numCellCenterDofs = gridGeometry->numCellCenterDofs();
const int numFaceDofs = gridGeometry->numFaceDofs();
std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
@@ -218,7 +164,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 2ae0f471cf90e2dd98db9792566a5cf992b72560..adcad9b279ebd6ddcb701f9ad15a6e4edd6c59cb 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;
@@ -143,9 +141,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,25 +156,25 @@ 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;
}
/*!
* \brief Returns the analytical solution of the problem at a given time and position.
* \param globalPos The global position
+ * \param time The current simulation time
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time) const
{
const Scalar x = globalPos[0];
const Scalar y = globalPos[1];
- const Scalar t = time_;
PrimaryVariables values;
- values[Indices::pressureIdx] = - 0.25 * std::exp(-10.0 * kinematicViscosity_ * M_PI * M_PI * t) * M_PI * M_PI * (4.0 * std::cos(2.0 * M_PI * x) + std::cos(4.0 * M_PI * y))*rho_;
- values[Indices::velocityXIdx] = - 2.0 * M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * t) * std::cos(M_PI * x) * std::sin(2.0 * M_PI * y);
- values[Indices::velocityYIdx] = M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * t) * std::sin(M_PI * x) * std::cos(2.0 * M_PI * y);
+ values[Indices::pressureIdx] = - 0.25 * std::exp(-10.0 * kinematicViscosity_ * M_PI * M_PI * time) * M_PI * M_PI * (4.0 * std::cos(2.0 * M_PI * x) + std::cos(4.0 * M_PI * y))*rho_;
+ values[Indices::velocityXIdx] = - 2.0 * M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * std::cos(M_PI * x) * std::sin(2.0 * M_PI * y);
+ values[Indices::velocityYIdx] = M_PI * std::exp(- 5.0 * kinematicViscosity_ * M_PI * M_PI * time) * std::sin(M_PI * x) * std::cos(2.0 * M_PI * y);
return values;
}
@@ -194,7 +192,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(), time_)[Indices::pressureIdx];
return priVars;
}
@@ -211,9 +209,9 @@ public:
PrimaryVariables initial(const SubControlVolumeFace& scvf) const
{
if (useVelocityAveragingForInitial_)
- return averagedVelocity_(scvf);
+ return averagedVelocity_(scvf, time_);
else
- return analyticalSolution(scvf.center());
+ return analyticalSolution(scvf.center(), time_);
}
// \}
@@ -227,7 +225,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 +234,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;
}
diff --git a/test/freeflow/navierstokes/channel/1d/main.cc b/test/freeflow/navierstokes/channel/1d/main.cc
index f1aa15c4465948b231fc8c564180bd9dd1a99977..fed25692a7e7d9be7eb9997e828264e810d0d21f 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 <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -95,9 +97,12 @@ int main(int argc, char** argv)
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");
+
+ NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem);
+ 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 19dc7edbf1fc739aa772e592435ccc11722da073..7f345c52b9f6def6bc1e65a2b1a13b59716c5a77 100644
--- a/test/freeflow/navierstokes/channel/1d/problem.hh
+++ b/test/freeflow/navierstokes/channel/1d/problem.hh
@@ -53,7 +53,6 @@ class NavierStokesAnalyticProblem : public NavierStokesProblem<TypeTag>
using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
@@ -67,13 +66,14 @@ class NavierStokesAnalyticProblem : public NavierStokesProblem<TypeTag>
using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
public:
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
NavierStokesAnalyticProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
printL2Error_ = getParam<bool>("Problem.PrintL2Error");
density_ = getParam<Scalar>("Component.LiquidDensity");
kinematicViscosity_ = getParam<Scalar>("Component.LiquidKinematicViscosity");
- createAnalyticalSolution_();
}
/*!
@@ -208,8 +208,9 @@ public:
* \brief Returns the analytical solution of the problem at a given position
*
* \param globalPos The global position
+ * \param time A parameter for consistent signatures. It is ignored here as this is a stationary test.
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
PrimaryVariables values;
values[Indices::pressureIdx] = p(globalPos);
@@ -285,76 +286,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 2c180790aa28daa3df6e7cca86cf110db7b0d98f..98d268c59f5fa0355ccd9873e5e5b447d5e02538 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 <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -115,11 +117,12 @@ int main(int argc, char** argv)
const bool isStationary = getParam<bool>("Problem.IsStationary", false);
+ NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem);
if (problem->hasAnalyticalSolution())
{
- 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(restartTime);
diff --git a/test/freeflow/navierstokes/channel/2d/problem.hh b/test/freeflow/navierstokes/channel/2d/problem.hh
index 96122cf3616465c8ef77827e3ef44570515ba6a9..89389521fc6164b808d9b0ba42e9f4a129d3f6d2 100644
--- a/test/freeflow/navierstokes/channel/2d/problem.hh
+++ b/test/freeflow/navierstokes/channel/2d/problem.hh
@@ -57,7 +57,6 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
@@ -66,8 +65,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>>;
@@ -78,6 +75,8 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
};
public:
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
ChannelTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
@@ -102,8 +101,6 @@ public:
const bool isStationary = getParam<bool>("Problem.IsStationary", false);
hasAnalyticalSolution_ = isStationary && useVelocityProfile_ && (outletCondition_ != OutletCondition::doNothing);
- if (hasAnalyticalSolution_)
- createAnalyticalSolution_();
}
/*!
@@ -273,8 +270,9 @@ public:
* \brief Return the analytical solution of the problem at a given position
*
* \param globalPos The global position
+ * \param time A parameter for consistent signatures. It is ignored here as this analytical solution is for a stationary version of the test only.
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
Scalar x = globalPos[0];
Scalar y = globalPos[1];
@@ -338,73 +336,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 +361,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/channel/pipe/problem.hh b/test/freeflow/navierstokes/channel/pipe/problem.hh
index 625a5389c671f42290440792de464b55531fccb1..833363066c14e12f253cddfe7d77e23d722eb3ad 100644
--- a/test/freeflow/navierstokes/channel/pipe/problem.hh
+++ b/test/freeflow/navierstokes/channel/pipe/problem.hh
@@ -117,7 +117,7 @@ public:
return values;
}
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
PrimaryVariables values(0.0);
diff --git a/test/freeflow/navierstokes/donea/main.cc b/test/freeflow/navierstokes/donea/main.cc
index 177bb4f4349c2c8549076732e2c8145479bae0f8..a777e1babc58e03ef7f29637d0e4eaee194f756b 100644
--- a/test/freeflow/navierstokes/donea/main.cc
+++ b/test/freeflow/navierstokes/donea/main.cc
@@ -45,6 +45,8 @@
#include "properties.hh"
+#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -97,9 +99,12 @@ int main(int argc, char** argv)
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");
+
+ NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem);
+ 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 e61ac2b385cef7bfb38f59c6498f7810879194b4..73631c676b7b431dd54a0a9b4948ae84ad2663f9 100644
--- a/test/freeflow/navierstokes/donea/problem.hh
+++ b/test/freeflow/navierstokes/donea/problem.hh
@@ -50,24 +50,22 @@ class DoneaTestProblem : public NavierStokesProblem<TypeTag>
using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
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:
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
DoneaTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
printL2Error_ = getParam<bool>("Problem.PrintL2Error");
- createAnalyticalSolution_();
mu_ = getParam<Scalar>("Component.LiquidKinematicViscosity", 1.0);
}
@@ -174,8 +172,9 @@ public:
* \brief Return the analytical solution of the problem at a given position
*
* \param globalPos The global position
+ * \param time A parameter for consistent signatures. It is ignored here as this is a stationary test
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
Scalar x = globalPos[0];
Scalar y = globalPos[1];
@@ -210,30 +209,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 +258,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 c9af40ea99d41d8dd953c7f255e6a75c3e753742..c42b7457bc7b5d3ef5614f0bf5409236cb90b56c 100644
--- a/test/freeflow/navierstokes/kovasznay/main.cc
+++ b/test/freeflow/navierstokes/kovasznay/main.cc
@@ -46,6 +46,8 @@
#include "properties.hh"
+#include <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
+
int main(int argc, char** argv)
{
using namespace Dumux;
@@ -125,9 +127,12 @@ int main(int argc, char** argv)
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");
+
+ NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem);
+ 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 748d6b4654ef2140cac386b946ffa47035fcb78a..a6bdb73ec2a93cc9cda534ffd4bff04b4eb40e58 100644
--- a/test/freeflow/navierstokes/kovasznay/problem.hh
+++ b/test/freeflow/navierstokes/kovasznay/problem.hh
@@ -53,21 +53,20 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag>
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
- using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
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>();
public:
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
KovasznayTestProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
@@ -78,8 +77,6 @@ public:
Scalar reynoldsNumber = 1.0 / kinematicViscosity_;
lambda_ = 0.5 * reynoldsNumber
- std::sqrt(reynoldsNumber * reynoldsNumber * 0.25 + 4.0 * M_PI * M_PI);
-
- createAnalyticalSolution_();
}
/*!
@@ -178,8 +175,9 @@ public:
* \brief Returns the analytical solution of the problem at a given position.
*
* \param globalPos The global position
+ * \param time A parameter for consistent signatures. It is ignored here as this is a stationary test
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
Scalar x = globalPos[0];
Scalar y = globalPos[1];
@@ -214,80 +212,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/l2error.hh b/test/freeflow/navierstokes/l2error.hh
index 9a5f866ebd7c73887f82bfb9019ce0410b26fb97..631df5a88f86f23ad5300e0df90290ad3c8a21e1 100644
--- a/test/freeflow/navierstokes/l2error.hh
+++ b/test/freeflow/navierstokes/l2error.hh
@@ -48,7 +48,7 @@ public:
* \param curSol Vector containing the current solution
*/
template<class Problem, class SolutionVector>
- static auto calculateL2Error(const Problem& problem, const SolutionVector& curSol)
+ static auto calculateL2Error(const Problem& problem, const SolutionVector& curSol, Scalar time = 0.0)
{
using GridGeometry = std::decay_t<decltype(problem.gridGeometry())>;
using Extrusion = Extrusion_t<GridGeometry>;
@@ -72,7 +72,7 @@ public:
// treat cell-center dofs
const auto dofIdxCellCenter = scv.dofIndex();
const auto& posCellCenter = scv.dofPosition();
- const auto analyticalSolutionCellCenter = problem.analyticalSolution(posCellCenter)[Indices::pressureIdx];
+ const auto analyticalSolutionCellCenter = problem.analyticalSolution(posCellCenter, time)[Indices::pressureIdx];
const auto numericalSolutionCellCenter = curSol[GridGeometry::cellCenterIdx()][dofIdxCellCenter][Indices::pressureIdx - ModelTraits::dim()];
sumError[Indices::pressureIdx] += squaredDiff_(analyticalSolutionCellCenter, numericalSolutionCellCenter) * Extrusion::volume(scv);
sumReference[Indices::pressureIdx] += analyticalSolutionCellCenter * analyticalSolutionCellCenter * Extrusion::volume(scv);
@@ -83,7 +83,7 @@ public:
{
const int dofIdxFace = scvf.dofIndex();
const int dirIdx = scvf.directionIndex();
- const auto analyticalSolutionFace = problem.analyticalSolution(scvf.center())[Indices::velocity(dirIdx)];
+ const auto analyticalSolutionFace = problem.analyticalSolution(scvf.center(), time)[Indices::velocity(dirIdx)];
const auto numericalSolutionFace = curSol[GridGeometry::faceIdx()][dofIdxFace][0];
directionIndex[dofIdxFace] = dirIdx;
errorVelocity[dofIdxFace] = squaredDiff_(analyticalSolutionFace, numericalSolutionFace);
diff --git a/test/freeflow/navierstokes/sincos/main.cc b/test/freeflow/navierstokes/sincos/main.cc
index 90b37841fcdc4de9f7243daa2785d91ab0e425d0..d358514427c4b089f2a18b6c731ab8330aaa9546 100644
--- a/test/freeflow/navierstokes/sincos/main.cc
+++ b/test/freeflow/navierstokes/sincos/main.cc
@@ -48,62 +48,14 @@
#include "properties.hh"
+#include <test/freeflow/navierstokes/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)
{
@@ -137,18 +89,17 @@ auto createSource(const Problem& problem)
return source;
}
-template<class Problem, class SolutionVector, class GridGeometry>
-void printL2Error(const Problem& problem, const SolutionVector& x, const GridGeometry& gridGeometry)
+template<class Problem, class SolutionVector, class GridGeometry, class Scalar = double>
+void printL2Error(const Problem& problem, const SolutionVector& x, const GridGeometry& gridGeometry, Scalar time = 0.0)
{
using namespace Dumux;
- using Scalar = double;
using TypeTag = Properties::TTag::SincosTest;
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 auto l2error = L2Error::calculateL2Error(*problem, x, time);
const int numCellCenterDofs = gridGeometry->numCellCenterDofs();
const int numFaceDofs = gridGeometry->numFaceDofs();
std::cout << std::setprecision(8) << "** L2 error (abs/rel) for "
@@ -237,10 +188,11 @@ 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");
+ NavierStokesAnalyticalSolutionVectors analyticalSolVectors(problem, 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");
@@ -272,7 +224,7 @@ int main(int argc, char** argv)
}
// write vtk output
- analyticalSolution = createAnalyticalSolution(timeLoop->time(), *problem);
+ analyticalSolVectors.update();
vtkWriter.write(1.0);
timer.stop();
@@ -296,13 +248,13 @@ int main(int argc, char** argv)
if (shouldPrintL2Error)
{
- printL2Error(problem, x, gridGeometry);
+ printL2Error(problem, x, gridGeometry, timeLoop->time()+timeLoop->timeStepSize());
}
// 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 d197933e7d236dcfe0edbac131906bdd099b4133..f78eb27e1732f44539f6806197ff73a3ef56baa7 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>;
@@ -187,16 +185,6 @@ public:
return values;
}
- /*!
- * \brief Returns the analytical solution of the problem at a given position.
- *
- * \param globalPos The global position
- */
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
- {
- return analyticalSolution(globalPos, time_+timeStepSize_);
- }
-
// \}
/*!
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/main.cc b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/main.cc
index d39d370e52ecb68522525cce4b72df06773d0f5a..89490287feb3036685288b449f1f772c192e42c4 100644
--- a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/main.cc
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/main.cc
@@ -47,60 +47,7 @@
#include "testcase.hh"
#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 Scalar, class Problem>
-auto createFreeFlowAnalyticalSolution(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<Scalar> 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] = 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 <test/freeflow/navierstokes/analyticalsolutionvectors.hh>
/*!
* \brief Creates analytical solution.
@@ -303,10 +250,12 @@ int main(int argc, char** argv)
using Scalar = typename Traits::Scalar;
StaggeredVtkOutputModule<FreeFlowGridVariables, decltype(freeFlowSol)> freeFlowVtkWriter(*freeFlowGridVariables, freeFlowSol, freeFlowProblem->name());
GetPropType<FreeFlowTypeTag, Properties::IOFields>::initOutputModule(freeFlowVtkWriter);
- const auto freeFlowAnalyticalSolution = createFreeFlowAnalyticalSolution<Scalar>(*freeFlowProblem);
- freeFlowVtkWriter.addField(std::get<0>(freeFlowAnalyticalSolution), "pressureExact");
- freeFlowVtkWriter.addField(std::get<1>(freeFlowAnalyticalSolution), "velocityExact");
- freeFlowVtkWriter.addFaceField(std::get<2>(freeFlowAnalyticalSolution), "faceVelocityExact");
+
+ NavierStokesAnalyticalSolutionVectors freeFlowAnalyticalSolVectors(freeFlowProblem);
+ freeFlowVtkWriter.addField(freeFlowAnalyticalSolVectors.getAnalyticalPressureSolution(), "pressureExact");
+ freeFlowVtkWriter.addField(freeFlowAnalyticalSolVectors.getAnalyticalVelocitySolution(), "velocityExact");
+ freeFlowVtkWriter.addFaceField(freeFlowAnalyticalSolVectors.getAnalyticalVelocitySolutionOnFace(), "faceVelocityExact");
+
freeFlowVtkWriter.write(0.0);
VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyProblem->name());
diff --git a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_stokes.hh b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_stokes.hh
index 42e954634e6c8a3e392e80efc86f1af8ce4c97ab..1c71a5355641ca9f069fb4660a212ec93df94072 100644
--- a/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_stokes.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p_1p/convergencetest/problem_stokes.hh
@@ -233,8 +233,9 @@ public:
* \brief Returns the analytical solution of the problem at a given position.
*
* \param globalPos The global position
+ * \param time A parameter for consistent signatures. It is ignored here as this is a stationary test
*/
- PrimaryVariables analyticalSolution(const GlobalPosition& globalPos) const
+ PrimaryVariables analyticalSolution(const GlobalPosition& globalPos, Scalar time = 0.0) const
{
switch (testCase_)
{