From 7e8cf7cc5c6a161a7f35cbf625f15f382b39454a Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 27 Oct 2011 15:12:06 +0000
Subject: [PATCH] moved unused headers from stable to devel
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6797 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../2p2cni/2p2cniboundaryvariables.hh | 271 ---------
dumux/boxmodels/pdelab/pdelabadapter.hh | 170 ------
.../pdelab/pdelabboxistlvectorbackend.hh | 143 -----
.../pdelab/pdelabboxlocaloperator.hh | 147 -----
dumux/common/pdelabpreconditioner.hh | 522 ------------------
dumux/material/binarycoefficients/h2_n2.hh | 99 ----
.../Mp/2poftadapter.hh | 83 ---
.../fluidsystems/simplefluidsystem.hh | 326 -----------
8 files changed, 1761 deletions(-)
delete mode 100644 dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh
delete mode 100644 dumux/boxmodels/pdelab/pdelabadapter.hh
delete mode 100644 dumux/boxmodels/pdelab/pdelabboxistlvectorbackend.hh
delete mode 100644 dumux/boxmodels/pdelab/pdelabboxlocaloperator.hh
delete mode 100644 dumux/common/pdelabpreconditioner.hh
delete mode 100644 dumux/material/binarycoefficients/h2_n2.hh
delete mode 100644 dumux/material/fluidmatrixinteractions/Mp/2poftadapter.hh
delete mode 100644 dumux/material/fluidsystems/simplefluidsystem.hh
diff --git a/dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh b/dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh
deleted file mode 100644
index 19a054e025..0000000000
--- a/dumux/boxmodels/2p2cni/2p2cniboundaryvariables.hh
+++ /dev/null
@@ -1,271 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2010 by Katherina Baber, Klaus Mosthaf *
- * 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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This file contains the data which is required to calculate
- * all fluxes of the fluid phase over the boundary of a finite volume.
- *
- * This means pressure and temperature gradients, phase densities at
- * the integration point of the boundary, etc.
- */
-#ifndef DUMUX_2P2CNI_BOUNDARY_VARIABLES_HH
-#define DUMUX_2P2CNI_BOUNDARY_VARIABLES_HH
-
-#include <dumux/common/math.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup TwoPTwoCNIModel
- * \brief This template class contains the data which is required to
- * calculate the fluxes of the fluid phases over the boundary of a
- * finite volume for the 2p2cni model.
- *
- * This means pressure and velocity gradients, phase density and viscosity at
- * the integration point of the boundary, etc.
- */
-template <class TypeTag>
-class TwoPTwoCNIBoundaryVariables
-{
- 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(Problem)) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPTwoCIndices)) Indices;
-
- enum {
- dim = GridView::dimension,
- };
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef Dune::FieldVector<Scalar, dim> VelocityVector;
- typedef Dune::FieldVector<Scalar, dim> ScalarGradient;
- typedef Dune::FieldMatrix<Scalar, dim, dim> VectorGradient;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
- typedef typename FVElementGeometry::SubControlVolume SCV;
- typedef typename FVElementGeometry::BoundaryFace BoundaryFace;
-
- enum {
- lPhaseIdx = Indices::lPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
-
- lCompIdx = Indices::lCompIdx,
- gCompIdx = Indices::gCompIdx,
-
- numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases))
- };
-
-public:
- TwoPTwoCNIBoundaryVariables(const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int boundaryFaceIdx,
- const ElementVolumeVariables &elemDat,
- int scvIdx)
- : fvElemGeom_(elemGeom), scvIdx_(scvIdx)
- {
- boundaryFace_ = &fvElemGeom_.boundaryFace[boundaryFaceIdx];
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- densityAtIP_[phaseIdx] = Scalar(0);
- molarDensityAtIP_[phaseIdx] = Scalar(0);
- potentialGrad_[phaseIdx] = Scalar(0);
- concentrationGrad_[phaseIdx] = Scalar(0);
- molarConcGrad_[phaseIdx] = Scalar(0);
- }
-
- calculateBoundaryValues_(problem, element, elemDat);
- };
-
- Scalar KmvpNormal(int phaseIdx) const
- { return KmvpNormal_[phaseIdx]; }
-
- /*!
- * \brief The binary diffusion coefficient for each fluid phase.
- */
- Scalar porousDiffCoeff(int phaseIdx) const
- { return porousDiffCoeff_[phaseIdx]; };
-
- /*!
- * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
- * point.
- */
- Scalar densityAtIP(int phaseIdx) const
- { return densityAtIP_[phaseIdx]; }
-
- /*!
- * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration
- * point.
- */
- Scalar molarDensityAtIP(int phaseIdx) const
- { return molarDensityAtIP_[phaseIdx]; }
-
- /*!
- * \brief The concentration gradient of a component in a phase.
- */
- const ScalarGradient &concentrationGrad(int phaseIdx) const
- { return concentrationGrad_[phaseIdx]; };
-
- /*!
- * \brief The molar concentration gradient of a component in a phase.
- */
- const ScalarGradient &molarConcGrad(int phaseIdx) const
- { return molarConcGrad_[phaseIdx]; };
-
- /*!
- * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
- * of the rock matrix over the sub-control volume's face in
- * direction of the face normal.
- */
- Scalar normalMatrixHeatFlux() const
- { return normalMatrixHeatFlux_; }
-
- const BoundaryFace& boundaryFace() const
- { return *boundaryFace_; }
-
-protected:
- void calculateBoundaryValues_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemDat)
- {
- ScalarGradient tmp(0.0);
-
- // calculate gradients and secondary variables at IPs of the boundary
- for (int idx = 0;
- idx < fvElemGeom_.numVertices;
- idx++) // loop over adjacent vertices
- {
- // FE gradient at vertex idx
- const ScalarGradient& feGrad = boundaryFace_->grad[idx];
-
- // compute sum of pressure gradients for each phase
- for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
- {
- // the pressure gradient
- tmp = feGrad;
- tmp *= elemDat[idx].pressure(phaseIdx);
- potentialGrad_[phaseIdx] += tmp;
- }
-
- // the concentration gradient of the non-wetting
- // component in the wetting phase
- tmp = feGrad;
- tmp *= elemDat[idx].fluidState().massFrac(lPhaseIdx, gCompIdx);
- concentrationGrad_[lPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemDat[idx].fluidState().moleFrac(lPhaseIdx, gCompIdx);
- molarConcGrad_[lPhaseIdx] += tmp;
-
- // the concentration gradient of the wetting component
- // in the non-wetting phase
- tmp = feGrad;
- tmp *= elemDat[idx].fluidState().massFrac(gPhaseIdx, lCompIdx);
- concentrationGrad_[gPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemDat[idx].fluidState().moleFrac(gPhaseIdx, lCompIdx);
- molarConcGrad_[gPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemDat[idx].temperature();
- temperatureGrad_ += tmp;
-
- for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
- {
- densityAtIP_[phaseIdx] += elemDat[idx].density(phaseIdx)*boundaryFace_->shapeValue[idx];
- molarDensityAtIP_[phaseIdx] += elemDat[idx].molarDensity(phaseIdx)*boundaryFace_->shapeValue[idx];
- }
- }
-
- for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
- {
- tmp = problem.gravity();
- tmp *= densityAtIP_[phaseIdx];
-
- potentialGrad_[phaseIdx] -= tmp;
-
- Scalar k = problem.spatialParameters().intrinsicPermeability(element, fvElemGeom_, scvIdx_);
- VectorGradient K(0);
- K[0][0] = K[1][1] = k;
- ScalarGradient Kmvp;
- K.mv(potentialGrad_[phaseIdx], Kmvp);
- KmvpNormal_[phaseIdx] = - (Kmvp*boundaryFace_->normal);
-
- const VolumeVariables &vertDat = elemDat[scvIdx_];
-
- if (vertDat.saturation(phaseIdx) <= 0)
- porousDiffCoeff_[phaseIdx] = 0.0;
- else
- {
- Scalar tau = 1.0/(vertDat.porosity()*vertDat.porosity())*
- pow(vertDat.porosity()*vertDat.saturation(phaseIdx), 7.0/3);
-
- porousDiffCoeff_[phaseIdx] = vertDat.porosity()*vertDat.saturation(phaseIdx)*tau*vertDat.diffCoeff(phaseIdx);
- }
- }
-
- // The spatial parameters calculates the actual heat flux vector
- problem.spatialParameters().boundaryMatrixHeatFlux(tmp,
- elemDat,
- temperatureGrad_,
- element,
- fvElemGeom_,
- scvIdx_);
- // project the heat flux vector on the face's normal vector
- normalMatrixHeatFlux_ = tmp*boundaryFace_->normal;
-
- }
-
- const FVElementGeometry &fvElemGeom_;
- const BoundaryFace *boundaryFace_;
-
- // gradients
- ScalarGradient potentialGrad_[numPhases];
- ScalarGradient concentrationGrad_[numPhases];
- ScalarGradient molarConcGrad_[numPhases];
-
- // density of each face at the integration point
- Scalar densityAtIP_[numPhases];
- Scalar molarDensityAtIP_[numPhases];
-
- // intrinsic permeability times pressure potential gradient
- // projected on the face normal
- Scalar KmvpNormal_[numPhases];
-
- // the diffusion coefficient for the porous medium
- Scalar porousDiffCoeff_[numPhases];
-
- ScalarGradient temperatureGrad_;
- Scalar normalMatrixHeatFlux_;
-
- int scvIdx_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/boxmodels/pdelab/pdelabadapter.hh b/dumux/boxmodels/pdelab/pdelabadapter.hh
deleted file mode 100644
index 9d97342355..0000000000
--- a/dumux/boxmodels/pdelab/pdelabadapter.hh
+++ /dev/null
@@ -1,170 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2011 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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This file provides wrappers which allow the Dumux box models to
- * be used with dune pdelab
- */
-#ifndef DUMUX_PDELAB_ADAPTER_HH
-#define DUMUX_PDELAB_ADAPTER_HH
-
-#if ! HAVE_DUNE_PDELAB
-#error "DUNE-PDELab must be available in order to include this file!"
-#endif
-
-#include <dune/pdelab/backend/istlmatrixbackend.hh>
-#include <dune/pdelab/finiteelementmap/q1fem.hh>
-#include <dune/pdelab/finiteelementmap/conformingconstraints.hh>
-#include <dune/pdelab/backend/istlvectorbackend.hh>
-#include "pdelabboxistlvectorbackend.hh"
-#include "pdelabboxlocaloperator.hh"
-
-namespace Dumux
-{
-namespace Properties
-{
-/*!
- * \ingroup ModelCoupling
- */
-// \{
-
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tag for box problems which ought to be used in conjunction
-//! with PDELab
-NEW_TYPE_TAG(BoxPDELab);
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-//! Specifies the host grid
-NEW_PROP_TAG(Grid);
-
-//! Specifies the scalar grid function space used for sub-problems
-NEW_PROP_TAG(ScalarGridFunctionSpace);
-
-//! Specifies the grid function space used for sub-problems
-NEW_PROP_TAG(GridFunctionSpace);
-
-//! Specifies the grid operator used for sub-problems
-NEW_PROP_TAG(GridOperator);
-
-//! Specifies the grid operator space used for sub-problems
-NEW_PROP_TAG(GridOperatorSpace);
-
-//! Specifies the type of the constraints
-NEW_PROP_TAG(Constraints);
-
-//! Specifies the type of the constraints transformation
-NEW_PROP_TAG(ConstraintsTrafo);
-
-//! Specifies the local finite element space
-NEW_PROP_TAG(LocalFEMSpace);
-
-//! Specifies the local operator
-NEW_PROP_TAG(LocalOperator);
-
-}
-}
-
-namespace Dumux
-{
-namespace Properties
-{
-// set the grid functions space for the sub-models
-SET_PROP(BoxPDELab, ScalarGridFunctionSpace)
-{private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Constraints)) Constraints;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
-public:
- typedef Dune::PDELab::GridFunctionSpace<GridView, FEM, Constraints, Dumux::PDELab::BoxISTLVectorBackend<TypeTag> > type;
-};
-
-// set the grid functions space for the sub-models
-SET_PROP(BoxPDELab, GridFunctionSpace)
-{private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ScalarGridFunctionSpace)) ScalarGridFunctionSpace;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
-public:
- typedef Dune::PDELab::PowerGridFunctionSpace<ScalarGridFunctionSpace, numEq, Dune::PDELab::GridFunctionSpaceBlockwiseMapper> type;
-};
-
-// set the grid function space for the sub-models
-SET_TYPE_PROP(BoxPDELab, Constraints, Dune::PDELab::NoConstraints);
-
-// set the grid functions space for the sub-models
-SET_PROP(BoxPDELab, ConstraintsTrafo)
-{private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
-public:
- typedef typename GridFunctionSpace::template ConstraintsContainer<Scalar>::Type type;
-};
-
-// set the local operator used for submodels
-SET_PROP(BoxPDELab, LocalOperator)
-{ typedef Dumux::PDELab::BoxLocalOperator<TypeTag> type; };
-
-// set the grid operator space used for submodels
-// DEPRECATED: use GridOperator instead
-SET_PROP(BoxPDELab, GridOperatorSpace)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
-
-#warning HACK: first line does not work. but why???
- //typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalOperator)) LocalOperator;
- typedef Dumux::PDELab::BoxLocalOperator<TypeTag> LocalOperator;
-
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
-
-public:
-
- typedef Dune::PDELab::GridOperatorSpace<GridFunctionSpace,
- GridFunctionSpace,
- LocalOperator,
- ConstraintsTrafo,
- ConstraintsTrafo,
- Dune::PDELab::ISTLBCRSMatrixBackend<numEq, numEq>,
- true
- > type;
-};
-
-//! use the local FEM space associated with cubes by default
-SET_PROP(BoxPDELab, 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;
-};
-
-}
-}
-
-#endif
diff --git a/dumux/boxmodels/pdelab/pdelabboxistlvectorbackend.hh b/dumux/boxmodels/pdelab/pdelabboxistlvectorbackend.hh
deleted file mode 100644
index a9cd347883..0000000000
--- a/dumux/boxmodels/pdelab/pdelabboxistlvectorbackend.hh
+++ /dev/null
@@ -1,143 +0,0 @@
-/*****************************************************************************
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This is basically a copy of PDELab's ISTLVectorBackend which
- * allows to set the global function space later
- */
-#ifndef DUMUX_BOX_SOLUTION_VECTOR_HH
-#define DUMUX_BOX_SOLUTION_VECTOR_HH
-
-#include <dumux/boxmodels/common/boxproperties.hh>
-
-#include <dune/common/fvector.hh>
-#include <dune/istl/bvector.hh>
-
-namespace Dumux {
-namespace Properties {
-NEW_PROP_TAG(NumEq);
-NEW_PROP_TAG(Model);
-};
-
-namespace PDELab {
-
-/*!
- * \brief ISTL backend for FunctionSpace
- *
- * This is basically a copy of PDELab's ISTLVectorBackend which allows
- * to set the global function space later
- *
- * \internal
- */
-template<class TypeTag>
-class BoxISTLVectorBackend
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
-
-public:
- //! \brief export the block size
- static const int BlockSize = GET_PROP_VALUE(TypeTag, PTAG(NumEq));
-
- /*!
- * \brief Container construction
- *
- * \internal
- */
- template<typename T, typename E>
- class VectorContainer : public Dune::BlockVector< Dune::FieldVector<E,BlockSize> >
- {
- public:
- typedef E ElementType;
- typedef Dune::BlockVector< Dune::FieldVector<E,BlockSize> > BaseT;
- typedef BoxISTLVectorBackend<TypeTag> Backend;
-
- VectorContainer()
- {}
- VectorContainer(const Model &model) : BaseT(model.numDofs()) {}
- VectorContainer(const Model &model, E val) : BaseT(model.numDofs(), val) {}
- VectorContainer (const T& t_) : BaseT(t_.globalSize()/BlockSize)
- {}
- VectorContainer (const T& t_, const E& e) : BaseT(t_.globalSize()/BlockSize)
- {
- BaseT::operator=(e);
- }
- VectorContainer& operator= (const E& e)
- {
- BaseT::operator=(e);
- return *this;
- }
-
- // for debugging and AMG access
- BaseT& base ()
- {
- return *this;
- }
-
- const BaseT& base () const
- {
- return *this;
- }
-
- template<typename X>
- void std_copy_to (std::vector<X>& x) const
- {
- size_t n = this->size()*BlockSize;
- x.resize(n);
- for (size_t i=0; i<n; i++)
- x[i] = (*this)[i/BlockSize][i%BlockSize];
- }
-
- template<typename X>
- void std_copy_from (const std::vector<X>& x)
- {
- size_t n = this->size()*BlockSize;
- x.resize(n);
- for (size_t i=0; i<n; i++)
- (*this)[i/BlockSize][i%BlockSize] = x[i];
- }
- };
-
- //! extract type of container element
- template<class C>
- struct Value
- {
- typedef typename C::field_type Type;
- };
-
- //! The size type
- typedef typename Dune::BlockVector< Dune::FieldVector<float,BlockSize> >::size_type size_type;
-
- // get const_reference to container element
- // note: this method does not depend on T!
- template<typename C>
- static const typename C::field_type& access (const C& c, size_type i)
- {
- return c[i/BlockSize][i%BlockSize];
- }
-
- // get non const_reference to container element
- // note: this method does not depend on T!
- template<typename C>
- static typename C::field_type& access (C& c, size_type i)
- {
- return c[i/BlockSize][i%BlockSize];
- }
-};
-
-} // namespace PDELab
-} // namespace Dune
-
-#endif
diff --git a/dumux/boxmodels/pdelab/pdelabboxlocaloperator.hh b/dumux/boxmodels/pdelab/pdelabboxlocaloperator.hh
deleted file mode 100644
index aa707dc4a2..0000000000
--- a/dumux/boxmodels/pdelab/pdelabboxlocaloperator.hh
+++ /dev/null
@@ -1,147 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief A local operator for PDELab which wraps the box models.
- */
-#ifndef DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
-#define DUMUX_PDELAB_BOX_LOCAL_OPERATOR_HH
-
-#if ! HAVE_DUNE_PDELAB
-#error "DUNE-PDELab must be available in order to include this file!"
-#endif
-
-#include<dune/pdelab/localoperator/pattern.hh>
-#include<dune/pdelab/localoperator/flags.hh>
-
-#include <dumux/boxmodels/common/boxproperties.hh>
-
-namespace Dumux {
-
-namespace PDELab {
-
-/*!
- * \brief A local operator for PDELab which wraps the box models.
- */
-template<class TypeTag>
-class BoxLocalOperator
- :
-// : public Dune::PDELab::NumericalJacobianApplyVolume<BoxLocalOperatorPDELab<TypeTag> >,
-//public Dune::PDELab::NumericalJacobianVolume<BoxLocalOperatorPDELab<TypeTag> >,
- public Dune::PDELab::FullVolumePattern,
- public Dune::PDELab::LocalOperatorDefaultFlags
-{
- // copying the local operator for PDELab is not a good idea
- BoxLocalOperator(const BoxLocalOperator &);
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
- enum{numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq))};
-
-public:
- // pattern assembly flags
- enum { doPatternVolume = true };
-
- // residual assembly flags
- enum { doAlphaVolume = true };
-
-
- /*!
- * \param model The physical model for the box scheme.
- */
- BoxLocalOperator(Model &model)
- : model_(model)
- {}
-
- /*!
- * \brief Volume integral depending on test and ansatz functions
- *
- * \tparam EG The entity geometry type from PDELab
- * \tparam LFSU The type of the local function space of the ansatz functions
- * \tparam X The type of the container for the coefficients for the ansatz functions
- * \tparam LFSV The type of the local function space of the test functions
- * \tparam R The range type (usually FieldVector<double>)
- *
- * \param eg The entity geometry object
- * \param lfsu The local function space object of the ansatz functions
- * \param x The object of the container for the coefficients for the ansatz functions
- * \param lfsv The local function space object of the test functions
- * \param r The object storing the volume integral
- */
- 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
- {
- typedef typename LFSU::Traits::SizeType size_type;
-
- //enum { Green = JacobianAssembler::Green };
- //if (model_.jacobianAssembler().elementColor(eg.entity()) == Green)
- // Green elements don't need to be reassembled
- // return;
-
- model_.localResidual().eval(eg.entity());
-
- int numVertices = x.size()/numEq;
- for (size_type comp = 0; comp < r.size(); comp++)
- r[comp] = model_.localResidual().residual(comp%numVertices)[comp/numVertices];
- }
-
- /*!
- * \brief Jacobian of volume term
- *
- * \tparam EG The entity geometry type from PDELab
- * \tparam LFSU The type of the local function space of the ansatz functions
- * \tparam X The type of the container for the coefficients for the ansatz functions
- * \tparam LFSV The type of the local function space of the test functions
- * \tparam R The range type (usually FieldVector<double>)
- *
- * \param eg The entity geometry object
- * \param lfsu The local function space object of the ansatz functions
- * \param x The object of the container for the coefficients for the ansatz functions
- * \param lfsv The local function space object of the test functions
- * \param mat The object containing the local jacobian matrix
- */
- template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
- void jacobian_volume (const EG& eg,
- const LFSU& lfsu,
- const X& x,
- const LFSV& lfsv,
- Dune::PDELab::LocalMatrix<R>& mat) const
- {
- typedef typename LFSU::Traits::SizeType size_type;
-
- model_.localJacobian().assemble(eg.entity());
-
- int numVertices = x.size()/numEq;
- for (size_type j=0; j<lfsu.size(); j++) {
- for (size_type i=0; i<lfsu.size(); i++) {
- mat(i,j) = (model_.localJacobian().mat(i%numVertices,j%numVertices))[i/numVertices][j/numVertices];
- }
- }
- }
-
-private:
- Model& model_;
-};
-
-} // namespace PDELab
-} // namespace Dumux
-
-#endif
diff --git a/dumux/common/pdelabpreconditioner.hh b/dumux/common/pdelabpreconditioner.hh
deleted file mode 100644
index 8297d7b492..0000000000
--- a/dumux/common/pdelabpreconditioner.hh
+++ /dev/null
@@ -1,522 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief A parallel solver for nonoverlapping grid partitions.
- */
-#ifndef DUMUX_PDELAB_PRECONDITIONER_HH
-#define DUMUX_PDELAB_PRECONDITIONER_HH
-
-#include <dumux/common/pardiso.hh>
-#include <dumux/common/propertysystem.hh>
-
-#include <dune/pdelab/backend/istlsolverbackend.hh>
-
-namespace Dumux {
-// forward declaration of property tags
-namespace Properties {
-NEW_PROP_TAG(Problem);
-NEW_PROP_TAG(Model);
-NEW_PROP_TAG(GridView);
-NEW_PROP_TAG(NumEq);
-NEW_PROP_TAG(Grid);
-NEW_PROP_TAG(Scalar);
-NEW_PROP_TAG(VertexMapper);
-NEW_PROP_TAG(GridOperatorSpace);
-NEW_PROP_TAG(JacobianMatrix);
-NEW_PROP_TAG(ReferenceElements);
-NEW_PROP_TAG(GridFunctionSpace);
-NEW_PROP_TAG(ConstraintsTrafo);
-};
-
-
-namespace PDELab {
-
-/*!
- * \brief exchanges matrix entries for parallel computations
- */
-template<class TypeTag>
-class Exchanger
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- 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(VertexMapper)) VertexMapper;
- 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 GridView::ctype CoordScalar;
-
- typedef typename Dune::GenericReferenceElements<CoordScalar, dim> ReferenceElements;
- typedef typename Dune::GenericReferenceElement<CoordScalar, dim> ReferenceElement;
-
-public:
-
- Exchanger(const Problem& problem)
- : gridView_(problem.gridView()), vertexMapper_(problem.vertexMapper()), borderIndices_(0)
- {
- gid2Index_.clear();
- index2GID_.clear();
-
- VertexIterator vertexEndIt = gridView_.template end<dim>();
- for (VertexIterator vertexIt = gridView_.template begin<dim>(); vertexIt != vertexEndIt; ++vertexIt)
- {
-// std::cout << gridView_.comm().rank() << ": node " << vertexMapper_.map(*vertexIt)
-// << " at (" << vertexIt->geometry().corner(0) << ") is of type "
-// << vertexIt->partitionType() << ", GID = "
-// << gridView_.grid().globalIdSet().id(*vertexIt) << std::endl;
- if (vertexIt->partitionType() == Dune::BorderEntity)
- {
- int localIdx = vertexMapper_.map(*vertexIt);
- IdType globalIdx = gridView_.grid().globalIdSet().id(*vertexIt);
-
- std::pair<IdType,int> g2iPair(globalIdx, localIdx);
- gid2Index_.insert(g2iPair);
-
- std::pair<int,IdType> i2gPair(localIdx, globalIdx);
- index2GID_.insert(i2gPair);
-
- borderIndices_.push_back(localIdx);
- }
- }
- }
-
- /*!
- * \brief matrix entry for the MatEntryExchange class
- */
- struct MatEntry
- {
- IdType first;
- BlockType second;
- MatEntry (const IdType& f, const BlockType& s) : first(f),second(s) {}
- MatEntry () {}
- };
-
- /*!
- * \brief A DataHandle class to exchange matrix entries
- */
- class MatEntryExchange
- : 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
- {
- return (codim==dim);
- }
-
- //! returns true if size per entity of given dim and codim is a constant
- bool fixedsize (int dim, int codim) const
- {
- return false;
- }
-
- /*! how many objects of type DataType have to be sent for a given entity
-
- Note: Only the sender side needs to know this size.
- */
- template<class EntityType>
- size_t size (EntityType& e) const
- {
-// std::cout << gridView_.comm().rank() << ": begin loop over vertices.\n";
-// VertexIterator vertexEndIt = gridView_.template end<dim>();
-// for (VertexIterator vertexIt = gridView_.template begin<dim>(); vertexIt != vertexEndIt; ++vertexIt)
-// {
-// std::cout << gridView_.comm().rank() << ": node " << vertexMapper_.map(*vertexIt)
-// << " at (" << vertexIt->geometry().corner(0) << ") is of type "
-// << vertexIt->partitionType() << ", GID = "
-// << gridView_.grid().globalIdSet().id(*vertexIt) << std::endl;
-// }
-// std::cout << gridView_.comm().rank() << ": end loop over vertices.\n";
-// std::cout.flush();
-// std::cout << gridView_.comm().rank() << ": node " << vertexMapper_.map(e)
-// << " on level " << e.level() << " at (" << e.geometry().corner(0) << ") is of type "
-// << e.partitionType() << ", GID = "
-// << gridView_.grid().globalIdSet().id(e) << std::endl;;
- int i = vertexMapper_.map(e);
- int n = 0;
- for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
- {
- // only count those entries corresponding to border entities
- typename std::map<int,IdType>::const_iterator it = index2GID_.find(j.index());
- if (it != index2GID_.end())
- n++;
- }
-// std::cout << gridView_.comm().rank() << ": node " << i << " has sending size " << n << std::endl;
- return n;
- }
-
- //! pack data from user to message buffer
- template<class MessageBuffer, class EntityType>
- void gather (MessageBuffer& buff, const EntityType& e) const
- {
- int i = vertexMapper_.map(e);
- for (ColIterator j = A_[i].begin(); j != A_[i].end(); ++j)
- {
- // only send those entries corresponding to border entities
- typename std::map<int,IdType>::const_iterator it=index2GID_.find(j.index());
- if (it != index2GID_.end()) {
- buff.write(MatEntry(it->second,*j));
-// std::cout << gridView_.comm().rank() << ": node " << i << " gathers (" << it->second
-// << ", " << *j << ") for j = " << j.index() << std::endl;
- }
- }
- }
-
- /*! unpack data from message buffer to user
-
- n is the number of objects sent by the sender
- */
- template<class MessageBuffer, class EntityType>
- void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
- {
- int i = vertexMapper_.map(e);
- for (size_t k = 0; k < n; k++)
- {
- MatEntry m;
- buff.read(m);
- // only add entries corresponding to border entities
- typename std::map<IdType,int>::const_iterator it = gid2Index_.find(m.first);
- if (it != gid2Index_.end())
- {
-// std::cout << gridView_.comm().rank() << ": node " << i << " adds " << m.second
-// << " to j = " << it->second << ", GID = " << m.first << std::endl;
- if (A_[i].find(it->second) != A_[i].end())
- A_[i][it->second] += m.second;
- }
- }
- }
-
- //! constructor
- MatEntryExchange (const GridView& gridView, const std::map<IdType,int>& g2i,
- const std::map<int,IdType>& i2g,
- const VertexMapper& vm,
- JacobianMatrix& A)
- : gridView_(gridView), gid2Index_(g2i), index2GID_(i2g), vertexMapper_(vm), A_(A)
- {}
-
-private:
- const GridView& gridView_;
- const std::map<IdType,int>& gid2Index_;
- const std::map<int,IdType>& index2GID_;
- const VertexMapper& vertexMapper_;
- JacobianMatrix& A_;
- };
-
- void sumEntries (JacobianMatrix& A)
- {
- if (gridView_.comm().size() > 1)
- {
- MatEntryExchange datahandle(gridView_, gid2Index_, index2GID_, vertexMapper_, A);
- gridView_.communicate(datahandle,
- Dune::InteriorBorder_InteriorBorder_Interface,
- Dune::ForwardCommunication);
- }
- }
-
- const std::vector<int>& borderIndices()
- {
- return borderIndices_;
- }
-
-private:
- const GridView& gridView_;
- std::map<IdType,int> gid2Index_;
- std::map<int,IdType> index2GID_;
- const VertexMapper& vertexMapper_;
- std::vector<int> borderIndices_;
-};
-
-/*!
- * \brief Wrapper for a sequential preconditioner
- */
-template<class CC, class GFS, class P>
-class NonoverlappingWrappedPreconditioner
- : public Dune::Preconditioner<typename P::domain_type,typename P::range_type>
-{
-public:
- //! \brief The domain type of the preconditioner.
- typedef typename P::domain_type domain_type;
- //! \brief The range type of the preconditioner.
- typedef typename P::range_type range_type;
-
- // define the category
- enum {
- //! \brief The category the preconditioner is part of.
- category=Dune::SolverCategory::nonoverlapping
- };
-
- //! Constructor.
- NonoverlappingWrappedPreconditioner (const GFS& gfs_, P& prec_, const CC& cc_,
- const std::vector<int>& borderIndices,
- const Dune::PDELab::ParallelISTLHelper<GFS>& helper_)
- : gfs(gfs_), prec(prec_), cc(cc_), borderIndices_(borderIndices), helper(helper_)
- {}
-
- /*!
- \brief \copybrief Dune::SeqPardiso::pre
-
- \copydetails Dune::SeqPardiso::pre
- */
- virtual void pre (domain_type& x, range_type& b)
- {
- prec.pre(x,b);
- }
-
- /*!
- \brief \copybrief Dune::SeqPardiso::apply(X&,const Y&)
-
- \copydetails Dune::SeqPardiso::apply(X&,const Y&)
- */
- virtual void apply (domain_type& v, const range_type& d)
- {
- range_type dd(d);
- set_constrained_dofs(cc,0.0,dd);
- prec.apply(v,dd);
-
- Dune::PDELab::AddDataHandle<GFS,domain_type> adddh(gfs,v);
- if (gfs.gridview().comm().size()>1)
- gfs.gridview().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
-
- for (int k = 0; k < borderIndices_.size(); k++)
- v[borderIndices_[k]] *= 0.5;
- }
-
- /*!
- \brief \copybrief Dune::SeqPardiso::post(X&)
-
- \copydetails Dune::SeqPardiso::post(X&)
- */
- virtual void post (domain_type& x)
- {
- prec.post(x);
- }
-
-private:
- const GFS& gfs;
- P& prec;
- const CC& cc;
- const std::vector<int>& borderIndices_;
- const Dune::PDELab::ParallelISTLHelper<GFS>& helper;
-};
-
-/*!
- * \brief backend for an ISTL parallel ILU preconditioned BiCGSTAB solver
- */
-template<class TypeTag>
-class ISTLBackend_NoOverlap_BCGS_ILU
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- 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:
- /*! \brief make a linear solver object
-
- \param[in] problem The Dumux problem
- \param[in] maxiter_ Maximum number of iterations to do
- \param[in] verbose_ Verbosity level
- */
- 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
-
- \param[in] v the given vector
- */
- template<class Vector>
- typename Vector::ElementType norm (const Vector& v) const
- {
- Vector x(v); // make a copy because it has to be made consistent
- typedef Dune::PDELab::NonoverlappingScalarProduct<GridFunctionSpace,Vector> PSP;
- PSP psp(gfs,phelper);
- psp.make_consistent(x);
- return psp.norm(x);
- }
-
- /*! \brief solve the given linear system
-
- \param[in] A the given matrix
- \param[out] z the solution vector to be computed
- \param[in] r right hand side
- \param[in] reduction to be achieved
- */
- template<class JacobianMatrix, class SolVector, class RhsVector>
- void apply(const JacobianMatrix& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
- {
- typedef Dune::SeqILU0<JacobianMatrix,SolVector,RhsVector> SeqPreCond;
- JacobianMatrix B(A);
- exchanger_.sumEntries(B);
- SeqPreCond seqPreCond(B, 1.0);
-
- 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 NonoverlappingWrappedPreconditioner<ConstraintsTrafo, GridFunctionSpace, SeqPreCond> ParPreCond;
- ParPreCond parPreCond(gfs, seqPreCond, constraintsTrafo_, exchanger_.borderIndices(), phelper);
-
- int verb=0;
- if (gfs.gridview().comm().rank()==0) verb=verbose;
- Dune::BiCGSTABSolver<SolVector> solver(pop,psp,parPreCond,reduction,maxiter,verb);
- Dune::InverseOperatorResult stat;
- solver.apply(z,r,stat);
- res.converged = stat.converged;
- res.iterations = stat.iterations;
- res.elapsed = stat.elapsed;
- res.reduction = stat.reduction;
- }
-
- /*! \brief Return access to result data */
- const Dune::PDELab::LinearSolverResult<double>& result() const
- {
- return res;
- }
-
-private:
- const GridFunctionSpace& gfs;
- PHELPER phelper;
- Dune::PDELab::LinearSolverResult<double> res;
- unsigned maxiter;
- int verbose;
- const ConstraintsTrafo& constraintsTrafo_;
- Exchanger<TypeTag> exchanger_;
-};
-
-/*!
- * \brief backend for an ISTL parallel Pardiso preconditioned loop solver
- */
-template<class TypeTag>
-class ISTLBackend_NoOverlap_Loop_Pardiso
-{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- 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(JacobianMatrix)) JacobianMatrix;
- typedef Dune::PDELab::ParallelISTLHelper<GridFunctionSpace> PHELPER;
-
-public:
- /*! \brief make a linear solver object
-
- \param[in] problem The Dumux problem
- \param[in] maxiter_ Maximum number of iterations to do
- \param[in] verbose_ Verbosity level
- */
- 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
-
- \param[in] v the given vector
- */
- template<class Vector>
- typename Vector::ElementType norm (const Vector& v) const
- {
- Vector x(v); // make a copy because it has to be made consistent
- typedef Dune::PDELab::NonoverlappingScalarProduct<GridFunctionSpace,Vector> PSP;
- PSP psp(gfs,phelper);
- psp.make_consistent(x);
- return psp.norm(x);
- }
-
- /*! \brief solve the given linear system
-
- \param[in] A the given matrix
- \param[out] z the solution vector to be computed
- \param[in] r right hand side
- \param[in] reduction to be achieved
- */
- template<class M, class SolVector, class RhsVector>
- void apply(M& A, SolVector& z, RhsVector& r, typename SolVector::ElementType reduction)
- {
- typedef Dune::SeqPardiso<JacobianMatrix,SolVector,RhsVector> SeqPreCond;
- JacobianMatrix B(A);
- exchanger_.sumEntries(B);
- SeqPreCond seqPreCond(B);
-
- 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 NonoverlappingWrappedPreconditioner<ConstraintsTrafo, GridFunctionSpace, SeqPreCond> ParPreCond;
- ParPreCond parPreCond(gfs, seqPreCond, constraintsTrafo_, exchanger_.borderIndices(), phelper);
-
- // typedef Dune::PDELab::NonoverlappingRichardson<GridFunctionSpace,SolVector,RhsVector> PRICH;
- // PRICH prich(gfs,phelper);
- int verb=0;
- if (gfs.gridview().comm().rank()==0) verb=verbose;
- Dune::LoopSolver<SolVector> solver(pop,psp,parPreCond,reduction,maxiter,verb);
-// Dune::BiCGSTABSolver<SolVector> solver(pop,psp,parPreCond,reduction,maxiter,verb);
- Dune::InverseOperatorResult stat;
- solver.apply(z,r,stat);
- res.converged = stat.converged;
- res.iterations = stat.iterations;
- res.elapsed = stat.elapsed;
- res.reduction = stat.reduction;
- }
-
- /*! \brief Return access to result data */
- const Dune::PDELab::LinearSolverResult<double>& result() const
- {
- return res;
- }
-
-private:
- const GridFunctionSpace& gfs;
- PHELPER phelper;
- Dune::PDELab::LinearSolverResult<double> res;
- unsigned maxiter;
- int verbose;
- const ConstraintsTrafo& constraintsTrafo_;
- Exchanger<TypeTag> exchanger_;
-};
-
-
-
-} // namespace PDELab
-} // namespace Dumux
-
-#endif
diff --git a/dumux/material/binarycoefficients/h2_n2.hh b/dumux/material/binarycoefficients/h2_n2.hh
deleted file mode 100644
index a06a897f0a..0000000000
--- a/dumux/material/binarycoefficients/h2_n2.hh
+++ /dev/null
@@ -1,99 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009 by Andreas Lauser *
- * Institute of Hydraulic Engineering *
- * University of Stuttgart, Germany *
- * email: <givenname>.<name>@iws.uni-stuttgart.de *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Binary coefficients for hydrogen and nitrogen.
- */
-#ifndef DUMUX_BINARY_COEFF_H2_N2_HH
-#define DUMUX_BINARY_COEFF_H2_N2_HH
-
-#include "henryiapws.hh"
-#include "fullermethod.hh"
-
-#include <dumux/material/components/h2.hh>
-#include <dumux/material/components/n2.hh>
-
-namespace Dumux
-{
-namespace BinaryCoeff
-{
-
-/*!
- * \ingroup Binarycoefficients
- * \brief Binary coefficients for hydrogen and nitrogen.
- */
-class H2_N2
-{
-public:
- /*!
- * \brief Henry coefficent \f$\mathrm{[N/m^2]}\f$ for molecular nitrogen in liquid hydrogen.
- *
- * \param temperature the temperature \f$\mathrm{[K]}\f$
- */
- template <class Scalar>
- static Scalar henry(Scalar temperature)
- {
- DUNE_THROW(Dune::NotImplemented, "henry coefficient for nitrogen in liquid hydrogen");
- };
-
- /*!
- * \brief Binary diffusion coefficent \f$\mathrm{[m^2/s]}\f$ for molecular hydrogen and nitrogen.
- *
- * This function estimates the diffusion coefficents in binary gases
- * using to the method proposed by Fuller. This method and is only
- * valid at "low" pressures.
- *
- * See: R. Reid, et al.: The Properties of Gases and Liquids, 4th
- * edition, McGraw-Hill, 1987, pp. 587-588
- * \param temperature the temperature \f$\mathrm{[K]}\f$
- * \param pressure the phase pressure \f$\mathrm{[Pa]}\f$
- */
- template <class Scalar>
- static Scalar gasDiffCoeff(Scalar temperature, Scalar pressure)
- {
- typedef Dumux::H2<Scalar> H2;
- typedef Dumux::N2<Scalar> N2;
-
- // atomic diffusion volumes
- const Scalar SigmaNu[2] = { 6.12 /* H2 */, 18.5 /* N2 */ };
- // molar masses [g/mol]
- const Scalar M[2] = { H2::molarMass()*1e3, N2::molarMass()*1e3 };
-
- return fullerMethod(M, SigmaNu, temperature, pressure);
- };
-
- /*!
- * \brief Diffusion coefficent \f$\mathrm{[m^2/s]}\f$ for molecular nitrogen in liquid hydrogen.
- *
- * \param temperature the temperature \f$\mathrm{[K]}\f$
- * \param pressure the phase pressure \f$\mathrm{[Pa]}\f$
- */
- template <class Scalar>
- static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
- {
- DUNE_THROW(Dune::NotImplemented, "diffusion coefficient for liquid nitrogen and hydrogen");
- };
-};
-
-}
-} // end namepace
-
-#endif
diff --git a/dumux/material/fluidmatrixinteractions/Mp/2poftadapter.hh b/dumux/material/fluidmatrixinteractions/Mp/2poftadapter.hh
deleted file mode 100644
index 035b6dfb65..0000000000
--- a/dumux/material/fluidmatrixinteractions/Mp/2poftadapter.hh
+++ /dev/null
@@ -1,83 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2010 by Philipp Nuske *
- * 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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file 2poftadapter.hh
- *
- * Makes the twophase capillary pressure-saturation relations
- * available under the M-phase API for material laws.
- *
- * Also use the temperature dependent version of the material laws.
- */
-#ifndef DUMUX_MP_2P_OFT_ADAPTER_HH
-#define DUMUX_MP_2P_OFT_ADAPTER_HH
-
-#include <algorithm>
-
-namespace Dumux
-{
-/*!
- * \ingroup material
- *
- * \brief Adapts the interface of the MpNc material law to the standard-Dumux material law.
- *
- * Also use the temperature dependent version of the material laws.
- */
-template <int wPhaseIdx, class TwoPLaw>
-class TwoPOfTAdapter
-{
- enum { nPhaseIdx = (wPhaseIdx == 0)?1:0 };
-
-public:
- typedef typename TwoPLaw::Params Params;
- typedef typename Params::Scalar Scalar;
- enum { numPhases = 2 };
-
- /*!
- * \brief The capillary pressure-saturation curve.
- */
- template <class pcContainerT, class FluidState>
- static void capillaryPressures(pcContainerT &pc,
- const Params ¶ms,
- const FluidState &fluidState)
- {
- // non-wetting phase gets the capillary pressure added
- pc[nPhaseIdx] = 0;
-
- // wetting phase does not get anything added
- pc[wPhaseIdx] = - TwoPLaw::pC(params, fluidState.saturation(wPhaseIdx), fluidState.temperature(wPhaseIdx));
- }
-
-
-
- /*!
- * \brief The relative permeability of all phases.
- */
- template <class krContainerT, class FluidState>
- static void relativePermeabilities(krContainerT &kr,
- const Params ¶ms,
- const FluidState &fluidState)
- {
- kr[wPhaseIdx] = TwoPLaw::krw(params, fluidState.saturation(wPhaseIdx));
- kr[nPhaseIdx] = TwoPLaw::krn(params, fluidState.saturation(wPhaseIdx));
- };
-};
-}
-
-#endif
diff --git a/dumux/material/fluidsystems/simplefluidsystem.hh b/dumux/material/fluidsystems/simplefluidsystem.hh
deleted file mode 100644
index 09310f6df9..0000000000
--- a/dumux/material/fluidsystems/simplefluidsystem.hh
+++ /dev/null
@@ -1,326 +0,0 @@
-/*****************************************************************************
- * 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. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Base class for all fluid systems which do not require
- * mutable parameters.
- */
-#ifndef DUMUX_SIMPLE_FLUID_SYSTEM_HH
-#define DUMUX_SIMPLE_FLUID_SYSTEM_HH
-
-#include <dumux/material/fluidstates/genericfluidstate.hh>
-
-#include <dumux/material/components/simpleh2o.hh>
-#include <dumux/material/components/h2o.hh>
-#include <dumux/material/components/n2.hh>
-#include <dumux/material/idealgas.hh>
-
-#include <dumux/material/binarycoefficients/h2o_n2.hh>
-
-#include <dumux/common/exceptions.hh>
-
-#include "basicmutableparameters.hh"
-
-namespace Dumux
-{
-/*!
- * \brief Base class for all fluid systems which do not require
- * mutable parameters.
- */
-template <class Scalar, class StaticParametersT, class Implementation>
-class SimpleFluidSystem : public StaticParametersT
-{
-public:
- typedef StaticParametersT StaticParameters;
- typedef SimpleMutableParameters<Scalar, StaticParameters> MutableParameters;
- typedef typename MutableParameters::FluidState FluidState;
-
-public:
- /*!
- * \brief Initialize the fluid system's static parameters
- */
- static void init()
- {}
-
- static Scalar computeMolarVolume(const MutableParameters &mutParams, int phaseIdx)
- { Implementation::computeMolarVolume_(mutParams.fluidState(), phaseIdx); }
-
- /*!
- * \brief Calculate the molar volume [m^3/mol] of a fluid phase.
- *
- * This is the method which does not require any mutable
- * parameters and can thus only gets a fluid state object as
- * argument.
- */
- static Scalar computeMolarVolume_(const FluidState &fs, int phaseIdx)
- {
- DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
- };
-
- static Scalar computePureFugacity(const MutableParameters ¶ms,
- int phaseIdx,
- int compIdx)
- {
- const FluidState &fs = params.fluidState();
- Scalar T = fs.temperature(phaseIdx);
- switch (phaseIdx) {
- case SP::lPhaseIdx:
- switch (compIdx) {
- case SP::H2OIdx: return SP::H2O::vaporPressure(T);
- case SP::N2Idx: return BinaryCoeff::H2O_N2::henry(T);
- };
- case SP::gPhaseIdx: {
- switch (compIdx) {
- case SP::H2OIdx:
- SP::H2O::vaporPressure(T);
- case SP::N2Idx:
- // ideal gas
- return fs.pressure(phaseIdx);
- };
- };
-
- default: DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
- }
- };
-
- /*!
- * \brief Calculate the fugacity of a component in a fluid phase
- * [Pa]
- *
- * The components chemical \f$mu_\kappa\f$ potential is connected
- * to the component's fugacity \f$f_\kappa\f$ by the relation
- *
- * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f]
- *
- * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
- * pressure and temperature.
- */
- static Scalar computeFugacity(const MutableParameters ¶ms,
- int phaseIdx,
- int compIdx)
- {
- const FluidState &fs = params.fluidState();
-
- Scalar x = fs.moleFrac(phaseIdx,compIdx);
- Scalar p = fs.pressure(phaseIdx);
- return x*p*computeFugacityCoeff(params, phaseIdx, compIdx);
- };
-
- /*!
- * \brief Calculate the fugacity coefficient [Pa] of an individual
- * component in a fluid phase
- *
- * The fugacity coefficient \f$\phi_\kappa\f$ is connected to the
- * fugacity \f$f_\kappa\f$ and the component's molarity
- * \f$x_\kappa\f$ by means of the relation
- *
- * \f[ f_\kappa = \phi_\kappa * x_{\kappa} \f]
- */
- static Scalar computeFugacityCoeff(const MutableParameters ¶ms,
- int phaseIdx,
- int compIdx)
- {
- assert(0 <= phaseIdx && phaseIdx <= SP::numPhases);
- assert(0 <= compIdx && compIdx <= SP::numComponents);
-
- const FluidState &fs = params.fluidState();
-
- Scalar T = fs.temperature(phaseIdx);
- Scalar p = fs.pressure(phaseIdx);
- switch (phaseIdx) {
- case SP::lPhaseIdx:
- switch (compIdx) {
- case SP::H2OIdx: return SP::H2O::vaporPressure(T)/p;
- case SP::N2Idx: return BinaryCoeff::H2O_N2::henry(T)/p;
- };
- case SP::gPhaseIdx:
- return 1.0; // ideal gas
- };
-
- DUNE_THROW(Dune::InvalidStateException, "Unhandled phase or component index");
- }
-
- /*!
- * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
- */
- static Scalar computeViscosity(const MutableParameters ¶ms,
- int phaseIdx)
- {
- assert(0 <= phaseIdx && phaseIdx <= SP::numPhases);
-
- const FluidState &fs = params.fluidState();
- Scalar T = fs.temperature(phaseIdx);
- Scalar p = fs.pressure(phaseIdx);
- switch (phaseIdx) {
- case SP::lPhaseIdx:
- // assume pure water for the liquid phase
- return SP::H2O::liquidViscosity(T, p);
- case SP::gPhaseIdx:
- // assume pure water for the gas phase
- return SP::N2::gasViscosity(T, p);
- }
-
- DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
- };
-
- /*!
- * \brief Calculate the binary molecular diffusion coefficient for
- * a component in a fluid phase [mol^2 * s / (kg*m^3)]
- *
- * Molecular diffusion of a compoent \f$\kappa\f$ is caused by a
- * gradient of the chemical potential and follows the law
- *
- * \f[ J = - D \grad mu_\kappa \f]
- *
- * where \f$\mu_\kappa\f$ is the component's chemical potential,
- * \f$D\f$ is the diffusion coefficient and \f$J\f$ is the
- * diffusive flux. \f$mu_\kappa\f$ is connected to the component's
- * fugacity \f$f_\kappa\f$ by the relation
- *
- * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f]
- *
- * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
- * pressure and temperature.
- */
- static Scalar computeDiffusionCoeff(const MutableParameters ¶ms,
- int phaseIdx,
- int compIdx)
- {
- // TODO!
- DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
- };
-
- /*!
- * \brief Given a phase's composition, temperature and pressure,
- * return the binary diffusion coefficient for components
- * \f$i\f$ and \f$j\f$ in this phase.
- */
- template <class FluidState>
- static Scalar binaryDiffCoeff(const MutableParameters ¶ms,
- int phaseIdx,
- int compIIdx,
- int compJIdx)
-
- {
- if (compIIdx > compJIdx)
- std::swap(compIIdx, compJIdx);
-
-#ifndef NDEBUG
- if (compIIdx == compJIdx ||
- phaseIdx > SP::numPhases - 1 ||
- compJIdx > SP::numComponents - 1)
- {
- DUNE_THROW(Dune::InvalidStateException,
- "Binary diffusion coefficient of components "
- << compIIdx << " and " << compJIdx
- << " in phase " << phaseIdx << " is undefined!\n");
- }
-#endif
-
- const FluidState &fs = params.fluidState();
- Scalar T = fs.temperature(phaseIdx);
- Scalar p = fs.pressure(phaseIdx);
-
- switch (phaseIdx) {
- case SP::lPhaseIdx:
- switch (compIIdx) {
- case SP::H2OIdx:
- switch (compJIdx) {
- case SP::N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
- }
- default:
- DUNE_THROW(Dune::InvalidStateException,
- "Binary diffusion coefficients of trace "
- "substances in liquid phase is undefined!\n");
- }
- case SP::gPhaseIdx:
- switch (compIIdx) {
- case SP::H2OIdx:
- switch (compJIdx) {
- case SP::N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(T, p);
- }
- }
- }
-
- DUNE_THROW(Dune::InvalidStateException,
- "Binary diffusion coefficient of components "
- << compIIdx << " and " << compJIdx
- << " in phase " << phaseIdx << " is undefined!\n");
- };
-
- /*!
- * \brief Given a phase's composition, temperature, pressure and
- * density, calculate its specific enthalpy [J/kg].
- */
- /*!
- * \todo This system neglects the contribution of gas-molecules in the liquid phase.
- * This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ...
- */
- static Scalar computeEnthalpy(const MutableParameters ¶ms,
- int phaseIdx)
- {
- const FluidState &fs = params.fluidState();
- Scalar T = fs.temperature(phaseIdx);
- Scalar p = fs.pressure(phaseIdx);
-
- if (phaseIdx == SP::lPhaseIdx) {
-#warning hack
- T = 300.0;
- // TODO: correct way to deal with the solutes???
- return
- SP::H2O::liquidEnthalpy(T, p);
- }
- else {
- // assume ideal gas
- Scalar XH2O = fs.massFrac(SP::lPhaseIdx, SP::H2OIdx);
- Scalar XN2 = fs.massFrac(SP::lPhaseIdx, SP::N2Idx);
- Scalar result = 0;
- result +=
- SP::H2O::gasEnthalpy(T, p) *
- fs.massFrac(SP::gPhaseIdx, SP::H2OIdx);
- result +=
- SP::N2::gasEnthalpy(T, p) *
- fs.massFrac(SP::gPhaseIdx, SP::N2Idx);
- return result;
- }
- };
-
- /*!
- * \brief Given a phase's composition, temperature, pressure and
- * density, calculate its specific internal energy [J/kg].
- */
- static Scalar internalEnergy(const MutableParameters ¶ms,
- int phaseIdx)
- {
- const FluidState &fs = params.fluidState();
- Scalar T = fs.temperature(phaseIdx);
- Scalar p = fs.pressure(phaseIdx);
- Scalar rho = fs.density(phaseIdx);
- return
- computeEnthalpy(params, phaseIdx) -
- p/rho;
-
- };
-};
-
-} // end namepace
-
-#endif
--
GitLab