From 642be0796eac930d4a8f34302ad4a32196e891e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 8 Mar 2021 16:00:55 +0100 Subject: [PATCH 01/32] [disc] introduce local context --- dumux/discretization/CMakeLists.txt | 1 + dumux/discretization/localcontext.hh | 82 ++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) create mode 100644 dumux/discretization/localcontext.hh diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt index 3b145a92fb..625fe8b3cc 100644 --- a/dumux/discretization/CMakeLists.txt +++ b/dumux/discretization/CMakeLists.txt @@ -20,6 +20,7 @@ functionspacebasis.hh fvgridvariables.hh fvproperties.hh gridvariables.hh +localcontext.hh localview.hh method.hh rotationpolicy.hh diff --git a/dumux/discretization/localcontext.hh b/dumux/discretization/localcontext.hh new file mode 100644 index 0000000000..f79742702c --- /dev/null +++ b/dumux/discretization/localcontext.hh @@ -0,0 +1,82 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Discretization + * \brief Class that contains the element-local (or element stencil-local) data + * required to evaluate the terms of discrete equations. + */ +#ifndef DUMUX_LOCAL_CONTEXT_HH +#define DUMUX_LOCAL_CONTEXT_HH + +namespace Dumux { +namespace Experimental { + +/*! + * \ingroup Discretization + * \brief Implementation of element-stencil-local contexts, which solely store + * the local geometry and primary/secondary variables. This implementation + * defines the minimum interface required for contexts in single-domain + * applications to work with the default assembly mechanism. + * \tparam EV The element-local view on the grid variables + */ +template +class DefaultLocalContext +{ + using GridVariables = typename EV::GridVariables; + using GridGeometry = typename GridVariables::GridGeometry; + +public: + using ElementGridGeometry = typename GridGeometry::LocalView; + using ElementVariables = EV; + + //! Constructor + DefaultLocalContext(const ElementGridGeometry& eg, + const ElementVariables& ev) + : elemGeom_(&eg) + , elemVars_(&ev) + {} + + //! Return the element-local view on the grid geometry + const ElementGridGeometry elementGridGeometry() const + { return *elemGeom_; } + + //! Return the element-local view on the grid variables + const ElementVariables& elementVariables() const + { return *elemVars_; } + +private: + const ElementGridGeometry* elemGeom_; + const ElementVariables* elemVars_; +}; + +/*! + * \ingroup Discretization + * \brief Factory function to create a context from local views. + */ +template +DefaultLocalContext +makeLocalContext(const typename EV::GridVariables::GridGeometry::LocalView& gglocalView, + const EV& elemVariables) +{ return {gglocalView, elemVariables}; } + +} // end namespace Experimental +} // end namespace Dumux + +#endif -- GitLab From b4259840ae06f5fc4a458354f812547c8cb45efa Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:02:20 +0200 Subject: [PATCH 02/32] [assembly] introduce fv operators --- dumux/assembly/fv/CMakeLists.txt | 3 + dumux/assembly/fv/operators.hh | 133 +++++++++++++++++++++++++++++++ 2 files changed, 136 insertions(+) create mode 100644 dumux/assembly/fv/CMakeLists.txt create mode 100644 dumux/assembly/fv/operators.hh diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt new file mode 100644 index 0000000000..72883ac6ae --- /dev/null +++ b/dumux/assembly/fv/CMakeLists.txt @@ -0,0 +1,3 @@ +install(FILES +operators.hh +DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv) diff --git a/dumux/assembly/fv/operators.hh b/dumux/assembly/fv/operators.hh new file mode 100644 index 0000000000..f2103086aa --- /dev/null +++ b/dumux/assembly/fv/operators.hh @@ -0,0 +1,133 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief The base class for the sub-control entity-local evaluation of + * the terms of equations in the context of finite-volume schemes + */ +#ifndef DUMUX_FV_OPERATORS_HH +#define DUMUX_FV_OPERATORS_HH + +#include +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup Assembly + * \brief The base class for the sub-control entity-local evaluation of + * the terms of equations in the context of finite-volume schemes + * \tparam LC the element-stencil-local data required to evaluate the terms + */ +template +class FVOperators +{ + // The grid geometry on which the scheme operates + using FVElementGeometry = typename LC::ElementGridGeometry; + using GridGeometry = typename FVElementGeometry::GridGeometry; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Extrusion = Extrusion_t; + + // The variables (primary/secondary) on the grid + using GridVariables = typename LC::ElementVariables::GridVariables; + using PrimaryVariables = typename GridVariables::PrimaryVariables; + +public: + //! export the local context on which this operates + using LocalContext = LC; + + //! export the type used to store scalar values for all equations + using NumEqVector = Dumux::NumEqVector; + + // export the types of the flux/source/storage terms + using FluxTerm = NumEqVector; + using SourceTerm = NumEqVector; + using StorageTerm = NumEqVector; + + /*! + * \name Model specific interfaces + * \note The following method are the model specific implementations of the + * operators appearing in the equation and should be overloaded by the + * implementations. + */ + // \{ + + /*! + * \brief Compute the storage term of the equations for the given sub-control volume + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scv The sub-control volume + * \note This must be overloaded by the implementation + */ + template + static StorageTerm storage(const Problem& problem, + const LocalContext& context, + const SubControlVolume& scv) + { DUNE_THROW(Dune::NotImplemented, "Storage operator not implemented!"); } + + /*! + * \brief Compute the flux term of the equations for the given sub-control volume face + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scvf The sub-control volume face for which the flux term is to be computed + * \note This must be overloaded by the implementation + */ + template + static FluxTerm flux(const Problem& problem, + const LocalContext& context, + const SubControlVolumeFace& scvf) + { DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!"); } + + /*! + * \brief Compute the source term of the equations for the given sub-control volume + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scv The sub-control volume for which the source term is to be computed + * \note This is a default implementation forwarding to interfaces in the problem + */ + template + static SourceTerm source(const Problem& problem, + const LocalContext& context, + const SubControlVolume& scv) + { + SourceTerm source(0.0); + + // TODO: The problem-interfaces should be adapted to receive context? + const auto& fvGeometry = context.elementGridGeometry(); + const auto& elemVolVars = context.elementVariables().elemVolVars(); + const auto& element = fvGeometry.element(); + + // add contributions from volume flux sources + source += problem.source(element, fvGeometry, elemVolVars, scv); + + // add contribution from possible point sources + source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv); + + // multiply with scv volume + source *= Extrusion::volume(scv)*elemVolVars[scv].extrusionFactor(); + + return source; + } +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 46346c46569c68c2e6a1d7f8ca19e00f74a76bb7 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:07:19 +0200 Subject: [PATCH 03/32] [nonisothermal] add fv-type operators --- .../nonisothermal/CMakeLists.txt | 1 + .../nonisothermal/operators.hh | 154 ++++++++++++++++++ 2 files changed, 155 insertions(+) create mode 100644 dumux/porousmediumflow/nonisothermal/operators.hh diff --git a/dumux/porousmediumflow/nonisothermal/CMakeLists.txt b/dumux/porousmediumflow/nonisothermal/CMakeLists.txt index 24dc71bd45..5db30722a0 100644 --- a/dumux/porousmediumflow/nonisothermal/CMakeLists.txt +++ b/dumux/porousmediumflow/nonisothermal/CMakeLists.txt @@ -3,5 +3,6 @@ indices.hh iofields.hh localresidual.hh model.hh +operators.hh volumevariables.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/nonisothermal) diff --git a/dumux/porousmediumflow/nonisothermal/operators.hh b/dumux/porousmediumflow/nonisothermal/operators.hh new file mode 100644 index 0000000000..1949e41397 --- /dev/null +++ b/dumux/porousmediumflow/nonisothermal/operators.hh @@ -0,0 +1,154 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup NIModel + * \brief Sub-control entity-local evaluation of the operators + * involved in the system of equations non-isothermal models + * that assume thermal equilibrium between all phases. + */ +#ifndef DUMUX_FV_NON_ISOTHERMAL_OPERATORS_HH +#define DUMUX_FV_NON_ISOTHERMAL_OPERATORS_HH + +namespace Dumux::Experimental { + +/*! + * \ingroup NIModel + * \brief Sub-control entity-local evaluation of the operators + * involved in the system of equations of non-isothermal models + * that assume thermal equilibrium between all phases. + * \tparam ModelTraits Model-specific traits. + * \tparam LocalContext Element-local context (geometry & primary/secondary variables) + */ +template +class FVNonIsothermalOperators +{ + // The variables required for the evaluation of the equation + using ElementVariables = typename LocalContext::ElementVariables; + + // The grid geometry on which the scheme operates + using GridGeometry = typename LocalContext::ElementGridGeometry::GridGeometry; + using SubControlVolume = typename GridGeometry::SubControlVolume; + +public: + + /*! + * \brief The energy storage in the fluid phase with index phaseIdx + * + * \param storage The mass of the component within the sub-control volume + * \param scv The sub-control volume + * \param context The element-local context (primary/secondary variables) + * \param phaseIdx The phase index + */ + template + static void addFluidPhaseStorage(NumEqVector& storage, + const SubControlVolume& scv, + const LocalContext& context, + int phaseIdx) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + { + const auto& volVars = context.elementVariables().elemVolVars()[scv]; + storage[ModelTraits::Indices::energyEqIdx] += volVars.porosity() + *volVars.density(phaseIdx) + *volVars.internalEnergy(phaseIdx) + *volVars.saturation(phaseIdx); + } + } + + /*! + * \brief The energy storage in the solid matrix + * + * \param storage The mass of the component within the sub-control volume + * \param scv The sub-control volume + * \param context The element-local context (primary/secondary variables) + */ + template + static void addSolidPhaseStorage(NumEqVector& storage, + const SubControlVolume& scv, + const LocalContext& context) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + { + const auto& volVars = context.elementVariables().elemVolVars()[scv]; + storage[ModelTraits::Indices::energyEqIdx] += volVars.temperature() + *volVars.solidHeatCapacity() + *volVars.solidDensity() + *(1.0 - volVars.porosity()); + } + } + + /*! + * \brief The advective phase energy fluxes + * + * \param flux The flux + * \param fluxVars The flux variables. + * \param phaseIdx The phase index + */ + template + static void addHeatConvectionFlux(NumEqVector& flux, + FluxVariables& fluxVars, + int phaseIdx) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + { + auto upwindTerm = [phaseIdx](const auto& volVars) + { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx)*volVars.enthalpy(phaseIdx); }; + + flux[ModelTraits::Indices::energyEqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm); + } + } + + /*! + * \brief The diffusive energy fluxes + * + * \param flux The flux + * \param fluxVars The flux variables. + */ + template + static void addHeatConductionFlux(NumEqVector& flux, + FluxVariables& fluxVars) + { + // do no-op in case energy balance is not active + if constexpr (ModelTraits::enableEnergyBalance()) + flux[ModelTraits::Indices::energyEqIdx] += fluxVars.heatConductionFlux(); + } + + /*! + * \brief heat transfer between the phases for nonequilibrium models + * + * \param source The source which ought to be simulated + * \param element An element which contains part of the control volume + * \param fvGeometry The finite-volume geometry + * \param context The element-local context (primary/secondary variables) + * \param scv The sub-control volume over which we integrate the source term + */ + template + static void addSourceEnergy(NumEqVector& source, + const LocalContext& context, + const SubControlVolume &scv) + {} +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From ba75f3adabd84a84260872034f5b2c294910e139 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:09:16 +0200 Subject: [PATCH 04/32] [immiscible] add fv-type operators --- .../immiscible/CMakeLists.txt | 1 + .../porousmediumflow/immiscible/operators.hh | 185 ++++++++++++++++++ 2 files changed, 186 insertions(+) create mode 100644 dumux/porousmediumflow/immiscible/operators.hh diff --git a/dumux/porousmediumflow/immiscible/CMakeLists.txt b/dumux/porousmediumflow/immiscible/CMakeLists.txt index 192e7955c6..ccb0268dc6 100644 --- a/dumux/porousmediumflow/immiscible/CMakeLists.txt +++ b/dumux/porousmediumflow/immiscible/CMakeLists.txt @@ -1,3 +1,4 @@ install(FILES localresidual.hh +operators.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/immiscible) diff --git a/dumux/porousmediumflow/immiscible/operators.hh b/dumux/porousmediumflow/immiscible/operators.hh new file mode 100644 index 0000000000..50f2a69fb8 --- /dev/null +++ b/dumux/porousmediumflow/immiscible/operators.hh @@ -0,0 +1,185 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup PorousmediumflowModels + * \brief Sub-control entity-local evaluation of the operators + * within in the systems of equations of n-phase immiscible models. + */ +#ifndef DUMUX_FV_IMMISCIBLE_OPERATORS_HH +#define DUMUX_FV_IMMISCIBLE_OPERATORS_HH + +#include +#include +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup PorousmediumflowModels + * \brief Sub-control entity-local evaluation of the operators + * within in the systems of equations of n-phase immiscible models. + * \tparam ModelTraits defines model-related types and variables (e.g. number of phases) + * \tparam FluxVariables the type that is responsible for computing the individual + * flux contributions, i.e., advective, diffusive, convective... + * \tparam LocalContext the element-stencil-local data required to evaluate the terms + * \tparam EnergyOperators optional template argument, specifying the class that + * handles the operators related to non-isothermal effects. + * The default energy operators are compatible with isothermal + * simulations. + */ +template> +class FVImmiscibleOperators +: public FVOperators +{ + using ParentType = FVOperators; + + // The variables required for the evaluation of the equation + using ElementVariables = typename LocalContext::ElementVariables; + using GridVariables = typename ElementVariables::GridVariables; + using VolumeVariables = typename GridVariables::VolumeVariables; + using ElementVolumeVariables = typename ElementVariables::ElementVolumeVariables; + using ElementFluxVariablesCache = typename ElementVariables::ElementFluxVariablesCache; + + // The grid geometry on which the scheme operates + using GridGeometry = typename GridVariables::GridGeometry; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Extrusion = Extrusion_t; + + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + + using Problem = std::decay_t().gridVolVars().problem())>; + using Indices = typename ModelTraits::Indices; + + static constexpr int numPhases = ModelTraits::numFluidPhases(); + static constexpr int conti0EqIdx = ModelTraits::Indices::conti0EqIdx; + static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();; + +public: + //! export the type used to store scalar values for all equations + using typename ParentType::NumEqVector; + + // export the types of the flux/source/storage terms + using typename ParentType::FluxTerm; + using typename ParentType::SourceTerm; + using typename ParentType::StorageTerm; + + /*! + * \brief Compute the storage term of the equations for the given sub-control volume + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scv The sub-control volume + */ + template + static StorageTerm storage(const Problem& problem, + const LocalContext& context, + const SubControlVolume& scv) + { + const auto& volVars = context.elementVariables().elemVolVars()[scv]; + + StorageTerm storage; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + auto eqIdx = conti0EqIdx + phaseIdx; + storage[eqIdx] = volVars.porosity() + * volVars.density(phaseIdx) + * volVars.saturation(phaseIdx); + + // The energy storage in the fluid phase with index phaseIdx + if constexpr (enableEnergyBalance) + EnergyOperators::fluidPhaseStorage(storage, scv, volVars, phaseIdx); + } + + // The energy storage in the solid matrix + if constexpr (enableEnergyBalance) + EnergyOperators::solidPhaseStorage(storage, scv, volVars); + + // multiply with volume + storage *= Extrusion::volume(scv)*volVars.extrusionFactor(); + + return storage; + } + + /*! + * \brief Compute the flux term of the equations for the given sub-control volume face + * \param problem The problem to be solved (could store additionally required quantities) + * \param context The element-stencil-local required data + * \param scvf The sub-control volume face for which the flux term is to be computed + * \note This must be overloaded by the implementation + */ + template + static FluxTerm flux(const Problem& problem, + const LocalContext& context, + const SubControlVolumeFace& scvf) + { + const auto& fvGeometry = context.elementGridGeometry(); + const auto& elemVolVars = context.elementVariables().elemVolVars(); + const auto& elemFluxVarsCache = context.elementVariables().elemFluxVarsCache(); + + FluxVariables fluxVars; + + // for box, the "inside" element is always the one the grid geom is bound to + if constexpr (GridGeometry::discMethod == DiscretizationMethod::box) + fluxVars.init(problem, fvGeometry.element(), fvGeometry, elemVolVars, scvf, elemFluxVarsCache); + else + { + const auto& gridGeom = fvGeometry.gridGeometry(); + const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); + const auto boundElementIdx = gridGeom.elementMapper().index(fvGeometry.element()); + + if (insideScv.elementIndex() == boundElementIdx) + fluxVars.init(problem, fvGeometry.element(), fvGeometry, elemVolVars, scvf, elemFluxVarsCache); + else + fluxVars.init(problem, gridGeom.element(insideScv.elementIndex()), + fvGeometry, elemVolVars, scvf, elemFluxVarsCache); + } + + + NumEqVector flux; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + { + // the physical quantities for which we perform upwinding + auto upwindTerm = [phaseIdx](const auto& volVars) + { return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); }; + + auto eqIdx = conti0EqIdx + phaseIdx; + flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindTerm); + + // Add advective phase energy fluxes. For isothermal model the contribution is zero. + if constexpr (enableEnergyBalance) + EnergyOperators::heatConvectionFlux(flux, fluxVars, phaseIdx); + } + + // Add diffusive energy fluxes. For isothermal model the contribution is zero. + if constexpr (enableEnergyBalance) + EnergyOperators::heatConductionFlux(flux, fluxVars); + + return flux; + } +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 1d7da0cd2512965ef27b5b7059236acd809cec4e Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:11:54 +0200 Subject: [PATCH 05/32] [assembly] add fv local operator --- dumux/assembly/fv/CMakeLists.txt | 1 + dumux/assembly/fv/localoperator.hh | 237 +++++++++++++++++++++++++++++ 2 files changed, 238 insertions(+) create mode 100644 dumux/assembly/fv/localoperator.hh diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt index 72883ac6ae..fb70764d7e 100644 --- a/dumux/assembly/fv/CMakeLists.txt +++ b/dumux/assembly/fv/CMakeLists.txt @@ -1,3 +1,4 @@ install(FILES +localoperator.hh operators.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv) diff --git a/dumux/assembly/fv/localoperator.hh b/dumux/assembly/fv/localoperator.hh new file mode 100644 index 0000000000..49f5ff3411 --- /dev/null +++ b/dumux/assembly/fv/localoperator.hh @@ -0,0 +1,237 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief An element-wise local operator for finite-volume schemes. + */ +#ifndef DUMUX_FV_LOCAL_OPERATOR_HH +#define DUMUX_FV_LOCAL_OPERATOR_HH + +#include + +#include +#include + +#include +#include +#include + +namespace Dumux::Experimental { + +/*! + * \ingroup Assembly + * \brief The element-wise local operator for finite volume schemes. + * This allows for element-wise evaluation of individual terms + * of the equations to be solved. + * \tparam OP The model-specific operators + */ +template +class FVLocalOperator +{ + using LC = typename OP::LocalContext; + using ElementVariables = typename LC::ElementVariables; + using GridVariables = typename ElementVariables::GridVariables; + using PrimaryVariables = typename GridVariables::PrimaryVariables; + using NumEqVector = Dumux::NumEqVector; + + using FVElementGeometry = typename LC::ElementGridGeometry; + using GridGeometry = typename FVElementGeometry::GridGeometry; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using Extrusion = Extrusion_t; + + static constexpr int numEq = NumEqVectorTraits::numEq(); + static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box; + +public: + //! export the expected local context type + using LocalContext = LC; + + //! export the underlying operators + using Operators = OP; + + //! vector type storing the operator values on all dofs of an element + //! TODO: Use ReservedBlockVector + using ElementResidualVector = Dune::BlockVector; + + //! Constructor from a local context + explicit FVLocalOperator(const LocalContext& context) + : context_(context) + {} + + /*! + * \name Main interface + * \note Methods used by the assembler to compute derivatives and residual + */ + // \{ + + /*! + * \brief Compute the terms of the local residual that do not appear in + * time derivatives. These are the sources and the fluxes. + */ + ElementResidualVector evalFluxesAndSources() const + { + auto result = getEmptyResidual(); + const auto& problem = context_.elementVariables().gridVariables().gridVolVars().problem(); + const auto& fvGeometry = context_.elementGridGeometry(); + + // source term + for (const auto& scv : scvs(fvGeometry)) + result[scv.localDofIndex()] -= Operators::source(problem, context_, scv); + + // flux term + for (const auto& scvf : scvfs(fvGeometry)) + addFlux_(result, scvf); + + return result; + } + + /*! + * \brief Compute the storage term, i.e. the term appearing in the time derivative. + */ + ElementResidualVector evalStorage() const + { + const auto& problem = context_.elementVariables().gridVariables().gridVolVars().problem(); + const auto& fvGeometry = context_.elementGridGeometry(); + + // TODO: Until now, FVLocalResidual defined storage as the entire + // time derivative. Now it means the term above the time derivative. + // We should think about the correct naming here... + // TODO: Should storage() NOT multiply with volume?? That was different until + // now but flux() also returns the flux multiplied with area so this should + // be more consistent + auto result = getEmptyResidual(); + for (const auto& scv : scvs(fvGeometry)) + result[scv.localDofIndex()] += Operators::storage(problem, context_, scv); + + return result; + } + + ElementResidualVector getEmptyResidual() const + { + ElementResidualVector res(context_.elementGridGeometry().numScv()); + res = 0.0; + return res; + } + + // \} + + /*! + * \name Interfaces for analytic Jacobian computation + */ + // \{ + + //! \todo TODO: Add interfaces. Or, should this be here at all!? + + //\} + + // \} + +protected: + + //! compute and add the flux across the given face to the container (cc schemes) + template = 0> + void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const + { + const auto& evv = context_.elementVariables().elemVolVars(); + const auto& problem = evv.gridVolVars().problem(); + + // TODO: Modify problem interfaces to receive context + const auto& fvGeometry = context_.elementGridGeometry(); + const auto& element = fvGeometry.element(); + const auto& efvc = context_.elementVariables().elemFluxVarsCache(); + + const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); + const auto localDofIdx = insideScv.localDofIndex(); + + if (!scvf.boundary()) + r[localDofIdx] += Operators::flux(problem, context_, scvf); + else + { + const auto& bcTypes = problem.boundaryTypes(element, scvf); + + // Dirichlet boundaries + if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) + r[localDofIdx] += Operators::flux(problem, context_, scvf); + + // Neumann and Robin ("solution dependent Neumann") boundary conditions + else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet()) + { + auto neumannFluxes = problem.neumann(element, fvGeometry, evv, efvc, scvf); + + // multiply neumann fluxes with the area and the extrusion factor + neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor(); + r[localDofIdx] += neumannFluxes; + } + + else + DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions for cell-centered schemes. " << + "Use pure boundary conditions by converting Dirichlet BCs to Robin BCs"); + } + } + + //! compute and add the flux across the given face to the container (box scheme) + template = 0> + void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const + { + const auto& evv = context_.elementVariables().elemVolVars(); + const auto& problem = evv.gridVolVars().problem(); + + // TODO: Modify problem interfaces to receive context + const auto& fvGeometry = context_.elementGridGeometry(); + const auto& element = fvGeometry.element(); + const auto& efvc = context_.elementVariables().elemFluxVarsCache(); + + // inner faces + if (!scvf.boundary()) + { + const auto flux = Operators::flux(problem, context_, scvf); + r[fvGeometry.scv(scvf.insideScvIdx()).localDofIndex()] += flux; + + if (scvf.numOutsideScvs() > 0) + r[fvGeometry.scv(scvf.outsideScvIdx()).localDofIndex()] -= flux; + } + + // boundary faces + else + { + const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); + const auto& bcTypes = problem.boundaryTypes(element, scv); + + // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions. + if (bcTypes.hasNeumann()) + { + const auto neumannFluxes = problem.neumann(element, fvGeometry, evv, efvc, scvf); + const auto area = Extrusion::area(scvf)*evv[scv].extrusionFactor(); + + // only add fluxes to equations for which Neumann is set + for (int eqIdx = 0; eqIdx < NumEqVector::size(); ++eqIdx) + if (bcTypes.isNeumann(eqIdx)) + r[scv.localDofIndex()][eqIdx] += neumannFluxes[eqIdx]*area; + } + } + } + +private: + const LocalContext& context_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 4e74341991281d0f4775a4618a6def10b9237ce9 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:12:51 +0200 Subject: [PATCH 06/32] [assembly] add new box local assembler --- dumux/assembly/fv/CMakeLists.txt | 1 + dumux/assembly/fv/boxlocalassembler.hh | 479 +++++++++++++++++++++++++ 2 files changed, 480 insertions(+) create mode 100644 dumux/assembly/fv/boxlocalassembler.hh diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt index fb70764d7e..40d6910f65 100644 --- a/dumux/assembly/fv/CMakeLists.txt +++ b/dumux/assembly/fv/CMakeLists.txt @@ -1,4 +1,5 @@ install(FILES +boxlocalassembler.hh localoperator.hh operators.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv) diff --git a/dumux/assembly/fv/boxlocalassembler.hh b/dumux/assembly/fv/boxlocalassembler.hh new file mode 100644 index 0000000000..9e49482ca1 --- /dev/null +++ b/dumux/assembly/fv/boxlocalassembler.hh @@ -0,0 +1,479 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief An assembler for Jacobian and residual contribution + * per element for the box scheme. + */ +#ifndef DUMUX_BOX_LOCAL_ASSEMBLER_HH +#define DUMUX_BOX_LOCAL_ASSEMBLER_HH + +#include +#include + +#include + +#include +#include + +#include +#include +#include + + +namespace Dumux::Experimental { + +/*! + * \ingroup Assembly + * \brief An assembler for Jacobian and residual contribution + * per element for the box scheme. + * \tparam LO The element-local operator + */ +template +class BoxLocalAssembler +{ + using LocalContext = typename LO::LocalContext; + using ElementVariables = typename LocalContext::ElementVariables; + using ElementGridGeometry = typename LocalContext::ElementGridGeometry; + + using GG = typename ElementGridGeometry::GridGeometry; + using GV = typename ElementVariables::GridVariables; + + using FVElementGeometry = typename GG::LocalView; + using SubControlVolume = typename GG::SubControlVolume; + using Element = typename GG::GridView::template Codim<0>::Entity; + + using PrimaryVariables = typename GV::PrimaryVariables; + using Scalar = typename GV::Scalar; + + static constexpr int numEq = NumEqVectorTraits::numEq; + +public: + + //! the parameters of a stage in time integration + using StageParams = MultiStageParams; + + //! export the grid variables type on which to operate + using GridVariables = GV; + + //! export the underlying local operator + using LocalOperator = LO; + + //! export the vector storing the residuals of all dofs of the element + using ElementResidualVector = typename LO::ElementResidualVector; + + /*! + * \brief Constructor for stationary problems. + */ + explicit BoxLocalAssembler(const Element& element, + GridVariables& gridVariables, + DiffMethod dm = DiffMethod::numeric) + : diffMethod_(dm) + , gridVariables_(gridVariables) + , fvGeometry_(localView(gridVariables.gridGeometry())) + , elemVariables_(localView(gridVariables)) + , prevElemVariables_() + , elementIsGhost_(element.partitionType() == Dune::GhostEntity) + , stageParams_(nullptr) + { + if (diffMethod_ != DiffMethod::numeric) + DUNE_THROW(Dune::NotImplemented, "Provided differentiation method"); + + fvGeometry_.bind(element); + elemVariables_.bind(element, fvGeometry_); + } + + /*! + * \brief Constructor for instationary problems. + * \note Using this constructor, we assemble one stage within + * a time integration step using multi-stage methods. + */ + explicit BoxLocalAssembler(const Element& element, + const std::vector& prevGridVars, + GridVariables& gridVariables, + std::shared_ptr stageParams, + DiffMethod dm = DiffMethod::numeric) + : diffMethod_(dm) + , gridVariables_(gridVariables) + , fvGeometry_(localView(gridVariables.gridGeometry())) + , elemVariables_(localView(gridVariables)) + , prevElemVariables_() + , elementIsGhost_(element.partitionType() == Dune::GhostEntity) + , stageParams_(stageParams) + { + if (diffMethod_ != DiffMethod::numeric) + DUNE_THROW(Dune::NotImplemented, "Provided differentiation method"); + + if (prevGridVars.size() != stageParams->size() - 1) + DUNE_THROW(Dune::InvalidStateException, "Grid Variables for all stages needed"); + + fvGeometry_.bind(element); + for (const auto& gridVars : prevGridVars) + { + prevElemVariables_.emplace_back(localView(gridVars)); + prevElemVariables_.back().bind(element, fvGeometry_); + } + + elemVariables_.bind(element, fvGeometry_); + } + + /*! + * \brief Computes the derivatives with respect to the given element and adds + * them to the global matrix. The element residual is written into the + * right hand side. + */ + template + void assembleJacobianAndResidual(JacobianMatrix& jac, + ResidualVector& res, + const PartialReassembler* partialReassembler = nullptr) + { + const auto& element = fvGeometry_.element(); + const auto eIdxGlobal = fvGeometry_.gridGeometry().elementMapper().index(element); + + if (partialReassembler && partialReassembler->elementColor(eIdxGlobal) == EntityColor::green) + { + const auto residual = evalLocalResidual(); + for (const auto& scv : scvs(fvGeometry_)) + res[scv.dofIndex()] += residual[scv.localDofIndex()]; + } + else if (!elementIsGhost_) + { + const auto residual = assembleJacobianAndResidual_(jac, partialReassembler); + for (const auto& scv : scvs(fvGeometry_)) + res[scv.dofIndex()] += residual[scv.localDofIndex()]; + } + else + { + static constexpr int dim = GG::GridView::dimension; + + for (const auto& scv : scvs(fvGeometry_)) + { + const auto vIdxLocal = scv.indexInElement(); + const auto& v = fvGeometry_.element().template subEntity(vIdxLocal); + + // do not change the non-ghost vertices + if (v.partitionType() == Dune::InteriorEntity || + v.partitionType() == Dune::BorderEntity) + continue; + + const auto dofIdx = scv.dofIndex(); + res[dofIdx] = 0; + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) + jac[dofIdx][dofIdx][pvIdx][pvIdx] = 1.0; + } + } + + auto applyDirichlet = [&] (const auto& scvI, + const auto& dirichletValues, + const auto eqIdx, + const auto pvIdx) + { + const auto& priVars = elemVariables_.elemVolVars()[scvI].priVars(); + res[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx]; + + auto& row = jac[scvI.dofIndex()]; + for (auto col = row.begin(); col != row.end(); ++col) + row[col.index()][eqIdx] = 0.0; + + jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0; + + // if a periodic dof has Dirichlet values also apply the same Dirichlet values to the other dof + if (fvGeometry_.gridGeometry().dofOnPeriodicBoundary(scvI.dofIndex())) + { + const auto periodicDof = fvGeometry_.gridGeometry().periodicallyMappedDof(scvI.dofIndex()); + res[periodicDof][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx]; + const auto end = jac[periodicDof].end(); + for (auto it = jac[periodicDof].begin(); it != end; ++it) + (*it) = periodicDof != it.index() ? 0.0 : 1.0; + } + }; + + enforceDirichletConstraints(applyDirichlet); + } + + /*! + * \brief Computes the derivatives with respect to the dofs of the given + * element and adds them to the given jacobian matrix. + */ + template + void assembleJacobian(JacobianMatrix& jac) + { + assembleJacobianAndResidual_(jac); + + auto applyDirichlet = [&] (const auto& scvI, + const auto& dirichletValues, + const auto eqIdx, + const auto pvIdx) + { + auto& row = jac[scvI.dofIndex()]; + for (auto col = row.begin(); col != row.end(); ++col) + row[col.index()][eqIdx] = 0.0; + + jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0; + }; + + enforceDirichletConstraints(applyDirichlet); + } + + //! Assemble the residual only + template + void assembleResidual(ResidualVector& res) + { + const auto residual = evalLocalResidual(); + for (const auto& scv : scvs(fvGeometry_)) + res[scv.dofIndex()] += residual[scv.localDofIndex()]; + + auto applyDirichlet = [&] (const auto& scvI, + const auto& dirichletValues, + const auto eqIdx, + const auto pvIdx) + { + const auto& priVars = elemVariables_.elemVolVars()[scvI].priVars(); + res[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx]; + }; + + enforceDirichletConstraints(applyDirichlet); + } + + //! Enforce Dirichlet constraints + template + void enforceDirichletConstraints(const ApplyFunction& applyDirichlet) + { + // enforce Dirichlet boundary conditions + evalDirichletBoundaries(applyDirichlet); + // TODO: take care of internal Dirichlet constraints (if enabled) + // enforceInternalDirichletConstraints(applyDirichlet); + } + + //! Evaluates Dirichlet boundaries + template< typename ApplyDirichletFunctionType > + void evalDirichletBoundaries(ApplyDirichletFunctionType applyDirichlet) + { + const auto& problem = gridVariables_.gridVolVars().problem(); + + for (const auto& scvI : scvs(fvGeometry_)) + { + const auto bcTypes = problem.boundaryTypes(fvGeometry_.element(), scvI); + if (bcTypes.hasDirichlet()) + { + const auto dirichletValues = problem.dirichlet(fvGeometry_.element(), scvI); + + // set the Dirichlet conditions in residual and jacobian + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (bcTypes.isDirichlet(eqIdx)) + { + const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx); + assert(0 <= pvIdx && pvIdx < numEq); + applyDirichlet(scvI, dirichletValues, eqIdx, pvIdx); + } + } + } + } + } + + //! Evaluate the complete local residual for the current element. + ElementResidualVector evalLocalResidual() const + { + if (isStationary()) + { + const auto context = makeLocalContext(fvGeometry_, elemVariables_); + LocalOperator localOperator(context); + return elementIsGhost_ ? localOperator.getEmptyResidual() + : localOperator.evalFluxesAndSources(); + } + else + { + ElementResidualVector residual(fvGeometry_.numScv()); + residual = 0.0; + + if (!elementIsGhost_) + { + // add the terms associated with previous stages + for (std::size_t k = 0; k < stageParams_->size()-1; ++k) + addStageTerms_(residual, k, + makeLocalContext(fvGeometry_, prevElemVariables_[k])); + + // add the terms of the current stage + addStageTerms_(residual, stageParams_->size()-1, + makeLocalContext(fvGeometry_, elemVariables_)); + } + + return residual; + } + } + + //! Return true if a stationary problem is assembled + bool isStationary() const + { return !stageParams_; } + +protected: + + //! add the terms of a stage to the current element residual + void addStageTerms_(ElementResidualVector& r, + std::size_t stageIdx, + const LocalContext& context) const + { + LocalOperator localOperator(context); + if (!stageParams_->skipTemporal(stageIdx)) + r.axpy(stageParams_->temporalWeight(stageIdx), + localOperator.evalStorage()); + if (!stageParams_->skipSpatial(stageIdx)) + r.axpy(stageParams_->spatialWeight(stageIdx), + localOperator.evalFluxesAndSources()); + } + + /*! + * \brief Computes the derivatives with respect to the dofs of the given + * element and adds them to the global matrix. + * \return The element residual at the current solution. + */ + template + ElementResidualVector assembleJacobianAndResidual_(JacobianMatrix& A, + const PartialReassembler* partialReassembler = nullptr) + { + // TODO: DO we need this to be constexpr? + if (diffMethod_ == DiffMethod::numeric) + return assembleJacobianAndResidualNumeric_(A, partialReassembler); + else + DUNE_THROW(Dune::NotImplemented, "Analytic assembler for box"); + } + + /*! + * \brief Computes the derivatives by means of numeric differentiation + * and adds them to the global matrix. + * \return The element residual at the current solution. + * \note This specialization is for the box scheme with numeric differentiation + */ + template + ElementResidualVector assembleJacobianAndResidualNumeric_(JacobianMatrix& A, + const PartialReassembler* partialReassembler = nullptr) + { + // alias for the variables of the current stage + auto& curVariables = elemVariables_; + auto& curElemVolVars = curVariables.elemVolVars(); + const auto& x = curVariables.gridVariables().dofs(); + const auto& element = fvGeometry_.element(); + const auto& problem = curElemVolVars.gridVolVars().problem(); + + const auto origResiduals = evalLocalResidual(); + const auto origElemSol = elementSolution(element, x, fvGeometry_.gridGeometry()); + auto elemSol = origElemSol; + + ////////////////////////////////////////////////////////////////////////////////////////////// + // Calculate derivatives of the residual of all dofs in element with respect to themselves. // + ////////////////////////////////////////////////////////////////////////////////////////////// + + ElementResidualVector partialDerivs(fvGeometry_.numScv()); + for (const auto& scvI : scvs(fvGeometry_)) + { + // dof index and corresponding actual pri vars + const auto globalI = scvI.dofIndex(); + const auto localI = scvI.localDofIndex(); + + const auto origCurVolVars = curElemVolVars[scvI]; + auto& curVolVars = getVolVarAcess_(curElemVolVars, gridVariables_.gridVolVars(), scvI); + + // calculate derivatives w.r.t to the privars at the dof at hand + for (int pvIdx = 0; pvIdx < numEq; pvIdx++) + { + partialDerivs = 0.0; + auto evalResiduals = [&](Scalar priVar) + { + // update the volume variables and compute element residual + elemSol[scvI.localDofIndex()][pvIdx] = priVar; + curVolVars.update(elemSol, problem, element, scvI); + return evalLocalResidual(); + }; + + static const NumericEpsilon numEps{problem.paramGroup()}; + static const int numDiffMethod = getParamFromGroup(problem.paramGroup(), + "Assembly.NumericDifferenceMethod"); + + // epsilon used for numeric differentiation + const auto eps = numEps(elemSol[localI][pvIdx], pvIdx); + + // derive the residuals numerically + NumericDifferentiation::partialDerivative(evalResiduals, elemSol[localI][pvIdx], partialDerivs, + origResiduals, eps, numDiffMethod); + + // TODO: Distinguish between implicit/explicit here. For explicit schemes, + // no entries between different scvs of an element are reserved. Thus, + // we currently get an error when using explicit schemes. + // TODO: Doesn't this mean we only have to evaluate the residual for a single + // scv instead of calling evalLocalResidual()? That computes the residuals + // and derivs for all other scvs of the element, too, which are never used. + // Note: this is the same in the current implementation of master. + // Should we try to optimize this for explicit schemes? Or adjust the Jacobian pattern? + // update the global stiffness matrix with the current partial derivatives + for (const auto& scvJ : scvs(fvGeometry_)) + { + const auto globalJ = scvJ.dofIndex(); + const auto localJ = scvJ.localDofIndex(); + + // don't add derivatives for green entities + if (!partialReassembler || partialReassembler->dofColor(globalJ) != EntityColor::green) + { + for (int eqIdx = 0; eqIdx < numEq; eqIdx++) + { + // A[i][col][eqIdx][pvIdx] is the rate of change of the + // the residual of equation 'eqIdx' at dof 'i' + // depending on the primary variable 'pvIdx' at dof 'col' + A[globalJ][globalI][eqIdx][pvIdx] += partialDerivs[localJ][eqIdx]; + } + } + } + + // restore the original element solution & volume variables + elemSol[localI][pvIdx] = origElemSol[localI][pvIdx]; + curVolVars = origCurVolVars; + + // TODO additional dof dependencies + } + } + + return origResiduals; + } + + //! Returns the volume variables of the local view (case of no caching) + template = 0> + auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv) + { return elemVolVars[scv]; } + + //! Returns the volume variables of the grid vol vars (case of global caching) + template = 0> + auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv) + { return gridVolVars.volVars(scv); } + + DiffMethod diffMethod_; //!< the type of differentiation method + GridVariables& gridVariables_; //!< reference to the grid variables + FVElementGeometry fvGeometry_; //!< element-local finite volume geometry + ElementVariables elemVariables_; //!< element variables of the current stage + std::vector prevElemVariables_; //!< element variables of prior stages + + bool elementIsGhost_; + std::shared_ptr stageParams_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 8277b6d25442661d0a16c265e3f3e75c6bdae56a Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 22 Oct 2020 17:14:08 +0200 Subject: [PATCH 07/32] [assembly] add scheme-agnostic assembler --- dumux/assembly/CMakeLists.txt | 1 + dumux/assembly/assembler.hh | 463 ++++++++++++++++++++++++++++++++++ 2 files changed, 464 insertions(+) create mode 100644 dumux/assembly/assembler.hh diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt index 2543a6631f..17a49d7b47 100644 --- a/dumux/assembly/CMakeLists.txt +++ b/dumux/assembly/CMakeLists.txt @@ -1,4 +1,5 @@ install(FILES +assembler.hh boxlocalassembler.hh boxlocalresidual.hh cclocalassembler.hh diff --git a/dumux/assembly/assembler.hh b/dumux/assembly/assembler.hh new file mode 100644 index 0000000000..fa79c7d9fa --- /dev/null +++ b/dumux/assembly/assembler.hh @@ -0,0 +1,463 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief Assembler class for residuals and jacobian matrices + * for grid-based numerical schemes. + */ +#ifndef DUMUX_ASSEMBLER_HH +#define DUMUX_ASSEMBLER_HH + +#include + +#include +#include +#include +#include + +#include +#include + +#include "jacobianpattern.hh" + +// TODO: Remove when linear algebra traits are available +#include +#include +#include +#include + +namespace Dumux::Experimental { + +//! Default types used for the linear system +template +struct DefaultLinearSystemTraits; + +//! Default linear system types for Dune::BlockVector +template +struct DefaultLinearSystemTraits> +{ +private: + static constexpr int numEq = NumEqVectorTraits::numEq; + using Scalar = typename PrimaryVariables::value_type; + using MatrixBlockType = Dune::FieldMatrix; + +public: + using ResidualVector = Dune::BlockVector; + using JacobianMatrix = Dune::BCRSMatrix; +}; + +/*! + * \ingroup Assembly + * \brief A linear system assembler (residual and Jacobian) for grid-based numerical schemes + * \tparam LA The local assembler (element-local assembly) + * \tparam LST The linear system traits (types used for jacobians and residuals) + */ +template> +class Assembler +{ + using GV = typename LA::GridVariables; + using GG = typename GV::GridGeometry; + + using GridView = typename GG::GridView; + using Element = typename GridView::template Codim<0>::Entity; + +public: + //! export the types used for the linear system + using JacobianMatrix = typename LST::JacobianMatrix; + using ResidualVector = typename LST::ResidualVector; + using ResidualType = ResidualVector; + + //! export the types underlying discrete solutions + using Scalar = typename GV::Scalar; + using SolutionVector = typename GV::SolutionVector; + + //! export the types involved in element-local assembly + using LocalAssembler = LA; + using LocalOperator = typename LA::LocalOperator; + + //! the variables representing an evaluation point + using GridVariables = GV; + using Variables = GridVariables; + + //! export the underlying grid geometry type + using GridGeometry = GG; + + //! export the type for parameters of a stage in time integration + using StageParams = MultiStageParams; + + /*! + * \brief The Constructor from a grid geometry. + * \param gridGeometry A grid geometry instance + * \param diffMethod The method used for computing partial derivatives + * \note This assembler class is, after construction, defined for a specific equation + * (defined by the LocalOperator inside the LocalAssembler) and a specific grid + * geometry - which defines the connectivity of the degrees of freedom of the + * underlying discretization scheme on a particular grid. The evaluation point, + * consisting of a particular solution/variables/parameters may vary, and therefore, + * an instance of the grid variables is passed to the assembly functions. + * \note This constructor (without time integration method) assumes hhat a stationary problem + * is assembled, and thus, an implicit jacobian pattern is used. + */ + Assembler(std::shared_ptr gridGeometry, + DiffMethod diffMethod = DiffMethod::numeric) + : gridGeometry_(gridGeometry) + , isImplicit_(true) + , diffMethod_(diffMethod) + {} + + /*! + * \brief The Constructor from a grid geometry. + * \param gridGeometry A grid geometry instance + * \param method the time integration method that will be used. + */ + template + Assembler(std::shared_ptr gridGeometry, + const TimeIntegrationMethod& method, + DiffMethod diffMethod = DiffMethod::numeric) + : gridGeometry_(gridGeometry) + , isImplicit_(method.implicit()) + , diffMethod_(diffMethod) + {} + + /*! + * \brief Assembles the Jacobian matrix and the residual around the given evaluation point + * which is determined by the grid variables, containing all quantities required + * to evaluate the equations to be assembled. + * \param gridVariables The variables corresponding to the given solution state + * \note We assume the grid geometry on which the grid variables are defined + * to be the same as the one used to instantiate this class. + * \todo TODO: The grid variables have to be non-const in order to compute + * the derivatives in the case of global caching. Could this be + * circumvented somehow? + */ + template + void assembleJacobianAndResidual(GridVariables& gridVariables, + const PartialReassembler* partialReassembler = nullptr) + { + assert(gridVariables.gridGeometry().numDofs() == gridGeometry().numDofs()); + + resetJacobian_(partialReassembler); + resetResidual_(); + + assemble_([&](const Element& element) + { + auto localAssembler = this->makeLocalAssembler_(element, gridVariables); + localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, partialReassembler); + }); + + // TODO: Put these into discretization-specific helpers? + enforceDirichletConstraints_(gridVariables, *jacobian_, *residual_); + enforceInternalConstraints_(gridVariables, *jacobian_, *residual_); + enforcePeriodicConstraints_(gridVariables, *jacobian_, *residual_); + } + + /*! + * \brief Assembles the Jacobian matrix of the discrete system of equations + * around a given state represented by the grid variables object. + */ + void assembleJacobian(GridVariables& gridVariables) + { + assert(gridVariables.gridGeometry().numDofs() == gridGeometry().numDofs()); + + resetJacobian_(); + + assemble_([&](const Element& element) + { + auto localAssembler = this->makeLocalAssembler_(element, gridVariables); + localAssembler.assembleJacobianAndResidual(*jacobian_); + }); + + // TODO: Put these into discretization-specific helpers? + enforceDirichletConstraints_(gridVariables, *jacobian_, *residual_); + enforceInternalConstraints_(gridVariables, *jacobian_, *residual_); + enforcePeriodicConstraints_(gridVariables, *jacobian_, *residual_); + } + + /*! + * \brief Assembles the residual for a given state represented by the provided + * grid variables object, using the internal residual vector to store the result. + */ + void assembleResidual(GridVariables& gridVariables) + { + resetResidual_(); + assembleResidual(*residual_, gridVariables); + } + + /*! + * \brief Assembles the residual for a given state represented by the provided + * grid variables object, using the provided residual vector to store the result. + */ + void assembleResidual(ResidualVector& r, GridVariables& gridVariables) const + { + assemble_([&](const Element& element) + { this->makeLocalAssembler_(element, gridVariables).assembleResidual(r); }); + } + + //! Will become obsolete once the new linear solvers are available + Scalar residualNorm(GridVariables& gridVars) const + { + ResidualType residual(numDofs()); + assembleResidual(residual, gridVars); + + // for box communicate the residual with the neighboring processes + if (GridGeometry::discMethod == DiscretizationMethod::box && gridView().comm().size() > 1) + { + using VertexMapper = typename GridGeometry::VertexMapper; + VectorCommDataHandleSum + sumResidualHandle(gridGeometry_->vertexMapper(), residual); + gridView().communicate(sumResidualHandle, + Dune::InteriorBorder_InteriorBorder_Interface, + Dune::ForwardCommunication); + } + + // calculate the square norm of the residual + Scalar result2 = residual.two_norm2(); + if (gridView().comm().size() > 1) + result2 = gridView().comm().sum(result2); + + using std::sqrt; + return sqrt(result2); + } + + + /*! + * \brief Tells the assembler which jacobian and residual to use. + * This also resizes the containers to the required sizes and sets the + * sparsity pattern of the jacobian matrix. + */ + void setLinearSystem(std::shared_ptr A, + std::shared_ptr r) + { + jacobian_ = A; + residual_ = r; + + // check and/or set the BCRS matrix's build mode + if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown) + jacobian_->setBuildMode(JacobianMatrix::random); + else if (jacobian_->buildMode() != JacobianMatrix::BuildMode::random) + DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment"); + + setJacobianPattern(); + setResidualSize(); + } + + /*! + * \brief The version without arguments uses the default constructor to create + * the jacobian and residual objects in this assembler if you don't need them outside this class + */ + void setLinearSystem() + { + jacobian_ = std::make_shared(); + jacobian_->setBuildMode(JacobianMatrix::random); + residual_ = std::make_shared(); + + setJacobianPattern(); + setResidualSize(); + } + + /*! + * \brief Resizes the jacobian and sets the jacobian' sparsity pattern. + */ + void setJacobianPattern() + { + // resize the jacobian and the residual + const auto numDofs = this->numDofs(); + jacobian_->setSize(numDofs, numDofs); + + // create occupation pattern of the jacobian + // TODO: Does the bool need to be at compile time in getJacobianPattern? + const auto occupationPattern = isImplicit_ ? getJacobianPattern(gridGeometry()) + : getJacobianPattern(gridGeometry()); + + // export pattern to jacobian + occupationPattern.exportIdx(*jacobian_); + } + + /*! + * \brief Prepare for a new stage within a time integration step. + * This caches the given grid variables, which are then used as a + * representation of the previous stage. Moreover, the given grid + * variables are then updated to the time level of the upcoming stage. + * \param gridVars the grid variables representing the previous stage + * \param params the parameters with the weights to be used in the upcoming stage + */ + void prepareStage(GridVariables& gridVars, + std::shared_ptr params) + { + stageParams_ = params; + const auto curStage = params->size() - 1; + + // we keep track of previous stages, they are needed for residual assembly + prevStageVariables_.push_back(gridVars); + + // Now we update the time level of the given grid variables + const auto t = params->timeAtStage(curStage); + const auto prevT = params->timeAtStage(0); + const auto dtFraction = params->timeStepFraction(curStage); + gridVars.updateTime(TimeLevel{t, prevT, dtFraction}); + } + + /*! + * \brief Remove traces from stages within a time integration step. + */ + void clearStages() + { + prevStageVariables_.clear(); + stageParams_ = nullptr; + } + + //! Resizes the residual + void setResidualSize() + { residual_->resize(numDofs()); } + + //! Returns the number of degrees of freedom + std::size_t numDofs() const + { return gridGeometry().numDofs(); } + + //! The global finite volume geometry + const GridGeometry& gridGeometry() const + { return *gridGeometry_; } + + //! The grid view + const GridView& gridView() const + { return gridGeometry().gridView(); } + + //! The jacobian matrix + JacobianMatrix& jacobian() + { return *jacobian_; } + + //! The residual vector (rhs) + ResidualVector& residual() + { return *residual_; } + +protected: + // make a local assembler instance + LocalAssembler makeLocalAssembler_(const Element& element, GridVariables& gridVars) const + { + // instationary assembly + if (!stageParams_) + return LocalAssembler{element, gridVars, diffMethod_}; + + // stationary assembly + return LocalAssembler{element, prevStageVariables_, gridVars, stageParams_, diffMethod_}; + } + + // reset the residual vector to 0.0 + void resetResidual_() + { + if (!residual_) + { + residual_ = std::make_shared(); + setResidualSize(); + } + + (*residual_) = 0.0; + } + + // reset the Jacobian matrix to 0.0 + template + void resetJacobian_(const PartialReassembler *partialReassembler = nullptr) + { + if (!jacobian_) + { + jacobian_ = std::make_shared(); + jacobian_->setBuildMode(JacobianMatrix::random); + setJacobianPattern(); + } + + if (partialReassembler) + partialReassembler->resetJacobian(*this); + else + *jacobian_ = 0.0; + } + + /*! + * \brief A method assembling something per element + * \note Handles exceptions for parallel runs + * \throws NumericalProblem on all processes if something throwed during assembly + */ + template + void assemble_(AssembleElementFunc&& assembleElement) const + { + // a state that will be checked on all processes + bool succeeded = false; + + // try assembling using the local assembly function + try + { + // let the local assembler add the element contributions + for (const auto& element : elements(gridView())) + assembleElement(element); + + // if we get here, everything worked well on this process + succeeded = true; + } + // throw exception if a problem ocurred + catch (NumericalProblem &e) + { + std::cout << "rank " << gridView().comm().rank() + << " caught an exception while assembling:" << e.what() + << "\n"; + succeeded = false; + } + + // make sure everything worked well on all processes + if (gridView().comm().size() > 1) + succeeded = gridView().comm().min(succeeded); + + // if not succeeded rethrow the error on all processes + if (!succeeded) + DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system"); + } + + void enforceDirichletConstraints_(const GridVariables& gridVars, JacobianMatrix& jac, SolutionVector& res) + { /*TODO: Implement / where to put this? Currently, local assemblers do this*/ } + + void enforceInternalConstraints_(const GridVariables& gridVars, JacobianMatrix& jac, SolutionVector& res) + { /*TODO: Implement / where to put this? Currently, local assemblers do this*/ } + + void enforcePeriodicConstraints_(const GridVariables& gridVars, JacobianMatrix& jac, SolutionVector& res) + { /*TODO: Implement / where to put this? Currently, local assemblers do this*/ } + +private: + //! the grid geometry on which it is assembled + std::shared_ptr gridGeometry_; + + //! states if an implicit of explicit scheme is used (affects jacobian pattern) + bool isImplicit_; + + //! the method to compute the derivatives + DiffMethod diffMethod_; + + //! shared pointers to the jacobian matrix and residual + std::shared_ptr jacobian_; + std::shared_ptr residual_; + + //! parameters containing information on the current stage of time integration + std::shared_ptr stageParams_ = nullptr; + + //! keep track of the states of previous stages within a time integration step + std::vector prevStageVariables_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 91c7a7b1b14fe5438d23052a1a545b10f8efa470 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 15:09:39 +0100 Subject: [PATCH 08/32] [experimental][gridvars] add helper boolean --- dumux/discretization/fvgridvariables.hh | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh index b4b078629f..d0ea810c50 100644 --- a/dumux/discretization/fvgridvariables.hh +++ b/dumux/discretization/fvgridvariables.hh @@ -173,7 +173,10 @@ private: // Experimental implementation of new grid variables layout // ////////////////////////////////////////////////////////////// +#include +#include #include + #include #include @@ -374,6 +377,23 @@ private: GridFluxVariablesCache gridFluxVarsCache_; //!< the flux variables cache }; + namespace Detail { + struct hasGridVolVars + { + template + auto operator()(const GV& gv) -> decltype(gv.gridVolVars()) {} + }; + } // end namespace Detail + + // helper bool to check if a grid variables type fulfills the new experimental interface + // can be used in the transition period in places where both interfaces should be supported + // Note that this doesn't check if the provided type actually models grid variables. Instead, + // it is assumed that that is known, and it is only checked if the new or old behaviour can be + // expected from it. + template + inline constexpr bool areExperimentalGridVars = + decltype(isValid(Detail::hasGridVolVars())(std::declval()))::value; + } // end namespace Dumux::Experimental #endif -- GitLab From e31392399ad1c43e2c6143675a5c669b9c784067 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 15:10:26 +0100 Subject: [PATCH 09/32] [newton] use non-const assembly arguments The background is that in the case of global caching, we need to pass the grid variables as non-const in order to be able to deflect them during numeric differentiation. It would be nice if we could find a way around it. --- dumux/nonlinear/newtonsolver.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh index 121e68cac2..e1fd60ef37 100644 --- a/dumux/nonlinear/newtonsolver.hh +++ b/dumux/nonlinear/newtonsolver.hh @@ -481,7 +481,7 @@ public: * * \param vars The current iteration's variables */ - virtual void assembleLinearSystem(const Variables& vars) + virtual void assembleLinearSystem(Variables& vars) { assembleLinearSystem_(this->assembler(), vars); @@ -879,7 +879,7 @@ protected: this->assembler().updateGridVariables(Backend::dofs(vars)); } - void computeResidualReduction_(const Variables& vars) + void computeResidualReduction_(Variables& vars) { // we assume that the assembler works on solution vectors // if it doesn't export the variables type @@ -1067,7 +1067,7 @@ private: //! assembleLinearSystem_ for assemblers that support partial reassembly template - auto assembleLinearSystem_(const A& assembler, const Variables& vars) + auto assembleLinearSystem_(const A& assembler, Variables& vars) -> typename std::enable_if_t { this->assembler().assembleJacobianAndResidual(vars, partialReassembler_.get()); @@ -1075,7 +1075,7 @@ private: //! assembleLinearSystem_ for assemblers that don't support partial reassembly template - auto assembleLinearSystem_(const A& assembler, const Variables& vars) + auto assembleLinearSystem_(const A& assembler, Variables& vars) -> typename std::enable_if_t { this->assembler().assembleJacobianAndResidual(vars); -- GitLab From d8cab3ddc0fcd5201b6c44d70c8c2c04cd29201a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 15:11:58 +0100 Subject: [PATCH 10/32] [io][vtk] introduce compatibility layer with new/old grid vars --- dumux/io/vtkoutputmodule.hh | 18 +++++++++++++++--- dumux/porousmediumflow/velocity.hh | 12 +++++++++++- dumux/porousmediumflow/velocityoutput.hh | 18 +++++++++++++++++- 3 files changed, 43 insertions(+), 5 deletions(-) diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh index 0805552aae..faec1248d4 100644 --- a/dumux/io/vtkoutputmodule.hh +++ b/dumux/io/vtkoutputmodule.hh @@ -43,7 +43,9 @@ #include #include + #include +#include #include "vtkfunction.hh" #include "velocityoutput.hh" @@ -378,8 +380,18 @@ public: } protected: + // extract the grid volume variables from the grid variables (experimental interface) + template, + std::enable_if_t = 0> + decltype(auto) gridVolVars() const { return gridVariables_.gridVolVars(); } + + // extract the grid volume variables from the grid variables (experimental interface) + template, + std::enable_if_t = 0> + decltype(auto) gridVolVars() const { return gridVariables_.curGridVolVars(); } + // some return functions for differing implementations to use - const auto& problem() const { return gridVariables_.curGridVolVars().problem(); } + const auto& problem() const { return gridVolVars().problem(); } const GridVariables& gridVariables() const { return gridVariables_; } const GridGeometry& gridGeometry() const { return gridVariables_.gridGeometry(); } const SolutionVector& sol() const { return sol_; } @@ -447,7 +459,7 @@ private: const auto eIdxGlobal = gridGeometry().elementMapper().index(element); auto fvGeometry = localView(gridGeometry()); - auto elemVolVars = localView(gridVariables_.curGridVolVars()); + auto elemVolVars = localView(gridVolVars()); // If velocity output is enabled we need to bind to the whole stencil // otherwise element-local data is sufficient @@ -634,7 +646,7 @@ private: const auto numCorners = element.subEntities(dim); auto fvGeometry = localView(gridGeometry()); - auto elemVolVars = localView(gridVariables_.curGridVolVars()); + auto elemVolVars = localView(gridVolVars()); // resize element-local data containers for (std::size_t i = 0; i < volVarScalarDataInfo_.size(); ++i) diff --git a/dumux/porousmediumflow/velocity.hh b/dumux/porousmediumflow/velocity.hh index e34bf588da..088e626cbb 100644 --- a/dumux/porousmediumflow/velocity.hh +++ b/dumux/porousmediumflow/velocity.hh @@ -36,6 +36,7 @@ #include #include #include +#include #include namespace Dumux { @@ -108,7 +109,7 @@ public: * \param gridVariables The grid variables */ PorousMediumFlowVelocity(const GridVariables& gridVariables) - : problem_(gridVariables.curGridVolVars().problem()) + : problem_(getProblem_(gridVariables)) , gridGeometry_(gridVariables.gridGeometry()) , gridVariables_(gridVariables) { @@ -460,6 +461,15 @@ private: } private: + + // extract problem from grid vars (experimental interface) + template, int> = 0> + const Problem& getProblem_(const GV& gv) { return gv.gridVolVars().problem(); } + + // extract problem from grid vars (standards interface) + template, int> = 0> + const Problem& getProblem_(const GV& gv) { return gv.curGridVolVars().problem(); } + const Problem& problem_; const GridGeometry& gridGeometry_; const GridVariables& gridVariables_; diff --git a/dumux/porousmediumflow/velocityoutput.hh b/dumux/porousmediumflow/velocityoutput.hh index a3319b14e8..227f2c4298 100644 --- a/dumux/porousmediumflow/velocityoutput.hh +++ b/dumux/porousmediumflow/velocityoutput.hh @@ -29,9 +29,11 @@ #include #include + #include #include #include +#include #include namespace Dumux { @@ -67,6 +69,15 @@ class PorousMediumFlowVelocityOutput : public VelocityOutput using Problem = typename GridVolumeVariables::Problem; using VelocityBackend = PorousMediumFlowVelocity; + struct hasCurGridVolVars + { + template + auto operator()(const GV& gv) -> decltype(gv.curGridVolVars()) {} + }; + + static constexpr auto isOldGVInterface = + decltype(isValid(hasCurGridVolVars())(std::declval()))::value; + public: using VelocityVector = typename ParentType::VelocityVector; @@ -78,7 +89,12 @@ public: PorousMediumFlowVelocityOutput(const GridVariables& gridVariables) { // check, if velocity output can be used (works only for cubes so far) - enableOutput_ = getParamFromGroup(gridVariables.curGridVolVars().problem().paramGroup(), "Vtk.AddVelocity"); + // compatibility layer with new and old-style grid variables + if constexpr (Experimental::areExperimentalGridVars) + enableOutput_ = getParamFromGroup(gridVariables.gridVolVars().problem().paramGroup(), "Vtk.AddVelocity"); + else + enableOutput_ = getParamFromGroup(gridVariables.curGridVolVars().problem().paramGroup(), "Vtk.AddVelocity"); + if (enableOutput_) velocityBackend = std::make_unique(gridVariables); } -- GitLab From 91aa77ffc469453b5e7fa1144620866b6021f0de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 15:13:17 +0100 Subject: [PATCH 11/32] [test][1p][compressible][box] use new time integration --- .../compressible/instationary/CMakeLists.txt | 10 +-- .../1p/compressible/instationary/assembler.hh | 54 --------------- .../instationary/gridvariables.hh | 68 ------------------- .../instationary/main_experimental.cc | 44 ++++++++---- .../compressible/instationary/properties.hh | 5 +- 5 files changed, 35 insertions(+), 146 deletions(-) delete mode 100644 test/porousmediumflow/1p/compressible/instationary/assembler.hh delete mode 100644 test/porousmediumflow/1p/compressible/instationary/gridvariables.hh diff --git a/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt b/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt index f824ba4a8f..05d17f3dcf 100644 --- a/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt +++ b/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt @@ -31,13 +31,13 @@ dumux_add_test(NAME test_1p_compressible_instationary_box ${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_box-00010.vtu --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_box params.input -Problem.Name test_1p_compressible_instationary_box") -dumux_add_test(NAME test_1p_compressible_instationary_tpfa_experimental +dumux_add_test(NAME test_1p_compressible_instationary_box_experimental LABELS porousmediumflow 1p experimental SOURCES main_experimental.cc - COMPILE_DEFINITIONS TYPETAG=OnePCompressibleTpfa + COMPILE_DEFINITIONS TYPETAG=OnePCompressibleBox COMPILE_DEFINITIONS EXPERIMENTAL=1 COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py CMD_ARGS --script fuzzy - --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_cc-reference.vtu - ${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_tpfa_experimental-00010.vtu - --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_tpfa_experimental params.input -Problem.Name test_1p_compressible_instationary_tpfa_experimental") + --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_box-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_box_experimental-00010.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_box_experimental params.input -Problem.Name test_1p_compressible_instationary_box_experimental") diff --git a/test/porousmediumflow/1p/compressible/instationary/assembler.hh b/test/porousmediumflow/1p/compressible/instationary/assembler.hh deleted file mode 100644 index 115b25570c..0000000000 --- a/test/porousmediumflow/1p/compressible/instationary/assembler.hh +++ /dev/null @@ -1,54 +0,0 @@ -// -*- 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 3 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 . * - *****************************************************************************/ -/*! - * \file - * \ingroup OnePTests - * \brief Dummy assembler that fulfills the new experimental assembly interface. - */ -#ifndef DUMUX_COMPRESSIBLE_ONEP_TEST_ASSEMBLER_HH -#define DUMUX_COMPRESSIBLE_ONEP_TEST_ASSEMBLER_HH - -namespace Dumux::OnePCompressibleTest { - -// Custom assembler to test assembly with grid variables, -// fulfilling the foreseen required interface -template -class GridVarsAssembler : public Assembler -{ -public: - using Assembler::Assembler; - using typename Assembler::GridVariables; - using typename Assembler::ResidualType; - - using Variables = GridVariables; - - void assembleJacobianAndResidual(const GridVariables& gridVars) - { Assembler::assembleJacobianAndResidual(gridVars.dofs()); } - - void assembleResidual(const GridVariables& gridVars) - { Assembler::assembleResidual(gridVars.dofs()); } - - //! Remove residualNorm once this is fixed in the solvers !2113 - auto residualNorm(const GridVariables& gridVars) - { return Assembler::residualNorm(gridVars.dofs()); } -}; - -} // end namespace Dumux::OnePCompressibleTest - -#endif diff --git a/test/porousmediumflow/1p/compressible/instationary/gridvariables.hh b/test/porousmediumflow/1p/compressible/instationary/gridvariables.hh deleted file mode 100644 index c9234e9c05..0000000000 --- a/test/porousmediumflow/1p/compressible/instationary/gridvariables.hh +++ /dev/null @@ -1,68 +0,0 @@ -// -*- 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 3 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 . * - *****************************************************************************/ -/*! - * \file - * \ingroup OnePTests - * \brief Wrapper around the current FVGridVariables to fulfill the layout - * of the new grid variables to test grid variables-based assembly. - */ -#ifndef DUMUX_COMPRESSIBLE_ONEP_TEST_GRID_VARIABLES_HH -#define DUMUX_COMPRESSIBLE_ONEP_TEST_GRID_VARIABLES_HH - -#include - -namespace Dumux::OnePCompressibleTest { - -template -class TestGridVariables -: public BaseGridVariables -, public Dumux::Experimental::GridVariables -{ - using ExperimentalBase = Dumux::Experimental::GridVariables; - -public: - // export some types to avoid ambiguity - using GridGeometry = GG; - using Scalar = typename BaseGridVariables::Scalar; - - template - TestGridVariables(std::shared_ptr problem, - std::shared_ptr gridGeometry, - const SolutionVector& x) - : BaseGridVariables(problem, gridGeometry) - , ExperimentalBase(gridGeometry, x) - { - BaseGridVariables::init(x); - } - - // update to a new solution - void update(const SolutionVector& x) - { - BaseGridVariables::update(x); - ExperimentalBase::update(x); - } - - // overload some functions to avoid ambiguity - decltype(auto) gridGeometry() const - { return ExperimentalBase::gridGeometry(); } -}; - -} // end namespace Dumux::OnePCompressibleTest - -#endif diff --git a/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc b/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc index 06f71a35b3..4b2281feb5 100644 --- a/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc +++ b/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc @@ -19,7 +19,7 @@ /*! * \file * \ingroup OnePTests - * \brief Test for the one-phase CC model + * \brief Test for the one-phase model */ #include #include @@ -31,6 +31,7 @@ #include #include +#include #include #include @@ -39,13 +40,17 @@ #include #include -#include +#include +#include + +#include +#include +#include #include #include #include "properties.hh" -#include "assembler.hh" int main(int argc, char** argv) { @@ -94,7 +99,6 @@ int main(int argc, char** argv) using SolutionVector = GetPropType; SolutionVector x(gridGeometry->numDofs()); problem->applyInitialSolution(x); - auto xOld = x; // the grid variables using GridVariables = GetPropType; @@ -118,10 +122,19 @@ int main(int argc, char** argv) auto timeLoop = std::make_shared>(0.0, dt, tEnd); timeLoop->setMaxTimeStepSize(maxDt); - // the assembler with time loop for instationary problem - using BaseAssembler = FVAssembler; - using Assembler = Dumux::OnePCompressibleTest::GridVarsAssembler; - auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); + // use implicit Euler for time integration + auto timeMethod = std::make_shared>(); + + // the assembler (we use the immiscible operators to define the system of equations) + using ModelTraits = GetPropType; + using FluxVariables = GetPropType; + using LocalContext = Experimental::LocalContext; + using Operators = Experimental::FVImmiscibleOperators; + + using LocalOperator = Experimental::FVLocalOperator; + using LocalAssembler = Experimental::BoxLocalAssembler; + using Assembler = Experimental::Assembler; + auto assembler = std::make_shared(gridGeometry, DiffMethod::numeric); // the linear solver using LinearSolver = ILU0BiCGSTABBackend; @@ -129,7 +142,11 @@ int main(int argc, char** argv) // the non-linear solver using NewtonSolver = Dumux::NewtonSolver; - NewtonSolver nonLinearSolver(assembler, linearSolver); + auto nonLinearSolver = std::make_shared(assembler, linearSolver); + + // the time stepper for time integration + using TimeStepper = MultiStageTimeStepper; + TimeStepper timeStepper(nonLinearSolver, timeMethod); // set some check points for the time loop timeLoop->setPeriodicCheckPoint(tEnd/10.0); @@ -137,11 +154,8 @@ int main(int argc, char** argv) // time loop timeLoop->start(); do { - // linearize & solve - nonLinearSolver.solve(*gridVariables, *timeLoop); - - // make the new solution the old solution - xOld = gridVariables->dofs(); + // do time integraiton + timeStepper.step(*gridVariables, timeLoop->time(), timeLoop->timeStepSize()); // advance to the time loop to the next step timeLoop->advanceTimeStep(); @@ -154,7 +168,7 @@ int main(int argc, char** argv) timeLoop->reportTimeStep(); // set new dt as suggested by the newton solver - timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); + timeLoop->setTimeStepSize(nonLinearSolver->suggestTimeStepSize(timeLoop->timeStepSize())); } while (!timeLoop->finished()); diff --git a/test/porousmediumflow/1p/compressible/instationary/properties.hh b/test/porousmediumflow/1p/compressible/instationary/properties.hh index 98a21fb4c9..f98db1bbd7 100644 --- a/test/porousmediumflow/1p/compressible/instationary/properties.hh +++ b/test/porousmediumflow/1p/compressible/instationary/properties.hh @@ -38,7 +38,6 @@ #include "problem.hh" #include "spatialparams.hh" -#include "gridvariables.hh" #ifndef EXPERIMENTAL #define EXPERIMENTAL 0 @@ -88,14 +87,12 @@ template struct GridVariables { private: - using GG = GetPropType; using GVV = GetPropType; using GFC = GetPropType; - using Base = Dumux::FVGridVariables; using X = GetPropType; public: - using type = Dumux::OnePCompressibleTest::TestGridVariables; + using type = Dumux::Experimental::FVGridVariables; }; #endif -- GitLab From b3ff6f71a1981fa9eafe07b6ef9c45fb118f4e4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 18:04:43 +0100 Subject: [PATCH 12/32] [fv][localop] add convenience function for scvf flux --- dumux/assembly/fv/localoperator.hh | 78 ++++++++++++++++++++++++------ 1 file changed, 62 insertions(+), 16 deletions(-) diff --git a/dumux/assembly/fv/localoperator.hh b/dumux/assembly/fv/localoperator.hh index 49f5ff3411..55483f5fe3 100644 --- a/dumux/assembly/fv/localoperator.hh +++ b/dumux/assembly/fv/localoperator.hh @@ -53,6 +53,7 @@ class FVLocalOperator using FVElementGeometry = typename LC::ElementGridGeometry; using GridGeometry = typename FVElementGeometry::GridGeometry; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; using Extrusion = Extrusion_t; @@ -143,20 +144,59 @@ public: // \} + /*! + * \brief Compute the flux across a single face embedded in the given element. + * \note This function behaves different than the flux function in the + * operators, as it checks if the face is on the boundary. If so, + * it returns the flux resulting from the user-specified conditions. + * \note This assumes that the given element and face are contained within + * the context that was used to instantiate this class. + */ + NumEqVector computeFlux(const Element& element, + const SubControlVolumeFace& scvf) const + { + const auto& evv = context_.elementVariables().elemVolVars(); + const auto& problem = evv.gridVolVars().problem(); + + // interior faces + if (!scvf.boundary()) + return Operators::flux(problem, context_, scvf); + + const auto& insideScv = context_.elementGridGeometry().scv(scvf.insideScvIdx()); + + // for cell-centered schemes, evaluate fluxes also across Dirichlet boundaries + if constexpr (!isBox) + { + const auto& bcTypes = problem.boundaryTypes(element, scvf); + if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) + return Operators::flux(problem, context_, scvf); + else if (bcTypes.hasNeumann() && bcTypes.hasDirichlet()) + DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions for cell-centered schemes. " << + "Use pure boundary conditions by converting Dirichlet BCs to Robin BCs."); + } + // for the box scheme, only pure Neumann boundaries are supported + else + { + if (!problem.boundaryTypes(element, insideScv).hasOnlyNeumann()) + DUNE_THROW(Dune::NotImplemented, "For the box scheme only Neumann boundary fluxes can be computed."); + } + + // compute and return the Neumann boundary fluxes + auto neumannFluxes = getUnscaledNeumannFluxes_(element, scvf); + neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor(); + return neumannFluxes; + } + protected: //! compute and add the flux across the given face to the container (cc schemes) template = 0> void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const { + const auto& fvGeometry = context_.elementGridGeometry(); const auto& evv = context_.elementVariables().elemVolVars(); const auto& problem = evv.gridVolVars().problem(); - // TODO: Modify problem interfaces to receive context - const auto& fvGeometry = context_.elementGridGeometry(); - const auto& element = fvGeometry.element(); - const auto& efvc = context_.elementVariables().elemFluxVarsCache(); - const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); const auto localDofIdx = insideScv.localDofIndex(); @@ -164,7 +204,7 @@ protected: r[localDofIdx] += Operators::flux(problem, context_, scvf); else { - const auto& bcTypes = problem.boundaryTypes(element, scvf); + const auto& bcTypes = problem.boundaryTypes(fvGeometry.element(), scvf); // Dirichlet boundaries if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann()) @@ -173,9 +213,7 @@ protected: // Neumann and Robin ("solution dependent Neumann") boundary conditions else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet()) { - auto neumannFluxes = problem.neumann(element, fvGeometry, evv, efvc, scvf); - - // multiply neumann fluxes with the area and the extrusion factor + auto neumannFluxes = getUnscaledNeumannFluxes_(fvGeometry.element(), scvf); neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor(); r[localDofIdx] += neumannFluxes; } @@ -190,14 +228,10 @@ protected: template = 0> void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const { + const auto& fvGeometry = context_.elementGridGeometry(); const auto& evv = context_.elementVariables().elemVolVars(); const auto& problem = evv.gridVolVars().problem(); - // TODO: Modify problem interfaces to receive context - const auto& fvGeometry = context_.elementGridGeometry(); - const auto& element = fvGeometry.element(); - const auto& efvc = context_.elementVariables().elemFluxVarsCache(); - // inner faces if (!scvf.boundary()) { @@ -212,12 +246,12 @@ protected: else { const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); - const auto& bcTypes = problem.boundaryTypes(element, scv); + const auto& bcTypes = problem.boundaryTypes(fvGeometry.element(), scv); // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions. if (bcTypes.hasNeumann()) { - const auto neumannFluxes = problem.neumann(element, fvGeometry, evv, efvc, scvf); + const auto neumannFluxes = getUnscaledNeumannFluxes_(fvGeometry.element(), scvf); const auto area = Extrusion::area(scvf)*evv[scv].extrusionFactor(); // only add fluxes to equations for which Neumann is set @@ -228,6 +262,18 @@ protected: } } + //! get the user-specified Neumann boundary flux + NumEqVector getUnscaledNeumannFluxes_(const Element& element, + const SubControlVolumeFace& scvf) const + { + // TODO: Modify problem interfaces to receive context + const auto& fvGeometry = context_.elementGridGeometry(); + const auto& evv = context_.elementVariables().elemVolVars(); + const auto& efvc = context_.elementVariables().elemFluxVarsCache(); + const auto& problem = evv.gridVolVars().problem(); + return problem.neumann(element, fvGeometry, evv, efvc, scvf); + } + private: const LocalContext& context_; }; -- GitLab From 239af063e97a0a4eb0e9d110379cc6fd37f6f8ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 15:36:59 +0100 Subject: [PATCH 13/32] [experimental][assembly] add cc local assembler --- dumux/assembly/fv/CMakeLists.txt | 1 + dumux/assembly/fv/cclocalassembler.hh | 477 ++++++++++++++++++++++++++ 2 files changed, 478 insertions(+) create mode 100644 dumux/assembly/fv/cclocalassembler.hh diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt index 40d6910f65..890727a926 100644 --- a/dumux/assembly/fv/CMakeLists.txt +++ b/dumux/assembly/fv/CMakeLists.txt @@ -1,5 +1,6 @@ install(FILES boxlocalassembler.hh +cclocalassembler.hh localoperator.hh operators.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv) diff --git a/dumux/assembly/fv/cclocalassembler.hh b/dumux/assembly/fv/cclocalassembler.hh new file mode 100644 index 0000000000..6210660368 --- /dev/null +++ b/dumux/assembly/fv/cclocalassembler.hh @@ -0,0 +1,477 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief An assembler for Jacobian and residual contribution + * per element for cell-centered finite volume schemes. + */ +#ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH +#define DUMUX_CC_LOCAL_ASSEMBLER_HH + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + + +namespace Dumux::Experimental { + +/*! + * \ingroup Assembly + * \brief An assembler for Jacobian and residual contribution + * per element for cell-centered finite volume schemes. + * \tparam LO The element-local operator + */ +template +class CCLocalAssembler +{ + using Operators = typename LO::Operators; + using LocalContext = typename LO::LocalContext; + using ElementVariables = typename LocalContext::ElementVariables; + using ElementGridGeometry = typename LocalContext::ElementGridGeometry; + + using GG = typename ElementGridGeometry::GridGeometry; + using GV = typename ElementVariables::GridVariables; + + using FVElementGeometry = typename GG::LocalView; + using SubControlVolume = typename GG::SubControlVolume; + using Element = typename GG::GridView::template Codim<0>::Entity; + + using PrimaryVariables = typename GV::PrimaryVariables; + using NumEqVector = Dumux::NumEqVector; + using Scalar = typename GV::Scalar; + + static constexpr int numEq = NumEqVectorTraits::numEq; + static constexpr int maxElementStencilSize = GG::maxElementStencilSize; + static constexpr bool enableGridFluxVarsCache = GV::GridFluxVariablesCache::cachingEnabled; + +public: + + //! the parameters of a stage in time integration + using StageParams = MultiStageParams; + + //! export the grid variables type on which to operate + using GridVariables = GV; + + //! export the underlying local operator + using LocalOperator = LO; + + //! export the vector storing the residuals of all dofs of the element + using ElementResidualVector = typename LO::ElementResidualVector; + + /*! + * \brief Constructor for stationary problems. + */ + explicit CCLocalAssembler(const Element& element, + GridVariables& gridVariables, + DiffMethod dm = DiffMethod::numeric) + : diffMethod_(dm) + , gridVariables_(gridVariables) + , fvGeometry_(localView(gridVariables.gridGeometry())) + , elemVariables_(localView(gridVariables)) + , prevElemVariables_() + , elementIsGhost_(element.partitionType() == Dune::GhostEntity) + , stageParams_(nullptr) + { + if (diffMethod_ != DiffMethod::numeric) + DUNE_THROW(Dune::NotImplemented, "Provided differentiation method"); + + fvGeometry_.bind(element); + elemVariables_.bind(element, fvGeometry_); + } + + /*! + * \brief Constructor for instationary problems. + * \note Using this constructor, we assemble one stage within + * a time integration step using multi-stage methods. + */ + explicit CCLocalAssembler(const Element& element, + const std::vector& prevGridVars, + GridVariables& gridVariables, + std::shared_ptr stageParams, + DiffMethod dm = DiffMethod::numeric) + : diffMethod_(dm) + , gridVariables_(gridVariables) + , fvGeometry_(localView(gridVariables.gridGeometry())) + , elemVariables_(localView(gridVariables)) + , prevElemVariables_() + , elementIsGhost_(element.partitionType() == Dune::GhostEntity) + , stageParams_(stageParams) + { + if (diffMethod_ != DiffMethod::numeric) + DUNE_THROW(Dune::NotImplemented, "Provided differentiation method"); + + if (prevGridVars.size() != stageParams->size() - 1) + DUNE_THROW(Dune::InvalidStateException, "Grid Variables for all stages needed"); + + fvGeometry_.bind(element); + for (const auto& gridVars : prevGridVars) + { + prevElemVariables_.emplace_back(localView(gridVars)); + prevElemVariables_.back().bind(element, fvGeometry_); + } + + elemVariables_.bind(element, fvGeometry_); + } + + /*! + * \brief Computes the derivatives with respect to the given element and adds + * them to the global matrix. The element residual is written into the + * right hand side. + */ + template + void assembleJacobianAndResidual(JacobianMatrix& jac, + ResidualVector& res, + const PartialReassembler* partialReassembler = nullptr) + { + const auto& element = fvGeometry_.element(); + const auto globalI = fvGeometry_.gridGeometry().elementMapper().index(element); + if (partialReassembler + && partialReassembler->elementColor(globalI) == EntityColor::green) + { + res[globalI] = evalLocalResidual()[0]; // forward to the internal implementation + } + else + { + res[globalI] = assembleJacobianAndResidual_(jac); // forward to the internal implementation + } + } + + + /*! + * \brief Computes the derivatives with respect to the dofs of the given + * element and adds them to the given jacobian matrix. + */ + template + void assembleJacobian(JacobianMatrix& jac) + { + assembleJacobianAndResidual_(jac); + } + + //! Assemble the residual only + template + void assembleResidual(ResidualVector& res) + { + const auto residual = evalLocalResidual(); + for (const auto& scv : scvs(fvGeometry_)) + res[scv.dofIndex()] += residual[scv.localDofIndex()]; + } + + //! Evaluate the complete local residual for the current element. + ElementResidualVector evalLocalResidual() const + { + if (isStationary()) + { + const auto context = makeLocalContext(fvGeometry_, elemVariables_); + LocalOperator localOperator(context); + return elementIsGhost_ ? localOperator.getEmptyResidual() + : localOperator.evalFluxesAndSources(); + } + else + { + ElementResidualVector residual(fvGeometry_.numScv()); + residual = 0.0; + + if (!elementIsGhost_) + { + // add the terms associated with previous stages + for (std::size_t k = 0; k < stageParams_->size()-1; ++k) + addStageTerms_(residual, k, + makeLocalContext(fvGeometry_, prevElemVariables_[k])); + + // add the terms of the current stage + addStageTerms_(residual, stageParams_->size()-1, + makeLocalContext(fvGeometry_, elemVariables_)); + } + + return residual; + } + } + + //! Return true if a stationary problem is assembled + bool isStationary() const + { return !stageParams_; } + +protected: + + //! add the terms of a stage to the current element residual + void addStageTerms_(ElementResidualVector& r, + std::size_t stageIdx, + const LocalContext& context) const + { + LocalOperator localOperator(context); + if (!stageParams_->skipTemporal(stageIdx)) + r.axpy(stageParams_->temporalWeight(stageIdx), + localOperator.evalStorage()); + if (!stageParams_->skipSpatial(stageIdx)) + r.axpy(stageParams_->spatialWeight(stageIdx), + localOperator.evalFluxesAndSources()); + } + + /*! + * \brief Computes the derivatives with respect to the dofs of the given + * element and adds them to the global matrix. + * \return The element residual at the current solution. + */ + template + NumEqVector assembleJacobianAndResidual_(JacobianMatrix& A, + const PartialReassembler* partialReassembler = nullptr) + { + // TODO: DO we need this to be constexpr? + if (diffMethod_ == DiffMethod::numeric) + return assembleJacobianAndResidualNumeric_(A, partialReassembler); + else + DUNE_THROW(Dune::NotImplemented, "Analytic assembler for box"); + } + + /*! + * \brief Computes the derivatives by means of numeric differentiation + * and adds them to the global matrix. + * \return The element residual at the current solution. + * \note This specialization is for the box scheme with numeric differentiation + */ + template + NumEqVector assembleJacobianAndResidualNumeric_(JacobianMatrix& A, + const PartialReassembler* partialReassembler = nullptr) + { + // alias for the variables of the current stage + auto& curVariables = elemVariables_; + auto& curElemVolVars = curVariables.elemVolVars(); + auto& curElemFluxVarsCache = curVariables.elemFluxVarsCache(); + const auto& x = curVariables.gridVariables().dofs(); + + // get stencil informations + const auto& element = fvGeometry_.element(); + const auto& gridGeometry = fvGeometry_.gridGeometry(); + const auto& connectivityMap = gridGeometry.connectivityMap(); + const auto globalI = gridGeometry.elementMapper().index(element); + const auto numNeighbors = connectivityMap[globalI].size(); + + // deduce the problem type + const auto& problem = curElemVolVars.gridVolVars().problem(); + using Problem = std::decay_t; + + // assemble the undeflected residual + using Residuals = ReservedBlockVector; + Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0; + origResiduals[0] = evalLocalResidual()[0]; + + // lambda for convenient evaluation of the fluxes across scvfs in the neighbors + auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf) + { + const auto context = makeLocalContext(fvGeometry_, elemVariables_); + return LocalOperator(context).computeFlux(neighbor, scvf); + }; + + // get the elements in which we need to evaluate the fluxes + // and calculate these in the undeflected state + Dune::ReservedVector neighborElements; + neighborElements.resize(numNeighbors); + + unsigned int j = 1; + for (const auto& dataJ : connectivityMap[globalI]) + { + neighborElements[j-1] = gridGeometry.element(dataJ.globalJ); + + if (neighborElements[j-1].partitionType() != Dune::GhostEntity) + for (const auto scvfIdx : dataJ.scvfsJ) + origResiduals[j] += evalNeighborFlux(neighborElements[j-1], + fvGeometry_.scvf(scvfIdx)); + + ++j; + } + + // reference to the element's scv (needed later) and corresponding vol vars + const auto& scv = fvGeometry_.scv(globalI); + auto& curVolVars = getVolVarAcess_(curElemVolVars, gridVariables_.gridVolVars(), scv); + + // save a copy of the original privars and vol vars in order + // to restore the original solution after deflection + const auto origPriVars = x[globalI]; + const auto origVolVars = curVolVars; + + // element solution container to be deflected + auto elemSol = elementSolution(element, x, gridGeometry); + + // derivatives in the neighbors with repect to the current elements + // in index 0 we save the derivative of the element residual with respect to it's own dofs + Residuals partialDerivs(numNeighbors + 1); + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) + { + partialDerivs = 0.0; + + auto evalResiduals = [&](Scalar priVar) + { + Residuals partialDerivsTmp(numNeighbors + 1); + partialDerivsTmp = 0.0; + + // update the volume variables and the flux var cache + elemSol[0][pvIdx] = priVar; + curVolVars.update(elemSol, problem, element, scv); + curElemFluxVarsCache.update(element, fvGeometry_, curElemVolVars); + if (enableGridFluxVarsCache) + gridVariables_.gridFluxVarsCache().updateElement(element, fvGeometry_, curElemVolVars); + + // calculate the residual with the deflected primary variables + partialDerivsTmp[0] = evalLocalResidual()[0]; + + // calculate the fluxes in the neighbors with the deflected primary variables + for (std::size_t k = 0; k < numNeighbors; ++k) + if (neighborElements[k].partitionType() != Dune::GhostEntity) + for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ) + partialDerivsTmp[k+1] += evalNeighborFlux(neighborElements[k], + fvGeometry_.scvf(scvfIdx)); + + return partialDerivsTmp; + }; + + // derive the residuals numerically + static const NumericEpsilon numEps{problem.paramGroup()}; + static const int numDiffMethod = getParamFromGroup(problem.paramGroup(), "Assembly.NumericDifferenceMethod"); + + const auto eps = numEps(elemSol[0][pvIdx], pvIdx); + NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], + partialDerivs, origResiduals, + eps, numDiffMethod); + + // Correct derivative for ghost elements, i.e. set a 1 for the derivative w.r.t. the + // current primary variable and a 0 elsewhere. As we always solve for a delta of the + // solution with repect to the initial one, this results in a delta of 0 for ghosts. + if (elementIsGhost_) + { + partialDerivs[0] = 0.0; + partialDerivs[0][pvIdx] = 1.0; + } + + // For instationary simulations, scale the coupling + // fluxes of the current stage correctly + if (stageParams_) + { + for (std::size_t k = 0; k < numNeighbors; ++k) + partialDerivs[k+1] *= stageParams_->spatialWeight(stageParams_->size()-1); + } + + // add the current partial derivatives to the global jacobian matrix + // no special treatment is needed if globalJ is a ghost because then derivatives have been assembled to 0 above + if constexpr (Problem::enableInternalDirichletConstraints()) + { + // check if own element has internal Dirichlet constraint + const auto internalDirichletConstraintsOwnElement = problem.hasInternalDirichletConstraint(element, scv); + const auto dirichletValues = problem.internalDirichlet(element, scv); + + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (internalDirichletConstraintsOwnElement[eqIdx]) + { + origResiduals[0][eqIdx] = origVolVars.priVars()[eqIdx] - dirichletValues[eqIdx]; + A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0; + } + else + A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx]; + } + + // off-diagonal entries + j = 1; + for (const auto& dataJ : connectivityMap[globalI]) + { + const auto& neighborElement = gridGeometry.element(dataJ.globalJ); + const auto& neighborScv = fvGeometry_.scv(dataJ.globalJ); + const auto internalDirichletConstraintsNeighbor = problem().hasInternalDirichletConstraint(neighborElement, neighborScv); + + for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) + { + if (internalDirichletConstraintsNeighbor[eqIdx]) + A[dataJ.globalJ][globalI][eqIdx][pvIdx] = 0.0; + else + A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j][eqIdx]; + } + + ++j; + } + } + else // no internal Dirichlet constraints specified + { + for (int eqIdx = 0; eqIdx < numEq; eqIdx++) + { + // the diagonal entries + A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx]; + + // off-diagonal entries + j = 1; + for (const auto& dataJ : connectivityMap[globalI]) + A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx]; + } + } + + // restore the original state of the scv's volume variables + curVolVars = origVolVars; + + // restore the current element solution + elemSol[0][pvIdx] = origPriVars[pvIdx]; + } + + // restore original state of the flux vars cache in case of global caching. + // This has to be done in order to guarantee that everything is in an undeflected + // state before the assembly of another element is called. In the case of local caching + // this is obsolete because the elemFluxVarsCache used here goes out of scope after this. + // We only have to do this for the last primary variable, for all others the flux var cache + // is updated with the correct element volume variables before residual evaluations + curElemFluxVarsCache.update(element, fvGeometry_, curElemVolVars); + if (enableGridFluxVarsCache) + gridVariables_.gridFluxVarsCache().updateElement(element, fvGeometry_, curElemVolVars); + + // return the original residual + return origResiduals[0]; + } + + //! Returns the volume variables of the local view (case of no caching) + template = 0> + auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv) + { return elemVolVars[scv]; } + + //! Returns the volume variables of the grid vol vars (case of global caching) + template = 0> + auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv) + { return gridVolVars.volVars(scv); } + + DiffMethod diffMethod_; //!< the type of differentiation method + GridVariables& gridVariables_; //!< reference to the grid variables + FVElementGeometry fvGeometry_; //!< element-local finite volume geometry + ElementVariables elemVariables_; //!< element variables of the current stage + std::vector prevElemVariables_; //!< element variables of prior stages + + bool elementIsGhost_; + std::shared_ptr stageParams_; +}; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 69a6cc584e707a5e47e19d77266bf969349eb09f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 15:37:21 +0100 Subject: [PATCH 14/32] [experimental][assembly] add convenience alias for local assembler --- dumux/assembly/CMakeLists.txt | 1 + dumux/assembly/localassembler.hh | 65 ++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+) create mode 100644 dumux/assembly/localassembler.hh diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt index 17a49d7b47..63a13d1e43 100644 --- a/dumux/assembly/CMakeLists.txt +++ b/dumux/assembly/CMakeLists.txt @@ -11,6 +11,7 @@ fvlocalassemblerbase.hh fvlocalresidual.hh initialsolution.hh jacobianpattern.hh +localassembler.hh numericepsilon.hh partialreassembler.hh staggeredfvassembler.hh diff --git a/dumux/assembly/localassembler.hh b/dumux/assembly/localassembler.hh new file mode 100644 index 0000000000..5e51dfd9f4 --- /dev/null +++ b/dumux/assembly/localassembler.hh @@ -0,0 +1,65 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief Convenience alias to select a local assembler based on the discretization scheme. + */ +#ifndef DUMUX_LOCAL_ASSEMBLER_HH +#define DUMUX_LOCAL_ASSEMBLER_HH + +#include +#include "fv/boxlocalassembler.hh" +#include "fv/cclocalassembler.hh" + +namespace Dumux::Experimental { +namespace Detail { + + template + struct LocalAssemblerChooser; + + template + struct LocalAssemblerChooser + { using type = BoxLocalAssembler; }; + + template + struct LocalAssemblerChooser + { using type = CCLocalAssembler; }; + + template + struct LocalAssemblerChooser + { using type = CCLocalAssembler; }; + + template + using LocalAssemblerType + = typename LocalAssemblerChooser::type; + +} // end namespace Detail + +/*! + * \ingroup Assembly + * \brief Helper alias to select the local assembler type from a local operator. + */ +template +using LocalAssembler = Detail::LocalAssemblerType; + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 0fe1aa6f966ace7dbc89a7d1f5fecba2fad4c867 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Tue, 9 Mar 2021 16:30:31 +0100 Subject: [PATCH 15/32] [test][1p][compressible][tpfa] use new time integration --- .../1p/compressible/instationary/CMakeLists.txt | 11 +++++++++++ .../1p/compressible/instationary/main_experimental.cc | 10 +++++----- 2 files changed, 16 insertions(+), 5 deletions(-) diff --git a/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt b/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt index 05d17f3dcf..ed0822004d 100644 --- a/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt +++ b/test/porousmediumflow/1p/compressible/instationary/CMakeLists.txt @@ -41,3 +41,14 @@ dumux_add_test(NAME test_1p_compressible_instationary_box_experimental --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_box-reference.vtu ${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_box_experimental-00010.vtu --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_box_experimental params.input -Problem.Name test_1p_compressible_instationary_box_experimental") + +dumux_add_test(NAME test_1p_compressible_instationary_tpfa_experimental + LABELS porousmediumflow 1p experimental + SOURCES main_experimental.cc + COMPILE_DEFINITIONS TYPETAG=OnePCompressibleTpfa + COMPILE_DEFINITIONS EXPERIMENTAL=1 + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_1p_cc-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_tpfa_experimental-00010.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_1p_compressible_instationary_tpfa_experimental params.input -Problem.Name test_1p_compressible_instationary_tpfa_experimental") diff --git a/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc b/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc index 4b2281feb5..f1c7d683aa 100644 --- a/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc +++ b/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc @@ -19,7 +19,7 @@ /*! * \file * \ingroup OnePTests - * \brief Test for the one-phase model + * \brief test for the one-phase model */ #include #include @@ -44,7 +44,7 @@ #include #include -#include +#include #include #include @@ -123,7 +123,7 @@ int main(int argc, char** argv) timeLoop->setMaxTimeStepSize(maxDt); // use implicit Euler for time integration - auto timeMethod = std::make_shared>(); + auto timeMethod = std::make_shared>(); // the assembler (we use the immiscible operators to define the system of equations) using ModelTraits = GetPropType; @@ -132,7 +132,7 @@ int main(int argc, char** argv) using Operators = Experimental::FVImmiscibleOperators; using LocalOperator = Experimental::FVLocalOperator; - using LocalAssembler = Experimental::BoxLocalAssembler; + using LocalAssembler = Experimental::LocalAssembler; using Assembler = Experimental::Assembler; auto assembler = std::make_shared(gridGeometry, DiffMethod::numeric); @@ -145,7 +145,7 @@ int main(int argc, char** argv) auto nonLinearSolver = std::make_shared(assembler, linearSolver); // the time stepper for time integration - using TimeStepper = MultiStageTimeStepper; + using TimeStepper = Experimental::MultiStageTimeStepper; TimeStepper timeStepper(nonLinearSolver, timeMethod); // set some check points for the time loop -- GitLab From 017b723d248b3291e8a2b1f71b3f05363316872c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 15 Mar 2021 18:42:08 +0100 Subject: [PATCH 16/32] [operators] hide context definition from user --- dumux/assembly/fv/operators.hh | 16 ++++++++++++++-- dumux/porousmediumflow/immiscible/operators.hh | 16 ++++++++++------ .../porousmediumflow/nonisothermal/operators.hh | 11 +++++++---- .../instationary/main_experimental.cc | 4 ++-- 4 files changed, 33 insertions(+), 14 deletions(-) diff --git a/dumux/assembly/fv/operators.hh b/dumux/assembly/fv/operators.hh index f2103086aa..7e07999ee3 100644 --- a/dumux/assembly/fv/operators.hh +++ b/dumux/assembly/fv/operators.hh @@ -27,18 +27,30 @@ #include #include +#include namespace Dumux::Experimental { +/*! + * \ingroup Assembly + * \brief Convenience alias to define the context finite-volume operators work on. + * \tparam EV The element-(stencil)-local variables + */ +template +using FVOperatorsContext = DefaultLocalContext; + /*! * \ingroup Assembly * \brief The base class for the sub-control entity-local evaluation of * the terms of equations in the context of finite-volume schemes - * \tparam LC the element-stencil-local data required to evaluate the terms + * \tparam EV The element-(stencil)-local variables */ -template +template class FVOperators { + // context type on which to operate + using LC = FVOperatorsContext; + // The grid geometry on which the scheme operates using FVElementGeometry = typename LC::ElementGridGeometry; using GridGeometry = typename FVElementGeometry::GridGeometry; diff --git a/dumux/porousmediumflow/immiscible/operators.hh b/dumux/porousmediumflow/immiscible/operators.hh index 50f2a69fb8..3c2cd9718c 100644 --- a/dumux/porousmediumflow/immiscible/operators.hh +++ b/dumux/porousmediumflow/immiscible/operators.hh @@ -38,7 +38,7 @@ namespace Dumux::Experimental { * \tparam ModelTraits defines model-related types and variables (e.g. number of phases) * \tparam FluxVariables the type that is responsible for computing the individual * flux contributions, i.e., advective, diffusive, convective... - * \tparam LocalContext the element-stencil-local data required to evaluate the terms + * \tparam ElemVars The element-(stencil)-local variables * \tparam EnergyOperators optional template argument, specifying the class that * handles the operators related to non-isothermal effects. * The default energy operators are compatible with isothermal @@ -46,15 +46,16 @@ namespace Dumux::Experimental { */ template> + class ElemVars, + class EnergyOperators = FVNonIsothermalOperators>> class FVImmiscibleOperators -: public FVOperators +: public FVOperators { - using ParentType = FVOperators; + using ParentType = FVOperators; + using LC = typename ParentType::LocalContext; // The variables required for the evaluation of the equation - using ElementVariables = typename LocalContext::ElementVariables; + using ElementVariables = typename LC::ElementVariables; using GridVariables = typename ElementVariables::GridVariables; using VolumeVariables = typename GridVariables::VolumeVariables; using ElementVolumeVariables = typename ElementVariables::ElementVolumeVariables; @@ -78,6 +79,9 @@ class FVImmiscibleOperators static constexpr bool enableEnergyBalance = ModelTraits::enableEnergyBalance();; public: + //! export the type of context on which this class operates + using LocalContext = LC; + //! export the type used to store scalar values for all equations using typename ParentType::NumEqVector; diff --git a/dumux/porousmediumflow/nonisothermal/operators.hh b/dumux/porousmediumflow/nonisothermal/operators.hh index 1949e41397..0f15b79bf1 100644 --- a/dumux/porousmediumflow/nonisothermal/operators.hh +++ b/dumux/porousmediumflow/nonisothermal/operators.hh @@ -34,20 +34,23 @@ namespace Dumux::Experimental { * involved in the system of equations of non-isothermal models * that assume thermal equilibrium between all phases. * \tparam ModelTraits Model-specific traits. - * \tparam LocalContext Element-local context (geometry & primary/secondary variables) + * \tparam LC Element-local context (geometry & primary/secondary variables) */ -template +template class FVNonIsothermalOperators { // The variables required for the evaluation of the equation - using ElementVariables = typename LocalContext::ElementVariables; + using ElementVariables = typename LC::ElementVariables; // The grid geometry on which the scheme operates - using GridGeometry = typename LocalContext::ElementGridGeometry::GridGeometry; + using GridGeometry = typename LC::ElementGridGeometry::GridGeometry; using SubControlVolume = typename GridGeometry::SubControlVolume; public: + //! export the type of context on which this class operates + using LocalContext = LC; + /*! * \brief The energy storage in the fluid phase with index phaseIdx * diff --git a/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc b/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc index f1c7d683aa..ec2ed71b0c 100644 --- a/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc +++ b/test/porousmediumflow/1p/compressible/instationary/main_experimental.cc @@ -128,8 +128,8 @@ int main(int argc, char** argv) // the assembler (we use the immiscible operators to define the system of equations) using ModelTraits = GetPropType; using FluxVariables = GetPropType; - using LocalContext = Experimental::LocalContext; - using Operators = Experimental::FVImmiscibleOperators; + using ElementVariables = typename GridVariables::LocalView; + using Operators = Experimental::FVImmiscibleOperators; using LocalOperator = Experimental::FVLocalOperator; using LocalAssembler = Experimental::LocalAssembler; -- GitLab From 203afe2c7eecd6ef3262d9bea52c85bd748ad17c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 3 Mar 2021 16:21:17 +0100 Subject: [PATCH 17/32] [disc] add solution state class --- dumux/discretization/CMakeLists.txt | 1 + dumux/discretization/solutionstate.hh | 161 ++++++++++++++++++++++++++ 2 files changed, 162 insertions(+) create mode 100644 dumux/discretization/solutionstate.hh diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt index 625fe8b3cc..9519315c93 100644 --- a/dumux/discretization/CMakeLists.txt +++ b/dumux/discretization/CMakeLists.txt @@ -28,6 +28,7 @@ rotationsymmetricgridgeometrytraits.hh rotationsymmetricscv.hh rotationsymmetricscvf.hh scvandscvfiterators.hh +solutionstate.hh staggered.hh subcontrolvolumebase.hh subcontrolvolumefacebase.hh diff --git a/dumux/discretization/solutionstate.hh b/dumux/discretization/solutionstate.hh new file mode 100644 index 0000000000..55aa657f88 --- /dev/null +++ b/dumux/discretization/solutionstate.hh @@ -0,0 +1,161 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Discretization + * \brief State classes to represent the global and element-local states of the + * primary variables of a numerical model, consisting of the spatial degrees + * of freedom and potentially a corresponding time level. + */ +#ifndef DUMUX_SOLUTION_STATE_HH +#define DUMUX_SOLUTION_STATE_HH + +#include +#include + +namespace Dumux::Experimental { + +namespace Concept { + +//! Concept for solution states +struct SolutionState +{ + template + auto require(const T& t) -> decltype( + t.dofs(), + t.timeLevel() + ); +}; + +//! Concept for element solution states +struct ElementSolutionState +{ + template + auto require(const T& t) -> decltype( + t.elementSolution(), + t.timeLevel() + ); +}; + +} // end namespace Concept + +/*! + * \ingroup Discretization + * \brief State class to represent the state of the primary variables of a + * numerical model, consisting of the spatial degrees of freedom and + * potentially a corresponding time level. + */ +template +class SolutionState +{ +public: + using SolutionVector = SV; + using TimeLevel = TL; + + //! Constructor + SolutionState(const SolutionVector* dofs, + const TimeLevel* timeLevel) + : dofs_(dofs) + , timeLevel_(timeLevel) + {} + + //! return the solution vector + const SolutionVector& dofs() const + { return *dofs_; } + + //! return the time level + const TimeLevel& timeLevel() const + { return *timeLevel_; } + +private: + const SolutionVector* dofs_; + const TimeLevel* timeLevel_; +}; + +/*! + * \ingroup Discretization + * \brief State class to represent the element-local state of the primary + * variables of a numerical model, consisting of the spatial degrees + * of freedom and potentially a corresponding time level. + */ +template +class ElementSolutionState +{ +public: + using ElementSolution = ES; + using TimeLevel = TL; + + //! Constructor + ElementSolutionState(ElementSolution&& es, + const TimeLevel* timeLevel) + : elemSol_(std::move(es)) + , timeLevel_(timeLevel) + {} + + //! return the element solution vector + const ElementSolution& elementSolution() const + { return elemSol_; } + + //! return the element solution vector + ElementSolution& elementSolution() + { return elemSol_; } + + //! return the time level + const TimeLevel& timeLevel() const + { return *timeLevel_; } + +private: + ElementSolution elemSol_; + const TimeLevel* timeLevel_; +}; + +/*! + * \ingroup Discretization + * \brief Convenience factory function to make an element solution state + * from a global solution state. + */ +template +auto makeElementSolutionState(const Element& element, + const SolutionState& state, + const GridGeometry& gridGeom) +{ + static_assert(Dune::models(), + "Provided type does not fulfill solution state interface"); + + auto elemSol = elementSolution(element, state.dofs(), gridGeom); + return ElementSolutionState{std::move(elemSol), &state.timeLevel()}; +} + +/*! + * \ingroup Discretization + * \brief Convenience factory function to make an element solution state + * from grid variables. + */ +template +auto makeElementSolutionState(const Element& element, + const GridVariables& gridVariables) +{ + const auto& gridGeom = gridVariables.gridGeometry(); + auto elemSol = elementSolution(element, gridVariables.dofs(), gridGeom); + return ElementSolutionState{std::move(elemSol), &gridVariables.timeLevel()}; +} + +} // end namespace Dumux::Experimental + +#endif -- GitLab From 68a898a1a8d9b186c887ce1120f9559813090410 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 24 Feb 2021 14:41:37 +0100 Subject: [PATCH 18/32] [experimental] fvproblem with new interfaces --- dumux/common/fvproblem.hh | 116 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index f306356947..32b3eb1b2a 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -601,6 +601,122 @@ private: PointSourceMap pointSourceMap_; }; +namespace Experimental { + +//! Experimental problem implementation compatible with new time +//! integration schemes and corresponding assembly. +template +class FVProblem : public Dumux::FVProblem +{ + using ParentType = Dumux::FVProblem; + + using GridGeometry = GetPropType; + using FVElementGeometry = typename GridGeometry::LocalView; + using SubControlVolume = typename GridGeometry::SubControlVolume; + using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + + using GridView = typename GridGeometry::GridView; + using Element = typename GridView::template Codim<0>::Entity; + + using PrimaryVariables = typename ParentType::Traits::PrimaryVariables; + using NumEqVector = typename ParentType::Traits::NumEqVector; + using Scalar = typename ParentType::Traits::Scalar; + + static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box; + +public: + // pull up constructor + using ParentType::ParentType; + + /*! + * \brief Specify the kind of boundary conditions used on a discrete + * entity on the boundary. + * + * \param element The finite element + * \param boundaryEntity The boundary entity (scv/scvf) + * \note In cell-centered schemes, boundaryEntity is a sub-control + * volume face (scvf). In the box scheme, a sub-control volume (scv). + */ + template + auto boundaryTypes(const Element& element, + const BoundaryEntity& boundaryEntity) const + { + if constexpr (isBox) + return this->asImp_().boundaryTypesAtPos(boundaryEntity.dofPosition()); + else + return this->asImp_().boundaryTypesAtPos(boundaryEntity.ipGlobal()); + } + + /*! + * \brief Evaluate the boundary conditions for a discrete entity on the boundary. + * + * \param element The finite element + * \param boundaryEntity The boundary entity (scv/scvf) + * \param elementSolution The element-local state of the discrete solution. + * \note In cell-centered schemes, boundaryEntity is a sub-control + * volume face (scvf). In the box scheme, a sub-control volume (scv). + */ + template + PrimaryVariables dirichlet(const Element& element, + const BoundaryEntity& boundaryEntity, + const ElementSolution& elementSolution) const + { + if constexpr (isBox) + return this->asImp_().dirichletAtPos(boundaryEntity.dofPosition()); + else + return this->asImp_().dirichletAtPos(boundaryEntity.ipGlobal()); + } + + /*! + * \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 + * + * \param context The element-local context + * \param scvf The sub control volume face + * \note The element-local context consists of an element and local + * geometric information, together with the primary/secondary + * variables in that local scope. Potentially, users can define + * the context to carry an additional object containing further + * locally required data. + */ + template + NumEqVector neumann(const LocalContext& context, + const SubControlVolumeFace& scvf) const + { return this->asImp_().neumannAtPos(scvf.ipGlobal()); } + + /*! + * \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 context The element-local context + * \param scv The sub control volume + * \note The element-local context consists of an element and local + * geometric information, together with the primary/secondary + * variables in that local scope. Potentially, users can define + * the context to carry an additional object containing further + * locally required data. + * + * For this method, the return 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$. + */ + template + NumEqVector source(const LocalContext& context, + const SubControlVolume& scv) const + { return this->asImp_().sourceAtPos(scv.center()); } + + //! TODO: Point sources! +}; + +} // end namespace Experimental } // end namespace Dumux #endif -- GitLab From f292b193789499abfac2599aa207c22e32158660 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Fri, 26 Feb 2021 15:44:36 +0100 Subject: [PATCH 19/32] [experimental][pmflow] introduce experimental problem --- dumux/porousmediumflow/problem.hh | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/dumux/porousmediumflow/problem.hh b/dumux/porousmediumflow/problem.hh index 801f846340..298e5b038d 100644 --- a/dumux/porousmediumflow/problem.hh +++ b/dumux/porousmediumflow/problem.hh @@ -34,10 +34,11 @@ namespace Dumux { * * TODO: derive from base problem property? */ -template -class PorousMediumFlowProblem : public FVProblem +template> +class PorousMediumFlowProblem : public BaseProblem { - using ParentType = FVProblem; + using ParentType = BaseProblem; using GridGeometry = GetPropType; using GridView = typename GetPropType::GridView; @@ -139,6 +140,15 @@ protected: std::shared_ptr spatialParams_; }; +namespace Experimental { + +//! Experimental problem implementation compatible with new time integration +//! schemes and corresponding assembly. +template +using PorousMediumFlowProblem + = Dumux::PorousMediumFlowProblem>; + +} // end namespace Experimental } // end namespace Dumux #endif -- GitLab From 25afc8c6ad5dac9048cd3c7d6af87becad1e63bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 10 Mar 2021 16:35:41 +0100 Subject: [PATCH 20/32] [assembly] add helpers for compatibility with experimental --- dumux/assembly/CMakeLists.txt | 1 + dumux/assembly/experimentalhelpers.hh | 135 ++++++++++++++++++++++++++ 2 files changed, 136 insertions(+) create mode 100644 dumux/assembly/experimentalhelpers.hh diff --git a/dumux/assembly/CMakeLists.txt b/dumux/assembly/CMakeLists.txt index 63a13d1e43..4288fe602d 100644 --- a/dumux/assembly/CMakeLists.txt +++ b/dumux/assembly/CMakeLists.txt @@ -6,6 +6,7 @@ cclocalassembler.hh cclocalresidual.hh diffmethod.hh entitycolor.hh +experimentalhelpers.hh fvassembler.hh fvlocalassemblerbase.hh fvlocalresidual.hh diff --git a/dumux/assembly/experimentalhelpers.hh b/dumux/assembly/experimentalhelpers.hh new file mode 100644 index 0000000000..4a7e5c859d --- /dev/null +++ b/dumux/assembly/experimentalhelpers.hh @@ -0,0 +1,135 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Assembly + * \brief Helper classes and functions to be used in places where compatibility + * between the current default and the new experimental assembly is seeked. + */ +#ifndef DUMUX_ASSEMBLY_EXPERIMENTAL_HELPERS_HH +#define DUMUX_ASSEMBLY_EXPERIMENTAL_HELPERS_HH + +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace Dumux::Experimental { +namespace CompatibilityHelpers { + +// implements necessary interfaces for compatibility +template +class GridVariablesFacade +{ + using P = typename ElemVolVars::GridVolumeVariables::Problem; +public: + using GridVolumeVariables = typename ElemVolVars::GridVolumeVariables; + using GridFluxVariablesCache = typename ElemFluxVarsCache::GridFluxVariablesCache; + using GridGeometry = typename ProblemTraits

::GridGeometry; +}; + +// implements necessary interfaces for compatibility +template +class ElemVariablesFacade +{ + const ElemVolVars* elemVolVars_; + const ElemFluxVarsCache* elemFluxVarsCache_; + +public: + + using GridVariables = GridVariablesFacade; + + ElemVariablesFacade(const ElemVolVars* elemVolVars, + const ElemFluxVarsCache* elemFluxVarsCache) + : elemVolVars_(elemVolVars) + , elemFluxVarsCache_(elemFluxVarsCache) + {} + + const ElemVolVars& elemVolVars() const { return *elemVolVars_; } + const ElemFluxVarsCache& elemFluxVarsCache() const { return *elemFluxVarsCache_; } +}; + +// helper function to construct an element variables facade from given local views +template +ElemVariablesFacade +makeElemVariablesFacade(const ElemVolVars& elemVolVars, + const ElemFluxVarsCache& elemFluxVarsCache) +{ return {&elemVolVars, &elemFluxVarsCache}; } + +// Get an element solution state from a global state. +// If a solution vector is passed (current default style), +// then an element solution is returned, whereas if an actual +// solution state is given (new concept), we return an elem sol state. +template +auto getElemSolState(const Element& element, + const SolState& solState, + const GridGeometry& gridGeometry) +{ + if constexpr (Dune::models()) + return makeElementSolutionState(element, solState, gridGeometry); + else + return elementSolution(element, solState, gridGeometry); +} + +// Get an element solution state from grid variables. +// If grid variables with the experimental layout are passed, +// we return an actual element solution state. Otherwise, +// we return an empty element solution. +template +auto getElemSolState(const Element& element, + const GridVariables& gridVariables) +{ + if constexpr (areExperimentalGridVars) + return makeElementSolutionState(element, gridVariables); + else + { + using DummySolVec = std::vector>; + using ElemSol = decltype(elementSolution(element, + std::declval(), + gridVariables.gridGeometry())); + return ElemSol{}; + } +} + +// get the Dirichlet values for a boundary entity from a user problem +template +decltype(auto) getDirichletValues(const Problem& problem, + const Element& element, + const BoundaryEntity& boundaryEntity, + const ElemSolState& elemSolState) +{ + // given element solution state models an actual state + if constexpr (Dune::models()) + return problem.dirichlet(element, boundaryEntity, elemSolState); + + // current default style (elemSolState = element solution, not passed to user interface) + else + return problem.dirichlet(element, boundaryEntity); +} + +} // end namespace CompatibilityHelpers +} // end namespace Dumux::Experimental + +#endif -- GitLab From 9cebf9a9de8ea9b3d3e070be6d06f6ef5f4ad765 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 10 Mar 2021 16:42:32 +0100 Subject: [PATCH 21/32] [pmflow] compatibility with experimental assembly --- dumux/porousmediumflow/volumevariables.hh | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/dumux/porousmediumflow/volumevariables.hh b/dumux/porousmediumflow/volumevariables.hh index bfca5829db..51452884b5 100644 --- a/dumux/porousmediumflow/volumevariables.hh +++ b/dumux/porousmediumflow/volumevariables.hh @@ -26,6 +26,8 @@ #ifndef DUMUX_POROUSMEDIUMFLOW_VOLUME_VARIABLES_HH #define DUMUX_POROUSMEDIUMFLOW_VOLUME_VARIABLES_HH +#include + namespace Dumux { /*! @@ -66,8 +68,19 @@ public: const Element& element, const Scv& scv) { - priVars_ = elemSol[scv.localDofIndex()]; - extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol); + // compatibility layer with new experimental assembly style (elem sol = element solution state) + if constexpr (Dune::models()) + { + const auto& elemSolState = elemSol; + priVars_ = elemSolState.elementSolution()[scv.localDofIndex()]; + extrusionFactor_ = problem.spatialParams().extrusionFactor(element, scv, elemSolState); + } + // current standard style (element solution-based) + else + { + priVars_ = elemSol[scv.localDofIndex()]; + extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol); + } } /*! -- GitLab From 40e51e77e1844785fb496fc028c9941edf9d6518 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 10 Mar 2021 16:43:01 +0100 Subject: [PATCH 22/32] [nonisothermal] compatibility with experimental assembly --- .../nonisothermal/volumevariables.hh | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/dumux/porousmediumflow/nonisothermal/volumevariables.hh b/dumux/porousmediumflow/nonisothermal/volumevariables.hh index a8594028e4..ebeaedddd3 100644 --- a/dumux/porousmediumflow/nonisothermal/volumevariables.hh +++ b/dumux/porousmediumflow/nonisothermal/volumevariables.hh @@ -176,10 +176,19 @@ public: FluidState& fluidState, SolidState& solidState) { + // compatibility with new experimental assembly style + const auto& priVars = [&elemSol, &scv] () + { + if constexpr (Dune::models()) + return elemSol.elementSolution()[scv.localDofIndex()]; + else + return elemSol[scv.localDofIndex()]; + } (); + if constexpr (fullThermalEquilibrium) { // retrieve temperature from solution vector, all phases have the same temperature - const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx]; + const Scalar T = priVars[temperatureIdx]; for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { fluidState.setTemperature(phaseIdx, T); @@ -192,7 +201,7 @@ public: // this means we have 1 temp for fluid phase, one for solid if constexpr (fluidThermalEquilibrium) { - const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx]; + const Scalar T = priVars[temperatureIdx]; for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { fluidState.setTemperature(phaseIdx, T); @@ -204,11 +213,11 @@ public: for(int phaseIdx=0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { // retrieve temperatures from solution vector, phases might have different temperature - const Scalar T = elemSol[scv.localDofIndex()][temperatureIdx + phaseIdx]; + const Scalar T = priVars[temperatureIdx + phaseIdx]; fluidState.setTemperature(phaseIdx, T); } } - const Scalar solidTemperature = elemSol[scv.localDofIndex()][temperatureIdx+numEnergyEq-1]; + const Scalar solidTemperature = priVars[temperatureIdx+numEnergyEq-1]; solidState.setTemperature(solidTemperature); } } -- GitLab From f79f696ae73d1f650a8c4dae190ff573e796dfc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 10 Mar 2021 16:44:14 +0100 Subject: [PATCH 23/32] [1p] compatibility with experimental assembly --- dumux/porousmediumflow/1p/volumevariables.hh | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/dumux/porousmediumflow/1p/volumevariables.hh b/dumux/porousmediumflow/1p/volumevariables.hh index 88da7efd0f..03469aa48a 100644 --- a/dumux/porousmediumflow/1p/volumevariables.hh +++ b/dumux/porousmediumflow/1p/volumevariables.hh @@ -25,8 +25,11 @@ #ifndef DUMUX_1P_VOLUME_VARIABLES_HH #define DUMUX_1P_VOLUME_VARIABLES_HH +#include + #include #include + #include #include @@ -51,6 +54,7 @@ class OnePVolumeVariables using Scalar = typename Traits::PrimaryVariables::value_type; using PermeabilityType = typename Traits::PermeabilityType; static constexpr int numFluidComps = ParentType::numFluidComponents(); + public: //! Export the underlying fluid system using FluidSystem = typename Traits::FluidSystem; @@ -109,10 +113,18 @@ public: FluidState& fluidState, SolidState& solidState) { + // compatibility with new experimental assembly style + const auto& priVars = [&elemSol, &scv] () + { + if constexpr (Dune::models()) + return elemSol.elementSolution()[scv.localDofIndex()]; + else + return elemSol[scv.localDofIndex()]; + } (); + EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState); fluidState.setSaturation(/*phaseIdx=*/0, 1.); - const auto& priVars = elemSol[scv.localDofIndex()]; fluidState.setPressure(/*phaseIdx=*/0, priVars[Indices::pressureIdx]); // saturation in a single phase is always 1 and thus redundant -- GitLab From c834045b42b5e2c59383b935b5778fd9aaf71d83 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 10 Mar 2021 16:46:14 +0100 Subject: [PATCH 24/32] [io] compatibility with experimental assembly --- dumux/io/vtkoutputmodule.hh | 37 ++++++++++++++++++++++++++---- dumux/porousmediumflow/velocity.hh | 19 ++++++++++++++- 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh index faec1248d4..530b5ed243 100644 --- a/dumux/io/vtkoutputmodule.hh +++ b/dumux/io/vtkoutputmodule.hh @@ -46,6 +46,7 @@ #include #include +#include #include "vtkfunction.hh" #include "velocityoutput.hh" @@ -466,12 +467,26 @@ private: if (velocityOutput_->enableOutput()) { fvGeometry.bind(element); - elemVolVars.bind(element, fvGeometry, sol_); + if constexpr (Experimental::areExperimentalGridVars) + { + const Experimental::SolutionState solState(&gridVariables_.dofs(), + &gridVariables_.timeLevel()); + elemVolVars.bind(element, fvGeometry, solState); + } + else + elemVolVars.bind(element, fvGeometry, sol_); } else { fvGeometry.bindElement(element); - elemVolVars.bindElement(element, fvGeometry, sol_); + if constexpr (Experimental::areExperimentalGridVars) + { + const Experimental::SolutionState solState(&gridVariables_.dofs(), + &gridVariables_.timeLevel()); + elemVolVars.bindElement(element, fvGeometry, solState); + } + else + elemVolVars.bindElement(element, fvGeometry, sol_); } if (!volVarScalarDataInfo_.empty() @@ -659,12 +674,26 @@ private: if (velocityOutput_->enableOutput()) { fvGeometry.bind(element); - elemVolVars.bind(element, fvGeometry, sol_); + if constexpr (Experimental::areExperimentalGridVars) + { + const Experimental::SolutionState solState(&gridVariables_.dofs(), + &gridVariables_.timeLevel()); + elemVolVars.bind(element, fvGeometry, solState); + } + else + elemVolVars.bind(element, fvGeometry, sol_); } else { fvGeometry.bindElement(element); - elemVolVars.bindElement(element, fvGeometry, sol_); + if constexpr (Experimental::areExperimentalGridVars) + { + const Experimental::SolutionState solState(&gridVariables_.dofs(), + &gridVariables_.timeLevel()); + elemVolVars.bindElement(element, fvGeometry, solState); + } + else + elemVolVars.bindElement(element, fvGeometry, sol_); } if (!volVarScalarDataInfo_.empty() diff --git a/dumux/porousmediumflow/velocity.hh b/dumux/porousmediumflow/velocity.hh index 088e626cbb..464783284d 100644 --- a/dumux/porousmediumflow/velocity.hh +++ b/dumux/porousmediumflow/velocity.hh @@ -39,6 +39,10 @@ #include #include +// for compatibility with experimental assembly +#include +#include + namespace Dumux { #ifndef DOXYGEN @@ -327,7 +331,20 @@ public: else { // check if we have Neumann no flow, we can just use 0 - const auto neumannFlux = problem_.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); + const auto neumannFlux = [&] () + { + // compatibility layer with new experimental assembly + if constexpr (Experimental::areExperimentalGridVars) + { + using namespace Experimental::CompatibilityHelpers; + const auto elemVarsFacade = makeElemVariablesFacade(elemVolVars, elemFluxVarsCache); + const auto context = makeLocalContext(fvGeometry, elemVarsFacade); + return problem_.neumann(context, scvf); + } + else + return problem_.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf); + } (); + using NumEqVector = std::decay_t; if (Dune::FloatCmp::eq(neumannFlux, NumEqVector(0.0), 1e-30)) scvfFluxes[scvfIndexInInside[localScvfIdx]] = 0; -- GitLab From 01783a8d2a20463cbf860bfb6e3053ffc3842d4a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Wed, 10 Mar 2021 16:47:03 +0100 Subject: [PATCH 25/32] [experimentalassembly] use solution state --- dumux/assembly/fv/boxlocalassembler.hh | 15 ++-- dumux/assembly/fv/cclocalassembler.hh | 9 +- dumux/assembly/fv/localoperator.hh | 18 +--- dumux/assembly/fv/operators.hh | 11 +-- .../box/elementvolumevariables.hh | 51 +++++++---- .../mpfa/elementvolumevariables.hh | 90 ++++++++++++------- .../tpfa/elementvolumevariables.hh | 66 ++++++++------ dumux/discretization/fvgridvariables.hh | 15 +++- .../solidstates/updatesolidvolumefractions.hh | 16 +++- .../1p/compressible/instationary/problem.hh | 16 +++- 10 files changed, 194 insertions(+), 113 deletions(-) diff --git a/dumux/assembly/fv/boxlocalassembler.hh b/dumux/assembly/fv/boxlocalassembler.hh index 9e49482ca1..a6bbb28ca5 100644 --- a/dumux/assembly/fv/boxlocalassembler.hh +++ b/dumux/assembly/fv/boxlocalassembler.hh @@ -35,6 +35,7 @@ #include #include +#include #include @@ -268,14 +269,16 @@ public: template< typename ApplyDirichletFunctionType > void evalDirichletBoundaries(ApplyDirichletFunctionType applyDirichlet) { + const auto& element = fvGeometry_.element(); const auto& problem = gridVariables_.gridVolVars().problem(); + const auto elemSolState = makeElementSolutionState(element, gridVariables_); for (const auto& scvI : scvs(fvGeometry_)) { - const auto bcTypes = problem.boundaryTypes(fvGeometry_.element(), scvI); + const auto bcTypes = problem.boundaryTypes(element, scvI); if (bcTypes.hasDirichlet()) { - const auto dirichletValues = problem.dirichlet(fvGeometry_.element(), scvI); + const auto dirichletValues = problem.dirichlet(element, scvI, elemSolState); // set the Dirichlet conditions in residual and jacobian for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) @@ -371,13 +374,13 @@ protected: // alias for the variables of the current stage auto& curVariables = elemVariables_; auto& curElemVolVars = curVariables.elemVolVars(); - const auto& x = curVariables.gridVariables().dofs(); const auto& element = fvGeometry_.element(); const auto& problem = curElemVolVars.gridVolVars().problem(); const auto origResiduals = evalLocalResidual(); - const auto origElemSol = elementSolution(element, x, fvGeometry_.gridGeometry()); - auto elemSol = origElemSol; + auto elemSolState = makeElementSolutionState(element, curVariables.gridVariables()); + auto& elemSol = elemSolState.elementSolution(); + const auto origElemSol = elemSol; ////////////////////////////////////////////////////////////////////////////////////////////// // Calculate derivatives of the residual of all dofs in element with respect to themselves. // @@ -401,7 +404,7 @@ protected: { // update the volume variables and compute element residual elemSol[scvI.localDofIndex()][pvIdx] = priVar; - curVolVars.update(elemSol, problem, element, scvI); + curVolVars.update(elemSolState, problem, element, scvI); return evalLocalResidual(); }; diff --git a/dumux/assembly/fv/cclocalassembler.hh b/dumux/assembly/fv/cclocalassembler.hh index 6210660368..75bd41587e 100644 --- a/dumux/assembly/fv/cclocalassembler.hh +++ b/dumux/assembly/fv/cclocalassembler.hh @@ -37,6 +37,7 @@ #include #include #include +#include #include @@ -318,8 +319,10 @@ protected: const auto origPriVars = x[globalI]; const auto origVolVars = curVolVars; - // element solution container to be deflected - auto elemSol = elementSolution(element, x, gridGeometry); + // element solution to be deflected + auto elemSolState = makeElementSolutionState(element, curVariables.gridVariables()); + auto& elemSol = elemSolState.elementSolution(); + const auto origElemSol = elemSolState.elementSolution(); // derivatives in the neighbors with repect to the current elements // in index 0 we save the derivative of the element residual with respect to it's own dofs @@ -335,7 +338,7 @@ protected: // update the volume variables and the flux var cache elemSol[0][pvIdx] = priVar; - curVolVars.update(elemSol, problem, element, scv); + curVolVars.update(elemSolState, problem, element, scv); curElemFluxVarsCache.update(element, fvGeometry_, curElemVolVars); if (enableGridFluxVarsCache) gridVariables_.gridFluxVarsCache().updateElement(element, fvGeometry_, curElemVolVars); diff --git a/dumux/assembly/fv/localoperator.hh b/dumux/assembly/fv/localoperator.hh index 55483f5fe3..55cf9e9ad3 100644 --- a/dumux/assembly/fv/localoperator.hh +++ b/dumux/assembly/fv/localoperator.hh @@ -182,7 +182,7 @@ public: } // compute and return the Neumann boundary fluxes - auto neumannFluxes = getUnscaledNeumannFluxes_(element, scvf); + auto neumannFluxes = problem.neumann(context_, scvf); neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor(); return neumannFluxes; } @@ -213,7 +213,7 @@ protected: // Neumann and Robin ("solution dependent Neumann") boundary conditions else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet()) { - auto neumannFluxes = getUnscaledNeumannFluxes_(fvGeometry.element(), scvf); + auto neumannFluxes = problem.neumann(context_, scvf); neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor(); r[localDofIdx] += neumannFluxes; } @@ -251,7 +251,7 @@ protected: // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions. if (bcTypes.hasNeumann()) { - const auto neumannFluxes = getUnscaledNeumannFluxes_(fvGeometry.element(), scvf); + const auto neumannFluxes = problem.neumann(context_, scvf); const auto area = Extrusion::area(scvf)*evv[scv].extrusionFactor(); // only add fluxes to equations for which Neumann is set @@ -262,18 +262,6 @@ protected: } } - //! get the user-specified Neumann boundary flux - NumEqVector getUnscaledNeumannFluxes_(const Element& element, - const SubControlVolumeFace& scvf) const - { - // TODO: Modify problem interfaces to receive context - const auto& fvGeometry = context_.elementGridGeometry(); - const auto& evv = context_.elementVariables().elemVolVars(); - const auto& efvc = context_.elementVariables().elemFluxVarsCache(); - const auto& problem = evv.gridVolVars().problem(); - return problem.neumann(element, fvGeometry, evv, efvc, scvf); - } - private: const LocalContext& context_; }; diff --git a/dumux/assembly/fv/operators.hh b/dumux/assembly/fv/operators.hh index 7e07999ee3..867f3e50c0 100644 --- a/dumux/assembly/fv/operators.hh +++ b/dumux/assembly/fv/operators.hh @@ -122,18 +122,15 @@ public: { SourceTerm source(0.0); - // TODO: The problem-interfaces should be adapted to receive context? - const auto& fvGeometry = context.elementGridGeometry(); - const auto& elemVolVars = context.elementVariables().elemVolVars(); - const auto& element = fvGeometry.element(); - // add contributions from volume flux sources - source += problem.source(element, fvGeometry, elemVolVars, scv); + source += problem.source(context, scv); // add contribution from possible point sources - source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv); + // TODO: Point sources + // source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv); // multiply with scv volume + const auto& elemVolVars = context.elementVariables().elemVolVars(); source *= Extrusion::volume(scv)*elemVolVars[scv].extrusionFactor(); return source; diff --git a/dumux/discretization/box/elementvolumevariables.hh b/dumux/discretization/box/elementvolumevariables.hh index 2ddb9591a6..3ebbdb9096 100644 --- a/dumux/discretization/box/elementvolumevariables.hh +++ b/dumux/discretization/box/elementvolumevariables.hh @@ -27,6 +27,7 @@ #include #include +#include namespace Dumux { @@ -68,19 +69,19 @@ public: // For compatibility reasons with the case of not storing the vol vars. // function to be called before assembling an element, preparing the vol vars within the stencil - template + template void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { - bindElement(element, fvGeometry, sol); + bindElement(element, fvGeometry, solState); } // function to prepare the vol vars within the element - template + template void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { eIdx_ = fvGeometry.gridGeometry().elementMapper().index(element); } @@ -114,27 +115,43 @@ public: : gridVolVarsPtr_(&gridVolVars) {} // specialization for box models, simply forwards to the bindElement method - template + template void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { - bindElement(element, fvGeometry, sol); + bindElement(element, fvGeometry, solState); } // specialization for box models - template + template void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { - // get the solution at the dofs of the element - auto elemSol = elementSolution(element, sol, fvGeometry.gridGeometry()); - - // resize volume variables to the required size - volumeVariables_.resize(fvGeometry.numScv()); - for (auto&& scv : scvs(fvGeometry)) - volumeVariables_[scv.indexInElement()].update(elemSol, gridVolVars().problem(), element, scv); + // compatibility layer with new experimental assembly style + if constexpr (Dune::models()) + { + const auto elemSolState = makeElementSolutionState(element, + solState, + fvGeometry.gridGeometry()); + + // resize volume variables to the required size + volumeVariables_.resize(fvGeometry.numScv()); + for (auto&& scv : scvs(fvGeometry)) + volumeVariables_[scv.indexInElement()].update(elemSolState, gridVolVars().problem(), element, scv); + } + // current standard style (solution state = solution vector) + else + { + // get the solution at the dofs of the element + auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry()); + + // resize volume variables to the required size + volumeVariables_.resize(fvGeometry.numScv()); + for (auto&& scv : scvs(fvGeometry)) + volumeVariables_[scv.indexInElement()].update(elemSol, gridVolVars().problem(), element, scv); + } } const VolumeVariables& operator [](std::size_t scvIdx) const diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index 83b6b3a93a..bfa6015505 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -29,6 +29,7 @@ #include #include +#include #include namespace Dumux { @@ -76,14 +77,16 @@ namespace CCMpfa { * \param element The element to which the finite volume geometry is bound * \param fvGeometry The element finite volume geometry * \param nodalIndexSet The dual grid index set around a node + * \param solState The state of the discrete solution */ - template + template void addBoundaryVolVarsAtNode(std::vector& volVars, std::vector& volVarIndices, const Problem& problem, const typename FVElemGeom::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElemGeom& fvGeometry, - const NodalIndexSet& nodalIndexSet) + const NodalIndexSet& nodalIndexSet, + const SolState& solState) { if (nodalIndexSet.numBoundaryScvfs() == 0) return; @@ -91,6 +94,11 @@ namespace CCMpfa { // index of the element the fvGeometry was bound to const auto boundElemIdx = fvGeometry.gridGeometry().elementMapper().index(element); + // compatibility layer with new assembly: get element-local state + // This causes a slight overhead when using the old assembly (elem sol not needed) + using namespace Experimental::CompatibilityHelpers; + auto elemSolState = getElemSolState(element, solState, fvGeometry.gridGeometry()); + // check each scvf in the index set for boundary presence for (auto scvfIdx : nodalIndexSet.gridScvfIndices()) { @@ -108,11 +116,15 @@ namespace CCMpfa { // boundaries the "outside" vol vars cannot be properly defined. if (bcTypes.hasOnlyDirichlet()) { + // update the element solution state with the Dirichlet values + if constexpr (Dune::models()) + elemSolState.elementSolution()[0] = getDirichletValues(problem, insideElement, ivScvf, elemSolState); + // old assembly style (elemSolState = elementSolution) + else + elemSolState[0] = getDirichletValues(problem, insideElement, ivScvf, elemSolState); + VolumeVariables dirichletVolVars; - dirichletVolVars.update(elementSolution(problem.dirichlet(insideElement, ivScvf)), - problem, - insideElement, - fvGeometry.scv(insideScvIdx)); + dirichletVolVars.update(elemSolState, problem, insideElement, fvGeometry.scv(insideScvIdx)); volVars.emplace_back(std::move(dirichletVolVars)); volVarIndices.push_back(ivScvf.outsideScvIdx()); @@ -130,16 +142,23 @@ namespace CCMpfa { * \param problem The problem containing the Dirichlet boundary conditions * \param element The element to which the finite volume geometry was bound * \param fvGeometry The element finite volume geometry + * \param solState The state of the discrete solution */ - template + template void addBoundaryVolVars(std::vector& volVars, std::vector& volVarIndices, const Problem& problem, const typename FVElemGeom::GridGeometry::GridView::template Codim<0>::Entity& element, - const FVElemGeom& fvGeometry) + const FVElemGeom& fvGeometry, + const SolState& solState) { const auto& gridGeometry = fvGeometry.gridGeometry(); + // compatibility layer with new assembly: get element-local state + // This causes a slight overhead when using old assembly (elem sol not needed) + using namespace Experimental::CompatibilityHelpers; + auto elemSolState = getElemSolState(element, solState, fvGeometry.gridGeometry()); + // treat the BCs inside the element if (fvGeometry.hasBoundaryScvf()) { @@ -155,11 +174,15 @@ namespace CCMpfa { // boundaries the "outside" vol vars cannot be properly defined. if (problem.boundaryTypes(element, scvf).hasOnlyDirichlet()) { + // update the element solution state with the Dirichlet values + if constexpr (Dune::models()) + elemSolState.elementSolution()[0] = getDirichletValues(problem, element, scvf, elemSolState); + // old assembly style (elemSolState = elementSolution) + else + elemSolState[0] = getDirichletValues(problem, element, scvf, elemSolState); + VolumeVariables dirichletVolVars; - dirichletVolVars.update(elementSolution(problem.dirichlet(element, scvf)), - problem, - element, - scvI); + dirichletVolVars.update(elemSolState, problem, element, scvI); volVars.emplace_back(std::move(dirichletVolVars)); volVarIndices.push_back(scvf.outsideScvIdx()); @@ -173,10 +196,10 @@ namespace CCMpfa { { if (!gridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex())) addBoundaryVolVarsAtNode( volVars, volVarIndices, problem, element, fvGeometry, - gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet() ); + gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet(), solState ); else addBoundaryVolVarsAtNode( volVars, volVarIndices, problem, element, fvGeometry, - gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet() ); + gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet(), solState ); } } } // end namespace CCMpfa @@ -228,10 +251,10 @@ public: } //! precompute all volume variables in a stencil of an element - bind Dirichlet vol vars in the stencil - template + template void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { clear(); @@ -241,7 +264,8 @@ public: { boundaryVolVars_.reserve(maxNumBoundaryVolVars); boundaryVolVarIndices_.reserve(maxNumBoundaryVolVars); - CCMpfa::addBoundaryVolVars(boundaryVolVars_, boundaryVolVarIndices_, gridVolVars().problem(), element, fvGeometry); + CCMpfa::addBoundaryVolVars(boundaryVolVars_, boundaryVolVarIndices_, + gridVolVars().problem(), element, fvGeometry, solState); } } @@ -299,10 +323,10 @@ public: : gridVolVarsPtr_(&gridVolVars) {} //! Prepares the volume variables within the element stencil - template + template void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { clear(); @@ -319,12 +343,13 @@ public: volumeVariables_.reserve(numVolVars+maxNumBoundaryVolVars); volVarIndices_.reserve(numVolVars+maxNumBoundaryVolVars); + // compatibility layer with new assembly: get element-local state + using namespace Experimental::CompatibilityHelpers; + const auto elemSolState = getElemSolState(element, solState, gridGeometry); + VolumeVariables volVars; const auto& scvI = fvGeometry.scv(globalI); - volVars.update(elementSolution(element, sol, gridGeometry), - problem, - element, - scvI); + volVars.update(elemSolState, problem, element, scvI); volVarIndices_.push_back(scvI.dofIndex()); volumeVariables_.emplace_back(std::move(volVars)); @@ -333,12 +358,10 @@ public: for (auto&& dataJ : assemblyMapI) { const auto& elementJ = gridGeometry.element(dataJ.globalJ); + const auto elemSolStateJ = getElemSolState(elementJ, solState, gridGeometry); const auto& scvJ = fvGeometry.scv(dataJ.globalJ); VolumeVariables volVarsJ; - volVarsJ.update(elementSolution(elementJ, sol, gridGeometry), - problem, - elementJ, - scvJ); + volVarsJ.update(elemSolStateJ, problem, elementJ, scvJ); volVarIndices_.push_back(scvJ.dofIndex()); volumeVariables_.emplace_back(std::move(volVarsJ)); @@ -346,7 +369,8 @@ public: // maybe prepare boundary volume variables if (maxNumBoundaryVolVars > 0) - CCMpfa::addBoundaryVolVars(volumeVariables_, volVarIndices_, problem, element, fvGeometry); + CCMpfa::addBoundaryVolVars(volumeVariables_, volVarIndices_, + problem, element, fvGeometry, solState); // //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends // //! on additional DOFs not included in the discretization schemes' occupation pattern @@ -373,10 +397,10 @@ public: } //! Prepares the volume variables of an element - template + template void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { clear(); @@ -385,9 +409,13 @@ public: volumeVariables_.resize(1); volVarIndices_.resize(1); + // compatibility layer with new assembly: get element-local state + using namespace Experimental::CompatibilityHelpers; + const auto elemSolState = getElemSolState(element, solState, gridGeometry); + // update the volume variables of the element const auto& scv = fvGeometry.scv(eIdx); - volumeVariables_[0].update(elementSolution(element, sol, gridGeometry), + volumeVariables_[0].update(elemSolState, gridVolVars().problem(), element, scv); diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh index e10d37ac5c..e64257ba6f 100644 --- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh @@ -28,6 +28,7 @@ #include #include +#include #include namespace Dumux { @@ -84,10 +85,10 @@ public: } //! precompute all boundary volume variables in a stencil of an element, the remaining ones are cached - template + template void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solutionState) { if (!fvGeometry.hasBoundaryScvf()) return; @@ -96,6 +97,11 @@ public: boundaryVolVarIndices_.reserve(fvGeometry.numScvf()); boundaryVolumeVariables_.reserve(fvGeometry.numScvf()); + // compatibility layer with new assembly: get element-local state + // This causes a slight overhead when using the old assembly (elem sol not needed) + using namespace Experimental::CompatibilityHelpers; + auto elemSolState = getElemSolState(element, solutionState, fvGeometry.gridGeometry()); + for (const auto& scvf : scvfs(fvGeometry)) { if (!scvf.boundary()) @@ -106,14 +112,16 @@ public: const auto bcTypes = problem.boundaryTypes(element, scvf); if (bcTypes.hasOnlyDirichlet()) { - const auto dirichletPriVars = elementSolution(problem.dirichlet(element, scvf)); - auto&& scvI = fvGeometry.scv(scvf.insideScvIdx()); - + // update the element solution state with the Dirichlet values + if constexpr (Dune::models()) + elemSolState.elementSolution()[0] = getDirichletValues(problem, element, scvf, elemSolState); + // old assembly style (elemSolState = elementSolution) + else + elemSolState[0] = getDirichletValues(problem, element, scvf, elemSolState); + + const auto& scvI = fvGeometry.scv(scvf.insideScvIdx()); VolumeVariables volVars; - volVars.update(dirichletPriVars, - problem, - element, - scvI); + volVars.update(elemSolState, problem, element, scvI); boundaryVolumeVariables_.emplace_back(std::move(volVars)); boundaryVolVarIndices_.push_back(scvf.outsideScvIdx()); @@ -174,10 +182,10 @@ public: : gridVolVarsPtr_(&gridVolVars) {} //! Prepares the volume variables within the element stencil - template + template void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { clear_(); @@ -193,8 +201,10 @@ public: int localIdx = 0; // update the volume variables of the element at hand - auto&& scvI = fvGeometry.scv(globalI); - volumeVariables_[localIdx].update(elementSolution(element, sol, gridGeometry), + using namespace Experimental::CompatibilityHelpers; + auto elemSolState = getElemSolState(element, solState, gridGeometry); + const auto& scvI = fvGeometry.scv(globalI); + volumeVariables_[localIdx].update(elemSolState, problem, element, scvI); @@ -206,7 +216,7 @@ public: { const auto& elementJ = gridGeometry.element(dataJ.globalJ); auto&& scvJ = fvGeometry.scv(dataJ.globalJ); - volumeVariables_[localIdx].update(elementSolution(elementJ, sol, gridGeometry), + volumeVariables_[localIdx].update(getElemSolState(elementJ, solState, gridGeometry), problem, elementJ, scvJ); @@ -227,18 +237,20 @@ public: const auto bcTypes = problem.boundaryTypes(element, scvf); if (bcTypes.hasOnlyDirichlet()) { - const auto dirichletPriVars = elementSolution(problem.dirichlet(element, scvf)); + // update the element solution state with the Dirichlet values + if constexpr (Dune::models()) + elemSolState.elementSolution()[0] = getDirichletValues(problem, element, scvf, elemSolState); + // old assembly style (elemSolState = elementSolution) + else + elemSolState[0] = getDirichletValues(problem, element, scvf, elemSolState); volumeVariables_.resize(localIdx+1); volVarIndices_.resize(localIdx+1); - volumeVariables_[localIdx].update(dirichletPriVars, - problem, - element, - scvI); - volVarIndices_[localIdx] = scvf.outsideScvIdx(); - ++localIdx; - } + volumeVariables_[localIdx].update(elemSolState, problem, element, scvI); + volVarIndices_[localIdx] = scvf.outsideScvIdx(); + ++localIdx; } + } } //! Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends @@ -263,10 +275,10 @@ public: // } } - template + template void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element, const FVElementGeometry& fvGeometry, - const SolutionVector& sol) + const SolutionState& solState) { clear_(); @@ -275,8 +287,10 @@ public: volVarIndices_.resize(1); // update the volume variables of the element - auto&& scv = fvGeometry.scv(eIdx); - volumeVariables_[0].update(elementSolution(element, sol, fvGeometry.gridGeometry()), + using namespace Experimental::CompatibilityHelpers; + const auto elemSolState = getElemSolState(element, solState, fvGeometry.gridGeometry()); + const auto& scv = fvGeometry.scv(eIdx); + volumeVariables_[0].update(elemSolState, gridVolVars().problem(), element, scv); diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh index d0ea810c50..e4b5a2d974 100644 --- a/dumux/discretization/fvgridvariables.hh +++ b/dumux/discretization/fvgridvariables.hh @@ -179,6 +179,7 @@ private: #include #include +#include namespace Dumux::Experimental { @@ -221,8 +222,9 @@ public: void bind(const Element& element, const FVElementGeometry& fvGeometry) { - const auto& x = gridVariables().dofs(); - elemVolVars_.bind(element, fvGeometry, x); + const Experimental::SolutionState solState(&gridVariables().dofs(), + &gridVariables().timeLevel()); + elemVolVars_.bind(element, fvGeometry, solState); elemFluxVarsCache_.bind(element, fvGeometry, elemVolVars_); } @@ -234,7 +236,13 @@ public: void bindElemVolVars(const Element& element, const FVElementGeometry& fvGeometry) { - elemVolVars_.bind(element, fvGeometry, gridVariables().dofs()); + Dumux::Experimental::SolutionState solState(&gridVariables().dofs(), + &gridVariables().timeLevel()); + + if (bindEntireStencil) + elemVolVars_.bind(element, fvGeometry, solState); + else + elemVolVars_.bindElement(element, fvGeometry, solState); // unbind flux variables cache elemFluxVarsCache_ = localView(gridVariables().gridFluxVarsCache()); @@ -377,6 +385,7 @@ private: GridFluxVariablesCache gridFluxVarsCache_; //!< the flux variables cache }; + // implementation details namespace Detail { struct hasGridVolVars { diff --git a/dumux/material/solidstates/updatesolidvolumefractions.hh b/dumux/material/solidstates/updatesolidvolumefractions.hh index 0e3a185341..2f327be0b9 100644 --- a/dumux/material/solidstates/updatesolidvolumefractions.hh +++ b/dumux/material/solidstates/updatesolidvolumefractions.hh @@ -19,11 +19,13 @@ /*! * \file * \ingroup SolidStates - * \brief Update the solid volume fractions (inert and reacitve) and set them in the solidstate + * \brief Update the solid volume fractions (inert and reactive) and set them in the solidstate */ #ifndef DUMUX_UPDATE_SOLID_VOLUME_FRACTION_HH #define DUMUX_UPDATE_SOLID_VOLUME_FRACTION_HH +#include + namespace Dumux { /*! @@ -50,8 +52,16 @@ void updateSolidVolumeFractions(const ElemSol& elemSol, if (!(solidState.isInert())) { - auto&& priVars = elemSol[scv.localDofIndex()]; - for (int sCompIdx = 0; sCompIdx < solidState.numComponents- solidState.numInertComponents; ++sCompIdx) + // compatibility with new experimental assembly style + const auto& priVars = [&] () + { + if constexpr (Dune::models()) + return elemSol.elementSolution()[scv.localDofIndex()]; + else + return elemSol[scv.localDofIndex()]; + } (); + + for (int sCompIdx = 0; sCompIdx < solidState.numComponents-solidState.numInertComponents; ++sCompIdx) solidState.setVolumeFraction(sCompIdx, priVars[solidVolFracOffset + sCompIdx]); } } diff --git a/test/porousmediumflow/1p/compressible/instationary/problem.hh b/test/porousmediumflow/1p/compressible/instationary/problem.hh index 561252815b..d8003bb8d4 100644 --- a/test/porousmediumflow/1p/compressible/instationary/problem.hh +++ b/test/porousmediumflow/1p/compressible/instationary/problem.hh @@ -33,7 +33,19 @@ #include #include + +#ifndef EXPERIMENTAL +#define EXPERIMENTAL false +#endif + +#if EXPERIMENTAL +#define BASEPROBLEM Experimental::PorousMediumFlowProblem +#else +#define BASEPROBLEM PorousMediumFlowProblem +#endif + namespace Dumux { + /*! * \ingroup OnePTests * \brief Test problem for the compressible one-phase model. @@ -42,9 +54,9 @@ namespace Dumux { * ./test_cc1pfv */ template -class OnePTestProblem : public PorousMediumFlowProblem +class OnePTestProblem : public BASEPROBLEM { - using ParentType = PorousMediumFlowProblem; + using ParentType = BASEPROBLEM; using GridView = typename GetPropType::GridView; using Element = typename GridView::template Codim<0>::Entity; using Scalar = GetPropType; -- GitLab From 36f42c136475eda852b3f69ab2a3c7cf8e0dbfe4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 11 Mar 2021 09:38:12 +0100 Subject: [PATCH 26/32] [assembly][fvoperators] add comment on quasi-newton --- dumux/assembly/fv/operators.hh | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/dumux/assembly/fv/operators.hh b/dumux/assembly/fv/operators.hh index 867f3e50c0..b25a46496c 100644 --- a/dumux/assembly/fv/operators.hh +++ b/dumux/assembly/fv/operators.hh @@ -88,6 +88,13 @@ public: * \param context The element-stencil-local required data * \param scv The sub-control volume * \note This must be overloaded by the implementation + * \note In cell-centered schemes, depending on the assembler or the policy used + * therein, the context may contain variables (e.g. from neighboring cells + * within the stencil) with respect to which the storage term is not + * differentiated. This means that if a Newton solver is used to solve + * the system of equations, one ends up with a quasi-Newton scheme. If + * this is undesireable, make sure to set up an assembler that takes into + * account the dependencies of the storage term w.r.t variables in the context. */ template static StorageTerm storage(const Problem& problem, @@ -114,6 +121,13 @@ public: * \param context The element-stencil-local required data * \param scv The sub-control volume for which the source term is to be computed * \note This is a default implementation forwarding to interfaces in the problem + * \note In cell-centered schemes, depending on the assembler or the policy used + * therein, the context may contain variables (e.g. from neighboring cells + * within the stencil) with respect to which the source term is not + * differentiated. This means that if a Newton solver is used to solve + * the system of equations, one ends up with a quasi-Newton scheme. If + * this is undesireable, make sure to set up an assembler that takes into + * account the dependencies of the storage term w.r.t variables in the context. */ template static SourceTerm source(const Problem& problem, -- GitLab From 95840553a4bb5a99ded5f7214c48e51a461406cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 11 Mar 2021 10:00:28 +0100 Subject: [PATCH 27/32] [problemtraits] add bool for transient dirichlet BCs --- dumux/common/typetraits/problem.hh | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/dumux/common/typetraits/problem.hh b/dumux/common/typetraits/problem.hh index 0131254abd..4de22d482f 100644 --- a/dumux/common/typetraits/problem.hh +++ b/dumux/common/typetraits/problem.hh @@ -25,6 +25,9 @@ #define DUMUX_TYPETRAITS_PROBLEM_HH #include + +#include +#include #include namespace Dumux { @@ -33,6 +36,29 @@ namespace Dumux { namespace Detail { template struct ProblemTraits; + +template +struct hasTransientDirichlet +{ +private: + using TimeLevel = Dumux::Experimental::TimeLevel; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; + using BoundaryEntity = std::conditional_t; +public: + template + auto operator()(const Problem& p) + -> decltype(p.dirichlet(std::declval(), + std::declval(), + std::declval())) + {} +}; + +template +inline constexpr bool hasTransientDirichletInterface + = decltype(isValid(hasTransientDirichlet())(std::declval()))::value; + } // end namespace Detail /*! @@ -44,6 +70,9 @@ struct ProblemTraits { using GridGeometry = std::decay_t().gridGeometry())>; using BoundaryTypes = typename Detail::template ProblemTraits::BoundaryTypes; + + static constexpr bool hasTransientDirichletInterface + = Detail::hasTransientDirichletInterface; }; } // end namespace Dumux -- GitLab From 900779338d979f17b5397c33412821cf4ff7dbf2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 11 Mar 2021 10:13:40 +0100 Subject: [PATCH 28/32] [experimentalassembly] use time level in dirichlet --- dumux/assembly/experimentalhelpers.hh | 7 ++++- dumux/assembly/fv/boxlocalassembler.hh | 37 +++++++++++++++++++------- dumux/common/fvproblem.hh | 19 ++++++++++--- 3 files changed, 49 insertions(+), 14 deletions(-) diff --git a/dumux/assembly/experimentalhelpers.hh b/dumux/assembly/experimentalhelpers.hh index 4a7e5c859d..3003eeae57 100644 --- a/dumux/assembly/experimentalhelpers.hh +++ b/dumux/assembly/experimentalhelpers.hh @@ -122,7 +122,12 @@ decltype(auto) getDirichletValues(const Problem& problem, { // given element solution state models an actual state if constexpr (Dune::models()) - return problem.dirichlet(element, boundaryEntity, elemSolState); + { + if constexpr (ProblemTraits::hasTransientDirichletInterface) + return problem.dirichlet(element, boundaryEntity, elemSolState.timeLevel()); + else + return problem.dirichlet(element, boundaryEntity); + } // current default style (elemSolState = element solution, not passed to user interface) else diff --git a/dumux/assembly/fv/boxlocalassembler.hh b/dumux/assembly/fv/boxlocalassembler.hh index a6bbb28ca5..e54bcc4536 100644 --- a/dumux/assembly/fv/boxlocalassembler.hh +++ b/dumux/assembly/fv/boxlocalassembler.hh @@ -34,6 +34,8 @@ #include #include +#include + #include #include #include @@ -271,14 +273,14 @@ public: { const auto& element = fvGeometry_.element(); const auto& problem = gridVariables_.gridVolVars().problem(); - const auto elemSolState = makeElementSolutionState(element, gridVariables_); for (const auto& scvI : scvs(fvGeometry_)) { const auto bcTypes = problem.boundaryTypes(element, scvI); if (bcTypes.hasDirichlet()) { - const auto dirichletValues = problem.dirichlet(element, scvI, elemSolState); + const auto dirichletValues = getDirichletValues_(problem, element, scvI, + gridVariables_.timeLevel()); // set the Dirichlet conditions in residual and jacobian for (int eqIdx = 0; eqIdx < numEq; ++eqIdx) @@ -457,15 +459,30 @@ protected: return origResiduals; } - //! Returns the volume variables of the local view (case of no caching) - template = 0> - auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv) - { return elemVolVars[scv]; } - //! Returns the volume variables of the grid vol vars (case of global caching) - template = 0> - auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv) - { return gridVolVars.volVars(scv); } + template + auto& getVolVarAcess_(ElemVolVars& elemVolVars, + GridVolVars& gridVolVars, + const SubControlVolume& scv) + { + if constexpr (GridVolVars::cachingEnabled) + return gridVolVars.volVars(scv); + else + return elemVolVars[scv]; + } + + //! get the user-defined Dirichlet boundary conditions + template + auto getDirichletValues_(const Problem& problem, + const Element& element, + const SubControlVolume& scv, + const TimeLevel& timeLevel) + { + if constexpr(ProblemTraits::hasTransientDirichletInterface) + return problem.dirichlet(element, scv, timeLevel); + else + return problem.dirichlet(element, scv); + } DiffMethod diffMethod_; //!< the type of differentiation method GridVariables& gridVariables_; //!< reference to the grid variables diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index 32b3eb1b2a..e4cc354f3a 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -652,14 +652,27 @@ public: * * \param element The finite element * \param boundaryEntity The boundary entity (scv/scvf) - * \param elementSolution The element-local state of the discrete solution. + * \param timeLevel The time level on which to evaluate the boundary conditions * \note In cell-centered schemes, boundaryEntity is a sub-control * volume face (scvf). In the box scheme, a sub-control volume (scv). */ - template + template PrimaryVariables dirichlet(const Element& element, const BoundaryEntity& boundaryEntity, - const ElementSolution& elementSolution) const + const TimeLevel& timeLevel) const + { return this->asImp_().dirichlet(element, boundaryEntity); } + + /*! + * \brief Evaluate the boundary conditions for a discrete entity on the boundary. + * + * \param element The finite element + * \param boundaryEntity The boundary entity (scv/scvf) + * \note In cell-centered schemes, boundaryEntity is a sub-control + * volume face (scvf). In the box scheme, a sub-control volume (scv). + */ + template + PrimaryVariables dirichlet(const Element& element, + const BoundaryEntity& boundaryEntity) const { if constexpr (isBox) return this->asImp_().dirichletAtPos(boundaryEntity.dofPosition()); -- GitLab From 0678fb55d0f856bf6fe528bc22732f89af052299 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 11 Mar 2021 13:14:56 +0100 Subject: [PATCH 29/32] [experimental] define&use extended element solution concept This gets rid of the experimental element solution state concept and instead, this is now represented by an extension of the elementSolution concept we already have. --- dumux/assembly/experimentalhelpers.hh | 45 +----- dumux/assembly/fv/boxlocalassembler.hh | 8 +- dumux/assembly/fv/cclocalassembler.hh | 8 +- .../box/elementvolumevariables.hh | 30 +--- .../mpfa/elementvolumevariables.hh | 49 ++----- .../tpfa/elementvolumevariables.hh | 39 ++---- dumux/discretization/fvgridvariables.hh | 7 +- dumux/discretization/solutionstate.hh | 128 ++++++++++-------- dumux/io/vtkoutputmodule.hh | 12 +- .../solidstates/updatesolidvolumefractions.hh | 12 +- dumux/porousmediumflow/1p/volumevariables.hh | 12 +- .../nonisothermal/volumevariables.hh | 9 +- dumux/porousmediumflow/volumevariables.hh | 17 +-- 13 files changed, 127 insertions(+), 249 deletions(-) diff --git a/dumux/assembly/experimentalhelpers.hh b/dumux/assembly/experimentalhelpers.hh index 3003eeae57..ccec82e842 100644 --- a/dumux/assembly/experimentalhelpers.hh +++ b/dumux/assembly/experimentalhelpers.hh @@ -78,53 +78,18 @@ makeElemVariablesFacade(const ElemVolVars& elemVolVars, const ElemFluxVarsCache& elemFluxVarsCache) { return {&elemVolVars, &elemFluxVarsCache}; } -// Get an element solution state from a global state. -// If a solution vector is passed (current default style), -// then an element solution is returned, whereas if an actual -// solution state is given (new concept), we return an elem sol state. -template -auto getElemSolState(const Element& element, - const SolState& solState, - const GridGeometry& gridGeometry) -{ - if constexpr (Dune::models()) - return makeElementSolutionState(element, solState, gridGeometry); - else - return elementSolution(element, solState, gridGeometry); -} - -// Get an element solution state from grid variables. -// If grid variables with the experimental layout are passed, -// we return an actual element solution state. Otherwise, -// we return an empty element solution. -template -auto getElemSolState(const Element& element, - const GridVariables& gridVariables) -{ - if constexpr (areExperimentalGridVars) - return makeElementSolutionState(element, gridVariables); - else - { - using DummySolVec = std::vector>; - using ElemSol = decltype(elementSolution(element, - std::declval(), - gridVariables.gridGeometry())); - return ElemSol{}; - } -} - // get the Dirichlet values for a boundary entity from a user problem -template +template decltype(auto) getDirichletValues(const Problem& problem, const Element& element, const BoundaryEntity& boundaryEntity, - const ElemSolState& elemSolState) + const ElemSol& elemSol) { - // given element solution state models an actual state - if constexpr (Dune::models()) + // given element solution models the experimental elemsol concept + if constexpr (Dune::models()) { if constexpr (ProblemTraits::hasTransientDirichletInterface) - return problem.dirichlet(element, boundaryEntity, elemSolState.timeLevel()); + return problem.dirichlet(element, boundaryEntity, elemSol.timeLevel()); else return problem.dirichlet(element, boundaryEntity); } diff --git a/dumux/assembly/fv/boxlocalassembler.hh b/dumux/assembly/fv/boxlocalassembler.hh index e54bcc4536..ea6e3178c2 100644 --- a/dumux/assembly/fv/boxlocalassembler.hh +++ b/dumux/assembly/fv/boxlocalassembler.hh @@ -380,9 +380,9 @@ protected: const auto& problem = curElemVolVars.gridVolVars().problem(); const auto origResiduals = evalLocalResidual(); - auto elemSolState = makeElementSolutionState(element, curVariables.gridVariables()); - auto& elemSol = elemSolState.elementSolution(); - const auto origElemSol = elemSol; + const auto solStateView = solutionStateView(curVariables.gridVariables()); + const auto origElemSol = elementSolution(element, solStateView, fvGeometry_.gridGeometry()); + auto elemSol = origElemSol; ////////////////////////////////////////////////////////////////////////////////////////////// // Calculate derivatives of the residual of all dofs in element with respect to themselves. // @@ -406,7 +406,7 @@ protected: { // update the volume variables and compute element residual elemSol[scvI.localDofIndex()][pvIdx] = priVar; - curVolVars.update(elemSolState, problem, element, scvI); + curVolVars.update(elemSol, problem, element, scvI); return evalLocalResidual(); }; diff --git a/dumux/assembly/fv/cclocalassembler.hh b/dumux/assembly/fv/cclocalassembler.hh index 75bd41587e..3703b16a0f 100644 --- a/dumux/assembly/fv/cclocalassembler.hh +++ b/dumux/assembly/fv/cclocalassembler.hh @@ -320,9 +320,9 @@ protected: const auto origVolVars = curVolVars; // element solution to be deflected - auto elemSolState = makeElementSolutionState(element, curVariables.gridVariables()); - auto& elemSol = elemSolState.elementSolution(); - const auto origElemSol = elemSolState.elementSolution(); + const auto solStateView = solutionStateView(curVariables.gridVariables()); + const auto origElemSol = elementSolution(element, solStateView, fvGeometry_.gridGeometry()); + auto elemSol = origElemSol; // derivatives in the neighbors with repect to the current elements // in index 0 we save the derivative of the element residual with respect to it's own dofs @@ -338,7 +338,7 @@ protected: // update the volume variables and the flux var cache elemSol[0][pvIdx] = priVar; - curVolVars.update(elemSolState, problem, element, scv); + curVolVars.update(elemSol, problem, element, scv); curElemFluxVarsCache.update(element, fvGeometry_, curElemVolVars); if (enableGridFluxVarsCache) gridVariables_.gridFluxVarsCache().updateElement(element, fvGeometry_, curElemVolVars); diff --git a/dumux/discretization/box/elementvolumevariables.hh b/dumux/discretization/box/elementvolumevariables.hh index 3ebbdb9096..e072f7ff66 100644 --- a/dumux/discretization/box/elementvolumevariables.hh +++ b/dumux/discretization/box/elementvolumevariables.hh @@ -27,7 +27,6 @@ #include #include -#include namespace Dumux { @@ -129,29 +128,12 @@ public: const FVElementGeometry& fvGeometry, const SolutionState& solState) { - // compatibility layer with new experimental assembly style - if constexpr (Dune::models()) - { - const auto elemSolState = makeElementSolutionState(element, - solState, - fvGeometry.gridGeometry()); - - // resize volume variables to the required size - volumeVariables_.resize(fvGeometry.numScv()); - for (auto&& scv : scvs(fvGeometry)) - volumeVariables_[scv.indexInElement()].update(elemSolState, gridVolVars().problem(), element, scv); - } - // current standard style (solution state = solution vector) - else - { - // get the solution at the dofs of the element - auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry()); - - // resize volume variables to the required size - volumeVariables_.resize(fvGeometry.numScv()); - for (auto&& scv : scvs(fvGeometry)) - volumeVariables_[scv.indexInElement()].update(elemSol, gridVolVars().problem(), element, scv); - } + const auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry()); + + // resize volume variables to the required size + volumeVariables_.resize(fvGeometry.numScv()); + for (auto&& scv : scvs(fvGeometry)) + volumeVariables_[scv.indexInElement()].update(elemSol, gridVolVars().problem(), element, scv); } const VolumeVariables& operator [](std::size_t scvIdx) const diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh index bfa6015505..1666053e04 100644 --- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh @@ -93,11 +93,7 @@ namespace CCMpfa { // index of the element the fvGeometry was bound to const auto boundElemIdx = fvGeometry.gridGeometry().elementMapper().index(element); - - // compatibility layer with new assembly: get element-local state - // This causes a slight overhead when using the old assembly (elem sol not needed) - using namespace Experimental::CompatibilityHelpers; - auto elemSolState = getElemSolState(element, solState, fvGeometry.gridGeometry()); + auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry()); // check each scvf in the index set for boundary presence for (auto scvfIdx : nodalIndexSet.gridScvfIndices()) @@ -116,15 +112,11 @@ namespace CCMpfa { // boundaries the "outside" vol vars cannot be properly defined. if (bcTypes.hasOnlyDirichlet()) { - // update the element solution state with the Dirichlet values - if constexpr (Dune::models()) - elemSolState.elementSolution()[0] = getDirichletValues(problem, insideElement, ivScvf, elemSolState); - // old assembly style (elemSolState = elementSolution) - else - elemSolState[0] = getDirichletValues(problem, insideElement, ivScvf, elemSolState); + using namespace Experimental::CompatibilityHelpers; + elemSol[0] = getDirichletValues(problem, insideElement, ivScvf, elemSol); VolumeVariables dirichletVolVars; - dirichletVolVars.update(elemSolState, problem, insideElement, fvGeometry.scv(insideScvIdx)); + dirichletVolVars.update(elemSol, problem, insideElement, fvGeometry.scv(insideScvIdx)); volVars.emplace_back(std::move(dirichletVolVars)); volVarIndices.push_back(ivScvf.outsideScvIdx()); @@ -153,11 +145,7 @@ namespace CCMpfa { const SolState& solState) { const auto& gridGeometry = fvGeometry.gridGeometry(); - - // compatibility layer with new assembly: get element-local state - // This causes a slight overhead when using old assembly (elem sol not needed) - using namespace Experimental::CompatibilityHelpers; - auto elemSolState = getElemSolState(element, solState, fvGeometry.gridGeometry()); + auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry()); // treat the BCs inside the element if (fvGeometry.hasBoundaryScvf()) @@ -174,15 +162,11 @@ namespace CCMpfa { // boundaries the "outside" vol vars cannot be properly defined. if (problem.boundaryTypes(element, scvf).hasOnlyDirichlet()) { - // update the element solution state with the Dirichlet values - if constexpr (Dune::models()) - elemSolState.elementSolution()[0] = getDirichletValues(problem, element, scvf, elemSolState); - // old assembly style (elemSolState = elementSolution) - else - elemSolState[0] = getDirichletValues(problem, element, scvf, elemSolState); + using namespace Experimental::CompatibilityHelpers; + elemSol[0] = getDirichletValues(problem, element, scvf, elemSol); VolumeVariables dirichletVolVars; - dirichletVolVars.update(elemSolState, problem, element, scvI); + dirichletVolVars.update(elemSol, problem, element, scvI); volVars.emplace_back(std::move(dirichletVolVars)); volVarIndices.push_back(scvf.outsideScvIdx()); @@ -343,13 +327,10 @@ public: volumeVariables_.reserve(numVolVars+maxNumBoundaryVolVars); volVarIndices_.reserve(numVolVars+maxNumBoundaryVolVars); - // compatibility layer with new assembly: get element-local state - using namespace Experimental::CompatibilityHelpers; - const auto elemSolState = getElemSolState(element, solState, gridGeometry); - VolumeVariables volVars; const auto& scvI = fvGeometry.scv(globalI); - volVars.update(elemSolState, problem, element, scvI); + const auto elemSol = elementSolution(element, solState, gridGeometry); + volVars.update(elemSol, problem, element, scvI); volVarIndices_.push_back(scvI.dofIndex()); volumeVariables_.emplace_back(std::move(volVars)); @@ -358,10 +339,10 @@ public: for (auto&& dataJ : assemblyMapI) { const auto& elementJ = gridGeometry.element(dataJ.globalJ); - const auto elemSolStateJ = getElemSolState(elementJ, solState, gridGeometry); + const auto elemSolJ = elementSolution(elementJ, solState, gridGeometry); const auto& scvJ = fvGeometry.scv(dataJ.globalJ); VolumeVariables volVarsJ; - volVarsJ.update(elemSolStateJ, problem, elementJ, scvJ); + volVarsJ.update(elemSolJ, problem, elementJ, scvJ); volVarIndices_.push_back(scvJ.dofIndex()); volumeVariables_.emplace_back(std::move(volVarsJ)); @@ -409,13 +390,9 @@ public: volumeVariables_.resize(1); volVarIndices_.resize(1); - // compatibility layer with new assembly: get element-local state - using namespace Experimental::CompatibilityHelpers; - const auto elemSolState = getElemSolState(element, solState, gridGeometry); - // update the volume variables of the element const auto& scv = fvGeometry.scv(eIdx); - volumeVariables_[0].update(elemSolState, + volumeVariables_[0].update(elementSolution(element, solState, gridGeometry), gridVolVars().problem(), element, scv); diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh index e64257ba6f..19da0f3774 100644 --- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh +++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh @@ -97,11 +97,7 @@ public: boundaryVolVarIndices_.reserve(fvGeometry.numScvf()); boundaryVolumeVariables_.reserve(fvGeometry.numScvf()); - // compatibility layer with new assembly: get element-local state - // This causes a slight overhead when using the old assembly (elem sol not needed) - using namespace Experimental::CompatibilityHelpers; - auto elemSolState = getElemSolState(element, solutionState, fvGeometry.gridGeometry()); - + auto elemSol = elementSolution(element, solutionState, fvGeometry.gridGeometry()); for (const auto& scvf : scvfs(fvGeometry)) { if (!scvf.boundary()) @@ -112,16 +108,11 @@ public: const auto bcTypes = problem.boundaryTypes(element, scvf); if (bcTypes.hasOnlyDirichlet()) { - // update the element solution state with the Dirichlet values - if constexpr (Dune::models()) - elemSolState.elementSolution()[0] = getDirichletValues(problem, element, scvf, elemSolState); - // old assembly style (elemSolState = elementSolution) - else - elemSolState[0] = getDirichletValues(problem, element, scvf, elemSolState); - + using namespace Experimental::CompatibilityHelpers; + elemSol[0] = getDirichletValues(problem, element, scvf, elemSol); const auto& scvI = fvGeometry.scv(scvf.insideScvIdx()); VolumeVariables volVars; - volVars.update(elemSolState, problem, element, scvI); + volVars.update(elemSol, problem, element, scvI); boundaryVolumeVariables_.emplace_back(std::move(volVars)); boundaryVolVarIndices_.push_back(scvf.outsideScvIdx()); @@ -201,10 +192,9 @@ public: int localIdx = 0; // update the volume variables of the element at hand - using namespace Experimental::CompatibilityHelpers; - auto elemSolState = getElemSolState(element, solState, gridGeometry); + auto elemSol = elementSolution(element, solState, gridGeometry); const auto& scvI = fvGeometry.scv(globalI); - volumeVariables_[localIdx].update(elemSolState, + volumeVariables_[localIdx].update(elemSol, problem, element, scvI); @@ -216,7 +206,7 @@ public: { const auto& elementJ = gridGeometry.element(dataJ.globalJ); auto&& scvJ = fvGeometry.scv(dataJ.globalJ); - volumeVariables_[localIdx].update(getElemSolState(elementJ, solState, gridGeometry), + volumeVariables_[localIdx].update(elementSolution(elementJ, solState, gridGeometry), problem, elementJ, scvJ); @@ -237,16 +227,12 @@ public: const auto bcTypes = problem.boundaryTypes(element, scvf); if (bcTypes.hasOnlyDirichlet()) { - // update the element solution state with the Dirichlet values - if constexpr (Dune::models()) - elemSolState.elementSolution()[0] = getDirichletValues(problem, element, scvf, elemSolState); - // old assembly style (elemSolState = elementSolution) - else - elemSolState[0] = getDirichletValues(problem, element, scvf, elemSolState); + using namespace Experimental::CompatibilityHelpers; + elemSol[0] = getDirichletValues(problem, element, scvf, elemSol); volumeVariables_.resize(localIdx+1); volVarIndices_.resize(localIdx+1); - volumeVariables_[localIdx].update(elemSolState, problem, element, scvI); + volumeVariables_[localIdx].update(elemSol, problem, element, scvI); volVarIndices_[localIdx] = scvf.outsideScvIdx(); ++localIdx; } @@ -287,10 +273,9 @@ public: volVarIndices_.resize(1); // update the volume variables of the element - using namespace Experimental::CompatibilityHelpers; - const auto elemSolState = getElemSolState(element, solState, fvGeometry.gridGeometry()); + const auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry()); const auto& scv = fvGeometry.scv(eIdx); - volumeVariables_[0].update(elemSolState, + volumeVariables_[0].update(elemSol, gridVolVars().problem(), element, scv); diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh index e4b5a2d974..7f1ac235e1 100644 --- a/dumux/discretization/fvgridvariables.hh +++ b/dumux/discretization/fvgridvariables.hh @@ -222,8 +222,7 @@ public: void bind(const Element& element, const FVElementGeometry& fvGeometry) { - const Experimental::SolutionState solState(&gridVariables().dofs(), - &gridVariables().timeLevel()); + const auto solState = solutionStateView(gridVariables()); elemVolVars_.bind(element, fvGeometry, solState); elemFluxVarsCache_.bind(element, fvGeometry, elemVolVars_); } @@ -236,9 +235,7 @@ public: void bindElemVolVars(const Element& element, const FVElementGeometry& fvGeometry) { - Dumux::Experimental::SolutionState solState(&gridVariables().dofs(), - &gridVariables().timeLevel()); - + const auto solState = solutionStateView(gridVariables()); if (bindEntireStencil) elemVolVars_.bind(element, fvGeometry, solState); else diff --git a/dumux/discretization/solutionstate.hh b/dumux/discretization/solutionstate.hh index 55aa657f88..e6ff94c689 100644 --- a/dumux/discretization/solutionstate.hh +++ b/dumux/discretization/solutionstate.hh @@ -26,6 +26,9 @@ #ifndef DUMUX_SOLUTION_STATE_HH #define DUMUX_SOLUTION_STATE_HH +#include +#include + #include #include @@ -43,12 +46,13 @@ struct SolutionState ); }; -//! Concept for element solution states -struct ElementSolutionState +//! Concept for element solutions +struct ElementSolution { template auto require(const T& t) -> decltype( - t.elementSolution(), + t.size(), + t[0], t.timeLevel() ); }; @@ -57,103 +61,115 @@ struct ElementSolutionState /*! * \ingroup Discretization - * \brief State class to represent the state of the primary variables of a - * numerical model, consisting of the spatial degrees of freedom and - * potentially a corresponding time level. + * \brief View on the state of a numerical solution, consisting of the spatial + * degrees of freedom and a corresponding time level. */ template -class SolutionState +class SolutionStateView { public: using SolutionVector = SV; using TimeLevel = TL; //! Constructor - SolutionState(const SolutionVector* dofs, - const TimeLevel* timeLevel) + template + SolutionStateView(SolutionVector&& dofs, + TimeLevel&& timeLevel) : dofs_(dofs) , timeLevel_(timeLevel) - {} + { + static_assert(std::is_lvalue_reference_v, + "Solution vector must be an lvalue reference"); + static_assert(std::is_lvalue_reference_v, + "Time level must be an lvalue reference"); + } + + //! Constructor from a solution state + template + SolutionStateView(const SolState& other) + : SolutionStateView(other.dofs(), other.timeLevel()) + { + static_assert(Dune::models(), + "Given type does not model SolutionState concept"); + } //! return the solution vector const SolutionVector& dofs() const - { return *dofs_; } + { return dofs_; } //! return the time level const TimeLevel& timeLevel() const - { return *timeLevel_; } + { return timeLevel_; } private: - const SolutionVector* dofs_; - const TimeLevel* timeLevel_; + const SolutionVector& dofs_; + const TimeLevel& timeLevel_; }; +/*! + * \ingroup Discretization + * \brief Make a solution state view from a solution state + */ +template +auto solutionStateView(const SolutionState& solState) +{ + static_assert(Dune::models(), + "Given type does not model SolutionState concept"); + using SV = std::decay_t; + using TL = std::decay_t; + return SolutionStateView(solState); +} + /*! * \ingroup Discretization * \brief State class to represent the element-local state of the primary * variables of a numerical model, consisting of the spatial degrees * of freedom and potentially a corresponding time level. + * \note Here, we extend the concept of elementSolution to additionally carry + * time information. Therefore, for now we inherit from the given element + * solution, extending its interface */ -template -class ElementSolutionState +template +class ElementSolution : public BaseElemSol { public: - using ElementSolution = ES; using TimeLevel = TL; //! Constructor - ElementSolutionState(ElementSolution&& es, - const TimeLevel* timeLevel) - : elemSol_(std::move(es)) + template, TL>, int> = 0> + ElementSolution(BaseElemSol&& base, + TheTimeLevel&& timeLevel) + : BaseElemSol(std::forward(base)) , timeLevel_(timeLevel) - {} - - //! return the element solution vector - const ElementSolution& elementSolution() const - { return elemSol_; } - - //! return the element solution vector - ElementSolution& elementSolution() - { return elemSol_; } + { + static_assert(std::is_lvalue_reference_v, + "Time level must be an lvalue reference"); + } //! return the time level const TimeLevel& timeLevel() const - { return *timeLevel_; } + { return timeLevel_; } private: - ElementSolution elemSol_; - const TimeLevel* timeLevel_; + const TimeLevel& timeLevel_; }; /*! * \ingroup Discretization - * \brief Convenience factory function to make an element solution state - * from a global solution state. + * \brief Make an element solution from a global solution state. */ -template -auto makeElementSolutionState(const Element& element, - const SolutionState& state, - const GridGeometry& gridGeom) +template(), int> = 0> +auto elementSolution(const Element& element, + const SolState& solState, + const GridGeometry& gridGeometry) { - static_assert(Dune::models(), - "Provided type does not fulfill solution state interface"); - - auto elemSol = elementSolution(element, state.dofs(), gridGeom); - return ElementSolutionState{std::move(elemSol), &state.timeLevel()}; -} + auto elemSol = elementSolution(element, solState.dofs(), gridGeometry); -/*! - * \ingroup Discretization - * \brief Convenience factory function to make an element solution state - * from grid variables. - */ -template -auto makeElementSolutionState(const Element& element, - const GridVariables& gridVariables) -{ - const auto& gridGeom = gridVariables.gridGeometry(); - auto elemSol = elementSolution(element, gridVariables.dofs(), gridGeom); - return ElementSolutionState{std::move(elemSol), &gridVariables.timeLevel()}; + using Base = std::decay_t; + using TL = std::decay_t; + return ElementSolution{std::move(elemSol), solState.timeLevel()}; } } // end namespace Dumux::Experimental diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh index 530b5ed243..be46ab7168 100644 --- a/dumux/io/vtkoutputmodule.hh +++ b/dumux/io/vtkoutputmodule.hh @@ -469,8 +469,7 @@ private: fvGeometry.bind(element); if constexpr (Experimental::areExperimentalGridVars) { - const Experimental::SolutionState solState(&gridVariables_.dofs(), - &gridVariables_.timeLevel()); + const auto solState = solutionStateView(gridVariables_); elemVolVars.bind(element, fvGeometry, solState); } else @@ -481,8 +480,7 @@ private: fvGeometry.bindElement(element); if constexpr (Experimental::areExperimentalGridVars) { - const Experimental::SolutionState solState(&gridVariables_.dofs(), - &gridVariables_.timeLevel()); + const auto solState = solutionStateView(gridVariables_); elemVolVars.bindElement(element, fvGeometry, solState); } else @@ -676,8 +674,7 @@ private: fvGeometry.bind(element); if constexpr (Experimental::areExperimentalGridVars) { - const Experimental::SolutionState solState(&gridVariables_.dofs(), - &gridVariables_.timeLevel()); + const auto solState = solutionStateView(gridVariables_); elemVolVars.bind(element, fvGeometry, solState); } else @@ -688,8 +685,7 @@ private: fvGeometry.bindElement(element); if constexpr (Experimental::areExperimentalGridVars) { - const Experimental::SolutionState solState(&gridVariables_.dofs(), - &gridVariables_.timeLevel()); + const auto solState = solutionStateView(gridVariables_); elemVolVars.bindElement(element, fvGeometry, solState); } else diff --git a/dumux/material/solidstates/updatesolidvolumefractions.hh b/dumux/material/solidstates/updatesolidvolumefractions.hh index 2f327be0b9..dd9082e65a 100644 --- a/dumux/material/solidstates/updatesolidvolumefractions.hh +++ b/dumux/material/solidstates/updatesolidvolumefractions.hh @@ -24,8 +24,6 @@ #ifndef DUMUX_UPDATE_SOLID_VOLUME_FRACTION_HH #define DUMUX_UPDATE_SOLID_VOLUME_FRACTION_HH -#include - namespace Dumux { /*! @@ -52,15 +50,7 @@ void updateSolidVolumeFractions(const ElemSol& elemSol, if (!(solidState.isInert())) { - // compatibility with new experimental assembly style - const auto& priVars = [&] () - { - if constexpr (Dune::models()) - return elemSol.elementSolution()[scv.localDofIndex()]; - else - return elemSol[scv.localDofIndex()]; - } (); - + const auto& priVars = elemSol[scv.localDofIndex()]; for (int sCompIdx = 0; sCompIdx < solidState.numComponents-solidState.numInertComponents; ++sCompIdx) solidState.setVolumeFraction(sCompIdx, priVars[solidVolFracOffset + sCompIdx]); } diff --git a/dumux/porousmediumflow/1p/volumevariables.hh b/dumux/porousmediumflow/1p/volumevariables.hh index 03469aa48a..dc0a76dce6 100644 --- a/dumux/porousmediumflow/1p/volumevariables.hh +++ b/dumux/porousmediumflow/1p/volumevariables.hh @@ -25,8 +25,6 @@ #ifndef DUMUX_1P_VOLUME_VARIABLES_HH #define DUMUX_1P_VOLUME_VARIABLES_HH -#include - #include #include @@ -113,18 +111,10 @@ public: FluidState& fluidState, SolidState& solidState) { - // compatibility with new experimental assembly style - const auto& priVars = [&elemSol, &scv] () - { - if constexpr (Dune::models()) - return elemSol.elementSolution()[scv.localDofIndex()]; - else - return elemSol[scv.localDofIndex()]; - } (); - EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState); fluidState.setSaturation(/*phaseIdx=*/0, 1.); + const auto& priVars = elemSol[scv.localDofIndex()]; fluidState.setPressure(/*phaseIdx=*/0, priVars[Indices::pressureIdx]); // saturation in a single phase is always 1 and thus redundant diff --git a/dumux/porousmediumflow/nonisothermal/volumevariables.hh b/dumux/porousmediumflow/nonisothermal/volumevariables.hh index ebeaedddd3..0c45856d29 100644 --- a/dumux/porousmediumflow/nonisothermal/volumevariables.hh +++ b/dumux/porousmediumflow/nonisothermal/volumevariables.hh @@ -176,14 +176,7 @@ public: FluidState& fluidState, SolidState& solidState) { - // compatibility with new experimental assembly style - const auto& priVars = [&elemSol, &scv] () - { - if constexpr (Dune::models()) - return elemSol.elementSolution()[scv.localDofIndex()]; - else - return elemSol[scv.localDofIndex()]; - } (); + const auto& priVars = elemSol[scv.localDofIndex()]; if constexpr (fullThermalEquilibrium) { diff --git a/dumux/porousmediumflow/volumevariables.hh b/dumux/porousmediumflow/volumevariables.hh index 51452884b5..bfca5829db 100644 --- a/dumux/porousmediumflow/volumevariables.hh +++ b/dumux/porousmediumflow/volumevariables.hh @@ -26,8 +26,6 @@ #ifndef DUMUX_POROUSMEDIUMFLOW_VOLUME_VARIABLES_HH #define DUMUX_POROUSMEDIUMFLOW_VOLUME_VARIABLES_HH -#include - namespace Dumux { /*! @@ -68,19 +66,8 @@ public: const Element& element, const Scv& scv) { - // compatibility layer with new experimental assembly style (elem sol = element solution state) - if constexpr (Dune::models()) - { - const auto& elemSolState = elemSol; - priVars_ = elemSolState.elementSolution()[scv.localDofIndex()]; - extrusionFactor_ = problem.spatialParams().extrusionFactor(element, scv, elemSolState); - } - // current standard style (element solution-based) - else - { - priVars_ = elemSol[scv.localDofIndex()]; - extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol); - } + priVars_ = elemSol[scv.localDofIndex()]; + extrusionFactor_ = problem.extrusionFactor(element, scv, elemSol); } /*! -- GitLab From a08e0d54d5ccee1fef0dc5a8dae6d454244c0b52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Thu, 11 Mar 2021 13:36:16 +0100 Subject: [PATCH 30/32] [experimental][fvgridvars] introduce bind policy --- dumux/discretization/fvgridvariables.hh | 46 +++++++++++++++---------- 1 file changed, 28 insertions(+), 18 deletions(-) diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh index 7f1ac235e1..c0f193b68e 100644 --- a/dumux/discretization/fvgridvariables.hh +++ b/dumux/discretization/fvgridvariables.hh @@ -183,6 +183,17 @@ private: namespace Dumux::Experimental { +/*! + * \ingroup Discretization + * \brief Policy for binding local views of finite volume grid variables. + */ +enum class FVGridVariablesBindPolicy +{ + all, + volVarsOnly, + volVarsOfElementOnly +}; + /*! * \ingroup Discretization * \brief Finite volume-specific local view on grid variables. @@ -220,31 +231,30 @@ public: * \param fvGeometry Local view on the grid geometry */ void bind(const Element& element, - const FVElementGeometry& fvGeometry) + const FVElementGeometry& fvGeometry, + FVGridVariablesBindPolicy policy = FVGridVariablesBindPolicy::all) { const auto solState = solutionStateView(gridVariables()); - elemVolVars_.bind(element, fvGeometry, solState); - elemFluxVarsCache_.bind(element, fvGeometry, elemVolVars_); - } - /*! - * \brief Bind only the volume variables local view to a grid element. - * \param element The grid element - * \param fvGeometry Local view on the grid geometry - */ - void bindElemVolVars(const Element& element, - const FVElementGeometry& fvGeometry) - { - const auto solState = solutionStateView(gridVariables()); - if (bindEntireStencil) + if (policy == FVGridVariablesBindPolicy::all) + { elemVolVars_.bind(element, fvGeometry, solState); + elemFluxVarsCache_.bind(element, fvGeometry, elemVolVars_); + } else - elemVolVars_.bindElement(element, fvGeometry, solState); - - // unbind flux variables cache - elemFluxVarsCache_ = localView(gridVariables().gridFluxVarsCache()); + { + // unbind flux variables cache + elemFluxVarsCache_ = localView(gridVariables().gridFluxVarsCache()); + if (policy == FVGridVariablesBindPolicy::volVarsOnly) + elemVolVars_.bind(element, fvGeometry, solState); + else if (policy == FVGridVariablesBindPolicy::volVarsOfElementOnly) + elemVolVars_.bindElement(element, fvGeometry, solState); + else + DUNE_THROW(Dune::NotImplemented, "Bind for given policy"); + } } + //! return reference to the elem vol vars const ElementVolumeVariables& elemVolVars() const { return elemVolVars_; } ElementVolumeVariables& elemVolVars() { return elemVolVars_; } -- GitLab From 963ef00ceac5c56dff320fdc1b92e6a2058478b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 15 Mar 2021 17:53:36 +0100 Subject: [PATCH 31/32] [disc] refactor solution state concept --- dumux/assembly/experimentalhelpers.hh | 4 +- dumux/assembly/fv/boxlocalassembler.hh | 5 +- dumux/assembly/fv/cclocalassembler.hh | 7 +- dumux/discretization/CMakeLists.txt | 3 +- dumux/discretization/concepts.hh | 62 +++++++++++ dumux/discretization/elementsolution.hh | 63 +++++++++++ dumux/discretization/fvgridvariables.hh | 4 +- ...{solutionstate.hh => solutionstateview.hh} | 102 ++---------------- dumux/io/vtkoutputmodule.hh | 10 +- 9 files changed, 152 insertions(+), 108 deletions(-) create mode 100644 dumux/discretization/concepts.hh rename dumux/discretization/{solutionstate.hh => solutionstateview.hh} (52%) diff --git a/dumux/assembly/experimentalhelpers.hh b/dumux/assembly/experimentalhelpers.hh index ccec82e842..720976d836 100644 --- a/dumux/assembly/experimentalhelpers.hh +++ b/dumux/assembly/experimentalhelpers.hh @@ -32,9 +32,7 @@ #include #include -#include -#include -#include +#include namespace Dumux::Experimental { namespace CompatibilityHelpers { diff --git a/dumux/assembly/fv/boxlocalassembler.hh b/dumux/assembly/fv/boxlocalassembler.hh index ea6e3178c2..633323c546 100644 --- a/dumux/assembly/fv/boxlocalassembler.hh +++ b/dumux/assembly/fv/boxlocalassembler.hh @@ -37,7 +37,8 @@ #include #include -#include +#include +#include #include @@ -380,7 +381,7 @@ protected: const auto& problem = curElemVolVars.gridVolVars().problem(); const auto origResiduals = evalLocalResidual(); - const auto solStateView = solutionStateView(curVariables.gridVariables()); + const Experimental::SolutionStateView solStateView(curVariables.gridVariables()); const auto origElemSol = elementSolution(element, solStateView, fvGeometry_.gridGeometry()); auto elemSol = origElemSol; diff --git a/dumux/assembly/fv/cclocalassembler.hh b/dumux/assembly/fv/cclocalassembler.hh index 3703b16a0f..396b732f7c 100644 --- a/dumux/assembly/fv/cclocalassembler.hh +++ b/dumux/assembly/fv/cclocalassembler.hh @@ -36,8 +36,11 @@ #include #include + #include -#include +#include +#include + #include @@ -320,7 +323,7 @@ protected: const auto origVolVars = curVolVars; // element solution to be deflected - const auto solStateView = solutionStateView(curVariables.gridVariables()); + const Experimental::SolutionStateView solStateView(curVariables.gridVariables()); const auto origElemSol = elementSolution(element, solStateView, fvGeometry_.gridGeometry()); auto elemSol = origElemSol; diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt index 9519315c93..82f9d4dcf9 100644 --- a/dumux/discretization/CMakeLists.txt +++ b/dumux/discretization/CMakeLists.txt @@ -11,6 +11,7 @@ box.hh ccmpfa.hh cctpfa.hh checkoverlapsize.hh +concepts.hh elementsolution.hh evalgradients.hh evalsolution.hh @@ -28,7 +29,7 @@ rotationsymmetricgridgeometrytraits.hh rotationsymmetricscv.hh rotationsymmetricscvf.hh scvandscvfiterators.hh -solutionstate.hh +solutionstateview.hh staggered.hh subcontrolvolumebase.hh subcontrolvolumefacebase.hh diff --git a/dumux/discretization/concepts.hh b/dumux/discretization/concepts.hh new file mode 100644 index 0000000000..65bdff436e --- /dev/null +++ b/dumux/discretization/concepts.hh @@ -0,0 +1,62 @@ +// -*- 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 3 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 . * + *****************************************************************************/ +/*! + * \file + * \ingroup Discretization + * \brief Concepts related to discretization methods. + */ +#ifndef DUMUX_DISCRETIZATION_CONCEPT_HH +#define DUMUX_DISCRETIZATION_CONCEPT_HH + +#include + +namespace Dumux{ + +// experimental concepts +namespace Experimental::Concept { + +//! Concept for states of discrete numerical solutions +struct SolutionState +{ + template + auto require(const T& t) -> decltype( + Dune::Concept::requireType(), + Dune::Concept::requireType(), + Dune::Concept::requireConvertible(t.dofs()), + Dune::Concept::requireConvertible(t.timeLevel()) + ); +}; + +//! Concept for element-local discrete numerical solutions +struct ElementSolution +{ + template + auto require(const T& t) -> decltype( + Dune::Concept::requireType(), + Dune::Concept::requireType(), + Dune::Concept::requireConvertible(t.size()), + Dune::Concept::requireConvertible(t[0]), + Dune::Concept::requireConvertible(t.timeLevel()) + ); +}; + +} // end namespace Experimental::Concept +} // end namespace Dumux + +#endif diff --git a/dumux/discretization/elementsolution.hh b/dumux/discretization/elementsolution.hh index 3f9d62f11b..3e5e8bae55 100644 --- a/dumux/discretization/elementsolution.hh +++ b/dumux/discretization/elementsolution.hh @@ -24,6 +24,11 @@ #ifndef DUMUX_DISCRETIZATION_ELEMENT_SOLUTION_HH #define DUMUX_DISCRETIZATION_ELEMENT_SOLUTION_HH +#include + +#include +#include + #include #include #include @@ -32,6 +37,64 @@ namespace Dumux { struct EmptyElementSolution {}; +namespace Experimental { + +/*! + * \ingroup Discretization + * \brief State class to represent the element-local state of the primary + * variables of a numerical model, consisting of the spatial degrees + * of freedom and potentially a corresponding time level. + * \note Here, we extend the concept of elementSolution to additionally carry + * time information. Therefore, for now we inherit from the given element + * solution, extending its interface + */ +template +class ElementSolution : public BaseElemSol +{ +public: + using TimeLevel = TL; + + //! Constructor + template, TL>, int> = 0> + ElementSolution(BaseElemSol&& base, + TheTimeLevel&& timeLevel) + : BaseElemSol(std::forward(base)) + , timeLevel_(timeLevel) + { + static_assert(std::is_lvalue_reference_v, + "Time level must be an lvalue reference"); + } + + //! return the time level + const TimeLevel& timeLevel() const + { return timeLevel_; } + +private: + const TimeLevel& timeLevel_; +}; + +// hide deduction guides from doxygen +#ifndef DOXYGEN +template ElementSolution(const B& b, const TL& t) -> ElementSolution; +template ElementSolution(B&& b, const TL& t) -> ElementSolution; +#endif // DOXYGEN + +/*! + * \ingroup Discretization + * \brief Make an element solution from a global solution state. + */ +template(), int> = 0> +auto elementSolution(const Element& element, + const SolState& solState, + const GridGeometry& gridGeometry) +{ + auto elemSol = elementSolution(element, solState.dofs(), gridGeometry); + return ElementSolution{std::move(elemSol), solState.timeLevel()}; +} + +} // end namespace Experimental } // end namespace Dumux #endif diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh index c0f193b68e..2f43626a5a 100644 --- a/dumux/discretization/fvgridvariables.hh +++ b/dumux/discretization/fvgridvariables.hh @@ -179,7 +179,7 @@ private: #include #include -#include +#include namespace Dumux::Experimental { @@ -234,7 +234,7 @@ public: const FVElementGeometry& fvGeometry, FVGridVariablesBindPolicy policy = FVGridVariablesBindPolicy::all) { - const auto solState = solutionStateView(gridVariables()); + const Experimental::SolutionStateView solState(gridVariables()); if (policy == FVGridVariablesBindPolicy::all) { diff --git a/dumux/discretization/solutionstate.hh b/dumux/discretization/solutionstateview.hh similarity index 52% rename from dumux/discretization/solutionstate.hh rename to dumux/discretization/solutionstateview.hh index e6ff94c689..cd619366bc 100644 --- a/dumux/discretization/solutionstate.hh +++ b/dumux/discretization/solutionstateview.hh @@ -23,42 +23,17 @@ * primary variables of a numerical model, consisting of the spatial degrees * of freedom and potentially a corresponding time level. */ -#ifndef DUMUX_SOLUTION_STATE_HH -#define DUMUX_SOLUTION_STATE_HH +#ifndef DUMUX_DISCRETIZATION_SOLUTION_STATE_HH +#define DUMUX_DISCRETIZATION_SOLUTION_STATE_HH #include #include #include -#include +#include namespace Dumux::Experimental { -namespace Concept { - -//! Concept for solution states -struct SolutionState -{ - template - auto require(const T& t) -> decltype( - t.dofs(), - t.timeLevel() - ); -}; - -//! Concept for element solutions -struct ElementSolution -{ - template - auto require(const T& t) -> decltype( - t.size(), - t[0], - t.timeLevel() - ); -}; - -} // end namespace Concept - /*! * \ingroup Discretization * \brief View on the state of a numerical solution, consisting of the spatial @@ -106,71 +81,12 @@ private: const TimeLevel& timeLevel_; }; -/*! - * \ingroup Discretization - * \brief Make a solution state view from a solution state - */ -template -auto solutionStateView(const SolutionState& solState) -{ - static_assert(Dune::models(), - "Given type does not model SolutionState concept"); - using SV = std::decay_t; - using TL = std::decay_t; - return SolutionStateView(solState); -} - -/*! - * \ingroup Discretization - * \brief State class to represent the element-local state of the primary - * variables of a numerical model, consisting of the spatial degrees - * of freedom and potentially a corresponding time level. - * \note Here, we extend the concept of elementSolution to additionally carry - * time information. Therefore, for now we inherit from the given element - * solution, extending its interface - */ -template -class ElementSolution : public BaseElemSol -{ -public: - using TimeLevel = TL; - - //! Constructor - template, TL>, int> = 0> - ElementSolution(BaseElemSol&& base, - TheTimeLevel&& timeLevel) - : BaseElemSol(std::forward(base)) - , timeLevel_(timeLevel) - { - static_assert(std::is_lvalue_reference_v, - "Time level must be an lvalue reference"); - } - - //! return the time level - const TimeLevel& timeLevel() const - { return timeLevel_; } - -private: - const TimeLevel& timeLevel_; -}; - -/*! - * \ingroup Discretization - * \brief Make an element solution from a global solution state. - */ -template(), int> = 0> -auto elementSolution(const Element& element, - const SolState& solState, - const GridGeometry& gridGeometry) -{ - auto elemSol = elementSolution(element, solState.dofs(), gridGeometry); - - using Base = std::decay_t; - using TL = std::decay_t; - return ElementSolution{std::move(elemSol), solState.timeLevel()}; -} +// hide deduction guides from doxygen +#ifndef DOXYGEN +template SolutionStateView(const SolState& b) + -> SolutionStateView; +#endif // DOXYGEN } // end namespace Dumux::Experimental diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh index be46ab7168..c1cb695f37 100644 --- a/dumux/io/vtkoutputmodule.hh +++ b/dumux/io/vtkoutputmodule.hh @@ -46,7 +46,7 @@ #include #include -#include +#include #include "vtkfunction.hh" #include "velocityoutput.hh" @@ -469,7 +469,7 @@ private: fvGeometry.bind(element); if constexpr (Experimental::areExperimentalGridVars) { - const auto solState = solutionStateView(gridVariables_); + const Experimental::SolutionStateView solState(gridVariables_); elemVolVars.bind(element, fvGeometry, solState); } else @@ -480,7 +480,7 @@ private: fvGeometry.bindElement(element); if constexpr (Experimental::areExperimentalGridVars) { - const auto solState = solutionStateView(gridVariables_); + const Experimental::SolutionStateView solState(gridVariables_); elemVolVars.bindElement(element, fvGeometry, solState); } else @@ -674,7 +674,7 @@ private: fvGeometry.bind(element); if constexpr (Experimental::areExperimentalGridVars) { - const auto solState = solutionStateView(gridVariables_); + const Experimental::SolutionStateView solState(gridVariables_); elemVolVars.bind(element, fvGeometry, solState); } else @@ -685,7 +685,7 @@ private: fvGeometry.bindElement(element); if constexpr (Experimental::areExperimentalGridVars) { - const auto solState = solutionStateView(gridVariables_); + const Experimental::SolutionStateView solState(gridVariables_); elemVolVars.bindElement(element, fvGeometry, solState); } else -- GitLab From 1db19b5c1afb5729c4b076899172a9ebc793e3f7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dennis=20Gl=C3=A4ser?= Date: Mon, 15 Mar 2021 18:02:16 +0100 Subject: [PATCH 32/32] [experimental][fvprob] add point source interface --- dumux/assembly/fv/operators.hh | 3 +-- dumux/common/fvproblem.hh | 23 ++++++++++++++++++++++- 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/dumux/assembly/fv/operators.hh b/dumux/assembly/fv/operators.hh index b25a46496c..0fbc70812b 100644 --- a/dumux/assembly/fv/operators.hh +++ b/dumux/assembly/fv/operators.hh @@ -140,8 +140,7 @@ public: source += problem.source(context, scv); // add contribution from possible point sources - // TODO: Point sources - // source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv); + source += problem.scvPointSources(context, scv); // multiply with scv volume const auto& elemVolVars = context.elementVariables().elemVolVars(); diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh index e4cc354f3a..ed4d05e119 100644 --- a/dumux/common/fvproblem.hh +++ b/dumux/common/fvproblem.hh @@ -726,7 +726,28 @@ public: const SubControlVolume& scv) const { return this->asImp_().sourceAtPos(scv.center()); } - //! TODO: Point sources! + /*! + * \brief Adds contribution of point sources for a specific sub control volume + * to the values. + * + * \param context The element-local context + * \param scv The sub control volume + * \note The element-local context consists of an element and local + * geometric information, together with the primary/secondary + * variables in that local scope. Potentially, users can define + * the context to carry an additional object containing further + * locally required data. + * + */ + template + NumEqVector scvPointSources(const LocalContext& context, + const SubControlVolume& scv) const + { + const auto& elemVolVars = context.elementVariables().elemVolVars(); + const auto& fvGeometry = context.elementGridGeometry(); + const auto& element = fvGeometry.element(); + return ParentType::scvPointSources(element, fvGeometry, elemVolVars, scv); + } }; } // end namespace Experimental -- GitLab