From 2885681ef1cd137934db1ec3213f1d949a5255fc Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 5 Dec 2012 12:49:09 +0000
Subject: [PATCH] implicit branch: add cell-centered treatment to the Richards
model, unify the test.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/branches/implicit@9782 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/richards/richardsmodel.hh | 85 ++++++++++++-------
.../richards/richardsnewtoncontroller.hh | 21 +++--
test/implicit/richards/Makefile.am | 5 +-
test/implicit/richards/richardslensproblem.hh | 11 ++-
test/implicit/richards/test_boxrichards.cc | 58 +++++++++++++
test/implicit/richards/test_boxrichards.input | 31 +++++++
.../{test_richards.cc => test_ccrichards.cc} | 2 +-
...t_richards.input => test_ccrichards.input} | 3 +
8 files changed, 173 insertions(+), 43 deletions(-)
create mode 100644 test/implicit/richards/test_boxrichards.cc
create mode 100644 test/implicit/richards/test_boxrichards.input
rename test/implicit/richards/{test_richards.cc => test_ccrichards.cc} (98%)
rename test/implicit/richards/{test_richards.input => test_ccrichards.input} (96%)
diff --git a/dumux/implicit/richards/richardsmodel.hh b/dumux/implicit/richards/richardsmodel.hh
index ff71be4f6f..8ed12767a9 100644
--- a/dumux/implicit/richards/richardsmodel.hh
+++ b/dumux/implicit/richards/richardsmodel.hh
@@ -119,6 +119,8 @@ class RichardsModel : public GET_PROP_TYPE(TypeTag, BaseModel)
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
enum { dim = GridView::dimension };
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief All relevant primary and secondary of a given
@@ -132,25 +134,27 @@ public:
{
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+ // get the number of degrees of freedom
+ unsigned numDofs = this->numDofs();
+
// create the required scalar fields
- unsigned numVertices = this->problem_().gridView().size(dim);
- ScalarField *pW = writer.allocateManagedBuffer(numVertices);
- ScalarField *pN = writer.allocateManagedBuffer(numVertices);
- ScalarField *pC = writer.allocateManagedBuffer(numVertices);
- ScalarField *Sw = writer.allocateManagedBuffer(numVertices);
- ScalarField *Sn = writer.allocateManagedBuffer(numVertices);
- ScalarField *rhoW = writer.allocateManagedBuffer(numVertices);
- ScalarField *rhoN = writer.allocateManagedBuffer(numVertices);
- ScalarField *mobW = writer.allocateManagedBuffer(numVertices);
- ScalarField *mobN = writer.allocateManagedBuffer(numVertices);
- ScalarField *poro = writer.allocateManagedBuffer(numVertices);
- ScalarField *Te = writer.allocateManagedBuffer(numVertices);
+ ScalarField *pW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *pN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *pC = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Sw = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Sn = writer.allocateManagedBuffer(numDofs);
+ ScalarField *rhoW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *rhoN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *mobW = writer.allocateManagedBuffer(numDofs);
+ ScalarField *mobN = writer.allocateManagedBuffer(numDofs);
+ ScalarField *poro = writer.allocateManagedBuffer(numDofs);
+ ScalarField *Te = writer.allocateManagedBuffer(numDofs);
unsigned numElements = this->gridView_().size(0);
ScalarField *rank =
writer.allocateManagedBuffer (numElements);
- FVElementGeometry fvElemGeom;
+ FVElementGeometry fvGeometry;
VolumeVariables volVars;
ElementIterator elemIt = this->gridView_().template begin<0>();
@@ -160,17 +164,21 @@ public:
int idx = this->problem_().model().elementMapper().map(*elemIt);
(*rank)[idx] = this->gridView_().comm().rank();
- fvElemGeom.update(this->gridView_(), *elemIt);
+ fvGeometry.update(this->gridView_(), *elemIt);
- int numVerts = elemIt->template count<dim> ();
- for (int i = 0; i < numVerts; ++i)
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
{
- int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ int globalIdx;
+ if (isBox)
+ globalIdx = this->vertexMapper().map(*elemIt, scvIdx, dim);
+ else
+ globalIdx = this->elementMapper().map(*elemIt);
+
volVars.update(sol[globalIdx],
this->problem_(),
*elemIt,
- fvElemGeom,
- i,
+ fvGeometry,
+ scvIdx,
false);
(*pW)[globalIdx] = volVars.pressure(wPhaseIdx);
@@ -187,17 +195,34 @@ public:
}
}
- writer.attachVertexData(*Sn, "Sn");
- writer.attachVertexData(*Sw, "Sw");
- writer.attachVertexData(*pN, "pn");
- writer.attachVertexData(*pW, "pw");
- writer.attachVertexData(*pC, "pc");
- writer.attachVertexData(*rhoW, "rhoW");
- writer.attachVertexData(*rhoN, "rhoN");
- writer.attachVertexData(*mobW, "mobW");
- writer.attachVertexData(*mobN, "mobN");
- writer.attachVertexData(*poro, "porosity");
- writer.attachVertexData(*Te, "temperature");
+ if (isBox) // vertex data
+ {
+ writer.attachVertexData(*Sn, "Sn");
+ writer.attachVertexData(*Sw, "Sw");
+ writer.attachVertexData(*pN, "pn");
+ writer.attachVertexData(*pW, "pw");
+ writer.attachVertexData(*pC, "pc");
+ writer.attachVertexData(*rhoW, "rhoW");
+ writer.attachVertexData(*rhoN, "rhoN");
+ writer.attachVertexData(*mobW, "mobW");
+ writer.attachVertexData(*mobN, "mobN");
+ writer.attachVertexData(*poro, "porosity");
+ writer.attachVertexData(*Te, "temperature");
+ }
+ else // cell data
+ {
+ writer.attachCellData(*Sn, "Sn");
+ writer.attachCellData(*Sw, "Sw");
+ writer.attachCellData(*pN, "pn");
+ writer.attachCellData(*pW, "pw");
+ writer.attachCellData(*pC, "pc");
+ writer.attachCellData(*rhoW, "rhoW");
+ writer.attachCellData(*rhoN, "rhoN");
+ writer.attachCellData(*mobW, "mobW");
+ writer.attachCellData(*mobN, "mobN");
+ writer.attachCellData(*poro, "porosity");
+ writer.attachCellData(*Te, "temperature");
+ }
writer.attachCellData(*rank, "process rank");
}
};
diff --git a/dumux/implicit/richards/richardsnewtoncontroller.hh b/dumux/implicit/richards/richardsnewtoncontroller.hh
index 635f6a6288..2f0b035d6d 100644
--- a/dumux/implicit/richards/richardsnewtoncontroller.hh
+++ b/dumux/implicit/richards/richardsnewtoncontroller.hh
@@ -54,6 +54,8 @@ class RichardsNewtonController : public NewtonController<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
enum { dim = GridView::dimension };
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
public:
/*!
* \brief Constructor
@@ -92,15 +94,20 @@ public:
const ElementIterator elemEndIt = gridView.template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
fvGeometry.update(gridView, *elemIt);
- for (int i = 0; i < fvGeometry.numSCV; ++i) {
- int globI = this->problem_().vertexMapper().map(*elemIt, i, dim);
+ for (int scvIdx = 0; scvIdx < fvGeometry.numSCV; ++scvIdx)
+ {
+ int globalIdx;
+ if (isBox)
+ globalIdx = this->problem_().vertexMapper().map(*elemIt, scvIdx, dim);
+ else
+ globalIdx = this->problem_().elementMapper().map(*elemIt);
// calculate the old wetting phase saturation
const SpatialParams &spatialParams = this->problem_().spatialParams();
- const MaterialLawParams &mp = spatialParams.materialLawParams(*elemIt, fvGeometry, i);
+ const MaterialLawParams &mp = spatialParams.materialLawParams(*elemIt, fvGeometry, scvIdx);
Scalar pcMin = MaterialLaw::pC(mp, 1.0);
- Scalar pW = uLastIter[globI][pwIdx];
- Scalar pN = std::max(this->problem_().referencePressure(*elemIt, fvGeometry, i),
+ Scalar pW = uLastIter[globalIdx][pwIdx];
+ Scalar pN = std::max(this->problem_().referencePressure(*elemIt, fvGeometry, scvIdx),
pW + pcMin);
Scalar pcOld = pN - pW;
Scalar SwOld = std::max<Scalar>(0.0, MaterialLaw::Sw(mp, pcOld));
@@ -111,9 +118,9 @@ public:
Scalar pwMax = pN - MaterialLaw::pC(mp, SwOld + 0.2);
// clamp the result
- pW = uCurrentIter[globI][pwIdx];
+ pW = uCurrentIter[globalIdx][pwIdx];
pW = std::max(pwMin, std::min(pW, pwMax));
- uCurrentIter[globI][pwIdx] = pW;
+ uCurrentIter[globalIdx][pwIdx] = pW;
}
}
diff --git a/test/implicit/richards/Makefile.am b/test/implicit/richards/Makefile.am
index 4f821cef5b..596324d827 100644
--- a/test/implicit/richards/Makefile.am
+++ b/test/implicit/richards/Makefile.am
@@ -1,8 +1,9 @@
-check_PROGRAMS = test_richards
+check_PROGRAMS = test_boxrichards test_ccrichards
noinst_HEADERS = *.hh
EXTRA_DIST=*reference*.vtu *.input grids/*.dgf CMakeLists.txt
-test_richards_SOURCES = test_richards.cc
+test_boxrichards_SOURCES = test_boxrichards.cc
+test_ccrichards_SOURCES = test_ccrichards.cc
include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/richards/richardslensproblem.hh b/test/implicit/richards/richardslensproblem.hh
index a33345db67..fb2f68d0da 100644
--- a/test/implicit/richards/richardslensproblem.hh
+++ b/test/implicit/richards/richardslensproblem.hh
@@ -45,7 +45,9 @@ class RichardsLensProblem;
//////////
namespace Properties
{
-NEW_TYPE_TAG(RichardsLensProblem, INHERITS_FROM(BoxRichards, RichardsLensSpatialParams));
+NEW_TYPE_TAG(RichardsLensProblem, INHERITS_FROM(Richards, RichardsLensSpatialParams));
+NEW_TYPE_TAG(RichardsLensBoxProblem, INHERITS_FROM(BoxModel, RichardsLensProblem));
+NEW_TYPE_TAG(RichardsLensCCProblem, INHERITS_FROM(CCModel, RichardsLensProblem));
// Use 2d YaspGrid
SET_TYPE_PROP(RichardsLensProblem, Grid, Dune::YaspGrid<2>);
@@ -156,6 +158,8 @@ public:
lensUpperRight_[1] = 3.0;
this->spatialParams().setLensCoords(lensLowerLeft_, lensUpperRight_);
+
+ name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
}
/*!
@@ -168,8 +172,8 @@ public:
*
* This is used as a prefix for files generated by the simulation.
*/
- const char *name() const
- { return "richardslens"; }
+ const std::string name() const
+ { return name_; }
/*!
* \brief Returns the temperature [K] within a finite volume
@@ -326,6 +330,7 @@ private:
GlobalPosition lensLowerLeft_;
GlobalPosition lensUpperRight_;
+ std::string name_;
};
} //end namespace
diff --git a/test/implicit/richards/test_boxrichards.cc b/test/implicit/richards/test_boxrichards.cc
new file mode 100644
index 0000000000..0e375eb705
--- /dev/null
+++ b/test/implicit/richards/test_boxrichards.cc
@@ -0,0 +1,58 @@
+// -*- 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 Richards box model.
+ */
+#include "config.h"
+#include "richardslensproblem.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 options for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(RichardsLensBoxProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/richards/test_boxrichards.input b/test/implicit/richards/test_boxrichards.input
new file mode 100644
index 0000000000..30bc118eb5
--- /dev/null
+++ b/test/implicit/richards/test_boxrichards.input
@@ -0,0 +1,31 @@
+###############################################################
+# Parameter file for test_richards.
+# Everything behind a '#' is a comment.
+# Type "./test_richards --help" for more information.
+###############################################################
+
+###############################################################
+# Mandatory arguments
+###############################################################
+
+[TimeManager]
+DtInitial = 100 # [s]
+TEnd = 3000 # [s]
+
+[Grid]
+File = ./grids/richardslens-24x16.dgf
+
+[Problem]
+Name = richardslensbox
+
+###############################################################
+# Simulation restart
+#
+# DuMux simulations can be restarted from *.drs files
+# Set Restart to the value of a specific file,
+# e.g.: 'Restart = 27184.1' for the restart file
+# name_time=27184.1_rank=0.drs
+# Please comment in the two lines below, if restart is desired.
+###############################################################
+# [TimeManager]
+# Restart = ...
diff --git a/test/implicit/richards/test_richards.cc b/test/implicit/richards/test_ccrichards.cc
similarity index 98%
rename from test/implicit/richards/test_richards.cc
rename to test/implicit/richards/test_ccrichards.cc
index b033c62fa7..349edf3981 100644
--- a/test/implicit/richards/test_richards.cc
+++ b/test/implicit/richards/test_ccrichards.cc
@@ -53,6 +53,6 @@ void usage(const char *progName, const std::string &errorMsg)
int main(int argc, char** argv)
{
- typedef TTAG(RichardsLensProblem) ProblemTypeTag;
+ typedef TTAG(RichardsLensCCProblem) ProblemTypeTag;
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
}
diff --git a/test/implicit/richards/test_richards.input b/test/implicit/richards/test_ccrichards.input
similarity index 96%
rename from test/implicit/richards/test_richards.input
rename to test/implicit/richards/test_ccrichards.input
index 8507fe30b1..348df3bbb2 100644
--- a/test/implicit/richards/test_richards.input
+++ b/test/implicit/richards/test_ccrichards.input
@@ -15,6 +15,9 @@ TEnd = 3000 # [s]
[Grid]
File = ./grids/richardslens-24x16.dgf
+[Problem]
+Name = richardslenscc
+
###############################################################
# Simulation restart
#
--
GitLab