Commit 9d58d440 authored by Kilian Weishaupt's avatar Kilian Weishaupt Committed by Timo Koch
Browse files

[staggeredGrid] Account for the possibility of having different kinds of priVars

* Cell-center privars might be related to two or even more equations
  while face priVars will probably always be related to only one equation
parent 922fc200
......@@ -100,17 +100,25 @@ public:
const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndexSelf()).velocity();
CellCenterPrimaryVariables flux(0.0);
const Scalar eps = 1e-6;
for(int direction = 0; direction < dim; ++direction)
{
if(scvf.unitOuterNormal()[direction] > eps && velocity > 0) // positive x-direction
return insideVolVars.density(0) * velocity;
if(scvf.unitOuterNormal()[direction] > eps && velocity > 0)
{
flux[0] = insideVolVars.density(0) * velocity;
return flux;
}
if(scvf.unitOuterNormal()[direction] < eps && velocity < 0)
return outsideVolVars.density(0) * velocity;
{
flux[0] = outsideVolVars.density(0) * velocity;
return flux;
}
}
return CellCenterPrimaryVariables(0.0);
return flux;
}
static CellCenterPrimaryVariables computeSourceForCellCenter(const Element &element,
......
......@@ -23,7 +23,7 @@
#ifndef DUMUX_NAVIERSTOKES_PROBLEM_HH
#define DUMUX_NAVIERSTOKES_PROBLEM_HH
#include <dumux/implicit/problem.hh>
#include <dumux/implicit/staggered/problem.hh>
#include "properties.hh"
......@@ -37,9 +37,9 @@ namespace Dumux
* This implements gravity (if desired) and a function returning the temperature.
*/
template<class TypeTag>
class NavierStokesProblem : public ImplicitProblem<TypeTag>
class NavierStokesProblem : public StaggeredProblem<TypeTag>
{
typedef ImplicitProblem<TypeTag> ParentType;
typedef StaggeredProblem<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
......
......@@ -47,7 +47,7 @@ class NavierStokesVolumeVariables : public StaggeredVolumeVariables<TypeTag>
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
......@@ -60,7 +60,7 @@ public:
/*!
* \copydoc ImplicitVolumeVariables::update
*/
void update(const PrimaryVariables &priVars,
void update(const CellCenterPrimaryVariables &priVars,
const Problem &problem,
const Element &element,
const SubControlVolume& scv)
......@@ -73,7 +73,7 @@ public:
/*!
* \copydoc ImplicitModel::completeFluidState
*/
static void completeFluidState(const PrimaryVariables& priVars,
static void completeFluidState(const CellCenterPrimaryVariables& priVars,
const Problem& problem,
const Element& element,
const SubControlVolume& scv,
......
......@@ -65,8 +65,12 @@ class StaggeredBaseModel : public ImplicitModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum { dim = GridView::dimension };
enum { dimWorld = GridView::dimensionworld };
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
......@@ -181,6 +185,28 @@ public:
}
/*!
* \brief Returns the maximum relative shift between two vectors of
* primary variables.
*
* \param priVars1 The first vector of primary variables
* \param priVars2 The second vector of primary variables
*/
template<class CellCenterOrFacePriVars>
Scalar relativeShiftAtDof(const CellCenterOrFacePriVars &priVars1,
const CellCenterOrFacePriVars &priVars2)
{
const auto numEq = priVars1.size();
Scalar result = 0.0;
for (int j = 0; j < numEq; ++j) {
Scalar eqErr = std::abs(priVars1[j] - priVars2[j]);
eqErr /= std::max<Scalar>(1.0, std::abs(priVars1[j] + priVars2[j])/2);
result = std::max(result, eqErr);
}
return result;
}
/*!
* \brief Called by the update() method if it was
* unsuccessful. This is primarily a hook which the actual
......@@ -235,6 +261,13 @@ public:
auto dofIdxGlobal = scv.dofIndex();
this->uCur_[cellCenterIdx][dofIdxGlobal] += initPriVars;
}
// loop over faces
for(auto&& scvf : scvfs(fvGeometry))
{
auto initPriVars = this->problem_().initialFaceValues(scvf);
this->uCur_[faceIdx][scvf.dofIndexSelf()] = initPriVars;
}
}
}
......
// -*- 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 Base class for all fully implicit problems
*/
#ifndef DUMUX_IMPLICIT_STAGGERED_PROBLEM_HH
#define DUMUX_IMPLICIT_STAGGERED_PROBLEM_HH
#include <dumux/implicit/problem.hh>
namespace Dumux
{
/*!
* \ingroup ImplicitBaseProblems
* \brief Base class for all fully implicit problems
*
* \note All quantities are specified assuming a threedimensional
* world. Problems discretized using 2D grids are assumed to be
* extruded by \f$1 m\f$ and 1D grids are assumed to have a
* cross section of \f$1m \times 1m\f$.
*/
template<class TypeTag>
class StaggeredProblem : public ImplicitProblem<TypeTag>
{
private:
using ParentType = ImplicitProblem<TypeTag>;
using Implementation = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
using VtkMultiWriter = typename GET_PROP_TYPE(TypeTag, VtkMultiWriter);
using NewtonMethod = typename GET_PROP_TYPE(TypeTag, NewtonMethod);
using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
using Model = typename GET_PROP_TYPE(TypeTag, Model);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
using ElementMapper = typename GET_PROP_TYPE(TypeTag, ElementMapper);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using PointSource = typename GET_PROP_TYPE(TypeTag, PointSource);
using PointSourceHelper = typename GET_PROP_TYPE(TypeTag, PointSourceHelper);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
using Element = typename GridView::template Codim<0>::Entity;
using Vertex = typename GridView::template Codim<dim>::Entity;
using Intersection = typename GridView::Intersection;
using CoordScalar = typename GridView::Grid::ctype;
using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
enum { dofCodim = isBox ? dim : 0 };
enum { adaptiveGrid = GET_PROP_VALUE(TypeTag, AdaptiveGrid) };
using GridAdaptModel = ImplicitGridAdapt<TypeTag, adaptiveGrid>;
using BoundingBoxTree = Dumux::BoundingBoxTree<GridView>;
// copying a problem is not a good idea
StaggeredProblem(const StaggeredProblem &);
public:
StaggeredProblem(TimeManager &timeManager, const GridView &gridView) : ParentType(timeManager, gridView)
{}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* control volume.
*
* \param values The dirichlet values for the primary variables
* \param scvFace the sub control volume face
*
* The method returns the boundary types information.
*/
CellCenterPrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
{
// forward it to the method which only takes the global coordinate
if (isBox)
{
DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for box method.");
}
else
return asImp_().dirichletAtPos(scvf.center());
}
CellCenterPrimaryVariables dirichlet(const Element &element, const SubControlVolume &scv) const
{
// forward it to the method which only takes the global coordinate
if (!isBox)
{
DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for cell-centered method.");
}
else
return asImp_().dirichletAtPos(scv.dofPosition());
}
/*!
* \brief Evaluate the boundary conditions for a dirichlet
* control volume.
*
* \param globalPos The position of the center of the finite volume
* for which the dirichlet condition ought to be
* set in global coordinates
*
* For this method, the \a values parameter stores primary variables.
*/
CellCenterPrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
// Throw an exception (there is no reasonable default value
// for Dirichlet conditions)
DUNE_THROW(Dune::InvalidStateException,
"The problem specifies that some boundary "
"segments are dirichlet, but does not provide "
"a dirichlet() method.");
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* This is the method for the case where the Neumann condition is
* potentially solution dependent and requires some quantities that
* are specific to the fully-implicit method.
*
* \param values The neumann values for the conservation equations in units of
* \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param intersection The intersection between element and boundary
* \param scvIdx The local subcontrolvolume index
* \param boundaryFaceIdx The index of the boundary face
* \param elemVolVars All volume variables for the element
*
* For this method, the \a values parameter stores the flux
* in normal direction of each phase. Negative values mean influx.
* E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
CellCenterPrimaryVariables neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolvars,
const SubControlVolumeFace& scvFace) const
{
// forward it to the interface with only the global position
return asImp_().neumannAtPos(scvFace.center());
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
* \param values The neumann values for the conservation equations in units of
* \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
* \param globalPos The position of the boundary face's integration point in global coordinates
*
* For this method, the \a values parameter stores the flux
* in normal direction of each phase. Negative values mean influx.
* E.g. for the mass balance that would be the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
CellCenterPrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
{
// Throw an exception (there is no reasonable default value
// for Neumann conditions)
DUNE_THROW(Dune::InvalidStateException,
"The problem specifies that some boundary "
"segments are neumann, but does not provide "
"a neumannAtPos() method.");
}
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* This is the method for the case where the source term is
* potentially solution dependent and requires some quantities that
* are specific to the fully-implicit method.
*
* \param values The source and sink values for the conservation equations in units of
* \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param elemVolVars All volume variables for the element
* \param scv The subcontrolvolume
*
* For this method, the \a values parameter stores the conserved quantity rate
* generated or annihilate per volume unit. Positive values mean
* that the conserved quantity is created, negative ones mean that it vanishes.
* E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
*/
CellCenterPrimaryVariables source(const Element &element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
// forward to solution independent, fully-implicit specific interface
return asImp_().sourceAtPos(scv.dofPosition());
}
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* \param values The source and sink values for the conservation equations in units of
* \f$ [ \textnormal{unit of conserved quantity} / (m^3 \cdot s )] \f$
* \param globalPos The position of the center of the finite volume
* for which the source term ought to be
* specified in global coordinates
*
* For this method, the \a values parameter stores the conserved quantity rate
* generated or annihilate per volume unit. Positive values mean
* that the conserved quantity is created, negative ones mean that it vanishes.
* E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
*/
CellCenterPrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
{
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide "
"a sourceAtPos() method.");
}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param values The initial values for the primary variables
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scvIdx The local subcontrolvolume index
*
* For this method, the \a values parameter stores primary
* variables.
*/
CellCenterPrimaryVariables initial(const SubControlVolume &scv) const
{
// forward to generic interface
return asImp_().initialAtPos(scv.dofPosition());
}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param values The initial values for the primary variables
* \param globalPos The position of the center of the finite volume
* for which the initial values ought to be
* set (in global coordinates)
*
* For this method, the \a values parameter stores primary variables.
*/
CellCenterPrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
// Throw an exception (there is no reasonable default value
// for initial values)
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide "
"a initialAtPos() method.");
}
/*!
* \brief Evaluate the initial value for a facet.
*
* \param values The initial values for the primary variables
* \param element The finite element
* \param fvGeometry The finite-volume geometry
* \param scvIdx The local subcontrolvolume index
*
* For this method, the \a values parameter stores primary
* variables.
*/
FacePrimaryVariables initialFaceValues(const SubControlVolumeFace &scvf) const
{
// forward to generic interface
return asImp_().initialFaceValueAtPos(scvf.center());
}
/*!
* \brief Evaluate the initial value for a facet.
*
* \param values The initial values for the primary variables
* \param globalPos The position of the center of the finite volume
* for which the initial values ought to be
* set (in global coordinates)
*
* For this method, the \a values parameter stores primary variables.
*/
FacePrimaryVariables initialFaceValueAtPos(const GlobalPosition &globalPos) const
{
// Throw an exception (there is no reasonable default value
// for initial values)
DUNE_THROW(Dune::InvalidStateException,
"The problem does not provide "
"a initialAtPos() method.");
}
protected:
//! Returns the implementation of the problem (i.e. static polymorphism)
Implementation &asImp_()
{ return *static_cast<Implementation *>(this); }
//! \copydoc asImp_()
const Implementation &asImp_() const
{ return *static_cast<const Implementation *>(this); }
};
} // namespace Dumux
#include <dumux/implicit/adaptive/gridadaptpropertydefaults.hh>
#endif
......@@ -63,6 +63,7 @@
#include "properties.hh"
#include "newtoncontroller.hh"
#include "model.hh"
#include "problem.hh"
namespace Dumux {
......
......@@ -132,6 +132,9 @@ class StaggeredTestProblem : public NavierStokesProblem<TypeTag>
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
public:
StaggeredTestProblem(TimeManager &timeManager, const GridView &gridView)
: ParentType(timeManager, gridView)
......@@ -171,9 +174,9 @@ public:
* \param values Stores the source values, acts as return value
* \param globalPos The global position
*/
PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
CellCenterPrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
{
return PrimaryVariables(0);
return CellCenterPrimaryVariables(0);
}
// \}
/*!
......@@ -210,9 +213,9 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
CellCenterPrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0);
CellCenterPrimaryVariables values(0);
values[pressureIdx] = 1.0e+5*(2.0 - globalPos[dimWorld-1]);
return values;
}
......@@ -225,9 +228,9 @@ public:
* in normal direction of each component. Negative values mean
* influx.
*/
PrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
CellCenterPrimaryVariables neumannAtPos(const GlobalPosition& globalPos) const
{
return PrimaryVariables(0);
return CellCenterPrimaryVariables(0);
}
// \}
......@@ -243,13 +246,33 @@ public:
* For this method, the \a priVars parameter stores primary
* variables.
*/
PrimaryVariables initial(const SubControlVolume& scv) const
CellCenterPrimaryVariables initial(const SubControlVolume& scv) const
{
PrimaryVariables priVars(0);
CellCenterPrimaryVariables priVars(0);
priVars[pressureIdx] = 1.0e+5;
return priVars;
}
/*!
* \brief Evaluate the initial value for a facet.
*
* \param values The initial values for the primary variables
* \param globalPos The position of the center of the finite volume
* for which the initial values ought to be
* set (in global coordinates)
*
* For this method, the \a values parameter stores primary variables.
*/
FacePrimaryVariables initialFaceValueAtPos(const GlobalPosition &globalPos) const
{
FacePrimaryVariables value(0.0);
if(globalPos[0] < 1e-6)
value[0] = 1.0;
return value;
}
// \}
private:
......
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