diff --git a/test/porousmediumflow/1p/implicit/1ptestproblem.hh b/test/porousmediumflow/1p/implicit/1ptestproblem.hh
index 7a91e85ebc34ac911424aadba0558d7dd8db0f19..9ad2ad894ae8546dc16bc6508f2016df6db00317 100644
--- a/test/porousmediumflow/1p/implicit/1ptestproblem.hh
+++ b/test/porousmediumflow/1p/implicit/1ptestproblem.hh
@@ -41,7 +41,7 @@ class OnePTestProblem;
namespace Properties
{
-NEW_TYPE_TAG(OnePTestProblem, INHERITS_FROM(OneP));
+NEW_TYPE_TAG(OnePTestProblem, INHERITS_FROM(OneP, OnePTestSpatialParams));
NEW_TYPE_TAG(OnePTestBoxProblem, INHERITS_FROM(BoxModel, OnePTestProblem));
NEW_TYPE_TAG(OnePTestCCProblem, INHERITS_FROM(CCModel, OnePTestProblem));
diff --git a/test/porousmediumflow/1p/implicit/1ptestspatialparams.hh b/test/porousmediumflow/1p/implicit/1ptestspatialparams.hh
index 127a4f02f38c8e0437226c51a2aeac564b3e6a36..14af3996deb05b2d1daccae6dfa346bff9a0b131 100644
--- a/test/porousmediumflow/1p/implicit/1ptestspatialparams.hh
+++ b/test/porousmediumflow/1p/implicit/1ptestspatialparams.hh
@@ -26,10 +26,25 @@
#define DUMUX_1P_TEST_SPATIALPARAMS_HH
#include <dumux/material/spatialparams/implicit1p.hh>
+#include <dumux/material/spatialparams/gstatrandomfield.hh>
namespace Dumux
{
+//forward declaration
+template<class TypeTag>
+class OnePTestSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(OnePTestSpatialParams);
+
+// Set properties of the porous medium
+NEW_PROP_TAG(SpatialParamsRandomField);
+SET_BOOL_PROP(OnePTestSpatialParams, SpatialParamsRandomField, false);
+}
+
/*!
* \ingroup OnePModel
* \ingroup ImplicitTestProblems
@@ -43,7 +58,9 @@ class OnePTestSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
typedef ImplicitSpatialParamsOneP<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef std::vector<Scalar> ScalarVector;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::IndexSet IndexSet;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
enum {
@@ -53,16 +70,23 @@ class OnePTestSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
public:
OnePTestSpatialParams(const GridView& gridView)
- : ParentType(gridView)
+ : ParentType(gridView),
+ randomPermeability_(gridView.size(dim), 0.0),
+ indexSet_(gridView.indexSet())
{
- permeability_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.Permeability);
+ randomField_ = GET_PARAM_FROM_GROUP(TypeTag, bool, SpatialParams, RandomField);
+ permeability_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.Permeability);
+ if(!randomField_)
permeabilityLens_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.PermeabilityLens);
+ else
+ initRandomField(gridView);
- lensLowerLeft_ = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, SpatialParams.LensLowerLeft);
- lensUpperRight_ = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, SpatialParams.LensUpperRight);
+ lensLowerLeft_ = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, SpatialParams.LensLowerLeft);
+ lensUpperRight_ = GET_RUNTIME_PARAM(TypeTag, GlobalPosition, SpatialParams.LensUpperRight);
}
/*!
@@ -80,7 +104,12 @@ public:
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
if (isInLens_(globalPos))
- return permeabilityLens_;
+ {
+ if(randomField_)
+ return randomPermeability_[indexSet_.index(element.template subEntity<dim> (0))];
+ else
+ return permeabilityLens_;
+ }
else
return permeability_;
}
@@ -96,6 +125,36 @@ public:
const int scvIdx) const
{ return 0.4; }
+ /*!
+ * \brief This method allows the generation of a statistical field using gstat
+ *
+ * \param gridView The GridView used by the problem
+ */
+ void initRandomField(const GridView& gridView)
+ {
+ std::string gStatControlFile = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Gstat, ControlFile);
+ std::string gStatInputFile = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Gstat, InputFile);
+ std::string outputFilePrefix = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Gstat, OutputFilePrefix);
+
+ // create random permeability object
+ GstatRandomField<GridView, Scalar> randomPermeabilityField(gridView,
+ gStatControlFile,
+ gStatInputFile,
+ outputFilePrefix + ".dat",
+ true);
+ randomPermeabilityField.createPowField();
+ randomPermeability_.resize(gridView.size(dim), 0.0);
+
+ ElementIterator eItEnd = gridView.template end<0> ();
+ for (ElementIterator eIt = gridView.template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ randomPermeability_[indexSet_.index(eIt->template subEntity<dim> (0/*scvIdx*/))]
+ = randomPermeabilityField.data(*eIt);;
+ }
+
+ randomPermeabilityField.writeVtk(outputFilePrefix, "absolute permeability");
+ }
+
private:
bool isInLens_(const GlobalPosition &globalPos) const
{
@@ -106,10 +165,14 @@ private:
return true;
}
+ bool randomField_;
GlobalPosition lensLowerLeft_;
GlobalPosition lensUpperRight_;
Scalar permeability_, permeabilityLens_;
+ ScalarVector randomPermeability_;
+
+ const IndexSet& indexSet_;
};
} // end namespace
#endif
diff --git a/test/porousmediumflow/1p/implicit/CMakeLists.txt b/test/porousmediumflow/1p/implicit/CMakeLists.txt
index ae6d385d03d498361d0c2c9813a02ac9b5ad6ed1..600bf56e508066708776e20d2438254744218af0 100644
--- a/test/porousmediumflow/1p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/CMakeLists.txt
@@ -1,6 +1,7 @@
add_subdirectory(pointsources)
add_input_file_links()
+add_gstat_file_links()
# isothermal tests
add_dumux_test(test_box1p test_box1p test_box1p.cc
@@ -31,6 +32,9 @@ add_dumux_test(test_cc1pwithamg test_cc1pwithamg test_cc1pwithamg.cc
${CMAKE_CURRENT_BINARY_DIR}/1ptestccwithamg-00002.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/test_cc1pwithamg")
+add_dumux_test(test_cc1pwithgstat test_cc1pwithgstat test_cc1pwithgstat.cc
+ ./test_cc1pwithamg)
+
# non-isothermal tests
add_dumux_test(test_box1pniconduction test_box1pniconduction test_box1pniconduction.cc
python ${CMAKE_SOURCE_DIR}/bin/runtest.py
diff --git a/test/porousmediumflow/1p/implicit/control.gstat b/test/porousmediumflow/1p/implicit/control.gstat
new file mode 100644
index 0000000000000000000000000000000000000000..8b55ee804e64071b3c184769fb03cdd40efdf608
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/control.gstat
@@ -0,0 +1,24 @@
+# For more info see the gstat manual on www.gstat.org !
+#
+# unconditional Gaussian simulation (random field)
+# "perm" is the name of the data field.
+data(perm): dummy, sk_mean=-9.5, max=100; # dummy = no data for conditioning
+ # sk_mean = mean value
+ # higher max -> simulation more accurate and slower
+#
+#Sph:spherical, Gau:Gauss
+variogram(perm): 1.0000 Gau(0.05, 90, 1.0); # variogram definition: the first value is the variance
+ # Then specification of the model
+ # Brackets: 1: correlation length in main direction
+ # 2: zenith specifying the main direction
+ # 3: azimuth specifying the main direction
+ # 4: rotation of domain around main direction
+ # 5: the correlation length ratios in second direction
+ # 6: the correlation length ratios in second direction
+#
+set zero=1e-10; # higher values stabilize the simulation
+method: gs; # gaussian simulation
+data(): 'gstatInput.txt', x=1, y=2; # specification of simulation points: filename and order of data
+set output='permeability.dat'; # specification of output location
+set n_uk = 500; # higher n_uk -> simulation more accurate and slower.
+#
diff --git a/test/porousmediumflow/1p/implicit/test_cc1pwithgstat.cc b/test/porousmediumflow/1p/implicit/test_cc1pwithgstat.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e84848f6733d7a9d9d9da155439b60ee0bd7eb08
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/test_cc1pwithgstat.cc
@@ -0,0 +1,73 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief test for the one-phase CC model on geostatistical random field
+ * generated by gstat
+ */
+#include <config.h>
+#include "1ptestproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.LowerLeft Lower left corner coordinates\n"
+ "\t-Grid.UpperRight Upper right corner coordinates\n"
+ "\t-Grid.Cells Number of cells in respective coordinate directions\n"
+ "\t-SpatialParams.RandomField Set true to use a random field generated by gstat \n"
+ "\t-SpatialParams.LensLowerLeft Coordinates of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRight Coordinates of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-Gstat.ControlFile Control file for generating geostatistical field using gstat \n"
+ "\t-Gstat.InputFile Input file generated for gstat \n"
+ "\t-Gstat.OutputFilePrefix Prefix of the output file generated by gstat and for the vtu\n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+#if HAVE_GSTAT
+ typedef TTAG(OnePTestCCProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+#else
+#warning External geostatistical module gstat needed to run this example.
+ std::cerr << "Test skipped, it needs gstat!" << std::endl;
+ return 77;
+#endif
+}
diff --git a/test/porousmediumflow/1p/implicit/test_cc1pwithgstat.input b/test/porousmediumflow/1p/implicit/test_cc1pwithgstat.input
new file mode 100644
index 0000000000000000000000000000000000000000..f3709218917357f49f837497c806cd2c2cfe4945
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/test_cc1pwithgstat.input
@@ -0,0 +1,22 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 10 # [s]
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 1 1
+Cells = 10 10
+
+[Problem]
+Name = 1ptestccwithgstat # name passed to the output routines
+
+[SpatialParams]
+RandomField = true
+LensLowerLeft = 0.2 0.2
+LensUpperRight = 0.8 0.8
+Permeability = 1e-10 # [m^2]
+
+[Gstat]
+ControlFile = control.gstat
+InputFile = gstatInput.txt
+OutputFilePrefix = permeability