Commit 5884e5e3 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

2p2cni box model: bring it up to stable standards

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4185 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 3fb309ae
......@@ -206,9 +206,9 @@ public:
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up =
this->curVolVars_(vars.upstreamIdx(phaseIdx));
const VolumeVariables &dn = this->curVolVars_()[vars.downstreamIdx(
phaseIdx)];
this->curVolVars_(vars.upstreamIdx(phaseIdx));
const VolumeVariables &dn =
this->curVolVars_(vars.downstreamIdx(phaseIdx));
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
......
......@@ -210,6 +210,17 @@ public:
*/
};
/*!
* \brief Returns the relative weight of a primary variable for
* calculating relative errors.
*/
Scalar primaryVarWeight(int vertIdx, int pvIdx) const
{
if (Indices::pressureIdx == pvIdx)
return 1e-5;
return 1;
}
/*!
* \brief Called by the BoxModel's update method.
*/
......
......@@ -146,7 +146,7 @@ SET_BOOL_PROP(BoxTwoPTwoC, EnableJacobianRecycling, true);
// enable partial reassembling by default
SET_BOOL_PROP(BoxTwoPTwoC, EnablePartialReassemble, true);
// enable time-step ramp up by default
SET_BOOL_PROP(BoxTwoPTwoC, EnableTimeStepRampUp, true);
SET_BOOL_PROP(BoxTwoPTwoC, EnableTimeStepRampUp, false);
// \}
}
......
......@@ -61,6 +61,8 @@ class TwoPTwoCNILocalResidual : public TwoPTwoCLocalResidual<TypeTag>
dimWorld = GridView::dimensionworld,
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
energyEqIdx = Indices::energyEqIdx,
temperatureIdx = Indices::temperatureIdx,
lPhaseIdx = Indices::lPhaseIdx,
......@@ -99,7 +101,7 @@ public:
const VolumeVariables &vertDat = elemDat[scvIdx];
// compute the energy storage
result[temperatureIdx] =
result[energyEqIdx] =
vertDat.porosity()*(vertDat.density(lPhaseIdx) *
vertDat.internalEnergy(lPhaseIdx) *
//vertDat.enthalpy(lPhaseIdx) *
......@@ -110,8 +112,7 @@ public:
//vertDat.enthalpy(gPhaseIdx) *
vertDat.saturation(gPhaseIdx))
+
vertDat.temperature() *
vertDat.heatCapacity();
vertDat.temperature()*vertDat.heatCapacity();
}
/*!
......@@ -128,13 +129,13 @@ public:
ParentType::computeAdvectiveFlux(flux, fluxData);
// advective heat flux in all phases
flux[temperatureIdx] = 0;
flux[energyEqIdx] = 0;
for (int phase = 0; phase < numPhases; ++phase) {
// vertex data of the upstream and the downstream vertices
const VolumeVariables &up = this->curVolVars_(fluxData.upstreamIdx(phase));
const VolumeVariables &dn = this->curVolVars_(fluxData.downstreamIdx(phase));
flux[temperatureIdx] +=
flux[energyEqIdx] +=
fluxData.KmvpNormal(phase) * (
mobilityUpwindAlpha * // upstream vertex
( up.density(phase) *
......
// $Id$
/*****************************************************************************
* Copyright (C) 2008-2009 by Melanie Darcis *
* Copyright (C) 2008-2009 by Klaus Mosthaf *
* Copyright (C) 2008-2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* Copyright (C) 2008,2009 by Melanie Darcis *
* Klaus Mosthaf, *
* Andreas Lauser, *
* Bernd Flemisch *
* Copyright (C) 2008-2009 by Bernd Flemisch *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
......
......@@ -74,8 +74,7 @@ class TwoPNILocalResidual : public TwoPLocalResidual<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
static const Scalar mobilityUpwindAlpha = GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
......@@ -120,29 +119,39 @@ public:
* This method is called by compute flux (base class)
*/
void computeAdvectiveFlux(PrimaryVariables &flux,
const FluxVariables &fluxData) const
const FluxVariables &fluxVars) const
{
// advective mass flux
ParentType::computeAdvectiveFlux(flux, fluxData);
ParentType::computeAdvectiveFlux(flux, fluxVars);
// advective heat flux in all phases
flux[energyEqIdx] = 0;
for (int phase = 0; phase < numPhases; ++phase) {
// vertex data of the upstream and the downstream vertices
const VolumeVariables &up = this->curVolVars_(fluxData.upstreamIdx(phase));
const VolumeVariables &dn = this->curVolVars_(fluxData.downstreamIdx(phase));
flux[energyEqIdx] +=
fluxData.KmvpNormal(phase) * (
mobilityUpwindAlpha * // upstream vertex
( up.density(phase) *
up.mobility(phase) *
up.enthalpy(phase))
+
(1 - mobilityUpwindAlpha) * // downstream vertex
( dn.density(phase) *
dn.mobility(phase) *
dn.enthalpy(phase)) );
Vector tmpVec;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// calculate the flux in the normal direction of the
// current sub control volume face
fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(phaseIdx),
tmpVec);
Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
// data attached to upstream and the downstream vertices
// of the current phase
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
// add advective energy flux in current phase
flux[energyEqIdx] +=
normalFlux
*
(( mobilityUpwindAlpha)*
up.density(phaseIdx)*
up.mobility(phaseIdx)*
up.enthalpy(phaseIdx)
+
(1 - mobilityUpwindAlpha)*
dn.density(phaseIdx)*
dn.mobility(phaseIdx)*
dn.enthalpy(phaseIdx));
}
}
......
......@@ -28,8 +28,6 @@
#define DUMUX_2PNI_MODEL_HH
#include <dumux/boxmodels/2p/2pmodel.hh>
#include "2pnilocalresidual.hh"
#include "2pniproblem.hh"
namespace Dumux {
......@@ -86,4 +84,6 @@ class TwoPNIModel: public TwoPModel<TypeTag>
}
#include "2pnipropertydefaults.hh"
#endif
......@@ -20,14 +20,10 @@
*
* \brief Defines the properties required for the non-isotherm two-phase BOX model.
*/
#ifndef DUMUX_2PNI_PROPERTIES_HH
#define DUMUX_2PNI_PROPERTIES_HH
#include <dumux/boxmodels/2p/2pproperties.hh>
#include "2pnivolumevariables.hh"
#include "2pnifluxvariables.hh"
namespace Dumux
{
......@@ -36,36 +32,6 @@ namespace Dumux
*/
// \{
////////////////////////////////
// forward declarations
////////////////////////////////
template<class TypeTag>
class TwoPNIModel;
template<class TypeTag>
class TwoPNILocalResidual;
template <class TypeTag>
class TwoPNIVolumeVariables;
template <class TypeTag>
class TwoPNIFluxVariables;
/*!
* \brief Enumerations for the non-isothermal two-phase model
*/
template <int PVOffset = 0>
class TwoPNIIndices : public TwoPIndices<PVOffset>
{
public:
static const int temperatureIdx = PVOffset + 2; //! The primary variable index for temperature
static const int energyEqIdx = PVOffset + 2; //! The equation index of the energy equation
};
////////////////////////////////
// properties
////////////////////////////////
namespace Properties
{
//////////////////////////////////////////////////////////////////
......@@ -80,33 +46,7 @@ NEW_TYPE_TAG(BoxTwoPNI, INHERITS_FROM(BoxTwoP));
//////////////////////////////////////////////////////////////////
NEW_PROP_TAG(TwoPNIIndices); //!< Enumerations for the non-isothermal 2p models
//////////////////////////////////////////////////////////////////
// Properties
//////////////////////////////////////////////////////////////////
SET_INT_PROP(BoxTwoPNI, NumEq, 3); //!< set the number of equations to 3
//! Use the 2pni local jacobian operator for the 2pni model
SET_TYPE_PROP(BoxTwoPNI,
LocalResidual,
TwoPNILocalResidual<TypeTag>);
//! the Model property
SET_TYPE_PROP(BoxTwoPNI, Model, TwoPNIModel<TypeTag>);
//! the VolumeVariables property
SET_TYPE_PROP(BoxTwoPNI, VolumeVariables, TwoPNIVolumeVariables<TypeTag>);
//! the FluxVariables property
SET_TYPE_PROP(BoxTwoPNI, FluxVariables, TwoPNIFluxVariables<TypeTag>);
//! The indices required by the non-isothermal two-phase model
SET_TYPE_PROP(BoxTwoPNI, TwoPIndices, TwoPNIIndices<0>);
SET_TYPE_PROP(BoxTwoPNI, TwoPNIIndices, TwoPNIIndices<0>);
}
}
#endif
......@@ -423,12 +423,15 @@ public:
N2::gasEnthalpy(temperature, pN2);
}
else {
Scalar pWater = fluidState.fugacity(H2OIdx);
Scalar pN2 = fluidState.fugacity(N2Idx);
Scalar cH2O = fluidState.concentration(gPhaseIdx, H2OIdx);
Scalar cN2 = fluidState.concentration(gPhaseIdx, N2Idx);
Scalar pH2O = H2O::gasPressure(temperature, cH2O*H2O::molarMass());
Scalar pN2 = N2::gasPressure(temperature, cN2*N2::molarMass());
Scalar result = 0;
result +=
H2O::gasEnthalpy(temperature, pWater) *
H2O::gasEnthalpy(temperature, pH2O) *
fluidState.massFrac(gPhaseIdx, H2OIdx);
result +=
N2::gasEnthalpy(temperature, pN2) *
......
......@@ -2,7 +2,7 @@ DGF
Interval
0 0 % first corner
40 40 % second corner
16 16 % 24 cells in x and 16 in y direction
64 64 % 24 cells in x and 16 in y direction
#
GridParameter
overlap 0
......
......@@ -82,7 +82,7 @@ SET_TYPE_PROP(WaterAirProblem,
SET_BOOL_PROP(WaterAirProblem, EnableGravity, true);
// Write newton convergence
SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, true);
SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, false);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment