diff --git a/dumux/discretization/box/subcontrolvolumeface.hh b/dumux/discretization/box/subcontrolvolumeface.hh
index 17a00578bf73e76c9e115d0d2991400da18b40fd..9eca2a0652c5053306c437733ce9f30412aae947 100644
--- a/dumux/discretization/box/subcontrolvolumeface.hh
+++ b/dumux/discretization/box/subcontrolvolumeface.hh
@@ -195,7 +195,7 @@ public:
return scvIndices_[1];
}
- //! The global index of this sub control volume face
+ //! The local index of this sub control volume face
GridIndexType index() const
{
return scvfIndex_;
diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
index bbe20095b9659c0e7b43a73a673987209371b2a4..916226f3ac35d970a215c1ae26e9e8cf90b73c88 100644
--- a/dumux/io/vtkoutputmodule.hh
+++ b/dumux/io/vtkoutputmodule.hh
@@ -81,7 +81,6 @@ class VtkOutputModule
using VolVarsVector = Dune::FieldVector<Scalar, dimWorld>;
static constexpr bool isBox = FVGridGeometry::discMethod == DiscretizationMethod::box;
- static constexpr int dofCodim = isBox ? dim : 0;
struct VolVarScalarDataInfo { std::function<Scalar(const VV&)> get; std::string name; };
struct VolVarVectorDataInfo { std::function<VolVarsVector(const VV&)> get; std::string name; };
@@ -258,7 +257,7 @@ private:
//! (1) Assemble all variable fields and add to writer
//////////////////////////////////////////////////////////////
- // instatiate the velocity output
+ // instantiate the velocity output
using VelocityVector = typename VelocityOutput::VelocityVector;
std::vector<VelocityVector> velocity(velocityOutput_->numPhases());
@@ -278,7 +277,7 @@ private:
|| addProcessRank)
{
const auto numCells = fvGridGeometry().gridView().size(0);
- const auto numDofs = numDofs_();
+ const auto numDofs = fvGridGeometry().numDofs();
// get fields for all volume variables
if (!volVarScalarDataInfo_.empty())
@@ -419,7 +418,7 @@ private:
writer_->clear();
}
- //! Assembles the fields and adds them to the writer (conforming output)
+ //! Assembles the fields and adds them to the writer (nonconforming output)
void writeNonConforming_(double time, Dune::VTK::OutputType type)
{
if(!isBox)
@@ -457,7 +456,7 @@ private:
|| addProcessRank)
{
const auto numCells = fvGridGeometry().gridView().size(0);
- const auto numDofs = numDofs_();
+ const auto numDofs = fvGridGeometry().numDofs();
// get fields for all volume variables
if (!volVarScalarDataInfo_.empty())
@@ -599,9 +598,6 @@ private:
template<class Vector, typename std::enable_if_t<!IsIndexable<decltype(std::declval<Vector>()[0])>::value, int> = 0>
std::size_t getNumberOfComponents_(const Vector& v) { return 1; }
- //! return the number of dofs, we only support vertex and cell data
- std::size_t numDofs_() const { return dofCodim == dim ? fvGridGeometry().vertexMapper().size() : fvGridGeometry().elementMapper().size(); }
-
const GridVariables& gridVariables_;
const SolutionVector& sol_;
diff --git a/dumux/porousmediumflow/velocityoutput.hh b/dumux/porousmediumflow/velocityoutput.hh
index 320cf6ff54454144c2cd3f06d3beb19352f27633..6df0c47e15a5fd29d1cc4f6016e1eba1c6bbe93f 100644
--- a/dumux/porousmediumflow/velocityoutput.hh
+++ b/dumux/porousmediumflow/velocityoutput.hh
@@ -182,15 +182,17 @@ public:
jacobianT1.mv(globalNormal, localNormal);
localNormal /= localNormal.two_norm();
- // insantiate the flux variables
+ // instantiate the flux variables
FluxVariables fluxVars;
fluxVars.init(problem_, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get the volume flux divided by the area of the
// subcontrolvolume face in the reference element
- // TODO: Divide by extrusion factor!!?
Scalar localArea = scvfReferenceArea_(geomType, scvf.index());
Scalar flux = fluxVars.advectiveFlux(phaseIdx, upwindTerm) / localArea;
+ flux /= problem_.extrusionFactor(element,
+ fvGeometry.scv(scvf.insideScvIdx()),
+ elementSolution(element, elemVolVars, fvGeometry));
// transform the volume flux into a velocity vector
Velocity tmpVelocity = localNormal;
@@ -230,17 +232,20 @@ public:
// first we extract the corner indices for each scv for the CIV method
// for network grids there might be multiple intersection with the same geometryInInside
- // we indentify those by the indexInInside for now (assumes conforming grids at branching facets)
+ // we identify those by the indexInInside for now (assumes conforming grids at branching facets)
// here we keep track of them
std::vector<bool> handledScvf;
- if (dim < dimWorld) handledScvf.resize(element.subEntities(1), false);
+ if (dim < dimWorld)
+ handledScvf.resize(element.subEntities(1), false);
// find the local face indices of the scvfs (for conforming meshes)
std::vector<unsigned int> scvfIndexInInside(fvGeometry.numScvf());
int localScvfIdx = 0;
for (const auto& intersection : intersections(fvGridGeometry_.gridView(), element))
{
- if (dim < dimWorld) if (handledScvf[intersection.indexInInside()]) continue;
+ if (dim < dimWorld)
+ if (handledScvf[intersection.indexInInside()])
+ continue;
if (intersection.neighbor() || intersection.boundary())
{
@@ -248,7 +253,8 @@ public:
scvfIndexInInside[localScvfIdx++] = intersection.indexInInside();
// for surface and network grids mark that we handled this face
- if (dim < dimWorld) handledScvf[intersection.indexInInside()] = true;
+ if (dim < dimWorld)
+ handledScvf[intersection.indexInInside()] = true;
}
}
@@ -297,8 +303,12 @@ public:
auto bcTypes = problemBoundaryTypes_(element, scvf);
if (bcTypes.hasNeumann())
{
+ // check if we have Neumann no flow, we can just use 0
+ const auto neumannFlux = problem_.neumann(element, fvGeometry, elemVolVars, scvf);
+ if (Dune::FloatCmp::eq<std::decay_t<decltype(neumannFlux)>, Dune::FloatCmp::CmpStyle::absolute>(neumannFlux, 0.0, 1e-30))
+ scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0;
// cubes
- if (dim == 1 || geomType.isCube())
+ else if (dim == 1 || geomType.isCube())
{
const auto fIdx = scvfIndexInInside[localScvfIdx];
const auto fIdxOpposite = fIdx%2 ? fIdx-1 : fIdx+1;
diff --git a/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt b/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt
index daa510481fba36bd9c362f5cad6951132263f0a6..30dae4f00b0ebc4df330b8e0f96e5d79926fff0b 100644
--- a/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/incompressible/CMakeLists.txt
@@ -1,14 +1,26 @@
dune_symlink_to_source_files(FILES "params.input")
+dune_symlink_to_source_files(FILES grids)
+
+# using tpfa and analytical Jacobian
+add_executable(test_1p_incompressible_tpfa_anadiff EXCLUDE_FROM_ALL main.cc)
+target_compile_definitions(test_1p_incompressible_tpfa_anadiff PUBLIC "TYPETAG=OnePIncompressibleTpfa" "NUMDIFFMETHOD=DiffMethod::analytic")
+
+# using mpfa and analytical Jacobian
+add_executable(test_1p_incompressible_mpfa_anadiff EXCLUDE_FROM_ALL main.cc)
+target_compile_definitions(test_1p_incompressible_mpfa_anadiff PUBLIC "TYPETAG=OnePIncompressibleMpfa" "NUMDIFFMETHOD=DiffMethod::analytic")
+
+# using box and analytical Jacobian
+add_executable(test_1p_incompressible_box_anadiff EXCLUDE_FROM_ALL main.cc)
+target_compile_definitions(test_1p_incompressible_box_anadiff PUBLIC "TYPETAG=OnePIncompressibleBox" "NUMDIFFMETHOD=DiffMethod::analytic")
# using tpfa and analytical Jacobian
dune_add_test(NAME test_1p_incompressible_tpfa
- SOURCES main.cc
- COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleTpfa NUMDIFFMETHOD=DiffMethod::analytic
+ TARGET test_1p_incompressible_tpfa_anadiff
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_1p_cc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_tpfa-00001.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_tpfa params.input -Problem.Name test_1p_incompressible_tpfa")
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_tpfa_anadiff params.input -Problem.Name test_1p_incompressible_tpfa")
# using tpfa, analytical Jacobian and quad precision
dune_add_test(NAME test_1p_incompressible_tpfa_quad
@@ -24,24 +36,21 @@ dune_add_test(NAME test_1p_incompressible_tpfa_quad
# using mpfa and analytical Jacobian
dune_add_test(NAME test_1p_incompressible_mpfa
- SOURCES main.cc
- COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleMpfa NUMDIFFMETHOD=DiffMethod::analytic
+ TARGET test_1p_incompressible_mpfa_anadiff
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_1p_cc-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_mpfa-00001.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_mpfa params.input -Problem.Name test_1p_incompressible_mpfa")
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_mpfa_anadiff params.input -Problem.Name test_1p_incompressible_mpfa")
# using box and analytical Jacobian
dune_add_test(NAME test_1p_incompressible_box
- SOURCES main.cc
- COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleBox NUMDIFFMETHOD=DiffMethod::analytic
+ TARGET test_1p_incompressible_box_anadiff
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_1p_box-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_box-00001.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_box params.input -Problem.Name test_1p_incompressible_box")
-
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_incompressible_box_anadiff params.input -Problem.Name test_1p_incompressible_box")
# using tpfa and numeric differentiation
dune_add_test(NAME test_1p_incompressible_tpfa_numdiff
@@ -67,4 +76,62 @@ dune_add_test(NAME test_1p_incompressible_box_numdiff
-Problem.Name test_1p_incompressible_box_numdiff
-Assembly.NumericDifference.PriVarMagnitude 1e5")
+# using tpfa and analytical Jacobian with extrusion factor
+dune_add_test(NAME test_1p_incompressible_tpfa_extrude
+ TARGET test_1p_incompressible_tpfa_anadiff
+ COMMAND ./test_1p_incompressible_tpfa_anadiff
+ CMD_ARGS params.input -Problem.Name test_1p_incompressible_tpfa_extrude
+ -Problem.ExtrusionFactor 10
+ -Vtk.AddVelocity 1
+ -Problem.CheckIsConstantVelocity true
+ -Problem.EnableGravity false)
+
+# using mpfa and analytical Jacobian with extrusion factor
+dune_add_test(NAME test_1p_incompressible_mpfa_extrude
+ TARGET test_1p_incompressible_mpfa_anadiff
+ COMMAND ./test_1p_incompressible_mpfa_anadiff
+ CMD_ARGS params.input -Problem.Name test_1p_incompressible_mpfa_extrude
+ -Problem.ExtrusionFactor 10
+ -Vtk.AddVelocity 1
+ -Problem.CheckIsConstantVelocity true
+ -Problem.EnableGravity false)
+
+# using box and analytical Jacobian with extrusion factor
+dune_add_test(NAME test_1p_incompressible_box_extrude
+ TARGET test_1p_incompressible_box_anadiff
+ COMMAND ./test_1p_incompressible_box_anadiff
+ CMD_ARGS params.input -Problem.Name test_1p_incompressible_box_extrude
+ -Problem.ExtrusionFactor 10
+ -Vtk.AddVelocity 1
+ -Problem.CheckIsConstantVelocity true
+ -Problem.EnableGravity false)
+
+# using box and analytical Jacobian with extrusion factor on distorted grid
+dune_add_test(NAME test_1p_incompressible_box_extrude_distorted
+ SOURCES main.cc
+ CMAKE_GUARD dune-uggrid_FOUND
+ COMPILE_DEFINITIONS GRIDTYPE=Dune::UGGrid<2>
+ COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleBox NUMDIFFMETHOD=DiffMethod::analytic
+ COMMAND ./test_1p_incompressible_box_extrude_distorted
+ CMD_ARGS params.input -Problem.Name test_1p_incompressible_box_extrude_distorted
+ -Problem.ExtrusionFactor 10
+ -Vtk.AddVelocity 1
+ -Problem.CheckIsConstantVelocity true
+ -Problem.EnableGravity false
+ -Grid.File ./grids/randomlydistorted.dgf)
+
+# using mpfa and analytical Jacobian with extrusion factor on distorted grid
+dune_add_test(NAME test_1p_incompressible_mpfa_extrude_distorted
+ SOURCES main.cc
+ CMAKE_GUARD dune-uggrid_FOUND
+ COMPILE_DEFINITIONS GRIDTYPE=Dune::UGGrid<2>
+ COMPILE_DEFINITIONS TYPETAG=OnePIncompressibleMpfa NUMDIFFMETHOD=DiffMethod::analytic
+ COMMAND ./test_1p_incompressible_mpfa_extrude_distorted
+ CMD_ARGS params.input -Problem.Name test_1p_incompressible_mpfa_extrude_distorted
+ -Problem.ExtrusionFactor 10
+ -Vtk.AddVelocity 1
+ -Problem.CheckIsConstantVelocity true
+ -Problem.EnableGravity false
+ -Grid.File ./grids/randomlydistorted.dgf)
+
set(CMAKE_BUILD_TYPE Release)
diff --git a/test/porousmediumflow/1p/implicit/incompressible/grids/randomlydistorted.dgf b/test/porousmediumflow/1p/implicit/incompressible/grids/randomlydistorted.dgf
new file mode 100644
index 0000000000000000000000000000000000000000..781f2bcdae5a184d58e3e09437b9797b6df101b7
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/incompressible/grids/randomlydistorted.dgf
@@ -0,0 +1,150 @@
+DGF
+Vertex
+0.0 0.0
+0.125 0.0
+0.25 0.0
+0.375 0.0
+0.5 0.0
+0.625 0.0
+0.75 0.0
+0.875 0.0
+1.0 0.0
+0.0 0.125
+0.0941385715747 0.138553971641
+0.221551511008 0.122209860023
+0.346853292055 0.129552124692
+0.482449274707 0.0983832818263
+0.651408372826 0.127916352327
+0.756728447773 0.137878324818
+0.846576088325 0.0965672249351
+1.0 0.125
+0.0 0.25
+0.131457007796 0.22608851401
+0.240456879904 0.271275285843
+0.345203098344 0.238509933903
+0.49793043015 0.222489717302
+0.653306271238 0.264447933308
+0.779921389529 0.237406544209
+0.867912433252 0.237206449115
+1.0 0.25
+0.0 0.375
+0.14939938624 0.3996354669
+0.225569987097 0.385474429209
+0.373490000873 0.405725027188
+0.525057845475 0.375418827149
+0.597759009065 0.351174934213
+0.731607015097 0.351029249509
+0.867855325892 0.347686472665
+1.0 0.375
+0.0 0.5
+0.099695749914 0.506732350839
+0.23976188556 0.478966982857
+0.369370419073 0.529651915591
+0.504659856464 0.482681252036
+0.605029156497 0.493909671035
+0.720928346604 0.529431790288
+0.888756786221 0.513641736326
+1.0 0.5
+0.0 0.625
+0.151921054497 0.64913187199
+0.270781382815 0.642285929151
+0.360323771837 0.633935978844
+0.475592672963 0.616708801221
+0.642406786794 0.622788518651
+0.763924855322 0.620843310615
+0.846036926139 0.612541711404
+1.0 0.625
+0.0 0.75
+0.106895114158 0.744668125056
+0.274954350162 0.740154282669
+0.361024506709 0.779420524896
+0.519113829969 0.73378760529
+0.654306310229 0.74963785855
+0.744256881762 0.76984643781
+0.864159923638 0.766042925423
+1.0 0.75
+0.0 0.875
+0.121631522308 0.87931798433
+0.23413160962 0.87565218627
+0.353753336001 0.870490662736
+0.523652797146 0.84803633498
+0.634049374115 0.904248570954
+0.750680324251 0.85884675384
+0.856851940426 0.889040325118
+1.0 0.875
+0.0 1.0
+0.125 1.0
+0.25 1.0
+0.375 1.0
+0.5 1.0
+0.625 1.0
+0.75 1.0
+0.875 1.0
+1.0 1.0
+#
+CUBE
+0 1 9 10
+1 2 10 11
+2 3 11 12
+3 4 12 13
+4 5 13 14
+5 6 14 15
+6 7 15 16
+7 8 16 17
+9 10 18 19
+10 11 19 20
+11 12 20 21
+12 13 21 22
+13 14 22 23
+14 15 23 24
+15 16 24 25
+16 17 25 26
+18 19 27 28
+19 20 28 29
+20 21 29 30
+21 22 30 31
+22 23 31 32
+23 24 32 33
+24 25 33 34
+25 26 34 35
+27 28 36 37
+28 29 37 38
+29 30 38 39
+30 31 39 40
+31 32 40 41
+32 33 41 42
+33 34 42 43
+34 35 43 44
+36 37 45 46
+37 38 46 47
+38 39 47 48
+39 40 48 49
+40 41 49 50
+41 42 50 51
+42 43 51 52
+43 44 52 53
+45 46 54 55
+46 47 55 56
+47 48 56 57
+48 49 57 58
+49 50 58 59
+50 51 59 60
+51 52 60 61
+52 53 61 62
+54 55 63 64
+55 56 64 65
+56 57 65 66
+57 58 66 67
+58 59 67 68
+59 60 68 69
+60 61 69 70
+61 62 70 71
+63 64 72 73
+64 65 73 74
+65 66 74 75
+66 67 75 76
+67 68 76 77
+68 69 77 78
+69 70 78 79
+70 71 79 80
+#
diff --git a/test/porousmediumflow/1p/implicit/incompressible/main.cc b/test/porousmediumflow/1p/implicit/incompressible/main.cc
index 830c99c4548b27dee670d99e1e3b86a6512b8fb4..290bf8af22974f260dc99c84cd53a7055bbdf668 100644
--- a/test/porousmediumflow/1p/implicit/incompressible/main.cc
+++ b/test/porousmediumflow/1p/implicit/incompressible/main.cc
@@ -29,6 +29,7 @@
// Support for quad precision has to be included before any other Dune module:
#include <dumux/common/quad.hh>
+#include <dune/common/float_cmp.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
@@ -72,8 +73,8 @@ int main(int argc, char** argv) try
//////////////////////////////////////////////////////////////////////
// try to create a grid (from the given grid file or the input file)
/////////////////////////////////////////////////////////////////////
-
- GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
+ using Grid = GetPropType<TypeTag, Properties::Grid>;
+ GridManager<Grid> gridManager;
gridManager.init();
// we compute on the leaf grid view
@@ -144,6 +145,46 @@ int main(int argc, char** argv) try
timer.stop();
+ const bool checkIsConstantVelocity = getParam<bool>("Problem.CheckIsConstantVelocity", false);
+ if(checkIsConstantVelocity)
+ {
+ // instantiate the velocity output
+ VelocityOutput velocityOutput(*gridVariables);
+ using VelocityVector = typename VelocityOutput::VelocityVector;
+ VelocityVector velocity;
+
+ constexpr bool isBox = FVGridGeometry::discMethod == Dumux::DiscretizationMethod::box;
+ constexpr int dimWorld = FVGridGeometry::GridView::dimensionworld;
+ const auto numCells = leafGridView.size(0);
+ const auto numDofs = fvGridGeometry->numDofs();
+ auto numVelocities = (isBox && dimWorld == 1) ? numCells : numDofs;
+
+ velocity.resize(numVelocities);
+
+ const auto exactVel = problem->velocity();
+
+ for (const auto& element : elements(leafGridView, Dune::Partitions::interior))
+ {
+ const auto eIdx = fvGridGeometry->elementMapper().index(element);
+
+ auto fvGeometry = localView(*fvGridGeometry);
+ auto elemVolVars = localView(gridVariables->curGridVolVars());
+
+ fvGeometry.bind(element);
+ elemVolVars.bind(element, fvGeometry, x);
+
+ velocityOutput.calculateVelocity(velocity, elemVolVars, fvGeometry, element, 0);
+
+ using Scalar = Grid::ctype;
+ // the y-component of the velocity should be exactly reproduced
+ // the x-component should be zero
+ // use a relative comparison for the y-component and an absolute one for the x-component
+ if(Dune::FloatCmp::ne(velocity[eIdx][dimWorld-1], exactVel[dimWorld-1], /*eps*/1e-8) ||
+ Dune::FloatCmp::ne<Scalar, Dune::FloatCmp::absolute>(velocity[eIdx][0], exactVel[0], /*eps*/1e-10))
+ DUNE_THROW(Dune::InvalidStateException, "Velocity is not exactly reproduced");
+ }
+ }
+
const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
if (mpiHelper.rank() == 0)
std::cout << "Simulation took " << timer.elapsed() << " seconds on "
diff --git a/test/porousmediumflow/1p/implicit/incompressible/params.input b/test/porousmediumflow/1p/implicit/incompressible/params.input
index ebf9ad8f55e047c86db93031d8d490fd57cc9168..1e734499f1c5ef5777d79d7619079cd3b01676dd 100644
--- a/test/porousmediumflow/1p/implicit/incompressible/params.input
+++ b/test/porousmediumflow/1p/implicit/incompressible/params.input
@@ -1,5 +1,6 @@
[Problem]
Name = incompressible
+ExtrusionFactor = 1
[Grid]
LowerLeft = 0 0
diff --git a/test/porousmediumflow/1p/implicit/incompressible/problem.hh b/test/porousmediumflow/1p/implicit/incompressible/problem.hh
index 4a505b6d2aa4622e73e2f85ec4f72480cc4d509c..a26e0cca7c2abd42795b246ad52b73c8c2efd63a 100644
--- a/test/porousmediumflow/1p/implicit/incompressible/problem.hh
+++ b/test/porousmediumflow/1p/implicit/incompressible/problem.hh
@@ -24,6 +24,9 @@
#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_PROBLEM_HH
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
#include <dune/grid/yaspgrid.hh>
#include <dumux/common/quad.hh>
@@ -40,6 +43,10 @@
#include "spatialparams.hh"
+#ifndef GRIDTYPE // default to yasp grid if not provided by CMake
+#define GRIDTYPE Dune::YaspGrid<2>
+#endif
+
namespace Dumux
{
// forward declarations
@@ -58,7 +65,7 @@ struct OnePIncompressibleBox { using InheritsFrom = std::tuple<OnePIncompressibl
// Set the grid type
template<class TypeTag>
-struct Grid<TypeTag, TTag::OnePIncompressible> { using type = Dune::YaspGrid<2>; };
+struct Grid<TypeTag, TTag::OnePIncompressible> { using type = GRIDTYPE; };
// Set the problem type
template<class TypeTag>
@@ -118,13 +125,26 @@ class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
static constexpr int dimWorld = GridView::dimensionworld;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using GlobalPosition = Dune::FieldVector<Scalar,dimWorld>;
public:
OnePTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
- : ParentType(fvGridGeometry) {}
+ : ParentType(fvGridGeometry), velocity_(0.0)
+ {
+ extrusionFactor_ = getParam<Scalar>("Problem.ExtrusionFactor");
+ Scalar permeability = getParam<Scalar>("SpatialParams.Permeability");
+ dp_dy_ = -1.0e+5;
+
+ const bool checkIsConstantVelocity = getParam<bool>("Problem.CheckIsConstantVelocity", false);
+ if(checkIsConstantVelocity)
+ {
+ velocity_[dimWorld-1] = -permeability * dp_dy_;
+ velocity_[dimWorld-1] /= FluidSystem::viscosity(temperature(), 1.0e5);
+ }
+ }
/*!
* \brief Specifies which kind of boundary condition should be
@@ -158,7 +178,8 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0);
- values[0] = 1.0e+5*(2.0 - globalPos[dimWorld-1]);
+ values[0] = 1.0e+5 + dp_dy_*(globalPos[dimWorld-1] - this->fvGridGeometry().bBoxMax()[dimWorld-1]);
+
return values;
}
@@ -173,6 +194,33 @@ public:
{
return 283.15; // 10°C
}
+
+ /*!
+ * \brief Return how much the domain is extruded at a given position.
+ *
+ * This means the factor by which a lower-dimensional
+ * entity needs to be expanded to get a full dimensional cell.
+ */
+ Scalar extrusionFactorAtPos(const GlobalPosition &globalPos) const
+ {
+ return extrusionFactor_;
+ }
+
+ /*!
+ * \brief Returns the velocity
+ *
+ * The velocity is given for the case of a linear pressure solution
+ * with constant permeablity and without gravity
+ */
+ const GlobalPosition velocity() const
+ {
+ return velocity_;
+ }
+
+private:
+ Scalar extrusionFactor_;
+ Scalar dp_dy_;
+ GlobalPosition velocity_;
};
} // end namespace Dumux
diff --git a/test/porousmediumflow/1p/implicit/incompressible/spatialparams.hh b/test/porousmediumflow/1p/implicit/incompressible/spatialparams.hh
index 2e36b043f1f586e52e867c300bbee447bf32d87f..9ebcaf99ebfd4a14d9af0cc21e3afbb21f8969cc 100644
--- a/test/porousmediumflow/1p/implicit/incompressible/spatialparams.hh
+++ b/test/porousmediumflow/1p/implicit/incompressible/spatialparams.hh
@@ -59,6 +59,8 @@ public:
lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
+
+ checkIsConstantVelocity_ = getParam<bool>("Problem.CheckIsConstantVelocity",false);
}
/*!
@@ -74,7 +76,7 @@ public:
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
- if (isInLens_(scv.dofPosition()))
+ if (isInLens_(scv.dofPosition()) && !checkIsConstantVelocity_)
return permeabilityLens_;
else
return permeability_;
@@ -104,6 +106,8 @@ private:
Scalar permeability_, permeabilityLens_;
static constexpr Scalar eps_ = 1.5e-7;
+
+ bool checkIsConstantVelocity_;
};
} // end namespace Dumux