diff --git a/dumux/implicit/cellcentered/stencils.hh b/dumux/implicit/cellcentered/stencils.hh new file mode 100644 index 0000000000000000000000000000000000000000..e4fcc3e503815bbcbfa0202d9229fda439ecd0bc --- /dev/null +++ b/dumux/implicit/cellcentered/stencils.hh @@ -0,0 +1,76 @@ +// -*- 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 Implements the notion of stencils + */ +#ifndef DUMUX_CC_STENCILS_HH +#define DUMUX_CC_STENCILS_HH + +#include <set> +#include <dumux/implicit/cellcentered/properties.hh> + +namespace Dumux +{ +/*! + * \ingroup ImplicitModel + * \brief Base class for a sub control volume, i.e a part of the control + * volume we are making the balance for. + */ +template<class TypeTag> +class CCStencils +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using GridView = typename GET_PROP_TYPE(TypeTag, GridView); + using IndexType = typename GridView::IndexSet::IndexType; + using Element = typename GridView::template Codim<0>::Entity; + // TODO a separate stencil class all stencils can derive from? + using Stencil = std::set<IndexType>; +public: + void update(const Problem& problem, const Element& element) + { + for (auto&& scvf : problem.model().fvGeometries(element).scvfs()) + { + auto&& fluxStencil = problem.model().fluxVars(scvf).stencil(); + elementStencil_.insert(fluxStencil.begin(), fluxStencil.end()); + } + neighborStencil_ = elementStencil_; + neighborStencil_.erase(problem.elementMapper().index(element)); + } + + //! The full element stencil (all element this element is interacting with) + const Stencil& elementStencil() const + { + return elementStencil_; + } + + //! //! The full element stencil without this element + const Stencil& neighborStencil() const + { + return neighborStencil_; + } + +private: + Stencil elementStencil_; + Stencil neighborStencil_; +}; + +} // end namespace + +#endif diff --git a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh index 138d21dc449098c1b4910d62abeb0cc963b0a8de..05e7edbdead719781a250f758ed1596f1324133f 100644 --- a/dumux/implicit/cellcentered/tpfa/propertydefaults.hh +++ b/dumux/implicit/cellcentered/tpfa/propertydefaults.hh @@ -28,10 +28,10 @@ #define DUMUX_CCTPFA_PROPERTY_DEFAULTS_HH #include <dumux/implicit/propertydefaults.hh> -#include <dumux/implicit/cellcentered/assembler.hh> #include <dumux/implicit/fvelementgeometry.hh> #include <dumux/implicit/cellcentered/tpfa/fvelementgeometryvector.hh> #include <dumux/implicit/cellcentered/elementboundarytypes.hh> +#include <dumux/implicit/cellcentered/stencils.hh> #include <dumux/implicit/cellcentered/localresidual.hh> #include <dumux/implicit/cellcentered/properties.hh> #include <dumux/porousmediumflow/implicit/cellcentered/tpfa/darcyfluxvariables.hh> @@ -55,11 +55,8 @@ SET_TYPE_PROP(CCTpfaModel, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMap //! Set the BaseLocalResidual to CCLocalResidual SET_TYPE_PROP(CCTpfaModel, BaseLocalResidual, Dumux::CCLocalResidual<TypeTag>); -//! An array of secondary variable containers -SET_TYPE_PROP(CCTpfaModel, ElementVolumeVariables, Dumux::CCElementVolumeVariables<TypeTag>); - -//! Assembler for the global jacobian matrix -SET_TYPE_PROP(CCTpfaModel, JacobianAssembler, Dumux::CCAssembler<TypeTag>); +//! The stencil container +SET_TYPE_PROP(CCTpfaModel, Stencils, Dumux::CCStencils<TypeTag>); //! indicate that this is no box discretization SET_BOOL_PROP(CCTpfaModel, ImplicitIsBox, false); diff --git a/dumux/implicit/model.hh b/dumux/implicit/model.hh index 397a34f28dff9415d51a5b9b49aadeaa295b0acc..7d0ab2430c162960a7290ab386f630b354c1d702 100644 --- a/dumux/implicit/model.hh +++ b/dumux/implicit/model.hh @@ -51,13 +51,14 @@ class ImplicitModel typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector; typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; typedef typename GET_PROP_TYPE(TypeTag, JacobianAssembler) JacobianAssembler; - typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariablesVector) VolumeVariablesVector; typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume; typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace; typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables; typedef typename GET_PROP_TYPE(TypeTag, FluxVariablesVector) FluxVariablesVector; + typedef typename GET_PROP_TYPE(TypeTag, Stencils) Stencils; + typedef typename GET_PROP_TYPE(TypeTag, StencilsVector) StencilsVector; enum { numEq = GET_PROP_VALUE(TypeTag, NumEq), @@ -109,14 +110,21 @@ public: if (isBox) boxVolume_.resize(numDofs); - localJacobian_.init(problem_()); - jacAsm_ = std::make_shared<JacobianAssembler>(); - jacAsm_->init(problem_()); - asImp_().applyInitialSolution_(); // update the volVars with the initial solution - volVarsVector_.update(problem_(), curSol()); + curVolVarsVector_.update(problem_(), curSol()); + + // update the flux vars (precompute transmissibilities) + fluxVarsVector_.update(problem_()); + + // update stencils + stencilsVector_.update(problem_()); + + // initialize assembler and create matrix + localJacobian_.init(problem_()); + jacAsm_ = std::make_shared<JacobianAssembler>(); + jacAsm_->init(problem_()); // also set the solution of the "previous" time step to the // initial solution. @@ -938,6 +946,9 @@ protected: // the flux variables vector FluxVariablesVector fluxVarsVector_; + // the stencils vector + StencilsVector stencilsVector_; + // the finite volume element geometries std::shared_ptr<FVElementGeometryVector> fvGeometries_; diff --git a/dumux/implicit/properties.hh b/dumux/implicit/properties.hh index 1d68a8f0cef87df51b52795948492bd09336a073..184177e1c2286a48da8ab8b9c66fd02bd98a3a81 100644 --- a/dumux/implicit/properties.hh +++ b/dumux/implicit/properties.hh @@ -86,6 +86,10 @@ NEW_PROP_TAG(DiffusionFluxVariables); //! The type that handles the calculation NEW_PROP_TAG(EnergyFluxes); //! specifies if the model consideres heat transport phenomena NEW_PROP_TAG(EnergyFluxVariables); //! The type for the heat flux calculation +// stencils +NEW_PROP_TAG(Stencils); //! The stencils container type +NEW_PROP_TAG(StencilsVector); //! The type of the global vector of stencils per element + // high level simulation control NEW_PROP_TAG(TimeManager); //!< Manages the simulation time NEW_PROP_TAG(NewtonMethod); //!< The type of the newton method diff --git a/dumux/implicit/propertydefaults.hh b/dumux/implicit/propertydefaults.hh index 6a605430bb73c8dbad6fd80b5c7395919dd77c04..f8647f59ea3e56b16c6eeb4131f561ea33ca63b0 100644 --- a/dumux/implicit/propertydefaults.hh +++ b/dumux/implicit/propertydefaults.hh @@ -37,12 +37,14 @@ #include "properties.hh" #include "model.hh" +#include "assembler.hh" #include "localjacobian.hh" #include "volumevariables.hh" #include "volumevariablesvector.hh" -#include "fluxvariables" +#include "fluxvariables.hh" #include "fluxvariablesvector.hh" #include "fvelementgeometry.hh" +#include "stencilsvector.hh" namespace Dumux { @@ -97,6 +99,9 @@ SET_TYPE_PROP(ImplicitBase, VolumeVariables, ImplicitVolumeVariables<TypeTag>); //! The global volume variables vector class SET_TYPE_PROP(ImplicitBase, VolumeVariablesVector, Dumux::VolumeVariablesVector<TypeTag>); +//! The global volume variables vector class +SET_TYPE_PROP(ImplicitBase, StencilsVector, Dumux::StencilsVector<TypeTag>); + //! Set darcy fluxes as the default fluxes to be considered SET_BOOL_PROP(ImplicitBase, DarcyFluxes, true); diff --git a/dumux/implicit/stencilsvector.hh b/dumux/implicit/stencilsvector.hh new file mode 100644 index 0000000000000000000000000000000000000000..83768e0d58cc80318effe1e1f5f7f7df2cedc3c6 --- /dev/null +++ b/dumux/implicit/stencilsvector.hh @@ -0,0 +1,54 @@ +// -*- 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 Implements the notion of stencils + */ +#ifndef DUMUX_STENCILSVECTOR_HH +#define DUMUX_STENCILSVECTOR_HH + +#include <dumux/implicit/properties.hh> + +namespace Dumux +{ +/*! + * \ingroup ImplicitModel + * \brief Base class for a sub control volume, i.e a part of the control + * volume we are making the balance for. + */ +template<class TypeTag> +class StencilsVector : public std::vector<typename GET_PROP_TYPE(TypeTag, Stencils)> +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + using Stencils = typename GET_PROP_TYPE(TypeTag, Stencils); +public: + void update(const Problem& problem) + { + this->resize(problem.gridView().size(0)); + for (const auto& element : elements(problem.gridView())) + { + auto eIdx = problem.elementMapper().index(element); + (*this)[eIdx].update(problem, element); + } + } +}; + +} // end namespace + +#endif