Commit 2885681e authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

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
parent e0e4f687
......@@ -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");
}
};
......
......@@ -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;
}
}
......
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
......@@ -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
......
// -*- 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);
}
###############################################################
# 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 = ...
......@@ -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);
}
......@@ -15,6 +15,9 @@ TEnd = 3000 # [s]
[Grid]
File = ./grids/richardslens-24x16.dgf
[Problem]
Name = richardslenscc
###############################################################
# Simulation restart
#
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment