Skip to content
GitLab
Commit 005b6d58 authored by Timo Koch's avatar Timo Koch
Browse files

[1p] Update 1p model to new structure

parent ed70dad8
Hide whitespace changes
Inline Side-by-side
dumux/porousmediumflow/1p/implicit/localresidual.hh
View file @ 005b6d58
......@@ -63,7 +63,8 @@ public:
// retrieve the upwind weight for the mass conservation equations. Use the value
// specified via the property system as default, and overwrite
// it by the run-time parameter from the Dune::ParameterTree
upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
}
/*!
......@@ -73,22 +74,20 @@ public:
* model.
*
* This function should not include the source and sink terms.
* \param storage The phase mass within the sub-control volume
* \param scvIdx The SCV (sub-control-volume) index
* \param usePrevSol Evaluate function with solution of current or previous time step
* \param scv The sub control volume
* \param volVars The current or previous volVars
*
* The volVars can be different to allow computing the implicit euler time derivative here
*/
void computeStorage(PrimaryVariables &storage, const int scvIdx, const bool usePrevSol) const
PrimaryVariables computeStorage(const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
// if flag usePrevSol is set, the solution from the previous
// time step is used, otherwise the current solution is
// used. The secondary variables are used accordingly. This
// is required to compute the derivative of the storage term
// using the implicit euler method.
const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
const VolumeVariables &volVars = elemVolVars[scvIdx];
PrimaryVariables storage(0);
// partial time derivative of the wetting phase mass
storage[conti0EqIdx] = volVars.density() * volVars.porosity();
storage[conti0EqIdx] = volVars.density() * volVars.porosity();
return storage;
}
......@@ -101,57 +100,26 @@ public:
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
void computeFlux(PrimaryVariables &flux, const int fIdx, const bool onBoundary=false) const
PrimaryVariables computeFlux(const SubControlVolumeFace& scvFace) const
{
FluxVariables fluxVars;
fluxVars.update(this->problem_(),
this->element_(),
this->fvGeometry_(),
fIdx,
this->curVolVars_(),
onBoundary);
asImp_()->computeAdvectiveFlux(flux, fluxVars);
asImp_()->computeDiffusiveFlux(flux, fluxVars);
}
const auto& fluxVars = this->model_().fluxVars(scvFace);
/*!
* \brief Evaluate the advective mass flux of all components over
* a face of a sub-control volume.
*
* \param flux The advective flux over the sub-control-volume face for each component
* \param fluxVars The flux variables at the current SCV
*/
void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(/*phaseIdx=*/0));
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(/*phaseIdx=*/0));
flux[conti0EqIdx] =
(( upwindWeight_)*up.density()
+
(1 - upwindWeight_)*dn.density())
*
fluxVars.volumeFlux(/*phaseIdx=*/0);
}
auto massWeight = massUpwindWeight_;
auto mobWeight = mobilityUpwindWeight_;
/*!
* \brief Adds the diffusive mass flux of all components over
* a face of a sub-control volume.
*
* \param flux The diffusive flux over the sub-control-volume face for each component
* \param fluxVars The flux variables at the current SCV
*/
void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
{
// diffusive fluxes
flux += 0.0;
auto upwindRule = [massWeight, mobWeight](const VolumeVariables& up, const VolumeVariables& dn)
{ return (up.density(0)*massWeight + dn.density(0)*(1-massWeight))
*(up.mobility(0)*mobWeight + dn.density(0)*(1-mobWeight)); };
auto flux = fluxVars.darcyFluxVars().computeFlux(0, upwindRule);
return flux;
}
/*!
* \brief Return the temperature given the solution vector of a
* finite volume.
*/
template <class PrimaryVariables>
Scalar temperature(const PrimaryVariables &priVars)
{ return this->problem_.temperature(); /* constant temperature */ }
......@@ -162,7 +130,8 @@ private:
const Implementation *asImp_() const
{ return static_cast<const Implementation *> (this); }
Scalar upwindWeight_;
Scalar massUpwindWeight_;
Scalar mobilityUpwindWeight_;
};
}
......
dumux/porousmediumflow/1p/implicit/model.hh
View file @ 005b6d58
......@@ -79,58 +79,30 @@ public:
MultiWriter &writer)
{
typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
typedef Dune::BlockVector<Dune::FieldVector<double, dimWorld> > VectorField;
// typedef Dune::BlockVector<Dune::FieldVector<double, dimWorld> > VectorField;
// create the required scalar fields
unsigned numDofs = this->numDofs();
ScalarField *p = writer.allocateManagedBuffer(numDofs);
VectorField *velocity = writer.template allocateManagedBuffer<double, dimWorld>(numDofs);
ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
if (velocityOutput.enableOutput())
{
// initialize velocity field
for (unsigned int i = 0; i < numDofs; ++i)
{
(*velocity)[i] = double(0);
}
}
unsigned numElements = this->gridView_().size(0);
ScalarField *rank = writer.allocateManagedBuffer(numElements);
for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
{
int eIdx = this->problem_().model().elementMapper().index(element);
int eIdx = this->elementMapper().index(element);
(*rank)[eIdx] = this->gridView_().comm().rank();
FVElementGeometry fvGeometry;
fvGeometry.update(this->gridView_(), element);
ElementVolumeVariables elemVolVars;
elemVolVars.update(this->problem_(),
element,
fvGeometry,
false /* oldSol? */);
for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
for (auto&& scv : this->fvGeometries(element).scvs())
{
int dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
(*p)[dofIdxGlobal] = elemVolVars[scvIdx].pressure();
}
const auto& spatialParams = this->problem_().spatialParams();
auto dofIdxGlobal = scv.dofIndex();
// velocity output
if (velocityOutput.enableOutput())
{
velocityOutput.calculateVelocity(*velocity, elemVolVars, fvGeometry, element, /*phaseIdx=*/0);
(*p)[dofIdxGlobal] = this->curVolVars(scv).pressure();
}
}
writer.attachDofData(*p, "p", isBox);
if (velocityOutput.enableOutput())
{
writer.attachDofData(*velocity, "velocity", isBox, dim);
}
writer.attachCellData(*rank, "process rank");
}
};
......
dumux/porousmediumflow/1p/implicit/propertydefaults.hh
View file @ 005b6d58
......@@ -40,7 +40,6 @@
#include <dumux/material/components/nullcomponent.hh>
#include <dumux/material/fluidsystems/1p.hh>
#include <dumux/material/spatialparams/implicit1p.hh>
#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
namespace Dumux
......@@ -63,9 +62,6 @@ SET_TYPE_PROP(OneP, Model, OnePModel<TypeTag>);
//! the VolumeVariables property
SET_TYPE_PROP(OneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
//! the FluxVariables property
SET_TYPE_PROP(OneP, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//! The indices required by the isothermal single-phase model
SET_TYPE_PROP(OneP, Indices, OnePIndices);
......@@ -133,9 +129,6 @@ SET_PROP(OnePNI, ThermalConductivityModel)
// set isothermal Model
SET_TYPE_PROP(OnePNI, IsothermalModel, OnePModel<TypeTag>);
//set isothermal FluxVariables
SET_TYPE_PROP(OnePNI, IsothermalFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
//set isothermal VolumeVariables
SET_TYPE_PROP(OnePNI, IsothermalVolumeVariables, OnePVolumeVariables<TypeTag>);
......
dumux/porousmediumflow/1p/implicit/volumevariables.hh
View file @ 005b6d58
......@@ -47,6 +47,7 @@ class OnePVolumeVariables : public ImplicitVolumeVariables<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
......@@ -64,20 +65,16 @@ public:
void update(const PrimaryVariables &priVars,
const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx,
const bool isOldSol)
const SubControlVolume& scv)
{
ParentType::update(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
ParentType::update(priVars, problem, element, scv);
completeFluidState(priVars, problem, element, fvGeometry, scvIdx, fluidState_);
completeFluidState(priVars, problem, element, scv, fluidState_);
// porosity
porosity_ = problem.spatialParams().porosity(element,
fvGeometry,
scvIdx);
porosity_ = problem.spatialParams().porosity(scv);
// energy related quantities not contained in the fluid state
asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
asImp_().updateEnergy_(priVars, problem, element, scv);
};
/*!
......@@ -86,12 +83,10 @@ public:
static void completeFluidState(const PrimaryVariables& priVars,
const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const int scvIdx,
const SubControlVolume& scv,
FluidState& fluidState)
{
Scalar t = Implementation::temperature_(priVars, problem, element,
fvGeometry, scvIdx);
Scalar t = Implementation::temperature_(priVars, problem, element, scv);
fluidState.setTemperature(t);
fluidState.setSaturation(/*phaseIdx=*/0, 1.);
......@@ -173,12 +168,11 @@ public:
protected:
static Scalar temperature_(const PrimaryVariables &priVars,
const Problem& problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx)
const Problem& problem,
const Element &element,
const SubControlVolume &scv)
{
return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
return problem.temperatureAtPos(scv.dofPosition());
}
template<class ParameterCache>
......@@ -195,9 +189,7 @@ protected:
void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx,
const bool isOldSol)
const SubControlVolume& scv)
{ }
FluidState fluidState_;
......
test/porousmediumflow/1p/implicit/1ptestproblem.hh
View file @ 005b6d58
......@@ -142,6 +142,7 @@ class OnePTestProblem : public ImplicitPorousMediaProblem<TypeTag>
typedef typename GridView::Intersection Intersection;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
......@@ -165,9 +166,9 @@ public:
*
* This is used as a prefix for files generated by the simulation.
*/
const char *name() const
std::string name() const
{
return name_.c_str();
return name_;
}
/*!
......@@ -184,10 +185,9 @@ public:
* \param values Stores the source values, acts as return value
* \param globalPos The global position
*/
void sourceAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
{
values = 0;
return PrimaryVariables(0);
}
// \}
/*!
......@@ -202,14 +202,17 @@ public:
* \param values The boundary types for the conservation equations
* \param globalPos The position of the center of the finite volume
*/
void boundaryTypesAtPos(BoundaryTypes &values,
const GlobalPosition &globalPos) const
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
double eps = 1.0e-6;
BoundaryTypes values;
Scalar eps = 1.0e-6;
if (globalPos[dimWorld-1] < eps || globalPos[dimWorld-1] > this->bBoxMax()[dimWorld-1] - eps)
values.setAllDirichlet();
else
values.setAllNeumann();
return values;
}
/*!
......@@ -221,10 +224,11 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
void dirichletAtPos(PrimaryVariables &values,
const GlobalPosition &globalPos) const
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0);
values[pressureIdx] = 1.0e+5*(2.0 - globalPos[dimWorld-1]);
return values;
}
/*!
......@@ -235,15 +239,9 @@ public:
* in normal direction of each component. Negative values mean
* influx.
*/
using ParentType::neumann;
void neumann(PrimaryVariables &priVars,
const Element &element,
const FVElementGeometry &fvGeometry,
const Intersection &intersection,
const int scvIdx,
const int boundaryFaceIdx) const
PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
{
priVars[conti0EqIdx] = 0;
return PrimaryVariables(0);
}
// \}
......@@ -259,12 +257,11 @@ public:
* For this method, the \a priVars parameter stores primary
* variables.
*/
void initial(PrimaryVariables &priVars,
const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
PrimaryVariables initial(const SubControlVolume& scv) const
{
PrimaryVariables priVars(0);
priVars[pressureIdx] = 1.0e+5;
return priVars;
}
// \}
......
test/porousmediumflow/1p/implicit/1ptestspatialparams.hh
View file @ 005b6d58
......@@ -59,9 +59,11 @@ class OnePTestSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef std::vector<Scalar> ScalarVector;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::IndexSet IndexSet;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
enum {
dim=GridView::dimension,
......@@ -97,13 +99,9 @@ public:
* \param scvIdx The index sub-control volume face where the
* intrinsic velocity ought to be calculated.
*/
Scalar intrinsicPermeability(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
Scalar intrinsicPermeability(const SubControlVolume &scv) const
{
const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
if (isInLens_(globalPos))
if (isInLens_(scv.dofPosition()))
{
if(randomField_)
return randomPermeability_[indexSet_.index(element.template subEntity<dim> (0))];
......@@ -120,9 +118,7 @@ public:
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume where
*/
Scalar porosity(const Element &element,
const FVElementGeometry &fvGeometry,
const int scvIdx) const
Scalar porosity(const SubControlVolume &scv) const
{ return 0.4; }
/*!
......@@ -160,7 +156,7 @@ private:
bool isInLens_(const GlobalPosition &globalPos) const
{
for (int i = 0; i < dimWorld; ++i) {
if (globalPos[i] < lensLowerLeft_[i] || globalPos[i] > lensUpperRight_[i])
if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
return false;
}
return true;
......@@ -174,7 +170,9 @@ private:
ScalarVector randomPermeability_;
const IndexSet& indexSet_;
static const Scalar eps_ = 1.5e-7;
};
} // end namespace
#endif
#endif
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