26#ifndef DUMUX_DONEA_TEST_PROBLEM_HH
27#define DUMUX_DONEA_TEST_PROBLEM_HH
29#include <dune/common/fmatrix.hh>
30#include <dune/grid/yaspgrid.hh>
38#include "../../l2error.hh"
42template <
class TypeTag>
43class NavierStokesAnalyticProblem;
52template<
class TypeTag>
60template<
class TypeTag>
61struct Grid<TypeTag, TTag::NavierStokesAnalytic> {
using type = Dune::YaspGrid<1>; };
64template<
class TypeTag>
67template<
class TypeTag>
70template<
class TypeTag>
72template<
class TypeTag>
75template<
class TypeTag>
76struct NormalizePressure<TypeTag, TTag::NavierStokesAnalytic> {
static constexpr bool value =
false; };
87template <
class TypeTag>
102 using Element =
typename GridGeometry::GridView::template Codim<0>::Entity;
103 using GlobalPosition =
typename Element::Geometry::GlobalCoordinate;
104 using DimVector = GlobalPosition;
105 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
109 : ParentType(gridGeometry), eps_(1e-6)
111 printL2Error_ = getParam<bool>(
"Problem.PrintL2Error");
112 density_ = getParam<Scalar>(
"Component.LiquidDensity");
113 kinematicViscosity_ = getParam<Scalar>(
"Component.LiquidKinematicViscosity");
114 createAnalyticalSolution_();
132 const auto l2error = L2Error::calculateL2Error(*
this, curSol);
133 const int numCellCenterDofs = this->gridGeometry().numCellCenterDofs();
134 const int numFaceDofs = this->gridGeometry().numFaceDofs();
135 std::cout << std::setprecision(8) <<
"** L2 error (abs/rel) for "
136 << std::setw(6) << numCellCenterDofs <<
" cc dofs and " << numFaceDofs <<
" face dofs (total: " << numCellCenterDofs + numFaceDofs <<
"): "
138 <<
"L2(p) = " << l2error.first[Indices::pressureIdx] <<
" / " << l2error.second[Indices::pressureIdx]
139 <<
" , L2(vx) = " << l2error.first[Indices::velocityXIdx] <<
" / " << l2error.second[Indices::velocityXIdx]
159 NumEqVector source(0.0);
162 for (
unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
164 source[Indices::conti0EqIdx] +=
dvdx(globalPos)[dimIdx][dimIdx];
166 source[Indices::conti0EqIdx] *= density_;
169 for (
unsigned int velIdx = 0; velIdx < dimWorld; ++velIdx)
171 for (
unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
175 source[Indices::velocity(velIdx)] += density_ *
dv2dx(globalPos)[velIdx][dimIdx];
178 source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_*
dvdx2(globalPos)[velIdx][dimIdx];
179 static const bool enableUnsymmetrizedVelocityGradient = getParam<bool>(
"FreeFlow.EnableUnsymmetrizedVelocityGradient",
false);
180 if (!enableUnsymmetrizedVelocityGradient)
181 source[Indices::velocity(velIdx)] -= density_ * kinematicViscosity_*
dvdx2(globalPos)[dimIdx][velIdx];
184 source[Indices::velocity(velIdx)] +=
dpdx(globalPos)[velIdx];
187 static const bool enableGravity = getParam<bool>(
"Problem.EnableGravity");
190 source[Indices::velocity(velIdx)] -= density_ * this->
gravity()[velIdx];
210 BoundaryTypes values;
213 values.setDirichlet(Indices::momentumXBalanceIdx);
227 const typename GridGeometry::LocalView& fvGeometry,
228 const typename GridGeometry::SubControlVolume& scv,
232 for (
const auto& scvf : scvfs(fvGeometry))
257 PrimaryVariables values;
258 values[Indices::pressureIdx] =
p(globalPos);
259 values[Indices::velocityXIdx] =
v(globalPos);
264 const DimVector
v(
const DimVector& globalPos)
const
267 v[0] = 2.0 * globalPos[0] * globalPos[0] * globalPos[0];
272 const DimMatrix
dvdx(
const DimVector& globalPos)
const
275 dvdx[0][0] = 6.0 * globalPos[0] * globalPos[0];
280 const DimMatrix
dv2dx(
const DimVector& globalPos)
const
283 for (
unsigned int velIdx = 0; velIdx < dimWorld; ++velIdx)
285 for (
unsigned int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
287 dv2dx[velIdx][dimIdx] =
dvdx(globalPos)[velIdx][dimIdx] *
v(globalPos)[dimIdx]
288 +
dvdx(globalPos)[dimIdx][dimIdx] *
v(globalPos)[velIdx];
295 const DimMatrix
dvdx2(
const DimVector& globalPos)
const
297 DimMatrix
dvdx2(0.0);
298 dvdx2[0][0] = 12.0 * globalPos[0];
303 const Scalar
p(
const DimVector& globalPos)
const
304 {
return 2.0 - 2.0 * globalPos[0]; }
307 const DimVector
dpdx(
const DimVector& globalPos)
const
336 return analyticalPressure_;
344 return analyticalVelocity_;
352 return analyticalVelocityOnFace_;
362 void createAnalyticalSolution_()
364 analyticalPressure_.resize(this->gridGeometry().numCellCenterDofs());
365 analyticalVelocity_.resize(this->gridGeometry().numCellCenterDofs());
366 analyticalVelocityOnFace_.resize(this->gridGeometry().numFaceDofs());
368 for (
const auto& element : elements(this->gridGeometry().gridView()))
370 auto fvGeometry =
localView(this->gridGeometry());
371 fvGeometry.bindElement(element);
372 for (
auto&& scv : scvs(fvGeometry))
374 auto ccDofIdx = scv.dofIndex();
375 auto ccDofPosition = scv.dofPosition();
379 for (
auto&& scvf : scvfs(fvGeometry))
381 const auto faceDofIdx = scvf.dofIndex();
382 const auto faceDofPosition = scvf.center();
383 const auto dirIdx = scvf.directionIndex();
385 analyticalVelocityOnFace_[faceDofIdx][dirIdx] = analyticalSolutionAtFace[Indices::velocity(dirIdx)];
388 analyticalPressure_[ccDofIdx] = analyticalSolutionAtCc[Indices::pressureIdx];
390 for(
int dirIdx = 0; dirIdx < ModelTraits::dim(); ++dirIdx)
391 analyticalVelocity_[ccDofIdx][dirIdx] = analyticalSolutionAtCc[Indices::velocity(dirIdx)];
399 Scalar kinematicViscosity_;
400 std::vector<Scalar> analyticalPressure_;
401 std::vector<GlobalPosition> analyticalVelocity_;
402 std::vector<GlobalPosition> analyticalVelocityOnFace_;
Setting constant fluid properties via the input file.
A liquid phase consisting of a single component.
GridCache::LocalView localView(const GridCache &gridCache)
Free function to get the local view of a grid cache object.
Definition: localview.hh:38
make the local view function available whenever we use the grid geometry
Definition: adapt.hh:29
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
The DUNE grid type.
Definition: common/properties.hh:57
Property to specify the type of a problem which has to be solved.
Definition: common/properties.hh:69
Definition: common/properties.hh:169
If disabled, the volume variables are not stored (reduces memory, but is slower)
Definition: common/properties.hh:178
specifies if data on flux vars should be saved (faster, but more memory consuming)
Definition: common/properties.hh:188
The type of the fluid system to use.
Definition: common/properties.hh:223
Returns whether to normalize the pressure term in the momentum balance or not.
Definition: common/properties.hh:346
Navier-Stokes problem base class.
Definition: dumux/freeflow/navierstokes/problem.hh:63
const GravityVector & gravity() const
Returns the acceleration due to gravity.
Definition: dumux/freeflow/navierstokes/problem.hh:134
bool enableInertiaTerms() const
Returns whether interia terms should be considered.
Definition: dumux/freeflow/navierstokes/problem.hh:140
A liquid phase consisting of a single component.
Definition: 1pliquid.hh:46
Test for the 1-D Navier-Stokes model with an analytical solution.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:89
const DimMatrix dvdx(const DimVector &globalPos) const
The velocity gradient.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:272
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
Returns the sources within the domain.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:157
NavierStokesAnalyticProblem(std::shared_ptr< const GridGeometry > gridGeometry)
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:108
const Scalar p(const DimVector &globalPos) const
The pressure.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:303
bool shouldWriteRestartFile() const
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:122
const DimMatrix dv2dx(const DimVector &globalPos) const
The gradient of the velocity squared (using product rule -> nothing to do here)
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:280
PrimaryVariables analyticalSolution(const GlobalPosition &globalPos) const
Returns the analytical solution of the problem at a given position.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:255
auto & getAnalyticalVelocitySolutionOnFace() const
Returns the analytical solution for the velocity at the faces.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:350
const DimVector dpdx(const DimVector &globalPos) const
The pressure gradient.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:307
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
Evaluates the initial value for a control volume.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:326
auto & getAnalyticalPressureSolution() const
Returns the analytical solution for the pressure.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:334
Scalar temperature() const
Returns the temperature within the domain in [K].
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:149
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
Specifies which kind of boundary condition should be used for which equation on a given boundary cont...
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:208
bool isDirichletCell(const Element &element, const typename GridGeometry::LocalView &fvGeometry, const typename GridGeometry::SubControlVolume &scv, int pvIdx) const
Returns whether a fixed Dirichlet value shall be used at a given cell.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:226
auto & getAnalyticalVelocitySolution() const
Returns the analytical solution for the velocity.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:342
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
Returns Dirichlet boundary values at a given position.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:244
const DimVector v(const DimVector &globalPos) const
The velocity.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:264
const DimMatrix dvdx2(const DimVector &globalPos) const
The gradient of the velocity gradient.
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:295
void printL2Error(const SolutionVector &curSol) const
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:127
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:48
std::tuple< NavierStokes, StaggeredFreeFlowModel > InheritsFrom
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:48
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:55
Dune::YaspGrid< 1 > type
Definition: test/freeflow/navierstokes/channel/1d/problem.hh:61
Routines to calculate the discrete L2 error.
Definition: l2error.hh:38
Defines a type tag and some properties for ree-flow models using the staggered scheme....
A single-phase, isothermal Navier-Stokes model.