Commit c80a076a authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

test/boxmodels/1p2c works in stable

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@3839 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 0f70a035
// $Id: 1p2cboxproblem.hh 3783 2010-06-24 11:33:53Z bernd $
// $Id: 1p2cboxproblem.hh 3838 2010-07-15 08:31:53Z bernd $
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
......@@ -23,7 +23,6 @@
#define DUMUX_1P2C_BOX_PROBLEM_HH
#include <dumux/boxmodels/boxscheme/boxproblem.hh>
#include <dumux/old_material/twophaserelations.hh>
namespace Dumux
{
......@@ -43,8 +42,7 @@ class OnePTwoCBoxProblem : public BoxProblem<TypeTag, Implementation>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
// material properties
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Fluid)) Fluid;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Soil)) Soil;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
enum {
dim = Grid::dimension,
......@@ -56,7 +54,8 @@ class OnePTwoCBoxProblem : public BoxProblem<TypeTag, Implementation>
public:
OnePTwoCBoxProblem(const GridView &gridView)
: ParentType(gridView),
gravity_(0)
gravity_(0),
spatialParams_(gridView)
{
if (GET_PROP_VALUE(TypeTag, PTAG(EnableGravity)))
gravity_[dim-1] = -9.81;
......@@ -84,23 +83,17 @@ public:
const GlobalPosition &gravity() const
{ return gravity_; }
/*!
* \brief Fluid properties of the wetting phase.
*/
const Fluid &fluid() const
{ return fluid_; }
/*!
* \brief Returns the soil properties object.
*/
Soil &soil()
{ return soil_; }
SpatialParameters &spatialParameters()
{ return spatialParams_; }
/*!
* \copydoc soil()
*/
const Soil &soil() const
{ return soil_; }
const SpatialParameters &spatialParameters() const
{ return spatialParams_; }
// \}
......@@ -113,11 +106,10 @@ private:
const Implementation *asImp_() const
{ return static_cast<const Implementation *>(this); }
// fluids and material properties
Fluid fluid_;
Soil soil_;
GlobalPosition gravity_;
// soil properties
SpatialParameters spatialParams_;
};
}
......
// $Id: 1p2cfluidstate.hh 3784 2010-06-24 13:43:57Z bernd $
/*****************************************************************************
* Copyright (C) 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 Calcultes the phase state from the primary variables in the
* 1p2c model.
*/
#ifndef DUMUX_1P2C_PHASE_STATE_HH
#define DUMUX_1P2C_PHASE_STATE_HH
#include "1p2cproperties.hh"
#include <dumux/material/fluidstate.hh>
namespace Dumux
{
/*!
* \brief Calcultes the phase state from the primary variables in the
* 1p2c model.
*/
template <class TypeTag>
class OnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)),
OnePTwoCFluidState<TypeTag> >
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
enum {
konti = Indices::konti,
transport = Indices::transport
};
public:
enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
enum { numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)) };
/*!
* \brief Update the phase state from the primary variables.
*/
void update(const PrimaryVarVector &primaryVars,
Scalar temperature)
{
Valgrind::CheckDefined(primaryVars);
temperature_ = temperature;
phasePressure_[0] = primaryVars[konti];
moleFraction_[0][konti] = 1.0 - primaryVars[transport];
moleFraction_[0][transport] = primaryVars[transport];
}
public:
/*!
* \brief Returns the molar fraction of a component in a fluid phase.
*/
Scalar moleFrac(int phaseIdx, int compIdx) const
{
return moleFraction_[phaseIdx][compIdx];
}
/*!
* \brief Returns the pressure of a fluid phase [Pa].
*/
Scalar phasePressure(int phaseIdx) const
{ return phasePressure_[phaseIdx]; }
/*!
* \brief Returns the temperature of the fluids [K].
*
* Note that we assume thermodynamic equilibrium, so all fluids
* and the rock matrix exhibit the same temperature.
*/
Scalar temperature() const
{ return temperature_; };
public:
Scalar moleFraction_[numPhases][numComponents];
Scalar phasePressure_[numPhases];
Scalar temperature_;
};
} // end namepace
#endif
// $Id: 1p2cfluxdata.hh 3357 2010-03-25 13:02:05Z lauser $
// $Id: 1p2cfluxdata.hh 3838 2010-07-15 08:31:53Z bernd $
/*****************************************************************************
* Copyright (C) 2009 by Karin Erbertseder *
* Copyright (C) 2009 by Andreas Lauser *
......@@ -104,7 +104,7 @@ private:
const VertexDataArray &elemDat)
{
GlobalPosition tmp;
if (!problem.soil().useTwoPointGradient(element, face->i, face->j)) {
if (!problem.spatialParameters().useTwoPointGradient(element, face->i, face->j)) {
// use finite-element gradients
tmp = 0.0;
for (int idx = 0;
......@@ -168,17 +168,14 @@ private:
const Element &element,
const VertexDataArray &elemDat)
{
// calculate the permeability tensor. TODO: this should be
// more flexible
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
const Tensor &Ki = problem.soil().K(insideSCV->global,
element,
insideSCV->local);
const Tensor &Kj = problem.soil().K(outsideSCV->global,
element,
outsideSCV->local);
Tensor K;
Dumux::harmonicMeanMatrix(K, Ki, Kj);
problem.spatialParameters().meanK(K,
problem.spatialParameters().intrinsicPermeability(element,
fvElemGeom,
face->i),
problem.spatialParameters().intrinsicPermeability(element,
fvElemGeom,
face->j));
// vector for the Darcy velocity
K.mv(pressureGrad, vDarcy); // vDarcy = K * grad p
......
// $Id: 1p2cproperties.hh 3736 2010-06-15 09:52:10Z lauser $
// $Id: 1p2cproperties.hh 3838 2010-07-15 08:31:53Z bernd $
/*****************************************************************************
* Copyright (C) 2009 by Karin Erbertseder *
* Copyright (C) 2009 by Andreas Lauser *
......@@ -25,6 +25,8 @@
#ifndef DUMUX_1P2C_PROPERTIES_HH
#define DUMUX_1P2C_PROPERTIES_HH
#include<dumux/boxmodels/boxscheme/boxproperties.hh>
namespace Dumux
{
////////////////////////////////
......@@ -72,8 +74,8 @@ NEW_TYPE_TAG(BoxOnePTwoC, INHERITS_FROM(BoxScheme));
NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
NEW_PROP_TAG(OnePTwoCIndices); //!< Enumerations for the 1p2c models
NEW_PROP_TAG(Soil); //!< The type of the soil properties object
NEW_PROP_TAG(Fluid); //!< The fluid for the single-phase problems
NEW_PROP_TAG(SpatialParameters); //!< The type of the soil
NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
NEW_PROP_TAG(UpwindAlpha); //!< The default value of the upwind parameter
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
......
// $Id: 1p2cvertexdata.hh 3784 2010-06-24 13:43:57Z bernd $
// $Id: 1p2cvertexdata.hh 3838 2010-07-15 08:31:53Z bernd $
/*****************************************************************************
* Copyright (C) 2009 by Karin Erbertseder *
* Copyright (C) 2009 by Andreas Lauser *
......@@ -24,7 +24,7 @@
#ifndef DUMUX_1P2C_VERTEX_DATA_HH
#define DUMUX_1P2C_VERTEX_DATA_HH
#include "1p2cproperties.hh"
#include "1p2cfluidstate.hh"
namespace Dumux
{
......@@ -39,6 +39,8 @@ class OnePTwoCVertexData
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(OnePTwoCIndices)) Indices;
typedef OnePTwoCFluidState<TypeTag> FluidState;
typedef typename GridView::template Codim<0>::Entity Element;
......@@ -49,6 +51,9 @@ class OnePTwoCVertexData
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
konti = Indices::konti,
transport = Indices::transport
};
typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
......@@ -57,9 +62,8 @@ class OnePTwoCVertexData
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
typedef typename SolutionTypes::PrimaryVarVector PrimaryVarVector;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
......@@ -75,25 +79,19 @@ public:
const Problem &problem,
bool isOldSol)
{
typedef Indices I;
// coordinates of the vertex
const GlobalPosition &global = element.geometry().corner(vertIdx);
const LocalPosition &local =
ReferenceElements::general(element.type()).position(vertIdx,
dim);
Scalar T = problem.temperature(element, elemGeom, vertIdx);
pressure = sol[I::konti];
molefraction = sol[I::transport];
porosity = problem.soil().porosity(global,element,local);
density = problem.fluid().density(T, pressure);
diffCoeff = problem.fluid().diffCoeff(T, pressure);
viscosity = problem.fluid().viscosity(T, pressure);
tortuosity = problem.soil().tortuosity(global,element,local);
dispersivity = problem.soil().dispersivity(global,element,local);
molarDensity = problem.fluid().molarDensity(T, pressure);
porosity = problem.spatialParameters().porosity(element, elemGeom, vertIdx);
tortuosity = problem.spatialParameters().tortuosity(element, elemGeom, vertIdx);
dispersivity = problem.spatialParameters().dispersivity(element, elemGeom, vertIdx);
Scalar temperature = problem.temperature(element, elemGeom, vertIdx);
fluidState_.update(sol, temperature);
pressure = fluidState_.phasePressure(konti);
molefraction = fluidState_.moleFrac(konti, transport);
density = FluidSystem::phaseDensity(konti, temperature, pressure, fluidState_);
molarDensity = FluidSystem::molarDensity(konti, temperature, pressure, fluidState_);
viscosity = FluidSystem::phaseViscosity(konti, temperature, pressure, fluidState_);
diffCoeff = FluidSystem::diffCoeff(konti, transport, konti, temperature, pressure, fluidState_);
}
Scalar porosity;
......@@ -105,6 +103,8 @@ public:
Scalar molefraction;
Scalar pressure;
Scalar molarDensity;
FluidState fluidState_;
};
}
......
// $Id: croperator.hh 3777 2010-06-24 06:46:46Z bernd $
// $Id: croperator.hh 3838 2010-07-15 08:31:53Z bernd $
/*****************************************************************************
* Copyright (C) 2007-2009 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
......
// $Id: mimeticgroundwater.hh 3732 2010-06-11 13:27:20Z bernd $
// $Id: mimeticgroundwater.hh 3834 2010-07-14 12:50:32Z bernd $
/*****************************************************************************
* Copyright (C) 2008-2009 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
......
// $Id: mimeticoperator.hh 3777 2010-06-24 06:46:46Z bernd $
// $Id: mimeticoperator.hh 3834 2010-07-14 12:50:32Z bernd $
/*****************************************************************************
* Copyright (C) 2007-2009 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
......
/*****************************************************************************
* Copyright (C) 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 A fluid system with one phase and an arbitrary number of components.
*/
#ifndef DUMUX_ISFLUID_TRAIL_SYSTEM_HH
#define DUMUX_ISFLUID_TRAIL_SYSTEM_HH
#include <dumux/common/propertysystem.hh>
#include <dumux/boxmodels/1p2c/1p2cproperties.hh>
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(Scalar);
};
/*!
* \brief A fluid system with one phase and an arbitrary number of components.
*/
template <class TypeTag, bool verbose=true>
class ISFluid_Trail_System
{
typedef ISFluid_Trail_System<TypeTag, verbose> ThisType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
enum {
isfluid = Indices::konti,
trail = Indices::transport
};
public:
static void init()
{}
/*!
* \brief Return the human readable name of a component
*/
static const char *componentName(int compIdx)
{
switch(compIdx)
{
case isfluid:
return "ISFluid";
case trail:
return "Trail";
default:
DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
}
}
/*!
* \brief Given all mole fractions in a phase, return the phase
* density [kg/m^3].
*/
template <class FluidState>
static Scalar phaseDensity(int phaseIdx,
Scalar temperature,
Scalar pressure,
const FluidState &fluidState)
{
if (phaseIdx == 0)
return 1.03e-9; // in [kg /microm^3]
// Umwandlung in Pascal notwendig -> density*10^6
// will man wieder mit kg/m^3 rechnen muss g auch wieder geändert werden
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief Given all mole fractions in a phase, return the phase
* density [kg/m^3].
*/
template <class FluidState>
static Scalar molarDensity(int phaseIdx,
Scalar temperature,
Scalar pressure,
const FluidState &fluidState)
{
if (phaseIdx == 0)
return 3.035e-16; // in [mol/microm^3] = 303.5 mol/m^3
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief Return the viscosity of a phase.
*/
template <class FluidState>
static Scalar phaseViscosity(int phaseIdx,
Scalar temperature,
Scalar pressure,
const FluidState &fluidState)
{
if (phaseIdx == 0)
return 0.00069152; // in [Pas]
DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
}
/*!
* \brief Given all mole fractions, return the diffusion
* coefficent of a component in a phase.
*/
template <class FluidState>
static Scalar diffCoeff(int phaseIdx,
int compIIdx,
int compJIdx,
Scalar temperature,
Scalar pressure,
const FluidState &fluidState)
{
return 0.088786695; // in [microm^2/s] = 3.7378e-12 m^2/s
}
};
} // end namepace
#endif
// $Id: tissue_soilproperties.hh 3783 2010-06-24 11:33:53Z bernd $
/*****************************************************************************
* Copyright (C) 2009 by Karin Erbertseder *
* 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_TISSUE_SOILPROPERTIES_HH
#define DUMUX_TISSUE_SOILPROPERTIES_HH
#include <dumux/old_material/property_baseclasses.hh>
namespace Dumux
{
/**
* \brief Definition of the properties of the human tissue
*
* \todo do not derive from Matrix2p, since this is the base class for
* two-phase models
*/
template<class Grid, class Scalar>
class TissueSoil: public Matrix2p<Grid, Scalar>
{
public:
typedef typename Grid::Traits::template Codim<0>::Entity Element;
typedef typename Grid::ctype CoordScalar;
enum
{ dim=Grid::dimension, numEq=2};
const Dune::FieldMatrix<CoordScalar,dim,dim> &K (const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos) const
{
if (isTumor_(globalPos))
return permTumor_;
else
return permTissue_;
}
double porosity(const Dune::FieldVector<CoordScalar,dim> &globalPos, const Element &element, const Dune::FieldVector<CoordScalar,dim> &localPos) const
{
if (isTumor_(globalPos))
return 0.31; // tumor tissue
else
return 0.13; // healthy tissue
}
double tortuosity (const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos) const
{
if (isTumor_(globalPos))
return 0.706; // tumor tissue
else
return 0.280; // healthy tissue
}
virtual double Sr_w(const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos, const double T = 283.15) const
{
return 0;
}
virtual double Sr_n(const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos, const double T = 283.15) const
{
return 0;
}
virtual std::vector<double> paramRelPerm(const Dune::FieldVector<CoordScalar,dim>& globalPos, const Element& element, const Dune::FieldVector<CoordScalar,dim>& localPos, const double T = 283.15) const
{
std::vector<double> param(0);
return param;
}
bool useTwoPointGradient(const Element &elem,
int vertexI,
int vertexJ) const
{
bool inTumor = isTumor_(elem.geometry().corner(0));
int n = elem.template count<dim>();
for (int i = 1; i < n; ++i) {
if ((inTumor && !isTumor_(elem.geometry().corner(i))) ||
(!inTumor && isTumor_(elem.geometry().corner(i))))
return true;
}
return false;
}
TissueSoil()
:Matrix2p<Grid,Scalar>()
{
permTumor_ = 0;
permTissue_ = 0;
for (int i = 0; i < dim; i++)
permTumor_[i][i] = 2.142e-11; //tumor tissue
for (int i = 0; i < dim; i++)
permTissue_[i][i] = 4.424e-12; //healthy tissue
}
private:
bool isTumor_(const Dune::FieldVector<CoordScalar,dim>& globalPos) const
{
if(globalPos[0] > 0.99 && globalPos[0] < 2.01 &&
globalPos[1] > 0.99 && globalPos[1] < 3.01)
return true;