Commit 1a3b5a7c authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/mpnc-on-next' into 'next'

Feature/mpnc on next

See merge request !666
parents 25688435 3547abd5
......@@ -181,6 +181,7 @@ NEW_PROP_TAG(OnlyGasPhaseCanDisappear); //!< reduces the phasestates to threePha
/////////////////////////////////////////////////////////////
// Properties used by the staggered-grid discretization method
/////////////////////////////////////////////////////////////
NEW_PROP_TAG(NumEqCellCenter); //!< The number of equations for cell-centered dofs
NEW_PROP_TAG(NumEqFace); //!< The number of equations for face dofs
NEW_PROP_TAG(CellCenterSolutionVector); //!< The solution vector type for cell-centered dofs
......@@ -196,6 +197,31 @@ NEW_PROP_TAG(BaseEpsilon); //!< A base epsilon for numer
NEW_PROP_TAG(FaceVariables); //!< Class containing local face-related data
NEW_PROP_TAG(BoundaryValues); //!< Class containing local boundary data
/////////////////////////////////////////////////////////////
// Properties used by the mpnc model
/////////////////////////////////////////////////////////////
NEW_PROP_TAG(PressureFormulation); //! the formulation of the pressure e.g most wetting first
/////////////////////////////////////////////////////////////
// Properties used by the nonequilibrium model
/////////////////////////////////////////////////////////////
NEW_PROP_TAG(EquilibriumLocalResidual);
NEW_PROP_TAG(EquilibriumIndices);
NEW_PROP_TAG(EquilibriumVtkOutputFields);
NEW_PROP_TAG(NumEqBalance);
NEW_PROP_TAG(EnableThermalNonEquilibrium);
NEW_PROP_TAG(EnableChemicalNonEquilibrium);
NEW_PROP_TAG(NumEnergyEqFluid);
NEW_PROP_TAG(NumEnergyEqSolid);
NEW_PROP_TAG(AwnSurface);
NEW_PROP_TAG(AwsSurface);
NEW_PROP_TAG(AnsSurface);
NEW_PROP_TAG(NusseltFormulation);
NEW_PROP_TAG(SherwoodFormulation);
} // end namespace Properties
} // end namespace Dumux
......
......@@ -18,81 +18,117 @@
*****************************************************************************/
/*!
* \file
*
* \brief This class contains the volume variables required for the
* modules which require the specific interfacial area between
* fluid phases.
* \brief This file contains the data which is required to calculate
* energy fluxes due to molecular diffusion with Fourier's law.
*/
#ifndef DUMUX_MPNC_VOLUME_VARIABLES_IA_HH
#define DUMUX_MPNC_VOLUME_VARIABLES_IA_HH
#ifndef DUMUX_DISCRETIZATION_BOX_FOURIERS_LAW_NONEQUILIBRIUM_HH
#define DUMUX_DISCRETIZATION_BOX_FOURIERS_LAW_NONEQUILIBRIUM_HH
#include <dumux/porousmediumflow/mpnc/implicit/properties.hh>
#include <dune/common/float_cmp.hh>
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/methods.hh>
namespace Dumux
{
/*!
* \brief This class contains the volume variables required for the
* modules which require the specific interfacial area between
* fluid phases.
*
* This is the specialization for the cases which do _not_ require
* specific interfacial area.
* \ingroup BoxFouriersLaw
* \brief Specialization of Fourier's Law for the box method for thermal nonequilibrium models.
*/
template <class TypeTag, bool enableKinetic, int numEnergyEquations>
class MPNCVolumeVariablesIA
template <class TypeTag>
class FouriersLawNonEquilibriumImplementation<TypeTag, DiscretizationMethods::Box>
{
static_assert(((numEnergyEquations < 1) && !enableKinetic),
"The kinetic energy modules need specific interfacial area "
"but no suitable specialization of the IA volume variables module "
"has been included.");
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
using ParameterCache = typename FluidSystem::ParameterCache;
using IndexType = typename GridView::IndexSet::IndexType;
using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
enum {dimWorld=GridView::dimensionworld};
enum { dim = GridView::dimension} ;
enum { dimWorld = GridView::dimensionworld} ;
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
enum { numEnergyEqFluid = GET_PROP_VALUE(TypeTag, NumEnergyEqFluid) };
enum {sPhaseIdx = FluidSystem::sPhaseIdx};
using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
/*!
* \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
*
* \param volVars The volume variables
* \param fluidState Container for all the secondary variables concerning the fluids
* \param paramCache Container for cache parameters
* \param priVars The primary Variables
* \param problem The problem
* \param element The finite element
* \param fvGeometry The finite-volume geometry in the fully implicit scheme
* \param scvIdx The index of the sub-control volumete element
*
*/
void update(const VolumeVariables & volVars,
const FluidState &fluidState,
const ParameterCache &paramCache,
const PrimaryVariables &priVars,
const Problem &problem,
const Element & element,
const FVElementGeometry & fvGeometry,
const unsigned int scvIdx)
static Scalar flux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const int phaseIdx,
const ElementFluxVariablesCache& elemFluxVarsCache)
{
}
// get inside and outside diffusion tensors and calculate the harmonic mean
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
const auto& insideVolVars = elemVolVars[insideScv];
const auto& outsideVolVars = elemVolVars[outsideScv];
Scalar insideLambda = 0.0;
Scalar outsideLambda = 0.0;
// effective diffusion tensors
if (phaseIdx != sPhaseIdx)
{
if (numEnergyEqFluid == 1)
{ //when only one energy equation for fluids is used, we need an effective law for that
insideLambda += ThermalConductivityModel::effectiveThermalConductivity(insideVolVars, problem.spatialParams(), element, fvGeometry, insideScv);
outsideLambda += ThermalConductivityModel::effectiveThermalConductivity(outsideVolVars, problem.spatialParams(), element, fvGeometry, outsideScv);
}
else
{
insideLambda += insideVolVars.fluidThermalConductivity(phaseIdx)*insideVolVars.saturation(phaseIdx)*insideVolVars.porosity();
outsideLambda += outsideVolVars.fluidThermalConductivity(phaseIdx)*outsideVolVars.saturation(phaseIdx)*outsideVolVars.porosity();
}
}
//solid phase
else
{
insideLambda += insideVolVars.solidThermalConductivity()*(1-insideVolVars.porosity());
outsideLambda +=outsideVolVars.solidThermalConductivity()*(1-outsideVolVars.porosity());
}
/*!
* \brief If running in valgrind this makes sure that all
* quantities in the volume variables are defined.
*/
void checkDefined() const
{ }
// scale by extrusion factor
insideLambda *= insideVolVars.extrusionFactor();
outsideLambda *= outsideVolVars.extrusionFactor();
// the resulting averaged diffusion tensor
const auto lambda = problem.spatialParams().harmonicMean(insideLambda, outsideLambda, scvf.unitOuterNormal());
// evaluate gradTemp at integration point
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
GlobalPosition gradTemp(0.0);
for (auto&& scv : scvs(fvGeometry))
{
// compute the temperature gradient with the shape functions
if (phaseIdx < numEnergyEqFluid)
gradTemp.axpy(elemVolVars[scv].temperature(phaseIdx), fluxVarsCache.gradN(scv.indexInElement()));
else
gradTemp.axpy(elemVolVars[scv].temperatureSolid(), fluxVarsCache.gradN(scv.indexInElement()));
}
// comute the heat conduction flux
return -1.0*vtmv(scvf.unitOuterNormal(), lambda, gradTemp)*scvf.area();
}
};
} // namespace Dumux
} // end namespace Dumux
#endif
......@@ -16,44 +16,32 @@
* 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 Test for the MpNc CC model.
* \brief This file contains the data which is required to calculate
* diffusive mass fluxes due to molecular diffusion with Fourier's law.
*/
#include <config.h>
#include "obstacleproblem.hh"
#include <dumux/common/start.hh>
#ifndef DUMUX_DISCRETIZATION_FOURIERS_LAW_NONEQUILIBRIUM_HH
#define DUMUX_DISCRETIZATION_FOURIERS_LAW_NONEQUILIBRIUM_HH
#include <dumux/discretization/methods.hh>
namespace Dumux
{
// forward declaration
template <class TypeTag, DiscretizationMethods Method>
class FouriersLawNonEquilibriumImplementation
{};
/*!
* \brief Provides an interface for customizing error messages associated with
* reading in parameters.
*
* \param progName The name of the program, that was tried to be started.
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
* \ingroup FouriersLaw
* \brief Evaluates the heat conduction flux according to Fouriers's law
*/
void usage(const char *progName, const std::string &errorMsg)
{
if (errorMsg.size() > 0) {
std::string errorMessageOut = "\nUsage: ";
errorMessageOut += progName;
errorMessageOut += " [options]\n";
errorMessageOut += errorMsg;
errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.File Name of the file containing the grid \n"
"\t definition in DGF format\n";
template <class TypeTag>
using FouriersLawNonEquilibrium = FouriersLawNonEquilibriumImplementation<TypeTag, GET_PROP_VALUE(TypeTag, DiscretizationMethod)>;
std::cout << errorMessageOut
<< "\n";
}
}
} // end namespace Dumux
int main(int argc, char** argv)
{
using ProblemTypeTag = TTAG(ObstacleCCProblem);
return Dumux::start<ProblemTypeTag>(argc, argv, usage);
}
#include <dumux/discretization/box/fourierslawnonequilibrium.hh>
#endif
......@@ -25,23 +25,25 @@
#define THERMALCONDUCTIVITY_SIMPLE_FLUID_LUMPING_HH
#include <algorithm>
#include <dumux/porousmediumflow/mpnc/implicit/properties.hh>
#include <dumux/common/properties.hh>
namespace Dumux
{
struct SimpleLumpingIndices
{
static const int wPhaseIdx = 0;
static const int nPhaseIdx = 1;
};
/*!
* \ingroup fluidmatrixinteractionslaws
* \brief Relation for the saturation-dependent effective thermal conductivity
* \todo This shouldn't depend on TypeTag!!
*/
template<class TypeTag>
template<class TypeTag, class Scalar, class Indices = SimpleLumpingIndices>
class ThermalConductivitySimpleFluidLumping
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { numEnergyEquations = GET_PROP_VALUE(TypeTag, NumEnergyEquations)};
enum { numEnergyEquationsFluid = GET_PROP_VALUE(TypeTag, NumEnergyEqFluid)};
public:
/*!
......@@ -55,12 +57,12 @@ public:
*
* \return effective thermal conductivity \f$\mathrm{[W/(m K)]}\f$
*/
template<class VolumeVariables, class SpatialParams, class Element, class FVGeometry>
template<class VolumeVariables, class SpatialParams, class Element, class FVGeometry, class SubControlVolume>
static Scalar effectiveThermalConductivity(const VolumeVariables& volVars,
const SpatialParams& spatialParams,
const Element& element,
const FVGeometry& fvGeometry,
int scvIdx)
SubControlVolume& scv)
{
Scalar sw = volVars.saturation(Indices::wPhaseIdx);
Scalar lambdaW = volVars.fluidThermalConductivity(Indices::wPhaseIdx);
......@@ -90,7 +92,7 @@ public:
const Scalar porosity,
const Scalar rhoSolid = 0.0 /*unused*/)
{
assert(numEnergyEquations != 3) ;
assert(numEnergyEquationsFluid != 2) ;
// Franz Lindner / Shi & Wang 2011
using std::max;
......@@ -100,7 +102,7 @@ public:
Scalar keff ;
if (numEnergyEquations == 2){ // solid dealed with individually (extra balance equation)
if (numEnergyEquationsFluid == 1){ // solid dealed with individually (extra balance equation)
keff = kfeff ;
}
else {
......
......@@ -23,8 +23,8 @@
* multi-phase, multi-component fluid system without using
* any assumptions.
*/
#ifndef DUMUX_GENERIC_FLUID_STATE_HH
#define DUMUX_GENERIC_FLUID_STATE_HH
#ifndef DUMUX_NONEQUILIBRIUM_FLUID_STATE_HH
#define DUMUX_NONEQUILIBRIUM_FLUID_STATE_HH
#include <cmath>
#include <algorithm>
......@@ -35,6 +35,7 @@
namespace Dumux
{
/*!
* \ingroup FluidStates
* \brief Represents all relevant thermodynamic quantities of a
......@@ -45,8 +46,8 @@ template <class Scalar, class FluidSystem>
class NonEquilibriumFluidState
{
public:
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
static constexpr int numPhases = FluidSystem::numPhases;
static constexpr int numComponents = FluidSystem::numComponents;
NonEquilibriumFluidState()
{
......@@ -116,7 +117,14 @@ public:
* @copydoc Dumux::CompositionalFluidState::fugacity()
*/
Scalar fugacity(int phaseIdx, int compIdx) const
{ return pressure_[phaseIdx]*fugacityCoefficient_[phaseIdx][compIdx]*moleFraction_[phaseIdx][compIdx]; }
{
return pressure_[phaseIdx]*fugacityCoefficient_[phaseIdx][compIdx]*moleFraction_[phaseIdx][compIdx];
}
Scalar fugacity(int compIdx) const
{
return fugacity(0, compIdx);
}
/*!
* @copydoc Dumux::CompositionalFluidState::molarVolume()
......@@ -140,8 +148,16 @@ public:
/*!
* @copydoc Dumux::CompositionalFluidState::temperature()
*/
Scalar temperature(int phaseIdx) const
{ return temperature_[phaseIdx]; }
Scalar temperature(const int phaseIdx) const
{ return temperature_[phaseIdx]; }
/*!
* \brief Get the equilibrium temperature \f$\mathrm{[K]}\f$ of the fluid phases.
*/
Scalar temperature() const
{
return temperatureEquil_ ;
}
/*!
* @copydoc Dumux::CompositionalFluidState::pressure()
......@@ -178,13 +194,17 @@ public:
* \brief Set the temperature \f$\mathrm{[K]}\f$ of a fluid phase
*/
void setTemperature(int phaseIdx, Scalar value)
{ temperature_[phaseIdx] = value; }
{
temperature_[phaseIdx] = value;
}
/*!
* \brief Set the temperature \f$\mathrm{[K]}\f$ of all fluid phases.
*/
void setTemperature(Scalar value)
{DUNE_THROW(Dune::NotImplemented, "This is a fluidstate for non-equilibrium, temperature in which phase?");}
{
temperatureEquil_ = value;
}
/*!
* \brief Set the fluid pressure of a phase \f$\mathrm{[Pa]}\f$
......@@ -319,7 +339,8 @@ protected:
Scalar density_[numPhases];
Scalar enthalpy_[numPhases];
Scalar viscosity_[numPhases];
Scalar temperature_[numPhases];
Scalar temperature_[numPhases]; //
Scalar temperatureEquil_;
};
} // end namespace Dumux
......
......@@ -18,3 +18,4 @@ add_subdirectory("richards")
add_subdirectory("richardsnc")
add_subdirectory("sequential")
add_subdirectory("tracer")
add_subdirectory("nonequilibrium")
......@@ -61,6 +61,7 @@ class PorousMediumFluxVariables : public FluxVariablesBase<TypeTag>
static constexpr bool enableAdvection = GET_PROP_VALUE(TypeTag, EnableAdvection);
static constexpr bool enableMolecularDiffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
static constexpr bool enableEnergyBalance = GET_PROP_VALUE(TypeTag, EnableEnergyBalance);
static constexpr bool enableThermalNonEquilibrium = GET_PROP_VALUE(TypeTag, EnableThermalNonEquilibrium);
public:
......@@ -115,6 +116,7 @@ public:
}
}
template <bool enable = !enableThermalNonEquilibrium, typename std::enable_if_t<enable, int> = 0>
Scalar heatConductionFlux() const
{
if (enableEnergyBalance)
......@@ -132,6 +134,19 @@ public:
}
}
template <bool enable = enableThermalNonEquilibrium, typename std::enable_if_t<enable, int> = 0>
Scalar heatConductionFlux(const int phaseIdx) const
{
return HeatConductionType::flux(this->problem(),
this->element(),
this->fvGeometry(),
this->elemVolVars(),
this->scvFace(),
phaseIdx,
this->elemFluxVarsCache());
}
private:
//! simple caching if advection flux is used twice with different upwind function
mutable std::bitset<numPhases> advFluxIsCached_;
......
add_subdirectory("implicit")
#install headers
install(FILES
indices.hh
localresidual.hh
model.hh
volumevariables.hh
vtkoutputfield.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/mpnc)
add_subdirectory("diffusion")
add_subdirectory("energy")
add_subdirectory("mass")
#install headers
install(FILES
fluxvariables.hh
indices.hh
localresidual.hh
model.hh
modelkinetic.hh
newtoncontroller.hh
properties.hh
propertieskinetic.hh
propertydefaults.hh
propertydefaultskinetic.hh
volumevariables.hh
volumevariablesia.hh
volumevariablesiakinetic.hh
vtkwritercommon.hh
vtkwriter.hh
vtkwritermodule.hh
velomodelnewtoncontroller.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/mpnc/implicit)
#install headers
install(FILES
diffusion.hh
fluxvariables.hh
volumevariables.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/mpnc/implicit/diffusion)
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* 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 parts to calculate the diffusive flux in
* the fully coupled MpNc model
*/
#ifndef DUMUX_MPNC_DIFFUSION_HH
#define DUMUX_MPNC_DIFFUSION_HH
#include <dune/common/float_cmp.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dumux/porousmediumflow/mpnc/implicit/properties.hh>
namespace Dumux {
/*!
* \brief Calculates the diffusive flux in the fully coupled MpNc model.
* Called from the mass module.
*/
template <class TypeTag, bool enableDiffusion>
class MPNCDiffusion
{
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents)};
enum { nPhaseIdx = FluidSystem::nPhaseIdx };
enum { wPhaseIdx = FluidSystem::wPhaseIdx };
enum { nCompIdx = FluidSystem::nCompIdx };
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using ComponentMatrix = Dune::FieldMatrix<Scalar, numComponents, numComponents>;
using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
public:
/*!
* \brief Calls the function for the diffusion in the gas and liquid phases, respectively.
*
* In the gas phase, the mutual influence of mole fractions can be considered.
*
* \param fluxes The Diffusive flux over the sub-control-volume face for each component
* \param phaseIdx The index of the phase we are calculating the diffusive flux for
* \param fluxVars The flux variables at the current sub control volume face
* \param molarDensity The (molar) density of the phase
*/
static void flux(ComponentVector &fluxes,
const unsigned int phaseIdx,
const Fl