Commit fad37139 authored by Benjamin Faigle's avatar Benjamin Faigle
Browse files

multiphysics on new decoupled framework and new data structure

consequently methods from standard model

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7372 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent eb277902
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2007-2008 by Jochen Fritz *
* 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 boundary condition flags for the decoupled 2p2c model
*/
#ifndef DUMUX_BOUNDARYCONDITIONS2P2C_HH
#define DUMUX_BOUNDARYCONDITIONS2P2C_HH
namespace Dumux
{
/**
* \ingroup IMPEC IMPETbc
* \brief Defines type of boundary conditions for 2p2c processes
*
* This is to distinguish BC types for 2p2c processes similar to
* the class Dumux::BoundaryConditions which distinguishes between
* dirichlet, process and neumann.
* BoundaryConditions for compositional transport can either be
* defined as saturations or as total concentrations (decoupled method).
* Each leads via pressure and constitutive relationships to the other
* and to the phase concentrations.
*/
struct BoundaryConditions2p2c
{
DUNE_DEPRECATED
enum Flags
{
saturation=1,
concentration=2
};
};
}
#endif
......@@ -69,7 +69,7 @@ private:
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents))
};
private:
protected:
Scalar mobility_[numPhases];
//numerical quantities
......@@ -83,7 +83,7 @@ private:
Scalar perimeter_;
FluidState fluidState_;
FluidState* fluidState_;
// FluxData fluxData_;
public:
......@@ -96,7 +96,7 @@ public:
mobility_({0.0, 0.0}),
numericalDensity_({0.0, 0.0}), volumeError_(0.), errorCorrection_(0.),
dv_dp_(0.), dv_({0.0, 0.0}), volumeDerivativesAvailable_(false),
perimeter_(0.)
perimeter_(0.),fluidState_(0)
{
}
......@@ -111,14 +111,22 @@ public:
// }
void setFluidState(FluidState& fluidState)
DUNE_DEPRECATED
{
fluidState_ = fluidState;
fluidState_ = &fluidState;
}
// const FluidState& fluidState() const
// {
// return fluidState_;
// }
const FluidState& fluidState() const
{
return *fluidState_;
}
FluidState& manipulateFluidState()
{
if(!fluidState_)
fluidState_ = new FluidState;
return *fluidState_;
}
const bool hasVolumeDerivatives() const
{ return volumeDerivativesAvailable_;}
......@@ -134,12 +142,12 @@ public:
////////////////////////////////////////////////////////////
Scalar pressure(int phaseIdx)
{
return fluidState_.pressure(phaseIdx);
return fluidState_->pressure(phaseIdx);
}
const Scalar& pressure(int phaseIdx) const
{
return fluidState_.pressure(phaseIdx);
return fluidState_->pressure(phaseIdx);
}
void setPressure(int phaseIdx, Scalar press)
......@@ -150,16 +158,16 @@ public:
//! Return saturation vector
void setTotalConcentration(int compIdx, Scalar value)
{
fluidState_.setMassConcentration(compIdx, value);
fluidState_->setMassConcentration(compIdx, value);
}
const Scalar totalConcentration(int compIdx) const
{
return fluidState_.massConcentration(compIdx);
return fluidState_->massConcentration(compIdx);
}
const Scalar massConcentration(int compIdx) const
{
return fluidState_.massConcentration(compIdx);
return fluidState_->massConcentration(compIdx);
}
//////////////////////////////////////////////////////////////
......@@ -253,30 +261,30 @@ public:
//! Return saturation vector
void setSaturation(int phaseIdx, Scalar value)
{
fluidState_.setSaturation(phaseIdx, value);
fluidState_->setSaturation(phaseIdx, value);
}
const Scalar saturation(int phaseIdx) const
{
return fluidState_.saturation(phaseIdx);
return fluidState_->saturation(phaseIdx);
}
//! Return density vector
void setViscosity(int phaseIdx, Scalar value)
{
fluidState_.setViscosity(phaseIdx, value);
fluidState_->setViscosity(phaseIdx, value);
}
const Scalar viscosity(int phaseIdx) const
{
return fluidState_.viscosity(phaseIdx);
return fluidState_->viscosity(phaseIdx);
}
//! Return capillary pressure vector
const Scalar capillaryPressure() const
{
return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx);
return fluidState_->pressure(nPhaseIdx) - fluidState_->pressure(wPhaseIdx);
}
void setCapillaryPressure(Scalar pc)
......@@ -287,30 +295,30 @@ public:
//! Return density vector
const Scalar density(int phaseIdx) const
{
return (fluidState_.density(phaseIdx));
return (fluidState_->density(phaseIdx));
}
//! Return density vector
const Scalar massFraction(int phaseIdx, int compIdx) const
{
return fluidState_.massFraction(phaseIdx, compIdx);
return fluidState_->massFraction(phaseIdx, compIdx);
}
//! Return density vector
const Scalar moleFraction(int phaseIdx, int compIdx) const
{
return fluidState_.moleFraction(phaseIdx, compIdx);
return fluidState_->moleFraction(phaseIdx, compIdx);
}
//! Return temperature
const Scalar temperature(int phaseIdx) const
{
return fluidState_.temperature(phaseIdx);
return fluidState_->temperature(phaseIdx);
}
//! Return phase phase mass fraction
const Scalar phaseMassFraction(int phaseIdx) const
{
return fluidState_.phaseMassFraction(phaseIdx);
return fluidState_->phaseMassFraction(phaseIdx);
}
};
......
/*****************************************************************************
* Copyright (C) 2011 by Benjamin Faigle *
* 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/>. *
*****************************************************************************/
#ifndef DUMUX_ELEMENTDATA2P2C_MULTYPHYSICS_HH
#define DUMUX_ELEMENTDATA2P2C_MULTYPHYSICS_HH
#include "2p2cproperties.hh"
#include "cellData2p2c.hh"
#include <dumux/decoupled/2p2c/pseudo1p2cfluidstate.hh>
/**
* @file
* @brief Class including the variables and data of discretized data of the constitutive relations for one element
* @author Markus Wolff
*/
namespace Dumux
{
/*!
* \ingroup IMPES
*/
//! Class including the variables and data of discretized data of the constitutive relations for one element.
/*! TODO: The variables of two-phase flow, which are one pressure and one saturation are stored in this class.
* Additionally, a velocity needed in the transport part of the decoupled two-phase flow is stored, as well as discretized data of constitutive relationships like
* mobilities, fractional flow functions and capillary pressure. Thus, they have to be callculated just once in every time step or every iteration step.
*
* @tparam TypeTag The Type Tag
1*/
template<class TypeTag>
class CellData2P2Cmultiphysics : public CellData2P2C<TypeTag>
{
private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
// typedef FluxData2P2C<TypeTag> FluxData;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
typedef PseudoOnePTwoCFluidState<TypeTag> SimpleFluidState;
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
enum
{
wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx,
wCompIdx = Indices::wPhaseIdx, nCompIdx = Indices::nPhaseIdx,
contiWEqIdx = Indices::contiWEqIdx, contiNEqIdx = Indices::contiNEqIdx
};
enum
{
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents))
};
enum
{
complex = 0, simple = 1
};
private:
int subdomain_;
int fluidStateType_;
SimpleFluidState* simpleFluidState_;
// FluxData fluxData_;
public:
//! Constructs a VariableClass object
/**
* @param gridView a DUNE gridview object corresponding to diffusion and transport equation
*/
CellData2P2Cmultiphysics() : CellData2P2C<TypeTag>(),
subdomain_(2), fluidStateType_(complex), simpleFluidState_(0)
{
}
// FluxData& setFluxData()
// {
// return fluxData_;
// }
//
// const FluxData& fluxData() const
// {
// return fluxData_;
// }
void setSimpleFluidState(SimpleFluidState& simpleFluidState)
{
assert (this->subdomain() != 2);
fluidStateType_ = simple;
simpleFluidState_ = &simpleFluidState;
}
SimpleFluidState& manipulateSimpleFluidState()
{
assert (this->subdomain() != 2);
fluidStateType_ = simple;
if(!simpleFluidState_)
simpleFluidState_ = new SimpleFluidState;
return *simpleFluidState_;
}
// const FluidState& fluidState() const
// {
// return *fluidState_;
// }
FluidState& manipulateFluidState()
{
assert(this->subdomain() == 2);
fluidStateType_ = complex;
if(!this->fluidState_)
this->fluidState_ = new FluidState;
return *this->fluidState_;
}
////////////////////////////////////////////////////////////
// functions returning primary variables
////////////////////////////////////////////////////////////
Scalar pressure(int phaseIdx)
{
if(fluidStateType_ == simple)
{
return simpleFluidState_->pressure(phaseIdx);
}
else
return this->fluidState_->pressure(phaseIdx);
}
const Scalar& pressure(int phaseIdx) const
{
if(fluidStateType_ == simple)
{
return simpleFluidState_->pressure(phaseIdx);
}
else
return this->fluidState_->pressure(phaseIdx);
}
//! Return saturation vector
void setTotalConcentration(int compIdx, Scalar value)
{
if(fluidStateType_ == simple)
return simpleFluidState_->setMassConcentration(compIdx, value);
else
return this->fluidState_->setMassConcentration(compIdx, value);
}
const Scalar totalConcentration(int compIdx) const
{
if(fluidStateType_ == simple)
return simpleFluidState_->massConcentration(compIdx);
else
return this->fluidState_->massConcentration(compIdx);
}
const Scalar massConcentration(int compIdx) const
{
if(fluidStateType_ == simple)
return simpleFluidState_->massConcentration(compIdx);
else
return this->fluidState_->massConcentration(compIdx);
}
//////////////////////////////////////////////////////////////
// functions returning the vectors of secondary variables
//////////////////////////////////////////////////////////////
//! Return subdomain information
int& subdomain()
{
return subdomain_;
}
const int& subdomain() const
{
return subdomain_;
}
/*** b) from fluidstate ***/
//! Return saturation vector
void setSaturation(int phaseIdx, Scalar value)
{
assert(this->subdomain() == 2);
this->fluidState_->setSaturation(phaseIdx, value);
}
const Scalar saturation(int phaseIdx) const
{
if(fluidStateType_ == simple)
{
return simpleFluidState_->saturation(phaseIdx);
}
else
return this->fluidState_->saturation(phaseIdx);
}
//! Return density vector
void setViscosity(int phaseIdx, Scalar value)
{
if(fluidStateType_ == simple)
{
assert(phaseIdx == simpleFluidState_->presentPhaseIdx());
simpleFluidState_->setViscosity(phaseIdx, value);
}
else
this->fluidState_->setViscosity(phaseIdx, value);
}
const Scalar viscosity(int phaseIdx) const
{
if(fluidStateType_ == simple)
{
if(phaseIdx != simpleFluidState_->presentPhaseIdx())
return 0.; // This should only happend for output
return simpleFluidState_->viscosity(phaseIdx);
}
else
return this->fluidState_->viscosity(phaseIdx);
}
//! Return capillary pressure vector
const Scalar capillaryPressure() const
{
if(fluidStateType_ == simple)
return 0.;
else
return this->fluidState_->pressure(nPhaseIdx) - this->fluidState_->pressure(wPhaseIdx);
}
//! Return density vector
const Scalar density(int phaseIdx) const
{
if(fluidStateType_ == simple)
{
return simpleFluidState_->density(phaseIdx);
}
else
return this->fluidState_->density(phaseIdx);
}
//! Return the mass fraction
const Scalar massFraction(int phaseIdx, int compIdx) const
{
if(fluidStateType_ == simple)
{
return simpleFluidState_->massFraction(phaseIdx, compIdx);
}
else
return this->fluidState_->massFraction(phaseIdx, compIdx);
}
//! Return the mole fraction
const Scalar moleFraction(int phaseIdx, int compIdx) const
{
if(fluidStateType_ == simple)
{
return simpleFluidState_->moleFraction(phaseIdx, compIdx);
}
else
return this->fluidState_->moleFraction(phaseIdx, compIdx);
}
//! Return temperature
const Scalar temperature(int phaseIdx) const
{
if(fluidStateType_ == simple)
{
return simpleFluidState_->temperature(phaseIdx);
}
else
return this->fluidState_->temperature(phaseIdx);
}
//! Return phase phase mass fraction
const Scalar phaseMassFraction(int phaseIdx) const
{
if(fluidStateType_ == simple)
{
if(phaseIdx == simpleFluidState_->presentPhaseIdx())
return 1.; // phase is present => nu = 1
else
return 0.;
}
else
return this->fluidState_->phaseMassFraction(phaseIdx);
}
};
}
#endif
......@@ -238,10 +238,8 @@ void FVPressure2P2C<TypeTag>::getStorage(Dune::FieldVector<Scalar, 2>& storageEn
int globalIdxI = problem().variables().index(elementI);
Scalar volume = elementI.geometry().volume();
storageEntry = 0.;
// determine maximum error to scale error-term
Scalar timestep_ = problem().timeManager().timeStepSize();
// maxErr /= timestep_;
// compressibility term
if (!first && timestep_ != 0)
......@@ -798,10 +796,6 @@ void FVPressure2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& en
template<class TypeTag>
void FVPressure2P2C<TypeTag>::updateMaterialLaws()
{
// instantiate a brandnew fluid state object
FluidState fluidState;
Scalar maxError = 0.;
// iterate through leaf grid an evaluate c0 at cell center
ElementIterator eItEnd = problem().gridView().template end<0> ();
......@@ -822,15 +816,15 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
template<class TypeTag>
void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& elementI)
{
// instantiate a brandnew fluid state object
FluidState fluidState;
// get global coordinate of cell center
GlobalPosition globalPos = elementI.geometry().center();
int globalIdx = problem().variables().index(elementI);
CellData& cellData = problem().variables().cellData(globalIdx);
// acess the fluid state and prepare for manipulation
FluidState& fluidState = cellData.manipulateFluidState();
Scalar temperature_ = problem().temperatureAtPos(globalPos);
// reset
cellData.reset();
......@@ -938,11 +932,6 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
if(iterout !=0)
Dune::dinfo << iterout << "times iteration of pc was applied at Idx " << globalIdx << ", pc delta still " << abs(oldPc-pc) << std::endl;
}
// secondary variables
// initialize saturation, capillary pressure
cellData.setFluidState(fluidState);
// initialize phase properties not stored in fluidstate
cellData.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, wPhaseIdx));
cellData.setViscosity(nPhaseIdx, FluidSystem::viscosity(fluidState, nPhaseIdx));
......
......@@ -46,7 +46,7 @@ class PseudoOnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTa
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, TwoPTwoCIndices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
public:
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
......@@ -73,16 +73,18 @@ public:
*
* \param Z1 Feed mass fraction \f$\mathrm{[-]}\f$
* \param press1p Pressure value for present phase \f$\mathrm{[Pa]}\f$
* \param satW Saturation of the wetting phase \f$\mathrm{[-]}\f$
* \param presentPhaseIdx Subdomain Index = Indication which phase is present
* \param temperature Temperature \f$\mathrm{[K]}\f$
*/
void update(const Scalar& Z1,const Scalar& press1p,const Scalar& satW, const Scalar& temperature)
void update(const Scalar& Z1,const Scalar& press1p,const int presentPhaseIdx, const Scalar& temperature)
{
Sw_ = satW;
pressure1p_=press1p;
temperature_ = temperature;
if (satW == 1.)