From b0db5d0e5d2686fea39012ecb761c07aeaa2d489 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Tue, 3 Aug 2010 17:11:29 +0000
Subject: [PATCH] big cleanup
timemanager API, restart, boxmodel are affected
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@3955 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
configure.ac | 3 +-
dumux/boxmodels/Makefile.am | 2 +-
dumux/boxmodels/boxscheme/Makefile.am | 4 -
dumux/boxmodels/boxscheme/boxjacobian.hh | 703 ------------------
dumux/boxmodels/boxscheme/boxproperties.hh | 310 --------
dumux/boxmodels/boxscheme/boxscheme.hh | 696 -----------------
.../boxmodels/boxscheme/fvelementgeometry.hh | 27 -
.../common/boxelementboundarytypes.hh | 110 +++
.../common/boxelementsecondaryvars.hh | 87 +++
.../boxfvelementgeometry.hh} | 77 +-
dumux/boxmodels/common/boxlocaljacobian.hh | 353 +++++++++
dumux/boxmodels/common/boxlocalresidual.hh | 456 ++++++++++++
dumux/boxmodels/common/boxmodel.hh | 587 +++++++++++++++
.../{boxscheme => common}/boxproblem.hh | 291 ++++----
dumux/boxmodels/common/boxproperties.hh | 361 +++++++++
.../pdelabboxassembler.hh} | 106 +--
.../pdelabboxistlvectorbackend.hh} | 4 +-
.../pdelabboxlocaloperator.hh} | 49 +-
dumux/boxmodels/pdelab/Makefile.am | 4 -
dumux/boxmodels/pdelab/boundarytypespdelab.hh | 49 --
.../pdelab/boxdirichletconstraints.hh | 235 ------
dumux/common/boundarytypes.hh | 92 +--
dumux/common/exceptions.hh | 2 +-
dumux/common/math.hh | 2 +-
.../pdelabpreconditioner.hh} | 124 +--
dumux/common/propertysystem.hh | 2 +-
dumux/common/spline.hh | 440 ++++++++++-
dumux/common/start.hh | 33 +-
dumux/common/timemanager.hh | 235 +++---
dumux/common/valgrind.hh | 2 +-
dumux/io/restart.hh | 76 +-
dumux/io/vtkmultiwriter.hh | 12 +-
.../binarycoefficients/fullermethod.hh | 2 +-
dumux/material/binarycoefficients/h2o_n2.hh | 2 +-
.../material/binarycoefficients/henryiapws.hh | 2 +-
dumux/material/components/component.hh | 8 +-
dumux/material/components/h2.hh | 2 +-
dumux/material/components/h2o.hh | 2 +-
dumux/material/components/iapws/common.hh | 2 +-
dumux/material/components/iapws/region1.hh | 2 +-
dumux/material/components/iapws/region2.hh | 2 +-
dumux/material/components/iapws/region4.hh | 2 +-
dumux/material/components/n2.hh | 2 +-
dumux/material/components/o2.hh | 2 +-
dumux/material/components/simplednapl.hh | 2 +-
dumux/material/components/simpleh2o.hh | 2 +-
.../material/components/tabulatedcomponent.hh | 2 +-
dumux/material/components/unit.hh | 2 +-
.../fluidmatrixinteractions/2p/brookscorey.hh | 2 +-
.../2p/brookscoreyparams.hh | 2 +-
.../fluidmatrixinteractions/2p/efftoabslaw.hh | 2 +-
.../2p/efftoabslawparams.hh | 2 +-
.../2p/linearmaterial.hh | 2 +-
.../2p/linearmaterialparams.hh | 2 +-
.../2p/regularizedbrookscorey.hh | 2 +-
.../2p/regularizedbrookscoreyparams.hh | 2 +-
.../2p/regularizedvangenuchten.hh | 2 +-
.../2p/regularizedvangenuchtenparams.hh | 2 +-
.../2p/vangenuchten.hh | 2 +-
.../2p/vangenuchtenparams.hh | 2 +-
dumux/material/fluidstate.hh | 6 +-
dumux/material/fluidsystems/2p_system.hh | 2 +-
dumux/material/fluidsystems/h2o_n2_system.hh | 242 ++----
dumux/material/idealgas.hh | 2 +-
dumux/material/settablephase.hh | 2 +-
.../spatialparameters/boxspatialparameters.hh | 48 +-
dumux/nonlinear/newtoncontroller.hh | 324 ++++----
dumux/nonlinear/newtonmethod.hh | 65 +-
68 files changed, 3271 insertions(+), 3016 deletions(-)
delete mode 100644 dumux/boxmodels/boxscheme/Makefile.am
delete mode 100644 dumux/boxmodels/boxscheme/boxjacobian.hh
delete mode 100644 dumux/boxmodels/boxscheme/boxproperties.hh
delete mode 100644 dumux/boxmodels/boxscheme/boxscheme.hh
delete mode 100644 dumux/boxmodels/boxscheme/fvelementgeometry.hh
create mode 100644 dumux/boxmodels/common/boxelementboundarytypes.hh
create mode 100644 dumux/boxmodels/common/boxelementsecondaryvars.hh
rename dumux/boxmodels/{boxscheme/fvelementgeometry-pdelab.hh => common/boxfvelementgeometry.hh} (93%)
create mode 100644 dumux/boxmodels/common/boxlocaljacobian.hh
create mode 100644 dumux/boxmodels/common/boxlocalresidual.hh
create mode 100644 dumux/boxmodels/common/boxmodel.hh
rename dumux/boxmodels/{boxscheme => common}/boxproblem.hh (60%)
create mode 100644 dumux/boxmodels/common/boxproperties.hh
rename dumux/boxmodels/{pdelab/assemblerpdelab.hh => common/pdelabboxassembler.hh} (76%)
rename dumux/boxmodels/{pdelab/istlvectorbackend.hh => common/pdelabboxistlvectorbackend.hh} (97%)
rename dumux/boxmodels/{pdelab/boxjacobianpdelab.hh => common/pdelabboxlocaloperator.hh} (68%)
delete mode 100644 dumux/boxmodels/pdelab/Makefile.am
delete mode 100644 dumux/boxmodels/pdelab/boundarytypespdelab.hh
delete mode 100644 dumux/boxmodels/pdelab/boxdirichletconstraints.hh
rename dumux/{boxmodels/pdelab/preconditionerpdelab.hh => common/pdelabpreconditioner.hh} (78%)
diff --git a/configure.ac b/configure.ac
index 02e9221cc7..55700aa1d4 100644
--- a/configure.ac
+++ b/configure.ac
@@ -95,8 +95,7 @@ AC_CONFIG_FILES([dumux.pc
dumux/boxmodels/2p2c/Makefile
dumux/boxmodels/2p2cni/Makefile
dumux/boxmodels/2pni/Makefile
- dumux/boxmodels/boxscheme/Makefile
- dumux/boxmodels/pdelab/Makefile
+ dumux/boxmodels/common/Makefile
dumux/boxmodels/richards/Makefile
dumux/common/Makefile
dumux/decoupled/Makefile
diff --git a/dumux/boxmodels/Makefile.am b/dumux/boxmodels/Makefile.am
index b00434af24..0fd56c7e68 100644
--- a/dumux/boxmodels/Makefile.am
+++ b/dumux/boxmodels/Makefile.am
@@ -1,4 +1,4 @@
-SUBDIRS = 1p 1p2c 2p 2p2c 2p2cni 2pni boxscheme pdelab richards
+SUBDIRS = common 1p 1p2c 2p 2p2c 2p2cni 2pni richards
boxmodelsdir = $(includedir)/dumux/boxmodels
diff --git a/dumux/boxmodels/boxscheme/Makefile.am b/dumux/boxmodels/boxscheme/Makefile.am
deleted file mode 100644
index 7030c59f01..0000000000
--- a/dumux/boxmodels/boxscheme/Makefile.am
+++ /dev/null
@@ -1,4 +0,0 @@
-boxschemedir = $(includedir)/dumux/boxmodels/boxscheme
-boxscheme_HEADERS = *.hh
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/boxmodels/boxscheme/boxjacobian.hh b/dumux/boxmodels/boxscheme/boxjacobian.hh
deleted file mode 100644
index 7823fe72b5..0000000000
--- a/dumux/boxmodels/boxscheme/boxjacobian.hh
+++ /dev/null
@@ -1,703 +0,0 @@
-// $Id: boxjacobian.hh 3784 2010-06-24 13:43:57Z bernd $
-/*****************************************************************************
- * Copyright (C) 2007 by Peter Bastian *
- * Institute of Parallel and Distributed System *
- * Department Simulation of Large Systems *
- * University of Stuttgart, Germany *
- * *
- * Copyright (C) 2008-2009 by Andreas Lauser *
- * Copyright (C) 2007-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
- * \brief Caculates the jacobian of models based on the box scheme element-wise.
- */
-#ifndef DUMUX_BOX_JACOBIAN_HH
-#define DUMUX_BOX_JACOBIAN_HH
-
-#include <dune/common/exceptions.hh>
-
-#include <dumux/common/valgrind.hh>
-#include <dumux/boxmodels/boxscheme/fvelementgeometry.hh>
-#include <dune/grid/common/genericreferenceelements.hh>
-
-#include <boost/format.hpp>
-
-#include <dune/common/fmatrix.hh>
-#include <dune/istl/matrix.hh>
-
-#include "boxproperties.hh"
-
-
-namespace Dumux
-{
-/*!
- * \ingroup BoxScheme
- * \brief Element-wise caculation of the jacobian matrix for models
- * based on the box scheme .
- *
- * \todo Please doc me more!
- */
-template<class TypeTag>
-class BoxJacobian
-{
-private:
- typedef BoxJacobian<TypeTag> ThisType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
-
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GridView::Grid::ctype CoordScalar;
-
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef typename Element::EntityPointer ElementPointer;
-
- typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
- typedef typename RefElemProp::Container ReferenceElements;
- typedef typename RefElemProp::ReferenceElement ReferenceElement;
-
- typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef typename Element::Geometry Geometry;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
-
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::DofEntityMapper DofEntityMapper;
- typedef typename SolutionTypes::VertexMapper VertexMapper;
- typedef typename SolutionTypes::ElementMapper ElementMapper;
- typedef typename SolutionTypes::SolutionVector SolutionVector;
- typedef typename SolutionTypes::SolutionOnElement SolutionOnElement;
- typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
- typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
- typedef typename SolutionTypes::JacobianAssembler JacobianAssembler;
-
-
- typedef typename GET_PROP(TypeTag, PTAG(VertexData))::type VertexData;
- typedef typename std::vector<VertexData> VertexDataArray;
-
- typedef std::vector<Dumux::BoundaryTypes<numEq> > BoundaryTypeArray;
- typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
- typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
-
-public:
- BoxJacobian(Problem &problem)
- : problem_(problem),
- gridView_(problem.gridView()),
- curElementPtr_(* ++gridView_.template begin<0>()),
- curElementGeom_(gridView_)
- {
- }
-
- /*!
- * \brief Compute the global residual right hand side
- * of an equation we would like to have zero.
- */
- void evalGlobalResidual(SolutionVector &residual)
- {
- residual = 0;
- SolutionOnElement tmpSol, tmpSolOld;
- SolutionOnElement localResid;
- localResid.resize(12);
- ElementIterator elemIt = gridView_.template begin<0>();
- const ElementIterator elemEndIt = gridView_.template end<0>();
- for (; elemIt != elemEndIt; ++elemIt) {
- this->setCurrentElement(*elemIt);
- this->restrictToElement(tmpSol, this->problem().model().curSol());
- this->restrictToElement(tmpSolOld, this->problem().model().prevSol());
- this->asImp_().setCurrentSolution(tmpSol);
- this->asImp_().setPreviousSolution(tmpSolOld);
-
- asImp_().evalLocalResidual(localResid);
-
- for (int i = 0; i < elemIt->template count<dim>(); ++i) {
- int globalI = this->problem().model().vertexMapper().map(*elemIt, i, dim);
- residual[globalI] += localResid[i];
- }
- };
- }
-
- /*!
- * \brief Assemble the linear system of equations for the
- * verts of a element, given a local solution 'localU'.
- */
- void assemble(const Element &element)
- {
- // set the current grid element
- asImp_().setCurrentElement(element);
-
- int numVertices = curElementGeom_.numVertices;
- SolutionOnElement localU(numVertices);
- restrictToElement(localU, problem_.model().curSol());
- asImp_().assemble_(element, localU);
- }
-
- /*!
- * \brief Compute the local residual, i.e. the right hand side
- * of an equation we would like to have zero.
- */
- void evalLocalResidual(SolutionOnElement &residual)
- {
- // reset residual
- for (int i = 0; i < curElementGeom_.numVertices; i++) {
- for (int j = 0; j < numEq; ++j)
- residual[i][j] = Scalar(0);
- }
-
- asImp_().evalFluxes_(residual);
- asImp_().evalVolumeTerms_(residual);
-
- // add the neumann fluxes
- asImp_().addNeumannFluxes_(residual);
- // set the defect of the equations used for dirichlet
- // conditions to 0
- for (int i = 0; i < curElementGeom_.numVertices; i++) {
- for (int j = 0; j < numEq; j++) {
- if (this->bctype[i].isDirichlet(j))
- residual[i][j] = 0;
- }
- }
-
-#if HAVE_VALGRIND
- for (int i=0; i < curElementGeom_.numVertices; i++)
- Valgrind::CheckDefined(residual[i]);
-#endif // HAVE_VALGRIND
- }
-
- /*!
- * \brief Restrict the global function 'globalFn' to the vertices
- * of the current element, save the result to 'dest'.
- */
- void restrictToElement(SolutionOnElement &dest,
- const SolutionVector &globalSol) const
- {
- const DofEntityMapper &dofMapper = this->problem_.model().dofEntityMapper();
- // we assert that the i-th shape function is
- // associated to the i-th vert of the element.
- int n = curElement_().template count<dim>();
- dest.resize(n);
- for (int i = 0; i < n; i++) {
- dest[i] = globalSol[dofMapper.map(curElement_(), i, dim)];
- }
- }
-
- /*!
- * \brief Set the current grid element.
- */
- void setCurrentElement(const Element &element)
- {
- if (curElementPtr_ != ElementPointer(element)) {
- // update the FV element geometry if the current element
- // is to be changed
- curElementPtr_ = element;
- curElementGeom_.update(curElement_());
-
- // resize the array for the local jacobian matrix to the
- // current number of degrees of freedom
- int n = curElementGeom_.numVertices;
- A.setSize(n, n);
- bctype.resize(n);
- updateBoundaryTypes_();
- }
- };
-
- /*!
- * \brief Initialize the static data of all elements. The current
- * solution is not yet valid when this method is
- * called. (updateStaticData() is called as soon the
- * current solution is valid.)
- *
- * This method should be overwritten by the child class if
- * necessary.
- */
- void initStaticData()
- { }
-
- /*!
- * \brief Update the static data of all elements with the current solution
- *
- * This method should be overwritten by the child class if
- * necessary.
- */
- void updateStaticData(SolutionVector &curSol, SolutionVector &oldSol)
- { };
-
- /*!
- * \brief Update the model specific vertex data of a whole
- * element.
- */
- void updateElementData_(VertexDataArray &dest, const SolutionOnElement &sol, bool isOldSol)
- {
- int n = sol.size();
- dest.resize(n);
-
-#ifdef ENABLE_VALGRIND
- for (int i = 0; i < dest.size(); ++i)
- Valgrind::SetUndefined(dest[i]);
-#endif // ENABLE_VALGRIND
-
- int numVertices = this->curElement_().template count<dim>();
- for (int vertIdx = 0; vertIdx < numVertices; vertIdx++) {
- dest[vertIdx].update(sol[vertIdx],
- curElement_(),
- curFvElementGeometry(),
- vertIdx,
- problem(),
- isOldSol);
- }
- }
-
-
- /*!
- * \brief This returns the finite volume geometric information of the current
- * element. <b>Use this method with great care!</b>
- */
- const FVElementGeometry &curFvElementGeometry() const
- { return curElementGeom_; }
-
- /*!
- * \brief Set current local solution
- */
- void setCurrentSolution(const SolutionOnElement &sol)
- {
- curElemDat_.resize(sol.size());
- asImp_().updateElementData_(curElemDat_, sol, false);
- }
-
- /*!
- * \brief Set local solution of the last time step
- */
- void setPreviousSolution(const SolutionOnElement &sol)
- {
- prevElemDat_.resize(sol.size());
- asImp_().updateElementData_(prevElemDat_, sol, true);
- }
-
- /*!
- * \brief Vary a single component of a single vert of the
- * local solution for the current element.
- *
- * This method is a optimization, since if varying a single
- * component at a degree of freedom not the whole element cache
- * needs to be recalculated. (Updating the element cache is very
- * expensive since material laws need to be evaluated.)
- */
- void deflectCurrentSolution(SolutionOnElement &curSol,
- int vertIdx,
- int eqIdx,
- Scalar value)
- {
- // stash away the orignial vertex data so that we do not need
- // to re calculate them if restoreCurSolution() is called.
- curVertexDataStash_ = curElemDat_[vertIdx];
-
- // recalculate the vertex data for the box which should be
- // changed
- curSol[vertIdx][eqIdx] = value;
-
- Valgrind::SetUndefined(curElemDat_[vertIdx]);
- curElemDat_[vertIdx].update(curSol[vertIdx],
- curElement_(),
- curFvElementGeometry(),
- vertIdx,
- problem(),
- false);
- Valgrind::CheckDefined(curElemDat_[vertIdx]);
- }
-
- /*!
- * \brief Restore the local jacobian to the state before
- * deflectCurSolution() was called.
- *
- * This only works if deflectSolution was only called with
- * (vert, component) as arguments.
- */
- void restoreCurrentSolution(SolutionOnElement &curSol,
- int vertIdx,
- int eqIdx,
- Scalar origValue)
- {
- curSol[vertIdx][eqIdx] = origValue;
- curElemDat_[vertIdx] = curVertexDataStash_;
- };
-
- /*!
- * \brief Returns a reference to the problem.
- */
- Problem &problem()
- { return problem_; };
-
- /*!
- * \brief Returns a reference to the problem.
- */
- const Problem &problem() const
- { return problem_; };
-
- /*!
- * \brief Returns a reference to the model.
- */
- Model &model()
- { return problem_.model(); };
-
- /*!
- * \brief Returns a reference to the model.
- */
- const Model &model() const
- { return problem_.model(); };
-
- /*!
- * \brief Returns a reference to the vertex mapper.
- */
- const VertexMapper &vertexMapper() const
- { return model().vertexMapper(); };
-
- /*!
- * \brief Returns a reference to the element mapper.
- */
- const ElementMapper &elementMapper() const
- { return model().elementMapper(); };
-
-
- const MatrixBlock &mat(int i, int j) const
- { return A[i][j]; }
-
- const VertexDataArray& currentElementData() const
- {
- return curElemDat_;
- }
-
-protected:
- const Element &curElement_() const
- { return *curElementPtr_; }
-
- // The problem we would like to solve
- Problem &problem_;
-
- // The grid view we are dealing with
- const GridView gridView_;
-
- ElementPointer curElementPtr_;
- FVElementGeometry curElementGeom_;
-
- // current and previous element data. (this is model specific.)
- VertexDataArray curElemDat_;
- VertexDataArray prevElemDat_;
-
- // temporary variable to store the variable vertex data
- VertexData curVertexDataStash_;
-
- void updateBoundaryTypes_()
- {
- Dune::GeometryType geoType = curElement_().geometry().type();
- const ReferenceElement &refElem = ReferenceElements::general(geoType);
-
- int numVerts = curElement_().template count<dim>();
- for (int i = 0; i < numVerts; ++i)
- this->bctype[i].reset();
-
- // evaluate boundary conditions
- IntersectionIterator isIt = gridView_.template ibegin(curElement_());
- const IntersectionIterator &endIt = gridView_.template iend(curElement_());
- for (; isIt != endIt; ++isIt)
- {
- // Ignore non- boundary faces.
- if (!isIt->boundary())
- continue;
-
- // Set the bctype for all vertices of the face
- int faceIdx = isIt->indexInInside();
- int numFaceVerts = refElem.size(faceIdx, 1, dim);
- for (int faceVertIdx = 0;
- faceVertIdx < numFaceVerts;
- faceVertIdx++)
- {
- int elemVertIdx = refElem.subEntity(faceIdx,
- 1,
- faceVertIdx,
- dim);
- int boundaryFaceIdx =
- curElementGeom_.boundaryFaceIndex(faceIdx,
- faceVertIdx);
- // set the boundary types
- problem_.boundaryTypes(this->bctype[elemVertIdx],
- curElement_(),
- curFvElementGeometry(),
- *isIt,
- elemVertIdx,
- boundaryFaceIdx);
- this->bctype[elemVertIdx].checkWellPosed();
- Valgrind::CheckDefined(this->bctype[elemVertIdx]);
- }
- }
- };
-
- void addNeumannFluxes_(SolutionOnElement &result)
- {
- Dune::GeometryType geoType = curElement_().geometry().type();
- const ReferenceElement &refElem = ReferenceElements::general(geoType);
-
- // evaluate boundary conditions for all intersections of
- // the current element
- IntersectionIterator isIt = gridView_.template ibegin(curElement_());
- const IntersectionIterator &endIt = gridView_.template iend(curElement_());
- for (; isIt != endIt; ++isIt)
- {
- // handle only faces on the boundary
- if (!isIt->boundary())
- continue;
-
- // Assemble the boundary for all vertices of the current
- // face
- int faceIdx = isIt->indexInInside();
- int numFaceVerts = refElem.size(faceIdx, 1, dim);
- for (int faceVertIdx = 0;
- faceVertIdx < numFaceVerts;
- ++faceVertIdx)
- {
- int elemVertIdx = refElem.subEntity(faceIdx,
- 1,
- faceVertIdx,
- dim);
-
- if (!this->bctype[elemVertIdx].hasNeumann())
- // the current boundary segment does not have any
- // equation where a neumann condition should be
- // applied
- continue;
-
- int boundaryFaceIdx =
- curElementGeom_.boundaryFaceIndex(faceIdx,
- faceVertIdx);
-
- // add the neuman fluxes of a single boundary segment
- addSingleNeumannSegment_(result[elemVertIdx],
- isIt,
- elemVertIdx,
- boundaryFaceIdx);
- }
- }
- }
-
- // handle boundary conditions for a single
- // sub-control volume face
- void addSingleNeumannSegment_(PrimaryVarVector &result,
- const IntersectionIterator &isIt,
- int scvIdx,
- int boundaryFaceIdx)
- {
- // temporary vector to store the neumann boundary fluxes
- PrimaryVarVector values(0.0);
-
- problem_.neumann(values,
- curElement_(),
- curElementGeom_,
- *isIt,
- scvIdx,
- boundaryFaceIdx);
- values *= curElementGeom_.boundaryFace[boundaryFaceIdx].area;
- Valgrind::CheckDefined(values);
-
- result += values;
- }
-
- void assemble_(const Element &element, SolutionOnElement& localU)
- {
- int numVertices = curElementGeom_.numVertices;
-
- // restrict the previous global solution to the current element
- SolutionOnElement localOldU(numVertices);
- restrictToElement(localOldU, problem_.model().prevSol());
-
- this->asImp_().setCurrentSolution(localU);
- this->asImp_().setPreviousSolution(localOldU);
-
- // approximate the local stiffness matrix numerically
- // TODO: do this analytically if possible
- SolutionOnElement partialStiffness(numVertices);
- for (int j = 0; j < numVertices; j++) {
- for (int pvIdx = 0; pvIdx < numEq; pvIdx++) {
- asImp_().assemblePartialStiffness_(partialStiffness,
- localU,
- j,
- pvIdx);
-
- // update the local stiffness matrix with the current partial
- // derivatives
- updateLocalStiffness_(j,
- pvIdx,
- partialStiffness);
- }
- }
- };
-
- void evalFluxes_(SolutionOnElement &residual)
- {
- // calculate the mass flux over the faces and subtract
- // it from the local rates
- for (int k = 0; k < curElementGeom_.numEdges; k++)
- {
- int i = curElementGeom_.subContVolFace[k].i;
- int j = curElementGeom_.subContVolFace[k].j;
-
- PrimaryVarVector flux;
- Valgrind::SetUndefined(flux);
- this->asImp_().computeFlux(flux, k);
- Valgrind::CheckDefined(flux);
-
- // subtract fluxes from the local mass rates of the
- // respective sub control volume adjacent to the face.
- for (int eq = 0; eq < numEq; ++ eq) {
- residual[i][eq] -= flux[eq];
- residual[j][eq] += flux[eq];
- }
- }
- }
-
- void evalVolumeTerms_(SolutionOnElement &residual)
- {
- // evaluate the volume terms (storage + source terms)
- for (int i=0; i < curElementGeom_.numVertices; i++)
- {
- PrimaryVarVector massContrib(0), tmp(0);
-
- // mass balance within the element. this is the
- // $\frac{m}{\partial t}$ term if using implicit
- // euler as time discretization.
- //
- // TODO (?): we might need a more explicit way for
- // doing the time discretization...
- this->asImp_().computeStorage(massContrib, i, false);
- this->asImp_().computeStorage(tmp, i, true);
-
- massContrib -= tmp;
- massContrib *=
- curElementGeom_.subContVol[i].volume
- /
- problem_.timeManager().timeStepSize();
-
- for (int j = 0; j < numEq; ++j)
- residual[i][j] += massContrib[j];
-
- // subtract the source term from the local rate
- PrimaryVarVector source;
- this->asImp_().computeSource(source, i);
- source *= curElementGeom_.subContVol[i].volume;
-
- for (int j = 0; j < numEq; ++j) {
- residual[i][j] -= source[j];
-
- // make sure that only defined quantities where used
- // to calculate the residual.
- Valgrind::CheckDefined(residual[i][j]);
- }
- }
- }
-
- /*!
- * \brief Update the stiffness matrix for all equations on all
- * vertices of the current element with the partial
- * derivative to the pvIdx's primary variable at the
- * vertexIdx-th vertex.
- *
- * This method can be overwritten by the implementation if a
- * better scheme than central differences ought to be used.
- */
- void assemblePartialStiffness_(SolutionOnElement &dest,
- SolutionOnElement &elemSol,
- int vertexIdx,
- int pvIdx)
- {
- Scalar eps = asImp_().numericEpsilon_(elemSol, vertexIdx, pvIdx);
- Scalar uJ = elemSol[vertexIdx][pvIdx];
-
- // vary the pvIdx-th primary variable at the element's j-th
- // vertex and calculate the residual, don't include the
- // boundary conditions
- asImp_().deflectCurrentSolution(elemSol, vertexIdx, pvIdx, uJ + eps);
- asImp_().evalLocalResidual(dest);
- asImp_().restoreCurrentSolution(elemSol, vertexIdx, pvIdx, uJ);
-
- asImp_().deflectCurrentSolution(elemSol, vertexIdx, pvIdx, uJ - eps);
- SolutionOnElement tmp(curElementGeom_.numVertices);
- asImp_().evalLocalResidual(tmp);
- asImp_().restoreCurrentSolution(elemSol, vertexIdx, pvIdx, uJ);
-
- // central differences
- dest -= tmp;
- dest /= 2*eps;
-
-#if HAVE_VALGRIND
- for (unsigned i = 0; i < dest.size(); ++i)
- Valgrind::CheckDefined(dest[i]);
-#endif
- }
-
- /*!
- * \brief Returns the epsilon value which is added and removed
- * from the current solution.
- *
- * \param elemSol The current solution on the element
- * \param vertexIdx The local index of the element's vertex for
- * which the local derivative ought to be calculated.
- * \param pvIdx The index of the primary variable which gets varied
- */
- Scalar numericEpsilon_(const SolutionOnElement &elemSol,
- int vertIdx,
- int pvIdx) const
- { return 1e-11*(std::abs(elemSol[vertIdx][pvIdx]) + 1); }
-
- /*!
- * \brief Updates the current local stiffness matrix with the
- * partial derivatives of all equations in regard to the
- * primary variable 'pvIdx' at vertex j .
- */
- void updateLocalStiffness_(int j,
- int pvIdx,
- const SolutionOnElement &stiffness)
- {
- for (int i = 0; i < curElementGeom_.numVertices; i++) {
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
- // A[i][j][eqIdx][pvIdx] is the approximate rate of
- // change of the residual of equation 'eqIdx' at
- // vertex 'i' depending on the primary variable
- // 'pvIdx' at vertex 'j'.
- this->A[i][j][eqIdx][pvIdx] = stiffness[i][eqIdx];
- Valgrind::CheckDefined(this->A[i][j][eqIdx][pvIdx]);
- }
- }
- }
-
- Implementation &asImp_()
- { return *static_cast<Implementation*>(this); }
-
- const Implementation &asImp_() const
- { return *static_cast<const Implementation*>(this); }
-
- LocalBlockMatrix A;
- BoundaryTypeArray bctype;
-};
-}
-
-#endif
diff --git a/dumux/boxmodels/boxscheme/boxproperties.hh b/dumux/boxmodels/boxscheme/boxproperties.hh
deleted file mode 100644
index 107a530981..0000000000
--- a/dumux/boxmodels/boxscheme/boxproperties.hh
+++ /dev/null
@@ -1,310 +0,0 @@
-// $Id: boxproperties.hh 3783 2010-06-24 11:33:53Z bernd $
-/*****************************************************************************
- * 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_BOX_PROPERTIES_HH
-#define DUMUX_BOX_PROPERTIES_HH
-
-#include <dumux/common/propertysystem.hh>
-
-#include <dumux/boxmodels/pdelab/istlvectorbackend.hh>
-
-/*!
- * \file
- * \brief Specify the shape functions, operator assemblers, etc
- * used for the BoxScheme.
- */
-namespace Dumux
-{
-namespace Properties
-{
-/*!
- * \addtogroup BoxScheme
- */
-// \{
-
-//////////////////////////////////////////////////////////////////
-// Type tags tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for models based on the box-scheme
-NEW_TYPE_TAG(BoxScheme);
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-//!< Property tag for scalar values
-NEW_PROP_TAG(Scalar);
-
-//! Property tag for types associated with the solution of the PDE.
-//! This means vectors of primary variables, solution functions on the
-//! grid, and elements, and shape functions.
-NEW_PROP_TAG(SolutionTypes);
-
-NEW_PROP_TAG(Grid); //!< The type of the DUNE grid
-NEW_PROP_TAG(GridView); //!< The type of the grid view
-
-NEW_PROP_TAG(ReferenceElements); //!< DUNE reference elements to be used
-NEW_PROP_TAG(FVElementGeometry); //! The type of the finite-volume geometry in the box scheme
-
-NEW_PROP_TAG(Problem); //!< The type of the problem
-NEW_PROP_TAG(Model); //!< The type of the discretization
-NEW_PROP_TAG(NumEq); //!< Number of equations in the system of PDEs
-NEW_PROP_TAG(LocalJacobian); //!< The type of the local jacobian operator
-
-#if HAVE_DUNE_PDELAB
-NEW_PROP_TAG(LocalFEMSpace); //!< The local finite element space used for the finite element interpolation
-NEW_PROP_TAG(PDELabTypes); //!< The types used from dune-pdelab
-#endif
-
-NEW_PROP_TAG(VertexData); //!< Data merging from constitutive relations defined on the vertices of the grid
-NEW_PROP_TAG(ElementData); //!< Data merging from constitutive relations defined on the elements of the grid
-NEW_PROP_TAG(FluxData); //!< Data required to calculate a flux over a face
-
-NEW_PROP_TAG(NewtonMethod); //!< The type of the newton method
-NEW_PROP_TAG(NewtonController); //!< The type of the newton controller
-}
-}
-
-#include <dumux/boxmodels/boxscheme/fvelementgeometry.hh>
-
-#include <dumux/nonlinear/newtonmethod.hh>
-#include <dumux/nonlinear/newtoncontroller.hh>
-
-#include <dumux/boxmodels/pdelab/assemblerpdelab.hh>
-#include <dumux/boxmodels/pdelab/boxdirichletconstraints.hh>
-
-#include <dumux/common/boundarytypes.hh>
-
-namespace Dumux {
-
-#if HAVE_DUNE_PDELAB
-template<typename TypeTag>
-class FVElementGeometry;
-#endif
-
-namespace Properties {
-//////////////////////////////////////////////////////////////////
-// Some defaults for very fundamental properties
-//////////////////////////////////////////////////////////////////
-
-//! Set the default type for scalar values to double
-SET_PROP_DEFAULT(Scalar)
-{ typedef double type; };
-
-//! Use the leaf grid view if not defined otherwise
-SET_PROP_DEFAULT(GridView)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
-
-public:
- typedef typename Grid::LeafGridView type;
-};
-
-/*!
- * \brief Specify the reference elements which we ought to use.
- *
- * We use Dune::ReferenceElements by default (-> old entity
- * numbering).
- *
- * TODO: Some specialization if the grid only supports one kind of
- * cells would be nice. this would be better fixed inside DUNE,
- * though. something like:
- * Dune::GenericReferenceElements<Dune::GeometryType<cube, dim> >
- */
-SET_PROP_DEFAULT(ReferenceElements)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
-
- typedef typename Grid::ctype CoordScalar;
- static const int dim = Grid::dimension;
-
-public:
-#if HAVE_DUNE_PDELAB
- typedef Dune::GenericReferenceElements<CoordScalar, dim> Container;
- typedef Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
-#else
- typedef Dune::ReferenceElements<CoordScalar, dim> Container;
- typedef Dune::ReferenceElement<CoordScalar, dim> ReferenceElement;
-#endif
-};
-
-//////////////////////////////////////////////////////////////////
-// Properties
-//////////////////////////////////////////////////////////////////
-
-//! Set the default for the FVElementGeometry
-SET_PROP(BoxScheme, FVElementGeometry)
-{
-#if HAVE_DUNE_PDELAB
- typedef Dumux::FVElementGeometry<TypeTag> type;
-
-#else
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
-public:
- typedef Dumux::FVElementGeometry<Grid> type;
-#endif
-};
-
-//! use the plain newton method for the box scheme by default
-SET_PROP(BoxScheme, NewtonMethod)
-{public:
- typedef Dumux::NewtonMethod<TypeTag> type;
-};
-
-//! use the plain newton controller for the box scheme by default
-SET_PROP(BoxScheme, NewtonController)
-{public:
- typedef Dumux::NewtonController<TypeTag> type;
-};
-
-/*!
- * \brief Specifies the types which are assoicated with a solution.
- *
- * This means shape functions, solution vectors, etc.
- */
-SET_PROP(BoxScheme, SolutionTypes)
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GridView::Grid Grid;
- typedef typename Grid::ctype CoordScalar;
-
- enum {
- dim = GridView::dimension,
- numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))
- };
-
- template<int dim>
- struct VertexLayout {
- bool contains (Dune::GeometryType gt) const
- { return gt.dim() == 0; }
- };
-
- template<int dim>
- struct ElementLayout {
- bool contains (Dune::GeometryType gt) const
- { return gt.dim() == dim; }
- };
-
- typedef typename GET_PROP(TypeTag, PTAG(PDELabTypes)) PDELabTypes;
- typedef typename PDELabTypes::GridFunctionSpace GridFunctionSpace;
- typedef typename PDELabTypes::GridOperatorSpace GridOperatorSpace;
-
-public:
- /*!
- * \brief Mapper for the grid view's vertices.
- */
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, VertexLayout> VertexMapper;
-
- /*!
- * \brief Mapper for the grid view's elements.
- */
- typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, ElementLayout> ElementMapper;
-
- /*!
- * \brief The type which maps an entity with an attached degree of
- * freedom to an index in the solution.
- */
- typedef VertexMapper DofEntityMapper;
-
- /*!
- * \brief The type of a solution at a fixed time.
- *
- * This is the representation of the solution function and defines
- * a primary variable vector at each degree of freedom.
- */
- typedef typename GridFunctionSpace::template VectorContainer<Scalar>::Type SolutionVector;
- typedef typename GridOperatorSpace::template MatrixContainer<Scalar>::Type Matrix;
-
- /*!
- * \brief A vector of primary variables.
- */
- typedef typename SolutionVector::block_type PrimaryVarVector;
-
- /*!
- * \brief The solution for a single finite element.
- */
- typedef Dune::BlockVector<PrimaryVarVector> SolutionOnElement;
-
- /*!
- * \brief Vector of boundary types at a degree of freedom.
- */
- typedef Dumux::BoundaryTypes<numEq> BoundaryTypeVector;
-
- /*!
- * \brief Assembler for the global jacobian matrix.
- */
-#ifdef HAVE_DUNE_PDELAB
- typedef AssemblerPDELab<TypeTag> JacobianAssembler;
-#else
- typedef Dune::P1OperatorAssembler<Grid, Scalar, GridView, GridView, numEq> JacobianAssembler;
-#endif
-};
-
-
-
-#ifdef HAVE_DUNE_PDELAB
-//! use the local FEM space associated with cubes by default
-SET_PROP(BoxScheme, 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;
-};
-
-SET_PROP(BoxScheme, PDELabTypes)
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
-public:
- //typedef typename Dune::PDELab::NonoverlappingConformingDirichletConstraints Constraints;
- typedef typename Dumux::NonoverlappingBoxDirichletConstraints<TypeTag> Constraints;
- typedef Dune::PDELab::GridFunctionSpace<GridView,
- FEM,
- Constraints,
- Dumux::PDELab::ISTLVectorBackend<TypeTag>
- > ScalarGridFunctionSpace;
- typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, numEq,
- Dune::PDELab::GridFunctionSpaceBlockwiseMapper> GridFunctionSpace;
- typedef typename GridFunctionSpace::template ConstraintsContainer<Scalar>::Type ConstraintsTrafo;
- typedef BoxJacobianPDELab<TypeTag> LocalOperator;
- typedef Dune::PDELab::GridOperatorSpace<GridFunctionSpace,
- GridFunctionSpace,
- LocalOperator,
- ConstraintsTrafo,
- ConstraintsTrafo,
- Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
- true
- > GridOperatorSpace;
-};
-#endif
-
-
-// \}
-
-}
-}
-
-#endif
diff --git a/dumux/boxmodels/boxscheme/boxscheme.hh b/dumux/boxmodels/boxscheme/boxscheme.hh
deleted file mode 100644
index d6f7234930..0000000000
--- a/dumux/boxmodels/boxscheme/boxscheme.hh
+++ /dev/null
@@ -1,696 +0,0 @@
-// $Id: boxscheme.hh 3784 2010-06-24 13:43:57Z bernd $
-/*****************************************************************************
- * Copyright (C) 2007 by Peter Bastian *
- * Institute of Parallel and Distributed System *
- * Department Simulation of Large Systems *
- * University of Stuttgart, Germany *
- * *
- * Copyright (C) 2008 by Andreas Lauser, 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_BOX_SCHEME_HH
-#define DUMUX_BOX_SCHEME_HH
-
-#include <dumux/boxmodels/boxscheme/boxjacobian.hh>
-#include <dumux/common/valgrind.hh>
-
-#include <dune/istl/operators.hh>
-#include <dune/grid/common/genericreferenceelements.hh>
-
-#include <boost/format.hpp>
-
-#include "boxproperties.hh"
-
-namespace Dumux
-{
-
-/*!
- * \defgroup BoxScheme Box-Scheme
- */
-/*!
- * \ingroup BoxScheme
- * \defgroup BoxProblems Box-Problems
- */
-/*!
- * \ingroup BoxScheme
- * \defgroup BoxModels Box-Models
- */
-
-
-/*!
- * \ingroup BoxScheme
- *
- * \brief The base class for the vertex centered finite volume
- * discretization scheme.
- */
-template<class TypeTag, class Implementation>
-class BoxScheme
-{
- typedef BoxScheme<TypeTag, Implementation> ThisType;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GridView::Grid::ctype CoordScalar;
-
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::DofEntityMapper DofEntityMapper;
- typedef typename SolutionTypes::VertexMapper VertexMapper;
- typedef typename SolutionTypes::ElementMapper ElementMapper;
- typedef typename SolutionTypes::SolutionVector SolutionVector;
- typedef typename SolutionTypes::SolutionOnElement SolutionOnElement;
- typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
- typedef typename SolutionTypes::BoundaryTypeVector BoundaryTypeVector;
- typedef typename SolutionTypes::JacobianAssembler JacobianAssembler;
-
- enum {
- numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
- dim = GridView::dimension
- };
-
- typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
- typedef typename RefElemProp::Container ReferenceElements;
- typedef typename RefElemProp::ReferenceElement ReferenceElement;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexData)) VertexData;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef typename GridView::template Codim<dim>::Entity Vertex;
- typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
-
-public:
- /*!
- * \brief This structure is Required to use models based on the BOX
- * scheme in conjunction with the newton Method.
- *
- * \todo change the newton method to the property system!
- */
- struct NewtonTraits {
- typedef typename ThisType::LocalJacobian LocalJacobian;
- typedef typename ThisType::SolutionVector Vector;
- typedef typename ThisType::JacobianAssembler JacobianAssembler;
- typedef typename ThisType::Scalar Scalar;
- typedef typename ThisType::GridView::Grid Grid;
- };
-
- /*!
- * \brief The constructor.
- */
- BoxScheme(Problem &prob)
- : problem_(prob),
- gridView_(prob.gridView()),
- dofEntityMapper_(gridView_),
- vertexMapper_(gridView_),
- elementMapper_(gridView_),
- localJacobian_(prob)
- {
- jacAsm_ = new JacobianAssembler(asImp_(), problem_);
-
- wasRestarted_ = false;
-
- uCur_.resize(asImp_().numDofs());
- uPrev_.resize(asImp_().numDofs());
- f_.resize(asImp_().numDofs());
-
- // check grid partitioning if we are parallel
-// assert((prob.gridView().comm().size() == 1) ||
-// (prob.gridView().overlapSize(0) > 0) ||
-// (prob.gridView().ghostSize(0) > 0));
- }
-
- /*!
- * \brief Apply the initial conditions to the model.
- */
- void initial()
- {
- if (!wasRestarted_) {
- this->localJacobian().initStaticData();
- applyInitialSolution_(uCur_);
- }
-
- applyDirichletBoundaries_(uCur_);
-
- // also set the solution of the "previous" time step to the
- // initial solution.
- uPrev_ = uCur_;
- }
-
- Scalar globalResidual(const SolutionVector &u, SolutionVector &tmp)
- {
- SolutionVector tmpU(asImp_(), 0.0);
- tmpU = uCur_;
- uCur_ = u;
- localJacobian_.evalGlobalResidual(tmp);
-
- Scalar result = tmp.two_norm();
- /*
- Scalar result = 0;
- for (int i = 0; i < (*tmp).size(); ++i) {
- for (int j = 0; j < numEq; ++j)
- result += std::abs((*tmp)[i][j]);
- }
- */
- uCur_ = tmpU;
- return result;
- };
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- const SolutionVector &curSol() const
- { return uCur_; }
-
- /*!
- * \brief Reference to the current solution as a block vector.
- */
- SolutionVector &curSol()
- { return uCur_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- const SolutionVector &prevSol() const
- { return uPrev_; }
-
- /*!
- * \brief Reference to the previous solution as a block vector.
- */
- SolutionVector &prevSol()
- { return uPrev_; }
-
- /*!
- * \brief Returns the operator assembler for the global jacobian of
- * the problem.
- */
- JacobianAssembler &jacobianAssembler()
- { return *jacAsm_; }
-
- /*!
- * \brief Returns the local jacobian which calculates the local
- * stiffness matrix for an arbitrary element.
- *
- * The local stiffness matrices of the element are used by
- * the jacobian assembler to produce a global linerization of the
- * problem.
- */
- LocalJacobian &localJacobian()
- { return localJacobian_; }
- /*!
- * \copydoc localJacobian()
- */
- const LocalJacobian &localJacobian() const
- { return localJacobian_; }
-
-
- /*!
- * \brief A reference to the problem on which the model is applied.
- */
- Problem &problem()
- { return problem_; }
- /*!
- * \copydoc problem()
- */
- const Problem &problem() const
- { return problem_; }
-
- /*!
- * \brief Reference to the grid view of the spatial domain.
- */
- const GridView &gridView() const
- { return gridView_; }
-
- /*!
- * \brief Try to progress the model to the next timestep.
- */
- void update(NewtonMethod &solver,
- NewtonController &controller)
- {
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i)
- Valgrind::CheckDefined(curSol()[i]);
-#endif // HAVE_VALGRIND
-
- asImp_().updateBegin();
-
- int numRetries = 0;
- while (true)
- {
- bool converged = solver.execute(controller);
- if (converged)
- break;
-
- ++numRetries;
- if (numRetries > 10) {
- problem_.updateFailed();
- asImp_().updateFailed();
- DUNE_THROW(Dune::MathError,
- "Newton solver didn't converge after 10 timestep divisions. dt="
- << problem_.timeManager().timeStepSize());
- }
-
- problem_.updateFailedTry();
- asImp_().updateFailedTry();
-
- if (verbose_())
- std::cout << boost::format("Newton didn't converge for rank %d. Retrying with timestep of %f\n")
- %gridView_.comm().rank()
- %problem_.timeManager().timeStepSize();
- }
-
- problem_.updateSuccessful();
- asImp_().updateSuccessful();
-
-#if HAVE_VALGRIND
- for (size_t i = 0; i < curSol().size(); ++i) {
- Valgrind::CheckDefined(curSol()[i]);
- }
-#endif // HAVE_VALGRIND
- }
-
-
- /*!
- * \brief Called by the update() method before it tries to
- * apply the newton method. This is primary a hook
- * which the actual model can overload.
- */
- void updateBegin()
- {
- applyDirichletBoundaries_(curSol());
- }
-
-
- /*!
- * \brief Called by the update() method if it was
- * successful. This is primary a hook which the actual
- * model can overload.
- */
- void updateSuccessful()
- {
- // make the current solution the previous one.
- uPrev_ = uCur_;
- };
-
- /*!
- * \brief Called by the update() method if a try was ultimately
- * unsuccessful. This is primary a hook which the
- * actual model can overload.
- */
- void updateFailed()
- {
- };
-
- /*!
- * \brief Called by the update() method if a try was
- * unsuccessful. This is primary a hook which the
- * actual model can overload.
- */
- void updateFailedTry()
- {
- // Reset the current solution to the one of the
- // previous time step so that we can start the next
- // update at a physically meaningful solution.
- uCur_ = uPrev_;
- applyDirichletBoundaries_(curSol());
- };
-
- /*!
- * \brief Calculate the global residual.
- *
- * The global deflection of the mass balance from zero.
- */
- void evalGlobalResidual(SolutionVector &globResidual)
- {
- globResidual = Scalar(0.0);
-
- // iterate through leaf grid
- ElementIterator it = gridView_.template begin<0>();
- const ElementIterator &eendit = gridView_.template end<0>();
- for (; it != eendit; ++it)
- {
- // tell the local jacobian which element it should
- // consider and evaluate the local residual for the
- // element. in order to do this we first have to
- // evaluate the element's local solutions for the
- // current and the last timestep.
- const Element& element = *it;
-
- const int numDofs = element.template count<dim>();
- SolutionOnElement localResidual(numDofs);
-
- SolutionOnElement localU(numDofs);
- SolutionOnElement localOldU(numDofs);
-
- localJacobian_.setCurrentElement(element);
- localJacobian_.restrictToElement(localU, curSol());
- localJacobian_.restrictToElement(localOldU, prevSol());
-
- localJacobian_.setCurrentSolution(localU);
- localJacobian_.setPreviousSolution(localOldU);
-
- localJacobian_.evalLocalResidual(localResidual);
-
- // loop over the element's shape functions, map their
- // associated degree of freedom to the corresponding
- // indices in the solution vector and add the element's
- // local residual at the index to the global residual at
- // this index.
- for(int dofIdx=0; dofIdx < numDofs; dofIdx++)
- {
- int globalIdx = dofEntityMapper().map(element, dofIdx, dim);
- (*globResidual)[globalIdx] += localResidual[dofIdx];
- }
- }
- }
-
- /*!
- * \brief Serializes the current state of the model.
- */
- template <class Restarter>
- void serialize(Restarter &res)
- { res.template serializeEntities<dim>(asImp_(), this->gridView_); }
-
- /*!
- * \brief Deserializes the state of the model.
- */
- template <class Restarter>
- void deserialize(Restarter &res)
- {
- res.template deserializeEntities<dim>(asImp_(), this->gridView_);
- wasRestarted_ = true;
- }
-
- /*!
- * \brief Write the current solution for a vertex to a restart
- * file.
- */
- void serializeEntity(std::ostream &outstream,
- const Vertex &vert)
- {
- int vertIdx = dofEntityMapper().map(vert);
-
- // write phase state
- if (!outstream.good()) {
- DUNE_THROW(Dune::IOError,
- "Could not serialize vertex "
- << vertIdx);
- }
-
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- outstream << curSol()[vertIdx][eqIdx] << " ";
- }
- };
-
- /*!
- * \brief Reads the current solution variables for a vertex from a
- * restart file.
- */
- void deserializeEntity(std::istream &instream,
- const Vertex &vert)
- {
- int vertIdx = dofEntityMapper().map(vert);
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
- if (!instream.good())
- DUNE_THROW(Dune::IOError,
- "Could not deserialize vertex "
- << vertIdx);
- instream >> curSol()[vertIdx][eqIdx];
- }
- };
-
- /*!
- * \brief Mapper for the entities where degrees of freedoms are
- * defined to indices.
- *
- * This usually means a mapper for vertices.
- */
- const DofEntityMapper &dofEntityMapper() const
- { return dofEntityMapper_; };
-
- /*!
- * \brief Returns the number of global degrees of freedoms (DOFs)
- */
- size_t numDofs() const
- { return gridView_.size(dim); }
-
- /*!
- * \brief Mapper for vertices to indices.
- */
- const VertexMapper &vertexMapper() const
- { return vertexMapper_; };
-
- /*!
- * \brief Mapper for elements to indices.
- */
- const ElementMapper &elementMapper() const
- { return elementMapper_; };
-
- void resetJacobianAssembler ()
- {
- delete jacAsm_;
-
- jacAsm_ = new JacobianAssembler(asImp_(), problem_);
- }
-
- ~BoxScheme()
- {
- delete jacAsm_;
- }
-
-protected:
- //! returns true iff the grid has an overlap
- bool hasOverlap_()
- { return gridView_.overlapSize(0) > 0; };
-
- void applyInitialSolution_(SolutionVector &u)
- {
- // first set the whole domain to zero. This is
- // necessary in order to also get a meaningful value
- // for ghost nodes (if we are running in parallel)
-// if (gridView_.comm().size() > 1) {
-// u = Scalar(0.0);
-// }
-
- // iterate through leaf grid and evaluate initial
- // condition at the center of each sub control volume
- //
- // TODO: the initial condition needs to be unique for
- // each vertex. we should think about the API...
- ElementIterator it = gridView_.template begin<0>();
- const ElementIterator &eendit = gridView_.template end<0>();
- for (; it != eendit; ++it)
- {
- // deal with the current element
- applyInitialSolutionElement_(u, *it);
- }
- };
-
- // apply the initial solition for a single element
- void applyInitialSolutionElement_(SolutionVector &u,
- const Element &element)
- {
- // HACK: set the current element for the local
- // solution in order to get an updated FVElementGeometry
- localJacobian_.setCurrentElement(element);
-
- // loop over all element vertices, i.e. sub control volumes
- int numScv = element.template count<dim>();
- for (int scvIdx = 0;
- scvIdx < numScv;
- scvIdx++)
- {
- // map the local vertex index to the global one
- int globalIdx = dofEntityMapper().map(element,
- scvIdx,
- dim);
-
- const FVElementGeometry &fvElemGeom
- = localJacobian_.curFvElementGeometry();
- // use the problem for actually doing the
- // dirty work of nailing down the initial
- // solution.
- this->problem_.initial(u[globalIdx],
- element,
- fvElemGeom,
- scvIdx);
- Valgrind::CheckDefined(u[globalIdx]);
- }
- }
-
- // apply dirichlet boundaries for the whole grid
- void applyDirichletBoundaries_(SolutionVector &u)
- {
- // set Dirichlet boundary conditions of the grid's
- // outer boundaries
- ElementIterator elementIt = gridView_.template begin<0>();
- const ElementIterator &elementEndIt = gridView_.template end<0>();
- for (; elementIt != elementEndIt; ++elementIt)
- {
- // ignore elements which are not on the boundary of
- // the domain
- if (!elementIt->hasBoundaryIntersections())
- continue;
-
- // set the current element of the local jacobian
- localJacobian_.setCurrentElement(*elementIt);
-
- // apply dirichlet boundary for the current element
- applyDirichletElement_(u, *elementIt);
- }
- };
-
- // apply dirichlet boundaries for a single element
- void applyDirichletElement_(SolutionVector &u,
- const Element &element)
- {
- Dune::GeometryType geoType = element.geometry().type();
- const ReferenceElement &refElem = ReferenceElements::general(geoType);
-
- // loop over all the element's surface patches
- IntersectionIterator isIt = gridView_.ibegin(element);
- const IntersectionIterator &isEndIt = gridView_.iend(element);
- for (; isIt != isEndIt; ++isIt) {
- // if the current intersection is not on the boundary,
- // we ignore it
- if (!isIt->boundary())
- continue;
-
- // Assemble the boundary for all vertices of the
- // current face
- int faceIdx = isIt->indexInInside();
- int numVerticesOfFace = refElem.size(faceIdx, 1, dim);
- for (int vertInFace = 0;
- vertInFace < numVerticesOfFace;
- vertInFace++)
- {
- // apply dirichlet boundaries for the current
- // sub-control volume face
- applyDirichletSCVF_(u, element, refElem, isIt, vertInFace);
- }
- }
- }
-
- // apply dirichlet boundaries for a single boundary
- // sub-control volume face of a finite volume cell.
- void applyDirichletSCVF_(SolutionVector &u,
- const Element &element,
- const ReferenceElement &refElem,
- const IntersectionIterator &isIt,
- int faceVertIdx)
- {
- // apply dirichlet boundaries but make sure
- // not to interfere with non-dirichlet
- // boundaries...
- const FVElementGeometry &fvElemGeom =
- localJacobian_.curFvElementGeometry();
-
- int elemVertIdx = refElem.subEntity(isIt->indexInInside(), 1,
- faceVertIdx, dim);
- int boundaryFaceIdx = fvElemGeom.boundaryFaceIndex(isIt->indexInInside(),
- faceVertIdx);
- int globalVertexIdx = dofEntityMapper().map(element,
- elemVertIdx, dim);
-
-
- PrimaryVarVector dirichletVal;
- BoundaryTypeVector boundaryTypes;
- problem_.boundaryTypes(boundaryTypes,
- element,
- fvElemGeom,
- *isIt,
- elemVertIdx,
- boundaryFaceIdx);
-
- bool dirichletEvaluated = false;
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- // ignore non-dirichlet boundary conditions
- if (!boundaryTypes.isDirichlet(eqIdx))
- { continue; }
-
- // make sure to evaluate the dirichlet boundary
- // conditions exactly once (and only if the boundary
- // type is actually dirichlet).
- if (!dirichletEvaluated)
- {
- dirichletEvaluated = true;
- problem_.dirichlet(dirichletVal,
- element,
- fvElemGeom,
- *isIt,
- elemVertIdx,
- boundaryFaceIdx);
- Valgrind::CheckDefined(dirichletVal);
- }
-
- // copy the dirichlet value for the current equation
- // to the global solution.
- //
- // TODO: we should propably use the sum weighted by the
- // sub-control volume instead of just overwriting
- // the previous values...
- u[globalVertexIdx][eqIdx] = dirichletVal[boundaryTypes.eqToDirichletIndex(eqIdx)];
- }
- }
-
- bool verbose_() const
- { return gridView_.comm().rank() == 0; };
-
- Implementation &asImp_()
- { return *static_cast<Implementation*>(this); }
- const Implementation &asImp_() const
- { return *static_cast<const Implementation*>(this); }
-
- // the problem we want to solve. defines the constitutive
- // relations, material laws, etc.
- Problem &problem_;
-
- // the grid view for which we need a solution
- const GridView gridView_;
-
- // mapper for the entities of a solution to their indices
- const DofEntityMapper dofEntityMapper_;
-
- // mapper for the vertices to indices
- const VertexMapper vertexMapper_;
-
- // mapper for the elements to indices
- const ElementMapper elementMapper_;
-
- // the solution we are looking for
-
- // calculates the local jacobian matrix for a given element
- LocalJacobian localJacobian_;
- // Linearizes the problem at the current time step using the
- // local jacobian
- JacobianAssembler *jacAsm_;
-
- // cur is the current solution, prev the solution of the previous
- // time step
- SolutionVector uCur_;
- SolutionVector uPrev_;
- // the right hand side
- SolutionVector f_;
-
- bool wasRestarted_;
-};
-}
-
-#endif
diff --git a/dumux/boxmodels/boxscheme/fvelementgeometry.hh b/dumux/boxmodels/boxscheme/fvelementgeometry.hh
deleted file mode 100644
index f96c0b51bb..0000000000
--- a/dumux/boxmodels/boxscheme/fvelementgeometry.hh
+++ /dev/null
@@ -1,27 +0,0 @@
-// $Id: fvelementgeometry.hh 3774 2010-06-24 06:22:44Z bernd $
-/*****************************************************************************
- * Copyright (C) 2008-2009 by Bernd Flemisch, 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_FVELEMENTGEOMETRY_HH
-#define DUMUX_FVELEMENTGEOMETRY_HH
-
-#if HAVE_DUNE_PDELAB
-#include "fvelementgeometry-pdelab.hh"
-#else
-#include "fvelementgeometry-disc.hh"
-#endif
-
-#endif
-
diff --git a/dumux/boxmodels/common/boxelementboundarytypes.hh b/dumux/boxmodels/common/boxelementboundarytypes.hh
new file mode 100644
index 0000000000..6c206e7342
--- /dev/null
+++ b/dumux/boxmodels/common/boxelementboundarytypes.hh
@@ -0,0 +1,110 @@
+// $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. *
+ *****************************************************************************/
+#ifndef DUMUX_BOX_ELEMENT_BOUNDARY_TYPES_HH
+#define DUMUX_BOX_ELEMENT_BOUNDARY_TYPES_HH
+
+#include "boxproperties.hh"
+
+#include <vector>
+#include <dumux/common/boundarytypes.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup BoxModel
+ *
+ * \brief This class stores an array of BoundaryTypes objects
+ */
+template<class TypeTag>
+class BoxElementBoundaryTypes : public std::vector<typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) >
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef std::vector<BoundaryTypes> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
+ typedef typename RefElemProp::Container ReferenceElements;
+ typedef typename RefElemProp::ReferenceElement ReferenceElement;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ enum { dim = GridView::dimension };
+
+public:
+ /*!
+ * \brief The constructor.
+ */
+ BoxElementBoundaryTypes()
+ { }
+
+ void update(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom)
+ {
+ Dune::GeometryType geoType = element.geometry().type();
+ const ReferenceElement &refElem = ReferenceElements::general(geoType);
+
+ int numVerts = element.template count<dim>();
+ this->resize(numVerts);
+ for (int i = 0; i < numVerts; ++i)
+ (*this)[i].reset();
+
+ // evaluate boundary conditions
+ IntersectionIterator isIt = problem.gridView().template ibegin(element);
+ const IntersectionIterator &endIt = problem.gridView().template iend(element);
+ for (; isIt != endIt; ++isIt) {
+ // Ignore non- boundary faces.
+ if (!isIt->boundary())
+ continue;
+
+ // Set the boundary type for all vertices of the face
+ int faceIdx = isIt->indexInInside();
+ int numFaceVerts = refElem.size(faceIdx, 1, dim);
+ for (int faceVertIdx = 0;
+ faceVertIdx < numFaceVerts;
+ faceVertIdx++)
+ {
+ int elemVertIdx = refElem.subEntity(faceIdx,
+ 1,
+ faceVertIdx,
+ dim);
+ int boundaryFaceIdx =
+ fvElemGeom.boundaryFaceIndex(faceIdx,
+ faceVertIdx);
+ // set the boundary types
+ problem.boundaryTypes((*this)[elemVertIdx],
+ element,
+ fvElemGeom,
+ *isIt,
+ elemVertIdx,
+ boundaryFaceIdx);
+ (*this)[elemVertIdx].checkWellPosed();
+ Valgrind::CheckDefined((*this)[elemVertIdx]);
+ }
+ }
+ };
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/boxmodels/common/boxelementsecondaryvars.hh b/dumux/boxmodels/common/boxelementsecondaryvars.hh
new file mode 100644
index 0000000000..1a1130bea7
--- /dev/null
+++ b/dumux/boxmodels/common/boxelementsecondaryvars.hh
@@ -0,0 +1,87 @@
+// $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. *
+ *****************************************************************************/
+#ifndef DUMUX_BOX_ELEMENT_SECONDARY_VARS_HH
+#define DUMUX_BOX_ELEMENT_SECONDARY_VARS_HH
+
+#include "boxproperties.hh"
+
+#include <vector>
+#include <dumux/common/boundarytypes.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup BoxModel
+ *
+ * \brief This class stores an array of SecondaryVars objects
+ */
+template<class TypeTag>
+class BoxElementSecondaryVars : public std::vector<typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) >
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) SecondaryVars;
+ typedef std::vector<SecondaryVars> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVarVector)) PrimaryVarVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ enum { dim = GridView::dimension };
+
+public:
+ /*!
+ * \brief The constructor.
+ */
+ BoxElementSecondaryVars()
+ { }
+
+ void update(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ bool oldSol)
+ {
+ const SolutionVector &globalSol =
+ oldSol?
+ problem.model().prevSol():
+ problem.model().curSol();
+ const VertexMapper &vertexMapper = problem.vertexMapper();
+ // we assert that the i-th shape function is
+ // associated to the i-th vert of the element.
+ int n = element.template count<dim>();
+ this->resize(n);
+ for (int i = 0; i < n; i++) {
+ const PrimaryVarVector &solI
+ = globalSol[vertexMapper.map(element, i, dim)];
+ (*this)[i].update(solI,
+ problem,
+ element,
+ fvElemGeom,
+ i,
+ oldSol);
+ }
+ };
+};
+
+} // namespace Dumux
+
+#endif
diff --git a/dumux/boxmodels/boxscheme/fvelementgeometry-pdelab.hh b/dumux/boxmodels/common/boxfvelementgeometry.hh
similarity index 93%
rename from dumux/boxmodels/boxscheme/fvelementgeometry-pdelab.hh
rename to dumux/boxmodels/common/boxfvelementgeometry.hh
index 69deffcd29..aaeed0b213 100644
--- a/dumux/boxmodels/boxscheme/fvelementgeometry-pdelab.hh
+++ b/dumux/boxmodels/common/boxfvelementgeometry.hh
@@ -1,4 +1,4 @@
-// $Id: fvelementgeometry-pdelab.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -13,13 +13,8 @@
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
-
-#ifndef DUMUX_FVELEMENTGEOMETRY_PDELAB_HH
-#define DUMUX_FVELEMENTGEOMETRY_PDELAB_HH
-
-#ifdef DUMUX_FVELEMENTGEOMETRY_DISC_HH
-#error "you can not use dune-disc and dune-pdelab at the same time!"
-#endif
+#ifndef DUMUX_BOX_FV_ELEMENTGEOMETRY_HH
+#define DUMUX_BOX_FV_ELEMENTGEOMETRY_HH
#include <dune/grid/common/intersectioniterator.hh>
#include <dune/grid/common/capabilities.hh>
@@ -37,7 +32,7 @@ NEW_PROP_TAG(Scalar);
/////////////////////
// HACK: Functions to initialize the subcontrol volume data
-// structures of FVElementGeomerty.
+// structures of BoxFVElementGeomerty.
//
// this is a work around to be able to to use specialization for
// doing this. it is required since it is not possible to
@@ -46,24 +41,24 @@ NEW_PROP_TAG(Scalar);
/** \todo Please doc me! */
-template <typename FVElementGeometry, int dim>
-class _FVElemGeomHelper
+template <typename BoxFVElementGeometry, int dim>
+class _BoxFVElemGeomHelper
{
public:
- static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
+ static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
- DUNE_THROW(Dune::NotImplemented, "_FVElemGeomHelper::fillSubContVolData dim = " << dim);
+ DUNE_THROW(Dune::NotImplemented, "_BoxFVElemGeomHelper::fillSubContVolData dim = " << dim);
};
};
/** \todo Please doc me! */
-template <typename FVElementGeometry>
-class _FVElemGeomHelper<FVElementGeometry, 1>
+template <typename BoxFVElementGeometry>
+class _BoxFVElemGeomHelper<BoxFVElementGeometry, 1>
{
public:
enum { dim = 1 };
- static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
+ static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
// 1D
eg.subContVol[0].volume = 0.5*eg.elementVolume;
@@ -73,13 +68,13 @@ public:
/** \todo Please doc me! */
-template <typename FVElementGeometry>
-class _FVElemGeomHelper<FVElementGeometry, 2>
+template <typename BoxFVElementGeometry>
+class _BoxFVElemGeomHelper<BoxFVElementGeometry, 2>
{
public:
enum { dim = 2 };
- static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
+ static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
switch (numVertices) {
case 3: // 2D, triangle
@@ -88,7 +83,7 @@ public:
eg.subContVol[2].volume = eg.elementVolume/3;
break;
case 4: // 2D, quadrilinear
- if (!Dune::Capabilities::IsUnstructured<typename FVElementGeometry::Grid>::v) {
+ if (!Dune::Capabilities::IsUnstructured<typename BoxFVElementGeometry::Grid>::v) {
eg.subContVol[0].volume = eg.elementVolume/4;
eg.subContVol[1].volume = eg.elementVolume/4;
eg.subContVol[2].volume = eg.elementVolume/4;
@@ -102,20 +97,20 @@ public:
}
break;
default:
- DUNE_THROW(Dune::NotImplemented, "_FVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
+ DUNE_THROW(Dune::NotImplemented, "_BoxFVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
}
}
};
/** \todo Please doc me! */
-template <typename FVElementGeometry>
-class _FVElemGeomHelper<FVElementGeometry, 3>
+template <typename BoxFVElementGeometry>
+class _BoxFVElemGeomHelper<BoxFVElementGeometry, 3>
{
public:
enum { dim = 3 };
- static void fillSubContVolData(FVElementGeometry &eg, int numVertices)
+ static void fillSubContVolData(BoxFVElementGeometry &eg, int numVertices)
{
switch (numVertices) {
case 4: // 3D, tetrahedron
@@ -278,7 +273,7 @@ public:
eg.edgeCoord[9]);
break;
default:
- DUNE_THROW(Dune::NotImplemented, "_FVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
+ DUNE_THROW(Dune::NotImplemented, "_BoxFVElemGeomHelper::fillSubContVolData dim = " << dim << ", numVertices = " << numVertices);
}
}
};
@@ -288,18 +283,18 @@ public:
/** \todo Please doc me! */
template<class TypeTag>
-class FVElementGeometry
+class BoxFVElementGeometry
{
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(LocalFEMSpace)) LocalFEMSpace;
- typedef FVElementGeometry<TypeTag> ThisType;
+ typedef BoxFVElementGeometry<TypeTag> ThisType;
/** \todo Please doc me! */
- friend class _FVElemGeomHelper<ThisType, Grid::dimension>;
+ friend class _BoxFVElemGeomHelper<ThisType, Grid::dimension>;
- typedef _FVElemGeomHelper<ThisType, Grid::dimension> FVElemGeomHelper;
+ typedef _BoxFVElemGeomHelper<ThisType, Grid::dimension> BoxFVElemGeomHelper;
enum{dim = Grid::dimension};
enum{maxNC = (dim < 3 ? 4 : 8)};
@@ -443,7 +438,7 @@ class FVElementGeometry
rightFace = edgeToFaceHex[1][k];
break;
default:
- DUNE_THROW(Dune::NotImplemented, "FVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
+ DUNE_THROW(Dune::NotImplemented, "BoxFVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
break;
}
}
@@ -525,7 +520,7 @@ class FVElementGeometry
rightEdge = faceAndVertexToRightEdgeHex[face][vert];
break;
default:
- DUNE_THROW(Dune::NotImplemented, "FVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
+ DUNE_THROW(Dune::NotImplemented, "BoxFVElementGeometry :: getFaceIndices for numVertices = " << numVertices);
break;
}
}
@@ -578,13 +573,12 @@ public:
int numFaces; //!< number of faces (0 in < 3D)
const LocalFEMSpace feMap_;
- const GridView& gridView_;
- FVElementGeometry(const GridView& gridView)
- : feMap_(), gridView_(gridView)
+ BoxFVElementGeometry()
+ : feMap_()
{}
- void update(const Element& e)
+ void update(const GridView& gridView, const Element& e)
{
const Geometry& geometry = e.geometry();
Dune::GeometryType gt = geometry.type();
@@ -623,9 +617,9 @@ public:
// fill sub control volume data use specialization for this
// \todo maybe it would be a good idea to migrate everything
// which is dependend of the grid's dimension to
- // _FVElemGeomHelper in order to benefit from more aggressive
+ // _BoxFVElemGeomHelper in order to benefit from more aggressive
// compiler optimizations...
- FVElemGeomHelper::fillSubContVolData(*this, numVertices);
+ BoxFVElemGeomHelper::fillSubContVolData(*this, numVertices);
// fill sub control volume face data:
for (int k = 0; k < numEdges; k++) { // begin loop over edges / sub control volume faces
@@ -700,8 +694,8 @@ public:
} // end loop over edges / sub control volume faces
// fill boundary face data:
- IntersectionIterator endit = gridView_.template iend(e);
- for (IntersectionIterator it = gridView_.template ibegin(e); it != endit; ++it)
+ IntersectionIterator endit = gridView.template iend(e);
+ for (IntersectionIterator it = gridView.template ibegin(e); it != endit; ++it)
if (it->boundary())
{
int face = it->indexInInside();
@@ -711,7 +705,7 @@ public:
int vertInElement = referenceElement.subEntity(face, 1, vertInFace, dim);
int bfIdx = boundaryFaceIndex(face, vertInFace);
subContVol[vertInElement].inner = false;
- switch (dim) {
+ switch ((short) dim) {
case 1:
boundaryFace[bfIdx].ipLocal = referenceElement.position(vertInElement, dim);
boundaryFace[bfIdx].area = 1.0;
@@ -735,7 +729,7 @@ public:
edgeCoord[rightEdge], faceCoord[face], edgeCoord[leftEdge]);
break;
default:
- DUNE_THROW(Dune::NotImplemented, "FVElementGeometry for dim = " << dim);
+ DUNE_THROW(Dune::NotImplemented, "BoxFVElementGeometry for dim = " << dim);
}
boundaryFace[bfIdx].ipGlobal = geometry.global(boundaryFace[bfIdx].ipLocal);
@@ -830,6 +824,5 @@ public:
}
-
#endif
diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh
new file mode 100644
index 0000000000..b1f7ce7138
--- /dev/null
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -0,0 +1,353 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2007 by Peter Bastian *
+ * Institute of Parallel and Distributed System *
+ * Department Simulation of Large Systems *
+ * University of Stuttgart, Germany *
+ * *
+ * Copyright (C) 2008-2009 by Andreas Lauser *
+ * Copyright (C) 2007-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
+ * \brief Caculates the jacobian of models based on the box scheme element-wise.
+ */
+#ifndef DUMUX_BOX_LOCAL_JACOBIAN_HH
+#define DUMUX_BOX_LOCAL_JACOBIAN_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/valgrind.hh>
+
+#include <dune/grid/common/genericreferenceelements.hh>
+
+#include <boost/format.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+#include "boxelementsecondaryvars.hh"
+#include "boxfvelementgeometry.hh"
+#include "boxlocalresidual.hh"
+
+#include "boxproperties.hh"
+
+
+
+namespace Dumux
+{
+/*!
+ * \ingroup BoxModel
+ * \brief Element-wise caculation of the jacobian matrix for models
+ * based on the box scheme .
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class BoxLocalJacobian
+{
+private:
+ typedef BoxLocalJacobian<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) LocalResidual;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GridView::Grid::ctype CoordScalar;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename Element::EntityPointer ElementPointer;
+
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
+ typedef typename RefElemProp::Container ReferenceElements;
+ typedef typename RefElemProp::ReferenceElement ReferenceElement;
+
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename Element::Geometry Geometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVarVector)) PrimaryVarVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) SecondaryVars;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSecondaryVars)) ElementSecondaryVars;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
+
+ typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
+ typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+
+public:
+ BoxLocalJacobian()
+ {
+ Valgrind::SetUndefined(problemPtr_);
+ }
+
+
+ void init(Problem &prob)
+ {
+ problemPtr_ = &prob;
+ localResidual_.init(prob);
+ // assume quadrilinears as elements with most vertices
+ A_.setSize(2<<dim, 2<<dim);
+ }
+
+ /*!
+ * \brief Assemble the linear system of equations for the
+ * verts of a element, given a local solution 'localU'.
+ */
+ void assemble(const Element &element)
+ {
+ // set the current grid element and update the element's
+ // finite volume geometry
+ elemPtr_ = &element;
+ fvElemGeom_.update(gridView_(), elem_());
+ bcTypes_.update(problem_(), elem_(), fvElemGeom_);
+
+ // this is pretty much a HACK because the internal state of
+ // the problem is not supposed to be changed during the
+ // evaluation of the residual. (Reasons: It is a violation of
+ // abstraction, makes everything more prone to errors and is
+ // not thread save.) The real solution are context objects!
+ problem_().updateCouplingParams(elem_());
+
+
+ int numVertices = fvElemGeom_.numVertices;
+
+ // update the secondary variables for the element at the last
+ // and the current time levels
+ prevSecVars_.update(problem_(),
+ elem_(),
+ fvElemGeom_,
+ true /* isOldSol? */);
+ curSecVars_.update(problem_(),
+ elem_(),
+ fvElemGeom_,
+ false /* isOldSol? */);
+
+ // calculate the local jacobian matrix
+ ElementSolutionVector partialDeriv(numVertices);
+ for (int j = 0; j < numVertices; j++) {
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++) {
+ asImp_().evalPartialDerivative_(partialDeriv,
+ j,
+ pvIdx);
+
+ // update the local stiffness matrix with the current partial
+ // derivatives
+ updateLocalJacobian_(j,
+ pvIdx,
+ partialDeriv);
+ }
+ }
+ }
+
+ /*!
+ * \brief Returns a reference to the local residual.
+ */
+ const LocalResidual &localResidual() const
+ { return localResidual_; }
+
+ /*!
+ * \brief Returns a reference to the local residual.
+ */
+ LocalResidual &localResidual()
+ { return localResidual_; }
+
+ /*!
+ * \brief Returns the jacobian of the equations at vertex i to the
+ * primary variables at vertex j.
+ */
+ const MatrixBlock &mat(int i, int j) const
+ { return A_[i][j]; }
+
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ const Problem &problem_() const
+ {
+ Valgrind::CheckDefined(problemPtr_);
+ return *problemPtr_;
+ };
+
+ /*!
+ * \brief Returns a reference to the grid view.
+ */
+ const GridView &gridView_() const
+ { return problem_().gridView(); }
+
+ /*!
+ * \brief Returns a reference to the element.
+ */
+ const Element &elem_() const
+ {
+ Valgrind::CheckDefined(elemPtr_);
+ return *elemPtr_;
+ };
+
+ /*!
+ * \brief Returns a reference to the model.
+ */
+ const Model &model_() const
+ { return problem_().model(); };
+
+ /*!
+ * \brief Returns a reference to the vertex mapper.
+ */
+ const VertexMapper &vertexMapper_() const
+ { return problem_().vertexMapper(); };
+
+ /*!
+ * \brief Compute the partial derivatives to a primary variable at
+ * an degree of freedom.
+ *
+ * This method can be overwritten by the implementation if a
+ * better scheme than central differences ought to be used.
+ */
+ void evalPartialDerivative_(ElementSolutionVector &dest,
+ int scvIdx,
+ int pvIdx)
+ {
+ int globalIdx = vertexMapper_().map(elem_(), scvIdx, dim);
+ PrimaryVarVector priVars(model_().curSol()[globalIdx]);
+ SecondaryVars origSecVars(curSecVars_[scvIdx]);
+
+ curSecVars_[scvIdx].setEvalPoint(&origSecVars);
+ Scalar eps = asImp_().numericEpsilon_(scvIdx, pvIdx);
+
+ // deflect primary variables
+ priVars[pvIdx] += eps;
+
+ // calculate the residual
+ curSecVars_[scvIdx].update(priVars,
+ problem_(),
+ elem_(),
+ fvElemGeom_,
+ scvIdx,
+ false);
+ localResidual().eval(elem_(),
+ fvElemGeom_,
+ prevSecVars_,
+ curSecVars_,
+ bcTypes_);
+
+ // store the residual
+ dest = localResidual().residual();
+
+ // deflect the primary variables
+ priVars[pvIdx] -= 2*eps;
+
+ // calculate residual again
+ curSecVars_[scvIdx].update(priVars,
+ problem_(),
+ elem_(),
+ fvElemGeom_,
+ scvIdx,
+ false);
+ localResidual().eval(elem_(),
+ fvElemGeom_,
+ prevSecVars_,
+ curSecVars_,
+ bcTypes_);
+
+ // central differences
+ dest -= localResidual().residual();
+ dest /= 2*eps;
+
+ // restore the orignal state of the element's secondary
+ // variables
+ curSecVars_[scvIdx] = origSecVars;
+
+#if HAVE_VALGRIND
+ for (unsigned i = 0; i < dest.size(); ++i)
+ Valgrind::CheckDefined(dest[i]);
+#endif
+ }
+
+ /*!
+ * \brief Returns the epsilon value which is added and removed
+ * from the current solution.
+ *
+ * \param elemSol The current solution on the element
+ * \param scvIdx The local index of the element's vertex for
+ * which the local derivative ought to be calculated.
+ * \param pvIdx The index of the primary variable which gets varied
+ */
+ Scalar numericEpsilon_(int scvIdx,
+ int pvIdx) const
+ {
+ Scalar pv = this->curSecVars_[scvIdx].primaryVars()[pvIdx];
+ return 1e-9*(std::abs(pv) + 1);
+ }
+
+ /*!
+ * \brief Updates the current local stiffness matrix with the
+ * partial derivatives of all equations in regard to the
+ * primary variable 'pvIdx' at vertex j .
+ */
+ void updateLocalJacobian_(int scvIdx,
+ int pvIdx,
+ const ElementSolutionVector &stiffness)
+ {
+ for (int i = 0; i < fvElemGeom_.numVertices; i++) {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++) {
+ // A[i][scvIdx][eqIdx][pvIdx] is the approximate rate
+ // of change of the residual of equation 'eqIdx' at
+ // vertex 'i' depending on the primary variable
+ // 'pvIdx' at vertex 'scvIdx'.
+ this->A_[i][scvIdx][eqIdx][pvIdx] = stiffness[i][eqIdx];
+ Valgrind::CheckDefined(this->A_[i][scvIdx][eqIdx][pvIdx]);
+ }
+ }
+ }
+
+ const Element *elemPtr_;
+
+ FVElementGeometry fvElemGeom_;
+ ElementBoundaryTypes bcTypes_;
+
+ // The problem we would like to solve
+ Problem *problemPtr_;
+
+ // secondary variables at the previous and at the current time
+ // levels
+ ElementSecondaryVars prevSecVars_;
+ ElementSecondaryVars curSecVars_;
+
+ LocalResidual localResidual_;
+ LocalBlockMatrix A_;
+};
+}
+
+#endif
diff --git a/dumux/boxmodels/common/boxlocalresidual.hh b/dumux/boxmodels/common/boxlocalresidual.hh
new file mode 100644
index 0000000000..bb8354ca16
--- /dev/null
+++ b/dumux/boxmodels/common/boxlocalresidual.hh
@@ -0,0 +1,456 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2007 by Peter Bastian *
+ * Institute of Parallel and Distributed System *
+ * Department Simulation of Large Systems *
+ * University of Stuttgart, Germany *
+ * *
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2007-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
+ * \brief Caculates the residual of models based on the box scheme element-wise.
+ */
+#ifndef DUMUX_BOX_LOCAL_RESIDUAL_HH
+#define DUMUX_BOX_LOCAL_RESIDUAL_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/valgrind.hh>
+#include <dune/grid/common/genericreferenceelements.hh>
+
+#include <boost/format.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+#include "boxproperties.hh"
+#include "boxfvelementgeometry.hh"
+#include "boxelementboundarytypes.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup BoxModel
+ * \brief Element-wise caculation of the residual matrix for models
+ * based on the box scheme .
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class BoxLocalResidual
+{
+private:
+ typedef BoxLocalResidual<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GridView::Grid::ctype CoordScalar;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
+ typedef typename RefElemProp::Container ReferenceElements;
+ typedef typename RefElemProp::ReferenceElement ReferenceElement;
+
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename Element::Geometry Geometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVarVector)) PrimaryVarVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) SecondaryVars;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSecondaryVars)) ElementSecondaryVars;
+
+ typedef Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
+ typedef Dune::Matrix<MatrixBlock> LocalBlockMatrix;
+
+public:
+ BoxLocalResidual()
+ { }
+
+ ~BoxLocalResidual()
+ { }
+
+ void init(Problem &prob)
+ { problemPtr_ = &prob; }
+
+ /*!
+ * \brief Compute the local residual, i.e. the deviation of the
+ * equations from zero.
+ */
+ void eval(const Element &element)
+ {
+ FVElementGeometry fvGeom;
+ fvGeom.update(gridView_(), element);
+ ElementSecondaryVars secVarsPrev, secVarsCur;
+ secVarsPrev.update(problem_(),
+ element,
+ fvGeom,
+ true /* oldSol? */);
+ secVarsCur.update(problem_(),
+ element,
+ fvGeom,
+ false /* oldSol? */);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(problem_(), element, fvGeom);
+
+ // this is pretty much a HACK because the internal state of
+ // the problem is not supposed to be changed during the
+ // evaluation of the residual. (Reasons: It is a violation of
+ // abstraction, makes everything more prone to errors and is
+ // not thread save.) The real solution are context objects!
+ problem_().updateCouplingParams(element);
+
+ asImp_().eval(element, fvGeom, secVarsPrev, secVarsCur, bcTypes);
+ }
+
+ /*!
+ * \brief Compute the flux term for the current solution.
+ */
+ void evalFluxes(const Element &element,
+ const ElementSecondaryVars &curSecVars)
+ {
+ FVElementGeometry fvGeom;
+ fvGeom.update(gridView_(), element);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(problem_(), element, fvGeom);
+
+ residual_.resize(fvGeom.numVertices);
+ residual_ = 0;
+
+ elemPtr_ = &element;
+ fvElemGeomPtr_ = &fvGeom;
+ bcTypesPtr_ = &bcTypes;
+ prevSecVarsPtr_ = 0;
+ curSecVarsPtr_ = &curSecVars;
+ asImp_().evalFluxes_();
+ }
+
+
+ /*!
+ * \brief Compute the local residual, i.e. the deviation of the
+ * equations from zero.
+ */
+ void eval(const Element &element,
+ const FVElementGeometry &fvGeom,
+ const ElementSecondaryVars &prevSecVars,
+ const ElementSecondaryVars &curSecVars,
+ const ElementBoundaryTypes &bcTypes)
+ {
+ //Valgrind::CheckDefined(fvGeom);
+ Valgrind::CheckDefined(prevSecVars);
+ Valgrind::CheckDefined(curSecVars);
+
+#if HAVE_VALGRIND
+ for (int i=0; i < fvGeom.numVertices; i++) {
+ Valgrind::CheckDefined(prevSecVars[i]);
+ Valgrind::CheckDefined(curSecVars[i]);
+ }
+#endif // HAVE_VALGRIND
+
+ elemPtr_ = &element;
+ fvElemGeomPtr_ = &fvGeom;
+ bcTypesPtr_ = &bcTypes;
+ prevSecVarsPtr_ = &prevSecVars;
+ curSecVarsPtr_ = &curSecVars;
+
+ // reset residual
+ int numVerts = fvElemGeom_().numVertices;
+ residual_.resize(numVerts);
+ residual_ = 0;
+
+ asImp_().evalFluxes_();
+ asImp_().evalVolumeTerms_();
+
+ // evaluate the boundary
+ asImp_().evalBoundary_();
+
+#if HAVE_VALGRIND
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ Valgrind::CheckDefined(residual_[i]);
+#endif // HAVE_VALGRIND
+ }
+
+ /*!
+ * \brief Returns the local residual for a given sub-control
+ * volume of the element.
+ */
+ const ElementSolutionVector &residual() const
+ { return residual_; }
+
+ /*!
+ * \brief Returns the local residual for a given sub-control
+ * volume of the element.
+ */
+ const PrimaryVarVector residual(int scvIdx) const
+ { return residual_[scvIdx]; }
+
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ void evalBoundary_()
+ {
+ Dune::GeometryType geoType = elem_().geometry().type();
+ const ReferenceElement &refElem = ReferenceElements::general(geoType);
+
+ // evaluate boundary conditions for all intersections of
+ // the current element
+ IntersectionIterator isIt = gridView_().template ibegin(elem_());
+ const IntersectionIterator &endIt = gridView_().template iend(elem_());
+ for (; isIt != endIt; ++isIt)
+ {
+ // handle only faces on the boundary
+ if (!isIt->boundary())
+ continue;
+
+ // Assemble the boundary for all vertices of the current
+ // face
+ int faceIdx = isIt->indexInInside();
+ int numFaceVerts = refElem.size(faceIdx, 1, dim);
+ for (int faceVertIdx = 0;
+ faceVertIdx < numFaceVerts;
+ ++faceVertIdx)
+ {
+ int elemVertIdx = refElem.subEntity(faceIdx,
+ 1,
+ faceVertIdx,
+ dim);
+
+ int boundaryFaceIdx =
+ fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
+
+ // add the neuman fluxes of a single boundary segment
+ evalBoundarySegment_(isIt,
+ elemVertIdx,
+ boundaryFaceIdx);
+ }
+ }
+ }
+
+ // handle boundary conditions for a single sub-control volume face
+ void evalBoundarySegment_(const IntersectionIterator &isIt,
+ int scvIdx,
+ int boundaryFaceIdx)
+ {
+ // temporary vector to store the neumann boundary fluxes
+ PrimaryVarVector values(0.0);
+ const BoundaryTypes &bcTypes = bcTypes_(scvIdx);
+
+ // deal with neumann boundaries
+ if (bcTypes.hasNeumann()) {
+ problem_().neumann(values,
+ elem_(),
+ fvElemGeom_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx);
+ values *= fvElemGeom_().boundaryFace[boundaryFaceIdx].area;
+ Valgrind::CheckDefined(values);
+ residual_[scvIdx] += values;
+ }
+
+ // deal with dirichlet boundaries
+ if (bcTypes.hasDirichlet()) {
+ problem_().dirichlet(values,
+ elem_(),
+ fvElemGeom_(),
+ *isIt,
+ scvIdx,
+ boundaryFaceIdx);
+
+ Valgrind::CheckDefined(values);
+ for (int i = 0; i < numEq; ++i) {
+ if (!bcTypes.isDirichlet(i))
+ continue;
+
+ int pvIdx = bcTypes.eqToDirichletIndex(i);
+ residual_[scvIdx][i] =
+ curPrimaryVars_(scvIdx)[pvIdx] - values[pvIdx];
+ }
+ }
+ }
+
+ void evalFluxes_()
+ {
+ // calculate the mass flux over the faces and subtract
+ // it from the local rates
+ for (int k = 0; k < fvElemGeom_().numEdges; k++)
+ {
+ int i = fvElemGeom_().subContVolFace[k].i;
+ int j = fvElemGeom_().subContVolFace[k].j;
+
+ PrimaryVarVector flux;
+ Valgrind::SetUndefined(flux);
+ this->asImp_().computeFlux(flux, k);
+ Valgrind::CheckDefined(flux);
+
+ // subtract fluxes from the local mass rates of the
+ // respective sub control volume adjacent to the face.
+ for (int eq = 0; eq < numEq; ++ eq) {
+ residual_[i][eq] -= flux[eq];
+ residual_[j][eq] += flux[eq];
+ }
+ }
+ }
+
+ void evalVolumeTerms_()
+ {
+ // evaluate the volume terms (storage + source terms)
+ for (int i=0; i < fvElemGeom_().numVertices; i++)
+ {
+ PrimaryVarVector massContrib(0), tmp(0);
+
+ // mass balance within the element. this is the
+ // $\frac{m}{\partial t}$ term if using implicit
+ // euler as time discretization.
+ //
+ // TODO (?): we might need a more explicit way for
+ // doing the time discretization...
+ this->asImp_().computeStorage(massContrib, i, false);
+ this->asImp_().computeStorage(tmp, i, true);
+
+ massContrib -= tmp;
+ massContrib *=
+ fvElemGeom_().subContVol[i].volume
+ /
+ problem_().timeManager().timeStepSize();
+
+ for (int j = 0; j < numEq; ++j)
+ residual_[i][j] += massContrib[j];
+
+ // subtract the source term from the local rate
+ PrimaryVarVector source;
+ this->asImp_().computeSource(source, i);
+ source *= fvElemGeom_().subContVol[i].volume;
+
+ for (int j = 0; j < numEq; ++j) {
+ residual_[i][j] -= source[j];
+
+ // make sure that only defined quantities where used
+ // to calculate the residual.
+ Valgrind::CheckDefined(residual_[i][j]);
+ }
+ }
+ }
+
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ const Problem &problem_() const
+ { return *problemPtr_; };
+
+ /*!
+ * \brief Returns a reference to the model.
+ */
+ const Model &model_() const
+ { return problem_().model(); };
+
+ /*!
+ * \brief Returns a reference to the vertex mapper.
+ */
+ const VertexMapper &vertexMapper_() const
+ { return problem_().vertexMapper(); };
+ const GridView &gridView_() const
+ { return problem_().gridView(); }
+
+ const Element &elem_() const
+ {
+ Valgrind::CheckDefined(elemPtr_);
+ return *elemPtr_;
+ }
+
+ const FVElementGeometry &fvElemGeom_() const
+ {
+ Valgrind::CheckDefined(fvElemGeomPtr_);
+ return *fvElemGeomPtr_;
+ }
+
+ const PrimaryVarVector &curPrimaryVars_(int i) const
+ {
+ return curSecVars_(i).primaryVars();
+ }
+
+ const ElementSecondaryVars &curSecVars_() const
+ {
+ Valgrind::CheckDefined(curSecVarsPtr_);
+ return *curSecVarsPtr_;
+ }
+ const SecondaryVars &curSecVars_(int i) const
+ {
+ return curSecVars_()[i];
+ }
+
+ const ElementSecondaryVars &prevSecVars_() const
+ {
+ Valgrind::CheckDefined(prevSecVarsPtr_);
+ return *prevSecVarsPtr_;
+ }
+ const SecondaryVars &prevSecVars_(int i) const
+ {
+ return prevSecVars_()[i];
+ }
+
+ const ElementBoundaryTypes &bcTypes_() const
+ {
+ Valgrind::CheckDefined(bcTypesPtr_);
+ return *bcTypesPtr_;
+ }
+ const BoundaryTypes &bcTypes_(int i) const
+ {
+ return bcTypes_()[i];
+ }
+
+protected:
+ ElementSolutionVector residual_;
+
+private:
+ // The problem we would like to solve
+ Problem *problemPtr_;
+
+ const Element *elemPtr_;
+ const FVElementGeometry *fvElemGeomPtr_;
+
+ // current and previous secondary variables for the element
+ const ElementSecondaryVars *prevSecVarsPtr_;
+ const ElementSecondaryVars *curSecVarsPtr_;
+
+ const ElementBoundaryTypes *bcTypesPtr_;
+};
+}
+
+#endif
diff --git a/dumux/boxmodels/common/boxmodel.hh b/dumux/boxmodels/common/boxmodel.hh
new file mode 100644
index 0000000000..9866385779
--- /dev/null
+++ b/dumux/boxmodels/common/boxmodel.hh
@@ -0,0 +1,587 @@
+// $Id$
+/*****************************************************************************
+ * Copyright (C) 2007 by Peter Bastian *
+ * Institute of Parallel and Distributed System *
+ * Department Simulation of Large Systems *
+ * University of Stuttgart, Germany *
+ * *
+ * Copyright (C) 2008 by Andreas Lauser, 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_BOX_SCHEME_HH
+#define DUMUX_BOX_SCHEME_HH
+
+#include <dumux/common/valgrind.hh>
+#include <dune/grid/common/genericreferenceelements.hh>
+
+#include <boost/format.hpp>
+
+#include "boxproperties.hh"
+
+#include "boxelementsecondaryvars.hh"
+#include "boxlocaljacobian.hh"
+#include "boxlocalresidual.hh"
+
+#include "pdelabboxlocaloperator.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \defgroup BoxModel Box-Scheme
+ */
+/*!
+ * \ingroup BoxModel
+ * \defgroup BoxProblems Box-Problems
+ */
+/*!
+ * \ingroup BoxModel
+ * \defgroup BoxModels Box-Models
+ */
+
+
+/*!
+ * \ingroup BoxModel
+ *
+ * \brief The base class for the vertex centered finite volume
+ * discretization scheme.
+ */
+template<class TypeTag>
+class BoxModel
+{
+ typedef BoxModel<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GridView::Grid::ctype CoordScalar;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(DofMapper)) DofMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVarVector)) PrimaryVarVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
+
+ enum {
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+ dim = GridView::dimension
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
+ typedef typename RefElemProp::Container ReferenceElements;
+ typedef typename RefElemProp::ReferenceElement ReferenceElement;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) LocalResidual;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SecondaryVars)) SecondaryVars;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+
+public:
+ /*!
+ * \brief The constructor.
+ */
+ BoxModel()
+ { }
+
+ ~BoxModel()
+ { delete jacAsm_; }
+
+ /*!
+ * \brief Apply the initial conditions to the model.
+ */
+ void init(Problem &prob)
+ {
+ problemPtr_ = &prob;
+
+ int nDofs = asImp_().numDofs();
+ uCur_.resize(nDofs);
+ uPrev_.resize(nDofs);
+ boxVolume_.resize(nDofs);
+
+ localJacobian_.init(problem_());
+ jacAsm_ = new JacobianAssembler();
+ jacAsm_->init(problem_());
+
+ applyInitialSolution_();
+
+ // also set the solution of the "previous" time step to the
+ // initial solution.
+ uPrev_ = uCur_;
+ }
+
+ /*!
+ * \brief Compute the global residual for an arbitrary solution
+ * vector.
+ */
+ Scalar globalResidual(SolutionVector &dest, const SolutionVector &u)
+ {
+ SolutionVector tmp(curSol());
+ curSol() = u;
+ Scalar res = globalResidual(dest);
+ curSol() = tmp;
+ return res;
+ }
+
+ /*!
+ * \brief Compute the global residual for the current solution
+ * vector.
+ */
+ Scalar globalResidual(SolutionVector &dest)
+ {
+ dest = 0;
+
+ ElementIterator elemIt = gridView_().template begin<0>();
+ const ElementIterator elemEndIt = gridView_().template end<0>();
+ for (; elemIt != elemEndIt; ++elemIt) {
+ localResidual().eval(*elemIt);
+
+ for (int i = 0; i < elemIt->template count<dim>(); ++i) {
+ int globalI = vertexMapper().map(*elemIt, i, dim);
+ dest[globalI] += localResidual().residual(i);
+ }
+ };
+
+ Scalar result = dest.two_norm();
+ /*
+ Scalar result = 0;
+ for (int i = 0; i < (*tmp).size(); ++i) {
+ for (int j = 0; j < numEq; ++j)
+ result += std::abs((*tmp)[i][j]);
+ }
+ */
+ return result;
+ }
+
+ /*!
+ * \brief Returns the volume of a given control volume.
+ */
+ Scalar boxVolume(int globalIdx) const
+ { return boxVolume_[globalIdx][0]; }
+
+ /*!
+ * \brief Reference to the current solution as a block vector.
+ */
+ const SolutionVector &curSol() const
+ { return uCur_; }
+
+ /*!
+ * \brief Reference to the current solution as a block vector.
+ */
+ SolutionVector &curSol()
+ { return uCur_; }
+
+ /*!
+ * \brief Reference to the previous solution as a block vector.
+ */
+ const SolutionVector &prevSol() const
+ { return uPrev_; }
+
+ /*!
+ * \brief Reference to the previous solution as a block vector.
+ */
+ SolutionVector &prevSol()
+ { return uPrev_; }
+
+ /*!
+ * \brief Returns the operator assembler for the global jacobian of
+ * the problem.
+ */
+ JacobianAssembler &jacobianAssembler()
+ { return *jacAsm_; }
+
+ /*!
+ * \brief Returns the local jacobian which calculates the local
+ * stiffness matrix for an arbitrary element.
+ *
+ * The local stiffness matrices of the element are used by
+ * the jacobian assembler to produce a global linerization of the
+ * problem.
+ */
+ LocalJacobian &localJacobian()
+ { return localJacobian_; }
+ /*!
+ * \copydoc localJacobian()
+ */
+ const LocalJacobian &localJacobian() const
+ { return localJacobian_; }
+
+ /*!
+ * \brief Returns the local residual function.
+ */
+ LocalResidual &localResidual()
+ { return localJacobian().localResidual(); }
+ /*!
+ * \copydoc localResidual()
+ */
+ const LocalResidual &localResidual() const
+ { return localJacobian().localResidual(); }
+
+ /*!
+ * \brief Try to progress the model to the next timestep.
+ */
+ bool update(NewtonMethod &solver,
+ NewtonController &controller)
+ {
+#if HAVE_VALGRIND
+ for (size_t i = 0; i < curSol().size(); ++i)
+ Valgrind::CheckDefined(curSol()[i]);
+#endif // HAVE_VALGRIND
+
+ asImp_().updateBegin();
+
+ bool converged = solver.execute(controller);
+ if (converged)
+ asImp_().updateSuccessful();
+ else
+ asImp_().updateFailed();
+
+#if HAVE_VALGRIND
+ for (size_t i = 0; i < curSol().size(); ++i) {
+ Valgrind::CheckDefined(curSol()[i]);
+ }
+#endif // HAVE_VALGRIND
+
+ return converged;
+ }
+
+
+ /*!
+ * \brief Called by the update() method before it tries to
+ * apply the newton method. This is primary a hook
+ * which the actual model can overload.
+ */
+ void updateBegin()
+ { }
+
+
+ /*!
+ * \brief Called by the update() method if it was
+ * successful. This is primary a hook which the actual
+ * model can overload.
+ */
+ void updateSuccessful()
+ {
+ // make the current solution the previous one.
+ uPrev_ = uCur_;
+ };
+
+ /*!
+ * \brief Called by the update() method if it was
+ * unsuccessful. This is primary a hook which the actual
+ * model can overload.
+ */
+ void updateFailed()
+ {
+ // Reset the current solution to the one of the
+ // previous time step so that we can start the next
+ // update at a physically meaningful solution.
+ uCur_ = uPrev_;
+ };
+
+ /*!
+ * \brief Serializes the current state of the model.
+ */
+ template <class Restarter>
+ void serialize(Restarter &res)
+ { res.template serializeEntities<dim>(asImp_(), this->gridView_()); }
+
+ /*!
+ * \brief Deserializes the state of the model.
+ */
+ template <class Restarter>
+ void deserialize(Restarter &res)
+ {
+ res.template deserializeEntities<dim>(asImp_(), this->gridView_());
+ prevSol() = curSol();
+ }
+
+ /*!
+ * \brief Write the current solution for a vertex to a restart
+ * file.
+ */
+ void serializeEntity(std::ostream &outstream,
+ const Vertex &vert)
+ {
+ int vertIdx = dofMapper().map(vert);
+
+ // write phase state
+ if (!outstream.good()) {
+ DUNE_THROW(Dune::IOError,
+ "Could not serialize vertex "
+ << vertIdx);
+ }
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+ outstream << curSol()[vertIdx][eqIdx] << " ";
+ }
+ };
+
+ /*!
+ * \brief Reads the current solution variables for a vertex from a
+ * restart file.
+ */
+ void deserializeEntity(std::istream &instream,
+ const Vertex &vert)
+ {
+ int vertIdx = dofMapper().map(vert);
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+ if (!instream.good())
+ DUNE_THROW(Dune::IOError,
+ "Could not deserialize vertex "
+ << vertIdx);
+ instream >> curSol()[vertIdx][eqIdx];
+ }
+ };
+
+ /*!
+ * \brief Returns the number of global degrees of freedoms (DOFs)
+ */
+ size_t numDofs() const
+ { return gridView_().size(dim); }
+
+ /*!
+ * \brief Mapper for the entities where degrees of freedoms are
+ * defined to indices.
+ *
+ * This usually means a mapper for vertices.
+ */
+ const DofMapper &dofMapper() const
+ { return problem_().vertexMapper(); };
+
+ /*!
+ * \brief Mapper for vertices to indices.
+ */
+ const VertexMapper &vertexMapper() const
+ { return problem_().vertexMapper(); };
+
+ /*!
+ * \brief Mapper for elements to indices.
+ */
+ const ElementMapper &elementMapper() const
+ { return problem_().elementMapper(); };
+
+ void resetJacobianAssembler ()
+ {
+ delete jacAsm_;
+ jacAsm_ = new JacobianAssembler;
+ jacAsm_->init(problem_());
+ }
+
+ /*!
+ * \brief Add the vector fields for analysing the convergence of
+ * the newton method to the a VTK multi writer.
+ *
+ * \param writer The VTK multi writer where the fields should be added.
+ * \param update The delte of the solution function before and after the Newton update
+ */
+ template <class MultiWriter>
+ void addConvergenceVtkFields(MultiWriter &writer,
+ const SolutionVector &u,
+ const SolutionVector &deltaU)
+ {
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+
+ SolutionVector globalResid(asImp_(), 0.0);
+ asImp_().globalResidual(globalResid, u);
+
+ // create the required scalar fields
+ unsigned numVertices = this->gridView_().size(dim);
+ //unsigned numElements = this->gridView_().size(0);
+
+ // global defect of the two auxiliary equations
+ ScalarField* def[numEq];
+ ScalarField* delta[numEq];
+ ScalarField* x[numEq];
+ for (int i = 0; i < numEq; ++i) {
+ x[i] = writer.template createField<Scalar, 1>(numVertices);
+ delta[i] = writer.template createField<Scalar, 1>(numVertices);
+ def[i] = writer.template createField<Scalar, 1>(numVertices);
+ }
+
+ VertexIterator vIt = this->gridView_().template begin<dim>();
+ VertexIterator vEndIt = this->gridView_().template end<dim>();
+ for (; vIt != vEndIt; ++ vIt)
+ {
+ int globalIdx = vertexMapper().map(*vIt);
+ for (int i = 0; i < numEq; ++i) {
+ (*x[i])[globalIdx] = u[globalIdx][i];
+ (*delta[i])[globalIdx] = - deltaU[globalIdx][i];
+ (*def[i])[globalIdx] = globalResid[globalIdx][i];
+ }
+ }
+
+ for (int i = 0; i < numEq; ++i) {
+ writer.addVertexData(x[i], (boost::format("x_%i")%i).str().c_str());
+ writer.addVertexData(delta[i], (boost::format("delta_%i")%i).str().c_str());
+ writer.addVertexData(def[i], (boost::format("defect_%i")%i).str().c_str());
+ }
+
+ asImp_().addOutputVtkFields(u, writer);
+ }
+
+ /*!
+ * \brief Add the quantities of a time step which ought to be
+ * written to disk.
+ *
+ * This should be overwritten by the acutal model if any secondary
+ * variables should be written out. Read: This should _always_ be
+ * overwritten by well behaved models!
+ *
+ * \param writer The VTK multi writer where the fields should be added.
+ */
+ template <class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol,
+ MultiWriter &writer)
+ {
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+
+ // create the required scalar fields
+ unsigned numVertices = this->gridView_().size(dim);
+
+ // global defect of the two auxiliary equations
+ ScalarField* x[numEq];
+ for (int i = 0; i < numEq; ++i) {
+ x[i] = writer.template createField<Scalar, 1>(numVertices);
+ }
+
+ VertexIterator vIt = this->gridView_().template begin<dim>();
+ VertexIterator vEndIt = this->gridView_().template end<dim>();
+ for (; vIt != vEndIt; ++ vIt)
+ {
+ int globalIdx = vertexMapper().map(*vIt);
+ for (int i = 0; i < numEq; ++i) {
+ (*x[i])[globalIdx] = sol[globalIdx][i];
+ }
+ }
+
+ for (int i = 0; i < numEq; ++i)
+ writer.addVertexData(x[i], (boost::format("primaryVar%i")%i).str().c_str());
+ }
+
+
+protected:
+ /*!
+ * \brief A reference to the problem on which the model is applied.
+ */
+ Problem &problem_()
+ { return *problemPtr_; }
+ /*!
+ * \copydoc problem_()
+ */
+ const Problem &problem_() const
+ { return *problemPtr_; }
+
+ /*!
+ * \brief Reference to the grid view of the spatial domain.
+ */
+ const GridView &gridView_() const
+ { return problem_().gridView(); }
+
+ LocalResidual &localResidual_()
+ { return localJacobian_.localResidual(); }
+
+ void applyInitialSolution_()
+ {
+ // first set the whole domain to zero
+ uCur_ = Scalar(0.0);
+ boxVolume_ = Scalar(0.0);
+
+ FVElementGeometry fvElemGeom;
+
+ // iterate through leaf grid and evaluate initial
+ // condition at the center of each sub control volume
+ //
+ // TODO: the initial condition needs to be unique for
+ // each vertex. we should think about the API...
+ ElementIterator it = gridView_().template begin<0>();
+ const ElementIterator &eendit = gridView_().template end<0>();
+ for (; it != eendit; ++it) {
+ // deal with the current element
+ fvElemGeom.update(gridView_(), *it);
+
+ // loop over all element vertices, i.e. sub control volumes
+ for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; scvIdx++)
+ {
+ // map the local vertex index to the global one
+ int globalIdx = vertexMapper().map(*it,
+ scvIdx,
+ dim);
+
+ // let the problem do the dirty work of nailing down
+ // the initial solution.
+ PrimaryVarVector initVal;
+ Valgrind::SetUndefined(initVal);
+ problem_().initial(initVal,
+ *it,
+ fvElemGeom,
+ scvIdx);
+ Valgrind::CheckDefined(initVal);
+
+ // add up the initial values of all sub-control
+ // volumes. If the initial values disagree for
+ // different sub control volumes, the initial value
+ // will be the arithmetic mean.
+ initVal *= fvElemGeom.subContVol[scvIdx].volume;
+ boxVolume_[globalIdx] += fvElemGeom.subContVol[scvIdx].volume;
+ uCur_[globalIdx] += initVal;
+ Valgrind::CheckDefined(uCur_[globalIdx]);
+ }
+ }
+ // divide all primary variables by the volume of their boxes
+ int n = gridView_().size(dim);
+ for (int i = 0; i < n; ++i) {
+ uCur_[i] /= boxVolume(i);
+ }
+ }
+
+ bool verbose_() const
+ { return gridView_().comm().rank() == 0; };
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ // the problem we want to solve. defines the constitutive
+ // relations, matxerial laws, etc.
+ Problem *problemPtr_;
+
+ // calculates the local jacobian matrix for a given element
+ LocalJacobian localJacobian_;
+ // Linearizes the problem at the current time step using the
+ // local jacobian
+ JacobianAssembler *jacAsm_;
+
+ // cur is the current iterative solution, prev the converged
+ // solution of the previous time step
+ SolutionVector uCur_;
+ SolutionVector uPrev_;
+
+ Dune::BlockVector<Dune::FieldVector<Scalar, 1> > boxVolume_;
+};
+}
+
+#endif
diff --git a/dumux/boxmodels/boxscheme/boxproblem.hh b/dumux/boxmodels/common/boxproblem.hh
similarity index 60%
rename from dumux/boxmodels/boxscheme/boxproblem.hh
rename to dumux/boxmodels/common/boxproblem.hh
index c2d52822c0..80dc594172 100644
--- a/dumux/boxmodels/boxscheme/boxproblem.hh
+++ b/dumux/boxmodels/common/boxproblem.hh
@@ -1,4 +1,4 @@
-// $Id: boxproblem.hh 3784 2010-06-24 13:43:57Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -30,43 +30,50 @@
namespace Dumux
{
/*!
- * \ingroup BoxScheme
- * \brief Base class for all problems which use the box scheme
+ * \ingroup BoxModel
+ * \brief Base class for all problems which use the box scheme
*
* \todo Please doc me more!
*/
-template<class TypeTag, class Implementation>
+template<class TypeTag>
class BoxProblem
{
private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- enum Episode {}; // the type of an episode of the simulation
- typedef Dumux::TimeManager<Episode> TimeManager;
+ typedef Dumux::VtkMultiWriter<GridView> VtkMultiWriter;
- typedef Dumux::VtkMultiWriter<GridView> VtkMultiWriter;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
- typedef typename GridView::Grid::ctype CoordScalar;
- typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
- typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ typedef typename GridView::Grid::ctype CoordScalar;
+ typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
public:
- BoxProblem(const GridView &gridView)
+ BoxProblem(TimeManager &timeManager, const GridView &gridView)
: gridView_(gridView),
bboxMin_(std::numeric_limits<double>::max()),
bboxMax_(-std::numeric_limits<double>::max()),
- resultWriter_(asImp_()->name())
+ elementMapper_(gridView),
+ vertexMapper_(gridView),
+ timeManager_(&timeManager),
+ resultWriter_(asImp_().name())
{
wasRestarted_ = false;
@@ -83,8 +90,9 @@ public:
bboxMin_[i] = gridView.comm().min(bboxMin_[i]);
bboxMax_[i] = gridView.comm().max(bboxMax_[i]);
}
- model_ = new Model(*asImp_());
- newtonMethod_ = new NewtonMethod(*model_);
+
+ model_ = new Model();
+ newtonMethod_ = new NewtonMethod(asImp_());
}
~BoxProblem()
@@ -93,28 +101,6 @@ public:
delete newtonMethod_;
};
- /*!
- * \name Simulation steering
- */
- // \{
-
- /*!
- * \brief Start the simulation procedure.
- *
- * This method is usually called by the main() function and simply
- * uses Dumux::TimeManager::runSimulation() to do the actual
- * work.
- */
- bool simulate(Scalar dtInitial, Scalar tEnd)
- {
- // set the initial time step and the time where the simulation ends
- timeManager_.setEndTime(tEnd);
- timeManager_.setTimeStepSize(dtInitial);
- timeManager_.runSimulation(*asImp_());
- return true;
- };
-
-
/*!
* \brief Called by the Dumux::TimeManager in order to
* initialize the problem.
@@ -122,16 +108,28 @@ public:
void init()
{
// set the initial condition of the model
- model().initial();
-
- // write the inital solution to disk
- asImp_()->writeCurrentResult_();
+ model().init(asImp_());
}
+ /*!
+ * \brief If model coupling is used, this updates the parameters
+ * required to calculate the coupling fluxes between the
+ * sub-models.
+ *
+ * By default it does nothing
+ */
+ void updateCouplingParams(const Element &element) const
+ {}
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
/*!
* \brief Called by the time manager before the time integration.
*/
- void timeStepBegin()
+ void preProcess()
{}
/*!
@@ -140,7 +138,26 @@ public:
*/
void timeIntegration()
{
- model().update(*newtonMethod_, newtonCtl_);
+ const int maxFails = 10;
+ for (int i = 0; i < maxFails; ++i) {
+ if (i > 0 && gridView().comm().rank() == 0)
+ std::cout << "Newton solver did not converge. Retrying with time step of "
+ << timeManager().timeStepSize() << "sec\n";
+
+ if (model_->update(*newtonMethod_, newtonCtl_))
+ return;
+
+ // update failed
+ Scalar dt = timeManager().timeStepSize();
+ Scalar nextDt = dt / 2;
+ timeManager().setTimeStepSize(nextDt);
+ }
+
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxFails
+ << " timestep divisions. dt="
+ << timeManager().timeStepSize());
}
/*!
@@ -150,39 +167,10 @@ public:
*/
Scalar nextTimeStepSize()
{
- Scalar dt = asImp_()->timeManager().timeStepSize();
+ Scalar dt = asImp_().timeManager().timeStepSize();
return newtonCtl_.suggestTimeStepSize(dt);
};
-
- /*!
- * \brief This method is called by the model if the update to the
- * next time step failed completely.
- */
- void updateSuccessful()
- {
- wasRestarted_ = false;
- asImp_()->writeCurrentResult_();
- };
-
- /*!
- * \brief This method is called by the model if the update to the
- * next time step failed completely.
- */
- void updateFailed()
- { };
-
- /*!
- * \brief This method is called by the model if the update to the
- * next time step failed with the curent time step size.
- */
- void updateFailedTry()
- {
- Scalar dt = asImp_()->timeManager().timeStepSize();
- Scalar nextDt = newtonCtl_.suggestTimeStepSize(dt);
- asImp_()->timeManager().setTimeStepSize(nextDt);
- };
-
/*!
* \brief Returns true if a restart file should be written to
* disk.
@@ -191,11 +179,11 @@ public:
* steps. This file is intented to be overwritten by the
* implementation.
*/
- bool shouldWriteRestartFile() const
+ bool doSerialize() const
{
return !restarted() &&
- timeManager().timeStepNum() > 0 &&
- (timeManager().timeStepNum() % 5 == 0);
+ timeManager().timeStepIndex() > 0 &&
+ (timeManager().timeStepIndex() % 10 == 0);
}
/*!
@@ -206,18 +194,14 @@ public:
* very time step. This file is intented to be overwritten by the
* implementation.
*/
- bool shouldWriteOutputFile() const
- { return !restarted(); }
+ bool doOutput() const
+ { return true; }
/*!
* \brief Called by the time manager after the time integration.
*/
- void timeStepEnd()
- {
- // write restart file if necessary
- if (asImp_()->shouldWriteRestartFile())
- serialize();
- }
+ void postProcess()
+ { }
// \}
@@ -265,17 +249,29 @@ public:
const GlobalPosition &bboxMax() const
{ return bboxMax_; }
+ /*!
+ * \brief Returns the mapper for vertices to indices.
+ */
+ const VertexMapper &vertexMapper() const
+ { return vertexMapper_; }
+
+ /*!
+ * \brief Returns the mapper for elements to indices.
+ */
+ const ElementMapper &elementMapper() const
+ { return elementMapper_; }
+
/*!
* \brief Returns TimeManager object used by the simulation
*/
TimeManager &timeManager()
- { return timeManager_; }
+ { return *timeManager_; }
/*!
* \copydoc timeManager()
*/
const TimeManager &timeManager() const
- { return timeManager_; }
+ { return *timeManager_; }
/*!
* \brief Returns numerical model used for the problem.
@@ -303,7 +299,7 @@ public:
{ return wasRestarted_; }
/*!
- * \brief This method writes the complete state of the problem
+ * \brief This method writes the complete state of the simulation
* to the harddisk.
*
* The file will start with the prefix returned by the name()
@@ -313,19 +309,47 @@ public:
*/
void serialize()
{
- typedef Dumux::Restart<GridView> Restarter;
-
+ typedef Dumux::Restart Restarter;
Restarter res;
- res.serializeBegin(gridView(),
- asImp_()->name(),
- timeManager_.time());
- std::cerr << "Serialize to file " << res.fileName() << "\n";
+ res.serializeBegin(asImp_());
+ std::cerr << "Serialize to file '" << res.fileName() << "'\n";
- timeManager_.serialize(res);
+ timeManager().serialize(res);
+ asImp_().serialize(res);
+ res.serializeEnd();
+ }
+
+ /*!
+ * \brief This method writes the complete state of the problem
+ * to the harddisk.
+ *
+ * The file will start with the prefix returned by the name()
+ * method, has the current time of the simulation clock in it's
+ * name and uses the extension <tt>.drs</tt>. (Dumux ReStart
+ * file.) See Dumux::Restart for details.
+ */
+ template <class Restarter>
+ void serialize(Restarter &res)
+ {
resultWriter_.serialize(res);
model().serialize(res);
+ }
- res.serializeEnd();
+ /*!
+ * \brief Load a previously saved state of the whole simulation
+ * from disk.
+ */
+ void restart(Scalar tRestart)
+ {
+ typedef Dumux::Restart Restarter;
+
+ Restarter res;
+
+ res.deserializeBegin(asImp_(), tRestart);
+ std::cout << "Deserialize from file '" << res.fileName() << "'\n";
+ timeManager().deserialize(res);
+ asImp_().deserialize(res);
+ res.deserializeEnd();
}
/*!
@@ -334,77 +358,72 @@ public:
*
* It is the inverse of the serialize() method.
*/
- void deserialize(double t)
+ template <class Restarter>
+ void deserialize(Restarter &res)
{
- typedef Dumux::Restart<GridView> Restarter;
-
- Restarter res;
- res.deserializeBegin(gridView(),
- asImp_()->name(),
- t);
- std::cerr << "Deserialize from file " << res.fileName() << "\n";
-
- timeManager_.deserialize(res);
resultWriter_.deserialize(res);
model().deserialize(res);
-
- res.deserializeEnd();
-
+
wasRestarted_ = true;
};
// \}
-protected:
- //! Returns the implementation of the problem (i.e. static polymorphism)
- Implementation *asImp_()
- { return static_cast<Implementation *>(this); }
-
- //! \copydoc asImp_()
- const Implementation *asImp_() const
- { return static_cast<const Implementation *>(this); }
-
- //! Write the fields current solution into an VTK output file.
- void writeCurrentResult_()
+ /*!
+ * \brief Write the relavant secondar variables of the current
+ * solution into an VTK output file.
+ */
+ void writeOutput()
{
// write the current result to disk
- if (asImp_()->shouldWriteOutputFile()) {
+ if (asImp_().doOutput()) {
if (gridView().comm().rank() == 0)
- std::cout << "Writing result file for current time step\n";
+ std::cout << "Writing result file for \"" << asImp_().name() << "\"\n";
// calculate the time _after_ the time was updated
- Scalar t = timeManager_.time() + timeManager_.timeStepSize();
- resultWriter_.beginTimestep(t,
- gridView());
- model().addOutputVtkFields(resultWriter_);
+ Scalar t = timeManager().time() + timeManager().timeStepSize();
+ resultWriter_.beginTimestep(t, gridView());
+ model().addOutputVtkFields(model().curSol(), resultWriter_);
resultWriter_.endTimestep();
}
}
+protected:
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation &asImp_()
+ { return *static_cast<Implementation *>(this); }
+
+ //! \copydoc asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation *>(this); }
+
private:
static std::string simname_; // a string for the name of the current simulation,
// which could be set by means of an program argument,
// for example.
- const GridView gridView_;
+ const GridView gridView_;
+
+ GlobalPosition bboxMin_;
+ GlobalPosition bboxMax_;
- GlobalPosition bboxMin_;
- GlobalPosition bboxMax_;
+ ElementMapper elementMapper_;
+ VertexMapper vertexMapper_;
- TimeManager timeManager_;
+ TimeManager *timeManager_;
- Model *model_;
+ Model *model_;
NewtonMethod *newtonMethod_;
NewtonController newtonCtl_;
- VtkMultiWriter resultWriter_;
+ VtkMultiWriter resultWriter_;
bool wasRestarted_;
};
// definition of the static class member simname_,
// which is necessary because it is of type string.
-template <class TypeTag, class Implementation>
-std::string BoxProblem<TypeTag, Implementation>::simname_="sim"; //initialized with default "sim"
+template <class TypeTag>
+std::string BoxProblem<TypeTag>::simname_="sim"; //initialized with default "sim"
}
diff --git a/dumux/boxmodels/common/boxproperties.hh b/dumux/boxmodels/common/boxproperties.hh
new file mode 100644
index 0000000000..eb84ea8966
--- /dev/null
+++ b/dumux/boxmodels/common/boxproperties.hh
@@ -0,0 +1,361 @@
+// $Id$
+/*****************************************************************************
+ * 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_BOX_PROPERTIES_HH
+#define DUMUX_BOX_PROPERTIES_HH
+
+#include <dumux/common/propertysystem.hh>
+
+/*!
+ * \file
+ * \brief Specify the shape functions, operator assemblers, etc
+ * used for the BoxModel.
+ */
+namespace Dumux
+{
+
+namespace Properties
+{
+/*!
+ * \addtogroup BoxModel
+ */
+// \{
+
+//////////////////////////////////////////////////////////////////
+// Type tags tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for models based on the box-scheme
+NEW_TYPE_TAG(BoxModel);
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+//!< Property tag for scalar vaslues
+NEW_PROP_TAG(Scalar);
+
+NEW_PROP_TAG(Grid); //!< The type of the DUNE grid
+NEW_PROP_TAG(GridView); //!< The type of the grid view
+
+NEW_PROP_TAG(ReferenceElements); //!< DUNE reference elements to be used
+NEW_PROP_TAG(FVElementGeometry); //! The type of the finite-volume geometry in the box scheme
+
+NEW_PROP_TAG(Problem); //!< The type of the problem
+NEW_PROP_TAG(Model); //!< The type of the discretization
+NEW_PROP_TAG(NumEq); //!< Number of equations in the system of PDEs
+NEW_PROP_TAG(LocalResidual); //!< The type of the local residual function
+NEW_PROP_TAG(LocalJacobian); //!< The type of the local jacobian operator
+
+NEW_PROP_TAG(JacobianAssembler); //!< Assembles the global jacobian matrix
+NEW_PROP_TAG(JacobianMatrix); //!< Type of the global jacobian matrix
+NEW_PROP_TAG(BoundaryTypes); //!< Stores the boundary types of a single degree of freedom
+NEW_PROP_TAG(ElementBoundaryTypes); //!< Stores the boundary types on an element
+
+NEW_PROP_TAG(PrimaryVarVector); //!< A vector of primary variables within a sub-control volume
+NEW_PROP_TAG(SolutionVector); //!< Vector containing all primary variable vector of the grid
+NEW_PROP_TAG(ElementSolutionVector); //!< A vector of primary variables within a sub-control volume
+
+NEW_PROP_TAG(SecondaryVars); //!< The secondary variables within a sub-control volume
+NEW_PROP_TAG(ElementSecondaryVars); //!< The secondary variables of all sub-control volumes in an element
+NEW_PROP_TAG(FluxVars); //!< Data required to calculate a flux over a face
+
+// high level simulation control
+NEW_PROP_TAG(TimeManager); //!< Manages the simulation time
+NEW_PROP_TAG(NewtonMethod); //!< The type of the newton method
+NEW_PROP_TAG(NewtonController); //!< The type of the newton controller
+
+// properties for the PDELab wrapper
+NEW_PROP_TAG(LocalFEMSpace); //!< The local finite element space used for the finite element interpolation
+NEW_PROP_TAG(ScalarGridFunctionSpace); //!< The used grid function space for a single finite element function
+NEW_PROP_TAG(GridFunctionSpace); //!< The used grid function space
+NEW_PROP_TAG(Constraints); //!< The constraints on the grid function space
+NEW_PROP_TAG(ConstraintsTrafo); //!< The type of PDELab's constraints transformation
+NEW_PROP_TAG(LocalOperator); //!< The type of the local operator used by PDELab
+NEW_PROP_TAG(GridOperatorSpace); //!< The used grid operator space
+
+// mappers from local to global indices
+NEW_PROP_TAG(VertexMapper);
+NEW_PROP_TAG(ElementMapper);
+NEW_PROP_TAG(DofMapper);
+}
+}
+
+
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+#include "boxfvelementgeometry.hh"
+#include "pdelabboxassembler.hh"
+#include "pdelabboxistlvectorbackend.hh"
+//#include <dumux/boxmodels/pdelab/boxdirichletconstraints.hh>
+
+#include <dumux/common/boundarytypes.hh>
+#include <dumux/common/timemanager.hh>
+
+namespace Dumux {
+
+template<typename TypeTag>
+class BoxFVElementGeometry;
+
+template<typename TypeTag>
+class BoxElementBoundaryTypes;
+
+template<typename TypeTag>
+class BoxElementSecondaryVars;
+
+namespace PDELab {
+template<typename TypeTag>
+class BoxLocalOperator;
+
+template<typename TypeTag>
+class BoxAssembler;
+
+template<typename TypeTag>
+class BoxISTLVectorBackend;
+};
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Some defaults for very fundamental properties
+//////////////////////////////////////////////////////////////////
+
+//! Set the default type for scalar values to double
+SET_PROP_DEFAULT(Scalar)
+{ typedef double type; };
+
+//! Set the default type for the time manager
+SET_PROP_DEFAULT(TimeManager)
+{ typedef Dumux::TimeManager<TypeTag> type; };
+
+/*!
+ * \brief Specify the reference elements which we ought to use.
+ *
+ * We use Dune::ReferenceElements by default (-> old entity
+ * numbering).
+ *
+ * TODO: Some specialization if the grid only supports one kind of
+ * cells would be nice. this would be better fixed inside DUNE,
+ * though. something like:
+ * Dune::GenericReferenceElements<Dune::GeometryType<cube, dim> >
+ */
+SET_PROP_DEFAULT(ReferenceElements)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+
+ typedef typename Grid::ctype CoordScalar;
+ static const int dim = Grid::dimension;
+
+public:
+ typedef Dune::GenericReferenceElements<CoordScalar, dim> Container;
+ typedef Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
+};
+
+//////////////////////////////////////////////////////////////////
+// Properties
+//////////////////////////////////////////////////////////////////
+
+//! Use the leaf grid view if not defined otherwise
+SET_PROP(BoxModel, GridView)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
+
+public:
+ typedef typename Grid::LeafGridView type;
+};
+
+//! Set the default for the FVElementGeometry
+SET_PROP(BoxModel, FVElementGeometry)
+{
+ typedef Dumux::BoxFVElementGeometry<TypeTag> type;
+};
+
+//! Set the default for the ElementBoundaryTypes
+SET_PROP(BoxModel, ElementBoundaryTypes)
+{ typedef Dumux::BoxElementBoundaryTypes<TypeTag> type; };
+
+//! use the plain newton method for the box scheme by default
+SET_PROP(BoxModel, NewtonMethod)
+{public:
+ typedef Dumux::NewtonMethod<TypeTag> type;
+};
+
+//! use the plain newton controller for the box scheme by default
+SET_PROP(BoxModel, NewtonController)
+{public:
+ typedef Dumux::NewtonController<TypeTag> type;
+};
+
+//! Mapper for the grid view's vertices.
+SET_PROP(BoxModel, VertexMapper)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ template<int dim>
+ struct VertexLayout {
+ bool contains (Dune::GeometryType gt) const
+ { return gt.dim() == 0; }
+ };
+public:
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, VertexLayout> type;
+};
+
+//! Mapper for the grid view's elements.
+SET_PROP(BoxModel, ElementMapper)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ template<int dim>
+ struct ElementLayout {
+ bool contains (Dune::GeometryType gt) const
+ { return gt.dim() == dim; }
+ };
+public:
+ typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, ElementLayout> type;
+};
+
+//! Mapper for the degrees of freedoms.
+SET_PROP(BoxModel, DofMapper)
+{ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) type; };
+
+
+/*!
+ * \brief The type of a solution for the whole grid at a fixed time.
+ */
+SET_PROP(BoxModel, SolutionVector)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+public:
+ typedef typename GridFunctionSpace::template VectorContainer<Scalar>::Type type;
+};
+
+/*!
+ * \brief The type of a solution for a whole element.
+ */
+SET_PROP(BoxModel, ElementSolutionVector)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVarVector)) PrimaryVarVector;
+public:
+ typedef Dune::BlockVector<PrimaryVarVector> type;
+};
+
+/*!
+ * \brief A vector of primary variables.
+ */
+SET_PROP(BoxModel, PrimaryVarVector)
+{ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector))::block_type type; };
+
+/*!
+ * \brief An array of secondary variable containers.
+ */
+SET_TYPE_PROP(BoxModel, ElementSecondaryVars, Dumux::BoxElementSecondaryVars<TypeTag>);
+
+/*!
+ * \brief Boundary types at a single degree of freedom.
+ */
+SET_PROP(BoxModel, BoundaryTypes)
+{ private:
+ enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
+public:
+ typedef Dumux::BoundaryTypes<numEq> type;
+};
+
+/*!
+ * \brief Assembler for the global jacobian matrix.
+ */
+SET_TYPE_PROP(BoxModel, JacobianAssembler, Dumux::PDELab::BoxAssembler<TypeTag>);
+
+//! Extract the type of a global jacobian matrix from the solution types
+SET_PROP(BoxModel, JacobianMatrix)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
+public:
+ typedef typename GridOperatorSpace::template MatrixContainer<Scalar>::Type type;
+};
+
+//! use the local FEM space associated with cubes by default
+SET_PROP(BoxModel, 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;
+};
+
+SET_PROP(BoxModel, GridFunctionSpace)
+{private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
+public:
+ //typedef MyBoxConstraints Constraints;
+ typedef Dune::PDELab::NoConstraints Constraints;
+
+ typedef Dune::PDELab::GridFunctionSpace<GridView, FEM, Constraints, Dumux::PDELab::BoxISTLVectorBackend<TypeTag> >
+ ScalarGridFunctionSpace;
+
+ typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, numEq, Dune::PDELab::GridFunctionSpaceBlockwiseMapper>
+ type;
+
+ typedef typename type::template ConstraintsContainer<Scalar>::Type
+ ConstraintsTrafo;
+};
+
+SET_PROP(BoxModel, ConstraintsTrafo)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::ConstraintsTrafo type; };
+SET_PROP(BoxModel, Constraints)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::Constraints type; };
+SET_PROP(BoxModel, ScalarGridFunctionSpace)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridFunctionSpace))::ScalarGridFunctionSpace type; };
+
+SET_PROP(BoxModel, GridOperatorSpace)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+ enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
+
+public:
+ typedef Dumux::PDELab::BoxLocalOperator<TypeTag> LocalOperator;
+
+ typedef Dune::PDELab::GridOperatorSpace<GridFunctionSpace,
+ GridFunctionSpace,
+ LocalOperator,
+ ConstraintsTrafo,
+ ConstraintsTrafo,
+ Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
+ true
+ > type;
+};
+
+
+SET_PROP(BoxModel, LocalOperator)
+{ typedef typename GET_PROP(TypeTag, PTAG(GridOperatorSpace))::LocalOperator type; };
+
+// \}
+
+}
+}
+
+#endif
diff --git a/dumux/boxmodels/pdelab/assemblerpdelab.hh b/dumux/boxmodels/common/pdelabboxassembler.hh
similarity index 76%
rename from dumux/boxmodels/pdelab/assemblerpdelab.hh
rename to dumux/boxmodels/common/pdelabboxassembler.hh
index 8ae27ddfa4..77c00e7187 100644
--- a/dumux/boxmodels/pdelab/assemblerpdelab.hh
+++ b/dumux/boxmodels/common/pdelabboxassembler.hh
@@ -1,4 +1,4 @@
-// $Id: assemblerpdelab.hh 3759 2010-06-21 16:59:10Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009-2010 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
@@ -13,26 +13,28 @@
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
-#ifdef HAVE_DUNE_PDELAB
-
-#ifndef DUMUX_ASSEMBLERPDELAB_HH
-#define DUMUX_ASSEMBLERPDELAB_HH
+#ifndef DUMUX_PDELAB_BOX_ASSEMBLER_HH
+#define DUMUX_PDELAB_BOX_ASSEMBLER_HH
#include<dune/pdelab/finiteelementmap/p1fem.hh>
#include<dune/pdelab/finiteelementmap/q1fem.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
+#include<dune/pdelab/function/const.hh>
#include<dune/pdelab/finiteelementmap/conformingconstraints.hh>
#include<dune/pdelab/backend/istlvectorbackend.hh>
#include<dune/pdelab/backend/istlmatrixbackend.hh>
-#include"boundarytypespdelab.hh"
-#include"boxjacobianpdelab.hh"
+//#include "pdelabboundarytypes.hh"
+#include "pdelabboxlocaloperator.hh"
+
+namespace Dumux {
+namespace PDELab {
+
template<class TypeTag>
-class AssemblerPDELab
+class BoxAssembler
{
-public:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
@@ -40,66 +42,73 @@ public:
enum{dim = GridView::dimension};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
- typedef typename GET_PROP(TypeTag, PTAG(PDELabTypes)) PDELabTypes;
- typedef typename PDELabTypes::Constraints Constraints;
- typedef typename PDELabTypes::ScalarGridFunctionSpace ScalarGridFunctionSpace;
- typedef typename PDELabTypes::GridFunctionSpace GridFunctionSpace;
- typedef typename PDELabTypes::ConstraintsTrafo ConstraintsTrafo;
- typedef typename PDELabTypes::LocalOperator LocalOperator;
- typedef typename PDELabTypes::GridOperatorSpace GridOperatorSpace;
- typedef BoundaryIndexHelperPDELab<TypeTag> BoundaryFunction;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Constraints)) Constraints;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ScalarGridFunctionSpace)) ScalarGridFunctionSpace;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalOperator)) LocalOperator;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
+
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::SolutionVector SolutionVector;
- typedef typename GridFunctionSpace::template VectorContainer<Scalar>::Type Vector;
- typedef typename GridOperatorSpace::template MatrixContainer<Scalar>::Type Matrix;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
+public:
+
+ typedef SolutionVector Vector;
+ typedef JacobianMatrix Matrix;
typedef Matrix RepresentationType;
- AssemblerPDELab(Model &model, Problem& problem)
- : problem_(problem)
+ BoxAssembler()
{
+ problem_ = 0;
fem_ = 0;
cn_ = 0;
scalarGridFunctionSpace_ = 0;
gridFunctionSpace_ = 0;
- bTypes_ = 0;
constraintsTrafo_ = 0;
localOperator_ = 0;
gridOperatorSpace_ = 0;
matrix_ = 0;
-
- fem_ = new FEM();
- cn_ = new Constraints(problem_);
- scalarGridFunctionSpace_ = new ScalarGridFunctionSpace(problem_.gridView(), *fem_, *cn_);
- gridFunctionSpace_ = new GridFunctionSpace(*scalarGridFunctionSpace_);
-
- cn_->compute_ghosts(*gridFunctionSpace_);
-
- bTypes_ = new BoundaryFunction();
- constraintsTrafo_ = new ConstraintsTrafo();
- Dune::PDELab::constraints(*bTypes_, *gridFunctionSpace_, *constraintsTrafo_, false);
-
- localOperator_ = new LocalOperator(model);
- gridOperatorSpace_ = new GridOperatorSpace(*gridFunctionSpace_, *constraintsTrafo_,
- *gridFunctionSpace_, *constraintsTrafo_, *localOperator_);
-
- matrix_ = new Matrix(*gridOperatorSpace_);
- *matrix_ = 0;
}
- ~AssemblerPDELab()
+ ~BoxAssembler()
{
delete matrix_;
delete gridOperatorSpace_;
delete localOperator_;
delete constraintsTrafo_;
- delete bTypes_;
delete gridFunctionSpace_;
delete scalarGridFunctionSpace_;
delete cn_;
delete fem_;
}
+ void init(Problem& problem)
+ {
+ problem_ = &problem;
+ fem_ = new FEM();
+ //cn_ = new Constraints(*problem_);
+ cn_ = new Constraints();
+ scalarGridFunctionSpace_ = new ScalarGridFunctionSpace(problem_->gridView(), *fem_, *cn_);
+ gridFunctionSpace_ = new GridFunctionSpace(*scalarGridFunctionSpace_);
+
+ //cn_->compute_ghosts(*gridFunctionSpace_);
+
+ //typedef BoundaryIndexHelper<TypeTag> BoundaryFunction;
+ //BoundaryFunction *bTypes = new BoundaryFunction();
+ constraintsTrafo_ = new ConstraintsTrafo();
+ //Dune::PDELab::constraints(*bTypes, *gridFunctionSpace_, *constraintsTrafo_, false);
+
+ localOperator_ = new LocalOperator(problem_->model());
+ gridOperatorSpace_ =
+ new GridOperatorSpace(*gridFunctionSpace_, *constraintsTrafo_,
+ *gridFunctionSpace_, *constraintsTrafo_, *localOperator_);
+
+ matrix_ = new Matrix(*gridOperatorSpace_);
+ *matrix_ = 0;
+ }
+
void assemble(SolutionVector &u)
{
*matrix_ = 0;
@@ -109,9 +118,6 @@ public:
residual_ = 0;
gridOperatorSpace_->residual(u, residual_);
- set_constrained_dofs(*constraintsTrafo_, 0.0, residual_);
- set_constrained_dofs(*constraintsTrafo_, 0.0, u);
-
#if 0
// rescale jacobian and right hand side to the largest
// entry on the main diagonal block matrix
@@ -188,12 +194,11 @@ public:
{ return residual_; }
private:
- Problem& problem_;
+ Problem *problem_;
Constraints *cn_;
FEM *fem_;
ScalarGridFunctionSpace *scalarGridFunctionSpace_;
GridFunctionSpace *gridFunctionSpace_;
- BoundaryFunction *bTypes_;
ConstraintsTrafo *constraintsTrafo_;
LocalOperator *localOperator_;
GridOperatorSpace *gridOperatorSpace_;
@@ -202,6 +207,7 @@ private:
SolutionVector residual_;
};
-#endif
+} // namespace PDELab
+} // namespace Dumux
#endif
diff --git a/dumux/boxmodels/pdelab/istlvectorbackend.hh b/dumux/boxmodels/common/pdelabboxistlvectorbackend.hh
similarity index 97%
rename from dumux/boxmodels/pdelab/istlvectorbackend.hh
rename to dumux/boxmodels/common/pdelabboxistlvectorbackend.hh
index dfb039705c..20fd2e9375 100644
--- a/dumux/boxmodels/pdelab/istlvectorbackend.hh
+++ b/dumux/boxmodels/common/pdelabboxistlvectorbackend.hh
@@ -22,7 +22,7 @@ namespace PDELab {
//! ISTL backend for FunctionSpace
template<class TypeTag>
-class ISTLVectorBackend
+class BoxISTLVectorBackend
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
@@ -37,7 +37,7 @@ public:
public:
typedef E ElementType;
typedef Dune::BlockVector< Dune::FieldVector<E,BlockSize> > BaseT;
- typedef ISTLVectorBackend<TypeTag> Backend;
+ typedef BoxISTLVectorBackend<TypeTag> Backend;
VectorContainer()
{}
diff --git a/dumux/boxmodels/pdelab/boxjacobianpdelab.hh b/dumux/boxmodels/common/pdelabboxlocaloperator.hh
similarity index 68%
rename from dumux/boxmodels/pdelab/boxjacobianpdelab.hh
rename to dumux/boxmodels/common/pdelabboxlocaloperator.hh
index ff8dc7b0d8..5e3a4b2aeb 100644
--- a/dumux/boxmodels/pdelab/boxjacobianpdelab.hh
+++ b/dumux/boxmodels/common/pdelabboxlocaloperator.hh
@@ -1,4 +1,4 @@
-// $Id: boxjacobianpdelab.hh 3736 2010-06-15 09:52:10Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009-2010 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
@@ -13,10 +13,8 @@
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
-#ifdef HAVE_DUNE_PDELAB
-
-#ifndef DUMUX_BOXJACOBIANPDELAB_HH
-#define DUMUX_BOXJACOBIANPDELAB_HH
+#ifndef DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
+#define DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
#include<vector>
#include<dune/common/fvector.hh>
@@ -27,11 +25,15 @@
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
+namespace Dumux {
+
+namespace PDELab {
+
template<class TypeTag>
-class BoxJacobianPDELab
+class BoxLocalOperator
:
-// : public Dune::PDELab::NumericalJacobianApplyVolume<BoxJacobianPDELab<TypeTag> >,
-//public Dune::PDELab::NumericalJacobianVolume<BoxJacobianPDELab<TypeTag> >,
+// : public Dune::PDELab::NumericalJacobianApplyVolume<BoxLocalOperatorPDELab<TypeTag> >,
+//public Dune::PDELab::NumericalJacobianVolume<BoxLocalOperatorPDELab<TypeTag> >,
public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags
{
@@ -42,39 +44,27 @@ public:
// residual assembly flags
enum { doAlphaVolume = true };
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::SolutionOnElement SolutionOnElement;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
- BoxJacobianPDELab (Model &model)
+ BoxLocalOperator(Model &model)
: model_(model)
{}
// volume integral depending on test and ansatz functions
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
- const LFSV& lfsv, R& r) const
+ const LFSV& lfsv, R& r) const
{
typedef typename LFSU::Traits::SizeType size_type;
- model_.localJacobian().setCurrentElement(eg.entity());
-
+ model_.localResidual().eval(eg.entity());
+
int numVertices = x.size()/numEq;
-
- SolutionOnElement localU(numVertices);
- model_.localJacobian().restrictToElement(localU, model_.curSol());
- model_.localJacobian().setCurrentSolution(localU);
-
- SolutionOnElement localUOld(numVertices);
- model_.localJacobian().restrictToElement(localUOld, model_.prevSol());
- model_.localJacobian().setPreviousSolution(localUOld);
-
- SolutionOnElement localResidual(numVertices);
- model_.localJacobian().evalLocalResidual(localResidual);
for (size_type comp = 0; comp < r.size(); comp++)
- r[comp] = localResidual[comp%numVertices][comp/numVertices];
+ r[comp] = model_.localResidual().residual(comp%numVertices)[comp/numVertices];
}
// jacobian of volume term
@@ -99,6 +89,7 @@ private:
Model& model_;
};
-#endif
+} // namespace PDELab
+} // namespace Dumux
#endif
diff --git a/dumux/boxmodels/pdelab/Makefile.am b/dumux/boxmodels/pdelab/Makefile.am
deleted file mode 100644
index 86f8ad6565..0000000000
--- a/dumux/boxmodels/pdelab/Makefile.am
+++ /dev/null
@@ -1,4 +0,0 @@
-pdelabdir = $(includedir)/dumux/boxmodels/pdelab
-pdelab_HEADERS = *.hh
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/boxmodels/pdelab/boundarytypespdelab.hh b/dumux/boxmodels/pdelab/boundarytypespdelab.hh
deleted file mode 100644
index ce75d911e9..0000000000
--- a/dumux/boxmodels/pdelab/boundarytypespdelab.hh
+++ /dev/null
@@ -1,49 +0,0 @@
-// $Id: boundarytypespdelab.hh 3736 2010-06-15 09:52:10Z lauser $
-/*****************************************************************************
- * Copyright (C) 2009-2010 by Bernd Flemisch, 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. *
- *****************************************************************************/
-#if HAVE_DUNE_PDELAB && !defined DUMUX_BOUNDARYTYPESPDELAB_HH
-#define DUMUX_BOUNDARYTYPESPDELAB_HH
-
-#include<dune/common/fvector.hh>
-#include<dune/pdelab/common/function.hh>
-
-template <class TypeTag>
-class BoundaryIndexHelperPDELab
-{
- enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
-public:
- template <int i>
- struct Child {
- struct Type {
- enum { isLeaf = true };
- enum { eqIdx = i };
- };
- };
-
- template <int i>
- const typename Child<i>::Type &getChild() const
- {
- static typename Child<i>::Type dummy;
- return dummy;
- };
-
- enum { isLeaf = false };
- enum { CHILDREN = numEq };
-
- BoundaryIndexHelperPDELab()
- {}
-};
-
-#endif
diff --git a/dumux/boxmodels/pdelab/boxdirichletconstraints.hh b/dumux/boxmodels/pdelab/boxdirichletconstraints.hh
deleted file mode 100644
index 6b8d6f0e1c..0000000000
--- a/dumux/boxmodels/pdelab/boxdirichletconstraints.hh
+++ /dev/null
@@ -1,235 +0,0 @@
-// $Id: boxdirichletconstraints.hh 3834 2010-07-14 12:50:32Z bernd $
-// -*- tab-width: 4; indent-tabs-mode: nil -*-
-/*****************************************************************************
- * Copyright (C) 2009-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. *
- *****************************************************************************/
-#ifndef DUMUX_BOXDIRICHLETCONSTRAINTS_HH
-#define DUMUX_BOXDIRICHLETCONSTRAINTS_HH
-
-#include <cstddef>
-
-#include<dune/common/exceptions.hh>
-#include<dune/grid/common/genericreferenceelements.hh>
-#include<dune/grid/common/grid.hh>
-#include<dune/common/geometrytype.hh>
-#include<dune/pdelab/common/geometrywrapper.hh>
-#include<dune/pdelab/finiteelementmap/conformingconstraints.hh>
-
-#include<dumux/common/boundarytypes.hh>
-#include<dumux/common/boundaryconditions.hh>
-
-namespace Dumux {
-//! Constraints construction
-// works in any dimension and on all element types
-template <class TypeTag>
-class BoxDirichletConstraints // : public Dune::PDELab::ConformingDirichletConstraints
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
-
- enum {numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
- enum {dim = GridView::dimension};
-
- typedef Dumux::BoundaryTypes<numEq> BoundaryTypeVector;
- Problem &problem_;
-
-public:
- enum { doBoundary = true };
- enum { doProcessor = false };
- enum { doSkeleton = false };
- enum { doVolume = false };
-
- BoxDirichletConstraints(Problem& problem)
- : problem_(problem)
- {}
-
- //! boundary constraints
- /**
- * \tparam F grid function returning boundary condition type
- * \tparam IG intersection geometry
- * \tparam LFS local function space
- * \tparam T TransformationType
- */
- template<typename F, typename I, typename LFS,
- typename T>
- void boundary (const F& f,
- const Dune::PDELab::IntersectionGeometry<I>& ig,
- const LFS& lfs,
- T& trafo) const
- {
- const ElementPointer& elementPointer = ig.inside();
- FVElementGeometry fvElemGeom(problem_.gridView());
-
- fvElemGeom.update(*elementPointer);
- BoundaryTypeVector bcTypes;
-
- //problem_.boundaryTypes(values, element, fvElemGeom, ig.intersection(), scvIdx, boundaryFaceIdx);
- //typename F::Traits::RangeType bctype;
-
- const int face = ig.indexInInside();
-
- // find all local indices of this face
- Dune::GeometryType gt = ig.inside()->type();
- typedef typename Dune::PDELab::IntersectionGeometry<I>::ctype DT;
- const int dim = Dune::PDELab::IntersectionGeometry<I>::Entity::Geometry::dimension;
- const Dune::GenericReferenceElement<DT,dim>& refelem = Dune::GenericReferenceElements<DT,dim>::general(gt);
-
- // empty map means Dirichlet constraint
- typename T::RowType empty;
-
- for (int faceVertIdx = 0; faceVertIdx < refelem.size(face, 1, dim); faceVertIdx ++){
- int elemVertIdx = refelem.subEntity(face, 1, faceVertIdx, dim);
- int boundaryFaceIdx = fvElemGeom.boundaryFaceIndex(face, faceVertIdx);
-
- bcTypes.reset();
- problem_.boundaryTypes(bcTypes, *elementPointer, fvElemGeom, ig.intersection(), elemVertIdx, boundaryFaceIdx);
- bcTypes.checkWellPosed();
-
- for (std::size_t i = 0; i < lfs.localFiniteElement().localCoefficients().size(); i++) {
- // The codim to which this dof is attached to
- unsigned int codim = lfs.localFiniteElement().localCoefficients().localKey(i).codim();
-
- if (codim!=dim)
- continue;
- if (lfs.localFiniteElement().localCoefficients().localKey(i).subEntity() != elemVertIdx)
- continue;
-
- if (bcTypes.isDirichlet(F::eqIdx)) {
- trafo[i] = empty;
- }
- }
- }
- }
-};
-
-// extend constraints class by processor boundary
-template <class TypeTag>
-class NonoverlappingBoxDirichletConstraints : public BoxDirichletConstraints<TypeTag>
-{
-public:
- enum { doVolume = true };
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
-
- NonoverlappingBoxDirichletConstraints(Problem& problem)
- : BoxDirichletConstraints<TypeTag>(problem)
- {}
-
-
- template<typename E, typename LFS, typename T>
- void volume (const Dune::PDELab::ElementGeometry<E>& eg, const LFS& lfs, T& trafo) const
- {
- // nothing to do for interior entities
- if (eg.entity().partitionType()==Dune::InteriorEntity)
- return;
-
- // empty map means Dirichlet constraint
- typename T::RowType empty;
-
- typedef typename LFS::Traits::GridFunctionSpaceType::Traits::BackendType B;
-
- // loop over all degrees of freedom and check if it is not owned by this processor
- for (size_t i=0; i<lfs.localFiniteElement().localCoefficients().size(); i++)
- {
- if (ghost_[lfs.globalIndex(i)]!=0)
- {
- trafo[i] = empty;
- }
- }
- }
-
- template<class GFS>
- void compute_ghosts (const GFS& gfs)
- {
- typedef typename GFS::template VectorContainer<int>::Type V;
- V ighost(gfs);
- Dune::PDELab::GhostDataHandle<GFS,V> gdh(gfs,ighost);
- if (gfs.gridview().comm().size()>1)
- gfs.gridview().communicate(gdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
- ighost.std_copy_to(ghost_);
- rank = gfs.gridview().comm().rank();
- }
-
- void print ()
- {
- std::cout << "/" << rank << "/ " << "ghost size="
- << ghost_.size() << std::endl;
- for (std::size_t i=0; i<ghost_.size(); i++)
- std::cout << "/" << rank << "/ " << "ghost[" << i << "]="
- << ghost_[i] << std::endl;
- }
-
-private:
- int rank;
- std::vector<int> ghost_;
-};
-
-// extend constraints class by processor boundary
-template <class TypeTag>
-class OverlappingBoxDirichletConstraints : public BoxDirichletConstraints<TypeTag>
-{
-public:
- enum { doProcessor = true };
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
-
- OverlappingBoxDirichletConstraints(Problem& problem)
- : BoxDirichletConstraints<TypeTag>(problem)
- {}
-
- // boundary constraints
- // IG : intersection geometry
- // LFS : local function space
- // T : TransformationType
- template<typename I, typename LFS, typename T>
- void processor (const Dune::PDELab::IntersectionGeometry<I>& ig,
- const LFS& lfs, T& trafo) const
- {
- // determine face
- const int face = ig.indexInInside();
-
- // find all local indices of this face
- Dune::GeometryType gt = ig.inside()->type();
- typedef typename Dune::PDELab::IntersectionGeometry<I>::ctype DT;
- const int dim = Dune::PDELab::IntersectionGeometry<I>::Entity::Geometry::dimension;
-
-
- const Dune::GenericReferenceElement<DT,dim>& refelem = Dune::GenericReferenceElements<DT,dim>::general(gt);
-
- // empty map means Dirichlet constraint
- typename T::RowType empty;
-
- // loop over all degrees of freedom and check if it is on given face
- for (size_t i=0; i<lfs.localFiniteElement().localCoefficients().size(); i++)
- {
- // The codim to which this dof is attached to
- unsigned int codim = lfs.localFiniteElement().localCoefficients().localKey(i).codim();
-
- if (codim==0) continue;
-
- for (int j=0; j<refelem.size(face,1,codim); j++)
- if (lfs.localFiniteElement().localCoefficients().localKey(i).subEntity()==refelem.subEntity(face,1,j,codim))
- trafo[i] = empty;
- }
- }
-};
-
-
-
-
-}
-
-#endif
diff --git a/dumux/common/boundarytypes.hh b/dumux/common/boundarytypes.hh
index ac29c8af89..520346b69c 100644
--- a/dumux/common/boundarytypes.hh
+++ b/dumux/common/boundarytypes.hh
@@ -1,4 +1,4 @@
-// $Id: boundarytypes.hh 3834 2010-07-14 12:50:32Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -26,8 +26,6 @@
#include <dune/common/mpihelper.hh>
#include <dumux/common/valgrind.hh>
-#include <dumux/common/boundaryconditions.hh>
-#include <dumux/common/boundaryconditions.hh>
namespace Dumux
{
@@ -35,42 +33,6 @@ namespace Dumux
template <int numEq>
class BoundaryTypes
{
- class BoundaryWrapper
- { public:
- BoundaryWrapper(BoundaryTypes *bt, int idx)
- : bt_(bt), idx_(idx) {};
-
- BoundaryWrapper &operator=(int bc) {
- switch (bc) {
- case Dumux::BoundaryConditions::neumann:
- bt_->boundaryInfo_[idx_].visited = 1;
- bt_->boundaryInfo_[idx_].isDirichlet = 0;
-
- Valgrind::SetDefined(bt_->boundaryInfo_[idx_]);
-
- break;
- case Dumux::BoundaryConditions::dirichlet:
- bt_->boundaryInfo_[idx_].visited = 1;
- bt_->boundaryInfo_[idx_].isDirichlet = 1;
-
- Valgrind::SetDefined(bt_->boundaryInfo_[idx_]);
-
- bt_->eq2pvIdx_[idx_] = idx_;
- bt_->pv2eqIdx_[idx_] = idx_;
- break;
- default:
- DUNE_THROW(Dune::InvalidStateException,
- "Only neumann and dirichlet "
- "conditions are supported!");
- }
- return *this;
- }
-
- private:
- BoundaryTypes *bt_;
- int idx_;
- };
-
public:
BoundaryTypes()
{ reset(); }
@@ -88,6 +50,7 @@ public:
boundaryInfo_[i].visited = 0;
boundaryInfo_[i].isDirichlet = 0;
+ boundaryInfo_[i].isNeumann = 0;
eq2pvIdx_[i] = i;
pv2eqIdx_[i] = i;
@@ -121,6 +84,7 @@ public:
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) {
boundaryInfo_[eqIdx].visited = 1;
boundaryInfo_[eqIdx].isDirichlet = 0;
+ boundaryInfo_[eqIdx].isNeumann = 1;
Valgrind::SetDefined(boundaryInfo_[eqIdx]);
}
@@ -134,6 +98,7 @@ public:
for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
boundaryInfo_[eqIdx].visited = 1;
boundaryInfo_[eqIdx].isDirichlet = 1;
+ boundaryInfo_[eqIdx].isNeumann = 0;
eq2pvIdx_[eqIdx] = eqIdx;
pv2eqIdx_[eqIdx] = eqIdx;
@@ -167,6 +132,7 @@ public:
{
boundaryInfo_[eqIdx].visited = 1;
boundaryInfo_[eqIdx].isDirichlet = 1;
+ boundaryInfo_[eqIdx].isNeumann = 0;
// update the equation <-> primary variable mapping
eq2pvIdx_[eqIdx] = pvIdx;
@@ -211,7 +177,7 @@ public:
* neumann condition.
*/
bool isNeumann(unsigned eqIdx) const
- { return !boundaryInfo_[eqIdx].isDirichlet; };
+ { return boundaryInfo_[eqIdx].isNeumann; };
/*!
* \brief Returns true if some equation is used to specify a
@@ -220,7 +186,7 @@ public:
bool hasNeumann() const
{
for (int i = 0; i < numEq; ++i)
- if (!boundaryInfo_[i].isDirichlet)
+ if (boundaryInfo_[i].isNeumann)
return true;
return false;
};
@@ -241,54 +207,12 @@ public:
unsigned eqToDirichletIndex(unsigned eqIdx) const
{ return eq2pvIdx_[eqIdx]; };
- /*!
- * \brief Allows to assign Dumux::BoundaryCondition to a single
- * primary varible or equation.
- *
- * In the case of Dirichlet conditions, the equation with the same
- * index is always disabled.
- */
- BoundaryWrapper operator[](int i) DUNE_DEPRECATED
- { return BoundaryWrapper(this, i); }
-
- /*!
- * \brief Allows to assign Dumux::BoundaryCondition to a all
- * primary varibles or equations.
- *
- * To get rid of the "deprecated warning" do the following:
- * boundaryTypeVector is now a class, which has the methods
- * setAllNeumann() and setAllDirichlet().
- *
- * deprecated System -> current System
- *
- * values=BoundaryConditions::neumann; -> values.setAllNeumann();
- *
- */
- BoundaryTypes &operator=(int bc) DUNE_DEPRECATED
- {
- switch (bc) {
- case Dumux::BoundaryConditions::neumann:
- reset();
- setAllNeumann();
- break;
- case Dumux::BoundaryConditions::dirichlet:
- reset();
- setAllDirichlet();
- break;
-
- default:
- DUNE_THROW(Dune::InvalidStateException,
- "Only neumann and dirichlet "
- "conditions are supported!");
- }
- return *this;
- }
-
private:
// this is a bitfield structure!
struct __packed__ {
unsigned char visited : 1;
unsigned char isDirichlet : 1;
+ unsigned char isNeumann : 1;
} boundaryInfo_[numEq];
unsigned char eq2pvIdx_[numEq];
diff --git a/dumux/common/exceptions.hh b/dumux/common/exceptions.hh
index 40e929052a..f4dfa6a561 100644
--- a/dumux/common/exceptions.hh
+++ b/dumux/common/exceptions.hh
@@ -1,4 +1,4 @@
-// $Id: exceptions.hh 3356 2010-03-25 13:01:37Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
diff --git a/dumux/common/math.hh b/dumux/common/math.hh
index a9ab956f43..86d71e4c6b 100644
--- a/dumux/common/math.hh
+++ b/dumux/common/math.hh
@@ -1,4 +1,4 @@
-// $Id: math.hh 3736 2010-06-15 09:52:10Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/boxmodels/pdelab/preconditionerpdelab.hh b/dumux/common/pdelabpreconditioner.hh
similarity index 78%
rename from dumux/boxmodels/pdelab/preconditionerpdelab.hh
rename to dumux/common/pdelabpreconditioner.hh
index e049226361..c8124867f7 100644
--- a/dumux/boxmodels/pdelab/preconditionerpdelab.hh
+++ b/dumux/common/pdelabpreconditioner.hh
@@ -1,4 +1,4 @@
-// $Id: preconditionerpdelab.hh 3736 2010-06-15 09:52:10Z lauser $: preconditionerpdelab.hh 3728 2010-06-10 15:44:39Z bernd $
+// $Id$: preconditionerpdelab.hh 3728 2010-06-10 15:44:39Z bernd $
/*****************************************************************************
* Copyright (C) 2009-2010 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
@@ -13,41 +13,41 @@
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
-#ifndef DUMUX_PRECONDITIONERPDELAB_HH
-#define DUMUX_PRECONDITIONERPDELAB_HH
+#ifndef DUMUX_PDELAB_PRECONDITIONER_HH
+#define DUMUX_PDELAB_PRECONDITIONER_HH
#include<dune/pdelab/backend/istlsolverbackend.hh>
-namespace Dune {
+namespace Dumux {
namespace PDELab {
template<class TypeTag>
class Exchanger
{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
dim = GridView::dimension
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Grid)) Grid;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::VertexMapper VertexMapper;
- typedef typename GET_PROP(TypeTag, PTAG(PDELabTypes)) PDELabTypes;
- typedef typename PDELabTypes::GridOperatorSpace GridOperatorSpace;
- typedef typename GridOperatorSpace::template MatrixContainer<Scalar>::Type Matrix;
- typedef typename Matrix::block_type BlockType;
- typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
+
+ typedef typename JacobianMatrix::block_type BlockType;
+ typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename Grid::Traits::GlobalIdSet IDS;
typedef typename IDS::IdType IdType;
typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
- typedef typename RefElemProp::Container ReferenceElements;
- typedef typename RefElemProp::ReferenceElement ReferenceElement;
+ typedef typename RefElemProp::Container ReferenceElements;
+ typedef typename RefElemProp::ReferenceElement ReferenceElement;
public:
- Exchanger(const Model& model)
- : gridView_(model.gridView()), vertexMapper_(model.vertexMapper()), borderIndices_(0)
+ Exchanger(const Problem& problem)
+ : gridView_(problem.gridView()), vertexMapper_(problem.vertexMapper()), borderIndices_(0)
{
gid2Index_.clear();
index2GID_.clear();
@@ -59,7 +59,7 @@ public:
// << " at (" << vertexIt->geometry().corner(0) << ") is of type "
// << vertexIt->partitionType() << ", GID = "
// << gridView_.grid().globalIdSet().id(*vertexIt) << std::endl;
- if (vertexIt->partitionType() == BorderEntity)
+ if (vertexIt->partitionType() == Dune::BorderEntity)
{
int localIdx = vertexMapper_.map(*vertexIt);
IdType globalIdx = gridView_.grid().globalIdSet().id(*vertexIt);
@@ -82,15 +82,16 @@ public:
MatEntry (const IdType& f, const BlockType& s) : first(f),second(s) {}
MatEntry () {}
};
-
+
// A DataHandle class to exchange matrix entries
class MatEntryExchange
- : public CommDataHandleIF<MatEntryExchange,MatEntry> {
- typedef typename Matrix::RowIterator RowIterator;
- typedef typename Matrix::ColIterator ColIterator;
-public:
- //! export type of data for message buffer
- typedef MatEntry DataType;
+ : public Dune::CommDataHandleIF<MatEntryExchange,MatEntry>
+ {
+ typedef typename JacobianMatrix::RowIterator RowIterator;
+ typedef typename JacobianMatrix::ColIterator ColIterator;
+ public:
+ //! export type of data for message buffer
+ typedef MatEntry DataType;
//! returns true if data for this codim should be communicated
bool contains (int dim, int codim) const
@@ -184,7 +185,7 @@ public:
MatEntryExchange (const GridView& gridView, const std::map<IdType,int>& g2i,
const std::map<int,IdType>& i2g,
const VertexMapper& vm,
- Matrix& A)
+ JacobianMatrix& A)
: gridView_(gridView), gid2Index_(g2i), index2GID_(i2g), vertexMapper_(vm), A_(A)
{}
@@ -193,15 +194,17 @@ private:
const std::map<IdType,int>& gid2Index_;
const std::map<int,IdType>& index2GID_;
const VertexMapper& vertexMapper_;
- Matrix& A_;
+ JacobianMatrix& A_;
};
- void sumEntries (Matrix& A)
+ void sumEntries (JacobianMatrix& A)
{
if (gridView_.comm().size() > 1)
{
MatEntryExchange datahandle(gridView_, gid2Index_, index2GID_, vertexMapper_, A);
- gridView_.communicate(datahandle, InteriorBorder_InteriorBorder_Interface, ForwardCommunication);
+ gridView_.communicate(datahandle,
+ Dune::InteriorBorder_InteriorBorder_Interface,
+ Dune::ForwardCommunication);
}
}
@@ -237,7 +240,8 @@ public:
//! Constructor.
NonoverlappingWrappedPreconditioner (const GFS& gfs_, P& prec_, const CC& cc_,
- const std::vector<int>& borderIndices, const ParallelISTLHelper<GFS>& helper_)
+ const std::vector<int>& borderIndices,
+ const Dune::PDELab::ParallelISTLHelper<GFS>& helper_)
: gfs(gfs_), prec(prec_), cc(cc_), borderIndices_(borderIndices), helper(helper_)
{}
@@ -285,16 +289,16 @@ private:
P& prec;
const CC& cc;
const std::vector<int>& borderIndices_;
- const ParallelISTLHelper<GFS>& helper;
+ const Dune::PDELab::ParallelISTLHelper<GFS>& helper;
};
template<class TypeTag>
class ISTLBackend_NoOverlap_BCGS_ILU
{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP(TypeTag, PTAG(PDELabTypes)) PDELabTypes;
- typedef typename PDELabTypes::GridFunctionSpace GridFunctionSpace;
- typedef typename PDELabTypes::ConstraintsTrafo ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
typedef Dune::PDELab::ParallelISTLHelper<GridFunctionSpace> PHELPER;
public:
@@ -304,10 +308,13 @@ public:
\param[in] maxiter maximum number of iterations to do
\param[in] verbose print messages if true
*/
- explicit ISTLBackend_NoOverlap_BCGS_ILU (Model& model, unsigned maxiter_=5000, int verbose_=1)
- : gfs(model.jacobianAssembler().gridFunctionSpace()), phelper(gfs),
- maxiter(maxiter_), verbose(verbose_), constraintsTrafo_(model.jacobianAssembler().constraintsTrafo()),
- exchanger_(model)
+ explicit ISTLBackend_NoOverlap_BCGS_ILU (Problem& problem, unsigned maxiter_=5000, int verbose_=1)
+ : gfs(problem.model().jacobianAssembler().gridFunctionSpace()),
+ phelper(gfs),
+ maxiter(maxiter_),
+ verbose(verbose_),
+ constraintsTrafo_(problem.model().jacobianAssembler().constraintsTrafo()),
+ exchanger_(problem)
{}
/*! \brief compute global norm of a vector
@@ -331,19 +338,19 @@ public:
\param[in] r right hand side
\param[in] reduction to be achieved
*/
- template<class Matrix, class SolVector, class RhsVector>
- void apply(Matrix& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
+ template<class JacobianMatrix, class SolVector, class RhsVector>
+ void apply(JacobianMatrix& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
{
- typedef Dune::SeqILU0<Matrix,SolVector,RhsVector> SeqPreCond;
- Matrix B(A);
+ typedef Dune::SeqILU0<JacobianMatrix,SolVector,RhsVector> SeqPreCond;
+ JacobianMatrix B(A);
exchanger_.sumEntries(B);
SeqPreCond seqPreCond(B, 0.9);
- typedef Dune::PDELab::NonoverlappingOperator<GridFunctionSpace,Matrix,SolVector,RhsVector> POP;
+ typedef Dune::PDELab::NonoverlappingOperator<GridFunctionSpace,JacobianMatrix,SolVector,RhsVector> POP;
POP pop(gfs,A,phelper);
typedef Dune::PDELab::NonoverlappingScalarProduct<GridFunctionSpace,SolVector> PSP;
PSP psp(gfs,phelper);
- typedef Dune::PDELab::NonoverlappingWrappedPreconditioner<ConstraintsTrafo, GridFunctionSpace, SeqPreCond> ParPreCond;
+ typedef NonoverlappingWrappedPreconditioner<ConstraintsTrafo, GridFunctionSpace, SeqPreCond> ParPreCond;
ParPreCond parPreCond(gfs, seqPreCond, constraintsTrafo_, exchanger_.borderIndices(), phelper);
int verb=0;
@@ -376,10 +383,10 @@ private:
template<class TypeTag>
class ISTLBackend_NoOverlap_Loop_Pardiso
{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP(TypeTag, PTAG(PDELabTypes)) PDELabTypes;
- typedef typename PDELabTypes::GridFunctionSpace GridFunctionSpace;
- typedef typename PDELabTypes::ConstraintsTrafo ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
typedef Dune::PDELab::ParallelISTLHelper<GridFunctionSpace> PHELPER;
public:
@@ -389,10 +396,13 @@ public:
\param[in] maxiter maximum number of iterations to do
\param[in] verbose print messages if true
*/
- explicit ISTLBackend_NoOverlap_Loop_Pardiso (Model& model, unsigned maxiter_=5000, int verbose_=1)
- : gfs(model.jacobianAssembler().gridFunctionSpace()), phelper(gfs),
- maxiter(maxiter_), verbose(verbose_), constraintsTrafo_(model.jacobianAssembler().constraintsTrafo()),
- exchanger_(model)
+ explicit ISTLBackend_NoOverlap_Loop_Pardiso (Problem& problem, unsigned maxiter_=5000, int verbose_=1)
+ : gfs(problem.model().jacobianAssembler().gridFunctionSpace()),
+ phelper(gfs),
+ maxiter(maxiter_),
+ verbose(verbose_),
+ constraintsTrafo_(problem.model().jacobianAssembler().constraintsTrafo()),
+ exchanger_(problem)
{}
/*! \brief compute global norm of a vector
@@ -416,19 +426,19 @@ public:
\param[in] r right hand side
\param[in] reduction to be achieved
*/
- template<class Matrix, class SolVector, class RhsVector>
- void apply(Matrix& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
+ template<class JacobianMatrix, class SolVector, class RhsVector>
+ void apply(JacobianMatrix& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
{
- typedef Dune::SeqPardiso<Matrix,SolVector,RhsVector> SeqPreCond;
- Matrix B(A);
+ typedef Dune::SeqPardiso<JacobianMatrix,SolVector,RhsVector> SeqPreCond;
+ JacobianMatrix B(A);
exchanger_.sumEntries(B);
SeqPreCond seqPreCond(B);
- typedef Dune::PDELab::NonoverlappingOperator<GridFunctionSpace,Matrix,SolVector,RhsVector> POP;
+ typedef Dune::PDELab::NonoverlappingOperator<GridFunctionSpace,JacobianMatrix,SolVector,RhsVector> POP;
POP pop(gfs,A,phelper);
typedef Dune::PDELab::NonoverlappingScalarProduct<GridFunctionSpace,SolVector> PSP;
PSP psp(gfs,phelper);
- typedef Dune::PDELab::NonoverlappingWrappedPreconditioner<ConstraintsTrafo, GridFunctionSpace, SeqPreCond> ParPreCond;
+ typedef NonoverlappingWrappedPreconditioner<ConstraintsTrafo, GridFunctionSpace, SeqPreCond> ParPreCond;
ParPreCond parPreCond(gfs, seqPreCond, constraintsTrafo_, exchanger_.borderIndices(), phelper);
// typedef Dune::PDELab::NonoverlappingRichardson<GridFunctionSpace,SolVector,RhsVector> PRICH;
diff --git a/dumux/common/propertysystem.hh b/dumux/common/propertysystem.hh
index 3c7219b93f..34e2f73f57 100644
--- a/dumux/common/propertysystem.hh
+++ b/dumux/common/propertysystem.hh
@@ -1,4 +1,4 @@
-// $Id: propertysystem.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/common/spline.hh b/dumux/common/spline.hh
index 74724156c9..0d719d77b0 100644
--- a/dumux/common/spline.hh
+++ b/dumux/common/spline.hh
@@ -1,4 +1,4 @@
-// $Id: spline.hh 3736 2010-06-15 09:52:10Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -31,6 +31,7 @@
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
+#include <dune/istl/btdmatrix.hh>
namespace Dumux
{
@@ -99,7 +100,7 @@ public:
const Scalar d_i =
6 / (h_(i) + h_(i+1))
*
- ( (y_[i+1] - y[i])/h_(i+1) - (y_[i] - y[i-1])/h_(i));
+ ( (y_[i+1] - y_[i])/h_(i+1) - (y_[i] - y_[i-1])/h_(i));
M[i][i-1] = mu_i;
M[i][i] = 2;
@@ -144,6 +145,18 @@ public:
bool applies(Scalar x) const
{ return x_[0] <= x && x <= x_[numSamples - 1]; };
+ /*!
+ * \brief Return the x value of the leftmost sampling point.
+ */
+ Scalar xMin() const
+ { return x_[0]; };
+
+ /*!
+ * \brief Return the x value of the rightmost sampling point.
+ */
+ Scalar xMax() const
+ { return x_[numSamples - 1]; };
+
/*!
* \brief Evaluate the spline at a given position.
*/
@@ -402,6 +415,416 @@ private:
FieldVector y_; // y positions of the sampling points
};
+/*
+ * \brief A 3rd order polynomial spline.
+
+ * This class implements a spline s(x) for which, given n (known at
+ * run-time) sampling points x_1, ..., x_n, the following holds:
+ *
+ * s(x_i) = y_i
+ * s'(x_1) = m_1
+ * s'(x_n) = m_n
+ *
+ * for any given boundary slopes m_1 and m_n.
+ */
+template<class ScalarT>
+class Spline<ScalarT, -1>
+{
+ typedef ScalarT Scalar;
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar,1> > BlockVector;
+ typedef Dune::BTDMatrix<Dune::FieldMatrix<Scalar, 1, 1> > BTDMatrix;
+
+public:
+ Spline()
+ { };
+
+ /*!
+ * \brief Set the sampling points and the boundary slopes of the
+ * spline function.
+ */
+ void set(int numSamples,
+ const Scalar *x,
+ const Scalar *y,
+ Scalar m0,
+ Scalar m1)
+ {
+ BTDMatrix M(numSamples);
+ BlockVector d(numSamples);
+ fillNatural_(M, d, numSamples, x, y);
+
+ int n = numSamples - 1;
+ // first row
+ M[0][1] = 1;
+ d[0] = 6/h_(1) * ( (y_[1] - y_[0])/h_(1) - m0);
+
+ // last row
+ M[n][n - 1] = 1;
+ d[n] =
+ 6/h_(n)
+ *
+ (m1 - (y_[n] - y_[n - 1])/h_(n));
+
+ // solve for the moments
+ M.solve(moments_, d);
+ }
+
+ /*!
+ * \brief Set the sampling points of the
+ * spline function.
+ *
+ * The second derivatives at the boundaries are 0 in this
+ * case. (i.e. this is a natural spline)
+ */
+ void set(int numSamples,
+ const Scalar *x,
+ const Scalar *y)
+ {
+ BTDMatrix M(numSamples);
+ BlockVector d(numSamples);
+ fillNatural_(M, d, numSamples, x, y);
+
+ // solve for the moments
+ M.solve(moments_, d);
+ }
+
+ /*!
+ * \brief Return true iff the given x is in range [x1, xn].
+ */
+ bool applies(Scalar x) const
+ {
+ return x_[0] <= x && x <= x_[numSamples() - 1];
+ };
+
+ /*!
+ * \brief Return the x value of the leftmost sampling point.
+ */
+ Scalar xMin() const
+ { return x_[0]; };
+
+ /*!
+ * \brief Return the x value of the rightmost sampling point.
+ */
+ Scalar xMax() const
+ { return x_[numSamples() - 1]; };
+
+ /*!
+ * \brief Returns the number of sampling points.
+ */
+ int numSamples() const
+ { return x_.size(); }
+
+ /*!
+ * \brief Evaluate the spline at a given position.
+ */
+ Scalar eval(Scalar x) const
+ {
+ assert(applies(x));
+
+ // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
+ // Springer, 2005, p. 109
+
+ int i = findIntervalIdx_(x);
+
+ Scalar h_i1 = h_(i + 1);
+ Scalar x_i = x - x_[i];
+ Scalar x_i1 = x_[i+1] - x;
+
+ Scalar A_i =
+ (y_[i+1] - y_[i])/h_i1
+ -
+ h_i1/6*(moments_[i+1] - moments_[i]);
+ Scalar B_i = y_[i] - moments_[i]* (h_i1*h_i1) / 6;
+
+ return
+ moments_[i]* x_i1*x_i1*x_i1 / (6 * h_i1)
+ +
+ moments_[i + 1]* x_i*x_i*x_i / (6 * h_i1)
+ +
+ A_i*x_i
+ +
+ B_i;
+ }
+
+ /*!
+ * \brief Evaluate the spline's derivative at a given position.
+ */
+ Scalar evalDerivative(Scalar x) const
+ {
+ assert(applies(x));
+
+ // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
+ // Springer, 2005, p. 109
+
+ int i = findIntervalIdx_(x);
+
+ Scalar h_i1 = h_(i + 1);
+ Scalar x_i = x - x_[i];
+ Scalar x_i1 = x_[i+1] - x;
+
+ Scalar A_i =
+ (y_[i+1] - y_[i])/h_i1
+ -
+ h_i1/6*(moments_[i+1] - moments_[i]);
+
+ return
+ -moments_[i] * x_i1*x_i1 / (2 * h_i1)
+ +
+ moments_[i + 1] * x_i*x_i / (2 * h_i1)
+ +
+ A_i;
+ }
+
+ /*!
+ * \brief Returns 1 if the spline is monotonically increasing, -1
+ * if the spline is mononously decreasing and 0 if the
+ * spline is not monotonous in the interval (x0, x1).
+ *
+ * In the corner case where the whole spline is flat, it returns
+ * 2.
+ */
+ int monotonic(Scalar xi0, Scalar xi1) const
+ {
+ assert(applies(xi0));
+ assert(applies(xi1));
+
+ Scalar x0 = std::min(xi0, xi1);
+ Scalar x1 = std::max(xi0, xi1);
+
+ // corner case where the whole spline is a constant
+ if (moments_[0] == 0 &&
+ moments_[1] == 0 &&
+ y_[0] == y_[1])
+ {
+ // actually the is monotonically increasing as well as
+ // monotonously decreasing
+ return 2;
+ }
+
+
+ int i = findIntervalIdx_(x0);
+ if (x_[i] <= x0 && x1 <= x_[i+1]) {
+ // the interval to check is completely included in a
+ // single segment
+ return monotonic_(i, x0, x1);
+ }
+
+ // make sure that the segments which are completly in the
+ // interval [x0, x1] all exhibit the same monotonicity.
+ int iEnd = findIntervalIdx_(x1);
+ int r = monotonic_(i, x0, x_[1]);
+ for (; i < iEnd - 1; ++i)
+ if (r != monotonic_(i, x_[i], x_[i + 1]))
+ return 0;
+
+ // check for the last segment
+ if (x_[iEnd] < x1 && r != monotonic_(iEnd, x_[iEnd], x1))
+ { return 0; }
+
+ return r;
+ }
+
+ /*!
+ * \brief Same as monotonic(x0, x1), but with the entire range of the
+ * spline as interval.
+ */
+ int monotonic() const
+ { return monotonic(x_[0], x_[numSamples() - 1]); }
+
+ /*!
+ * \brief Prints k tuples of the format (x, y, dx/dy, isMonotonic)
+ * to stdout.
+ *
+ * If the spline does not apply for parts of [x0, x1] it is
+ * extrapolated using a straight line. The result can be inspected
+ * using the following commands:
+ *
+ ----------- snip -----------
+ ./yourProgramm > 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 "Monotonic"
+ ----------- snap -----------
+ */
+ void printCSV(Scalar xi0, Scalar xi1, int k) const
+ {
+ Scalar x0 = std::min(xi0, xi1);
+ Scalar x1 = std::max(xi0, xi1);
+ const int n = numSamples() - 1;
+ for (int i = 0; i <= k; ++i) {
+ double x = i*(x1 - x0)/k + x0;
+ double x_p1 = x + (x1 - x0)/k;
+ double y;
+ double dy_dx;
+ double mono = 1;
+ if (!applies(x)) {
+ if (x < 0) {
+ dy_dx = evalDerivative(x_[0]);
+ y = (x - x_[0])*dy_dx + y_[0];
+ mono = (dy_dx>0)?1:-1;
+ }
+ else if (x > x_[n]) {
+ dy_dx = evalDerivative(x_[n]);
+ y = (x - x_[n])*dy_dx + y_[n];
+ mono = (dy_dx>0)?1:-1;
+ }
+ else {
+ std::cerr << "ooops: " << x << "\n";
+ exit(1);
+ }
+ }
+ else {
+ y = eval(x);
+ dy_dx = evalDerivative(x);
+ mono = monotonic(std::max<Scalar>(x_[0][0], x), std::min<Scalar>(x_[n][0], x_p1));
+ }
+
+ std::cout << x << " " << y << " " << dy_dx << " " << mono << "\n";
+ }
+ }
+
+private:
+ // fill the system of equations for a natrural spline
+ void fillNatural_(BTDMatrix &M,
+ BlockVector &d,
+ int numSamples,
+ const Scalar *x,
+ const Scalar *y)
+ {
+ x_.resize(numSamples);
+ y_.resize(numSamples);
+ moments_.resize(numSamples);
+
+ // copy sample points, make sure that the first x value is
+ // smaller than the last one
+ for (int i = 0; i < numSamples; ++i) {
+ int idx = 0;
+ if (x[0] > x[numSamples - 1])
+ idx = numSamples - i - 1;
+ x_[i] = x[idx];
+ y_[i] = y[idx];
+ }
+
+ M = 0;
+ d = 0;
+
+ // See: J. Stoer: "Numerische Mathematik 1", 9th edition,
+ // Springer, 2005, p. 111
+ const int n = numSamples - 1;
+
+ // second to next to last rows
+ for (int i = 1; i < n; ++i) {
+ const Scalar lambda_i = h_(i + 1) / (h_(i) + h_(i+1));
+ const Scalar mu_i = 1 - lambda_i;
+ const Scalar d_i =
+ 6 / (h_(i) + h_(i+1))
+ *
+ ( (y_[i+1] - y_[i])/h_(i+1) - (y_[i] - y_[i-1])/h_(i));
+
+ M[i][i-1] = mu_i;
+ M[i][i] = 2;
+ M[i][i + 1] = lambda_i;
+ d[i] = d_i;
+ };
+
+ // first row
+ M[0][0] = 2;
+
+ // last row
+ M[n][n] = 2;
+ }
+
+ // returns the monotonicality of an interval of a spline segment
+ int monotonic_(int i, Scalar x0, Scalar x1) const
+ {
+ // shift the interval so that it is consistent with the
+ // definitions by Stoer
+ x0 = x0 - x_[i];
+ x1 = x1 - x_[i];
+
+ Scalar a = a_(i);
+ Scalar b = b_(i);
+ Scalar c = c_(i);
+
+ Scalar disc = 4*b*b - 12*a*c;
+ if (disc < 0) {
+ // discriminant is smaller than 0, i.e. the segment does
+ // not exhibit any extrema.
+ return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
+ }
+ disc = std::sqrt(disc);
+ Scalar xE1 = (-2*b + disc)/(6*a);
+ Scalar xE2 = (-2*b - disc)/(6*a);
+
+ if (disc == 0) {
+ // saddle point -> no extrema
+ if (xE1 == x0)
+ // make sure that we're not picking the saddle point
+ // to determine whether we're monotonically increasing
+ // or decreasing
+ x0 = x1;
+ return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
+ };
+ if ((x0 < xE1 && xE1 < x1) ||
+ (x0 < xE2 && xE2 < x1))
+ {
+ // there is an extremum in the range (x0, x1)
+ return 0;
+ }
+ // no extremum in range (x0, x1)
+ x0 = (x0 + x1)/2; // pick point in the middle of the interval
+ // to avoid extrema on the boundaries
+ return (x0*(x0*3*a + 2*b) + c > 0) ? 1 : -1;
+ };
+
+ Scalar h_(int i) const
+ { return x_[i] - x_[i-1]; }
+
+ int findIntervalIdx_(Scalar x) const
+ {
+ // bisection
+ int iLow = 0;
+ int iHigh = numSamples() - 1;
+
+ while (iLow + 1 < iHigh) {
+ int i = (iLow + iHigh) / 2;
+ if (x_[i] > x)
+ iHigh = i;
+ else
+ iLow = i;
+ };
+ return iLow;
+ };
+
+ // returns the coefficient in front of the x^3 term. In Stoer this
+ // is delta.
+ Scalar a_(int i) const
+ { return (moments_[i+1] - moments_[i])/(6*h_(i+1)); }
+
+ // returns the coefficient in front of the x^2 term In Stoer this
+ // is gamma.
+ Scalar b_(int i) const
+ { return moments_[i]/2; }
+
+ // returns the coefficient in front of the x term. In Stoer this
+ // is beta.
+ Scalar c_(int i) const
+ { return (y_[i+1] - y_[i])/h_(i + 1) - h_(i+1)/6*(2*moments_[i] + moments_[i+1]); }
+
+ // returns the coefficient in front of the constant term. In Stoer this
+ // is alpha.
+ Scalar d_(int i) const
+ { return y_[i]; }
+
+ const BlockVector &toBlockVector_(const Scalar *v) const
+ { return *reinterpret_cast<const BlockVector*>(v); };
+
+ BlockVector moments_; // moments
+ BlockVector x_; // x positions of the sampling points
+ BlockVector y_; // y positions of the sampling points
+};
+
/*
* \brief A 3rd order polynomial p(x) = a x^3 + b x^2 + c x + d for
* which, given two distinct points x1 and x2, the following holds:
@@ -520,6 +943,19 @@ public:
bool applies(Scalar x) const
{ return x1_ <= x && x <= x2_; };
+
+ /*!
+ * \brief Return the x value of the leftmost sampling point.
+ */
+ Scalar xMin() const
+ { return x1_; };
+
+ /*!
+ * \brief Return the x value of the rightmost sampling point.
+ */
+ Scalar xMax() const
+ { return x2_; };
+
/*!
* \brief Evaluate the polynomial at a given position.
*/
diff --git a/dumux/common/start.hh b/dumux/common/start.hh
index 0ef34fbfb4..a64a69a810 100644
--- a/dumux/common/start.hh
+++ b/dumux/common/start.hh
@@ -1,4 +1,4 @@
-// $Id: start.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -52,9 +52,11 @@ int startFromDGF(int argc, char **argv)
#ifdef NDEBUG
try {
#endif
+
typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(Grid)) Grid;
typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(ProblemTypeTag, PTAG(TimeManager)) TimeManager;
typedef Dune::GridPtr<Grid> GridPointer;
// initialize MPI, finalize is done automatically on exit
@@ -65,24 +67,24 @@ int startFromDGF(int argc, char **argv)
printUsage(argv[0]);
// deal with the restart stuff
- int argPos = 1;
+ int argIdx = 1;
bool restart = false;
double restartTime = 0;
- if (std::string("--restart") == argv[argPos]) {
+ if (std::string("--restart") == argv[argIdx]) {
restart = true;
- ++argPos;
+ ++argIdx;
- std::istringstream(argv[argPos++]) >> restartTime;
+ std::istringstream(argv[argIdx++]) >> restartTime;
}
- if (argc - argPos != 3) {
+ if (argc - argIdx != 3) {
printUsage(argv[0]);
}
double tEnd, dt;
- const char *dgfFileName = argv[argPos++];
- std::istringstream(argv[argPos++]) >> tEnd;
- std::istringstream(argv[argPos++]) >> dt;
+ const char *dgfFileName = argv[argIdx++];
+ std::istringstream(argv[argIdx++]) >> tEnd;
+ std::istringstream(argv[argIdx++]) >> dt;
// create grid
// -> load the grid from file
@@ -90,15 +92,14 @@ int startFromDGF(int argc, char **argv)
Dune::gridinfo(*gridPtr);
// instantiate and run the concrete problem
- Problem problem(gridPtr->leafView());
-
+ TimeManager timeManager;
+ Problem problem(timeManager, gridPtr->leafView());
+ timeManager.init(problem, 0, dt, tEnd, !restart);
if (restart)
- problem.deserialize(restartTime);
-
- if (!problem.simulate(dt, tEnd))
- return 2;
-
+ problem.restart(restartTime);
+ timeManager.run();
return 0;
+
#ifdef NDEBUG
}
catch (Dune::Exception &e) {
diff --git a/dumux/common/timemanager.hh b/dumux/common/timemanager.hh
index 1de0366569..5f06653f75 100644
--- a/dumux/common/timemanager.hh
+++ b/dumux/common/timemanager.hh
@@ -1,4 +1,4 @@
-// $Id: timemanager.hh 3736 2010-06-15 09:52:10Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -25,11 +25,14 @@
#include <dune/common/timer.hh>
#include <dune/common/mpihelper.hh>
+#include <dumux/io/restart.hh>
+
namespace Dumux
{
/*!
* \addtogroup SimControl Simulation Supervision
*/
+
/*!
* \ingroup SimControl
* \brief Simplify the handling of time dependent problems.
@@ -50,16 +53,15 @@ namespace Dumux
* \todo Change the time manager to the property system (?)
* \todo Remove the episode identifier stuff
*/
-template <class EpisodeIdentiferT>
+template <class TypeTag>
class TimeManager
{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
public:
- typedef EpisodeIdentiferT EpisodeIdentifer;
- TimeManager(double endTime = 1e100,
- bool verbose = true)
+ TimeManager(bool verbose = true)
{
- wasRestarted_ = false;
verbose_ =
verbose &&
Dune::MPIHelper::getCollectiveCommunication().rank() == 0;
@@ -68,59 +70,81 @@ public:
episodeStartTime_ = 0;
time_ = 0.0;
- endTime_ = endTime;
- wasRestarted_ = false;
+ endTime_ = -1e100;
- stepSize_ = 1.0;
- stepNum_ = 0;
+ timeStepSize_ = 1.0;
+ timeStepIdx_ = 0;
finished_ = false;
- wasRestarted_ = false;
-
episodeLength_ = 1e100;
+
+ if (verbose_)
+ std::cout <<
+ "Welcome aboard DuMuX airlines. Please fasten your seatbelts! Emergency exits are near the time integration.\n";
}
/*!
- * \name Simulated time and time step management
- * @{
+ * \brief Initialize the model and problem and write the initial
+ * condition to disk.
*/
+ void init(Problem &problem,
+ Scalar tStart,
+ Scalar dtInitial,
+ Scalar tEnd,
+ bool writeInitialSol = true)
+ {
+ std::cout << "Initializing '" << problem_->name() << "'\n";
+
+ problem_ = &problem;
+ time_ = tStart;
+ timeStepSize_ = dtInitial;
+ endTime_ = tEnd;
+
+ // initialize the problem
+ problem_->init();
+
+ // initialize the problem
+ if (writeInitialSol) {
+ time_ -= timeStepSize_;
+ problem_->writeOutput();
+ time_ += timeStepSize_;
+ }
+ }
/*!
- * \brief Set the current simulated time, don't change the
- * current time step number.
+ * \name Simulated time and time step management
+ * @{
*/
- void setTime(double t)
- { time_ = t; }
/*!
- * \brief Set the current simulated time and the time step
- * number.
+ * \brief Set the current simulated time, don't change the current
+ * time step index.
*/
- void setTime(double t, int stepNum)
- { time_ = t; stepNum_ = stepNum; }
+ void setTime(Scalar t)
+ { time_ = t; }
/*!
- * \brief Let one time step size pass.
+ * \brief Set the current simulated time and the time step index.
*/
- void proceed()
- { time_ += stepSize_; ++stepNum_; }
+ void setTime(Scalar t, int stepIdx)
+ { time_ = t; timeStepIdx_ = stepIdx; }
/*!
* \brief Return the current simulated time.
*/
- double time() const
+ Scalar time() const
{ return time_; }
/*!
* \brief Returns the number of (simulated) seconds which the simulation runs.
*/
- double endTime() const
+ Scalar endTime() const
{ return endTime_; }
/*!
* \brief Set the time of simulated seconds at which the simulation runs.
*/
- void setEndTime(double val)
+ void setEndTime(Scalar val)
{ endTime_ = val; }
/*!
@@ -130,26 +154,22 @@ public:
* episode, the timeStep() method will take care that the
* step size won't exceed the episode, though.
*/
- void setTimeStepSize(double stepSize)
- {
- stepSize_ = stepSize;
- }
+ void setTimeStepSize(Scalar stepSize)
+ { timeStepSize_ = stepSize; }
/*!
* \brief Returns a suggested timestep length so that we don't
* miss the beginning of the next episode.
*/
- double timeStepSize() const
- {
- return stepSize_;
- }
+ Scalar timeStepSize() const
+ { return timeStepSize_; }
/*!
* \brief Returns number of time steps which have been
* executed since t=0.
*/
- int timeStepNum() const
- { return stepNum_; }
+ int timeStepIndex() const
+ { return timeStepIdx_; }
/*!
* \brief Specify whether the simulation is finished
@@ -180,30 +200,12 @@ public:
* \param tStart Time when the episode began
* \param len Length of the episode
*/
- void startNextEpisode(const EpisodeIdentifer &id,
- double tStart,
- double len)
+ void startNextEpisode(Scalar tStart,
+ Scalar len)
{
++ episodeIndex_;
episodeStartTime_ = tStart;
episodeLength_ = len;
-
- episode_ = id;
- }
-
- /*!
- * \brief Change the current episode of the simulation
- * assuming that the episode starts at the current
- * time.
- */
- void startNextEpisode(const EpisodeIdentifer &id,
- double len = 1e100)
- {
- ++ episodeIndex_;
- episodeStartTime_ = time_;
- episodeLength_ = len;
-
- episode_ = id;
}
/*!
@@ -212,25 +214,13 @@ public:
*
* \param len Length of the episode, infinite if not specified.
*/
- void startNextEpisode(double len = 1e100)
+ void startNextEpisode(Scalar len = 1e100)
{
++ episodeIndex_;
episodeStartTime_ = time_;
episodeLength_ = len;
}
- /*!
- * \brief Returns the identifier of the current episode
- */
- const EpisodeIdentifer &episode() const
- { return episode_; }
-
- /*!
- * \brief Returns the identifier of the current episode
- */
- EpisodeIdentifer &episode()
- { return episode_; }
-
/*!
* \brief Returns the index of the current episode.
*
@@ -243,14 +233,14 @@ public:
* \brief Returns the absolute time when the current episode
* started.
*/
- double episodeStartTime() const
+ Scalar episodeStartTime() const
{ return episodeStartTime_; }
/*!
* \brief Returns the length of the current episode in
* simulated time.
*/
- double episodeLength() const
+ Scalar episodeLength() const
{ return std::min(episodeLength_, endTime_ - episodeStartTime_); }
/*!
@@ -263,7 +253,7 @@ public:
/*!
* \brief Aligns dt to the episode boundary if t+dt exceeds the current episode.
*/
- double episodeMaxTimeStepSize() const
+ Scalar episodeMaxTimeStepSize() const
{
// if the current episode is over and the simulation
// wants to give it some extra time, we will return
@@ -289,49 +279,46 @@ public:
* This method makes sure that time steps sizes are aligned to
* episode boundaries, amongst other stuff.
*/
- template <class Problem>
- void runSimulation(Problem &problem)
+ void run()
{
- if (verbose_)
- std::cout <<
- "Welcome aboard DuMuX airlines. Please fasten your seatbelts! Emergency exits are near the time integration.\n";
-
Dune::Timer timer;
timer.reset();
- // initialize them model and write the initial
- // condition to disk
- double dtInitial = stepSize_;
- stepSize_ = 0.0;
- problem.init();
- stepSize_ = dtInitial;
-
// do the time steps
while (!finished())
{
- problem.timeStepBegin();
-
- // execute the time integration (i.e. Runge-Kutta
- // or Euler).
- problem.timeIntegration();
+ // pre-process the current solution
+ problem_->preProcess();
+
+ // execute the time integration scheme
+ problem_->timeIntegration();
+
+ // post-process the current solution
+ problem_->postProcess();
+
+ // write the result to disk
+ if (problem_->doOutput())
+ problem_->writeOutput();
// advance the simulated time by the current time step
// size
- double dt = timeStepSize();
- proceed();
-
- problem.timeStepEnd();
+ Scalar dt = timeStepSize();
+ time_ += timeStepSize_;
+ ++timeStepIdx_;
+ if (problem_->doSerialize())
+ problem_->serialize();
+
// notify the problem that the timestep is done and ask it
// for a suggestion for the next timestep size
- double nextDt =
- std::min(problem.nextTimeStepSize(),
+ Scalar nextDt =
+ std::min(problem_->nextTimeStepSize(),
episodeMaxTimeStepSize());
if (verbose_) {
std::cout <<
- boost::format("Timestep %d done. CPUt=%.4g, t=%.4g, StepSize=%.4g, NextStepSize=%.4g\n")
- %timeStepNum()%timer.elapsed()%time()%dt%nextDt;
+ boost::format("Time step %d done. CPU time:%.4g, time:%.4g, last step size:%.4g, next step size:%.4g\n")
+ %timeStepIndex()%timer.elapsed()%time()%dt%nextDt;
}
// set the time step size for the next step
@@ -339,9 +326,10 @@ public:
}
if (verbose_)
- std::cout <<
- boost::format("Simulation took %.3f seconds. Hopefully you enjoyed simulating with us. We hope that you also simulate with us next time.\n")
- %timer.elapsed();
+ std::cout << "Simulation took " << timer.elapsed() <<" seconds.\n"
+ << "We hope that you enjoyed simulating with us\n"
+ << "and that you will chose us next time, too.\n";
+
}
/*!
@@ -354,11 +342,12 @@ public:
template <class Restarter>
void serialize(Restarter &res)
{
- res.serializeSection("TimeManager");
+ res.serializeSectionBegin("TimeManager");
res.serializeStream() << episodeIndex_ << " "
<< episodeStartTime_ << " "
<< time_ << " "
- << stepNum_ << "\n";
+ << timeStepIdx_ << " ";
+ res.serializeSectionEnd();
}
/*!
@@ -367,16 +356,12 @@ public:
template <class Restarter>
void deserialize(Restarter &res)
{
- res.deserializeSection("TimeManager");
+ res.deserializeSectionBegin("TimeManager");
res.deserializeStream() >> episodeIndex_
>> episodeStartTime_
>> time_
- >> stepNum_;
-
-
- std::string dummy;
- std::getline(res.deserializeStream(), dummy);
- wasRestarted_ = true;
+ >> timeStepIdx_;
+ res.deserializeSectionEnd();
}
/*
@@ -384,20 +369,18 @@ public:
*/
private:
-
- int episodeIndex_;
- double episodeStartTime_;
- double episodeLength_;
- EpisodeIdentifer episode_;
-
- double time_;
- double endTime_;
-
- double stepSize_;
- int stepNum_;
- bool finished_;
- bool verbose_;
- bool wasRestarted_;
+ Problem *problem_;
+ int episodeIndex_;
+ Scalar episodeStartTime_;
+ Scalar episodeLength_;
+
+ Scalar time_;
+ Scalar endTime_;
+
+ Scalar timeStepSize_;
+ int timeStepIdx_;
+ bool finished_;
+ bool verbose_;
};
}
diff --git a/dumux/common/valgrind.hh b/dumux/common/valgrind.hh
index e77a2823e5..818f79a0f5 100644
--- a/dumux/common/valgrind.hh
+++ b/dumux/common/valgrind.hh
@@ -1,4 +1,4 @@
-// $Id: valgrind.hh 3732 2010-06-11 13:27:20Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/io/restart.hh b/dumux/io/restart.hh
index 4f9ae534ac..cf3c4a71cd 100644
--- a/dumux/io/restart.hh
+++ b/dumux/io/restart.hh
@@ -1,4 +1,4 @@
-// $Id: restart.hh 3736 2010-06-15 09:52:10Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -32,16 +32,14 @@
#include <fstream>
#include <sstream>
-
namespace Dumux {
-
/*!
* \brief Load or save a state of a model to/from the harddisk.
*/
-template <class GridView>
class Restart {
//! \brief Create a magic cookie for restart files, so that it is
//! unlikely to load a restart file for an incorrectly.
+ template <class GridView>
static const std::string magicRestartCookie_(const GridView &gridView)
{
const std::string gridName = "blubb"; //gridView.grid().name();
@@ -70,6 +68,7 @@ class Restart {
}
//! \brief Return the restart file name.
+ template <class GridView>
static const std::string restartFileName_(const GridView &gridView,
const std::string &simName,
double t)
@@ -92,18 +91,20 @@ public:
/*!
* \brief Write the current state of the model to disk.
*/
- void serializeBegin(const GridView &gridView,
- const std::string &simName,
- double t)
+ template <class Problem>
+ void serializeBegin(Problem &problem)
{
- const std::string magicCookie = magicRestartCookie_(gridView);
- fileName_ = restartFileName_(gridView, simName, t);
+ const std::string magicCookie = magicRestartCookie_(problem.gridView());
+ fileName_ = restartFileName_(problem.gridView(),
+ problem.name(),
+ problem.timeManager().time());
// open output file and write magic cookie
outStream_.open(fileName_.c_str());
outStream_.precision(20);
- serializeSection(magicCookie);
+ serializeSectionBegin(magicCookie);
+ serializeSectionEnd();
}
/*!
@@ -117,22 +118,28 @@ public:
/*!
* \brief Start a new section in the serialized output.
*/
- void serializeSection(const std::string &cookie)
+ void serializeSectionBegin(const std::string &cookie)
{
outStream_ << cookie << "\n";
}
+ /*!
+ * \brief End of a section in the serialized output.
+ */
+ void serializeSectionEnd()
+ { outStream_ << "\n"; }
+
/*!
* \brief Serialize all leaf entities of a codim in a gridView.
*
* The actual work is done by Serializer::serialize(Entity)
*/
- template <int codim, class Serializer>
+ template <int codim, class Serializer, class GridView>
void serializeEntities(Serializer &serializer,
const GridView &gridView)
{
std::string cookie = (boost::format("Entities: Codim %d")%codim).str();
- serializeSection(cookie);
+ serializeSectionBegin(cookie);
// write element data
typedef typename GridView::template Codim<codim>::Iterator Iterator;
@@ -143,6 +150,8 @@ public:
serializer.serializeEntity(outStream_, *it);
outStream_ << "\n";
};
+
+ serializeSectionEnd();
}
/*!
@@ -158,11 +167,12 @@ public:
* \brief Start reading a restart file at a certain simulated
* time.
*/
- void deserializeBegin(const GridView &grid,
- const std::string &simName,
- double t)
+ template <class Problem>
+ void deserializeBegin(Problem &problem, double t)
{
- fileName_ = restartFileName_(grid, simName, t);
+ fileName_ = restartFileName_(problem.gridView(),
+ problem.name(),
+ t);
// open input file and read magic cookie
inStream_.open(fileName_.c_str());
@@ -184,8 +194,11 @@ public:
}
inStream_.seekg(0, std::ios::beg);
- const std::string magicCookie = magicRestartCookie_(grid);
- deserializeSection(magicCookie);
+ const std::string magicCookie =
+ magicRestartCookie_(problem.gridView());
+
+ deserializeSectionBegin(magicCookie);
+ deserializeSectionEnd();
}
/*!
@@ -200,7 +213,7 @@ public:
/*!
* \brief Start reading a new section of the restart file.
*/
- void deserializeSection(const std::string &cookie)
+ void deserializeSectionBegin(const std::string &cookie)
{
if (!inStream_.good())
DUNE_THROW(Dune::IOError,
@@ -212,18 +225,33 @@ public:
"Could not start section '" << cookie << "'");
};
+ /*!
+ * \brief End of a section in the serialized output.
+ */
+ void deserializeSectionEnd()
+ {
+ std::string dummy;
+ std::getline(inStream_, dummy);
+ for (int i = 0; i < dummy.length(); ++i) {
+ if (!std::isspace(dummy[i])) {
+ DUNE_THROW(Dune::InvalidStateException,
+ "Encountered unread values while deserializing");
+ };
+ }
+ }
+
/*!
* \brief Deserialize all leaf entities of a codim in a grid.
*
* The actual work is done by Deserializer::deserialize(Entity)
*/
- template <int codim, class Deserializer>
+ template <int codim, class Deserializer, class GridView>
void deserializeEntities(Deserializer &deserializer,
const GridView &gridView)
{
std::string cookie = (boost::format("Entities: Codim %d")%codim).str();
- deserializeSection(cookie);
+ deserializeSectionBegin(cookie);
std::string curLine;
@@ -241,6 +269,8 @@ public:
std::istringstream curLineStream(curLine);
deserializer.deserializeEntity(curLineStream, *it);
};
+
+ deserializeSectionEnd();
}
/*!
@@ -263,7 +293,7 @@ public:
private:
- std::string fileName_;
+ std::string fileName_;
std::ifstream inStream_;
std::ofstream outStream_;
};
diff --git a/dumux/io/vtkmultiwriter.hh b/dumux/io/vtkmultiwriter.hh
index d73b8fea0e..6f939970fd 100644
--- a/dumux/io/vtkmultiwriter.hh
+++ b/dumux/io/vtkmultiwriter.hh
@@ -1,4 +1,4 @@
-// $Id: vtkmultiwriter.hh 3736 2010-06-15 09:52:10Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser *
@@ -287,8 +287,8 @@ public:
template <class Restarter>
void serialize(Restarter &res)
{
- res.serializeSection("VTKMultiWriter");
- res.serializeStream() << writerNum_ - 1 << "\n";
+ res.serializeSectionBegin("VTKMultiWriter");
+ res.serializeStream() << writerNum_ << "\n";
if (commRank_ == 0) {
size_t fileLen = 0;
@@ -311,6 +311,8 @@ public:
delete[] tmp;
}
}
+
+ res.serializeSectionEnd();
}
/*!
@@ -321,7 +323,7 @@ public:
{
wasRestarted_ = true;
- res.deserializeSection("VTKMultiWriter");
+ res.deserializeSectionBegin("VTKMultiWriter");
res.deserializeStream() >> writerNum_;
std::string dummy;
@@ -346,6 +348,8 @@ public:
multiFile_.seekp(filePos);
}
+
+ res.deserializeSectionEnd();
}
diff --git a/dumux/material/binarycoefficients/fullermethod.hh b/dumux/material/binarycoefficients/fullermethod.hh
index b7895d135b..decea899c9 100644
--- a/dumux/material/binarycoefficients/fullermethod.hh
+++ b/dumux/material/binarycoefficients/fullermethod.hh
@@ -1,4 +1,4 @@
-// $Id: fullermethod.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/binarycoefficients/h2o_n2.hh b/dumux/material/binarycoefficients/h2o_n2.hh
index aa152aee9a..682abd2a98 100644
--- a/dumux/material/binarycoefficients/h2o_n2.hh
+++ b/dumux/material/binarycoefficients/h2o_n2.hh
@@ -1,4 +1,4 @@
-// $Id: h2o_n2.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/binarycoefficients/henryiapws.hh b/dumux/material/binarycoefficients/henryiapws.hh
index f0baa26e8a..df4c30afdc 100644
--- a/dumux/material/binarycoefficients/henryiapws.hh
+++ b/dumux/material/binarycoefficients/henryiapws.hh
@@ -1,4 +1,4 @@
-// $Id: henryiapws.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/component.hh b/dumux/material/components/component.hh
index 466d3737ee..9ae42e9edc 100644
--- a/dumux/material/components/component.hh
+++ b/dumux/material/components/component.hh
@@ -1,4 +1,4 @@
-// $Id: component.hh 3824 2010-07-13 13:30:02Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
@@ -21,8 +21,6 @@
#ifndef DUMUX_COMPONENT_HH
#define DUMUX_COMPONENT_HH
-#include <dune/common/exceptions.hh>
-
namespace Dumux
{
@@ -33,6 +31,10 @@ template <class Scalar, class Implementation>
class Component
{
public:
+ static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
+ Scalar pressMin, Scalar pressMax, unsigned nPress)
+ { Dune::dwarn << "No init routine defined - make shure that this is not necessary!" << std::endl; }
+
/*!
* \brief A human readable name for the compoent.
*/
diff --git a/dumux/material/components/h2.hh b/dumux/material/components/h2.hh
index c736b782e1..31d707bb74 100644
--- a/dumux/material/components/h2.hh
+++ b/dumux/material/components/h2.hh
@@ -106,7 +106,7 @@ public:
/*!
* \brief The density [kg/m^3] of H2 at a given pressure and temperature.
*/
- static Scalar density(Scalar temperature, Scalar pressure)
+ static Scalar gasDensity(Scalar temperature, Scalar pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), temperature, pressure);
diff --git a/dumux/material/components/h2o.hh b/dumux/material/components/h2o.hh
index a197c7285a..53dd670ad5 100644
--- a/dumux/material/components/h2o.hh
+++ b/dumux/material/components/h2o.hh
@@ -1,4 +1,4 @@
-// $Id: h2o.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/iapws/common.hh b/dumux/material/components/iapws/common.hh
index 659b7a59b1..991cca533e 100644
--- a/dumux/material/components/iapws/common.hh
+++ b/dumux/material/components/iapws/common.hh
@@ -1,4 +1,4 @@
-// $Id: common.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009-2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/iapws/region1.hh b/dumux/material/components/iapws/region1.hh
index 528db976f2..cbcb8be2fe 100644
--- a/dumux/material/components/iapws/region1.hh
+++ b/dumux/material/components/iapws/region1.hh
@@ -1,4 +1,4 @@
-// $Id: region1.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009-2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/iapws/region2.hh b/dumux/material/components/iapws/region2.hh
index c8fb50d44f..f045998814 100644
--- a/dumux/material/components/iapws/region2.hh
+++ b/dumux/material/components/iapws/region2.hh
@@ -1,4 +1,4 @@
-// $Id: region2.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009-2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/iapws/region4.hh b/dumux/material/components/iapws/region4.hh
index e8c5066c5c..c527846b33 100644
--- a/dumux/material/components/iapws/region4.hh
+++ b/dumux/material/components/iapws/region4.hh
@@ -1,4 +1,4 @@
-// $Id: region4.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009-2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/n2.hh b/dumux/material/components/n2.hh
index 1b9ad9662c..c40943810d 100644
--- a/dumux/material/components/n2.hh
+++ b/dumux/material/components/n2.hh
@@ -1,4 +1,4 @@
-// $Id: n2.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/o2.hh b/dumux/material/components/o2.hh
index 7fa102bf5a..6cef5df49a 100644
--- a/dumux/material/components/o2.hh
+++ b/dumux/material/components/o2.hh
@@ -116,7 +116,7 @@ public:
*
* \todo density liquid oxygen
*/
- static Scalar density(Scalar temperature, Scalar pressure)
+ static Scalar gasDensity(Scalar temperature, Scalar pressure)
{
// Assume an ideal gas
return IdealGas::density(molarMass(), temperature, pressure);
diff --git a/dumux/material/components/simplednapl.hh b/dumux/material/components/simplednapl.hh
index a00a7d34f5..2ec18d3b36 100644
--- a/dumux/material/components/simplednapl.hh
+++ b/dumux/material/components/simplednapl.hh
@@ -1,4 +1,4 @@
-// $Id: simplednapl.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/simpleh2o.hh b/dumux/material/components/simpleh2o.hh
index ad8fd2aa64..601b209932 100644
--- a/dumux/material/components/simpleh2o.hh
+++ b/dumux/material/components/simpleh2o.hh
@@ -1,4 +1,4 @@
-// $Id: simpleh2o.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/tabulatedcomponent.hh b/dumux/material/components/tabulatedcomponent.hh
index 588cf6d14b..7cf871f744 100644
--- a/dumux/material/components/tabulatedcomponent.hh
+++ b/dumux/material/components/tabulatedcomponent.hh
@@ -1,4 +1,4 @@
-// $Id: tabulatedcomponent.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/components/unit.hh b/dumux/material/components/unit.hh
index 31b3d4afb4..333b290fd6 100644
--- a/dumux/material/components/unit.hh
+++ b/dumux/material/components/unit.hh
@@ -1,4 +1,4 @@
-// $Id: unit.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
index 5cd1d5d53d..2bf6fa2f0d 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscorey.hh
@@ -1,4 +1,4 @@
-// $Id: brookscorey.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh b/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh
index 3bf6a58395..ae4bf6cdde 100644
--- a/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh
@@ -1,4 +1,4 @@
-// $Id: brookscoreyparams.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
index 96f4c92eb0..c8369dbed6 100644
--- a/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh
@@ -1,4 +1,4 @@
-// $Id: efftoabslaw.hh 3811 2010-07-05 12:59:17Z pnuske $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh b/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh
index a1a7049e62..4ef0655952 100644
--- a/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/efftoabslawparams.hh
@@ -1,4 +1,4 @@
-// $Id: efftoabslawparams.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
index f8b6dd43aa..7e96ef2b33 100644
--- a/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh
@@ -1,4 +1,4 @@
-// $Id: linearmaterial.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh b/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh
index f3a9fe506f..b0c95985a5 100644
--- a/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/linearmaterialparams.hh
@@ -1,4 +1,4 @@
-// $Id: linearmaterialparams.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
index 51d7952bec..e5b4b664b5 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh
@@ -1,4 +1,4 @@
-// $Id: regularizedbrookscorey.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh
index 767a47fb09..c65b2869e3 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedbrookscoreyparams.hh
@@ -1,4 +1,4 @@
-// $Id: regularizedbrookscoreyparams.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
index 072df298db..80233098b2 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh
@@ -1,4 +1,4 @@
-// $Id: regularizedvangenuchten.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh
index bbac9ad8d8..c65ff9f7a1 100644
--- a/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchtenparams.hh
@@ -1,4 +1,4 @@
-// $Id: regularizedvangenuchtenparams.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
index 45102ffe65..88603f05c3 100644
--- a/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh
@@ -1,4 +1,4 @@
-// $Id: vangenuchten.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh b/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh
index 25181813b7..0812da7910 100644
--- a/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/vangenuchtenparams.hh
@@ -1,4 +1,4 @@
-// $Id: vangenuchtenparams.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidstate.hh b/dumux/material/fluidstate.hh
index e496438f9e..1d9dc25998 100644
--- a/dumux/material/fluidstate.hh
+++ b/dumux/material/fluidstate.hh
@@ -1,4 +1,4 @@
-// $Id: fluidstate.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
@@ -117,8 +117,8 @@ public:
*
* Unit: [Pa] = [N/m^2]
*/
- Scalar partialPressure(int componentIdx) const
- { DUNE_THROW(Dune::NotImplemented, "FluidState::partialPressure()"); }
+ Scalar fugacity(int componentIdx) const
+ { DUNE_THROW(Dune::NotImplemented, "FluidState::fugacity()"); }
/*!
* \brief Return the total pressure of the gas phase.
diff --git a/dumux/material/fluidsystems/2p_system.hh b/dumux/material/fluidsystems/2p_system.hh
index e1dab713ba..eea985f1ed 100644
--- a/dumux/material/fluidsystems/2p_system.hh
+++ b/dumux/material/fluidsystems/2p_system.hh
@@ -1,4 +1,4 @@
-// $Id: 2p_system.hh 3784 2010-06-24 13:43:57Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/fluidsystems/h2o_n2_system.hh b/dumux/material/fluidsystems/h2o_n2_system.hh
index a4154846c5..f1eb18b80c 100644
--- a/dumux/material/fluidsystems/h2o_n2_system.hh
+++ b/dumux/material/fluidsystems/h2o_n2_system.hh
@@ -1,4 +1,4 @@
-// $Id: h2o_n2_system.hh 3824 2010-07-13 13:30:02Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
@@ -91,57 +91,6 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
}
- /*!
- * \brief Given the gas phase's composition, temperature and
- * pressure, compute the partial presures of all components
- * in [Pa] and assign it to the FluidState.
- *
- * This is required for models which cannot calculate the the
- * partial pressures of the components in the gas phase from the
- * degasPressure(). To use this method, the FluidState must have a
- * setPartialPressure(componentIndex, pressure) method.
- */
- template <class FluidState>
- static void computePartialPressures(Scalar temperature,
- Scalar pg,
- FluidState &fluidState)
- {
- Scalar X1 = fluidState.massFrac(gPhaseIdx, H2OIdx);
-
- // We use the newton method for this. For the initial value we
- // assume all components to be an ideal gas
- Scalar pH2O = fluidState.moleFrac(gPhaseIdx, H2OIdx) * pg;
- Scalar eps = pg*1e-9;
-
- Scalar deltaP = pH2O;
- Valgrind::CheckDefined(pH2O);
- Valgrind::CheckDefined(deltaP);
- for (int i = 0; i < 5 && std::abs(deltaP/pg) > 1e-9; ++i) {
- Scalar f =
- H2O::gasDensity(temperature, pH2O)*(X1 - 1) +
- X1*N2::gasDensity(temperature, pg - pH2O);
-
- Scalar df_dp;
- df_dp =
- H2O::gasDensity(temperature, pH2O + eps)*(X1 - 1) +
- X1*N2::gasDensity(temperature, pg - (pH2O + eps));
- df_dp -=
- H2O::gasDensity(temperature, pH2O - eps)*(X1 - 1) +
- X1*N2::gasDensity(temperature, pg - (pH2O - eps));
- df_dp /=
- 2*eps;
-
- deltaP = - f/df_dp;
-
- pH2O += deltaP;
- Valgrind::CheckDefined(pH2O);
- Valgrind::CheckDefined(deltaP);
- }
-
- fluidState.setPartialPressure(H2OIdx, pH2O);
- fluidState.setPartialPressure(N2Idx, pg - pH2O);
- };
-
/*!
* \brief Given a phase's composition, temperature, pressure, and
* the partial pressures of all components, return its
@@ -153,16 +102,26 @@ public:
Scalar pressure,
const FluidState &fluidState)
{
- if (phaseIdx == lPhaseIdx)
- return liquidPhaseDensity_(temperature,
- pressure,
- fluidState.moleFrac(lPhaseIdx, H2OIdx),
- fluidState.moleFrac(lPhaseIdx, N2Idx));
- else if (phaseIdx == gPhaseIdx)
- return gasPhaseDensity_(temperature,
- pressure,
- fluidState.moleFrac(gPhaseIdx, H2OIdx),
- fluidState.moleFrac(gPhaseIdx, N2Idx));
+ if (phaseIdx == lPhaseIdx) {
+ // See: Ochs 2008
+ // \todo: proper citation
+ Scalar rholH2O = H2O::liquidDensity(temperature, pressure);
+ Scalar clH2O = rholH2O/H2O::molarMass();
+
+ // this assumes each nitrogen molecule displaces exactly one
+ // water molecule in the liquid
+ return
+ clH2O*(H2O::molarMass()*fluidState.moleFrac(lPhaseIdx, H2OIdx)
+ +
+ N2::molarMass()*fluidState.moleFrac(lPhaseIdx, N2Idx));
+ }
+ else if (phaseIdx == gPhaseIdx) {
+ Scalar fugH2O = fluidState.fugacity(H2OIdx);
+ Scalar fugN2 = fluidState.fugacity(N2Idx);
+ return
+ H2O::gasDensity(temperature, fugH2O) +
+ N2::gasDensity(temperature, fugN2);
+ }
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
@@ -337,23 +296,46 @@ public:
}
/*!
- * \brief Returns the fugacity coefficient of a component in a
+ * \brief Returns the activity coefficient of a component in a
* phase.
*
- * For solutions with only traces in a solvent this boils down to
- * the inverse Henry constant for the solutes and the partial
- * pressure for the solvent.
+ * We define the activity coefficent \f$\gamma_{\alpha,\kappa}\f$
+ * of component \f$\kappa\f$ by the following equation:
+ * \f[ f_\kappa = p_\alpha \gamma_{\alpha,\kappa} \f]
+ * where \f$f_\kappa\f$ is the component's fugacity and \f$p_\alpha\f$
+ * is the phase' pressure
+ *
+ * For liquids with very low miscibility this boils down to the
+ * inverse Henry constant for the solutes and the partial pressure
+ * for the solvent.
+ *
+ * For ideal gases this is equivalent to the gas pressure, in real
+ * gases it is the gas pressure times the component's fugacity
+ * coefficient.
*/
template <class FluidState>
- static Scalar fugacityCoeff(int phaseIdx,
+ static Scalar activityCoeff(int phaseIdx,
int compIdx,
Scalar temperature,
Scalar pressure,
const FluidState &state)
{
- if (phaseIdx != lPhaseIdx)
- DUNE_THROW(Dune::NotImplemented,
- "Fugacities in the gas phase");
+ if (phaseIdx == gPhaseIdx) {
+ return pressure;
+ Scalar fugH2O = std::max(1e-3, state.fugacity(H2OIdx));
+ Scalar fugN2 = std::max(1e-3, state.fugacity(N2Idx));
+ Scalar cH2O = H2O::gasDensity(temperature, fugH2O) / H2O::molarMass();
+ Scalar cN2 = N2::gasDensity(temperature, fugN2) / N2::molarMass();
+
+ Scalar alpha = (fugH2O + fugN2)/pressure;
+
+ if (compIdx == H2OIdx)
+ return fugH2O/(alpha*cH2O/(cH2O + cN2));
+ else if (compIdx == N2Idx)
+ return fugN2/(alpha*cN2/(cH2O + cN2));
+
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
switch (compIdx) {
case H2OIdx: return H2O::vaporPressure(temperature);
@@ -362,59 +344,6 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
}
- /*!
- * \brief Given a component's pressure and temperature, return its
- * density in a phase [kg/m^3].
- */
- static Scalar componentDensity(int phaseIdx,
- int compIdx,
- Scalar temperature,
- Scalar pressure)
- {
- if (phaseIdx == lPhaseIdx) {
- if (compIdx == H2OIdx)
- return H2O::liquidDensity(temperature, pressure);
- else if (compIdx == N2Idx)
- return N2::liquidDensity(temperature, pressure);
- DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
- }
- else if (phaseIdx == gPhaseIdx) {
- if (compIdx == H2OIdx) {
- return H2O::gasDensity(temperature, pressure);
- }
- else if (compIdx == N2Idx)
- return N2::gasDensity(temperature, pressure);
- DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
- }
- DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
- };
-
- /*!
- * \brief Given a component's density and temperature, return the
- * corresponding pressure in a phase [Pa].
- */
- static Scalar componentPressure(int phaseIdx,
- int compIdx,
- Scalar temperature,
- Scalar density)
- {
- if (phaseIdx == lPhaseIdx) {
- if (compIdx == H2OIdx)
- return H2O::liquidPressure(temperature, density);
- else if (compIdx == N2Idx)
- return N2::liquidPressure(temperature, density);
- DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
- }
- else if (phaseIdx == gPhaseIdx) {
- if (compIdx == H2OIdx)
- return H2O::gasPressure(temperature, density);
- else if (compIdx == N2Idx)
- return N2::gasPressure(temperature, density);
- DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
- }
- DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
- };
-
/*!
* \brief Given a phase's composition, temperature and pressure,
* return the binary diffusion coefficent for components
@@ -496,8 +425,8 @@ public:
N2::gasEnthalpy(temperature, pN2);
}
else {
- Scalar pWater = fluidState.partialPressure(H2OIdx);
- Scalar pN2 = fluidState.partialPressure(N2Idx);
+ Scalar pWater = fluidState.fugacity(H2OIdx);
+ Scalar pN2 = fluidState.fugacity(N2Idx);
Scalar result = 0;
result +=
@@ -533,7 +462,7 @@ private:
// \todo: proper citation
Scalar rholH2O = H2O::liquidDensity(T, pl);
Scalar clH2O = rholH2O/H2O::molarMass();
-
+
// this assumes each nitrogen molecule displaces exactly one
// water molecule in the liquid
return
@@ -542,69 +471,10 @@ private:
xlN2*N2::molarMass());
}
- // defect of gas density
- static inline Scalar gasDefect_(Scalar pH2O, Scalar T, Scalar pg, Scalar xgH2O, Scalar xgN2)
- {
- // this assumes that sum of the fugacities of the individual
- // components add up to the gas pressure
- return
- H2O::gasDensity(T, pH2O)/H2O::molarMass()*(xgH2O - 1) +
- xgH2O*N2::gasDensity(T, pg - pH2O)/N2::molarMass();
- }
-
static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgH2O, Scalar xgN2)
{
- // assume ideal gas for the initial condition
- Scalar pH2O = pg*xgH2O;
- Scalar delta = 1e100;
-
- // newton method. makes sure that the total pressure of the
- // gas phase is the sum of the individual gas phase fugacities
- // of the individual components
- for (int i = 0; i < 10; ++i) {
- Scalar f;
- Scalar df_dpH2O;
-
- f = gasDefect_(pH2O, T, pg, xgH2O, xgN2);
- Scalar eps = (std::abs(pg) + 1)*1e-9;
- df_dpH2O = gasDefect_(pH2O + eps, T, pg, xgH2O, xgN2);
- df_dpH2O -= gasDefect_(pH2O - eps, T, pg, xgH2O, xgN2);
- df_dpH2O /= 2*eps;
-
- delta = - f/df_dpH2O;
- pH2O += delta;
-
- if (std::abs(delta) < 1e-10*pg) {
- return H2O::gasDensity(T, pH2O) + N2::gasDensity(T, pg - pH2O);
- }
- };
- DUNE_THROW(NumericalProblem,
- "Can not calculate the gas density at "
- << " T=" << T
- << " pg=" << pg
- << " xgH2O=" << xgH2O
- << " xgN2=" << xgN2
- << "\n");
+ return H2O::gasDensity(T, pg*xgH2O) + N2::gasDensity(T, pg*xgN2);
};
-
- template <class FluidState>
- static void checkConsistentGasDensity_(Scalar rhoToTest,
- Scalar pg,
- const FluidState &fluidState)
- {
-#ifndef NDEBUG
- Scalar rho = gasPhaseDensity_(fluidState.temperature(),
- pg,
- fluidState.moleFrac(gPhaseIdx, H2OIdx),
- fluidState.moleFrac(gPhaseIdx, N2Idx));
- static const Scalar eps = 1e-8;
- if (std::abs(rhoToTest/rho - 1.0) > eps) {
- DUNE_THROW(Dune::InvalidStateException,
- "Density of gas phase is inconsistent: rho1/rho2 = "
- << rhoToTest/rho);
- }
-#endif
- }
};
} // end namepace
diff --git a/dumux/material/idealgas.hh b/dumux/material/idealgas.hh
index 9bc012277c..5592c17b5c 100644
--- a/dumux/material/idealgas.hh
+++ b/dumux/material/idealgas.hh
@@ -1,4 +1,4 @@
-// $Id: idealgas.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/settablephase.hh b/dumux/material/settablephase.hh
index 29320b5784..621ce608c9 100644
--- a/dumux/material/settablephase.hh
+++ b/dumux/material/settablephase.hh
@@ -1,4 +1,4 @@
-// $Id: settablephase.hh 3777 2010-06-24 06:46:46Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2010 by Andreas Lauser
* Institute of Hydraulic Engineering *
diff --git a/dumux/material/spatialparameters/boxspatialparameters.hh b/dumux/material/spatialparameters/boxspatialparameters.hh
index f15cdd6b3a..3a1a4b97bc 100644
--- a/dumux/material/spatialparameters/boxspatialparameters.hh
+++ b/dumux/material/spatialparameters/boxspatialparameters.hh
@@ -1,4 +1,4 @@
-// $Id: boxspatialparameters.hh 3783 2010-06-24 11:33:53Z bernd $
+// $Id$
/*****************************************************************************
* Copyright (C) 2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
@@ -39,10 +39,14 @@ class BoxSpatialParameters
{
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(SpatialParameters)) Implementation;
enum {
dimWorld = GridView::dimensionworld
};
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GridView::ctype CoordScalar;
typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
@@ -54,6 +58,41 @@ public:
~BoxSpatialParameters()
{}
+ /*!
+ * \brief Returns the factor by which the volume of a sub control
+ * volume needs to be multiplied in order to get cubic
+ * meters.
+ *
+ * By default that's just 1.0
+ */
+ Scalar extrusionFactorScv(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvIdx) const
+ { return 1.0; }
+
+ /*!
+ * \brief Returns the factor by which the area of a sub control
+ * volume face needs to be multiplied in order to get
+ * square meters.
+ *
+ * By default it is the arithmetic mean of the extrusion factor of
+ * the face's two sub-control volumes.
+ */
+ Scalar extrusionFactorScvf(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ int scvfIdx) const
+ {
+ return
+ 0.5 *
+ (asImp_().extrusionFactorScv(element,
+ fvElemGeom,
+ fvElemGeom.subContVolFace[scvfIdx].i)
+ +
+ asImp_().extrusionFactorScv(element,
+ fvElemGeom,
+ fvElemGeom.subContVolFace[scvfIdx].j));
+ }
+
/*!
* \brief Averages the intrinsic permeability.
*/
@@ -82,6 +121,13 @@ public:
for (int j = 0; j < dimWorld; ++j)
result[i][j] = harmonicMean(K1[i][j], K2[i][j]);
}
+
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
};
} // namespace Dumux
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index fcc436adc3..0029169ed7 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -1,25 +1,25 @@
-// $Id: newtoncontroller.hh 3826 2010-07-14 07:03:41Z bernd $
+// $Id$
/****************************************************************************
-* Copyright (C) 2008-2010 by Andreas Lauser *
-* Copyright (C) 2008-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. *
-*****************************************************************************/
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-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. *
+ *****************************************************************************/
/*!
-* \file
-* \brief Reference implementation of a newton controller solver.
-*
-* Usually for most cases this controller should be sufficient.
-*/
+ * \file
+ * \brief Reference implementation of a newton controller solver.
+ *
+ * Usually for most cases this controller should be sufficient.
+ */
#ifndef DUMUX_NEWTON_CONTROLLER_HH
#define DUMUX_NEWTON_CONTROLLER_HH
@@ -40,7 +40,7 @@
#include <dumux/common/pardiso.hh>
-#include <dumux/boxmodels/pdelab/preconditionerpdelab.hh>
+#include <dumux/common/pdelabpreconditioner.hh>
#include <dune/pdelab/backend/istlsolverbackend.hh>
@@ -86,12 +86,11 @@ struct NewtonConvergenceWriter
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::SolutionVector SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef Dumux::VtkMultiWriter<GridView> VtkMultiWriter;
NewtonConvergenceWriter(NewtonController &ctl)
- : ctl_(ctl)
+ : ctl_(ctl)
{
timeStepNum_ = 0;
iteration_ = 0;
@@ -111,19 +110,17 @@ struct NewtonConvergenceWriter
{
++ iteration_;
vtkMultiWriter_->beginTimestep(timeStepNum_ + iteration_ / 100.0,
- gv);
+ gv);
};
void writeFields(const SolutionVector &uOld,
const SolutionVector &deltaU)
{
- ctl_.model().localJacobian().addConvergenceVtkFields(*vtkMultiWriter_, uOld, deltaU);
+ ctl_.method().model().addConvergenceVtkFields(*vtkMultiWriter_, uOld, deltaU);
};
void endIteration()
- {
- vtkMultiWriter_->endTimestep();
- };
+ { vtkMultiWriter_->endTimestep(); };
void endTimestep()
{
@@ -143,9 +140,8 @@ struct NewtonConvergenceWriter<TypeTag, false>
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::SolutionVector SolutionVector;
typedef Dumux::VtkMultiWriter<GridView> VtkMultiWriter;
NewtonConvergenceWriter(NewtonController &ctl)
@@ -158,7 +154,7 @@ struct NewtonConvergenceWriter<TypeTag, false>
{ };
void writeFields(const SolutionVector &uOld,
- const SolutionVector &deltaU)
+ const SolutionVector &deltaU)
{ };
void endIteration()
@@ -169,36 +165,34 @@ struct NewtonConvergenceWriter<TypeTag, false>
};
/*!
-* \brief The reference implementation of a newton controller.
-*
-* If you want to specialize only some methods but are happy with
-* the defaults of the reference controller, derive your
-* controller from this class and simply overload the required
-* methods.
-*/
+ * \brief The reference implementation of a newton controller.
+ *
+ * If you want to specialize only some methods but are happy with
+ * the defaults of the reference controller, derive your
+ * controller from this class and simply overload the required
+ * methods.
+ */
template <class TypeTag>
class NewtonController
{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
- typedef typename GET_PROP(TypeTag, PTAG(PDELabTypes)) PDELabTypes;
- typedef typename PDELabTypes::GridFunctionSpace GridFunctionSpace;
- typedef typename PDELabTypes::ConstraintsTrafo ConstraintsTrafo;
-
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::SolutionVector SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef NewtonConvergenceWriter<TypeTag, GET_PROP_VALUE(TypeTag, PTAG(NewtonWriteConvergence))> ConvergenceWriter;
public:
NewtonController()
- : endIterMsgStream_(std::ostringstream::out),
- convergenceWriter_(asImp_())
+ : endIterMsgStream_(std::ostringstream::out),
+ convergenceWriter_(asImp_())
{
numSteps_ = 0;
@@ -233,8 +227,8 @@ public:
{ maxSteps_ = maxSteps; }
/*!
- * \brief Returns true if another iteration should be done.
- */
+ * \brief Returns true if another iteration should be done.
+ */
bool newtonProceed(SolutionVector &u)
{
if (numSteps_ < 2)
@@ -250,7 +244,7 @@ public:
return false; // we are below the desired tolerance
Scalar tmp = asImp_().physicalness_(u);
- curPhysicalness_ = model().gridView().comm().min(tmp);
+ curPhysicalness_ = gridView_().comm().min(tmp);
curPhysicalness_ = std::min(curPhysicalness_, Scalar(1.0));
@@ -293,21 +287,21 @@ public:
}
/*!
- * \brief Returns true iff the error of the solution is below the
- * tolerance.
- */
+ * \brief Returns true iff the error of the solution is below the
+ * tolerance.
+ */
bool newtonConverged() const
{
return (error_ <= tolerance_) && (curPhysicalness_ >= 1.0);
}
/*!
- * \brief Called before the newton method is applied to an
- * non-linear system of equations.
- */
- void newtonBegin(NewtonMethod *method, SolutionVector &u)
+ * \brief Called before the newton method is applied to an
+ * non-linear system of equations.
+ */
+ void newtonBegin(NewtonMethod &method, SolutionVector &u)
{
- method_ = method;
+ method_ = &method;
numSteps_ = 0;
probationCount_ = 0;
maxPhysicalness_ = 0;
@@ -317,23 +311,23 @@ public:
}
/*!
- * \brief Indidicates the beginning of a newton iteration.
- */
+ * \brief Indidicates the beginning of a newton iteration.
+ */
void newtonBeginStep()
{ lastError_ = error_; }
/*!
- * \brief Returns the number of steps done since newtonBegin() was
- * called.
- */
+ * \brief Returns the number of steps done since newtonBegin() was
+ * called.
+ */
int newtonNumSteps()
{ return numSteps_; }
/*!
- * \brief Update the error of the solution compared to the
- * previous iteration.
- */
+ * \brief Update the error of the solution compared to the
+ * previous iteration.
+ */
void newtonUpdateRelError(const SolutionVector &uOld,
const SolutionVector &deltaU)
{
@@ -352,23 +346,23 @@ public:
/ std::max(std::abs(uOld[i][j]), Scalar(1e-4));
if (tmp > error_)
{
- idxI = i;
- idxJ = j;
+ idxI = i;
+ idxJ = j;
+ error_ = tmp;
}
- error_ = std::max(error_, tmp);
}
}
- error_ = model().gridView().comm().max(error_);
+ error_ = gridView_().comm().max(error_);
}
/*!
- * \brief Solve the linear system of equations \f$ \mathbf{A}x - b
- * = 0\f$.
- *
- * Throws Dumux::NumericalProblem if the linear solver didn't
- * converge.
- */
+ * \brief Solve the linear system of equations \f$ \mathbf{A}x - b
+ * = 0\f$.
+ *
+ * Throws Dumux::NumericalProblem if the linear solver didn't
+ * converge.
+ */
template <class Matrix, class Vector>
void newtonSolveLinear(Matrix &A,
Vector &u,
@@ -382,21 +376,21 @@ public:
Scalar residReduction = 1e-6;
try {
- solveLinear_(A, u, b, residReduction);
+ solveLinear_(A, u, b, residReduction);
- // make sure all processes converged
- int converged = 1;
- model().gridView().comm().min(converged);
+ // make sure all processes converged
+ int converged = 1;
+ gridView_().comm().min(converged);
- if (!converged) {
- DUNE_THROW(NumericalProblem,
- "A process threw MatrixBlockError");
- }
+ if (!converged) {
+ DUNE_THROW(NumericalProblem,
+ "A process threw MatrixBlockError");
+ }
}
catch (Dune::MatrixBlockError e) {
// make sure all processes converged
int converged = 0;
- model().gridView().comm().min(converged);
+ gridView_().comm().min(converged);
Dumux::NumericalProblem p;
std::string msg;
@@ -408,20 +402,20 @@ public:
}
/*!
- * \brief Update the current solution function with a delta vector.
- *
- * The error estimates required for the newtonConverged() and
- * newtonProceed() methods should be updated here.
- *
- * Different update strategies, such as line search and chopped
- * updates can be implemented. The default behaviour is just to
- * subtract deltaU from uOld.
- *
- * \param deltaU The delta as calculated from solving the linear
- * system of equations. This parameter also stores
- * the updated solution.
- * \param uOld The solution of the last iteration
- */
+ * \brief Update the current solution function with a delta vector.
+ *
+ * The error estimates required for the newtonConverged() and
+ * newtonProceed() methods should be updated here.
+ *
+ * Different update strategies, such as line search and chopped
+ * updates can be implemented. The default behaviour is just to
+ * subtract deltaU from uOld.
+ *
+ * \param deltaU The delta as calculated from solving the linear
+ * system of equations. This parameter also stores
+ * the updated solution.
+ * \param uOld The solution of the last iteration
+ */
void newtonUpdate(SolutionVector &deltaU, const SolutionVector &uOld)
{
writeConvergence_(uOld, deltaU);
@@ -433,52 +427,52 @@ public:
}
/*!
- * \brief Indicates that one newton iteration was finished.
- */
+ * \brief Indicates that one newton iteration was finished.
+ */
void newtonEndStep(SolutionVector &u, SolutionVector &uOld)
{
++numSteps_;
curPhysicalness_ = asImp_().physicalness_(u);
- if (this->method().verbose())
+ if (verbose())
std::cout << "\rNewton iteration " << numSteps_ << " done: "
- << "error=" << error_ << endIterMsg().str() << "\n";
+ << "error=" << error_ << endIterMsg().str() << "\n";
endIterMsgStream_.str("");
}
/*!
- * \brief Indicates that we're done solving the non-linear system of equations.
- */
+ * \brief Indicates that we're done solving the non-linear system of equations.
+ */
void newtonEnd()
{
convergenceWriter_.endTimestep();
}
/*!
- * \brief Called if the newton method broke down.
- *
- * This method is called _after_ newtonEnd()
- */
+ * \brief Called if the newton method broke down.
+ *
+ * This method is called _after_ newtonEnd()
+ */
void newtonFail()
{
numSteps_ = targetSteps_*2;
}
/*!
- * \brief Called when the newton method was sucessful.
- *
- * This method is called _after_ newtonEnd()
- */
+ * \brief Called when the newton method was sucessful.
+ *
+ * This method is called _after_ newtonEnd()
+ */
void newtonSucceed()
{ }
/*!
- * \brief Suggest a new time stepsize based on the old time step size.
- *
- * The default behaviour is to suggest the old time step size
- * scaled by the ratio between the target iterations and the
- * iterations required to actually solve the last time step.
- */
+ * \brief Suggest a new time stepsize based on the old time step size.
+ *
+ * The default behaviour is to suggest the old time step size
+ * scaled by the ratio between the target iterations and the
+ * iterations required to actually solve the last time step.
+ */
Scalar suggestTimeStepSize(Scalar oldTimeStep) const
{
// be agressive reducing the timestep size but
@@ -493,42 +487,66 @@ public:
else {
/*Scalar percent = (Scalar(1))/targetSteps_;
return oldTimeStep*(1 + percent);
- */
+ */
Scalar percent = ((Scalar) targetSteps_ - numSteps_)/targetSteps_;
return oldTimeStep*(1.0 + percent/1.2);
}
}
/*!
- * \brief Returns a reference to the current newton method
- * which is controlled by this controller.
- */
+ * \brief Returns a reference to the current newton method
+ * which is controlled by this controller.
+ */
NewtonMethod &method()
{ return *method_; }
/*!
- * \brief Returns a reference to the current newton method
- * which is controlled by this controller.
- */
+ * \brief Returns a reference to the current newton method
+ * which is controlled by this controller.
+ */
const NewtonMethod &method() const
{ return *method_; }
+ std::ostringstream &endIterMsg()
+ { return endIterMsgStream_; }
+
/*!
- * \brief Returns a reference to the current numeric model.
- */
- Model &model()
- { return method_->model(); }
+ * \brief Returns true iff the newton method ought to be chatty.
+ */
+ bool verbose() const
+ { return gridView_().comm().rank() == 0; }
+protected:
/*!
- * \brief Returns a reference to the current numeric model.
- */
- const Model &model() const
- { return method_->model(); }
+ * \brief Returns a reference to the grid view.
+ */
+ const GridView &gridView_() const
+ { return problem_().gridView(); }
- std::ostringstream &endIterMsg()
- { return endIterMsgStream_; }
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ Problem &problem_()
+ { return method_->problem(); }
+
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ const Problem &problem_() const
+ { return method_->problem(); }
+
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ Model &model_()
+ { return problem_().model(); }
+
+ /*!
+ * \brief Returns a reference to the problem.
+ */
+ const Model &model_() const
+ { return problem_().model(); }
-protected:
// returns the actual implementation for the cotroller we do
// it this way in order to allow "poor man's virtual methods",
// i.e. methods of subclasses which can be called by the base
@@ -541,7 +559,7 @@ protected:
void writeConvergence_(const SolutionVector &uOld,
const SolutionVector &deltaU)
{
- convergenceWriter_.beginIteration(this->model().gridView());
+ convergenceWriter_.beginIteration(this->gridView_());
convergenceWriter_.writeFields(uOld, deltaU);
convergenceWriter_.endIteration();
};
@@ -556,25 +574,25 @@ protected:
Vector &x = u;
int verbosity = GET_PROP_VALUE(TypeTag, PTAG(NewtonLinearSolverVerbosity));
- if (model().gridView().grid().comm().rank() != 0)
+ if (gridView_().comm().rank() != 0)
verbosity = 0;
#if HAVE_PARDISO
- typedef Dune::PDELab::ISTLBackend_NoOverlap_Loop_Pardiso<TypeTag> Solver;
- Solver solver(model(), 500, verbosity);
+ typedef Dumux::PDELab::ISTLBackend_NoOverlap_Loop_Pardiso<TypeTag> Solver;
+ Solver solver(problem_(), 500, verbosity);
#else // !HAVE_PARDISO
#if HAVE_MPI
// typedef Dune::PDELab::ISTLBackend_NOVLP_BCGS_NOPREC<GridFunctionSpace> Solver;
-// Solver solver(model().jacobianAssembler().gridFunctionSpace(), 50000, verbosity);
- typedef Dune::PDELab::ISTLBackend_NoOverlap_BCGS_ILU<TypeTag> Solver;
- Solver solver(model(), 500, verbosity);
+// Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 50000, verbosity);
+ typedef Dumux::PDELab::ISTLBackend_NoOverlap_BCGS_ILU<TypeTag> Solver;
+ Solver solver(problem_(), 500, verbosity);
#else
- typedef Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR Solver;
+ typedef Dumux::PDELab::ISTLBackend_SEQ_BCGS_SSOR Solver;
Solver solver(500, verbosity);
#endif // HAVE_MPI
#endif // HAVE_PARDISO
- // Solver solver(model().jacobianAssembler().gridFunctionSpace(), 500, verbosity);
+ // Solver solver(model_().jacobianAssembler().gridFunctionSpace(), 500, verbosity);
solver.apply(A, x, b, residReduction);
if (!solver.result().converged)
@@ -585,7 +603,7 @@ protected:
// somewhere. this should never happen but for some strange
// reason it happens anyway.
Scalar xNorm2 = x.two_norm2();
- model().gridView().grid().comm().sum(xNorm2);
+ gridView_().comm().sum(xNorm2);
if (std::isnan(xNorm2) || !std::isfinite(xNorm2))
DUNE_THROW(Dumux::NumericalProblem,
"The linear solver produced a NaN or inf somewhere.");
@@ -629,8 +647,6 @@ protected:
// actual number of steps done so far
int numSteps_;
};
-
-}
-
+} // namespace Dumux
#endif
diff --git a/dumux/nonlinear/newtonmethod.hh b/dumux/nonlinear/newtonmethod.hh
index b8a300c162..3ff95481b7 100644
--- a/dumux/nonlinear/newtonmethod.hh
+++ b/dumux/nonlinear/newtonmethod.hh
@@ -1,4 +1,4 @@
-// $Id: newtonmethod.hh 3738 2010-06-15 14:01:09Z lauser $
+// $Id$
/*****************************************************************************
* Copyright (C) 2008 by Andreas Lauser, Bernd Flemisch *
* Institute of Hydraulic Engineering *
@@ -23,10 +23,6 @@
#ifndef DUMUX_NEWTONMETHOD_HH
#define DUMUX_NEWTONMETHOD_HH
-#ifdef DUMUX_NEWTONMETHOD_DEPRECATED_HH
-# error "Never use the old and the new newton method in one program!"
-#endif
-
#include <limits>
#include <dumux/common/exceptions.hh>
@@ -44,54 +40,59 @@ template <class TypeTag>
class NewtonMethod
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalJacobian)) LocalJacobian;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) NewtonController;
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::SolutionVector SolutionVector;
- typedef typename SolutionTypes::JacobianAssembler JacobianAssembler;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianAssembler)) JacobianAssembler;
public:
- NewtonMethod(Model &model)
- : model_(model)
+ NewtonMethod(Problem &problem)
+ : problem_(problem)
{ }
/*!
- * \brief Returns a reference to the current numeric model.
+ * \brief Returns a reference to the current numeric problem.
*/
- Model &model()
- { return model_; }
+ Problem &problem()
+ { return problem_; }
/*!
- * \brief Returns a reference to the current numeric model.
+ * \brief Returns a reference to the current numeric problem.
*/
- const Model &model() const
- { return model_; }
+ const Problem &problem() const
+ { return problem_; }
+
+ /*!
+ * \brief Returns a reference to the numeric model.
+ */
+ Model &model()
+ { return problem().model(); }
/*!
- * \brief Returns true iff the newton method should be verbose
- * about what it is going on;
+ * \brief Returns a reference to the numeric model.
*/
- bool verbose() const
- { return (model().gridView().comm().rank() == 0); }
+ const Model &model() const
+ { return problem().model(); }
+
/*!
* \brief Run the newton method. The controller is responsible
* for all the strategic decisions.
*/
- template <class NewtonController>
bool execute(NewtonController &ctl)
{
try {
return execute_(ctl);
}
catch (const Dune::ISTLError &e) {
- if (verbose())
+ if (ctl.verbose())
std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
ctl.newtonFail();
return false;
}
catch (const Dumux::NumericalProblem &e) {
- if (verbose())
+ if (ctl.verbose())
std::cout << "Newton: Caught exception: \"" << e.what() << "\"\n";
ctl.newtonFail();
return false;
@@ -99,15 +100,14 @@ public:
}
protected:
- template <class NewtonController>
bool execute_(NewtonController &ctl)
{
// TODO (?): u shouldn't be hard coded to the model
- SolutionVector &u = model_.curSol();
- JacobianAssembler &jacobianAsm = model_.jacobianAssembler();
+ SolutionVector &u = model().curSol();
+ JacobianAssembler &jacobianAsm = model().jacobianAssembler();
// tell the controller that we begin solving
- ctl.newtonBegin(this, u);
+ ctl.newtonBegin(*this, u);
// execute the method as long as the controller thinks
// that we should do another iteration
@@ -120,7 +120,7 @@ protected:
// make the current solution to the old one
uOld_ = u;
- if (verbose()) {
+ if (ctl.verbose()) {
std::cout << "Assembling global jacobian";
std::cout.flush();
}
@@ -128,7 +128,7 @@ protected:
jacobianAsm.assemble(u);
// solve the resultuing linear equation system
- if (verbose()) {
+ if (ctl.verbose()) {
std::cout << "\rSolve Mx = r ";
std::cout.flush();
}
@@ -161,14 +161,13 @@ protected:
return true;
}
-
private:
SolutionVector uOld_;
SolutionVector residual_;
- Model &model_;
- bool verbose_;
+ Problem &problem_;
};
+
}
#endif
--
GitLab