From 233f3355afea832f19740f6b3248ced4c35fc8c7 Mon Sep 17 00:00:00 2001
From: Timo Koch <timokoch@math.uio.no>
Date: Sat, 18 Mar 2023 16:57:15 +0100
Subject: [PATCH] [assembler][md] Use unified subdomain CVFE assembler
---
dumux/multidomain/fvassembler.hh | 22 +-
...bler.hh => subdomaincvfelocalassembler.hh} | 52 +--
.../subdomainfcdiamondlocalassembler.hh | 418 ------------------
.../subdomainpq1bubblelocalassembler.hh | 418 ------------------
4 files changed, 30 insertions(+), 880 deletions(-)
rename dumux/multidomain/{subdomainboxlocalassembler.hh => subdomaincvfelocalassembler.hh} (91%)
delete mode 100644 dumux/multidomain/subdomainfcdiamondlocalassembler.hh
delete mode 100644 dumux/multidomain/subdomainpq1bubblelocalassembler.hh
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index d2cac7ef7b..dee7c820ee 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -45,11 +45,9 @@
#include "couplingjacobianpattern.hh"
#include "subdomaincclocalassembler.hh"
-#include "subdomainboxlocalassembler.hh"
+#include "subdomaincvfelocalassembler.hh"
#include "subdomainstaggeredlocalassembler.hh"
#include "subdomainfclocalassembler.hh"
-#include "subdomainfcdiamondlocalassembler.hh"
-#include "subdomainpq1bubblelocalassembler.hh"
#include <dumux/discretization/method.hh>
@@ -153,10 +151,10 @@ private:
using type = SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
};
- template<std::size_t id>
- struct SubDomainAssemblerType<DiscretizationMethods::Box, id>
+ template<std::size_t id, class DM>
+ struct SubDomainAssemblerType<DiscretizationMethods::CVFE<DM>, id>
{
- using type = SubDomainBoxLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
+ using type = SubDomainCVFELocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
};
template<std::size_t id>
@@ -171,18 +169,6 @@ private:
using type = SubDomainFaceCenteredLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
};
- template<std::size_t id>
- struct SubDomainAssemblerType<DiscretizationMethods::FCDiamond, id>
- {
- using type = SubDomainFaceCenteredDiamondLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
- };
-
- template<std::size_t id>
- struct SubDomainAssemblerType<DiscretizationMethods::PQ1Bubble, id>
- {
- using type = SubDomainPQ1BubbleLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
- };
-
template<std::size_t id>
using SubDomainAssembler = typename SubDomainAssemblerType<typename GridGeometry<id>::DiscretizationMethod, id>::type;
diff --git a/dumux/multidomain/subdomainboxlocalassembler.hh b/dumux/multidomain/subdomaincvfelocalassembler.hh
similarity index 91%
rename from dumux/multidomain/subdomainboxlocalassembler.hh
rename to dumux/multidomain/subdomaincvfelocalassembler.hh
index bd33d1c237..4ae95d50c6 100644
--- a/dumux/multidomain/subdomainboxlocalassembler.hh
+++ b/dumux/multidomain/subdomaincvfelocalassembler.hh
@@ -19,12 +19,12 @@
/*!
* \file
* \ingroup Assembly
- * \ingroup BoxDiscretization
+ * \ingroup CVFEDiscretization
* \ingroup MultiDomain
- * \brief An assembler for Jacobian and residual contribution per element (box methods) for multidomain problems
+ * \brief An assembler for Jacobian and residual contribution per element (CVFE methods) for multidomain problems
*/
-#ifndef DUMUX_MULTIDOMAIN_BOX_LOCAL_ASSEMBLER_HH
-#define DUMUX_MULTIDOMAIN_BOX_LOCAL_ASSEMBLER_HH
+#ifndef DUMUX_MULTIDOMAIN_SUBDOMAIN_CVFE_LOCAL_ASSEMBLER_HH
+#define DUMUX_MULTIDOMAIN_SUBDOMAIN_CVFE_LOCAL_ASSEMBLER_HH
#include <dune/common/reservedvector.hh>
#include <dune/common/indices.hh>
@@ -46,9 +46,9 @@ namespace Dumux {
/*!
* \ingroup Assembly
- * \ingroup BoxDiscretization
+ * \ingroup CVFEDiscretization
* \ingroup MultiDomain
- * \brief A base class for all box local assemblers
+ * \brief A base class for all CVFE subdomain local assemblers
* \tparam id the id of the sub domain
* \tparam TypeTag the TypeTag
* \tparam Assembler the assembler type
@@ -56,7 +56,7 @@ namespace Dumux {
* \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
*/
template<std::size_t id, class TypeTag, class Assembler, class Implementation, DiffMethod dm, bool implicit>
-class SubDomainBoxLocalAssemblerBase : public CVFELocalAssembler<TypeTag, Assembler, dm, implicit, Implementation>
+class SubDomainCVFELocalAssemblerBase : public CVFELocalAssembler<TypeTag, Assembler, dm, implicit, Implementation>
{
using ParentType = CVFELocalAssembler<TypeTag, Assembler, dm, implicit, Implementation>;
@@ -93,7 +93,7 @@ public:
using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
// the constructor
- explicit SubDomainBoxLocalAssemblerBase(
+ explicit SubDomainCVFELocalAssemblerBase(
const Assembler& assembler,
const Element& element,
const SolutionVector& curSol,
@@ -189,7 +189,7 @@ public:
{
// get some references for convenience
const auto& element = this->element();
- const auto& curSol = this->curSol()[domainId];
+ const auto& curSol = this->curSol(domainId);
auto&& fvGeometry = this->fvGeometry();
auto&& curElemVolVars = this->curElemVolVars();
auto&& elemFluxVarsCache = this->elemFluxVarsCache();
@@ -238,9 +238,9 @@ private:
/*!
* \ingroup Assembly
- * \ingroup BoxDiscretization
+ * \ingroup CVFEDiscretization
* \ingroup MultiDomain
- * \brief The box scheme multidomain local assembler
+ * \brief The CVFE scheme multidomain local assembler
* \tparam id the id of the sub domain
* \tparam TypeTag the TypeTag
* \tparam Assembler the assembler type
@@ -248,21 +248,21 @@ private:
* \tparam implicit whether the assembler is explicit or implicit in time
*/
template<std::size_t id, class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true>
-class SubDomainBoxLocalAssembler;
+class SubDomainCVFELocalAssembler;
/*!
* \ingroup Assembly
- * \ingroup BoxDiscretization
+ * \ingroup CVFEDiscretization
* \ingroup MultiDomain
- * \brief Box scheme multi domain local assembler using numeric differentiation and implicit time discretization
+ * \brief CVFE scheme multi domain local assembler using numeric differentiation and implicit time discretization
*/
template<std::size_t id, class TypeTag, class Assembler>
-class SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
-: public SubDomainBoxLocalAssemblerBase<id, TypeTag, Assembler,
- SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, true>, DiffMethod::numeric, true >
+class SubDomainCVFELocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
+: public SubDomainCVFELocalAssemblerBase<id, TypeTag, Assembler,
+ SubDomainCVFELocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, true>, DiffMethod::numeric, true >
{
- using ThisType = SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>;
- using ParentType = SubDomainBoxLocalAssemblerBase<id, TypeTag, Assembler, ThisType, DiffMethod::numeric, /*implicit=*/true>;
+ using ThisType = SubDomainCVFELocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>;
+ using ParentType = SubDomainCVFELocalAssemblerBase<id, TypeTag, Assembler, ThisType, DiffMethod::numeric, /*implicit=*/true>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
@@ -426,17 +426,17 @@ public:
/*!
* \ingroup Assembly
- * \ingroup BoxDiscretization
+ * \ingroup CVFEDiscretization
* \ingroup MultiDomain
- * \brief Box scheme multi domain local assembler using numeric differentiation and explicit time discretization
+ * \brief CVFE scheme multi domain local assembler using numeric differentiation and explicit time discretization
*/
template<std::size_t id, class TypeTag, class Assembler>
-class SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>
-: public SubDomainBoxLocalAssemblerBase<id, TypeTag, Assembler,
- SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, false>, DiffMethod::numeric, false >
+class SubDomainCVFELocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>
+: public SubDomainCVFELocalAssemblerBase<id, TypeTag, Assembler,
+ SubDomainCVFELocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, false>, DiffMethod::numeric, false >
{
- using ThisType = SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>;
- using ParentType = SubDomainBoxLocalAssemblerBase<id, TypeTag, Assembler, ThisType, DiffMethod::numeric, /*implicit=*/false>;
+ using ThisType = SubDomainCVFELocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>;
+ using ParentType = SubDomainCVFELocalAssemblerBase<id, TypeTag, Assembler, ThisType, DiffMethod::numeric, /*implicit=*/false>;
public:
using ParentType::ParentType;
diff --git a/dumux/multidomain/subdomainfcdiamondlocalassembler.hh b/dumux/multidomain/subdomainfcdiamondlocalassembler.hh
deleted file mode 100644
index bd5c3ef7d6..0000000000
--- a/dumux/multidomain/subdomainfcdiamondlocalassembler.hh
+++ /dev/null
@@ -1,418 +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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Assembly
- * \ingroup DiamondDiscretization
- * \ingroup MultiDomain
- * \brief An assembler for Jacobian and residual contribution per element (face-centered diamond methods) for multidomain problems
- */
-#ifndef DUMUX_MULTIDOMAIN_FACECENTERED_DIAMOND_LOCAL_ASSEMBLER_HH
-#define DUMUX_MULTIDOMAIN_FACECENTERED_DIAMOND_LOCAL_ASSEMBLER_HH
-
-#include <dune/common/indices.hh>
-#include <dune/common/hybridutilities.hh>
-#include <dune/grid/common/gridenums.hh> // for GhostEntity
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-#include <dumux/common/numericdifferentiation.hh>
-#include <dumux/assembly/numericepsilon.hh>
-#include <dumux/assembly/diffmethod.hh>
-#include <dumux/assembly/cvfelocalassembler.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Assembly
- * \ingroup DiamondDiscretization
- * \ingroup MultiDomain
- * \brief A base class for all face-centered staggered local assemblers
- * \tparam id the id of the sub domain
- * \tparam TypeTag the TypeTag
- * \tparam Assembler the assembler type
- * \tparam Implementation the actual implementation type
- * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
- */
-template<std::size_t id, class TypeTag, class Assembler, class Implementation, DiffMethod dm, bool implicit>
-class SubDomainFaceCenteredDiamondLocalAssemblerBase : public CVFELocalAssembler<TypeTag, Assembler, dm, implicit, Implementation>
-{
- using ParentType = CVFELocalAssembler<TypeTag, Assembler, dm, implicit, Implementation>;
-
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using SolutionVector = typename Assembler::SolutionVector;
-
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
- using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
- using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
- using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
- using Scalar = typename GridVariables::Scalar;
-
- using GridGeometry = typename GridVariables::GridGeometry;
- using FVElementGeometry = typename GridGeometry::LocalView;
- using SubControlVolume = typename GridGeometry::SubControlVolume;
- using GridView = typename GridGeometry::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
-
- using CouplingManager = typename Assembler::CouplingManager;
-
- static constexpr auto numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
-
-public:
- //! export the domain id of this sub-domain
- static constexpr auto domainId = typename Dune::index_constant<id>();
- //! pull up constructor of parent class
- using ParentType::ParentType;
- //! export element residual vector type
- using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
-
- // the constructor
- explicit SubDomainFaceCenteredDiamondLocalAssemblerBase(
- const Assembler& assembler,
- const Element& element,
- const SolutionVector& curSol,
- CouplingManager& couplingManager
- )
- : ParentType(assembler,
- element,
- curSol,
- localView(assembler.gridGeometry(domainId)),
- localView(assembler.gridVariables(domainId).curGridVolVars()),
- localView(assembler.gridVariables(domainId).prevGridVolVars()),
- localView(assembler.gridVariables(domainId).gridFluxVarsCache()),
- assembler.localResidual(domainId),
- (element.partitionType() == Dune::GhostEntity))
- , couplingManager_(couplingManager)
- {}
-
- /*!
- * \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<class JacobianMatrixRow, class SubResidual, class GridVariablesTuple>
- void assembleJacobianAndResidual(JacobianMatrixRow& jacRow, SubResidual& res, GridVariablesTuple& gridVariables)
- {
- auto assembleCouplingBlocks = [&](const auto& residual)
- {
- // assemble the coupling blocks
- using namespace Dune::Hybrid;
- forEach(integralRange(Dune::Hybrid::size(jacRow)), [&](auto&& i)
- {
- if (i != id)
- this->assembleJacobianCoupling(i, jacRow, residual, gridVariables);
- });
- };
-
- // the coupled model does not support partial reassembly yet
- const DefaultPartialReassembler* noReassembler = nullptr;
- ParentType::assembleJacobianAndResidual(jacRow[domainId], res, *std::get<domainId>(gridVariables), noReassembler, assembleCouplingBlocks);
- }
-
- /*!
- * \brief Assemble the entries in a coupling block of the jacobian.
- * There is no coupling block between a domain and itself.
- */
- template<std::size_t otherId, class JacRow, class GridVariables,
- typename std::enable_if_t<(otherId == id), int> = 0>
- void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow,
- const ElementResidualVector& res, GridVariables& gridVariables)
- {}
-
- /*!
- * \brief Assemble the entries in a coupling block of the jacobian.
- */
- template<std::size_t otherId, class JacRow, class GridVariables,
- typename std::enable_if_t<(otherId != id), int> = 0>
- void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow,
- const ElementResidualVector& res, GridVariables& gridVariables)
- {
- this->asImp_().assembleJacobianCoupling(domainJ, jacRow[domainJ], res, *std::get<domainId>(gridVariables));
- }
-
- /*!
- * \brief Evaluates the local source term for an element and given element volume variables
- */
- ElementResidualVector evalLocalSourceResidual(const Element& element, const ElementVolumeVariables& elemVolVars) const
- {
- // initialize the residual vector for all scvs in this element
- ElementResidualVector residual(this->fvGeometry().numScv());
-
- // evaluate the volume terms (storage + source terms)
- // forward to the local residual specialized for the discretization methods
- for (auto&& scv : scvs(this->fvGeometry()))
- {
- const auto& curVolVars = elemVolVars[scv];
- auto source = this->localResidual().computeSource(problem(), element, this->fvGeometry(), elemVolVars, scv);
- source *= -scv.volume()*curVolVars.extrusionFactor();
- residual[scv.localDofIndex()] = std::move(source);
- }
-
- return residual;
- }
-
- /*!
- * \brief Evaluates the local source term depending on time discretization scheme
- */
- ElementResidualVector evalLocalSourceResidual(const Element& neighbor) const
- { return this->evalLocalSourceResidual(neighbor, implicit ? this->curElemVolVars() : this->prevElemVolVars()); }
-
- /*!
- * \brief Prepares all local views necessary for local assembly.
- */
- void bindLocalViews()
- {
- // get some references for convenience
- const auto& element = this->element();
- const auto& curSol = this->curSol(domainId);
- auto&& fvGeometry = this->fvGeometry();
- auto&& curElemVolVars = this->curElemVolVars();
- auto&& elemFluxVarsCache = this->elemFluxVarsCache();
-
- // bind the caches
- couplingManager_.bindCouplingContext(domainId, element, this->assembler());
- fvGeometry.bind(element);
-
- if (implicit)
- {
- curElemVolVars.bind(element, fvGeometry, curSol);
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
- if (!this->assembler().isStationaryProblem())
- this->prevElemVolVars().bindElement(element, fvGeometry, this->assembler().prevSol()[domainId]);
- }
- else
- {
- auto& prevElemVolVars = this->prevElemVolVars();
- const auto& prevSol = this->assembler().prevSol()[domainId];
-
- curElemVolVars.bindElement(element, fvGeometry, curSol);
- prevElemVolVars.bind(element, fvGeometry, prevSol);
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
- }
-
- this->elemBcTypes().update(problem(), this->element(), this->fvGeometry());
- }
-
- //! return reference to the underlying problem
- template<std::size_t i = domainId>
- const Problem& problem(Dune::index_constant<i> dId = domainId) const
- { return this->assembler().problem(domainId); }
-
- //! return reference to the underlying problem
- template<std::size_t i = domainId>
- const auto& curSol(Dune::index_constant<i> dId = domainId) const
- { return ParentType::curSol()[dId]; }
-
- //! return reference to the coupling manager
- CouplingManager& couplingManager()
- { return couplingManager_; }
-
-private:
- CouplingManager& couplingManager_; //!< the coupling manager
-};
-
-/*!
- * \ingroup Assembly
- * \ingroup DiamondDiscretization
- * \ingroup MultiDomain
- * \brief The face-centered staggered scheme multidomain local assembler
- * \tparam id the id of the sub domain
- * \tparam TypeTag the TypeTag
- * \tparam Assembler the assembler type
- * \tparam DM the numeric differentiation method
- * \tparam implicit whether the assembler is explicit or implicit in time
- */
-template<std::size_t id, class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true>
-class SubDomainFaceCenteredDiamondLocalAssembler;
-
-/*!
- * \ingroup Assembly
- * \ingroup DiamondDiscretization
- * \ingroup MultiDomain
- * \brief Face-centered staggered scheme multi domain local assembler using numeric differentiation and implicit time discretization
- */
-template<std::size_t id, class TypeTag, class Assembler>
-class SubDomainFaceCenteredDiamondLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
-: public SubDomainFaceCenteredDiamondLocalAssemblerBase<id, TypeTag, Assembler,
- SubDomainFaceCenteredDiamondLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, true>, DiffMethod::numeric, /*implicit=*/true>
-{
- using ThisType = SubDomainFaceCenteredDiamondLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>;
- using ParentType = SubDomainFaceCenteredDiamondLocalAssemblerBase<id, TypeTag, Assembler, ThisType, DiffMethod::numeric, /*implicit=*/true>;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
-
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using FVElementGeometry = typename GridGeometry::LocalView;
- using GridView = typename GridGeometry::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
-
- enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
- enum { dim = GridView::dimension };
-
- static constexpr bool enableGridFluxVarsCache = GetPropType<TypeTag, Properties::GridVariables>::GridFluxVariablesCache::cachingEnabled;
- static constexpr bool enableGridVolVarsCache = GetPropType<TypeTag, Properties::GridVariables>::GridVolumeVariables::cachingEnabled;
- static constexpr auto domainI = Dune::index_constant<id>();
-
-public:
- using ParentType::ParentType;
- //! export element residual vector type
- using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
-
- /*!
- * \brief Update the coupling context for coupled models.
- */
- template<class ElemSol>
- void maybeUpdateCouplingContext(const SubControlVolume& scv, ElemSol& elemSol, const int pvIdx)
- {
- this->couplingManager().updateCouplingContext(domainI, *this, domainI, scv.dofIndex(), elemSol[scv.localDofIndex()], pvIdx);
- }
-
- /*!
- * \brief Update the additional domain derivatives for coupled models.
- */
- template<class JacobianMatrixDiagBlock, class GridVariables>
- void maybeEvalAdditionalDomainDerivatives(const ElementResidualVector& origResiduals, JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
- {
- this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origResiduals, A, gridVariables);
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- *
- * \return The element residual at the current solution.
- */
- template<std::size_t otherId, class JacobianBlock, class GridVariables>
- void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
- const ElementResidualVector& res, GridVariables& gridVariables)
- {
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
- // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
-
- // get some aliases for convenience
- const auto& element = this->element();
- const auto& fvGeometry = this->fvGeometry();
- auto&& curElemVolVars = this->curElemVolVars();
- auto&& elemFluxVarsCache = this->elemFluxVarsCache();
-
- // convenience lambda for call to update self
- auto updateCoupledVariables = [&] ()
- {
- // Update ourself after the context has been modified. Depending on the
- // type of caching, other objects might have to be updated. All ifs can be optimized away.
- if constexpr (enableGridFluxVarsCache)
- {
- if constexpr (enableGridVolVarsCache)
- this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
- else
- this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache());
- }
- else
- {
- if constexpr (enableGridVolVarsCache)
- this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache);
- else
- this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache);
- }
- };
-
- const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ);
- const auto& curSolJ = this->curSol(domainJ);
- for (const auto globalJ : stencil)
- {
- // undeflected privars and privars to be deflected
- const auto origPriVarsJ = curSolJ[globalJ];
- auto priVarsJ = origPriVarsJ;
-
- // the undeflected coupling residual
- const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ);
-
- for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx)
- {
- auto evalCouplingResidual = [&](Scalar priVar)
- {
- priVarsJ[pvIdx] = priVar;
- this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
- updateCoupledVariables();
- return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ);
- };
-
- // derive the residuals numerically
- ElementResidualVector partialDerivs(element.subEntities(1));
-
- const auto& paramGroup = this->assembler().problem(domainJ).paramGroup();
- static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod");
- static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup);
-
- NumericDifferentiation::partialDerivative(evalCouplingResidual, origPriVarsJ[pvIdx], partialDerivs, origResidual,
- epsCoupl(origPriVarsJ[pvIdx], pvIdx), numDiffMethod);
-
- // update the global stiffness matrix with the current partial derivatives
- for (const auto& scv : scvs(fvGeometry))
- {
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- // A[i][col][eqIdx][pvIdx] is the rate of change of
- // the residual of equation 'eqIdx' at dof 'i'
- // depending on the primary variable 'pvIdx' at dof
- // 'col'.
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] += partialDerivs[scv.localDofIndex()][eqIdx];
-
- // If the dof is coupled by a Dirichlet condition,
- // set the derived value only once (i.e. overwrite existing values).
- if (this->elemBcTypes().hasDirichlet())
- {
- const auto bcTypes = this->elemBcTypes().get(fvGeometry, scv);
- if (bcTypes.isCouplingDirichlet(eqIdx))
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = partialDerivs[scv.localDofIndex()][eqIdx];
- else if (bcTypes.isDirichlet(eqIdx))
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0;
- }
-
- // enforce internal Dirichlet constraints
- if constexpr (Problem::enableInternalDirichletConstraints())
- {
- const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
- if (internalDirichletConstraints[eqIdx])
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0;
- }
- }
- }
-
- // restore the current element solution
- priVarsJ[pvIdx] = origPriVarsJ[pvIdx];
-
- // restore the undeflected state of the coupling context
- this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
- }
-
- // Restore original state of the flux vars cache and/or vol vars.
- // This has to be done in case they depend on variables of domainJ before
- // we continue with the numeric derivative w.r.t the next globalJ. Otherwise,
- // the next "origResidual" will be incorrect.
- updateCoupledVariables();
- }
- }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/multidomain/subdomainpq1bubblelocalassembler.hh b/dumux/multidomain/subdomainpq1bubblelocalassembler.hh
deleted file mode 100644
index f76f5730a1..0000000000
--- a/dumux/multidomain/subdomainpq1bubblelocalassembler.hh
+++ /dev/null
@@ -1,418 +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 <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Assembly
- * \ingroup PQ1BubbleDiscretization
- * \ingroup MultiDomain
- * \brief An assembler for Jacobian and residual contribution per element for multidomain problems
- */
-#ifndef DUMUX_MULTIDOMAIN_PQ1BUBBLE_SUBDOMAIN_LOCAL_ASSEMBLER_HH
-#define DUMUX_MULTIDOMAIN_PQ1BUBBLE_SUBDOMAIN_LOCAL_ASSEMBLER_HH
-
-#include <dune/common/indices.hh>
-#include <dune/common/hybridutilities.hh>
-#include <dune/grid/common/gridenums.hh> // for GhostEntity
-
-#include <dumux/common/properties.hh>
-#include <dumux/common/parameters.hh>
-#include <dumux/common/numericdifferentiation.hh>
-#include <dumux/assembly/numericepsilon.hh>
-#include <dumux/assembly/diffmethod.hh>
-#include <dumux/assembly/cvfelocalassembler.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup Assembly
- * \ingroup PQ1BubbleDiscretization
- * \ingroup MultiDomain
- * \brief A base class for PQ1Bubble local assemblers
- * \tparam id the id of the sub domain
- * \tparam TypeTag the TypeTag
- * \tparam Assembler the assembler type
- * \tparam Implementation the actual implementation type
- * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
- */
-template<std::size_t id, class TypeTag, class Assembler, class Implementation, DiffMethod dm, bool implicit>
-class SubDomainPQ1BubbleLocalAssemblerBase : public CVFELocalAssembler<TypeTag, Assembler, dm, implicit, Implementation>
-{
- using ParentType = CVFELocalAssembler<TypeTag, Assembler, dm, implicit, Implementation>;
-
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using SolutionVector = typename Assembler::SolutionVector;
-
- using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
- using GridVolumeVariables = typename GridVariables::GridVolumeVariables;
- using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
- using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
- using Scalar = typename GridVariables::Scalar;
-
- using GridGeometry = typename GridVariables::GridGeometry;
- using FVElementGeometry = typename GridGeometry::LocalView;
- using SubControlVolume = typename GridGeometry::SubControlVolume;
- using GridView = typename GridGeometry::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
-
- using CouplingManager = typename Assembler::CouplingManager;
-
- static constexpr auto numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
-
-public:
- //! export the domain id of this sub-domain
- static constexpr auto domainId = typename Dune::index_constant<id>();
- //! pull up constructor of parent class
- using ParentType::ParentType;
- //! export element residual vector type
- using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
-
- // the constructor
- explicit SubDomainPQ1BubbleLocalAssemblerBase(
- const Assembler& assembler,
- const Element& element,
- const SolutionVector& curSol,
- CouplingManager& couplingManager
- )
- : ParentType(assembler,
- element,
- curSol,
- localView(assembler.gridGeometry(domainId)),
- localView(assembler.gridVariables(domainId).curGridVolVars()),
- localView(assembler.gridVariables(domainId).prevGridVolVars()),
- localView(assembler.gridVariables(domainId).gridFluxVarsCache()),
- assembler.localResidual(domainId),
- (element.partitionType() == Dune::GhostEntity))
- , couplingManager_(couplingManager)
- {}
-
- /*!
- * \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<class JacobianMatrixRow, class SubResidualVector, class GridVariablesTuple>
- void assembleJacobianAndResidual(JacobianMatrixRow& jacRow, SubResidualVector& res, GridVariablesTuple& gridVariables)
- {
- auto assembleCouplingBlocks = [&](const auto& residual)
- {
- // assemble the coupling blocks
- using namespace Dune::Hybrid;
- forEach(integralRange(Dune::Hybrid::size(jacRow)), [&](auto&& i)
- {
- if (i != id)
- this->assembleJacobianCoupling(i, jacRow, residual, gridVariables);
- });
- };
-
- // the coupled model does not support partial reassembly yet
- const DefaultPartialReassembler* noReassembler = nullptr;
- ParentType::assembleJacobianAndResidual(jacRow[domainId], res, *std::get<domainId>(gridVariables), noReassembler, assembleCouplingBlocks);
- }
-
- /*!
- * \brief Assemble the entries in a coupling block of the jacobian.
- * There is no coupling block between a domain and itself.
- */
- template<std::size_t otherId, class JacRow, class GridVariables,
- typename std::enable_if_t<(otherId == id), int> = 0>
- void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow,
- const ElementResidualVector& res, GridVariables& gridVariables)
- {}
-
- /*!
- * \brief Assemble the entries in a coupling block of the jacobian.
- */
- template<std::size_t otherId, class JacRow, class GridVariables,
- typename std::enable_if_t<(otherId != id), int> = 0>
- void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacRow& jacRow,
- const ElementResidualVector& res, GridVariables& gridVariables)
- {
- this->asImp_().assembleJacobianCoupling(domainJ, jacRow[domainJ], res, *std::get<domainId>(gridVariables));
- }
-
- /*!
- * \brief Evaluates the local source term for an element and given element volume variables
- */
- ElementResidualVector evalLocalSourceResidual(const Element& element, const ElementVolumeVariables& elemVolVars) const
- {
- // initialize the residual vector for all scvs in this element
- ElementResidualVector residual(this->fvGeometry().numScv());
-
- // evaluate the volume terms (storage + source terms)
- // forward to the local residual specialized for the discretization methods
- for (auto&& scv : scvs(this->fvGeometry()))
- {
- const auto& curVolVars = elemVolVars[scv];
- auto source = this->localResidual().computeSource(problem(), element, this->fvGeometry(), elemVolVars, scv);
- source *= -scv.volume()*curVolVars.extrusionFactor();
- residual[scv.localDofIndex()] = std::move(source);
- }
-
- return residual;
- }
-
- /*!
- * \brief Evaluates the local source term depending on time discretization scheme
- */
- ElementResidualVector evalLocalSourceResidual(const Element& neighbor) const
- { return this->evalLocalSourceResidual(neighbor, implicit ? this->curElemVolVars() : this->prevElemVolVars()); }
-
- /*!
- * \brief Prepares all local views necessary for local assembly.
- */
- void bindLocalViews()
- {
- // get some references for convenience
- const auto& element = this->element();
- const auto& curSol = this->curSol(domainId);
- auto&& fvGeometry = this->fvGeometry();
- auto&& curElemVolVars = this->curElemVolVars();
- auto&& elemFluxVarsCache = this->elemFluxVarsCache();
-
- // bind the caches
- couplingManager_.bindCouplingContext(domainId, element, this->assembler());
- fvGeometry.bind(element);
-
- if (implicit)
- {
- curElemVolVars.bind(element, fvGeometry, curSol);
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
- if (!this->assembler().isStationaryProblem())
- this->prevElemVolVars().bindElement(element, fvGeometry, this->assembler().prevSol()[domainId]);
- }
- else
- {
- auto& prevElemVolVars = this->prevElemVolVars();
- const auto& prevSol = this->assembler().prevSol()[domainId];
-
- curElemVolVars.bindElement(element, fvGeometry, curSol);
- prevElemVolVars.bind(element, fvGeometry, prevSol);
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
- }
-
- this->elemBcTypes().update(problem(), this->element(), this->fvGeometry());
- }
-
- //! return reference to the underlying problem
- template<std::size_t i = domainId>
- const Problem& problem(Dune::index_constant<i> dId = domainId) const
- { return this->assembler().problem(domainId); }
-
- //! return reference to the underlying problem
- template<std::size_t i = domainId>
- const auto& curSol(Dune::index_constant<i> dId = domainId) const
- { return ParentType::curSol()[dId]; }
-
- //! return reference to the coupling manager
- CouplingManager& couplingManager()
- { return couplingManager_; }
-
-private:
- CouplingManager& couplingManager_; //!< the coupling manager
-};
-
-/*!
- * \ingroup Assembly
- * \ingroup PQ1BubbleDiscretization
- * \ingroup MultiDomain
- * \brief The PQ1Bubble scheme multidomain local assembler
- * \tparam id the id of the sub domain
- * \tparam TypeTag the TypeTag
- * \tparam Assembler the assembler type
- * \tparam DM the numeric differentiation method
- * \tparam implicit whether the assembler is explicit or implicit in time
- */
-template<std::size_t id, class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true>
-class SubDomainPQ1BubbleLocalAssembler;
-
-/*!
- * \ingroup Assembly
- * \ingroup PQ1BubbleDiscretization
- * \ingroup MultiDomain
- * \brief Control-volume fe staggered scheme multi domain local assembler using numeric differentiation and implicit time discretization
- */
-template<std::size_t id, class TypeTag, class Assembler>
-class SubDomainPQ1BubbleLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
-: public SubDomainPQ1BubbleLocalAssemblerBase<id, TypeTag, Assembler,
- SubDomainPQ1BubbleLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, true>, DiffMethod::numeric, /*implicit=*/true>
-{
- using ThisType = SubDomainPQ1BubbleLocalAssembler<id, TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>;
- using ParentType = SubDomainPQ1BubbleLocalAssemblerBase<id, TypeTag, Assembler, ThisType, DiffMethod::numeric, /*implicit=*/true>;
- using Scalar = GetPropType<TypeTag, Properties::Scalar>;
- using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
-
- using Problem = GetPropType<TypeTag, Properties::Problem>;
- using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
- using FVElementGeometry = typename GridGeometry::LocalView;
- using GridView = typename GridGeometry::GridView;
- using Element = typename GridView::template Codim<0>::Entity;
- using SubControlVolume = typename FVElementGeometry::SubControlVolume;
-
- enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
- enum { dim = GridView::dimension };
-
- static constexpr bool enableGridFluxVarsCache = GetPropType<TypeTag, Properties::GridVariables>::GridFluxVariablesCache::cachingEnabled;
- static constexpr bool enableGridVolVarsCache = GetPropType<TypeTag, Properties::GridVariables>::GridVolumeVariables::cachingEnabled;
- static constexpr auto domainI = Dune::index_constant<id>();
-
-public:
- using ParentType::ParentType;
- //! export element residual vector type
- using ElementResidualVector = typename ParentType::LocalResidual::ElementResidualVector;
-
- /*!
- * \brief Update the coupling context for coupled models.
- */
- template<class ElemSol>
- void maybeUpdateCouplingContext(const SubControlVolume& scv, ElemSol& elemSol, const int pvIdx)
- {
- this->couplingManager().updateCouplingContext(domainI, *this, domainI, scv.dofIndex(), elemSol[scv.localDofIndex()], pvIdx);
- }
-
- /*!
- * \brief Update the additional domain derivatives for coupled models.
- */
- template<class JacobianMatrixDiagBlock, class GridVariables>
- void maybeEvalAdditionalDomainDerivatives(const ElementResidualVector& origResiduals, JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
- {
- this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origResiduals, A, gridVariables);
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- *
- * \return The element residual at the current solution.
- */
- template<std::size_t otherId, class JacobianBlock, class GridVariables>
- void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
- const ElementResidualVector& res, GridVariables& gridVariables)
- {
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
- // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
-
- // get some aliases for convenience
- const auto& element = this->element();
- const auto& fvGeometry = this->fvGeometry();
- auto&& curElemVolVars = this->curElemVolVars();
- auto&& elemFluxVarsCache = this->elemFluxVarsCache();
-
- // convenience lambda for call to update self
- auto updateCoupledVariables = [&] ()
- {
- // Update ourself after the context has been modified. Depending on the
- // type of caching, other objects might have to be updated. All ifs can be optimized away.
- if constexpr (enableGridFluxVarsCache)
- {
- if constexpr (enableGridVolVarsCache)
- this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
- else
- this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache());
- }
- else
- {
- if constexpr (enableGridVolVarsCache)
- this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache);
- else
- this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache);
- }
- };
-
- const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ);
- const auto& curSolJ = this->curSol(domainJ);
- for (const auto globalJ : stencil)
- {
- // undeflected privars and privars to be deflected
- const auto origPriVarsJ = curSolJ[globalJ];
- auto priVarsJ = origPriVarsJ;
-
- // the undeflected coupling residual
- const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ);
-
- for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx)
- {
- auto evalCouplingResidual = [&](Scalar priVar)
- {
- priVarsJ[pvIdx] = priVar;
- this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
- updateCoupledVariables();
- return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ);
- };
-
- // derive the residuals numerically
- ElementResidualVector partialDerivs(fvGeometry.numScv());
-
- const auto& paramGroup = this->assembler().problem(domainJ).paramGroup();
- static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod");
- static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup);
-
- NumericDifferentiation::partialDerivative(evalCouplingResidual, origPriVarsJ[pvIdx], partialDerivs, origResidual,
- epsCoupl(origPriVarsJ[pvIdx], pvIdx), numDiffMethod);
-
- // update the global stiffness matrix with the current partial derivatives
- for (const auto& scv : scvs(fvGeometry))
- {
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- // A[i][col][eqIdx][pvIdx] is the rate of change of
- // the residual of equation 'eqIdx' at dof 'i'
- // depending on the primary variable 'pvIdx' at dof
- // 'col'.
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] += partialDerivs[scv.localDofIndex()][eqIdx];
-
- // If the dof is coupled by a Dirichlet condition,
- // set the derived value only once (i.e. overwrite existing values).
- if (this->elemBcTypes().hasDirichlet())
- {
- const auto bcTypes = this->elemBcTypes().get(fvGeometry, scv);
- if (bcTypes.isCouplingDirichlet(eqIdx))
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = partialDerivs[scv.localDofIndex()][eqIdx];
- else if (bcTypes.isDirichlet(eqIdx))
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0;
- }
-
- // enforce internal Dirichlet constraints
- if constexpr (Problem::enableInternalDirichletConstraints())
- {
- const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
- if (internalDirichletConstraints[eqIdx])
- A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0;
- }
- }
- }
-
- // restore the current element solution
- priVarsJ[pvIdx] = origPriVarsJ[pvIdx];
-
- // restore the undeflected state of the coupling context
- this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
- }
-
- // Restore original state of the flux vars cache and/or vol vars.
- // This has to be done in case they depend on variables of domainJ before
- // we continue with the numeric derivative w.r.t the next globalJ. Otherwise,
- // the next "origResidual" will be incorrect.
- updateCoupledVariables();
- }
- }
-};
-
-} // end namespace Dumux
-
-#endif
--
GitLab