From c4e37299bf70a5b2afa523df73cd8bfaa6c92d2c Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Wed, 14 Dec 2016 23:17:40 +0100
Subject: [PATCH] [test][1p] Add 1d3d network test cc. TODO fix box.
---
.../1p/implicit/CMakeLists.txt | 23 +-
.../1p/implicit/test_box1p1d3d.cc | 69 ----
.../1p/implicit/test_box1p1d3d.input | 19 -
.../1p/implicit/test_box1pnetwork1d3d.cc | 38 ++
.../1p/implicit/test_box1pnetwork1d3d.input | 11 +
.../1p/implicit/test_cc1p1d3d.cc | 68 ----
.../1p/implicit/test_cc1p1d3d.input | 19 -
.../1p/implicit/test_cc1pnetwork1d3d.cc | 38 ++
.../1p/implicit/test_cc1pnetwork1d3d.input | 11 +
.../1p/implicit/tubesconvergencetest.py | 62 ++++
.../1p/implicit/tubesproblem.hh | 324 ++++++++++++++++++
.../1p/implicit/tubesspatialparams.hh | 104 ++++++
12 files changed, 600 insertions(+), 186 deletions(-)
delete mode 100644 test/porousmediumflow/1p/implicit/test_box1p1d3d.cc
delete mode 100644 test/porousmediumflow/1p/implicit/test_box1p1d3d.input
create mode 100644 test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.cc
create mode 100644 test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.input
delete mode 100644 test/porousmediumflow/1p/implicit/test_cc1p1d3d.cc
delete mode 100644 test/porousmediumflow/1p/implicit/test_cc1p1d3d.input
create mode 100644 test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.cc
create mode 100644 test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.input
create mode 100755 test/porousmediumflow/1p/implicit/tubesconvergencetest.py
create mode 100644 test/porousmediumflow/1p/implicit/tubesproblem.hh
create mode 100644 test/porousmediumflow/1p/implicit/tubesspatialparams.hh
diff --git a/test/porousmediumflow/1p/implicit/CMakeLists.txt b/test/porousmediumflow/1p/implicit/CMakeLists.txt
index 2fdc06850d..93b23efd35 100644
--- a/test/porousmediumflow/1p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/CMakeLists.txt
@@ -3,7 +3,11 @@ add_subdirectory(pointsources)
add_input_file_links()
add_gstat_file_links()
-add_executable(test_mpfa1p test_mpfa1p.cc)
+dune_symlink_to_source_files(FILES "tubesconvergencetest.py")
+
+set(CMAKE_BUILD_TYPE Debug)
+
+add_executable(test_mpfa1p EXCLUDE_FROM_ALL test_mpfa1p.cc)
# isothermal tests
add_dumux_test(test_box1p test_box1p test_box1p.cc
@@ -70,12 +74,12 @@ add_dumux_test(test_cc1pniconvection test_cc1pniconvection test_cc1pniconvection
--command "${CMAKE_CURRENT_BINARY_DIR}/test_cc1pniconvection")
# dim < dimWorld tests with foamgrid
-add_dumux_test(test_box1p1d3d test_box1p1d3d test_box1p1d3d.cc
+add_dumux_test(test_box1pnetwork1d3d test_box1pnetwork1d3d test_box1pnetwork1d3d.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
--script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/1ptestbox1d3d-reference.vtp
- ${CMAKE_CURRENT_BINARY_DIR}/1ptestbox1d3d-00001.vtp
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_box1p1d3d")
+ --files ${CMAKE_SOURCE_DIR}/test/references/1pboxnetwork1d3d-reference.vtp
+ ${CMAKE_CURRENT_BINARY_DIR}/1pboxnetwork1d3d-00001.vtp
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_box1pnetwork1d3d")
add_dumux_test(test_box1pfracture2d3d test_box1pfracture2d3d test_box1pfracture2d3d.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
@@ -84,12 +88,9 @@ add_dumux_test(test_box1pfracture2d3d test_box1pfracture2d3d test_box1pfracture2
${CMAKE_CURRENT_BINARY_DIR}/1pboxfracture2d3d-00001.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_box1pfracture2d3d")
-add_dumux_test(test_cc1p1d3d test_cc1p1d3d test_cc1p1d3d.cc
- python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/1ptestcc1d3d-reference.vtp
- ${CMAKE_CURRENT_BINARY_DIR}/1ptestcc1d3d-00001.vtp
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_cc1p1d3d")
+dune_add_test(SOURCES test_cc1pnetwork1d3d.cc
+ NAME test_cc1pnetwork1d3d
+ COMMAND ./tubesconvergencetest.py test_cc1pnetwork1d3d)
add_dumux_test(test_cc1pfracture2d3d test_cc1pfracture2d3d test_cc1pfracture2d3d.cc
python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
diff --git a/test/porousmediumflow/1p/implicit/test_box1p1d3d.cc b/test/porousmediumflow/1p/implicit/test_box1p1d3d.cc
deleted file mode 100644
index 9afff531b9..0000000000
--- a/test/porousmediumflow/1p/implicit/test_box1p1d3d.cc
+++ /dev/null
@@ -1,69 +0,0 @@
-// -*- 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 2 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
- *
- * \brief test for the one-phase CC model
- */
-#include <config.h>
-#include "1ptestproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \brief Provides an interface for customizing error messages associated with
- * reading in parameters.
- *
- * \param progName The name of the program, that was tried to be started.
- * \param errorMsg The error message that was issued by the start function.
- * Comprises the thing that went wrong and a general help message.
- */
-void usage(const char *progName, const std::string &errorMsg)
-{
- if (errorMsg.size() > 0) {
- std::string errorMessageOut = "\nUsage: ";
- errorMessageOut += progName;
- errorMessageOut += " [options]\n";
- errorMessageOut += errorMsg;
- errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
- "\t-TimeManager.TEnd End of the simulation [s] \n"
- "\t-TimeManager.DtInitial Initial timestep size [s] \n"
- "\t-Grid.LowerLeft Lower left corner coordinates\n"
- "\t-Grid.UpperRight Upper right corner coordinates\n"
- "\t-Grid.Cells Number of cells in respective coordinate directions\n"
- "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
- "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
- "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
- "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
-
- std::cout << errorMessageOut
- << "\n";
- }
-}
-
-int main(int argc, char** argv)
-{
-#if HAVE_DUNE_FOAMGRID
- typedef TTAG(OnePOneDThreeDTestBoxProblem) ProblemTypeTag;
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
-#else
-#warning External grid module dune-foamgrid needed to run this example.
- std::cerr << "Test skipped, it needs dune-foamgrid!" << std::endl;
- return 77;
-#endif
-}
diff --git a/test/porousmediumflow/1p/implicit/test_box1p1d3d.input b/test/porousmediumflow/1p/implicit/test_box1p1d3d.input
deleted file mode 100644
index 8c972d99f6..0000000000
--- a/test/porousmediumflow/1p/implicit/test_box1p1d3d.input
+++ /dev/null
@@ -1,19 +0,0 @@
-[TimeManager]
-DtInitial = 1 # [s]
-TEnd = 1 # [s]
-
-[Grid]
-LowerLeft = 0 0 0
-UpperRight = 1 1 1
-Cells = 4
-
-[Problem]
-Name = 1ptestbox1d3d # name passed to the output routines
-
-[SpatialParams]
-LensLowerLeft = 0.2 0.2 0.2
-LensUpperRight = 0.8 0.8 0.8
-
-Permeability = 1e-10 # [m^2]
-PermeabilityLens = 1e-12 # [m^2]
-
diff --git a/test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.cc b/test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.cc
new file mode 100644
index 0000000000..0d698cbb6c
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.cc
@@ -0,0 +1,38 @@
+// -*- 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 2 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
+ *
+ * \brief test for the one-phase CC model
+ */
+#include <config.h>
+#include "tubesproblem.hh"
+#include <dumux/common/start.hh>
+
+int main(int argc, char** argv)
+{
+#if HAVE_DUNE_FOAMGRID
+ typedef TTAG(TubesTestBoxProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, [](const char *, const std::string &){});
+#else
+#warning External grid module dune-foamgrid needed to run this example.
+ std::cerr << "Test skipped, it needs dune-foamgrid!" << std::endl;
+ return 77;
+#endif
+}
diff --git a/test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.input b/test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.input
new file mode 100644
index 0000000000..3fc22439d7
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/test_box1pnetwork1d3d.input
@@ -0,0 +1,11 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 1 # [s]
+
+[Grid]
+File = ./grids/network1d3d.msh
+
+[Problem]
+Name = 1pccnetwork1d3d
+LiquidDensity = 1050 # blood density
+LiquidKinematicViscosity = 4.0e-3 # blood viscosity
diff --git a/test/porousmediumflow/1p/implicit/test_cc1p1d3d.cc b/test/porousmediumflow/1p/implicit/test_cc1p1d3d.cc
deleted file mode 100644
index fd20762ad8..0000000000
--- a/test/porousmediumflow/1p/implicit/test_cc1p1d3d.cc
+++ /dev/null
@@ -1,68 +0,0 @@
-// -*- 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 2 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
- *
- * \brief test for the one-phase CC model
- */
-#include <config.h>
-#include "1ptestproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \brief Provides an interface for customizing error messages associated with
- * reading in parameters.
- *
- * \param progName The name of the program, that was tried to be started.
- * \param errorMsg The error message that was issued by the start function.
- * Comprises the thing that went wrong and a general help message.
- */
-void usage(const char *progName, const std::string &errorMsg)
-{
- if (errorMsg.size() > 0) {
- std::string errorMessageOut = "\nUsage: ";
- errorMessageOut += progName;
- errorMessageOut += " [options]\n";
- errorMessageOut += errorMsg;
- errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
- "\t-TimeManager.TEnd End of the simulation [s] \n"
- "\t-TimeManager.DtInitial Initial timestep size [s] \n"
- "\t-Grid.File Name of the file containing the grid \n"
- "\t definition in DGF format\n"
- "\t-SpatialParams.LensLowerLeftX coordinates of the lower left corner of the lens [m] \n"
- "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
- "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
- "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
-
- std::cout << errorMessageOut
- << "\n";
- }
-}
-
-int main(int argc, char** argv)
-{
-#if HAVE_DUNE_FOAMGRID
- typedef TTAG(OnePOneDThreeDTestCCProblem) ProblemTypeTag;
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
-#else
-#warning External grid module dune-foamgrid needed to run this example.
- std::cerr << "Test skipped, it needs dune-foamgrid!" << std::endl;
- return 77;
-#endif
-}
diff --git a/test/porousmediumflow/1p/implicit/test_cc1p1d3d.input b/test/porousmediumflow/1p/implicit/test_cc1p1d3d.input
deleted file mode 100644
index cb05c034f6..0000000000
--- a/test/porousmediumflow/1p/implicit/test_cc1p1d3d.input
+++ /dev/null
@@ -1,19 +0,0 @@
-[TimeManager]
-DtInitial = 1 # [s]
-TEnd = 1 # [s]
-
-[Grid]
-LowerLeft = 0 0 0
-UpperRight = 1 1 1
-Cells = 4
-
-[Problem]
-Name = 1ptestcc1d3d # name passed to the output routines
-
-[SpatialParams]
-LensLowerLeft = 0.2 0.2 0.2
-LensUpperRight = 0.8 0.8 0.8
-
-Permeability = 1e-10 # [m^2]
-PermeabilityLens = 1e-12 # [m^2]
-
diff --git a/test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.cc b/test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.cc
new file mode 100644
index 0000000000..21253b0753
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.cc
@@ -0,0 +1,38 @@
+// -*- 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 2 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
+ *
+ * \brief test for the one-phase CC model
+ */
+#include <config.h>
+#include "tubesproblem.hh"
+#include <dumux/common/start.hh>
+
+int main(int argc, char** argv)
+{
+#if HAVE_DUNE_FOAMGRID
+ typedef TTAG(TubesTestCCTpfaProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, [](const char *, const std::string &){});
+#else
+#warning External grid module dune-foamgrid needed to run this example.
+ std::cerr << "Test skipped, it needs dune-foamgrid!" << std::endl;
+ return 77;
+#endif
+}
diff --git a/test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.input b/test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.input
new file mode 100644
index 0000000000..3fc22439d7
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/test_cc1pnetwork1d3d.input
@@ -0,0 +1,11 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 1 # [s]
+
+[Grid]
+File = ./grids/network1d3d.msh
+
+[Problem]
+Name = 1pccnetwork1d3d
+LiquidDensity = 1050 # blood density
+LiquidKinematicViscosity = 4.0e-3 # blood viscosity
diff --git a/test/porousmediumflow/1p/implicit/tubesconvergencetest.py b/test/porousmediumflow/1p/implicit/tubesconvergencetest.py
new file mode 100755
index 0000000000..8a2723208a
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/tubesconvergencetest.py
@@ -0,0 +1,62 @@
+#!/usr/bin/env python
+
+from math import *
+import subprocess
+import sys
+
+if len(sys.argv) != 2:
+ sys.stderr.write('Please provide a single argument <testname> to the script\n')
+ sys.exit(1)
+
+testname = str(sys.argv[1])
+
+# remove the old log file
+subprocess.call(['rm', testname + '.log'])
+print("Removed old log file ({})!".format(testname + '.log'))
+
+# do the runs with different refinement
+for i in [0, 1, 2, 3, 4, 5]:
+ subprocess.call(['./' + testname, '-Grid.Refinement', str(i),
+ '-Problem.Name', testname])
+
+# check the rates and append them to the log file
+logfile = open(testname + '.log', "r+")
+
+error = []
+hmax = []
+for line in logfile:
+ line = line.strip("\n")
+ line = line.strip("\[ConvergenceTest\]")
+ line = line.split()
+ error.append(float(line[2]))
+ hmax.append(float(line[5]))
+
+results = []
+logfile.truncate(0)
+logfile.write("n\thmax\t\terror\t\trate\n")
+logfile.write("-"*50 + "\n")
+for i in range(len(error)-1):
+ if isnan(error[i]) or isinf(error[i]):
+ continue
+ if not (error[i] < 1e-12 or error[i] < 1e-12):
+ rate = (log(error[i])-log(error[i+1]))/(log(hmax[i])-log(hmax[i+1]))
+ message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, hmax[i], error[i], rate)
+ logfile.write(message)
+ results.append(rate)
+ else:
+ logfile.write("error: exact solution!?")
+i = len(error)-1
+message = "{}\t{:0.4e}\t{:0.4e}\n\n\n".format(i, hmax[i], error[i])
+logfile.write(message)
+
+logfile.close()
+print("\nComputed the following convergence rates for {}:\n".format(testname))
+subprocess.call(['cat', testname + '.log'])
+
+# check the rates, we expect rates around 2
+for r in results:
+ if int(round(r)) is not 2:
+ sys.stderr.write("*"*70 + "\n", "The convergence rates were not close enough to 2! Test failed.\n", "*"*70 + "\n")
+ sys.exit(1)
+
+sys.exit(0)
diff --git a/test/porousmediumflow/1p/implicit/tubesproblem.hh b/test/porousmediumflow/1p/implicit/tubesproblem.hh
new file mode 100644
index 0000000000..be741ac146
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/tubesproblem.hh
@@ -0,0 +1,324 @@
+// -*- 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 2 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
+ *
+ * \brief A test problem for the 1p model. A pipe system with circular cross-section
+ * and a branching point embedded in a three-dimensional world
+ */
+#ifndef DUMUX_ONEP_TUBES_TEST_PROBLEM_HH
+#define DUMUX_ONEP_TUBES_TEST_PROBLEM_HH
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/1p/implicit/model.hh>
+#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/material/components/constant.hh>
+
+#include "tubesspatialparams.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class TubesTestProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TubesTestProblem, INHERITS_FROM(OneP));
+NEW_TYPE_TAG(TubesTestCCTpfaProblem, INHERITS_FROM(CCTpfaModel, TubesTestProblem));
+NEW_TYPE_TAG(TubesTestBoxProblem, INHERITS_FROM(BoxModel, TubesTestProblem));
+
+// Set the grid type
+SET_TYPE_PROP(TubesTestProblem, Grid, Dune::FoamGrid<1, 3>);
+
+// Set the problem property
+SET_TYPE_PROP(TubesTestProblem, Problem, TubesTestProblem<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(TubesTestProblem, SpatialParams, TubesTestSpatialParams<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(TubesTestProblem, LinearSolver, UMFPackBackend<TypeTag>);
+
+// the fluid
+SET_PROP(TubesTestProblem, Fluid)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+public:
+ using type = FluidSystems::LiquidPhase<Scalar, Constant<TypeTag, Scalar>>;
+};
+}
+
+/*!
+ * \ingroup OnePModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief A test problem for the 1p model. A pipe system with circular cross-section
+ * and a branching point embedded in a three-dimensional world
+ */
+template <class TypeTag>
+class TubesTestProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+ using ParentType = ImplicitPorousMediaProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+
+ // Grid and world dimension
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ enum {
+ // indices of the primary variables
+ conti0EqIdx = Indices::conti0EqIdx,
+ pressureIdx = Indices::pressureIdx
+ };
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
+public:
+ TubesTestProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+
+ //get hMax_ of the grid
+ hMax_ = 0.0;
+ for (const auto& element : elements(this->gridView()))
+ hMax_ = std::max(element.geometry().volume(), hMax_);
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const std::string& name() const
+ {
+ return name_;
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 273.15 + 37.0; } // Body temperature
+
+ /*!
+ * \brief Return how much the domain is extruded at a given sub-control volume.
+ *
+ * This means the factor by which a lower-dimensional (1D or 2D)
+ * entity needs to be expanded to get a full dimensional cell. The
+ * default is 1.0 which means that 1D problems are actually
+ * thought as pipes with a cross section of 1 m^2 and 2D problems
+ * are assumed to extend 1 m to the back.
+ */
+ Scalar boxExtrusionFactor(const Element &element,
+ const SubControlVolume &scv) const
+ {
+ const auto radius = this->spatialParams().radius(scv);
+ return M_PI*radius*radius;
+ }
+
+ // \}
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param globalPos The position of the center of the finite volume
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes bcTypes;
+ bcTypes.setAllDirichlet();
+ return bcTypes;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * control volume.
+ *
+ * \param globalPos The center of the finite volume which ought to be set.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a priVars parameter stores the mass flux
+ * in normal direction of each component. Negative values mean
+ * influx.
+ */
+ PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * This is the method for the case where the source term is
+ * potentially solution dependent and requires some quantities that
+ * are specific to the fully-implicit method.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param scv The sub control volume
+ *
+ * For this method, the \a values parameter stores the conserved quantity rate
+ * generated or annihilate per volume unit. Positive values mean
+ * that the conserved quantity is created, negative ones mean that it vanishes.
+ * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
+ */
+ PrimaryVariables source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ PrimaryVariables source(0.0);
+
+ const auto& globalPos = scv.center();
+ const auto& volVars = elemVolVars[scv];
+ const auto K = this->spatialParams().intrinsicPermeability(scv, elemVolVars[scv]);
+
+ using std::sin;
+ if (globalPos[2] > 0.5 - eps_)
+ source[conti0EqIdx] = K/volVars.viscosity()*volVars.density()
+ *4.0*4.0*M_PI*M_PI*sin(4.0*M_PI*globalPos[2]);
+ else
+ // the "/3.0" stems from the coordindate transformations on the lower branches
+ source[conti0EqIdx] = K/volVars.viscosity()*volVars.density()
+ *4.0*4.0*M_PI*M_PI*sin(4.0*M_PI*globalPos[2])/3.0;
+
+ return source;
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initial(const SubControlVolume& scv) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ // \}
+
+ void postTimeStep()
+ {
+ const auto& solution = this->model().curSol();
+
+ // calculate the discrete L2-Norm
+ Scalar LTwoNorm = 0.0;
+
+ // get the Gaussian quadrature rule for intervals
+ const auto& quad = Dune::QuadratureRules<Scalar, dim>::rule(Dune::GeometryType(1), 1);
+
+ for (const auto& element : elements(this->gridView()))
+ {
+ const unsigned int eIdx = this->elementMapper().index(element);
+ const auto geometry = element.geometry();
+ for(auto&& qp : quad)
+ {
+ const auto globalPos = geometry.global(qp.position());
+ const Scalar pe = exactPressure_(globalPos);
+ const Scalar integrationElement = geometry.integrationElement(qp.position());
+ Scalar p = 0.0;
+ if (!isBox)
+ p = solution[eIdx][pressureIdx];
+ else
+ {
+ // do interpolation with ansatz functions
+ std::vector<Dune::FieldVector<Scalar, 1> > shapeValues;
+ const auto& localFiniteElement = feCache_.get(geometry.type());
+ localFiniteElement.localBasis().evaluateFunction(qp.position(), shapeValues);
+ for (unsigned int i = 0; i < shapeValues.size(); ++i)
+ p += shapeValues[i]*solution[this->model().dofMapper().subIndex(element, i, dim)][pressureIdx];
+ }
+ LTwoNorm += (p - pe)*(p - pe)*qp.weight()*integrationElement;
+ }
+ }
+ LTwoNorm = std::sqrt(LTwoNorm);
+
+ // write the norm into a log file
+ logFile_.open(this->name() + ".log", std::ios::app);
+ logFile_ << "[ConvergenceTest] L2-norm(pressure) = " << LTwoNorm << " hMax = " << hMax_ << std::endl;
+ logFile_.close();
+ }
+
+private:
+
+ Scalar exactPressure_(const GlobalPosition& globalPos)
+ {
+ using std::sin;
+ return sin(4.0*M_PI*globalPos[2]);
+ }
+
+ static constexpr Scalar eps_ = 1e-8;
+ std::string name_;
+
+ Scalar hMax_;
+
+ std::ofstream logFile_;
+ typename Dune::PQkLocalFiniteElementCache<Scalar, Scalar, dim, 1> feCache_;
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/1p/implicit/tubesspatialparams.hh b/test/porousmediumflow/1p/implicit/tubesspatialparams.hh
new file mode 100644
index 0000000000..c7ab856603
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/tubesspatialparams.hh
@@ -0,0 +1,104 @@
+// -*- 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 2 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
+ *
+ * \brief A test problem for the 1p model. A pipe system with circular cross-section
+ * and a branching point embedded in a three-dimensional world
+ */
+#ifndef DUMUX_ONEP_TUBES_TEST_SPATIALPARAMS_HH
+#define DUMUX_ONEP_TUBES_TEST_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/implicit1p.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup OnePModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief A test problem for the 1p model. A pipe system with circular cross-section
+ * and a branching point embedded in a three-dimensional world
+ */
+template<class TypeTag>
+class TubesTestSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
+{
+ using ParentType = ImplicitSpatialParamsOneP<TypeTag>;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+
+public:
+ TubesTestSpatialParams(const Problem& problem, const GridView& gridView)
+ : ParentType(problem, gridView)
+ {
+ radius_ = 1.0;
+
+ using std::sqrt;
+ radiusMain_ = sqrt(sqrt(4.0/sqrt(3.0)));
+ }
+
+ /*!
+ * \brief Return the radius of the circular pipe for the current sub-control volume in [m].
+ *
+ * \param scv The sub control volume
+ */
+ Scalar radius(const SubControlVolume &scv) const
+ {
+ if(scv.center()[2] > 0.5 - eps_)
+ return radiusMain_;
+ else
+ return radius_;
+ }
+
+ /*!
+ * \brief Function for defining the intrinsic (absolute) permeability.
+ *
+ * \param scv The sub control volume
+ * \param volVars The volume variables in the sub control volumes
+ * \return the intrinsic permeability
+ */
+ Scalar intrinsicPermeability (const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ const Scalar radius = this->radius(scv);
+ const Scalar gamma = 2; // quadratic velocity profile (Poiseuille flow)
+ return radius*radius/(2*(2+gamma));
+ }
+
+ /*!
+ * \brief Returns the porosity \f$[-]\f$
+ *
+ * \param scv The sub control volume
+ * \return the porosity
+ */
+ Scalar porosity(const SubControlVolume &scv) const
+ { return 1.0; }
+
+private:
+ Scalar radius_, radiusMain_;
+ static constexpr Scalar eps_ = 1e-8;
+};
+
+} // end namespace Dumux
+
+#endif
--
GitLab