From 0fd50e1c27e73fac059ed99a7cec9ade9bd2de90 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 14 Jul 2010 08:57:38 +0000
Subject: [PATCH] test folder completed
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@3832 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
CMake/Modules/DumuxMacros.cmake | 183 ++++++++
CMake/Modules/FindALUGrid.cmake | 26 ++
CMake/Modules/FindAlberta.cmake | 20 +
CMake/Modules/FindDUNE_common.cmake | 18 +
CMake/Modules/FindDUNE_grid.cmake | 16 +
CMake/Modules/FindDUNE_istl.cmake | 17 +
CMake/Modules/FindDUNE_localfunctions.cmake | 17 +
CMake/Modules/FindDUNE_mux.cmake | 25 ++
CMake/Modules/FindDUNE_pdelab.cmake | 17 +
CMake/Modules/FindMETIS.cmake | 24 ++
CMake/Modules/FindMPI.cmake | 319 ++++++++++++++
CMake/Modules/FindUG.cmake | 60 +++
configure.ac | 126 +++++-
dumux/Makefile.am | 3 -
dumux/dumux/Makefile.am | 4 -
dumux/dumux/dumux.hh | 6 -
test/CMakeLists.txt | 7 +
test/Makefile.am | 4 +
test/boxmodels/1p/1pproblem.hh | 308 ++++++++++++++
test/boxmodels/1p/1pspatialparameters.hh | 75 ++++
test/boxmodels/1p/CMakeLists.txt | 18 +
test/boxmodels/1p/Makefile.am | 7 +
test/boxmodels/1p/test_1p.cc | 201 +++++++++
test/boxmodels/1p2c/CMakeLists.txt | 17 +
test/boxmodels/1p2c/Makefile.am | 5 +
test/boxmodels/1p2c/test_1p2c.cc | 95 +++++
test/boxmodels/1p2c/tissue_soilproperties.hh | 121 ++++++
test/boxmodels/1p2c/tissue_tumor_problem.hh | 319 ++++++++++++++
test/boxmodels/2p/CMakeLists.txt | 18 +
test/boxmodels/2p/Makefile.am | 5 +
test/boxmodels/2p/lensproblem.hh | 394 ++++++++++++++++++
test/boxmodels/2p/lensspatialparameters.hh | 164 ++++++++
test/boxmodels/2p/test_2p.cc | 224 ++++++++++
test/boxmodels/2p2c/CMakeLists.txt | 18 +
test/boxmodels/2p2c/Makefile.am | 5 +
test/boxmodels/2p2c/injectionproblem.hh | 321 ++++++++++++++
.../2p2c/injectionspatialparameters.hh | 243 +++++++++++
test/boxmodels/2p2c/test_2p2c.cc | 24 ++
test/boxmodels/2p2cni/CMakeLists.txt | 18 +
test/boxmodels/2p2cni/Makefile.am | 5 +
test/boxmodels/2p2cni/test_2p2cni.cc | 75 ++++
test/boxmodels/2p2cni/waterairproblem.hh | 353 ++++++++++++++++
.../2p2cni/waterairspatialparameters.hh | 247 +++++++++++
test/boxmodels/2pni/CMakeLists.txt | 17 +
test/boxmodels/2pni/Makefile.am | 5 +
test/boxmodels/2pni/injectionproblem2pni.hh | 362 ++++++++++++++++
test/boxmodels/2pni/test_2pni.cc | 26 ++
test/boxmodels/CMakeLists.txt | 11 +
test/boxmodels/Makefile.am | 4 +
test/boxmodels/richards/CMakeLists.txt | 23 +
test/boxmodels/richards/Makefile.am | 6 +
test/boxmodels/richards/richardsproblem.hh | 293 +++++++++++++
test/boxmodels/richards/richardssoil.hh | 94 +++++
test/boxmodels/richards/test_richards.cc | 93 +++++
test/common/CMakeLists.txt | 6 +
test/common/Makefile.am | 4 +
test/common/pardiso/CMakeLists.txt | 18 +
test/common/pardiso/Makefile.am | 18 +
test/common/pardiso/test_pardiso.cc | 198 +++++++++
test/common/propertysystem/CMakeLists.txt | 5 +
test/common/propertysystem/Makefile.am | 12 +
.../propertysystem/test_propertysystem.cc | 183 ++++++++
test/common/spline/Makefile.am | 5 +
test/common/spline/test_spline.cc | 69 +++
test/decoupled/1p/Makefile.am | 40 ++
test/decoupled/1p/test_diffusion.cc | 128 ++++++
test/decoupled/1p/test_diffusion_problem.hh | 294 +++++++++++++
.../1p/test_diffusion_spatialparams.hh | 100 +++++
test/decoupled/2p/Makefile.am | 7 +
test/decoupled/2p/test_2p.cc | 113 +++++
test/decoupled/2p/test_2p_problem.hh | 282 +++++++++++++
test/decoupled/2p/test_2p_spatialparams.hh | 108 +++++
test/decoupled/Makefile.am | 7 +
73 files changed, 6679 insertions(+), 24 deletions(-)
create mode 100644 CMake/Modules/DumuxMacros.cmake
create mode 100644 CMake/Modules/FindALUGrid.cmake
create mode 100644 CMake/Modules/FindAlberta.cmake
create mode 100644 CMake/Modules/FindDUNE_common.cmake
create mode 100644 CMake/Modules/FindDUNE_grid.cmake
create mode 100644 CMake/Modules/FindDUNE_istl.cmake
create mode 100644 CMake/Modules/FindDUNE_localfunctions.cmake
create mode 100644 CMake/Modules/FindDUNE_mux.cmake
create mode 100644 CMake/Modules/FindDUNE_pdelab.cmake
create mode 100644 CMake/Modules/FindMETIS.cmake
create mode 100644 CMake/Modules/FindMPI.cmake
create mode 100644 CMake/Modules/FindUG.cmake
delete mode 100644 dumux/Makefile.am
delete mode 100644 dumux/dumux/Makefile.am
delete mode 100644 dumux/dumux/dumux.hh
create mode 100644 test/CMakeLists.txt
create mode 100644 test/Makefile.am
create mode 100644 test/boxmodels/1p/1pproblem.hh
create mode 100644 test/boxmodels/1p/1pspatialparameters.hh
create mode 100644 test/boxmodels/1p/CMakeLists.txt
create mode 100644 test/boxmodels/1p/Makefile.am
create mode 100644 test/boxmodels/1p/test_1p.cc
create mode 100644 test/boxmodels/1p2c/CMakeLists.txt
create mode 100644 test/boxmodels/1p2c/Makefile.am
create mode 100644 test/boxmodels/1p2c/test_1p2c.cc
create mode 100644 test/boxmodels/1p2c/tissue_soilproperties.hh
create mode 100644 test/boxmodels/1p2c/tissue_tumor_problem.hh
create mode 100644 test/boxmodels/2p/CMakeLists.txt
create mode 100644 test/boxmodels/2p/Makefile.am
create mode 100644 test/boxmodels/2p/lensproblem.hh
create mode 100644 test/boxmodels/2p/lensspatialparameters.hh
create mode 100644 test/boxmodels/2p/test_2p.cc
create mode 100644 test/boxmodels/2p2c/CMakeLists.txt
create mode 100644 test/boxmodels/2p2c/Makefile.am
create mode 100644 test/boxmodels/2p2c/injectionproblem.hh
create mode 100644 test/boxmodels/2p2c/injectionspatialparameters.hh
create mode 100644 test/boxmodels/2p2c/test_2p2c.cc
create mode 100644 test/boxmodels/2p2cni/CMakeLists.txt
create mode 100644 test/boxmodels/2p2cni/Makefile.am
create mode 100644 test/boxmodels/2p2cni/test_2p2cni.cc
create mode 100644 test/boxmodels/2p2cni/waterairproblem.hh
create mode 100644 test/boxmodels/2p2cni/waterairspatialparameters.hh
create mode 100644 test/boxmodels/2pni/CMakeLists.txt
create mode 100644 test/boxmodels/2pni/Makefile.am
create mode 100644 test/boxmodels/2pni/injectionproblem2pni.hh
create mode 100644 test/boxmodels/2pni/test_2pni.cc
create mode 100644 test/boxmodels/CMakeLists.txt
create mode 100644 test/boxmodels/Makefile.am
create mode 100644 test/boxmodels/richards/CMakeLists.txt
create mode 100644 test/boxmodels/richards/Makefile.am
create mode 100644 test/boxmodels/richards/richardsproblem.hh
create mode 100644 test/boxmodels/richards/richardssoil.hh
create mode 100644 test/boxmodels/richards/test_richards.cc
create mode 100644 test/common/CMakeLists.txt
create mode 100644 test/common/Makefile.am
create mode 100644 test/common/pardiso/CMakeLists.txt
create mode 100644 test/common/pardiso/Makefile.am
create mode 100644 test/common/pardiso/test_pardiso.cc
create mode 100644 test/common/propertysystem/CMakeLists.txt
create mode 100644 test/common/propertysystem/Makefile.am
create mode 100644 test/common/propertysystem/test_propertysystem.cc
create mode 100644 test/common/spline/Makefile.am
create mode 100644 test/common/spline/test_spline.cc
create mode 100644 test/decoupled/1p/Makefile.am
create mode 100644 test/decoupled/1p/test_diffusion.cc
create mode 100644 test/decoupled/1p/test_diffusion_problem.hh
create mode 100644 test/decoupled/1p/test_diffusion_spatialparams.hh
create mode 100644 test/decoupled/2p/Makefile.am
create mode 100644 test/decoupled/2p/test_2p.cc
create mode 100644 test/decoupled/2p/test_2p_problem.hh
create mode 100644 test/decoupled/2p/test_2p_spatialparams.hh
create mode 100644 test/decoupled/Makefile.am
diff --git a/CMake/Modules/DumuxMacros.cmake b/CMake/Modules/DumuxMacros.cmake
new file mode 100644
index 0000000000..589878cb22
--- /dev/null
+++ b/CMake/Modules/DumuxMacros.cmake
@@ -0,0 +1,183 @@
+#############################################################
+# This sets up the DumuxMacros for the current CMake module.
+# Call DumuxSweep at the end of your CMake module in order to
+# not pullute the namespace with unused variables
+#############################################################
+macro(DumuxSetup
+ CMakeModuleName
+ ModuleName
+ Framework)
+ #############
+ # Set some internal variables which are used within
+ # the current CMake Module
+ set(DumuxModule ${CMakeModuleName})
+ set(DumuxModuleName ${ModuleName})
+ set(DumuxFramework ${Framework})
+
+ set(DumuxLibsFound 1)
+ set(DumuxLibraryNames)
+ set(DumuxFound 0)
+
+ set(DumuxPathMessage
+"Set the ${DumuxModule}_DIR cmake cache entry to the directory
+where the ${DumuxModuleName} libraries reside. Alternatively you can set
+the ${DumuxFramework}_DIR entry where all ${DumuxFramework} sub-modules have been compiled.")
+
+ # Base path to look for libraries and includes
+ if(${DumuxModule}_DIR)
+ list(APPEND DumuxModulePath ${${DumuxModule}_DIR})
+ endif(${DumuxModule}_DIR)
+ if(${DumuxFramework}_DIR)
+ list(APPEND DumuxModulePath "${${DumuxFramework}_DIR}/${DumuxModuleName}")
+ endif(${DumuxFramework}_DIR)
+
+ # Path to look for includes (->DumuxIncludePath) and libraries (-> DumuxLibraryPath)
+ foreach(tmp ${DumuxModulePath})
+ list(APPEND DumuxIncludePath "${tmp}" "${tmp}/include")
+ list(APPEND DumuxLibraryPath "${tmp}" "${tmp}/lib")
+ endforeach(tmp)
+ list(APPEND DumuxIncludePath "/usr/include" "/usr/local/include")
+ list(APPEND DumuxLibraryPath "/usr/lib" "/usr/local/lib")
+
+ set(DumuxLibraries)
+ set(DumuxFailedLibraries)
+endmacro(DumuxSetup)
+
+#############################################################
+# This adds some additional paths to the location where
+# includes and libraries are searched
+#############################################################
+macro(DumuxAddPathSuffixes
+ IncludeSuffixes
+ LibSuffixes)
+ foreach(tmp ${DumuxModulePath})
+ # deal with the user defined library locations
+ foreach(foo ${LibSuffixes})
+ list(APPEND DumuxLibraryPath "${tmp}/${foo}")
+ endforeach(foo)
+
+ # deal with the user defined include locations
+ foreach(foo ${IncludeSuffixes})
+ list(APPEND DumuxIncludePath "${tmp}/${foo}")
+ endforeach(foo)
+ endforeach(tmp)
+endmacro(DumuxAddPathSuffixes)
+#############################################################
+# Find a given library using some reasonable default
+# search paths. Sets Dumux${LibName}_LIBRARY to the location
+# where the library was found and extends the DumuxLibraries
+# variable.
+#############################################################
+macro(DumuxFindLibrary LibName)
+ set(Lib ${DumuxModule}_${LibName}_LIBRARY)
+
+ find_library(${Lib}
+ ${LibName}
+ PATHS ${DumuxLibraryPath}
+ PATH_SUFFIXES ".libs")
+
+ if(${Lib})
+ list(APPEND DumuxLibraries ${${Lib}})
+ list(APPEND DumuxLibraryNames ${LibName})
+ else(${Lib})
+ list(APPEND DumuxFailedLibraries ${LibName})
+ endif(${Lib})
+endmacro(DumuxFindLibrary)
+
+#############################################################
+# Find a given header file using some reasonable default
+# search paths.
+#############################################################
+macro(DumuxFindExtraIncludeDir VarName HeaderName)
+ set(Inc ${DumuxModule}_${VarName}_INCLUDE_DIR)
+ find_path(${Inc}
+ ${HeaderName}
+ PATHS ${DumuxIncludePath})
+ if(${Inc})
+ list(APPEND ${DumuxModule}_INCLUDE_DIRS ${${Inc}})
+ list(APPEND DumuxIncludes ${Inc})
+ else(${Inc})
+ list(APPEND DumuxFailedIncludes ${HeaderName})
+ endif(${Inc})
+endmacro(DumuxFindExtraIncludeDir)
+
+macro(DumuxFindIncludeDir HeaderName)
+ set(Inc ${DumuxModule}_INCLUDE_DIR)
+ find_path(${Inc}
+ ${HeaderName}
+ PATHS ${DumuxIncludePath})
+ if(${Inc})
+ list(APPEND ${DumuxModule}_INCLUDE_DIRS "${${Inc}}")
+ list(APPEND DumuxIncludes ${Inc})
+ else(${Inc})
+ list(APPEND DumuxFailedIncludes ${HeaderName})
+ endif(${Inc})
+endmacro(DumuxFindIncludeDir)
+
+macro(DumuxFindIncludeBaseDir HeaderName DirSuffix)
+ set(Inc ${DumuxModule}_INCLUDE_DIR)
+ find_path(${Inc}
+ ${HeaderName}
+ PATHS ${DumuxIncludePath})
+ if(${Inc})
+ list(APPEND ${DumuxModule}_INCLUDE_DIRS "${${Inc}}/${DirSuffix}")
+ list(APPEND DumuxIncludes ${Inc})
+ else(${Inc})
+ list(APPEND DumuxFailedIncludes ${HeaderName})
+ endif(${Inc})
+endmacro(DumuxFindIncludeBaseDir)
+
+#############################################################
+# Make sure the required libraries were found
+#############################################################
+macro(DumuxRequiredLibsFound)
+ set(DumuxLibsFound 1)
+ set(DumuxFailedLibsMessage "Could not find the required libraries ")
+ foreach(curLib ${ARGN})
+ set(curLibFound 0)
+ foreach(tmp ${DumuxLibraryNames})
+ if (tmp STREQUAL ${curLib})
+ set(curLibFound 1)
+ endif (tmp STREQUAL ${curLib})
+ endforeach(tmp)
+
+ if (NOT curLibFound)
+ set(DumuxLibsFound 0)
+ set(DumuxFailedLibsMessage "${DumuxFailedLibsMessage} '${curLib}'")
+ endif(NOT curLibFound)
+ endforeach(curLib)
+endmacro(DumuxRequiredLibsFound)
+
+#############################################################
+# Make sure the required libraries were found
+#############################################################
+macro(DumuxIncludeDirsFound)
+endmacro(DumuxIncludeDirsFound)
+
+#############################################################
+# Make sure everything required was found
+#############################################################
+macro(DumuxCheckFound)
+ # Set the global macros
+ set(DumuxFound 0)
+
+ if(DumuxLibsFound AND ${DumuxModule}_INCLUDE_DIR)
+ set(DumuxFound 1)
+ endif(DumuxLibsFound AND ${DumuxModule}_INCLUDE_DIR)
+ set(${DumuxModule}_FOUND ${DumuxFound})
+ set(${DumuxModule}_LIBRARIES ${DumuxLibraries})
+
+ # print status message if requested
+ if(NOT ${DumuxModule}_FIND_QUIETLY AND DumuxFound)
+ message(STATUS "Found ${DumuxModule}")
+ endif(NOT ${DumuxModule}_FIND_QUIETLY AND DumuxFound)
+
+ if(NOT DumuxFound AND ${DumuxModule}_FIND_REQUIRED)
+ if (DumuxLibsFound)
+ message(FATAL_ERROR "${DumuxPathMessage}")
+ else (DumuxLibsFound)
+ message(FATAL_ERROR "${DumuxPathMessage}
+${DumuxFailedLibsMessage}")
+ endif( DumuxLibsFound)
+ endif(NOT DumuxFound AND ${DumuxModule}_FIND_REQUIRED)
+endmacro(DumuxCheckFound)
diff --git a/CMake/Modules/FindALUGrid.cmake b/CMake/Modules/FindALUGrid.cmake
new file mode 100644
index 0000000000..efa2ad7b2f
--- /dev/null
+++ b/CMake/Modules/FindALUGrid.cmake
@@ -0,0 +1,26 @@
+# -*-cmake-*-
+# - Try to find the UG grid manager
+# Once done this will define:
+# ALUGrid_FOUND - system has dune-grid
+# UG_INCLUDE_DIR - incude paths to use dune-grid
+# UG_LIBRARIES - Link these to use dune-grid
+Include(DumuxMacros)
+
+DumuxSetup("ALUGrid" "ALUGrid" "ALUGrid")
+
+set(MyIncludeSuffixes
+ "include/serial"
+ "include/parallel"
+ "include/duneinterface")
+
+DumuxAddPathSuffixes("${MyIncludeSuffixes}" "")
+
+DumuxFindIncludeDir("alugrid_2d.h")
+DumuxFindExtraIncludeDir("ALU_SERIAL" "serialize.h")
+DumuxFindExtraIncludeDir("ALU_PARALLEL" "gitter_pll_impl.h")
+DumuxFindExtraIncludeDir("ALU_DUNE" "gitter_dune_impl.h")
+DumuxFindLibrary("alugrid")
+
+DumuxRequiredLibsFound("alugrid")
+DumuxIncludeDirsFound()
+DumuxCheckFound()
diff --git a/CMake/Modules/FindAlberta.cmake b/CMake/Modules/FindAlberta.cmake
new file mode 100644
index 0000000000..c4c779fa52
--- /dev/null
+++ b/CMake/Modules/FindAlberta.cmake
@@ -0,0 +1,20 @@
+# -*-cmake-*-
+# - Try to find the Alberta grid manager
+# Once done this will define:
+# Alberta_FOUND - system has dune-grid
+# Alberta_INCLUDE_DIR - incude paths to use dune-grid
+# Alberta_LIBRARIES - Link these to use dune-grid
+Include(DumuxMacros)
+
+DumuxSetup("Alberta" "Alberta" "Alberta")
+
+#DumuxAddPathSuffixes("${MyIncludeSuffixes}" "")
+
+DumuxFindIncludeDir("alberta.h")
+DumuxFindLibrary("ALBERTA22_0")
+DumuxFindLibrary("ALBERTA22_1")
+DumuxFindLibrary("alberta_util")
+
+DumuxRequiredLibsFound("ALBERTA22_0" "ALBERTA22_1" "alberta_util")
+DumuxIncludeDirsFound()
+DumuxCheckFound()
diff --git a/CMake/Modules/FindDUNE_common.cmake b/CMake/Modules/FindDUNE_common.cmake
new file mode 100644
index 0000000000..15c56705db
--- /dev/null
+++ b/CMake/Modules/FindDUNE_common.cmake
@@ -0,0 +1,18 @@
+# -*-cmake-*-
+# - Try to find tje DUNE common library
+# Once done this will define:
+# DUNE_common_FOUND - system has dune-common
+# DUNE_common_INCLUDE_DIR - incude paths to use dune-common
+# DUNE_common_LIBRARIES - Link these to use dune-common
+INCLUDE(DumuxMacros)
+
+DumuxSetup("DUNE_common" "dune-common" "DUNE")
+
+DumuxFindIncludeDir("dune/common/misc.hh")
+DumuxFindLibrary("dunecommon")
+
+DumuxRequiredLibsFound("dunecommon")
+DumuxIncludeDirsFound()
+DumuxCheckFound()
+
+
diff --git a/CMake/Modules/FindDUNE_grid.cmake b/CMake/Modules/FindDUNE_grid.cmake
new file mode 100644
index 0000000000..475cc544bc
--- /dev/null
+++ b/CMake/Modules/FindDUNE_grid.cmake
@@ -0,0 +1,16 @@
+# -*-cmake-*-
+# - Try to find tje DUNE grid library
+# Once done this will define:
+# Dune_grid_FOUND - system has dune-grid
+# Dune_grid_INCLUDE_DIR - incude paths to use dune-grid
+# Dune_grid_LIBRARIES - Link these to use dune-grid
+INCLUDE(DumuxMacros)
+
+DumuxSetup("DUNE_grid" "dune-grid" "DUNE")
+
+DumuxFindIncludeDir("dune/grid/sgrid.hh")
+DumuxFindLibrary("dunegrid")
+
+DumuxRequiredLibsFound("dunegrid")
+DumuxIncludeDirsFound()
+DumuxCheckFound()
diff --git a/CMake/Modules/FindDUNE_istl.cmake b/CMake/Modules/FindDUNE_istl.cmake
new file mode 100644
index 0000000000..32ea2a28b6
--- /dev/null
+++ b/CMake/Modules/FindDUNE_istl.cmake
@@ -0,0 +1,17 @@
+# -*-cmake-*-
+# - Try to find tje DUNE istl library
+# Once done this will define:
+# DUNE_istl_FOUND - system has dune-istl
+# DUNE_istl_INCLUDE_DIR - incude paths to use dune-istl
+# DUNE_istl_LIBRARIES - Link these to use dune-istl
+INCLUDE(DumuxMacros)
+
+DumuxSetup("DUNE_istl" "dune-istl" "DUNE")
+
+DumuxFindIncludeDir("dune/istl/io.hh")
+
+DumuxRequiredLibsFound()
+DumuxIncludeDirsFound()
+DumuxCheckFound()
+
+
diff --git a/CMake/Modules/FindDUNE_localfunctions.cmake b/CMake/Modules/FindDUNE_localfunctions.cmake
new file mode 100644
index 0000000000..05227725b4
--- /dev/null
+++ b/CMake/Modules/FindDUNE_localfunctions.cmake
@@ -0,0 +1,17 @@
+# -*-cmake-*-
+# - Try to find the DUNE localfunctions library
+# Once done this will define:
+# DUNE_localfunctions_FOUND - system has dune-localfunctions
+# DUNE_localfunctions_INCLUDE_DIR - incude paths to use dune-localfunctions
+# DUNE_localfunctions_LIBRARIES - Link these to use dune-localfunctions
+INCLUDE(DumuxMacros)
+
+DumuxSetup("DUNE_localfunctions" "dune-localfunctions" "DUNE")
+
+DumuxFindIncludeDir("dune/localfunctions/lagrange.hh")
+
+DumuxRequiredLibsFound()
+DumuxIncludeDirsFound()
+DumuxCheckFound()
+
+
diff --git a/CMake/Modules/FindDUNE_mux.cmake b/CMake/Modules/FindDUNE_mux.cmake
new file mode 100644
index 0000000000..1134ec0bb2
--- /dev/null
+++ b/CMake/Modules/FindDUNE_mux.cmake
@@ -0,0 +1,25 @@
+# -*-cmake-*-
+# Try to find the DUMUX library
+# Once done this will define:
+# DUNE_mux_FOUND - system has dune-mux
+# DUNE_mux_INCLUDE_DIR - incude paths to use dune-mux
+# DUNE_mux_LIBRARIES - Link these to use dune-mux
+INCLUDE(DumuxMacros)
+
+DumuxSetup("DUNE_mux" "dune-mux" "DUNE")
+
+set(MyLibSuffixes
+ "dumux"
+ "dumux/shapefunctions"
+ "dumux/onedinndgrid")
+DumuxAddPathSuffixes("" "${MyLibSuffixes}" )
+
+#DumuxFindIncludeBaseDir("dumux/material/twophaserelations.hh" "..")
+DumuxFindIncludeDir("dumux/material/twophaserelations.hh")
+
+DumuxFindLibrary("dumux")
+DumuxRequiredLibsFound("dumux")
+DumuxIncludeDirsFound()
+DumuxCheckFound()
+
+
diff --git a/CMake/Modules/FindDUNE_pdelab.cmake b/CMake/Modules/FindDUNE_pdelab.cmake
new file mode 100644
index 0000000000..e66cb021f4
--- /dev/null
+++ b/CMake/Modules/FindDUNE_pdelab.cmake
@@ -0,0 +1,17 @@
+# -*-cmake-*-
+# - Try to find the DUNE pdelab library
+# Once done this will define:
+# DUNE_pdelab_FOUND - system has dune-pdelab
+# DUNE_pdelab_INCLUDE_DIR - incude paths to use dune-pdelab
+# DUNE_pdelab_LIBRARIES - Link these to use dune-pdelab
+INCLUDE(DumuxMacros)
+
+DumuxSetup("DUNE_pdelab" "dune-pdelab" "DUNE")
+
+DumuxFindIncludeDir("dune/pdelab/common/function.hh")
+
+DumuxRequiredLibsFound()
+DumuxIncludeDirsFound()
+DumuxCheckFound()
+
+
diff --git a/CMake/Modules/FindMETIS.cmake b/CMake/Modules/FindMETIS.cmake
new file mode 100644
index 0000000000..bf4db205f6
--- /dev/null
+++ b/CMake/Modules/FindMETIS.cmake
@@ -0,0 +1,24 @@
+# -*-cmake-*-
+# - Try to find the libMETIS graph partioning library
+# Once done this will define:
+# METIS_FOUND - system has the libMETIS graph partioning library
+# METIS_INCLUDE_DIR - incude paths to use libMETIS
+# METIS_LIBRARIES - Link these to use libMETIS
+Include(DumuxMacros)
+
+DumuxSetup("METIS" "METIS" "METIS")
+
+set(MyIncludeSuffixes "METISLib")
+#set(MyLibSuffixes "build/Linux-x86_64/" "build/Linux-i686/" "GKlib/builds/Linux-x86_64/" "GKlib/builds/Linux-i686/")
+
+#DumuxAddPathSuffixes("${MyIncludeSuffixes}" "${MyLibSuffixes}")
+DumuxAddPathSuffixes("${MyIncludeSuffixes}" "")
+
+DumuxFindIncludeDir("metis.h")
+DumuxFindLibrary("metis")
+#DumuxFindLibrary("GKlib")
+
+#DumuxRequiredLibsFound("metis" "GKlib")
+DumuxRequiredLibsFound("metis")
+DumuxIncludeDirsFound()
+DumuxCheckFound()
diff --git a/CMake/Modules/FindMPI.cmake b/CMake/Modules/FindMPI.cmake
new file mode 100644
index 0000000000..64d8f3410c
--- /dev/null
+++ b/CMake/Modules/FindMPI.cmake
@@ -0,0 +1,319 @@
+# - Message Passing Interface (MPI) module.
+#
+# The Message Passing Interface (MPI) is a library used to write
+# high-performance parallel applications that use message passing, and
+# is typically deployed on a cluster. MPI is a standard interface
+# (defined by the MPI forum) for which many implementations are
+# available. All of these implementations have somewhat different
+# compilation approaches (different include paths, libraries to link
+# against, etc.), and this module tries to smooth out those differences.
+#
+# This module will set the following variables:
+# MPI_FOUND TRUE if we have found MPI
+# MPI_COMPILE_FLAGS Compilation flags for MPI programs
+# MPI_INCLUDE_PATH Include path(s) for MPI header
+# MPI_LINK_FLAGS Linking flags for MPI programs
+# MPI_LIBRARY First MPI library to link against (cached)
+# MPI_EXTRA_LIBRARY Extra MPI libraries to link against (cached)
+# MPI_LIBRARIES All libraries to link MPI programs against
+# MPIEXEC Executable for running MPI programs
+# MPIEXEC_NUMPROC_FLAG Flag to pass to MPIEXEC before giving it the
+# number of processors to run on
+# MPIEXEC_PREFLAGS Flags to pass to MPIEXEC directly before the
+# executable to run.
+# MPIEXEC_POSTFLAGS Flags to pass to MPIEXEC after all other flags.
+#
+# This module will attempt to auto-detect these settings, first by
+# looking for a MPI compiler, which many MPI implementations provide
+# as a pass-through to the native compiler to simplify the compilation
+# of MPI programs. The MPI compiler is stored in the cache variable
+# MPI_COMPILER, and will attempt to look for commonly-named drivers
+# mpic++, mpicxx, mpiCC, or mpicc. If the compiler driver is found and
+# recognized, it will be used to set all of the module variables. To
+# skip this auto-detection, set MPI_LIBRARY and MPI_INCLUDE_PATH in
+# the CMake cache.
+#
+# If no compiler driver is found or the compiler driver is not
+# recognized, this module will then search for common include paths
+# and library names to try to detect MPI.
+#
+# If CMake initially finds a different MPI than was intended, and you
+# want to use the MPI compiler auto-detection for a different MPI
+# implementation, set MPI_COMPILER to the MPI compiler driver you want
+# to use (e.g., mpicxx) and then set MPI_LIBRARY to the string
+# MPI_LIBRARY-NOTFOUND. When you re-configure, auto-detection of MPI
+# will run again with the newly-specified MPI_COMPILER.
+#
+# When using MPIEXEC to execute MPI applications, you should typically
+# use all of the MPIEXEC flags as follows:
+# ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} PROCS ${MPIEXEC_PREFLAGS} EXECUTABLE
+# ${MPIEXEC_POSTFLAGS} ARGS
+# where PROCS is the number of processors on which to execute the program,
+# EXECUTABLE is the MPI program, and ARGS are the arguments to pass to the
+# MPI program.
+
+# Try to find the MPI driver program
+find_program(MPI_COMPILER
+ NAMES mpicxx mpiCC mpicc mpic++
+ DOC "MPI compiler. Used only to detect MPI compilation flags.")
+mark_as_advanced(MPI_COMPILER)
+
+find_program(MPIEXEC
+ NAMES mpiexec mpirun lamexec
+ DOC "Executable for running MPI programs.")
+
+set(MPIEXEC_NUMPROC_FLAG "-np" CACHE STRING "Flag used by MPI to specify the number of processes for MPIEXEC; the next option will be the number of processes.")
+set(MPIEXEC_PREFLAGS "" CACHE STRING "These flags will be directly before the executable that is being run by MPIEXEC.")
+set(MPIEXEC_POSTFLAGS "" CACHE STRING "These flags will come after all flags given to MPIEXEC.")
+set(MPIEXEC_MAX_NUMPROCS "2" CACHE STRING "Maximum number of processors available to run MPI applications.")
+mark_as_advanced(MPIEXEC MPIEXEC_NUMPROC_FLAG MPIEXEC_PREFLAGS
+ MPIEXEC_POSTFLAGS MPIEXEC_MAX_NUMPROCS)
+
+if (MPI_INCLUDE_PATH AND MPI_LIBRARY)
+ # Do nothing: we already have MPI_INCLUDE_PATH and MPI_LIBRARY in
+ # the cache, and we don't want to override those settings.
+elseif (MPI_COMPILER)
+ # Check whether the -showme:compile option works. This indicates
+ # that we have either Open MPI or a newer version of LAM-MPI, and
+ # implies that -showme:link will also work.
+ exec_program(${MPI_COMPILER}
+ ARGS -showme:compile
+ OUTPUT_VARIABLE MPI_COMPILE_CMDLINE
+ RETURN_VALUE MPI_COMPILER_RETURN)
+
+ if (MPI_COMPILER_RETURN EQUAL 0)
+ # If we appear to have -showme:compile, then we should also have
+ # -showme:link. Try it.
+ exec_program(${MPI_COMPILER}
+ ARGS -showme:link
+ OUTPUT_VARIABLE MPI_LINK_CMDLINE
+ RETURN_VALUE MPI_COMPILER_RETURN)
+
+ # Note that we probably have -showme:incdirs and -showme:libdirs
+ # as well.
+ set(MPI_COMPILER_MAY_HAVE_INCLIBDIRS TRUE)
+ endif (MPI_COMPILER_RETURN EQUAL 0)
+
+ if (MPI_COMPILER_RETURN EQUAL 0)
+ # Do nothing: we have our command lines now
+ else (MPI_COMPILER_RETURN EQUAL 0)
+ # Older versions of LAM-MPI have "-showme". Try it.
+ exec_program(${MPI_COMPILER}
+ ARGS -showme
+ OUTPUT_VARIABLE MPI_COMPILE_CMDLINE
+ RETURN_VALUE MPI_COMPILER_RETURN)
+ endif (MPI_COMPILER_RETURN EQUAL 0)
+
+ if (MPI_COMPILER_RETURN EQUAL 0)
+ # Do nothing: we have our command lines now
+ else (MPI_COMPILER_RETURN EQUAL 0)
+ # MPICH uses "-compile-info". Try it.
+ exec_program(${MPI_COMPILER}
+ ARGS -compile-info
+ OUTPUT_VARIABLE MPI_COMPILE_CMDLINE
+ RETURN_VALUE MPI_COMPILER_RETURN)
+ endif (MPI_COMPILER_RETURN EQUAL 0)
+
+ if (MPI_COMPILER_RETURN EQUAL 0)
+ # We have our command lines, but we might need to copy
+ # MPI_COMPILE_CMDLINE into MPI_LINK_CMDLINE, if the underlying
+ if (NOT MPI_LINK_CMDLINE)
+ # MPICH uses "-link-info". Try it.
+ exec_program(${MPI_COMPILER}
+ ARGS -link-info
+ OUTPUT_VARIABLE MPI_LINK_CMDLINE
+ RETURN_VALUE MPI_COMPILER_RETURN)
+
+ if (NOT MPI_COMPILER_RETURN EQUAL 0)
+ message(STATUS "Unable to determine MPI from MPI driver ${MPI_COMPILER}")
+ endif (NOT MPI_COMPILER_RETURN EQUAL 0)
+ endif (NOT MPI_LINK_CMDLINE)
+ else (MPI_COMPILER_RETURN EQUAL 0)
+ message(STATUS "Unable to determine MPI from MPI driver ${MPI_COMPILER}")
+ endif (MPI_COMPILER_RETURN EQUAL 0)
+endif (MPI_INCLUDE_PATH AND MPI_LIBRARY)
+
+if (MPI_INCLUDE_PATH AND MPI_LIBRARY)
+ # Do nothing: we already have MPI_INCLUDE_PATH and MPI_LIBRARY in
+ # the cache, and we don't want to override those settings.
+elseif (MPI_COMPILE_CMDLINE)
+ # Extract compile flags from the compile command line.
+ string(REGEX MATCHALL "-D([^\" ]+|\"[^\"]+\")" MPI_ALL_COMPILE_FLAGS "${MPI_COMPILE_CMDLINE}")
+ set(MPI_COMPILE_FLAGS_WORK)
+ foreach(FLAG ${MPI_ALL_COMPILE_FLAGS})
+ if (MPI_COMPILE_FLAGS_WORK)
+ set(MPI_COMPILE_FLAGS_WORK "${MPI_COMPILE_FLAGS_WORK} ${FLAG}")
+ else(MPI_COMPILE_FLAGS_WORK)
+ set(MPI_COMPILE_FLAGS_WORK ${FLAG})
+ endif(MPI_COMPILE_FLAGS_WORK)
+ endforeach(FLAG)
+
+ # Extract include paths from compile command line
+ string(REGEX MATCHALL "-I([^\" ]+|\"[^\"]+\")" MPI_ALL_INCLUDE_PATHS "${MPI_COMPILE_CMDLINE}")
+ set(MPI_INCLUDE_PATH_WORK)
+ foreach(IPATH ${MPI_ALL_INCLUDE_PATHS})
+ string(REGEX REPLACE "^-I" "" IPATH ${IPATH})
+ string(REGEX REPLACE "//" "/" IPATH ${IPATH})
+ list(APPEND MPI_INCLUDE_PATH_WORK ${IPATH})
+ endforeach(IPATH)
+
+ if (NOT MPI_INCLUDE_PATH_WORK)
+ if (MPI_COMPILER_MAY_HAVE_INCLIBDIRS)
+ # The compile command line didn't have any include paths on it,
+ # but we may have -showme:incdirs. Use it.
+ exec_program(${MPI_COMPILER}
+ ARGS -showme:incdirs
+ OUTPUT_VARIABLE MPI_INCLUDE_PATH_WORK
+ RETURN_VALUE MPI_COMPILER_RETURN)
+ separate_arguments(MPI_INCLUDE_PATH_WORK)
+ endif (MPI_COMPILER_MAY_HAVE_INCLIBDIRS)
+ endif (NOT MPI_INCLUDE_PATH_WORK)
+
+ if (NOT MPI_INCLUDE_PATH_WORK)
+ # If all else fails, just search for mpi.h in the normal include
+ # paths.
+ find_path(MPI_INCLUDE_PATH mpi.h)
+ set(MPI_INCLUDE_PATH_WORK ${MPI_INCLUDE_PATH})
+ endif (NOT MPI_INCLUDE_PATH_WORK)
+
+ # Extract linker paths from the link command line
+ string(REGEX MATCHALL "-L([^\" ]+|\"[^\"]+\")" MPI_ALL_LINK_PATHS "${MPI_LINK_CMDLINE}")
+ set(MPI_LINK_PATH)
+
+ foreach(LPATH ${MPI_ALL_LINK_PATHS})
+ string(REGEX REPLACE "^-L" "" LPATH ${LPATH})
+ string(REGEX REPLACE "//" "/" LPATH ${LPATH})
+ list(APPEND MPI_LINK_PATH ${LPATH})
+ endforeach(LPATH)
+
+ if (NOT MPI_LINK_PATH)
+ if (MPI_COMPILER_MAY_HAVE_INCLIBDIRS)
+ # The compile command line didn't have any linking paths on it,
+ # but we may have -showme:libdirs. Use it.
+ exec_program(${MPI_COMPILER}
+ ARGS -showme:libdirs
+ OUTPUT_VARIABLE MPI_LINK_PATH
+ RETURN_VALUE MPI_COMPILER_RETURN)
+ separate_arguments(MPI_LINK_PATH)
+ endif (MPI_COMPILER_MAY_HAVE_INCLIBDIRS)
+ endif (NOT MPI_LINK_PATH)
+
+ # Extract linker flags from the link command line
+ string(REGEX MATCHALL "-Wl,([^\" ]+|\"[^\"]+\")" MPI_ALL_LINK_FLAGS "${MPI_LINK_CMDLINE}")
+ set(MPI_LINK_FLAGS_WORK)
+ foreach(FLAG ${MPI_ALL_LINK_FLAGS})
+ if (MPI_LINK_FLAGS_WORK)
+ set(MPI_LINK_FLAGS_WORK "${MPI_LINK_FLAGS_WORK} ${FLAG}")
+ else(MPI_LINK_FLAGS_WORK)
+ set(MPI_LINK_FLAGS_WORK ${FLAG})
+ endif(MPI_LINK_FLAGS_WORK)
+ endforeach(FLAG)
+
+
+ # Extract the set of libraries to link against from the link command
+ # line.
+ string(REGEX MATCHALL "[ ]-l([^\" ]+|\"[^\"]+\")" MPI_LIBNAMES "${MPI_LINK_CMDLINE}")
+ string(REPLACE " " "" MPI_LIBNAMES "${MPI_LIBNAMES}")
+
+ # Determine full path names for all of the libraries that one needs
+ # to link against in an MPI program
+ set(MPI_LIBRARIES)
+ foreach(LIB ${MPI_LIBNAMES})
+ string(REGEX REPLACE "^-l" "" LIB ${LIB})
+ set(MPI_LIB "MPI_LIB-NOTFOUND" CACHE FILEPATH "Cleared" FORCE)
+ find_library(MPI_LIB ${LIB} HINTS ${MPI_LINK_PATH} PATHS "/usr/lib64" "/lib64")
+ if (MPI_LIB)
+ list(APPEND MPI_LIBRARIES ${MPI_LIB})
+ else (MPI_LIB)
+ message(SEND_ERROR "Unable to find MPI library ${LIB}")
+ endif (MPI_LIB)
+ endforeach(LIB)
+ set(MPI_LIB "MPI_LIB-NOTFOUND" CACHE INTERNAL "Scratch variable for MPI detection" FORCE)
+
+ # Chop MPI_LIBRARIES into the old-style MPI_LIBRARY and
+ # MPI_EXTRA_LIBRARY.
+ list(LENGTH MPI_LIBRARIES MPI_NUMLIBS)
+ list(LENGTH MPI_LIBNAMES MPI_NUMLIBS_EXPECTED)
+ if (MPI_NUMLIBS EQUAL MPI_NUMLIBS_EXPECTED)
+ list(GET MPI_LIBRARIES 0 MPI_LIBRARY_WORK)
+ set(MPI_LIBRARY ${MPI_LIBRARY_WORK} CACHE FILEPATH "MPI library to link against" FORCE)
+ else (MPI_NUMLIBS EQUAL MPI_NUMLIBS_EXPECTED)
+ set(MPI_LIBRARY "MPI_LIBRARY-NOTFOUND" CACHE FILEPATH "MPI library to link against" FORCE)
+ endif (MPI_NUMLIBS EQUAL MPI_NUMLIBS_EXPECTED)
+ if (MPI_NUMLIBS GREATER 1)
+ set(MPI_EXTRA_LIBRARY_WORK ${MPI_LIBRARIES})
+ list(REMOVE_AT MPI_EXTRA_LIBRARY_WORK 0)
+ set(MPI_EXTRA_LIBRARY ${MPI_EXTRA_LIBRARY_WORK} CACHE STRING "Extra MPI libraries to link against" FORCE)
+ else (MPI_NUMLIBS GREATER 1)
+ set(MPI_EXTRA_LIBRARY "MPI_EXTRA_LIBRARY-NOTFOUND" CACHE STRING "Extra MPI libraries to link against" FORCE)
+ endif (MPI_NUMLIBS GREATER 1)
+
+ # Set up all of the appropriate cache entries
+ set(MPI_COMPILE_FLAGS ${MPI_COMPILE_FLAGS_WORK} CACHE STRING "MPI compilation flags" FORCE)
+ set(MPI_INCLUDE_PATH ${MPI_INCLUDE_PATH_WORK} CACHE STRING "MPI include path" FORCE)
+ set(MPI_LINK_FLAGS ${MPI_LINK_FLAGS_WORK} CACHE STRING "MPI linking flags" FORCE)
+else (MPI_COMPILE_CMDLINE)
+ find_path(MPI_INCLUDE_PATH mpi.h
+ /usr/local/include
+ /usr/include
+ /usr/include/mpi
+ /usr/local/mpi/include
+ "C:/Program Files/MPICH/SDK/Include"
+ "$ENV{SystemDrive}/Program Files/MPICH2/include"
+ "$ENV{SystemDrive}/Program Files/Microsoft Compute Cluster Pack/Include"
+ )
+
+ # Decide between 32-bit and 64-bit libraries for Microsoft's MPI
+ if (CMAKE_CL_64)
+ set(MS_MPI_ARCH_DIR amd64)
+ else (CMAKE_CL_64)
+ set(MS_MPI_ARCH_DIR i386)
+ endif (CMAKE_CL_64)
+
+ find_library(MPI_LIBRARY
+ NAMES mpi mpich msmpi
+ PATHS /usr/lib /usr/local/lib /usr/local/mpi/lib
+ "C:/Program Files/MPICH/SDK/Lib"
+ "$ENV{SystemDrive}/Program Files/MPICH/SDK/Lib"
+ "$ENV{SystemDrive}/Program Files/Microsoft Compute Cluster Pack/Lib/${MS_MPI_ARCH_DIR}"
+ )
+ find_library(MPI_LIBRARY
+ NAMES mpich2
+ PATHS
+ "$ENV{SystemDrive}/Program Files/MPICH2/Lib")
+
+ find_library(MPI_EXTRA_LIBRARY
+ NAMES mpi++
+ PATHS /usr/lib /usr/local/lib /usr/local/mpi/lib
+ "C:/Program Files/MPICH/SDK/Lib"
+ DOC "Extra MPI libraries to link against.")
+
+ set(MPI_COMPILE_FLAGS "" CACHE STRING "MPI compilation flags")
+ set(MPI_LINK_FLAGS "" CACHE STRING "MPI linking flags")
+endif (MPI_INCLUDE_PATH AND MPI_LIBRARY)
+
+# on BlueGene/L the MPI lib is named libmpich.rts.a, there also these additional libs are required
+if("${MPI_LIBRARY}" MATCHES "mpich.rts")
+ set(MPI_EXTRA_LIBRARY ${MPI_EXTRA_LIBRARY} msglayer.rts devices.rts rts.rts devices.rts)
+ set(MPI_LIBRARY ${MPI_LIBRARY} msglayer.rts devices.rts rts.rts devices.rts)
+endif("${MPI_LIBRARY}" MATCHES "mpich.rts")
+
+# Set up extra variables to conform to
+if (MPI_EXTRA_LIBRARY)
+ set(MPI_LIBRARIES ${MPI_LIBRARY} ${MPI_EXTRA_LIBRARY})
+else (MPI_EXTRA_LIBRARY)
+ set(MPI_LIBRARIES ${MPI_LIBRARY})
+endif (MPI_EXTRA_LIBRARY)
+
+if (MPI_INCLUDE_PATH AND MPI_LIBRARY)
+ set(MPI_FOUND TRUE)
+else (MPI_INCLUDE_PATH AND MPI_LIBRARY)
+ set(MPI_FOUND FALSE)
+endif (MPI_INCLUDE_PATH AND MPI_LIBRARY)
+
+include(FindPackageHandleStandardArgs)
+# handle the QUIETLY and REQUIRED arguments
+find_package_handle_standard_args(MPI DEFAULT_MSG MPI_LIBRARY MPI_INCLUDE_PATH)
+
+mark_as_advanced(MPI_INCLUDE_PATH MPI_COMPILE_FLAGS MPI_LINK_FLAGS MPI_LIBRARY
+ MPI_EXTRA_LIBRARY)
diff --git a/CMake/Modules/FindUG.cmake b/CMake/Modules/FindUG.cmake
new file mode 100644
index 0000000000..e964a4a8b1
--- /dev/null
+++ b/CMake/Modules/FindUG.cmake
@@ -0,0 +1,60 @@
+# -*-cmake-*-
+# - Try to find the UG grid manager
+# Once done this will define:
+# UG_FOUND - system has dune-grid
+# UG_INCLUDE_DIR - incude paths to use dune-grid
+# UG_LIBRARIES - Link these to use dune-grid
+Include(DumuxMacros)
+
+
+DumuxSetup("UG" "UG" "UG")
+
+set(MyIncludeSuffixes
+ "gm"
+ "np")
+set(MyLibSuffixes
+ "np"
+ "np/field"
+ "np/udm"
+ "np/procs"
+ "np/amglib"
+ "np/algebra"
+ "dev"
+ "dev/rif"
+ "dev/sif"
+ "dev/ps"
+ "dev/meta"
+ "dev/xif"
+ "dev/ppm"
+ "dom/lgm/ngin2d"
+ "dom/lgm/ngin"
+ "dom/lgm"
+ "dom/std"
+ "parallel/dddif"
+ "parallel/util"
+ "graphics"
+ "graphics/uggraph"
+ "graphics/grape"
+ "low"
+ "gm"
+ "gm/gg2"
+ "gm/gg3"
+ "ui")
+set(MyUgLibs
+ "ugS2"
+ "ugS3"
+ "devS"
+)
+
+DumuxAddPathSuffixes("${MyIncludeSuffixes}" "${MyLibSuffixes}" )
+
+DumuxFindIncludeDir("ugm.h")
+DumuxFindExtraIncludeDir("GM_H" "gm.h")
+
+foreach(tmp ${MyUgLibs})
+ DumuxFindLibrary(${tmp})
+endforeach(tmp)
+
+DumuxRequiredLibsFound(${MyUgLibs})
+DumuxIncludeDirsFound()
+DumuxCheckFound()
diff --git a/configure.ac b/configure.ac
index 28fe85daf5..9e30a764ca 100644
--- a/configure.ac
+++ b/configure.ac
@@ -3,31 +3,135 @@
AC_PREREQ(2.50)
DUNE_AC_INIT # gets module version from dune.module file
AM_INIT_AUTOMAKE
-AC_CONFIG_SRCDIR([src/dumux.cc])
+AC_CONFIG_SRCDIR([test/decoupled/1p/test_diffusion.cc])
AM_CONFIG_HEADER([config.h])
+AC_CHECK_HEADER([valgrind/memcheck.h],
+ [HAVE_VALGRIND_H="1"],
+ AC_MSG_WARN([valgrind/memcheck.h not found]))
+if test x"$HAVE_VALGRIND_H" = "x1"; then
+ AC_DEFINE(HAVE_VALGRIND,1,[Define whether the valgrind header files for client requests are present.])
+fi
# we need no more than the standard DE-stuff
# this module depends on dune-pdelab
# this implies checking for [dune-common], [dune-grid], [dune-localfunctions], [dune-istl], [dune-pdelab]
DUNE_CHECK_ALL
+# the Boost c++ template libraries
+AX_BOOST_BASE([1.33.1])
+
+AC_DEFUN([SET_PARDISO], [
+AC_PREREQ(2.50)
+AC_REQUIRE([AC_F77_LIBRARY_LDFLAGS])
+acx_pardiso_ok=no
+
+AC_ARG_WITH(pardiso,
+ [AC_HELP_STRING([--with-pardiso=<lib>], [use PARDISO library <lib>])])
+case $with_pardiso in
+ yes | "") ;;
+ no) acx_pardiso_ok=disable ;;
+ -* | */* | *.a | *.so | *.so.* | *.o) PARDISO_LIBS="$with_pardiso" ;;
+ *) PARDISO_LIBS="-l$with_pardiso" ;;
+esac
+
+# Get fortran linker names of PARDISO functions to check for.
+AC_F77_FUNC(pardisoinit)
+
+acx_pardiso_save_LIBS="$LIBS"
+LIBS="$BLAS_LIBS $LIBS $FLIBS"
+
+# First, check PARDISO_LIBS environment variable
+if test $acx_pardiso_ok = no; then
+if test "x$PARDISO_LIBS" != x; then
+ save_LIBS="$LIBS"; LIBS="$PARDISO_LIBS $LIBS"
+ AC_MSG_CHECKING([for $pardisoinit in $PARDISO_LIBS])
+ AC_TRY_LINK_FUNC($pardisoinit, [acx_pardiso_ok=yes], [PARDISO_LIBS=""])
+ AC_MSG_RESULT($acx_pardiso_ok)
+ LIBS="$save_LIBS"
+fi
+fi
+
+AC_SUBST(PARDISO_LIBS)
+
+LIBS="$acx_pardiso_save_LIBS"
+
+# Finally, execute ACTION-IF-FOUND/ACTION-IF-NOT-FOUND:
+if test x"$acx_pardiso_ok" = xyes; then
+ ifelse([$1],,AC_DEFINE(HAVE_PARDISO,1,[Define if you have a PARDISO library.]),[$1])
+ :
+else
+ acx_pardiso_ok=no
+ $2
+fi
+])dnl SET_PARDISO
+
+SET_PARDISO
+
# implicitly set the Dune-flags everywhere
AC_SUBST(AM_CPPFLAGS, $DUNE_CPPFLAGS)
AC_SUBST(AM_LDFLAGS, $DUNE_LDFLAGS)
LIBS="$DUNE_LIBS"
-AC_CONFIG_FILES([
- Makefile
- src/Makefile
- doc/Makefile
- doc/doxygen/Makefile
- doc/doxygen/Doxyfile
- dune/Makefile
- dune/dumux/Makefile
- m4/Makefile
- dumux.pc
+AC_CONFIG_FILES([dumux.pc
+ Makefile
+ doc/Makefile
+ doc/doxygen/Makefile
+ doc/handbook/Makefile
+ dumux/Makefile
+ dumux/boxmodels/Makefile
+ dumux/boxmodels/1p/Makefile
+ dumux/boxmodels/1p2c/Makefile
+ dumux/boxmodels/2p/Makefile
+ dumux/boxmodels/2p2c/Makefile
+ dumux/boxmodels/2p2cni/Makefile
+ dumux/boxmodels/2pNc/Makefile
+ dumux/boxmodels/2pni/Makefile
+ dumux/boxmodels/boxscheme/Makefile
+ dumux/boxmodels/pdelab/Makefile
+ dumux/boxmodels/richards/Makefile
+ dumux/common/Makefile
+ dumux/decoupled/Makefile
+ dumux/decoupled/2p/Makefile
+ dumux/decoupled/2p/diffusion/Makefile
+ dumux/decoupled/2p/diffusion/fv/Makefile
+ dumux/decoupled/2p/diffusion/fvmpfa/Makefile
+ dumux/decoupled/2p/diffusion/mimetic/Makefile
+ dumux/decoupled/2p/impes/Makefile
+ dumux/decoupled/2p/transport/Makefile
+ dumux/decoupled/2p/transport/fv/Makefile
+ dumux/decoupled/common/Makefile
+ dumux/material/Makefile
+ dumux/material/binarycoefficients/Makefile
+ dumux/material/components/Makefile
+ dumux/material/fluidmatrixinteractions/Makefile
+ dumux/material/fluidsystems/Makefile
+ dumux/material/spatialparameters/Makefile
+ dumux/nonlinear/Makefile
+ lib/Makefile
+ m4/Makefile
+ test/Makefile
+ test/boxmodels/Makefile
+ test/boxmodels/1p/Makefile
+ test/boxmodels/2p/Makefile
+ test/boxmodels/2pni/Makefile
+ test/boxmodels/1p2c/Makefile
+ test/boxmodels/2p2c/Makefile
+ test/boxmodels/2p2cni/Makefile
+ test/boxmodels/2pNc/Makefile
+ test/boxmodels/richards/Makefile
+ test/common/Makefile
+ test/common/pardiso/Makefile
+ test/common/propertysystem/Makefile
+ test/common/spline/Makefile
+ test/decoupled/Makefile
+ test/decoupled/1p/Makefile
+ test/decoupled/2p/Makefile
+ tutorial/Makefile
])
+
AC_OUTPUT
# finally print the summary information
DUNE_SUMMARY_ALL
+echo "Pardiso..........: $acx_pardiso_ok"
+echo
diff --git a/dumux/Makefile.am b/dumux/Makefile.am
deleted file mode 100644
index 4e5a986abe..0000000000
--- a/dumux/Makefile.am
+++ /dev/null
@@ -1,3 +0,0 @@
-SUBDIRS = dumux
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/dumux/Makefile.am b/dumux/dumux/Makefile.am
deleted file mode 100644
index 17e94c15d0..0000000000
--- a/dumux/dumux/Makefile.am
+++ /dev/null
@@ -1,4 +0,0 @@
-dumuxincludedir = $(includedir)/dune/dumux
-dumuxinclude_HEADERS = dumux.hh
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/dumux/dumux.hh b/dumux/dumux/dumux.hh
deleted file mode 100644
index 945dc06eb5..0000000000
--- a/dumux/dumux/dumux.hh
+++ /dev/null
@@ -1,6 +0,0 @@
-#ifndef DUNE_dumux.hh
-#define DUNE_dumux.hh
-
-// add your classes here
-
-#endif // DUNE_dumux.hh
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
new file mode 100644
index 0000000000..10cc887605
--- /dev/null
+++ b/test/CMakeLists.txt
@@ -0,0 +1,7 @@
+# directories which are commented out do not compile at the moment
+# (they also don't compile if using the old build system)!
+
+add_subdirectory("boxmodels")
+add_subdirectory("common")
+#add_subdirectory("decoupled")
+
diff --git a/test/Makefile.am b/test/Makefile.am
new file mode 100644
index 0000000000..19434ebd63
--- /dev/null
+++ b/test/Makefile.am
@@ -0,0 +1,4 @@
+SUBDIRS = boxmodels common decoupled freeflow io misc modelcoupling multiscale onedinndgrid pardiso .
+
+include $(top_srcdir)/am/global-rules
+
diff --git a/test/boxmodels/1p/1pproblem.hh b/test/boxmodels/1p/1pproblem.hh
new file mode 100644
index 0000000000..c18d10436f
--- /dev/null
+++ b/test/boxmodels/1p/1pproblem.hh
@@ -0,0 +1,308 @@
+// $Id: 1pproblem.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Onur Dogan *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_1PPROBLEM_HH
+#define DUMUX_1PPROBLEM_HH
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+#include <dune/grid/io/file/gmshreader.hh>
+
+#include <dumux/boxmodels/1p/1pboxmodel.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+
+#include "1pspatialparameters.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class OnePTestProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(OnePTestProblem, INHERITS_FROM(BoxOneP));
+
+SET_PROP(OnePTestProblem, Fluid)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
+};
+
+// Set the grid type
+#if HAVE_UG
+SET_PROP(OnePTestProblem, Grid) { typedef Dune::UGGrid<3> type; };
+#else
+//SET_PROP(OnePTestProblem, Grid) { typedef Dune::SGrid<3, 3> type; };
+SET_TYPE_PROP(OnePTestProblem, Grid, Dune::YaspGrid<3>);
+#endif
+
+#ifdef HAVE_DUNE_PDELAB
+SET_PROP(OnePTestProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif
+
+SET_PROP(OnePTestProblem, NewtonLinearSolverVerbosity)
+{public:
+ static const int value = 1;
+};
+
+// Set the problem property
+SET_PROP(OnePTestProblem, Problem)
+{
+ typedef Dumux::OnePTestProblem<TTAG(OnePTestProblem)> type;
+};
+
+// Set the wetting phase
+//SET_TYPE_PROP(OnePTestProblem, Fluid, Dumux::Air);
+
+//! Use the leaf grid view if not defined otherwise
+//SET_PROP(OnePTestProblem, GridView)
+//{
+//private:
+// typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+//
+//public:
+// typedef typename Grid::LevelGridView type;
+//};
+
+// Set the soil properties
+// Set the soil properties
+SET_PROP(OnePTestProblem, SpatialParameters)
+{
+ typedef Dumux::OnePSpatialParameters<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(OnePTestProblem, EnableGravity, true);
+}
+
+/*!
+ * \ingroup OnePBoxProblems
+ * \brief Air flow in porous media
+ *
+ * The domain is box shaped. All sides are closed (Neumann 0 boundary)
+ * except the top and bottom boundaries (Dirichlet), where air is
+ * flowing from bottom to top.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_1p ./grids/1p_2d.dgf 10 0.01</tt>
+ * where start simulation time = 0.01 second, end simulation time = 10 seconds
+ * The same file can be also used for 3d simulation but you need to change line
+ * <tt>typedef Dune::SGrid<2,2> type;</tt> to
+ * <tt>typedef Dune::SGrid<3,3> type;</tt> and use <tt>1p_3d.dgf</tt> grid.
+ */
+template <class TypeTag = TTAG(OnePTestProblem) >
+class OnePTestProblem : public OnePBoxProblem<TypeTag, OnePTestProblem<TypeTag> >
+{
+ typedef OnePTestProblem<TypeTag> ThisType;
+ typedef OnePBoxProblem<TypeTag, ThisType> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePIndices)) Indices;
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ // indices of the primary variables
+ pressureIdx = Indices::pressureIdx,
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
+ typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ OnePTestProblem(const GridView &gridView)
+ : ParentType(gridView)
+ {
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ {
+#ifdef HAVE_DUNE_PDELAB
+ return "1ptest_pdelab";
+#else
+ return "1ptest_disc";
+#endif
+ }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 36 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 10; // 10°C
+ };
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ void boundaryTypes(BoundaryTypeVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ double eps = 1.0e-3;
+ const GlobalPosition &globalPos
+ = fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal;
+
+ if (globalPos[dim-1] < eps || globalPos[dim-1] > 1.0 - eps)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ double eps = 1.0e-3;
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if (globalPos[dim-1] < eps) {
+ values[pressureIdx] = 2.0e+5;
+ }
+ else if (globalPos[dim-1] > 1.0 -eps) {
+ values[pressureIdx] = 1.0e+5;
+ }
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each component. Negative values mean
+ * influx.
+ */
+ void neumann(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ // const GlobalPosition &globalPos = fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal;
+
+ values[pressureIdx] = 0;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ void source(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ values[pressureIdx] = 1.0e+5 + 9.81*1.23*(20-globalPos[dim-1]);
+ }
+
+ // \}
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/1p/1pspatialparameters.hh b/test/boxmodels/1p/1pspatialparameters.hh
new file mode 100644
index 0000000000..42cd504380
--- /dev/null
+++ b/test/boxmodels/1p/1pspatialparameters.hh
@@ -0,0 +1,75 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_1PSPATIALPARAMETERS_HH
+#define DUMUX_1PSPATIALPARAMETERS_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class OnePSpatialParameters : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GridView::ctype CoordScalar;
+
+ enum {
+ dimWorld=GridView::dimensionworld,
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+public:
+
+ OnePSpatialParameters(const GridView& gridView)
+ : ParentType(gridView)
+ {}
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index sub-control volume face where the
+ * intrinsic velocity ought to be calculated.
+ */
+ Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 1e-10;
+ }
+
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 0.4;
+ }
+
+private:
+};
+
+} // end namespace
+#endif
+
diff --git a/test/boxmodels/1p/CMakeLists.txt b/test/boxmodels/1p/CMakeLists.txt
new file mode 100644
index 0000000000..b2dc4c2f8c
--- /dev/null
+++ b/test/boxmodels/1p/CMakeLists.txt
@@ -0,0 +1,18 @@
+add_definitions(-DYASPGRID -DENABLE_ALBERTA -DENABLE_ALUGRID -DENABLE_UG)
+
+# build target for the onephase-onecomponent test problem
+ADD_EXECUTABLE("test_1p" test_1p.cc)
+TARGET_LINK_LIBRARIES("test_1p" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_1p"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/1p/Makefile.am b/test/boxmodels/1p/Makefile.am
new file mode 100644
index 0000000000..98dc4d7666
--- /dev/null
+++ b/test/boxmodels/1p/Makefile.am
@@ -0,0 +1,7 @@
+bin_PROGRAMS = test_1p
+
+test_1p_SOURCES = test_1p.cc
+test_1p_CXXFLAGS = $(DUNEMPICPPFLAGS)
+test_1p_LDADD = $(DUNEMPILDFLAGS)
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/1p/test_1p.cc b/test/boxmodels/1p/test_1p.cc
new file mode 100644
index 0000000000..389989bc5c
--- /dev/null
+++ b/test/boxmodels/1p/test_1p.cc
@@ -0,0 +1,201 @@
+// $Id: test_1p.cc 3779 2010-06-24 07:07:56Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Onur Dogan *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "1pproblem.hh"
+#include <dune/common/exceptions.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+template <class Grid, class Scalar>
+class CreateGrid
+{
+};
+
+#if HAVE_UG
+template <class Scalar>
+class CreateGrid<Dune::UGGrid<2>, Scalar>
+{
+public:
+ static Dune::UGGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ Dune::UGGrid<2> *grid = new Dune::UGGrid<2>;
+ Dune::GridFactory<Dune::UGGrid<2> > factory(grid);
+ for (int i=0; i<=cellRes[0]; i++) {
+ for (int j=0; j<=cellRes[1]; j++) {
+ Dune::FieldVector<double,2> pos;
+ pos[0] = upperRight[0]*double(i)/cellRes[0];
+ pos[1] = upperRight[1]*double(j)/cellRes[1];
+ factory.insertVertex(pos);
+ }
+ }
+
+ for (int i=0; i<cellRes[0]; i++) {
+ for (int j=0; j<cellRes[1]; j++) {
+ std::vector<unsigned int> v(4);
+
+ int i0 = i*(cellRes[1]+1) + j;
+ int i1 = i*(cellRes[1]+1) + j+1;
+ int i2 = (i+1)*(cellRes[1]+1) + j;
+ int i3 = (i+1)*(cellRes[1]+1) + j+1;
+
+ v[0] = i0;
+ v[1] = i1;
+ v[2] = i2;
+ v[3] = i3;
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,2), v);
+ }
+ }
+
+ factory.createGrid();
+ grid->loadBalance();
+ return grid;
+ }
+};
+
+template <class Scalar>
+class CreateGrid<Dune::UGGrid<3>, Scalar>
+{
+public:
+ static Dune::UGGrid<3> *create(const Dune::FieldVector<Scalar, 3> &upperRight,
+ const Dune::FieldVector<int, 3> &cellRes)
+ {
+ Dune::UGGrid<3> *grid = new Dune::UGGrid<3>;
+ Dune::GridFactory<Dune::UGGrid<3> > factory(grid);
+ for (int k=0; k<=cellRes[2]; k++) {
+ for (int j=0; j<=cellRes[1]; j++) {
+ for (int i=0; i<=cellRes[0]; i++) {
+ Dune::FieldVector<double,3> pos;
+ pos[0] = upperRight[0]*double(i)/cellRes[0];
+ pos[1] = upperRight[1]*double(j)/cellRes[1];
+ pos[2] = upperRight[2]*double(k)/cellRes[2];
+ factory.insertVertex(pos);
+ }
+ }
+ }
+
+ for (int k=0; k<cellRes[2]; k++) {
+ for (int j=0; j<cellRes[1]; j++) {
+ for (int i=0; i<cellRes[0]; i++) {
+ std::vector<unsigned int> v(8);
+
+ v[0] = k*(cellRes[0]+1)*(cellRes[1]+1) + j*(cellRes[0]+1) + i;
+ v[1] = k*(cellRes[0]+1)*(cellRes[1]+1) + j*(cellRes[0]+1) + i+1;
+ v[2] = k*(cellRes[0]+1)*(cellRes[1]+1) + (j+1)*(cellRes[0]+1) + i;
+ v[3] = k*(cellRes[0]+1)*(cellRes[1]+1) + (j+1)*(cellRes[0]+1) + i+1;
+ v[4] = (k+1)*(cellRes[0]+1)*(cellRes[1]+1) + j*(cellRes[0]+1) + i;
+ v[5] = (k+1)*(cellRes[0]+1)*(cellRes[1]+1) + j*(cellRes[0]+1) + i+1;
+ v[6] = (k+1)*(cellRes[0]+1)*(cellRes[1]+1) + (j+1)*(cellRes[0]+1) + i;
+ v[7] = (k+1)*(cellRes[0]+1)*(cellRes[1]+1) + (j+1)*(cellRes[0]+1) + i+1;
+
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,3), v);
+ }
+ }
+ }
+
+ factory.createGrid();
+ grid->loadBalance();
+ return grid;
+ }
+};
+#endif
+
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] gridFile.dgf tEnd dt\n")%progname;
+ exit(1);
+};
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(OnePTestProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ enum {dim = Grid::dimensionworld};
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+ typedef Dune::GridPtr<Grid> GridPointer;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ // parse the command line arguments for the program
+ if (argc < 4)
+ usage(argv[0]);
+
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+
+ if (argc - argPos != 3) {
+ usage(argv[0]);
+ }
+
+ double tEnd, dt;
+ const char *dgfFileName = argv[argPos++];
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> dt;
+
+ // create grid
+
+ // -> load the grid from file
+ GridPointer gridPtr = GridPointer(dgfFileName);
+ (*gridPtr).loadBalance();
+ Dune::gridinfo(*gridPtr);
+
+// Grid grid;
+// Dune::GridFactory<Grid> factory(&grid);
+// Dune::GmshReader<Grid>::read(factory,"../../../dune-pdelab-howto/examples/grids/cube1045.msh",true,true);
+// factory.createGrid();
+// grid.loadBalance();
+
+// GlobalPosition lowerLeft(0.0);
+// GlobalPosition upperRight(1.0);
+// Dune::FieldVector<int,dim> res(8); // cell resolution
+// std::auto_ptr<Grid> grid(CreateGrid<Grid, Scalar>::create(upperRight, res));
+// Dune::gridinfo(*grid);
+
+ // instantiate and run the concrete problem
+// Problem problem(gridPtr->levelView(gridPtr->maxLevel()));
+ Problem problem(gridPtr->leafView());
+ if (!problem.simulate(dt, tEnd))
+ return 2;
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ }
+
+ return 3;
+}
diff --git a/test/boxmodels/1p2c/CMakeLists.txt b/test/boxmodels/1p2c/CMakeLists.txt
new file mode 100644
index 0000000000..156f2f8347
--- /dev/null
+++ b/test/boxmodels/1p2c/CMakeLists.txt
@@ -0,0 +1,17 @@
+add_definitions(-DYASPGRID -DENABLE_ALBERTA -DENABLE_ALUGRID -DENABLE_UG)
+
+ADD_EXECUTABLE("test_1p2c" test_1p2c.cc)
+TARGET_LINK_LIBRARIES("test_1p2c" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_1p2c"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/1p2c/Makefile.am b/test/boxmodels/1p2c/Makefile.am
new file mode 100644
index 0000000000..fc07940176
--- /dev/null
+++ b/test/boxmodels/1p2c/Makefile.am
@@ -0,0 +1,5 @@
+bin_PROGRAMS = test_1p2c
+
+test_1p2c_SOURCES = test_1p2c.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/1p2c/test_1p2c.cc b/test/boxmodels/1p2c/test_1p2c.cc
new file mode 100644
index 0000000000..431f9f1bec
--- /dev/null
+++ b/test/boxmodels/1p2c/test_1p2c.cc
@@ -0,0 +1,95 @@
+// $Id: test_1p2c.cc 3779 2010-06-24 07:07:56Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Karin Erbertseder *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "tissue_tumor_problem.hh"
+
+#include <dune/common/exceptions.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] gridFile.dgf tEnd dt\n")%progname;
+ exit(1);
+};
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(TissueTumorProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+ typedef Dune::GridPtr<Grid> GridPointer;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ // parse the command line arguments for the program
+ if (argc < 4)
+ usage(argv[0]);
+
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+
+ if (argc - argPos != 3) {
+ usage(argv[0]);
+ }
+
+ double tEnd, dt;
+ const char *dgfFileName = argv[argPos++];
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> dt;
+
+ // create grid
+
+ // -> load the grid from file
+ GridPointer gridPtr = GridPointer(dgfFileName);
+ (*gridPtr).loadBalance();
+ Dune::gridinfo(*gridPtr);
+
+ // instantiate and run the concrete problem
+ Problem problem(gridPtr->leafView());
+ if (!problem.simulate(dt, tEnd))
+ return 2;
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ }
+
+ return 3;
+}
diff --git a/test/boxmodels/1p2c/tissue_soilproperties.hh b/test/boxmodels/1p2c/tissue_soilproperties.hh
new file mode 100644
index 0000000000..c70c6cd266
--- /dev/null
+++ b/test/boxmodels/1p2c/tissue_soilproperties.hh
@@ -0,0 +1,121 @@
+// $Id: tissue_soilproperties.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Karin Erbertseder *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_TISSUE_SOILPROPERTIES_HH
+#define DUMUX_TISSUE_SOILPROPERTIES_HH
+
+#include <dumux/old_material/property_baseclasses.hh>
+
+namespace Dumux
+{
+
+/**
+ * \brief Definition of the properties of the human tissue
+ *
+ * \todo do not derive from Matrix2p, since this is the base class for
+ * two-phase models
+ */
+template<class Grid, class Scalar>
+class TissueSoil: public Matrix2p<Grid, Scalar>
+{
+public:
+ typedef typename Grid::Traits::template Codim<0>::Entity Element;
+ typedef typename Grid::ctype CoordScalar;
+ enum
+ { dim=Grid::dimension, numEq=2};
+
+ const Dune::FieldMatrix<CoordScalar,dim,dim> &K (const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos) const
+ {
+ if (isTumor_(globalPos))
+ return permTumor_;
+ else
+ return permTissue_;
+ }
+
+ double porosity(const Dune::FieldVector<CoordScalar,dim> &globalPos, const Element &element, const Dune::FieldVector<CoordScalar,dim> &localPos) const
+ {
+ if (isTumor_(globalPos))
+ return 0.31; // tumor tissue
+ else
+ return 0.13; // healthy tissue
+ }
+
+ double tortuosity (const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos) const
+ {
+ if (isTumor_(globalPos))
+ return 0.706; // tumor tissue
+ else
+ return 0.280; // healthy tissue
+ }
+
+ virtual double Sr_w(const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos, const double T = 283.15) const
+ {
+ return 0;
+ }
+
+ virtual double Sr_n(const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos, const double T = 283.15) const
+ {
+ return 0;
+ }
+
+ virtual std::vector<double> paramRelPerm(const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos, const double T = 283.15) const
+ {
+ std::vector<double> param(0);
+ return param;
+ }
+
+ bool useTwoPointGradient(const Element &elem,
+ int vertexI,
+ int vertexJ) const
+ {
+ bool inTumor = isTumor_(elem.geometry().corner(0));
+ int n = elem.template count<dim>();
+ for (int i = 1; i < n; ++i) {
+ if ((inTumor && !isTumor_(elem.geometry().corner(i))) ||
+ (!inTumor && isTumor_(elem.geometry().corner(i))))
+ return true;
+ }
+ return false;
+ }
+
+ TissueSoil()
+ :Matrix2p<Grid,Scalar>()
+ {
+ permTumor_ = 0;
+ permTissue_ = 0;
+ for (int i = 0; i < dim; i++)
+ permTumor_[i][i] = 2.142e-11; //tumor tissue
+ for (int i = 0; i < dim; i++)
+ permTissue_[i][i] = 4.424e-12; //healthy tissue
+ }
+
+private:
+ bool isTumor_(const Dune::FieldVector<CoordScalar,dim>& globalPos) const
+ {
+ if(globalPos[0] > 0.99 && globalPos[0] < 2.01 &&
+ globalPos[1] > 0.99 && globalPos[1] < 3.01)
+ return true;
+ return false;
+ }
+
+ Dune::FieldMatrix<CoordScalar,dim,dim> permTumor_;
+ Dune::FieldMatrix<CoordScalar,dim,dim> permTissue_;
+};
+
+} // end namespace
+
+#endif
diff --git a/test/boxmodels/1p2c/tissue_tumor_problem.hh b/test/boxmodels/1p2c/tissue_tumor_problem.hh
new file mode 100644
index 0000000000..1a447ea232
--- /dev/null
+++ b/test/boxmodels/1p2c/tissue_tumor_problem.hh
@@ -0,0 +1,319 @@
+// $Id: tissue_tumor_problem.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Karin Erbertseder *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Copyright (C) 2008 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/**
+ * \file
+ * \brief Definition of a problem, where the distribution of a therapeutic agent
+ * within pulmonary tissue is described
+ * \author Karin Erbertseder, Bernd Flemisch
+ */
+#ifndef DUMUX_TISSUE_TUMOR_PROBLEM_HH
+#define DUMUX_TISSUE_TUMOR_PROBLEM_HH
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/old_material/fluids/interstitialfluid_trail.hh>
+#include <dumux/boxmodels/1p2c/1p2cboxmodel.hh>
+
+#include "tissue_soilproperties.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class TissueTumorProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(TissueTumorProblem, INHERITS_FROM(BoxOnePTwoC));
+
+// Set the grid type
+SET_PROP(TissueTumorProblem, Grid)
+{
+#if HAVE_UG
+ typedef Dune::UGGrid<2> type;
+#else
+ typedef Dune::SGrid<2, 2> type;
+ //typedef Dune::YaspGrid<2> type;
+#endif
+};
+
+#ifdef HAVE_DUNE_PDELAB
+SET_PROP(TissueTumorProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif
+
+
+// Set the problem property
+SET_PROP(TissueTumorProblem, Problem)
+{
+ typedef Dumux::TissueTumorProblem<TTAG(TissueTumorProblem)> type;
+};
+
+// Set the wetting phase
+SET_TYPE_PROP(TissueTumorProblem, Fluid, Dumux::InterstitialFluid_Trail);
+
+// Set the soil properties
+SET_PROP(TissueTumorProblem, Soil)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::TissueSoil<Grid, Scalar> type;
+};
+
+// Disable gravity
+SET_BOOL_PROP(TissueTumorProblem, EnableGravity, false);
+}
+
+
+/**
+ * \brief Definition of a problem, where the distribution of a therapeutic agent
+ * within pulmonary tissue is described
+ *
+ * The model domain is 22 mm long in x-direction and in y-direction with a discretization length of 0.1
+ * mm. The tumour area is located in the middle of the model domain. The diameter of the tumour is
+ * assumed to be 2 mm.
+ *
+ * The intercapillary distance is in the range of 0.1 mm. So the distance between the grid nodes is
+ * equal to the intercapillary distance. It is assumed that at each node within the model domain the
+ * transition of the therapeutic agent from the blood capillary into the tissue can take place.
+ * Based on this assumption, the initial conditions are adapted. The mole fraction of the dissolved
+ * therapeutic agent x is set to 1.1249 e-8 within the normal pulmonary tissue. Within the tumour the
+ * mole fraction of dissolved therapeutic agent is set to zero due to the assumption of a blood vessel
+ * free tumour. As initial condition for the pressure p, an interstitial fluid pressure of -1067 Pa is
+ * assumed. This value corresponds to the interstitial fluid pressure in healthy pulmonary tissue.
+ *
+ * All four sides of the model domain are described by Dirichlet boundary conditions. The primary
+ * variable p is instantiated with the value -1067 Pa and the primary variable x with the value
+ * 1.1249 e-8.
+ *
+ * The pressure field of the tumour is generated by an additional source term of 1.98 e-9 l/h,
+ * that is set in the center of the tumour region.
+ */
+
+template <class TypeTag = TTAG(TissueTumorProblem) >
+class TissueTumorProblem : public OnePTwoCBoxProblem<TypeTag, TissueTumorProblem<TypeTag> >
+{
+ typedef TissueTumorProblem<TypeTag> ThisType;
+ typedef OnePTwoCBoxProblem<TypeTag, ThisType> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ // indices of the primary variables
+ konti = Indices::konti,
+ transport = Indices::transport,
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
+ typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ TissueTumorProblem(const GridView &gridView)
+ : ParentType(gridView)
+ {
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "tissue"; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 36 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 36; // in [K]
+ };
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ void boundaryTypes(BoundaryTypeVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ values.setAllDirichlet();
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each component. Negative values mean
+ * influx.
+ */
+ void neumann(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+ values = 0;
+
+ //int globalIdx = this->model().vertexMapper().map(element, scvIdx, dim);
+
+ //Scalar lambda = (globalPos[1])/height_;
+ if (globalPos[0] < eps_ ) {
+ values[konti] = -3.8676e-8;
+ //values[transport] = -4.35064e-10;
+ //Robin-Boundary
+ //values[transport] = (*this->model().curSolFunction())[globalIdx][transport];
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ void source(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ values = Scalar(0.0);
+
+ if(globalPos[0]>10 && globalPos[0]<12 && globalPos[1]>10 && globalPos[1]<12)
+ values[0]= 1.5e-6;
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos);
+ }
+
+ // \}
+
+private:
+ // the internal method for the initial condition
+ void initial_(PrimaryVarVector &values,
+ const GlobalPosition &globalPos) const
+ {
+
+ values[konti] = -1067; //initial condition for the pressure
+ values[transport] = 1.1249e-8; //initial condition for the molefraction
+
+ }
+
+ static const Scalar eps_ = 1e-6;
+}; //end namespace
+}
+#endif
diff --git a/test/boxmodels/2p/CMakeLists.txt b/test/boxmodels/2p/CMakeLists.txt
new file mode 100644
index 0000000000..09bec21429
--- /dev/null
+++ b/test/boxmodels/2p/CMakeLists.txt
@@ -0,0 +1,18 @@
+add_definitions(-DYASPGRID -DGRIDDIM=2 -DENABLE_UG)
+
+# build target for the simple twophase lens problem
+ADD_EXECUTABLE("test_2p" test_2p.cc)
+TARGET_LINK_LIBRARIES("test_2p" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_2p"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/2p/Makefile.am b/test/boxmodels/2p/Makefile.am
new file mode 100644
index 0000000000..0c7a5a900f
--- /dev/null
+++ b/test/boxmodels/2p/Makefile.am
@@ -0,0 +1,5 @@
+bin_PROGRAMS = test_2p
+
+test_2p_SOURCES = test_2p.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/2p/lensproblem.hh b/test/boxmodels/2p/lensproblem.hh
new file mode 100644
index 0000000000..30369c17c6
--- /dev/null
+++ b/test/boxmodels/2p/lensproblem.hh
@@ -0,0 +1,394 @@
+// $Id: lensproblem.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_LENSPROBLEM_HH
+#define DUMUX_LENSPROBLEM_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/simplednapl.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+
+#include <dumux/boxmodels/2p/2pboxmodel.hh>
+
+#include "lensspatialparameters.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class LensProblem;
+
+//////////
+// Specify the properties for the lens problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(LensProblem, INHERITS_FROM(BoxTwoP));
+
+// Set the grid type
+SET_PROP(LensProblem, Grid)
+{
+#if HAVE_UG
+ typedef Dune::UGGrid<2> type;
+#else
+ typedef Dune::YaspGrid<2> type;
+ //typedef Dune::SGrid<2, 2> type;
+#endif
+};
+
+#ifdef HAVE_DUNE_PDELAB
+SET_PROP(LensProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+#if CUBES
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+#else
+ typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+#endif
+};
+#endif
+
+// Set the problem property
+SET_PROP(LensProblem, Problem)
+{
+ typedef Dumux::LensProblem<TypeTag> type;
+};
+
+// Set the wetting phase
+SET_PROP(LensProblem, WettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(LensProblem, NonwettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleDNAPL<Scalar> > type;
+};
+
+// Set the soil properties
+SET_PROP(LensProblem, SpatialParameters)
+{
+ typedef Dumux::LensSpatialParameters<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(LensProblem, EnableGravity, true);
+}
+
+/*!
+ * \ingroup TwoPBoxProblems
+ * \brief Soil decontamination problem where DNAPL infiltrates a fully
+ * water saturated medium.
+ *
+ * The domain is sized 6m times 4m and features a rectangular lens
+ * with low permeablility which spans from (1 m , 2 m) to (4 m, 3 m)
+ * and is surrounded by a medium with higher permability.
+ *
+ * On the top and the bottom of the domain neumann boundary conditions
+ * are used, while dirichlet conditions apply on the left and right
+ * boundaries.
+ *
+ * DNAPL is injected at the top boundary from 3m to 4m at a rate of
+ * 0.04 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * The dirichlet boundaries on the left boundary is the hydrostatic
+ * pressure scaled by a factor of 1.125, while on the right side it is
+ * just the hydrostatic pressure. The DNAPL saturation on both sides
+ * is zero.
+ *
+ * This problem uses the \ref TwoPBoxModel.
+ *
+ * This problem should typically simulated until \f$t_{\text{end}} =
+ * 50\,000\;s\f$ is reached. A good choice for the initial time step size
+ * is \f$t_{\text{inital}} = 1\,000\;s\f$.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_2p 50000 1000</tt>
+ */
+template <class TypeTag >
+class LensProblem : public TwoPBoxProblem<TypeTag,
+ LensProblem<TypeTag> >
+{
+ typedef LensProblem<TypeTag> ThisType;
+ typedef TwoPBoxProblem<TypeTag, ThisType> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef TwoPFluidState<TypeTag> FluidState;
+
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ // primary variable indices
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+ pwIdx = Indices::pwIdx,
+ SnIdx = Indices::SnIdx,
+
+ // equation indices
+ contiWEqIdx = Indices::contiWEqIdx,
+ contiNEqIdx = Indices::contiNEqIdx,
+
+ // phase indices
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
+ typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ LensProblem(const GridView &gridView,
+ const GlobalPosition &lensLowerLeft,
+ const GlobalPosition &lensUpperRight)
+ : ParentType(gridView)
+ {
+ this->spatialParameters().setLensCoords(lensLowerLeft, lensUpperRight);
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "lens"; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 10; // -> 10°C
+ };
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ void boundaryTypes(BoundaryTypeVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) {
+ values.setAllDirichlet();
+ }
+ else {
+ values.setAllNeumann();
+ //values.setDirichlet(SnIdx, pwIdx);
+ //values.setNeumann(SnIdx);
+ }
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ Scalar densityW = FluidSystem::componentDensity(wPhaseIdx,
+ wPhaseIdx,
+ temperature(element, fvElemGeom, scvIdx),
+ 1e5);
+
+ if (onLeftBoundary_(globalPos))
+ {
+ Scalar height = this->bboxMax()[1] - this->bboxMin()[1];
+ Scalar depth = this->bboxMax()[1] - globalPos[1];
+ Scalar alpha = (1 + 1.5/height);
+
+ // hydrostatic pressure scaled by alpha
+ values[pwIdx] = - alpha*densityW*this->gravity()[1]*depth;
+ values[SnIdx] = 0.0;
+ }
+ else if (onRightBoundary_(globalPos))
+ {
+ Scalar depth = this->bboxMax()[1] - globalPos[1];
+
+ // hydrostatic pressure
+ values[pwIdx] = -densityW*this->gravity()[1]*depth;
+ values[SnIdx] = 0.0;
+ }
+ else
+ values = 0.0;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+
+ values = 0.0;
+ if (onInlet_(globalPos)) {
+ values[contiNEqIdx] = -0.04; // kg / (m * s)
+ }
+ }
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * generated or annihilate per volume unit. Positive values mean
+ * that mass is created, negative ones mean that it vanishes.
+ */
+ void source(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &,
+ int subControlVolumeIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ // no DNAPL, some random pressure
+ values[pwIdx] = 0.0;
+ values[SnIdx] = 0.0;
+ }
+ // \}
+
+private:
+
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] < this->bboxMin()[0] + eps_;
+ }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] > this->bboxMax()[0] - eps_;
+ }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] < this->bboxMin()[1] + eps_;
+ }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] > this->bboxMax()[1] - eps_;
+ }
+
+ bool onInlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar width = this->bboxMax()[0] - this->bboxMin()[0];
+ Scalar lambda = (this->bboxMax()[0] - globalPos[0])/width;
+ return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
+ }
+
+ static const Scalar eps_ = 3e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/2p/lensspatialparameters.hh b/test/boxmodels/2p/lensspatialparameters.hh
new file mode 100644
index 0000000000..c1c44d07b9
--- /dev/null
+++ b/test/boxmodels/2p/lensspatialparameters.hh
@@ -0,0 +1,164 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Markus Wolff *
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_LENSSPATIALPARAMETERS_HH
+#define DUMUX_LENSSPATIALPARAMETERS_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+/**
+ * @file
+ * @brief Class for defining spatial parameters
+ * @author Bernd Flemisch, Klaus Mosthaf, Markus Wolff
+ */
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class LensSpatialParameters : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ // define the material law which is parameterized by effective
+ // saturations
+ typedef RegularizedVanGenuchten<Scalar> EffectiveLaw;
+
+public:
+ // define the material law parameterized by absolute saturations
+ typedef EffToAbsLaw<EffectiveLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ LensSpatialParameters(const GridView& gridView)
+ : ParentType(gridView)
+ {
+ // residual saturations
+ lensMaterialParams_.setSwr(0.18);
+ lensMaterialParams_.setSnr(0.0);
+ outerMaterialParams_.setSwr(0.05);
+ outerMaterialParams_.setSnr(0.0);
+
+// // parameters for the Van Genuchten law
+// // alpha and n
+ lensMaterialParams_.setVgAlpha(0.0037);
+ outerMaterialParams_.setVgAlpha(0.00045);
+ lensMaterialParams_.setVgN(4.7);
+ outerMaterialParams_.setVgN(7.3);
+
+ // parameters for the linear law
+ // minimum and maximum pressures
+// lensMaterialParams_.setEntryPC(0);
+// outerMaterialParams_.setEntryPC(0);
+// lensMaterialParams_.setMaxPC(0);
+// outerMaterialParams_.setMaxPC(0);
+
+ lensK_ = 1e-13;
+ outerK_ = 5e-10;
+ }
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvIdx The index sub-control volume face where the
+ * intrinsic velocity ought to be calculated.
+ */
+ Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+ if (isInLens_(globalPos))
+ return lensK_;
+ return outerK_;
+ }
+
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ { return 0.4; }
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvElemGeom.subContVol[scvIdx].global;
+
+ if (isInLens_(globalPos))
+ return lensMaterialParams_;
+ return outerMaterialParams_;
+ }
+
+
+ //! Set the bounding box of the fine-sand lens
+ void setLensCoords(const GlobalPosition& lensLowerLeft,
+ const GlobalPosition& lensUpperRight)
+ {
+ lensLowerLeft_ = lensLowerLeft;
+ lensUpperRight_ = lensUpperRight;
+ }
+
+private:
+ bool isInLens_(const GlobalPosition &pos) const
+ {
+ for (int i = 0; i < dim; ++i) {
+ if (pos[i] < lensLowerLeft_[i] || pos[i] > lensUpperRight_[i])
+ return false;
+ }
+ return true;
+ }
+
+ GlobalPosition lensLowerLeft_;
+ GlobalPosition lensUpperRight_;
+
+ Scalar lensK_;
+ Scalar outerK_;
+ MaterialLawParams lensMaterialParams_;
+ MaterialLawParams outerMaterialParams_;
+};
+
+} // end namespace
+#endif
+
diff --git a/test/boxmodels/2p/test_2p.cc b/test/boxmodels/2p/test_2p.cc
new file mode 100644
index 0000000000..53d64b3b16
--- /dev/null
+++ b/test/boxmodels/2p/test_2p.cc
@@ -0,0 +1,224 @@
+// $Id: test_2p.cc 3779 2010-06-24 07:07:56Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#define CUBES 1
+
+#include "lensproblem.hh"
+
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+////////////////////////
+// helper class for grid instantiation
+////////////////////////
+template <class Grid, class Scalar>
+class CreateGrid
+{
+};
+
+#if HAVE_UG
+template <class Scalar>
+class CreateGrid<Dune::UGGrid<2>, Scalar>
+{
+public:
+ static Dune::UGGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ Dune::UGGrid<2> *grid = new Dune::UGGrid<2>;
+ Dune::GridFactory<Dune::UGGrid<2> > factory(grid);
+ for (int i=0; i<=cellRes[0]; i++) {
+ for (int j=0; j<=cellRes[1]; j++) {
+ Dune::FieldVector<double,2> pos;
+ pos[0] = upperRight[0]*double(i)/cellRes[0];
+ pos[1] = upperRight[1]*double(j)/cellRes[1];
+ factory.insertVertex(pos);
+ }
+ }
+
+ for (int i=0; i<cellRes[0]; i++) {
+ for (int j=0; j<cellRes[1]; j++) {
+#if CUBES
+ std::vector<unsigned int> v(4);
+#else
+ std::vector<unsigned int> v(3);
+#endif
+
+ int i0 = i*(cellRes[1]+1) + j;
+ int i1 = i*(cellRes[1]+1) + j+1;
+ int i2 = (i+1)*(cellRes[1]+1) + j;
+ int i3 = (i+1)*(cellRes[1]+1) + j+1;
+
+#if CUBES
+ v[0] = i0;
+ v[1] = i1;
+ v[2] = i2;
+ v[3] = i3;
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,2), v);
+#else
+ v[0] = i0;
+ v[1] = i1;
+ v[2] = i2;
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,2), v);
+
+ v[0] = i1;
+ v[1] = i2;
+ v[2] = i3;
+ factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,2), v);
+#endif
+ }
+ }
+
+ factory.createGrid();
+ grid->loadBalance();
+ return grid;
+ }
+};
+#endif
+
+template <class Scalar>
+class CreateGrid<Dune::YaspGrid<2>, Scalar>
+{
+public:
+ static Dune::YaspGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ return new Dune::YaspGrid<2>(
+#ifdef HAVE_MPI
+ Dune::MPIHelper::getCommunicator(),
+#endif
+ upperRight, // upper right
+ cellRes, // number of cells
+ Dune::FieldVector<bool,2>(false), // periodic
+ 4); // overlap
+ };
+};
+
+template <class Scalar>
+class CreateGrid<Dune::SGrid<2, 2>, Scalar>
+{
+public:
+ static Dune::SGrid<2, 2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
+ const Dune::FieldVector<int, 2> &cellRes)
+ {
+ return new Dune::SGrid<2,2>(cellRes, // number of cells
+ Dune::FieldVector<Scalar, 2>(0.0), // lower left
+ upperRight); // upper right
+
+ };
+};
+
+////////////////////////
+// the main function
+////////////////////////
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] tEnd dt\n")%progname;
+ exit(1);
+}
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(LensProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+ if (argc < 3)
+ usage(argv[0]);
+
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+
+ if (argc - argPos != 2) {
+ usage(argv[0]);
+ }
+
+ // read the initial time step and the end time
+ double tEnd, dt;
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> dt;
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ GlobalPosition lowerLeft(0.0);
+ GlobalPosition upperRight;
+ Dune::FieldVector<int,dim> res; // cell resolution
+ upperRight[0] = 6.0;
+ upperRight[1] = 4.0;
+ res[0] = 48;
+ res[1] = 32;
+
+ std::auto_ptr<Grid> grid(CreateGrid<Grid, Scalar>::create(upperRight, res));
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ // specify dimensions of the low-permeable lens
+ GlobalPosition lowerLeftLens, upperRightLens;
+ lowerLeftLens[0] = 1.0;
+ lowerLeftLens[1] = 2.0;
+ upperRightLens[0] = 4.0;
+ upperRightLens[1] = 3.0;
+ Problem problem(grid->leafView(), lowerLeftLens, upperRightLens);
+
+ // load restart file if necessarry
+ if (restart)
+ problem.deserialize(restartTime);
+
+ // run the simulation
+ if (!problem.simulate(dt, tEnd))
+ return 2;
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ throw;
+ }
+
+ return 3;
+}
diff --git a/test/boxmodels/2p2c/CMakeLists.txt b/test/boxmodels/2p2c/CMakeLists.txt
new file mode 100644
index 0000000000..ad33076f53
--- /dev/null
+++ b/test/boxmodels/2p2c/CMakeLists.txt
@@ -0,0 +1,18 @@
+add_definitions(-DYASPGRID -DENABLE_ALBERTA -DENABLE_ALUGRID -DENABLE_UG)
+
+# build target for the twophase-twocomponent test problem
+ADD_EXECUTABLE("test_2p2c" test_2p2c.cc)
+TARGET_LINK_LIBRARIES("test_2p2c" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_2p2c"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/2p2c/Makefile.am b/test/boxmodels/2p2c/Makefile.am
new file mode 100644
index 0000000000..fe6d530fef
--- /dev/null
+++ b/test/boxmodels/2p2c/Makefile.am
@@ -0,0 +1,5 @@
+bin_PROGRAMS = test_2p2c
+
+test_2p2c_SOURCES = test_2p2c.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/2p2c/injectionproblem.hh b/test/boxmodels/2p2c/injectionproblem.hh
new file mode 100644
index 0000000000..1074765554
--- /dev/null
+++ b/test/boxmodels/2p2c/injectionproblem.hh
@@ -0,0 +1,321 @@
+// $Id: injectionproblem.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/**
+ * @file
+ * \ingroup TwoPTwoCBoxProblems
+ * @brief Definition of a problem, where air is injected under a low permeable layer
+ * @author Klaus Mosthaf, Andreas Lauser, Bernd Flemisch
+ */
+#ifndef DUMUX_INJECTIONPROBLEM_HH
+#define DUMUX_INJECTIONPROBLEM_HH
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/boxmodels/2p2c/2p2cboxmodel.hh>
+
+#include <dumux/material/fluidsystems/h2o_n2_system.hh>
+
+//#include <dumux/material/fluidsystems/brine_co2_system.hh>
+//#include <appl/co2/ifp/ifpco2tables.hh>
+
+#include "injectionspatialparameters.hh"
+
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class InjectionProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(BoxTwoPTwoC));
+
+// Set the grid type
+SET_PROP(InjectionProblem, Grid)
+{
+ typedef Dune::SGrid<2,2> type;
+};
+
+#ifdef HAVE_DUNE_PDELAB
+SET_PROP(InjectionProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif
+
+// Set the problem property
+SET_PROP(InjectionProblem, Problem)
+{
+ typedef Dumux::InjectionProblem<TTAG(InjectionProblem)> type;
+};
+
+// Set fluid configuration
+SET_PROP(InjectionProblem,
+ FluidSystem)
+{
+ //typedef Dumux::Brine_CO2_System<TypeTag, Dumux::IFP::CO2Tables> type;
+ typedef Dumux::H2O_N2_System<TypeTag> type;
+};
+
+// Set the soil properties
+SET_TYPE_PROP(InjectionProblem,
+ SpatialParameters,
+ Dumux::InjectionSpatialParameters<TypeTag>);
+
+// Enable gravity
+SET_BOOL_PROP(InjectionProblem, EnableGravity, true);
+
+// Enable gravity
+SET_INT_PROP(InjectionProblem, NewtonLinearSolverVerbosity, 0);
+}
+
+
+/*!
+ * \ingroup TwoPBoxProblems
+ * \brief Problem where air is injected under a low permeable layer in a depth of 800m.
+ *
+ * The domain is sized 60m times 40m and consists of two layers, a moderately
+ * permeable soil (\f$ K=10e-12\f$) for \f$ y>22m\f$ and one with a lower permeablility (\f$ K=10e-13\f$)
+ * in the rest of the domain.
+ *
+ * Air enters a water-filled aquifer, which is situated 800m below sea level, at the right boundary
+ * (\f$ 5m<y<15m\f$) and migrates upwards due to buoyancy. It accumulates and
+ * partially enters the lower permeable aquitard.
+ * This problem uses the \ref TwoPTwoCBoxModel.
+ */
+template <class TypeTag = TTAG(InjectionProblem) >
+class InjectionProblem : public TwoPTwoCBoxProblem<TypeTag, InjectionProblem<TypeTag> >
+{
+ typedef InjectionProblem<TypeTag> ThisType;
+ typedef TwoPTwoCBoxProblem<TypeTag, ThisType> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ contiLEqIdx = Indices::contiLEqIdx,
+ contiGEqIdx = Indices::contiGEqIdx,
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
+ typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ InjectionProblem(const GridView &gridView)
+ : ParentType(gridView)
+ {
+ // initialize the tables of the fluid system
+ FluidSystem::init();
+ }
+
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "injection"; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return temperature_;
+ };
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ void boundaryTypes(BoundaryTypeVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if (globalPos[0] < eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each component. Negative values mean
+ * influx.
+ */
+ void neumann(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ values = 0;
+ if (globalPos[1] < 15 && globalPos[1] > 5) {
+ values[contiGEqIdx] = -1e-3; // kg/(s*m^2)
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ void source(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ */
+ int initialPhasePresence(const Vertex &vert,
+ int &globalIdx,
+ const GlobalPosition &globalPos) const
+ { return Indices::lPhaseOnly; }
+
+ // \}
+
+private:
+ // the internal method for the initial condition
+ void initial_(PrimaryVarVector &values,
+ const GlobalPosition &globalPos) const
+ {
+ Scalar densityW = FluidSystem::H2O::liquidDensity(temperature_, 1e5);
+
+ values[Indices::plIdx] = 1e5 - densityW*this->gravity()[1]*(depthBOR_ - globalPos[1]);
+ values[Indices::SgOrXIdx] =
+ values[Indices::plIdx]*0.95/
+ BinaryCoeff::H2O_N2::henry(temperature_);
+ }
+
+ static const Scalar temperature_ = 273.15 + 40; // [K]
+ static const Scalar depthBOR_ = 2700.0; // [m]
+ static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/2p2c/injectionspatialparameters.hh b/test/boxmodels/2p2c/injectionspatialparameters.hh
new file mode 100644
index 0000000000..f0551db1d7
--- /dev/null
+++ b/test/boxmodels/2p2c/injectionspatialparameters.hh
@@ -0,0 +1,243 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_INJECTION_SPATIAL_PARAMETERS_HH
+#define DUMUX_INJECTION_SPATIAL_PARAMETERS_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/**
+ * \brief Definition of the soil properties for the injection problem
+ *
+ */
+template<class TypeTag>
+class InjectionSpatialParameters : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ enum {
+ lPhaseIdx = FluidSystem::wPhaseIdx,
+ gPhaseIdx = FluidSystem::nPhaseIdx,
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> Vector;
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::SolutionVector SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexData)) VertexData;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxData)) FluxData;
+ typedef std::vector<VertexData> VertexDataArray;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
+ //typedef LinearMaterial<Scalar> EffMaterialLaw;
+public:
+ typedef EffToAbsLaw<EffMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ InjectionSpatialParameters(const GridView &gv)
+ : ParentType(gv)
+ {
+ layerBottom_ = 22.0;
+
+ // intrinsic permeabilities
+ fineK_ = 1e-13;
+ coarseK_ = 1e-12;
+
+ // porosities
+ finePorosity_ = 0.3;
+ coarsePorosity_ = 0.3;
+
+ // residual saturations
+ fineMaterialParams_.setSwr(0.2);
+ fineMaterialParams_.setSnr(0.0);
+ coarseMaterialParams_.setSwr(0.2);
+ coarseMaterialParams_.setSnr(0.0);
+
+ // parameters for the Brooks-Corey law
+ fineMaterialParams_.setPe(1e4);
+ coarseMaterialParams_.setPe(1e4);
+ fineMaterialParams_.setAlpha(2.0);
+ coarseMaterialParams_.setAlpha(2.0);
+ }
+
+ ~InjectionSpatialParameters()
+ {}
+
+
+ /*!
+ * \brief Update the spatial parameters with the flow solution
+ * after a timestep.
+ *
+ * \TODO call interface
+ */
+ void update(const SolutionVector &globalSolution)
+ {
+ };
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvfIdx The index sub-control volume face where the
+ * intrinsic velocity ought to be calculated.
+ */
+ const Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineK_;
+ return coarseK_;
+ }
+
+ /*!
+ * \brief Define the porosity \f$[-]\f$ of the soil
+ *
+ * \param vDat The data defined on the sub-control volume
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the porosity needs to be defined
+ */
+ double porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return finePorosity_;
+ return coarsePorosity_;
+ }
+
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineMaterialParams_;
+ return coarseMaterialParams_;
+ }
+
+ /*!
+ * \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the heat capacity needs to be defined
+ */
+ double heatCapacity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return
+ 790 // specific heat capacity of granite [J / (kg K)]
+ * 2700 // density of granite [kg/m^3]
+ * (1 - porosity(element, fvElemGeom, scvIdx));
+ }
+
+ /*!
+ * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
+ * rock matrix based on the temperature gradient \f$[K / m]\f$
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param heatFlux The result vector
+ * \param tempGrad The temperature gradient
+ * \param element The current finite element
+ * \param fvElemGeom The finite volume geometry of the current element
+ * \param scvfIdx The local index of the sub-control volume face where
+ * the matrix heat flux should be calculated
+ */
+ void matrixHeatFlux(Vector &heatFlux,
+ const FluxData &fluxDat,
+ const VertexDataArray &vDat,
+ const Vector &tempGrad,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvfIdx) const
+ {
+ static const Scalar lWater = 0.6;
+ static const Scalar lGranite = 2.8;
+
+ // arithmetic mean of the liquid saturation and the porosity
+ const int i = fvElemGeom.subContVolFace[scvfIdx].i;
+ const int j = fvElemGeom.subContVolFace[scvfIdx].j;
+ Scalar Sl = std::max(0.0, (vDat[i].saturation(lPhaseIdx) +
+ vDat[j].saturation(lPhaseIdx)) / 2);
+ Scalar poro = (porosity(element, fvElemGeom, i) +
+ porosity(element, fvElemGeom, j)) / 2;
+
+ Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
+ Scalar ldry = pow(lGranite, (1-poro));
+
+ // the heat conductivity of the matrix. in general this is a
+ // tensorial value, but we assume isotropic heat conductivity.
+ Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
+
+ // the matrix heat flux is the negative temperature gradient
+ // times the heat conductivity.
+ heatFlux = tempGrad;
+ heatFlux *= -heatCond;
+ }
+
+private:
+ bool isFineMaterial_(const GlobalPosition &pos) const
+ { return pos[dim-1] > layerBottom_; };
+
+ Scalar fineK_;
+ Scalar coarseK_;
+ Scalar layerBottom_;
+
+ Scalar finePorosity_;
+ Scalar coarsePorosity_;
+
+ MaterialLawParams fineMaterialParams_;
+ MaterialLawParams coarseMaterialParams_;
+};
+
+}
+
+#endif
diff --git a/test/boxmodels/2p2c/test_2p2c.cc b/test/boxmodels/2p2c/test_2p2c.cc
new file mode 100644
index 0000000000..c52c0d4d88
--- /dev/null
+++ b/test/boxmodels/2p2c/test_2p2c.cc
@@ -0,0 +1,24 @@
+// $Id: test_2p2c.cc 3779 2010-06-24 07:07:56Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+#include "injectionproblem.hh"
+#include <dumux/common/start.hh>
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(InjectionProblem) ProblemTypeTag;
+ return Dumux::startFromDGF<ProblemTypeTag>(argc, argv);
+}
diff --git a/test/boxmodels/2p2cni/CMakeLists.txt b/test/boxmodels/2p2cni/CMakeLists.txt
new file mode 100644
index 0000000000..c4f52acea6
--- /dev/null
+++ b/test/boxmodels/2p2cni/CMakeLists.txt
@@ -0,0 +1,18 @@
+add_definitions(-DENABLE_UG)
+
+# build target for the 2p2cni test problem
+ADD_EXECUTABLE("test_2p2cni" test_2p2cni.cc)
+TARGET_LINK_LIBRARIES("test_2p2cni" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_2p2cni"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/2p2cni/Makefile.am b/test/boxmodels/2p2cni/Makefile.am
new file mode 100644
index 0000000000..90364a3ade
--- /dev/null
+++ b/test/boxmodels/2p2cni/Makefile.am
@@ -0,0 +1,5 @@
+bin_PROGRAMS = test_2p2cni
+
+test_2p2cni_SOURCES = test_2p2cni.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/2p2cni/test_2p2cni.cc b/test/boxmodels/2p2cni/test_2p2cni.cc
new file mode 100644
index 0000000000..4b4f64a07c
--- /dev/null
+++ b/test/boxmodels/2p2cni/test_2p2cni.cc
@@ -0,0 +1,75 @@
+// $Id: test_2p2cni.cc 3779 2010-06-24 07:07:56Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Melanie Darcis *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "waterairproblem.hh"
+
+#include <dune/common/exceptions.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef GET_PROP_TYPE(TTAG(WaterAirProblem), PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TTAG(WaterAirProblem), PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TTAG(WaterAirProblem), PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+ typedef Dune::GridPtr<Grid> GridPointer;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ // parse the command line arguments for the program
+ if (argc != 4) {
+ std::cout << boost::format("usage: %s grid tEnd dt\n")%argv[0];
+ return 1;
+ }
+ double tEnd, dt;
+ const char *dgfFileName = argv[1];
+ std::istringstream(argv[2]) >> tEnd;
+ std::istringstream(argv[3]) >> dt;
+
+ // create grid
+
+ // -> load the grid from file
+ GridPointer gridPtr = GridPointer(dgfFileName);
+ (*gridPtr).loadBalance();
+ Dune::gridinfo(*gridPtr);
+
+
+ // instantiate and run the concrete problem
+ Problem problem(gridPtr->leafView());
+ if (!problem.simulate(dt, tEnd))
+ return 2;
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ }
+ return 3;
+}
diff --git a/test/boxmodels/2p2cni/waterairproblem.hh b/test/boxmodels/2p2cni/waterairproblem.hh
new file mode 100644
index 0000000000..824ebb8b38
--- /dev/null
+++ b/test/boxmodels/2p2cni/waterairproblem.hh
@@ -0,0 +1,353 @@
+// $Id: waterairproblem.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Klaus Mosthaf *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_WATERAIRPROBLEM_HH
+#define DUMUX_WATERAIRPROBLEM_HH
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/material/fluidsystems/h2o_n2_system.hh>
+
+#include <dumux/boxmodels/2p2cni/2p2cniboxmodel.hh>
+
+#include "waterairspatialparameters.hh"
+
+#define ISOTHERMAL 0
+
+/*!
+ * \ingroup BoxProblems
+ * \brief TwoPTwoCNIBoxProblems Non-isothermal two-phase two-component box problems
+ */
+
+namespace Dumux
+{
+template <class TypeTag>
+class WaterAirProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(WaterAirProblem, INHERITS_FROM(BoxTwoPTwoCNI));
+
+// Set the grid type
+SET_PROP(WaterAirProblem, Grid)
+{
+ typedef Dune::SGrid<2,2> type;
+};
+
+#ifdef HAVE_DUNE_PDELAB
+SET_PROP(WaterAirProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif
+
+// Set the problem property
+SET_PROP(WaterAirProblem, Problem)
+{
+ typedef Dumux::WaterAirProblem<TTAG(WaterAirProblem)> type;
+};
+
+// Set the wetting phase
+SET_TYPE_PROP(WaterAirProblem, FluidSystem, Dumux::H2O_N2_System<TypeTag>);
+
+// Set the soil properties
+SET_TYPE_PROP(WaterAirProblem,
+ SpatialParameters,
+ Dumux::WaterAirSpatialParameters<TypeTag>);
+
+// Enable gravity
+SET_BOOL_PROP(WaterAirProblem, EnableGravity, true);
+
+// Write newton convergence
+SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, true);
+}
+
+
+/*!
+ * \ingroup TwoPTwoCNIBoxProblems
+ * \brief Nonisothermal gas injection problem where a gas (e.g. air) is injected into a fully
+ * water saturated medium. During buoyancy driven upward migration the gas
+ * passes a high temperature area.
+ *
+ * The domain is sized 40 m times 40 m. The rectangular area with the increased temperature (380 K)
+ * starts at (20 m, 5 m) and ends at (30 m, 35 m).
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top and on the bottom of the domain, while dirichlet conditions
+ * apply on the left and the right boundary.
+ * For the energy conservation equation dirichlet boundary conditions are applied
+ * on all boundaries.
+ *
+ * Gas is injected at the bottom boundary from 15 m to 25 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure, a gas saturation of zero and
+ * a geothermal temperature gradient of 0.03 K/m are applied.
+ *
+ * This problem uses the \ref TwoPTwoCNIBoxModel.
+ *
+ * This problem should typically be simulated for 300000 s.
+ * A good choice for the initial time step size is 1000 s.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_2p2cni ./grids/test_2p2cni.dgf 300000 1000</tt>
+ * */
+template <class TypeTag = TTAG(WaterAirProblem) >
+class WaterAirProblem : public TwoPTwoCNIBoxProblem<TypeTag, WaterAirProblem<TypeTag> >
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GridView::Grid Grid;
+
+ typedef WaterAirProblem<TypeTag> ThisType;
+ typedef TwoPTwoCNIBoxProblem<TypeTag, ThisType> ParentType;
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx,
+#if !ISOTHERMAL
+ temperatureIdx = Indices::temperatureIdx,
+ energyEqIdx = Indices::energyEqIdx,
+#endif
+
+ // Phase State
+ lPhaseOnly = Indices::lPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ bothPhases = Indices::bothPhases,
+
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
+ typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ WaterAirProblem(const GridView &gridView)
+ : ParentType(gridView)
+ {
+ FluidSystem::init();
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "waterair"; }
+
+#if ISOTHERMAL
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 10; // -> 10°C
+ };
+#endif
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ void boundaryTypes(BoundaryTypeVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if(globalPos[0] > 40 - eps_ || globalPos[0] < eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+
+#if !ISOTHERMAL
+ values.setDirichlet(temperatureIdx, energyEqIdx);
+#endif
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = element.geometry().corner(scvIdx);
+ // const LocalPosition &localPos
+ // = DomainTraits::referenceElement(element.geometry().type()).position(dim,scvIdx);
+
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each component. Negative values mean
+ * influx.
+ */
+ void neumann(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ values = 0;
+
+ // negative values for injection
+ if (globalPos[0] > 15 && globalPos[0] < 25 &&
+ globalPos[1] < eps_)
+ {
+ values[Indices::contiGEqIdx] = -1e-3;
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ void source(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ initial_(values, globalPos);
+
+#if !ISOTHERMAL
+ if (globalPos[0] > 20 && globalPos[0] < 30 && globalPos[1] < 30)
+ values[temperatureIdx] = 380;
+#endif
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ */
+ int initialPhasePresence(const Vertex &vert,
+ int &globalIdx,
+ const GlobalPosition &globalPos) const
+ {
+ return lPhaseOnly;
+ }
+
+private:
+ // internal method for the initial condition (reused for the
+ // dirichlet conditions!)
+ void initial_(PrimaryVarVector &values,
+ const GlobalPosition &globalPos) const
+ {
+ Scalar densityW = 1000.0;
+ values[pressureIdx] = 1e5 + (depthBOR_ - globalPos[1])*densityW*9.81;
+ values[switchIdx] = 0.0;
+#if !ISOTHERMAL
+ values[temperatureIdx] = 283.0 + (depthBOR_ - globalPos[1])*0.03;
+#endif
+ }
+
+ static const Scalar depthBOR_ = 1000.0; // [m]
+ static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/2p2cni/waterairspatialparameters.hh b/test/boxmodels/2p2cni/waterairspatialparameters.hh
new file mode 100644
index 0000000000..1874a5a35f
--- /dev/null
+++ b/test/boxmodels/2p2cni/waterairspatialparameters.hh
@@ -0,0 +1,247 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_WATER_AIR_SPATIAL_PARAMETERS_HH
+#define DUMUX_WATER_AIR_SPATIAL_PARAMETERS_HH
+
+#include <dumux/material/spatialparameters/boxspatialparameters.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/**
+ * \brief Definition of the soil properties for the water-air problem
+ *
+ */
+template<class TypeTag>
+class WaterAirSpatialParameters : public BoxSpatialParameters<TypeTag>
+{
+ typedef BoxSpatialParameters<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
+ enum {
+ lPhaseIdx = Indices::lPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> Vector;
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::SolutionVector SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexData)) VertexData;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxData)) FluxData;
+ typedef std::vector<VertexData> VertexDataArray;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef RegularizedBrooksCorey<Scalar> EffMaterialLaw;
+public:
+ typedef EffToAbsLaw<EffMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ WaterAirSpatialParameters(const GridView &gv)
+ : ParentType(gv)
+ {
+ layerBottom_ = 22.0;
+
+ // intrinsic permeabilities
+ fineK_ = 1e-13;
+ coarseK_ = 1e-12;
+
+ // porosities
+ finePorosity_ = 0.3;
+ coarsePorosity_ = 0.3;
+
+ // residual saturations
+ fineMaterialParams_.setSwr(0.2);
+ fineMaterialParams_.setSnr(0.0);
+ coarseMaterialParams_.setSwr(0.2);
+ coarseMaterialParams_.setSnr(0.0);
+
+ // parameters for the Brooks-Corey law
+ fineMaterialParams_.setPe(1e4);
+ coarseMaterialParams_.setPe(1e4);
+ fineMaterialParams_.setAlpha(2.0);
+ coarseMaterialParams_.setAlpha(2.0);
+ }
+
+ ~WaterAirSpatialParameters()
+ {}
+
+
+ /*!
+ * \brief Update the spatial parameters with the flow solution
+ * after a timestep.
+ *
+ * \TODO call parameters
+ */
+ void update(const SolutionVector &globalSolution)
+ {
+ };
+
+ /*!
+ * \brief Apply the intrinsic permeability tensor to a pressure
+ * potential gradient.
+ *
+ * \param element The current finite element
+ * \param fvElemGeom The current finite volume geometry of the element
+ * \param scvfIdx The index sub-control volume face where the
+ * intrinsic velocity ought to be calculated.
+ */
+ const Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineK_;
+ return coarseK_;
+ }
+
+ /*!
+ * \brief Define the porosity \f$[-]\f$ of the soil
+ *
+ * \param vDat The data defined on the sub-control volume
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the porosity needs to be defined
+ */
+ double porosity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return finePorosity_;
+ else
+ return coarsePorosity_;
+ }
+
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
+ if (isFineMaterial_(pos))
+ return fineMaterialParams_;
+ else
+ return coarseMaterialParams_;
+ }
+
+ /*!
+ * \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvElemGeom The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume where
+ * the heat capacity needs to be defined
+ */
+ double heatCapacity(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return
+ 790 // specific heat capacity of granite [J / (kg K)]
+ * 2700 // density of granite [kg/m^3]
+ * (1 - porosity(element, fvElemGeom, scvIdx));
+ }
+
+ /*!
+ * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
+ * rock matrix based on the temperature gradient \f$[K / m]\f$
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param heatFlux The result vector
+ * \param tempGrad The temperature gradient
+ * \param element The current finite element
+ * \param fvElemGeom The finite volume geometry of the current element
+ * \param scvfIdx The local index of the sub-control volume face where
+ * the matrix heat flux should be calculated
+ */
+ void matrixHeatFlux(Vector &heatFlux,
+ const FluxData &fluxDat,
+ const VertexDataArray &vDat,
+ const Vector &tempGrad,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvfIdx) const
+ {
+ static const Scalar lWater = 0.6;
+ static const Scalar lGranite = 2.8;
+
+ // arithmetic mean of the liquid saturation and the porosity
+ const int i = fvElemGeom.subContVolFace[scvfIdx].i;
+ const int j = fvElemGeom.subContVolFace[scvfIdx].j;
+ Scalar Sl = std::max(0.0, (vDat[i].saturation(lPhaseIdx) +
+ vDat[j].saturation(lPhaseIdx)) / 2);
+ Scalar poro = (porosity(element, fvElemGeom, i) +
+ porosity(element, fvElemGeom, j)) / 2;
+
+ Scalar lsat = pow(lGranite, (1-poro)) * pow(lWater, poro);
+ Scalar ldry = pow(lGranite, (1-poro));
+
+ // the heat conductivity of the matrix. in general this is a
+ // tensorial value, but we assume isotropic heat conductivity.
+ Scalar heatCond = ldry + sqrt(Sl) * (ldry - lsat);
+
+ // the matrix heat flux is the negative temperature gradient
+ // times the heat conductivity.
+ heatFlux = tempGrad;
+ heatFlux *= -heatCond;
+ }
+
+private:
+ bool isFineMaterial_(const GlobalPosition &pos) const
+ { return pos[dim-1] > layerBottom_; };
+
+ Scalar fineK_;
+ Scalar coarseK_;
+ Scalar layerBottom_;
+
+ Scalar finePorosity_;
+ Scalar coarsePorosity_;
+
+ MaterialLawParams fineMaterialParams_;
+ MaterialLawParams coarseMaterialParams_;
+};
+
+}
+
+#endif
diff --git a/test/boxmodels/2pni/CMakeLists.txt b/test/boxmodels/2pni/CMakeLists.txt
new file mode 100644
index 0000000000..9870500082
--- /dev/null
+++ b/test/boxmodels/2pni/CMakeLists.txt
@@ -0,0 +1,17 @@
+add_definitions(-DYASPGRID -DENABLE_ALBERTA -DENABLE_ALUGRID -DENABLE_UG)
+
+ADD_EXECUTABLE("test_2pni" test_2pni.cc)
+TARGET_LINK_LIBRARIES("test_2pni" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_2pni"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
diff --git a/test/boxmodels/2pni/Makefile.am b/test/boxmodels/2pni/Makefile.am
new file mode 100644
index 0000000000..075f3a5861
--- /dev/null
+++ b/test/boxmodels/2pni/Makefile.am
@@ -0,0 +1,5 @@
+bin_PROGRAMS = test_2pni
+
+test_2pni_SOURCES = test_2pni.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/2pni/injectionproblem2pni.hh b/test/boxmodels/2pni/injectionproblem2pni.hh
new file mode 100644
index 0000000000..9a6807df59
--- /dev/null
+++ b/test/boxmodels/2pni/injectionproblem2pni.hh
@@ -0,0 +1,362 @@
+// $Id: injectionproblem2pni.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Melanie Darcis *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_INJECTIONPROBLEM2PNI_HH
+#define DUMUX_INJECTIONPROBLEM2PNI_HH
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/boxmodels/2pni/2pniboxmodel.hh>
+
+#include <dumux/material/fluidsystems/h2o_n2_system.hh>
+
+// use the same spatial parameters as the injection problem of the
+// 2p2c test program
+#include "../2p2c/injectionspatialparameters.hh"
+
+#define ISOTHERMAL 0
+
+namespace Dumux {
+
+template <class TypeTag>
+class InjectionProblem2PNI;
+
+namespace Properties
+{
+#if !ISOTHERMAL
+NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(BoxTwoPNI));
+#else
+NEW_TYPE_TAG(InjectionProblem2PNI, INHERITS_FROM(BoxTwoP));
+#endif
+
+// Set the grid type
+SET_PROP(InjectionProblem2PNI, Grid)
+{
+#if HAVE_UG
+ typedef Dune::UGGrid<2> type;
+#else
+ typedef Dune::SGrid<2, 2> type;
+ //typedef Dune::YaspGrid<2> type;
+#endif
+};
+
+SET_PROP(InjectionProblem2PNI, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+
+// Set the problem property
+SET_PROP(InjectionProblem2PNI, Problem)
+{
+ typedef Dumux::InjectionProblem2PNI<TypeTag> type;
+};
+
+// Set the soil properties. we use the same spatial parameters as the
+// 2p2c injection problem
+SET_PROP(InjectionProblem2PNI, SpatialParameters)
+{ typedef InjectionSpatialParameters<TypeTag> type; };
+
+#if 1
+// Use the same fluid system as the 2p2c injection problem
+SET_PROP(InjectionProblem2PNI, FluidSystem)
+{
+ typedef H2O_N2_System<TypeTag> type;
+};
+#else
+// Set the wetting phase
+SET_PROP(InjectionProblem2PNI, WettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(InjectionProblem2PNI, NonwettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::GasPhase<Scalar, Dumux::N2<Scalar> > type;
+};
+#endif
+
+// Enable gravity
+SET_BOOL_PROP(InjectionProblem2PNI, EnableGravity, true);
+
+// write convergence behaviour to disk?
+SET_BOOL_PROP(InjectionProblem2PNI, NewtonWriteConvergence, true);
+}
+
+/*!
+ * \ingroup TwoPNIBoxProblems
+ * \brief Nonisothermal gas injection problem where a gas (e.g. air) is injected into a fully
+ * water saturated medium. During buoyancy driven upward migration the gas
+ * passes a high temperature area.
+ *
+ * The domain is sized 40 m times 40 m. The rectangular area with the increased temperature (380 K)
+ * starts at (20 m, 5 m) and ends at (30 m, 35 m)
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top and on the bottom of the domain, while dirichlet conditions
+ * apply on the left and the right boundary.
+ * For the energy conservation equation dirichlet boundary conditions are applied
+ * on all boundaries.
+ *
+ * Gas is injected at the bottom boundary from 15 m to 25 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure, a gas saturation of zero and
+ * a geothermal temperature gradient of 0.03 K/m are applied.
+ *
+ * This problem uses the \ref TwoPNIBoxModel.
+ *
+ * This problem should typically be simulated for 300000 seconds.
+ * A good choice for the initial time step size is 1000 seconds.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_2pni ./grids/test_2pni.dgf 300000 1000</tt>
+ */
+template<class TypeTag>
+class InjectionProblem2PNI
+#if !ISOTHERMAL
+ : public TwoPNIBoxProblem<TypeTag,
+ InjectionProblem2PNI<TypeTag> >
+#else
+ : public TwoPNIBoxProblem<TypeTag,
+ InjectionProblem2PNI<TypeTag> >
+#endif
+{
+ typedef InjectionProblem2PNI<TypeTag> ThisType;
+ typedef TwoPNIBoxProblem<TypeTag, ThisType> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+#if ISOTHERMAL
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+#else
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPNIIndices)) Indices;
+#endif
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+ pressureIdx = Indices::pressureIdx,
+ saturationIdx = Indices::saturationIdx,
+
+ contiWEqIdx = Indices::contiWEqIdx,
+ contiNEqIdx = Indices::contiNEqIdx,
+
+#if !ISOTHERMAL
+ temperatureIdx = Indices::temperatureIdx,
+ energyEqIdx = Indices::energyEqIdx,
+#endif
+
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
+ typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ InjectionProblem2PNI(const GridView &grid)
+ : ParentType(grid)
+ {
+ // initialize the tables of the fluid system
+ FluidSystem::init();
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "injection2pni"; }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ void boundaryTypes(BoundaryTypeVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if (globalPos[0] < eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+
+#if !ISOTHERMAL
+ // set a dirichlet value for the temperature, use the energy
+ // equation to set the value
+ values.setDirichlet(temperatureIdx, energyEqIdx);
+#endif
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ Scalar densityW = 1000.0;
+ values[pressureIdx] = 1e5 + (depthBOR_ - globalPos[1])*densityW*9.81;
+ values[saturationIdx] = 0.0;
+#if !ISOTHERMAL
+ values[temperatureIdx] = 283.0 + (depthBOR_ - globalPos[1])*0.03;
+#endif
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ values = 0;
+ if (globalPos[1] < 15 && globalPos[1] > 5) {
+ // inject air. negative values mean injection
+ values[contiNEqIdx] = -1e-3; // kg/(s*m^2)
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+#if ISOTHERMAL
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 30 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ return 273.15 + 30; // [K]
+ };
+#endif
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate mass
+ * generated or annihilate per volume unit. Positive values mean
+ * that mass is created, negative ones mean that it vanishes.
+ */
+ void source(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &,
+ int subControlVolumeIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ Scalar densityW = 1000.0;
+ values[pressureIdx] = 1e5 + (depthBOR_ - globalPos[1])*densityW*9.81;
+ values[saturationIdx] = 0.0;
+
+#if !ISOTHERMAL
+ values[temperatureIdx] = 283.0 + (depthBOR_ - globalPos[1])*0.03;
+ if (globalPos[0] > 20 && globalPos[0] < 30 && globalPos[1] > 5 && globalPos[1] < 35)
+ values[temperatureIdx] = 380;
+#endif // !ISOTHERMAL
+ }
+ // \}
+
+private:
+ static const Scalar depthBOR_ = 2700.0; // [m]
+ static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/2pni/test_2pni.cc b/test/boxmodels/2pni/test_2pni.cc
new file mode 100644
index 0000000000..d4bf35bcce
--- /dev/null
+++ b/test/boxmodels/2pni/test_2pni.cc
@@ -0,0 +1,26 @@
+// $Id: test_2pni.cc 3779 2010-06-24 07:07:56Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Melanie Darcis *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+#include "injectionproblem2pni.hh"
+#include <dumux/common/start.hh>
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(InjectionProblem2PNI) ProblemTypeTag;
+ return Dumux::startFromDGF<ProblemTypeTag>(argc, argv);
+}
diff --git a/test/boxmodels/CMakeLists.txt b/test/boxmodels/CMakeLists.txt
new file mode 100644
index 0000000000..4b46945fbd
--- /dev/null
+++ b/test/boxmodels/CMakeLists.txt
@@ -0,0 +1,11 @@
+# directories which are commented out do not compile at the moment
+# (they also don't compile if using the old build system)!
+
+add_subdirectory("1p")
+add_subdirectory("1p2c")
+add_subdirectory("2p")
+add_subdirectory("2pni")
+add_subdirectory("2p2c")
+add_subdirectory("2p2cni")
+add_subdirectory("richards")
+
diff --git a/test/boxmodels/Makefile.am b/test/boxmodels/Makefile.am
new file mode 100644
index 0000000000..a6f99f6075
--- /dev/null
+++ b/test/boxmodels/Makefile.am
@@ -0,0 +1,4 @@
+SUBDIRS = 1p 1p2c 2p 2p2c 2p2cia 2p2cni 2pFract 2pia 2pMINC 2pNc 2pni doublecontinuum doublecontinuum_ctTSM linearelasticity minc2p2cni richards .
+
+include $(top_srcdir)/am/global-rules
+
diff --git a/test/boxmodels/richards/CMakeLists.txt b/test/boxmodels/richards/CMakeLists.txt
new file mode 100644
index 0000000000..cc51863d3b
--- /dev/null
+++ b/test/boxmodels/richards/CMakeLists.txt
@@ -0,0 +1,23 @@
+add_definitions(-DYASPGRID -DENABLE_ALBERTA -DENABLE_ALUGRID -DENABLE_UG)
+
+if (MPI_FOUND)
+ # tell UG that the model is parallelized
+ add_definitions(-DModelP)
+endif (MPI_FOUND)
+
+# build target for the twophase-twocomponent test problem
+ADD_EXECUTABLE("test_richards" test_richards.cc)
+TARGET_LINK_LIBRARIES("test_richards" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_richards"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/boxmodels/richards/Makefile.am b/test/boxmodels/richards/Makefile.am
new file mode 100644
index 0000000000..68e7b8d55b
--- /dev/null
+++ b/test/boxmodels/richards/Makefile.am
@@ -0,0 +1,6 @@
+bin_PROGRAMS = test_richards lens_richards
+
+test_richards_SOURCES = test_richards.cc
+lens_richards_SOURCES = lens_richards.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/boxmodels/richards/richardsproblem.hh b/test/boxmodels/richards/richardsproblem.hh
new file mode 100644
index 0000000000..ebe369bf24
--- /dev/null
+++ b/test/boxmodels/richards/richardsproblem.hh
@@ -0,0 +1,293 @@
+// $Id: richardsproblem.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Onur Dogan *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_RICHARDSPROBLEM_HH
+#define DUMUX_RICHARDSPROBLEM_HH
+
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
+#include <dune/grid/io/file/dgfparser/dgfs.hh>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/old_material/fluids/water.hh>
+#include <dumux/old_material/fluids/air.hh>
+
+#include <dumux/boxmodels/richards/richardsboxmodel.hh>
+
+#include "richardssoil.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class RichardsTestProblem;
+
+//////////
+// Specify the properties for the lens problem
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(RichardsTestProblem, INHERITS_FROM(BoxRichards));
+
+// Set the grid type
+// Set the grid type
+#if HAVE_UG
+SET_TYPE_PROP(RichardsTestProblem, Grid, Dune::UGGrid<2>);
+#else
+SET_PROP(RichardsTestProblem, Grid) { typedef Dune::SGrid<2, 2> type; };
+//SET_TYPE_PROP(RichardsTestProblem, Grid, Dune::YaspGrid<2>);
+#endif
+
+#ifdef HAVE_DUNE_PDELAB
+SET_PROP(RichardsTestProblem, LocalFEMSpace)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{dim = GridView::dimension};
+
+public:
+ typedef Dune::PDELab::Q1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for cubes
+// typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
+};
+#endif
+
+// Set the problem property
+SET_PROP(RichardsTestProblem, Problem)
+{
+ typedef Dumux::RichardsTestProblem<TTAG(RichardsTestProblem)> type;
+};
+
+// Set the wetting phase
+SET_TYPE_PROP(RichardsTestProblem, WettingPhase, Dumux::Water);
+
+// Set the non-wetting phase
+SET_TYPE_PROP(RichardsTestProblem, NonwettingPhase, Dumux::Air);
+
+// Set the soil properties
+SET_PROP(RichardsTestProblem, Soil)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::RichardsSoil<Grid, Scalar> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(RichardsTestProblem, EnableGravity, true);
+}
+
+/*!
+ * \ingroup RichardsBoxProblems
+ * \brief Base class for defining an instance of a Richard`s problem, where water is infiltrating into an initially unsaturated zone.
+ *
+ * The domain is box shaped. All sides are closed (Neumann 0 boundary) except the top boundary(Dirichlet),
+ * where water is infiltrating into an initially unsaturated zone. Linear capillary pressure-saturation relationship is used.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_richards ./grids/richards_2d.dgf 10000 1</tt>
+ *
+ * where start simulation time = 1 second, end simulation time = 10000 seconds.
+ *
+ * The same file can be also used for 3d simulation but you need to change line
+ * "typedef Dune::UGGrid<2> type;" with "typedef Dune::UGGrid<3> type;" and use richards_3d.dgf grid
+ */
+template <class TypeTag = TTAG(RichardsTestProblem) >
+class RichardsTestProblem : public RichardsBoxProblem<TypeTag,
+ RichardsTestProblem<TypeTag> >
+{
+ typedef RichardsTestProblem<TypeTag> ThisType;
+ typedef RichardsBoxProblem<TypeTag, ThisType> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(RichardsIndices)) Indices;
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ // copy some indices for convenience
+ pW = Indices::pW,
+
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
+ typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ RichardsTestProblem(const GridView &gridView)
+ : ParentType(gridView)
+ {
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ const char *name() const
+ { return "richardstest"; }
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature(const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx) const
+ {
+ return 283.15; // -> 10°C
+ };
+
+ /*!
+ * \brief Returns the reference pressure of the nonwetting phase.
+ *
+ * This problem assumes a pressure of 1 bar.
+ */
+ Scalar pNreference() const
+ {
+ return 1.0e+5; // reference non-wetting phase pressure [Pa] used for viscosity and density calculations
+ };
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ */
+ void boundaryTypes(BoundaryTypeVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos
+ = fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal;
+
+ if (globalPos[dim-1] > this->bboxMax()[dim-1] - eps)
+ // dirichlet on the lower boundary
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a dirichlet
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ void dirichlet(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+
+ if (globalPos[dim-1] > this->bboxMax()[dim-1] - eps)
+ values[pW] = 1.0e+5*0.99;
+ else
+ values[pW] = 0.0;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ void neumann(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ {
+ // const GlobalPosition &globalPos = fvElemGeom.boundaryFace[boundaryFaceIdx].ipGlobal;
+ values[pW] = 0;
+ }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a values parameter stores the rate
+ * wetting phase mass is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ void source(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ values = Scalar(0.0);
+ }
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * For this method, the \a values parameter stores primary
+ * variables.
+ */
+ void initial(PrimaryVarVector &values,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ {
+ // const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
+ values[pW] = 1.0e+5 - 5.0e+4;
+ }
+ // \}
+
+private:
+ static const Scalar eps = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/boxmodels/richards/richardssoil.hh b/test/boxmodels/richards/richardssoil.hh
new file mode 100644
index 0000000000..da01fdd148
--- /dev/null
+++ b/test/boxmodels/richards/richardssoil.hh
@@ -0,0 +1,94 @@
+// $Id: richardssoil.hh 3784 2010-06-24 13:43:57Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Onur Dogan *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_RICHARDS_SOIL_HH
+#define DUMUX_RICHARDS_SOIL_HH
+
+#include <dumux/old_material/property_baseclasses.hh>
+#include <dumux/old_material/relperm_pc_law.hh>
+#include <dumux/old_material/matrixproperties.hh>
+#include <dumux/old_material/twophaserelations.hh>
+
+namespace Dumux
+{
+
+/**
+ * Defines soil properties Richards problem.
+ */
+template<class Grid, class ScalarT>
+class RichardsSoil : public Matrix2p<Grid, ScalarT>
+{
+public:
+ typedef typename Grid::Traits::template Codim<0>::Entity Element;
+ typedef ScalarT Scalar;
+ typedef typename Grid::ctype CoordScalar;
+ enum {dim=Grid::dimension, dimWorld=Grid::dimensionworld};
+
+ typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+
+ RichardsSoil():Matrix2p<Grid,Scalar>()
+ {
+ Kout_ = 0.;
+ for(int i = 0; i < dim; i++)
+ Kout_[i][i] = 1e-11;
+ }
+
+ ~RichardsSoil()
+ {}
+
+ const Dune::FieldMatrix<CoordScalar,dim,dim> &K (const GlobalPosition &x, const Element& e, const LocalPosition &xi) const
+ {
+ return Kout_;
+ }
+
+ double porosity(const GlobalPosition &x, const Element& e, const LocalPosition &xi) const
+ {
+ return 0.4;
+ }
+
+ double Sr_w(const GlobalPosition &x, const Element& e, const LocalPosition &xi, const double T) const
+ {
+ return 0.05;
+ }
+
+ double Sr_n(const GlobalPosition &x, const Element& e, const LocalPosition &xi, const double T) const
+ {
+ return 0.0;
+ }
+
+ std::vector<double> paramRelPerm(const GlobalPosition &x, const Element& e, const LocalPosition &xi, const double T) const
+ {
+ // example for Brooks-Corey parameters
+ std::vector<double> param(2);
+ param[0] = 0; // pCMin
+ param[1] = 5.0e+4; // pCMax
+
+ return param;
+ }
+
+ typename Matrix2p<Grid,Scalar>::modelFlag relPermFlag(const GlobalPosition &x, const Element& e, const LocalPosition &xi) const
+ {
+ return Matrix2p<Grid,Scalar>::linear;
+ }
+
+private:
+ Dune::FieldMatrix<Scalar,dim,dim> Kout_;
+};
+
+} // end namespace
+
+#endif
diff --git a/test/boxmodels/richards/test_richards.cc b/test/boxmodels/richards/test_richards.cc
new file mode 100644
index 0000000000..8dc39569f1
--- /dev/null
+++ b/test/boxmodels/richards/test_richards.cc
@@ -0,0 +1,93 @@
+// $Id: test_richards.cc 3780 2010-06-24 07:39:50Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2009 by Onur Dogan *
+ * Copyright (C) 2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "richardsproblem.hh"
+#include <dune/common/exceptions.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] gridFile.dgf tEnd dt\n")%progname;
+ exit(1);
+};
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(RichardsTestProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+ typedef Dune::GridPtr<Grid> GridPointer;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ // parse the command line arguments for the program
+ if (argc < 4)
+ usage(argv[0]);
+
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+
+ if (argc - argPos != 3) {
+ usage(argv[0]);
+ }
+
+ double tEnd, dt;
+ const char *dgfFileName = argv[argPos++];
+ std::istringstream(argv[argPos++]) >> tEnd;
+ std::istringstream(argv[argPos++]) >> dt;
+
+ // create grid
+
+ // -> load the grid from file
+ GridPointer gridPtr = GridPointer(dgfFileName);
+ (*gridPtr).loadBalance();
+ Dune::gridinfo(*gridPtr);
+
+ // instantiate and run the concrete problem
+ Problem problem(gridPtr->leafView());
+ if (!problem.simulate(dt, tEnd))
+ return 2;
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ }
+
+ return 3;
+}
diff --git a/test/common/CMakeLists.txt b/test/common/CMakeLists.txt
new file mode 100644
index 0000000000..1639200866
--- /dev/null
+++ b/test/common/CMakeLists.txt
@@ -0,0 +1,6 @@
+# directories which are commented out do not compile at the moment
+# (they also don't compile if using the old build system)!
+
+add_subdirectory("propertysystem")
+#add_subdirectory("spline")
+
diff --git a/test/common/Makefile.am b/test/common/Makefile.am
new file mode 100644
index 0000000000..a381c552af
--- /dev/null
+++ b/test/common/Makefile.am
@@ -0,0 +1,4 @@
+SUBDIRS = pardiso propertysystem spline .
+
+include $(top_srcdir)/am/global-rules
+
diff --git a/test/common/pardiso/CMakeLists.txt b/test/common/pardiso/CMakeLists.txt
new file mode 100644
index 0000000000..cc1c834081
--- /dev/null
+++ b/test/common/pardiso/CMakeLists.txt
@@ -0,0 +1,18 @@
+add_definitions(-DYASPGRID -DGRIDDIM=2 -DENABLE_UG)
+
+# add build targets
+ADD_EXECUTABLE("test_box3d" test_box3d.cc)
+TARGET_LINK_LIBRARIES("test_box3d" ${DumuxLinkLibraries})
+
+# add required libraries and includes to the build flags
+LINK_DIRECTORIES(${DumuxLinkDirectories})
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
+
+# make sure the grids are present in the build directory
+add_custom_command(TARGET "test_box3d"
+ POST_BUILD
+ COMMAND ${CMAKE_COMMAND} -E
+ copy_directory
+ "${CMAKE_CURRENT_SOURCE_DIR}/grids"
+ "${CMAKE_CURRENT_BINARY_DIR}/grids")
+
diff --git a/test/common/pardiso/Makefile.am b/test/common/pardiso/Makefile.am
new file mode 100644
index 0000000000..5d3e07ec30
--- /dev/null
+++ b/test/common/pardiso/Makefile.am
@@ -0,0 +1,18 @@
+# tests where program to build and program to run are equal
+NORMALTESTS = test_pardiso test_matrixroot
+
+# list of tests to run
+#TESTS = $(NORMALTESTS)
+
+# programs just to build when "make check" is used
+check_PROGRAMS = $(NORMALTESTS)
+
+dist_noinst_DATA = test_pardiso.cc
+
+test_pardiso_SOURCES = test_pardiso.cc
+
+test_matrixroot_SOURCES = test_matrixroot.cc
+test_matrixroot_CXXFLAGS=-I/usr/include/gsl
+test_matrixroot_LDADD=-lgsl
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/common/pardiso/test_pardiso.cc b/test/common/pardiso/test_pardiso.cc
new file mode 100644
index 0000000000..880026cf55
--- /dev/null
+++ b/test/common/pardiso/test_pardiso.cc
@@ -0,0 +1,198 @@
+#include "config.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "config.h"
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/io.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include "dumux/common/pardiso.hh"
+
+int main(int argc, char** argv)
+{
+ try
+ {
+#ifdef HAVE_PARDISO
+ /* Matrix data. */
+ int n = 8;
+ int ia[ 9] = { 0, 4, 7, 9, 11, 12, 15, 17, 20 };
+ int ja[20] = { 0, 2, 5, 6,
+ 1, 2, 4,
+ 2, 7,
+ 3, 6,
+ 4,
+ 2, 5, 7,
+ 1, 6,
+ 2, 6, 7 };
+ double a[20] = { 7.0, 1.0, 2.0, 7.0,
+ -4.0, 8.0, 2.0,
+ 1.0, 5.0,
+ 7.0, 9.0,
+ -4.0,
+ 7.0, 3.0, 8.0,
+ 1.0, 11.0,
+ -3.0, 2.0, 5.0 };
+
+ int nnz = ia[n];
+ int mtype = 11; /* Real unsymmetric matrix */
+
+ /* RHS and solution vectors. */
+ double b[8], x[8];
+ int nrhs = 1; /* Number of right hand sides. */
+
+ /* Internal solver memory pointer pt, */
+ /* 32-bit: int pt[64]; 64-bit: long int pt[64] */
+ /* or void *pt[64] should be OK on both architectures */
+ void *pt[64];
+
+ /* Pardiso control parameters. */
+ int iparm[64];
+ int maxfct, mnum, phase, error, msglvl;
+
+ /* Number of processors. */
+ //int num_procs;
+
+ /* Auxiliary variables. */
+ //char *var;
+ int i;
+
+ double ddum; /* Double dummy */
+ int idum; /* Integer dummy. */
+
+ /* -------------------------------------------------------------------- */
+ /* .. Setup Pardiso control parameters. */
+ /* -------------------------------------------------------------------- */
+
+ F77_FUN(pardisoinit) (pt, &mtype, iparm);
+
+ iparm[2] = 1;
+
+ maxfct = 1; /* Maximum number of numerical factorizations. */
+ mnum = 1; /* Which factorization to use. */
+
+ msglvl = 0; /* Print statistical information */
+ error = 0; /* Initialize error flag */
+
+ /* -------------------------------------------------------------------- */
+ /* .. Convert matrix from 0-based C-notation to Fortran 1-based */
+ /* notation. */
+ /* -------------------------------------------------------------------- */
+ for (i = 0; i < n+1; i++) {
+ ia[i] += 1;
+ }
+ for (i = 0; i < nnz; i++) {
+ ja[i] += 1;
+ }
+
+ phase = 13;
+
+ iparm[7] = 1; /* Max numbers of iterative refinement steps. */
+
+ /* Set right hand side to one. */
+ for (i = 0; i < n; i++) {
+ b[i] = 1;
+ }
+
+ F77_FUN(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
+ &n, a, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, b, x, &error);
+
+ if (error != 0) {
+ printf("\nERROR during solution: %d", error);
+ exit(3);
+ }
+
+ printf("\nSolve completed ... ");
+ printf("\nThe solution of the system is: ");
+ for (i = 0; i < n; i++) {
+ printf("\n x [%d] = % f", i, x[i] );
+ }
+ printf ("\n\n");
+
+ /* -------------------------------------------------------------------- */
+ /* .. Termination and release of memory. */
+ /* -------------------------------------------------------------------- */
+ phase = -1; /* Release internal memory. */
+
+ F77_FUN(pardiso) (pt, &maxfct, &mnum, &mtype, &phase,
+ &n, &ddum, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, &ddum, &ddum, &error);
+
+ typedef Dune::FieldMatrix<double,1,1> M;
+ Dune::BCRSMatrix<M> B(8,8,Dune::BCRSMatrix<M>::random);
+
+ // initially set row size for each row
+ B.setrowsize(0,4);
+ B.setrowsize(1,3);
+ B.setrowsize(2,2);
+ B.setrowsize(3,2);
+ B.setrowsize(4,1);
+ B.setrowsize(5,3);
+ B.setrowsize(6,2);
+ B.setrowsize(7,3);
+
+ // finalize row setup phase
+ B.endrowsizes();
+
+ // add column entries to rows
+ B.addindex(0,0); B.addindex(0,2); B.addindex(0,5); B.addindex(0,6);
+ B.addindex(1,1); B.addindex(1,2); B.addindex(1,4);
+ B.addindex(2,2); B.addindex(2,7);
+ B.addindex(3,3); B.addindex(3,6);
+ B.addindex(4,4);
+ B.addindex(5,2); B.addindex(5,5); B.addindex(5,7);
+ B.addindex(6,1); B.addindex(6,6);
+ B.addindex(7,2); B.addindex(7,6); B.addindex(7,7);
+
+ // finalize column setup phase
+ B.endindices();
+
+ // set entries using the random access operator
+ B[0][0] = 7; B[0][2] = 1; B[0][5] = 2; B[0][6] = 7;
+ B[1][1] = -4; B[1][2] = 8; B[1][4] = 2;
+ B[2][2] = 1; B[2][7] = 5;
+ B[3][3] = 7; B[3][6] = 9;
+ B[4][4] = -4;
+ B[5][2] = 7; B[5][5] = 3; B[5][7] = 8;
+ B[6][1] = 1; B[6][6] = 11;
+ B[7][2] = -3; B[7][6] = 2; B[7][7] = 5;
+
+ //printmatrix(std::cout, B, "matrix B", "row", 9, 1);
+
+ typedef Dune::FieldVector<double, 1> VB;
+ typedef Dune::BlockVector<VB> Vector;
+ typedef Dune::BCRSMatrix<M> Matrix;
+ Dune::MatrixAdapter<Matrix,Vector,Vector> op(B); // make linear operator from A
+ Dune::SeqPardiso<Matrix,Vector,Vector> pardiso(B); // preconditioner object
+ Dune::LoopSolver<Vector> loop(op, pardiso, 1E-14, 2, 1); // an inverse operator
+
+ Vector f(n);
+ f = 1;
+ Vector y(n);
+ y = 0;
+ Dune::InverseOperatorResult r;
+ loop.apply(y, f, r);
+
+ std::cout << "\nSolve completed ... ";
+ std::cout << "\nThe solution of the system is: ";
+ for (i = 0; i < n; i++) {
+ std::cout << "\n x [" << i << "] = " << y[i];
+ }
+ std::cout << "\n";
+
+#else
+ DUNE_THROW(Dune::NotImplemented, "no Pardiso library available, reconfigure with correct --with-pardiso options");
+#endif
+
+
+ return 0;
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ }
+}
diff --git a/test/common/propertysystem/CMakeLists.txt b/test/common/propertysystem/CMakeLists.txt
new file mode 100644
index 0000000000..b6b5f391f9
--- /dev/null
+++ b/test/common/propertysystem/CMakeLists.txt
@@ -0,0 +1,5 @@
+# build the test for the property system
+ADD_EXECUTABLE("test_properties" test_properties.cc)
+
+# add required libraries and includes to the build flags
+INCLUDE_DIRECTORIES(${DumuxIncludeDirectories})
diff --git a/test/common/propertysystem/Makefile.am b/test/common/propertysystem/Makefile.am
new file mode 100644
index 0000000000..a9bcf8a215
--- /dev/null
+++ b/test/common/propertysystem/Makefile.am
@@ -0,0 +1,12 @@
+# tests where program to build and program to run are equal
+NORMALTESTS = test_properties
+
+# list of tests to run
+#TESTS = $(NORMALTESTS)
+
+# programs just to build when "make check" is used
+check_PROGRAMS = $(NORMALTESTS)
+
+test_properties_SOURCES = test_properties.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/common/propertysystem/test_propertysystem.cc b/test/common/propertysystem/test_propertysystem.cc
new file mode 100644
index 0000000000..aace3f60e0
--- /dev/null
+++ b/test/common/propertysystem/test_propertysystem.cc
@@ -0,0 +1,183 @@
+// $Id: test_propertysystem.cc 3831 2010-07-14 08:37:48Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: andreas.lauser _at_ iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+
+/*!
+ * \brief This file tests the properties system.
+ *
+ * We define a few type tags and property tags, then we attach values
+ * to (TypeTag, PropertyTag) tuples and finally we use them in the
+ * main function and print some diagnostic messages.
+ */
+
+#include <dumux/common/propertysystem.hh>
+
+namespace Dumux {
+namespace Properties {
+
+///////////////////
+// Define some hierarchy of type tags:
+//
+// CompactCar --- Sedan -_
+// \.
+// +- Pickup ---_
+// / \.
+// Truck ----------------^ \.
+// Tank ---------------------------------- HummerH1
+///////////////////
+NEW_TYPE_TAG(CompactCar);
+NEW_TYPE_TAG(Truck);
+NEW_TYPE_TAG(Tank);
+
+NEW_TYPE_TAG(Sedan, INHERITS_FROM(CompactCar));
+NEW_TYPE_TAG(Pickup, INHERITS_FROM(Sedan, Truck));
+
+NEW_TYPE_TAG(HummerH1, INHERITS_FROM(Sedan, Pickup, Tank));
+
+///////////////////
+// Define the property tags:
+// TopSpeed, NumSeats, CanonCaliber, GasUsage, AutomaticTransmission, Payload
+///////////////////
+NEW_PROP_TAG(TopSpeed); // [km/h]
+NEW_PROP_TAG(NumSeats); // []
+NEW_PROP_TAG(CanonCaliber); // [mm]
+NEW_PROP_TAG(GasUsage); // [l/100km]
+NEW_PROP_TAG(AutomaticTransmission); // true/false
+NEW_PROP_TAG(Payload); // [t]
+
+///////////////////
+// Make the AutomaticTransmission default to false
+SET_PROP_DEFAULT(AutomaticTransmission)
+{ static const bool value = false; };
+
+///////////////////
+// Define some values for the properties on the type tags:
+//
+// (CompactCar, TopSpeed) = GasUsage*35
+// (CompactCar, NumSeats) = 5
+// (CompactCar, GasUsage) = 4
+//
+// (Truck, TopSpeed) = 100
+// (Truck, NumSeats) = 2
+// (Truck, GasUsage) = 12
+// (Truck, Payload) = 35
+//
+// (Tank, TopSpeed) = 60
+// (Tank, GasUsage) = 65
+// (Tank, CanonCaliber) = 120
+//
+// (Sedan, GasUsage) = 7
+// (Sedan, AutomaticTransmission) = true
+//
+// (Pickup, TopSpeed) = 120
+// (Pickup, Payload) = 5
+//
+// (HummmerH1, TopSpeed) = (Pickup, TopSpeed)
+///////////////////
+
+SET_INT_PROP(CompactCar, TopSpeed, GET_PROP_VALUE(TypeTag, PTAG(GasUsage)) * 35);
+SET_INT_PROP(CompactCar, NumSeats, 5);
+SET_INT_PROP(CompactCar, GasUsage, 4);
+
+SET_INT_PROP(Truck, TopSpeed, 100);
+SET_INT_PROP(Truck, NumSeats, 2);
+SET_INT_PROP(Truck, GasUsage, 12);
+SET_INT_PROP(Truck, Payload, 35);
+
+SET_INT_PROP(Tank, TopSpeed, 60);
+SET_INT_PROP(Tank, GasUsage, 65);
+SET_INT_PROP(Tank, CanonCaliber, 120);
+
+SET_INT_PROP(Sedan, GasUsage, 7);
+SET_BOOL_PROP(Sedan, AutomaticTransmission, true);
+
+SET_INT_PROP(Pickup, TopSpeed, 120);
+SET_INT_PROP(Pickup, Payload, 5);
+
+SET_INT_PROP(HummerH1, TopSpeed, GET_PROP_VALUE(TTAG(Pickup), PTAG(TopSpeed)));
+
+///////////////////
+// Unmount the canon from the Hummer
+UNSET_PROP(HummerH1, CanonCaliber);
+
+}; // namespace Properties
+}; // namespace Dumux
+
+
+int main()
+{
+ // print all properties for all type tags
+ std::cout << "---------------------------------------\n";
+ std::cout << "-- Property values\n";
+ std::cout << "---------------------------------------\n";
+
+ std::cout << "---------- Values for CompactCar ----------\n";
+
+ std::cout << "(CompactCar, TopSpeed) = " << GET_PROP_VALUE(TTAG(CompactCar), PTAG(TopSpeed)) << "\n";
+ std::cout << "(CompactCar, NumSeats) = " << GET_PROP_VALUE(TTAG(CompactCar), PTAG(NumSeats)) << "\n";
+ std::cout << "(CompactCar, GasUsage) = " << GET_PROP_VALUE(TTAG(CompactCar), PTAG(GasUsage)) << "\n";
+ std::cout << "(CompactCar, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(CompactCar), PTAG(AutomaticTransmission)) << "\n";
+
+ std::cout << "---------- Values for Truck ----------\n";
+
+ std::cout << "(Truck, TopSpeed) = " << GET_PROP_VALUE(TTAG(Truck), PTAG(TopSpeed)) << "\n";
+ std::cout << "(Truck, NumSeats) = " << GET_PROP_VALUE(TTAG(Truck), PTAG(NumSeats)) << "\n";
+ std::cout << "(Truck, GasUsage) = " << GET_PROP_VALUE(TTAG(Truck), PTAG(GasUsage)) << "\n";
+ std::cout << "(Truck, Payload) = " << GET_PROP_VALUE(TTAG(Truck), PTAG(Payload)) << "\n";
+ std::cout << "(Truck, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Truck), PTAG(AutomaticTransmission)) << "\n";
+
+ std::cout << "---------- Values for Tank ----------\n";
+
+ std::cout << "(Tank, TopSpeed) = " << GET_PROP_VALUE(TTAG(Tank), PTAG(TopSpeed)) << "\n";
+ std::cout << "(Tank, GasUsage) = " << GET_PROP_VALUE(TTAG(Tank), PTAG(GasUsage)) << "\n";
+ std::cout << "(Tank, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Tank), PTAG(AutomaticTransmission)) << "\n";
+ std::cout << "(Tank, CanonCaliber) = " << GET_PROP_VALUE(TTAG(Tank), PTAG(CanonCaliber)) << "\n";
+
+ std::cout << "---------- Values for Sedan ----------\n";
+
+ std::cout << "(Sedan, TopSpeed) = " << GET_PROP_VALUE(TTAG(Sedan), PTAG(TopSpeed)) << "\n";
+ std::cout << "(Sedan, NumSeats) = " << GET_PROP_VALUE(TTAG(Sedan), PTAG(NumSeats)) << "\n";
+ std::cout << "(Sedan, GasUsage) = " << GET_PROP_VALUE(TTAG(Sedan), PTAG(GasUsage)) << "\n";
+ std::cout << "(Sedan, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Sedan), PTAG(AutomaticTransmission)) << "\n";
+
+ std::cout << "---------- Values for Pickup ----------\n";
+ std::cout << "(Pickup, TopSpeed) = " << GET_PROP_VALUE(TTAG(Pickup), PTAG(TopSpeed)) << "\n";
+ std::cout << "(Pickup, NumSeats) = " << GET_PROP_VALUE(TTAG(Pickup), PTAG(NumSeats)) << "\n";
+ std::cout << "(Pickup, GasUsage) = " << GET_PROP_VALUE(TTAG(Pickup), PTAG(GasUsage)) << "\n";
+ std::cout << "(Pickup, Payload) = " << GET_PROP_VALUE(TTAG(Pickup), PTAG(Payload)) << "\n";
+ std::cout << "(Pickup, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(Pickup), PTAG(AutomaticTransmission)) << "\n";
+
+ std::cout << "---------- Values for HummerH1 ----------\n";
+ std::cout << "(HummerH1, TopSpeed) = " << GET_PROP_VALUE(TTAG(HummerH1), PTAG(TopSpeed)) << "\n";
+ std::cout << "(HummerH1, NumSeats) = " << GET_PROP_VALUE(TTAG(HummerH1), PTAG(NumSeats)) << "\n";
+ std::cout << "(HummerH1, GasUsage) = " << GET_PROP_VALUE(TTAG(HummerH1), PTAG(GasUsage)) << "\n";
+ std::cout << "(HummerH1, Payload) = " << GET_PROP_VALUE(TTAG(HummerH1), PTAG(Payload)) << "\n";
+ std::cout << "(HummerH1, AutomaticTransmission) = " << GET_PROP_VALUE(TTAG(HummerH1), PTAG(AutomaticTransmission)) << "\n";
+ // CanonCaliber is explcitly unset for the Hummer -> this would not compile:
+ // std::cout << "(HummerH1, CanonCaliber) = " << GET_PROP_VALUE(TTAG(HummerH1), PTAG(CanonCaliber)) << "\n";
+
+ std::cout << "\n";
+ std::cout << "---------------------------------------\n";
+ std::cout << "-- Diagnostic messages\n";
+ std::cout << "---------------------------------------\n";
+ std::cout << "---- Message for (HummerH1, CanonCaliber) ---\n"
+ << PROP_DIAGNOSTIC(TTAG(HummerH1), PTAG(CanonCaliber));
+ std::cout << "---- Message for (HummerH1, GasUsage) ---\n"
+ << PROP_DIAGNOSTIC(TTAG(HummerH1), PTAG(GasUsage));
+ std::cout << "---- Message for (HummerH1, AutomaticTransmission) ---\n"
+ << PROP_DIAGNOSTIC(TTAG(HummerH1), PTAG(AutomaticTransmission));
+
+ return 0;
+};
diff --git a/test/common/spline/Makefile.am b/test/common/spline/Makefile.am
new file mode 100644
index 0000000000..380ee3411a
--- /dev/null
+++ b/test/common/spline/Makefile.am
@@ -0,0 +1,5 @@
+bin_PROGRAMS = test_spline
+
+test_spline_SOURCES = test_spline.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/common/spline/test_spline.cc b/test/common/spline/test_spline.cc
new file mode 100644
index 0000000000..41ebc4ce3a
--- /dev/null
+++ b/test/common/spline/test_spline.cc
@@ -0,0 +1,69 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2010 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+/*!
+ * \brief This is a program to test the polynomial spline interpolation.
+ *
+ * It just prints some function to stdout. You can look at the result
+ * using the following commands:
+ *
+----------- snip -----------
+./test_spline > spline.csv
+gnuplot
+
+gnuplot> plot "spline.csv" using 1:2 w l ti "Curve", \
+ "spline.csv" using 1:3 w l ti "Derivative", \
+ "spline.csv" using 1:4 w p ti "Monotonical"
+----------- snap -----------
+
+
+*/
+#include "config.h"
+
+#include <dumux/common/spline.hh>
+
+#include <dune/common/fvector.hh>
+
+template <int numSamples>
+void plot(bool reallyPlot)
+{
+ if (!reallyPlot)
+ return;
+
+ const int n = numSamples - 1;
+ typedef Dune::FieldVector<double, numSamples> FV;
+
+ double x_[] = { 0, 5, 7.5, 8.75, 9.375 };
+ double y_[] = { 10, 0, 10, 0, 10 };
+ double m1 = 10;
+ double m2 = -10;
+ FV &xs = *reinterpret_cast<FV*>(x_);
+ FV &ys = *reinterpret_cast<FV*>(y_);
+
+ Dumux::Spline<double, numSamples> sp(xs, ys, m1, m2);
+
+ sp.printCSV(-0.1*(x_[n] - x_[0]) + x_[0],
+ 0.1*(x_[n] - x_[0]) + x_[n],
+ 1000);
+
+ std::cerr << "Spline is monotonic: " << sp.monotonic(x_[0], x_[n]) << "\n";
+};
+
+int main(int argc, char** argv)
+{
+ plot<2>(false);
+ plot<5>(true);
+ return 0;
+}
diff --git a/test/decoupled/1p/Makefile.am b/test/decoupled/1p/Makefile.am
new file mode 100644
index 0000000000..6214399ca7
--- /dev/null
+++ b/test/decoupled/1p/Makefile.am
@@ -0,0 +1,40 @@
+# tests where program to build and program to run are equal
+NORMALTESTS = test_diffusion test_diffusion3d test_mpfa \
+ test_equilibrated test_gradedmesh test_finesubdomainmesh \
+ test_cpgrid test_normals test_dgf
+
+# list of tests to run
+#TESTS = $(NORMALTESTS)
+
+# programs just to build when "make check" is used
+check_PROGRAMS = $(NORMALTESTS)
+
+test_diffusion_SOURCES = test_diffusion.cc
+
+test_mpfa_SOURCES = test_mpfa.cc
+
+test_equilibrated_SOURCES = test_equilibrated.cc
+
+test_gradedmesh_SOURCES = test_gradedmesh.cc
+
+test_finesubdomainmesh_SOURCES = test_finesubdomainmesh.cc
+
+test_dgf_SOURCES = test_dgf.cc
+
+test_diffusion3d_SOURCES = test_diffusion3d.cc
+test_diffusion3d_CXXFLAGS = -I../../../dune-cornerpoint
+test_diffusion3d_LDADD = ../../../dune-cornerpoint/lib/libdunecornerpoint.la \
+ ../../../dune-cornerpoint/grid/preprocess/libpreprocess.la
+
+test_cpgrid_SOURCES = test_cpgrid.cc
+test_cpgrid_CXXFLAGS = -I../../../dune-cornerpoint
+test_cpgrid_LDADD = ../../../dune-cornerpoint/lib/libdunecornerpoint.la \
+ ../../../dune-cornerpoint/grid/preprocess/libpreprocess.la \
+ -lblas -llapack
+
+test_normals_SOURCES = test_normals.cc
+test_normals_CXXFLAGS = -I../../../dune-cornerpoint
+test_normals_LDADD = ../../../dune-cornerpoint/lib/libdunecornerpoint.la \
+ ../../../dune-cornerpoint/grid/preprocess/libpreprocess.la
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/decoupled/1p/test_diffusion.cc b/test/decoupled/1p/test_diffusion.cc
new file mode 100644
index 0000000000..330a7a08bc
--- /dev/null
+++ b/test/decoupled/1p/test_diffusion.cc
@@ -0,0 +1,128 @@
+// $Id: test_diffusion.cc 3780 2010-06-24 07:39:50Z bernd $
+/*****************************************************************************
+ * Copyright (C) 20010 by Markus Wolff *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+#include <iostream>
+#include <boost/format.hpp>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/mpihelper.hh>
+#include <dune/grid/common/gridinfo.hh>
+
+#include "test_diffusion_problem.hh"
+#include "benchmarkresult.hh"
+
+////////////////////////
+// the main function
+////////////////////////
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s #refine [delta]\n")%progname;
+ exit(1);
+}
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(DiffusionTestProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ static const int dim = Grid::dimension;
+ typedef Dune::FieldVector<Scalar, dim> GlobalPosition;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+ if (argc != 2 && argc != 3)
+ usage(argv[0]);
+
+ int numRefine;
+ std::istringstream(argv[1]) >> numRefine;
+
+ double delta = 1e-3;
+ if (argc == 3)
+ std::istringstream(argv[2]) >> delta;
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ Dune::FieldVector<int,dim> N(1);
+ GlobalPosition L(0.0);
+ GlobalPosition H(1.0);
+ Grid grid(N,L,H);
+ grid.globalRefine(numRefine);
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+ Dune::Timer timer;
+ bool consecutiveNumbering = true;
+
+ typedef GET_PROP_TYPE(TTAG(FVVelocity2PTestProblem), PTAG(Problem)) FVProblem;
+ FVProblem fvProblem(grid.leafView(), delta);
+ timer.reset();
+ fvProblem.init();
+ fvProblem.model().calculateVelocity();
+ double fvTime = timer.elapsed();
+ Dumux::ResultEvaluation fvResult;
+ fvResult.evaluate(grid.leafView(), fvProblem, fvProblem.variables().pressure(), fvProblem.variables().velocity(), consecutiveNumbering);
+
+ typedef GET_PROP_TYPE(TTAG(FVMPFAOVelocity2PTestProblem), PTAG(Problem)) MPFAOProblem;
+ MPFAOProblem mpfaProblem(grid.leafView(), delta);
+ timer.reset();
+ mpfaProblem.init();
+ mpfaProblem.model().calculateVelocity();
+ double mpfaTime = timer.elapsed();
+ Dumux::ResultEvaluation mpfaResult;
+ mpfaResult.evaluate(grid.leafView(), mpfaProblem, mpfaProblem.variables().pressure(), mpfaProblem.variables().velocity(), consecutiveNumbering);
+
+ typedef GET_PROP_TYPE(TTAG(MimeticPressure2PTestProblem), PTAG(Problem)) MimeticProblem;
+ MimeticProblem mimeticProblem(grid.leafView(), delta);
+ timer.reset();
+ mimeticProblem.init();
+ mimeticProblem.model().calculateVelocity();
+ double mimeticTime = timer.elapsed();
+ Dumux::ResultEvaluation mimeticResult;
+ mimeticResult.evaluate(grid.leafView(), mimeticProblem, mimeticProblem.variables().pressure(), mimeticProblem.variables().velocity(), consecutiveNumbering);
+
+ std::cout.setf(std::ios_base::scientific, std::ios_base::floatfield);
+ std::cout.precision(2);
+ std::cout << "\t error press \t error grad\t sumflux\t erflm\t\t uMin\t\t uMax\t\t time" << std::endl;
+ std::cout << "2pfa\t " << fvResult.relativeL2Error << "\t " << fvResult.ergrad << "\t " << fvResult.sumflux
+ << "\t " << fvResult.erflm << "\t " << fvResult.uMin << "\t " << fvResult.uMax << "\t " << fvTime << std::endl;
+ std::cout << "mpfa-o\t " << mpfaResult.relativeL2Error << "\t " << mpfaResult.ergrad << "\t " << mpfaResult.sumflux
+ << "\t " << mpfaResult.erflm << "\t " << mpfaResult.uMin << "\t " << mpfaResult.uMax << "\t " << mpfaTime << std::endl;
+ std::cout << "mimetic\t " << mimeticResult.relativeL2Error << "\t " << mimeticResult.ergrad << "\t " << mimeticResult.sumflux
+ << "\t " << mimeticResult.erflm << "\t " << mimeticResult.uMin << "\t " << mimeticResult.uMax << "\t " << mimeticTime << std::endl;
+
+
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ throw;
+ }
+
+ return 3;
+}
diff --git a/test/decoupled/1p/test_diffusion_problem.hh b/test/decoupled/1p/test_diffusion_problem.hh
new file mode 100644
index 0000000000..bf6c9842b1
--- /dev/null
+++ b/test/decoupled/1p/test_diffusion_problem.hh
@@ -0,0 +1,294 @@
+// $Id: test_diffusion_problem.hh 3655 2010-05-26 17:13:50Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_TEST_2P_PROBLEM_HH
+#define DUMUX_TEST_2P_PROBLEM_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+
+#include <dumux/material/components/unit.hh>
+
+#include <dumux/decoupled/2p/diffusion/fvmpfa/mpfaproperties.hh>
+#include <dumux/decoupled/2p/diffusion/diffusionproblem2p.hh>
+#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
+//#include <dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocities2p_upwind.hh>
+#include <dumux/decoupled/2p/diffusion/mimetic/mimeticpressure2p.hh>
+
+#include "test_diffusion_spatialparams.hh"
+
+namespace Dumux
+{
+
+struct FileNames
+{
+ enum
+ {
+ TPFAName = 0, MPFAName = 1, MimeticName = 2
+ };
+};
+
+template<class TypeTag>
+class TestDiffusionProblem;
+
+//////////
+// Specify the properties
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(DiffusionTestProblem, INHERITS_FROM(DecoupledTwoP, MPFAProperties))
+;
+
+// Set the grid type
+SET_PROP(DiffusionTestProblem, Grid)
+{// typedef Dune::YaspGrid<2> type;
+typedef Dune::SGrid<2, 2> type;
+};
+
+SET_PROP(DiffusionTestProblem, PressurePreconditioner)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureCoefficientMatrix)) Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PressureRHSVector)) Vector;
+public:
+ typedef Dune::SeqPardiso<Matrix, Vector, Vector> type;
+};
+
+//SET_INT_PROP(DiffusionTestProblem, VelocityFormulation,
+// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW);
+
+SET_INT_PROP(DiffusionTestProblem, PressureFormulation,
+ GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureGlobal);
+
+// Set the wetting phase
+SET_PROP(DiffusionTestProblem, WettingPhase)
+{
+private:
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(DiffusionTestProblem, NonwettingPhase)
+{
+private:
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+typedef Dumux::LiquidPhase<Scalar, Dumux::Unit<Scalar> > type;
+};
+
+// Set the soil properties
+SET_PROP(DiffusionTestProblem, SpatialParameters)
+{
+private:
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+typedef Dumux::TestDiffusionSpatialParams<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(DiffusionTestProblem, EnableGravity, false);
+
+NEW_PROP_TAG( FileName );
+
+// set the types for the 2PFA FV method
+NEW_TYPE_TAG(FVVelocity2PTestProblem, INHERITS_FROM(DiffusionTestProblem));
+SET_INT_PROP(FVVelocity2PTestProblem, FileName, FileNames::TPFAName);
+SET_INT_PROP(FVVelocity2PTestProblem, IterationNumberPreconditioner, 0);
+SET_TYPE_PROP(FVVelocity2PTestProblem, Model, Dumux::FVVelocity2P<TTAG(FVVelocity2PTestProblem)>);
+SET_TYPE_PROP(FVVelocity2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(FVVelocity2PTestProblem)>);
+
+// set the types for the MPFA-O FV method
+NEW_TYPE_TAG(FVMPFAOVelocity2PTestProblem, INHERITS_FROM(DiffusionTestProblem));
+SET_INT_PROP(FVMPFAOVelocity2PTestProblem, FileName, FileNames::MPFAName);
+SET_INT_PROP(FVMPFAOVelocity2PTestProblem, IterationNumberPreconditioner, 1);
+//SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, Model, Dumux::FVMPFAOVelocity2P<TypeTag>);
+SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, Model, Dumux::FVMPFAOVelocities2P<TypeTag>);
+SET_TYPE_PROP(FVMPFAOVelocity2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(FVMPFAOVelocity2PTestProblem)>);
+
+// set the types for the mimetic FD method
+NEW_TYPE_TAG(MimeticPressure2PTestProblem, INHERITS_FROM(DiffusionTestProblem));
+SET_INT_PROP(MimeticPressure2PTestProblem, FileName, FileNames::MimeticName);
+SET_INT_PROP(MimeticPressure2PTestProblem, IterationNumberPreconditioner, 1);
+SET_TYPE_PROP(MimeticPressure2PTestProblem, Model, Dumux::MimeticPressure2P<TTAG(MimeticPressure2PTestProblem)>);
+SET_TYPE_PROP(MimeticPressure2PTestProblem, Problem, Dumux::TestDiffusionProblem<TTAG(MimeticPressure2PTestProblem)>);
+}
+
+/*!
+ * \ingroup DecoupledProblems
+ */
+template<class TypeTag = TTAG(DiffusionTestProblem)>
+class TestDiffusionProblem: public DiffusionProblem2P<TypeTag, TestDiffusionProblem<TypeTag> >
+{
+typedef TestDiffusionProblem<TypeTag> ThisType;
+typedef DiffusionProblem2P<TypeTag, ThisType> ParentType;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+enum
+{
+dim = GridView::dimension, dimWorld = GridView::dimensionworld
+};
+
+enum
+{
+wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+};
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+typedef typename GridView::Traits::template Codim<0>::Entity Element;
+typedef typename GridView::Intersection Intersection;
+typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+
+public:
+TestDiffusionProblem(const GridView &gridView, const double delta = 1.0) :
+ParentType(gridView), delta_(delta)
+{
+this->variables().saturation() = 1.0;
+this->spatialParameters().setDelta(delta_);
+}
+
+/*!
+ * \name Problem parameters
+ */
+// \{
+
+/*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+const char *name() const
+{
+switch (GET_PROP_VALUE(TypeTag, PTAG(FileName)))
+{
+ case FileNames::TPFAName:
+ return "fvdiffusion";
+ case FileNames::MPFAName:
+ return "fvmpfaodiffusion";
+ case FileNames::MimeticName:
+ return "mimeticdiffusion";
+}
+return "test_diffusion";
+}
+
+bool shouldWriteRestartFile() const
+{
+return false;
+}
+
+/*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+{
+return 273.15 + 10; // -> 10°C
+}
+
+// \}
+
+Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+{
+return 1e5; // -> 10°C
+}
+
+std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element)
+{
+double pi = 4.0*atan(1.0);
+double rt = globalPos[0]*globalPos[0]+globalPos[1]*globalPos[1];
+double ux = pi*cos(pi*globalPos[0])*sin(pi*globalPos[1]);
+double uy = pi*cos(pi*globalPos[1])*sin(pi*globalPos[0]);
+double kxx = (delta_*globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])/rt;
+double kxy = -(1.0 - delta_)*globalPos[0]*globalPos[1]/rt;
+double kyy = (globalPos[0]*globalPos[0] + delta_*globalPos[1]*globalPos[1])/rt;
+double f0 = sin(pi*globalPos[0])*sin(pi*globalPos[1])*pi*pi*(1.0 + delta_)*(globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])
++ cos(pi*globalPos[0])*sin(pi*globalPos[1])*pi*(1.0 - 3.0*delta_)*globalPos[0]
++ cos(pi*globalPos[1])*sin(pi*globalPos[0])*pi*(1.0 - 3.0*delta_)*globalPos[1]
++ cos(pi*globalPos[1])*cos(pi*globalPos[0])*2.0*pi*pi*(1.0 - delta_)*globalPos[0]*globalPos[1];
+
+std::vector<double> result(2, 0.0);
+result[wPhaseIdx]=(f0 + 2.0*(globalPos[0]*(kxx*ux + kxy*uy) + globalPos[1]*(kxy*ux + kyy*uy)))/rt;
+
+return (result);
+}
+
+typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+return BoundaryConditions::dirichlet;
+}
+
+Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+return (exact(globalPos));
+}
+
+typename BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+return BoundaryConditions::dirichlet;
+}
+
+Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+return 1.0;
+}
+
+std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+std::vector<Scalar> neumannFlux(2, 0.0);
+
+return neumannFlux;
+}
+
+Scalar exact (const GlobalPosition& globalPos) const
+{
+double pi = 4.0*atan(1.0);
+
+return (sin(pi*globalPos[0])*sin(pi*globalPos[1]));
+}
+
+Dune::FieldVector<Scalar,dim> exactGrad (const GlobalPosition& globalPos) const
+{
+Dune::FieldVector<Scalar,dim> grad(0);
+double pi = 4.0*atan(1.0);
+grad[0] = pi*cos(pi*globalPos[0])*sin(pi*globalPos[1]);
+grad[1] = pi*cos(pi*globalPos[1])*sin(pi*globalPos[0]);
+
+return grad;
+}
+
+private:
+double delta_;
+};
+} //end namespace
+
+#endif
diff --git a/test/decoupled/1p/test_diffusion_spatialparams.hh b/test/decoupled/1p/test_diffusion_spatialparams.hh
new file mode 100644
index 0000000000..f532c83535
--- /dev/null
+++ b/test/decoupled/1p/test_diffusion_spatialparams.hh
@@ -0,0 +1,100 @@
+// $Id: test_2p_spatialparams.hh 3456 2010-04-09 12:11:51Z mwolff $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef TEST_DIFFUSION_SPATIALPARAMETERS_HH
+#define TEST_DIFFUSION_SPATIALPARAMETERS_HH
+
+
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class TestDiffusionSpatialParams
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+
+ enum
+ {dim=Grid::dimension, dimWorld=Grid::dimensionworld, numEq=1};
+ typedef typename Grid::Traits::template Codim<0>::Entity Element;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
+ typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
+
+ typedef LinearMaterial<Scalar> RawMaterialLaw;
+public:
+ typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ void update (Scalar saturationW, const Element& element)
+ {
+
+ }
+
+ const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const
+ {
+ double rt = globalPos[0]*globalPos[0]+globalPos[1]*globalPos[1];
+ permeability_[0][0] = (delta_*globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1])/rt;
+ permeability_[0][1] = permeability_[1][0] = -(1.0 - delta_)*globalPos[0]*globalPos[1]/rt;
+ permeability_[1][1] = (globalPos[0]*globalPos[0] + delta_*globalPos[1]*globalPos[1])/rt;
+
+ return permeability_;
+ }
+
+ double porosity(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 0.2;
+ }
+
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
+ {
+ return materialLawParams_;
+ }
+
+ void setDelta(const double delta)
+ {
+ delta_ = delta;
+ }
+
+ TestDiffusionSpatialParams(const GridView& gridView)
+ : permeability_(0)
+ {
+ // residual saturations
+ materialLawParams_.setSwr(0.0);
+ materialLawParams_.setSnr(0.0);
+
+ // parameters for the linear entry pressure function
+ materialLawParams_.setEntryPC(0);
+ materialLawParams_.setMaxPC(0);
+ }
+
+private:
+ MaterialLawParams materialLawParams_;
+ mutable FieldMatrix permeability_;
+ double delta_;
+};
+
+} // end namespace
+#endif
diff --git a/test/decoupled/2p/Makefile.am b/test/decoupled/2p/Makefile.am
new file mode 100644
index 0000000000..8678110a1f
--- /dev/null
+++ b/test/decoupled/2p/Makefile.am
@@ -0,0 +1,7 @@
+bin_PROGRAMS = test_2p test_2pinjection
+
+test_2p_SOURCES = test_2p.cc
+
+test_2pinjection_SOURCES = test_2pinjection.cc
+
+include $(top_srcdir)/am/global-rules
diff --git a/test/decoupled/2p/test_2p.cc b/test/decoupled/2p/test_2p.cc
new file mode 100644
index 0000000000..070a70f509
--- /dev/null
+++ b/test/decoupled/2p/test_2p.cc
@@ -0,0 +1,113 @@
+// $Id: test_2p.cc 3732 2010-06-11 13:27:20Z bernd $
+/*****************************************************************************
+ * Copyright (C) 20010 by Markus Wolff *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#include "config.h"
+
+#include "test_2p_problem.hh"
+
+#include <dune/grid/common/gridinfo.hh>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/mpihelper.hh>
+
+#include <iostream>
+#include <boost/format.hpp>
+
+
+////////////////////////
+// the main function
+////////////////////////
+void usage(const char *progname)
+{
+ std::cout << boost::format("usage: %s [--restart restartTime] tEnd\n")%progname;
+ exit(1);
+}
+
+int main(int argc, char** argv)
+{
+ try {
+ typedef TTAG(TwoPTestProblem) TypeTag;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef Dune::FieldVector<Scalar, Grid::dimensionworld> GlobalPosition;
+
+ static const int dim = Grid::dimension;
+
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper::instance(argc, argv);
+
+ ////////////////////////////////////////////////////////////
+ // parse the command line arguments
+ ////////////////////////////////////////////////////////////
+ if (argc < 2)
+ usage(argv[0]);
+
+ // deal with the restart stuff
+ int argPos = 1;
+ bool restart = false;
+ double restartTime = 0;
+ if (std::string("--restart") == argv[argPos]) {
+ restart = true;
+ ++argPos;
+
+ std::istringstream(argv[argPos++]) >> restartTime;
+ }
+
+ if (argc - argPos != 1) {
+ usage(argv[0]);
+ }
+
+ // read the initial time step and the end time
+ double tEnd, dt;
+ std::istringstream(argv[argPos++]) >> tEnd;
+ dt = tEnd;
+
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ Dune::FieldVector<int,dim> N(6); N[0] = 30;
+ Dune::FieldVector<double ,dim> L(0);
+ Dune::FieldVector<double,dim> H(60); H[0] = 300;
+ Grid grid(N,L,H);
+
+ ////////////////////////////////////////////////////////////
+ // instantiate and run the concrete problem
+ ////////////////////////////////////////////////////////////
+
+ Problem problem(grid.leafView(), L, H);
+
+ // load restart file if necessarry
+ if (restart)
+ problem.deserialize(restartTime);
+
+ // run the simulation
+ if (!problem.simulate(dt, tEnd))
+ return 2;
+
+ return 0;
+ }
+ catch (Dune::Exception &e) {
+ std::cerr << "Dune reported error: " << e << std::endl;
+ }
+ catch (...) {
+ std::cerr << "Unknown exception thrown!\n";
+ throw;
+ }
+
+ return 3;
+}
diff --git a/test/decoupled/2p/test_2p_problem.hh b/test/decoupled/2p/test_2p_problem.hh
new file mode 100644
index 0000000000..7574bc341b
--- /dev/null
+++ b/test/decoupled/2p/test_2p_problem.hh
@@ -0,0 +1,282 @@
+// $Id: test_2p_problem.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2007-2008 by Klaus Mosthaf *
+ * Copyright (C) 2007-2008 by Bernd Flemisch *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef DUMUX_TEST_2P_PROBLEM_HH
+#define DUMUX_TEST_2P_PROBLEM_HH
+
+#if HAVE_UG
+#include <dune/grid/uggrid.hh>
+#endif
+
+#include <dune/grid/yaspgrid.hh>
+#include <dune/grid/sgrid.hh>
+
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/oil.hh>
+
+#include <dumux/decoupled/2p/impes/impesproblem2p.hh>
+#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
+#include <dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh>
+#include <dumux/decoupled/2p/transport/fv/fvsaturation2p.hh>
+#include <dumux/decoupled/2p/transport/fv/capillarydiffusion.hh>
+#include <dumux/decoupled/2p/transport/fv/gravitypart.hh>
+
+#include "test_2p_spatialparams.hh"
+
+namespace Dumux
+{
+
+template<class TypeTag>
+class Test2PProblem;
+
+//////////
+// Specify the properties
+//////////
+namespace Properties
+{
+NEW_TYPE_TAG(TwoPTestProblem, INHERITS_FROM(DecoupledTwoP, MPFAProperties, Transport));
+
+// Set the grid type
+SET_PROP(TwoPTestProblem, Grid)
+{
+ // typedef Dune::YaspGrid<2> type;
+ typedef Dune::SGrid<2, 2> type;
+};
+
+// Set the problem property
+SET_PROP(TwoPTestProblem, Problem)
+{
+public:
+ typedef Dumux::Test2PProblem<TTAG(TwoPTestProblem)> type;
+};
+
+// Set the model properties
+SET_PROP(TwoPTestProblem, SaturationModel)
+{
+ typedef Dumux::FVSaturation2P<TTAG(TwoPTestProblem)> type;
+};
+SET_TYPE_PROP(TwoPTestProblem, DiffusivePart, Dumux::CapillaryDiffusion<TypeTag>);
+SET_TYPE_PROP(TwoPTestProblem, ConvectivePart, Dumux::GravityPart<TypeTag>);
+
+SET_PROP(TwoPTestProblem, PressureModel)
+{
+ typedef Dumux::FVVelocity2P<TTAG(TwoPTestProblem)> type;
+// typedef Dumux::FVMPFAOVelocity2P<TTAG(TwoPTestProblem)> type;
+};
+
+//SET_INT_PROP(TwoPTestProblem, VelocityFormulation,
+// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW);
+
+//SET_INT_PROP(TwoPTestProblem, PressureFormulation,
+// GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::pressureGlobal);
+
+// Set the wetting phase
+SET_PROP(TwoPTestProblem, WettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
+};
+
+// Set the non-wetting phase
+SET_PROP(TwoPTestProblem, NonwettingPhase)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::SimpleH2O<Scalar> > type;
+};
+
+// Set the soil properties
+SET_PROP(TwoPTestProblem, SpatialParameters)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+public:
+ typedef Dumux::Test2PSpatialParams<TypeTag> type;
+};
+
+// Enable gravity
+SET_BOOL_PROP(TwoPTestProblem, EnableGravity, false);
+
+SET_SCALAR_PROP(TwoPTestProblem, CFLFactor, 0.95);
+}
+
+/*!
+ * \ingroup DecoupledProblems
+ */
+template<class TypeTag = TTAG(TwoPTestProblem)>
+class Test2PProblem: public IMPESProblem2P<TypeTag, Test2PProblem<TypeTag> >
+{
+typedef Test2PProblem<TypeTag> ThisType;
+typedef IMPESProblem2P<TypeTag, ThisType> ParentType;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+enum
+{
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+};
+
+enum
+{
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+};
+
+typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+typedef typename GridView::Traits::template Codim<0>::Entity Element;
+typedef typename GridView::Intersection Intersection;
+typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+
+public:
+Test2PProblem(const GridView &gridView, const GlobalPosition lowerLeft = 0, const GlobalPosition upperRight = 0) :
+ParentType(gridView), lowerLeft_(lowerLeft), upperRight_(upperRight)
+{
+}
+
+/*!
+ * \name Problem parameters
+ */
+// \{
+
+/*!
+ * \brief The problem name.
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+const char *name() const
+{
+ return "test2p";
+}
+
+bool shouldWriteRestartFile() const
+{
+ return false;
+}
+
+/*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+Scalar temperature(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 273.15 + 10; // -> 10°C
+}
+
+// \}
+
+
+Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 1e5; // -> 10°C
+}
+
+std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element)
+{
+ return std::vector<Scalar>(2, 0.0);
+}
+
+typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ if ((globalPos[0] < eps_))
+ return BoundaryConditions::dirichlet;
+ // all other boundaries
+ return BoundaryConditions::neumann;
+}
+
+BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ if (globalPos[0] > (upperRight_[0] - eps_) || globalPos[0] < eps_)
+// if (globalPos[0] < eps_)
+ return Dumux::BoundaryConditions::dirichlet;
+ else
+ return Dumux::BoundaryConditions::neumann;
+}
+
+Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ if (globalPos[0] < eps_)
+ {
+ if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)))
+ {
+ const Element& element = *(intersection.inside());
+
+ Scalar pRef = referencePressure(globalPos, element);
+ Scalar temp = temperature(globalPos, element);
+ Scalar sat = this->variables().satElement(element);
+
+ FluidState fluidState;
+ fluidState.update(sat, pRef, pRef, temp);
+ return (2e5 + (upperRight_[dim-1] - globalPos[dim-1]) * FluidSystem::phaseDensity(wPhaseIdx, temp, pRef, fluidState) * this->gravity().two_norm());
+ }
+ else
+ return 2e5;
+ }
+ // all other boundaries
+ return 2e5;
+}
+
+Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ if (globalPos[0] < eps_)
+ return 0.8;
+ // all other boundaries
+ return 0.2;
+}
+
+std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+{
+ std::vector<Scalar> neumannFlux(2, 0.0);
+ if (globalPos[0] > upperRight_[0] - eps_)
+ {
+ neumannFlux[nPhaseIdx] = 3e-4;
+ }
+ return neumannFlux;
+}
+
+Scalar neumannSat(const GlobalPosition& globalPos, const Intersection& intersection, Scalar factor) const
+{
+ if (globalPos[0] > upperRight_[0] - eps_)
+ return factor;
+ return 0;
+}
+
+Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
+{
+ return 0.2;
+}
+
+private:
+GlobalPosition lowerLeft_;
+GlobalPosition upperRight_;
+
+static const Scalar eps_ = 1e-6;
+};
+} //end namespace
+
+#endif
diff --git a/test/decoupled/2p/test_2p_spatialparams.hh b/test/decoupled/2p/test_2p_spatialparams.hh
new file mode 100644
index 0000000000..effd44dd19
--- /dev/null
+++ b/test/decoupled/2p/test_2p_spatialparams.hh
@@ -0,0 +1,108 @@
+// $Id: test_2p_spatialparams.hh 3783 2010-06-24 11:33:53Z bernd $
+/*****************************************************************************
+ * Copyright (C) 2008-2009 by Markus Wolff *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * 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, as long as this copyright notice *
+ * is included in its original form. *
+ * *
+ * This program is distributed WITHOUT ANY WARRANTY. *
+ *****************************************************************************/
+#ifndef TEST_2P_SPATIALPARAMETERS_HH
+#define TEST_2P_SPATIALPARAMETERS_HH
+
+
+//#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux
+{
+
+/** \todo Please doc me! */
+
+template<class TypeTag>
+class Test2PSpatialParams
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+
+ enum
+ {dim=Grid::dimension, dimWorld=Grid::dimensionworld, numEq=1};
+ typedef typename Grid::Traits::template Codim<0>::Entity Element;
+
+ typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldVector<CoordScalar, dim> LocalPosition;
+ typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
+
+// typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw;
+ typedef LinearMaterial<Scalar> RawMaterialLaw;
+public:
+ typedef EffToAbsLaw<RawMaterialLaw> MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ void update (Scalar saturationW, const Element& element)
+ {
+
+ }
+
+ const FieldMatrix& intrinsicPermeability (const GlobalPosition& globalPos, const Element& element) const
+ {
+ return constPermeability_;
+ }
+
+ double porosity(const GlobalPosition& globalPos, const Element& element) const
+ {
+ return 0.2;
+ }
+
+
+ // return the brooks-corey context depending on the position
+ const MaterialLawParams& materialLawParams(const GlobalPosition& globalPos, const Element &element) const
+ {
+ return materialLawParams_;
+ }
+
+
+ Test2PSpatialParams(const GridView& gridView)
+ : constPermeability_(0)
+ {
+ // residual saturations
+ materialLawParams_.setSwr(0.2);
+ materialLawParams_.setSnr(0.2);
+
+// // parameters for the Brooks-Corey Law
+// // entry pressures
+// materialLawParams_.setPe(0);
+//
+// // Brooks-Corey shape parameters
+// materialLawParams_.setAlpha(2);
+
+ // parameters for the linear
+ // entry pressures function
+ materialLawParams_.setEntryPC(0);
+ materialLawParams_.setMaxPC(0);
+
+ for(int i = 0; i < dim; i++)
+ {
+ constPermeability_[i][i] = 1e-7;
+ }
+
+ }
+
+private:
+ MaterialLawParams materialLawParams_;
+ FieldMatrix constPermeability_;
+
+};
+
+} // end namespace
+#endif
diff --git a/test/decoupled/Makefile.am b/test/decoupled/Makefile.am
new file mode 100644
index 0000000000..20af332225
--- /dev/null
+++ b/test/decoupled/Makefile.am
@@ -0,0 +1,7 @@
+SUBDIRS = . \
+ 2p \
+ 2p2c \
+ 2p2cni
+
+include $(top_srcdir)/am/global-rules
+
--
GitLab