diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh index 5045cea0455864883214698db93df741d890ba5f..0d03535433eb016d82deaff4df8e19ad5bad990a 100644 --- a/dumux/freeflow/staggered/localresidual.hh +++ b/dumux/freeflow/staggered/localresidual.hh @@ -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, diff --git a/dumux/freeflow/staggered/problem.hh b/dumux/freeflow/staggered/problem.hh index cdd81818ef5b74dba3b66539aaf5790ac4898a04..5ea3896999d3f4f99743458a65a02ea0ab1ff3c5 100644 --- a/dumux/freeflow/staggered/problem.hh +++ b/dumux/freeflow/staggered/problem.hh @@ -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; diff --git a/dumux/freeflow/staggered/volumevariables.hh b/dumux/freeflow/staggered/volumevariables.hh index 497c97fbab9d34d705854184da13e51ea3502c60..39b84dd9956ccbab41807b79695cc134ba3a1087 100644 --- a/dumux/freeflow/staggered/volumevariables.hh +++ b/dumux/freeflow/staggered/volumevariables.hh @@ -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, diff --git a/dumux/implicit/staggered/model.hh b/dumux/implicit/staggered/model.hh index 0adec2a091a2974bb1a59146c4d75b830d359c47..5bf97ce7f67ea940a3235999a0c13c4dfb94072e 100644 --- a/dumux/implicit/staggered/model.hh +++ b/dumux/implicit/staggered/model.hh @@ -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; + } } } diff --git a/dumux/implicit/staggered/problem.hh b/dumux/implicit/staggered/problem.hh new file mode 100644 index 0000000000000000000000000000000000000000..9ab98e22919626d3ddf4141dfa9462aa6b76b7da --- /dev/null +++ b/dumux/implicit/staggered/problem.hh @@ -0,0 +1,341 @@ +// -*- 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 diff --git a/dumux/implicit/staggered/propertydefaults.hh b/dumux/implicit/staggered/propertydefaults.hh index 4af4ce7198d6e458a9bc4b5375e23884a2c1ee46..be6dc795d346eb11a64c3db76678038ce3578a4f 100644 --- a/dumux/implicit/staggered/propertydefaults.hh +++ b/dumux/implicit/staggered/propertydefaults.hh @@ -63,6 +63,7 @@ #include "properties.hh" #include "newtoncontroller.hh" #include "model.hh" +#include "problem.hh" namespace Dumux { diff --git a/test/freeflow/staggered/staggeredtestproblem.hh b/test/freeflow/staggered/staggeredtestproblem.hh index 56d8789ed0a5fe5a79db69302d47da03cba18152..8eaf56a13cecb4a0ba9cee0db7668b7fee7c7bde 100644 --- a/test/freeflow/staggered/staggeredtestproblem.hh +++ b/test/freeflow/staggered/staggeredtestproblem.hh @@ -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: