From 17422c49aa9eedbb5ef254f7662e19db10aaf436 Mon Sep 17 00:00:00 2001
From: Anna Mareike Kostelecky <anmako96@web.de>
Date: Wed, 2 Oct 2024 14:18:51 +0200
Subject: [PATCH] [multidomain][freeflow] add init function with prevSol, when
used as binary coupling manager
Co-authored-by: Timo Koch <timokoch@uio.no>
---
.../freeflowporenetwork/couplingmanager.hh | 87 +++++++++++++------
.../couplingmanager_base.hh | 70 +++++++++++----
.../freeflow/couplingmanager_staggered.hh | 53 ++++++++---
3 files changed, 155 insertions(+), 55 deletions(-)
diff --git a/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh b/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh
index 2386f22da2..2b4f338e09 100644
--- a/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh
+++ b/dumux/multidomain/boundary/freeflowporenetwork/couplingmanager.hh
@@ -173,6 +173,10 @@ class FreeFlowPoreNetworkCouplingManager
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
using SolutionVector = typename MDTraits::SolutionVector;
+ template<std::size_t id>
+ using SubSolutionVector
+ = std::decay_t<decltype(std::declval<SolutionVector>()[Dune::index_constant<id>()])>;
+
using CouplingConditions = FreeFlowPoreNetworkCouplingConditions<MDTraits, FreeFlowPoreNetworkCouplingManager<MDTraits>>;
using CouplingMapper = StaggeredFreeFlowPoreNetworkCouplingMapper;
@@ -195,38 +199,36 @@ public:
GridVarsTuple&& gridVarsTuple,
const SolutionVector& curSol)
{
- this->updateSolution(curSol); // generic coupling manager stores tuple of shared_ptr
-
- auto couplingMapper = std::make_shared<CouplingMapper>();
- couplingMapper->update(freeFlowMomentumProblem->gridGeometry(),
- freeFlowMassProblem->gridGeometry(),
- poreNetworkProblem->gridGeometry()
- );
+ // initialize sub coupling manager that are not stationary or transient problem specific
+ this->init_(freeFlowMomentumProblem, freeFlowMassProblem, poreNetworkProblem, std::forward<GridVarsTuple>(gridVarsTuple), curSol);
- // initialize the binary sub coupling managers
- typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage ffSolVecTuple;
- std::get<0>(ffSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
- std::get<1>(ffSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
+ // initialize stationary-specific sub coupling manager for free-flow
+ using FFSol = typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage;
this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).init(
freeFlowMomentumProblem, freeFlowMassProblem,
std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<freeFlowMassIndex>(gridVarsTuple)),
- ffSolVecTuple
+ FFSol{ std::get<freeFlowMomentumIndex>(this->curSol()), std::get<freeFlowMassIndex>(this->curSol()) }
);
+ }
- typename SubCouplingManager<freeFlowMassIndex, poreNetworkIndex>::SolutionVectorStorage ffMassPmSolVecTuple;
- std::get<0>(ffMassPmSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
- std::get<1>(ffMassPmSolVecTuple) = std::get<poreNetworkIndex>(this->curSol());
- this->subCouplingManager(freeFlowMassIndex, poreNetworkIndex).init(
- freeFlowMassProblem, poreNetworkProblem, couplingMapper, ffMassPmSolVecTuple
- );
+ template<class GridVarsTuple>
+ void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
+ std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
+ std::shared_ptr<Problem<poreNetworkIndex>> poreNetworkProblem,
+ GridVarsTuple&& gridVarsTuple,
+ const SolutionVector& curSol,
+ const SolutionVector& prevSol)
+ {
+ // initialize sub coupling manager that are not stationary or transient problem specific
+ this->init_(freeFlowMomentumProblem, freeFlowMassProblem, poreNetworkProblem, std::forward<GridVarsTuple>(gridVarsTuple), curSol);
- typename SubCouplingManager<freeFlowMomentumIndex, poreNetworkIndex>::SolutionVectorStorage ffMomentumPmSolVecTuple;
- std::get<0>(ffMomentumPmSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
- std::get<1>(ffMomentumPmSolVecTuple) = std::get<poreNetworkIndex>(this->curSol());
- this->subCouplingManager(freeFlowMomentumIndex, poreNetworkIndex).init(
- freeFlowMomentumProblem, poreNetworkProblem,
- std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<poreNetworkIndex>(gridVarsTuple)),
- couplingMapper, ffMomentumPmSolVecTuple
+ using FFSol = typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage;
+ using FFPrevSol = std::tuple<const SubSolutionVector<freeFlowMomentumIndex>*, const SubSolutionVector<freeFlowMassIndex>*>;
+ this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).init(
+ freeFlowMomentumProblem, freeFlowMassProblem,
+ std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<freeFlowMassIndex>(gridVarsTuple)),
+ FFSol{ std::get<freeFlowMomentumIndex>(this->curSol()), std::get<freeFlowMassIndex>(this->curSol()) },
+ FFPrevSol{ &prevSol[freeFlowMomentumIndex], &prevSol[freeFlowMassIndex] }
);
}
@@ -476,6 +478,41 @@ public:
return cm.couplingStencil(ii, elementI, scvI, jj);
});
}
+
+private:
+ /*
+ * \brief Initializes sub-coupling managers for stationary and transient problems
+ */
+ template<class GridVarsTuple>
+ void init_(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
+ std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
+ std::shared_ptr<Problem<poreNetworkIndex>> poreNetworkProblem,
+ GridVarsTuple&& gridVarsTuple,
+ const SolutionVector& curSol)
+ {
+ this->updateSolution(curSol); // generic coupling manager stores tuple of shared_ptr
+
+ auto couplingMapper = std::make_shared<CouplingMapper>();
+ couplingMapper->update(freeFlowMomentumProblem->gridGeometry(),
+ freeFlowMassProblem->gridGeometry(),
+ poreNetworkProblem->gridGeometry()
+ );
+
+ // initialize the binary sub coupling managers for stationary and transient problems
+ using FFMassPNMSol = typename SubCouplingManager<freeFlowMassIndex, poreNetworkIndex>::SolutionVectorStorage;
+ this->subCouplingManager(freeFlowMassIndex, poreNetworkIndex).init(
+ freeFlowMassProblem, poreNetworkProblem, couplingMapper,
+ FFMassPNMSol{ std::get<freeFlowMassIndex>(this->curSol()), std::get<poreNetworkIndex>(this->curSol()) }
+ );
+
+ using FFMomPNMSol = typename SubCouplingManager<freeFlowMomentumIndex, poreNetworkIndex>::SolutionVectorStorage;
+ this->subCouplingManager(freeFlowMomentumIndex, poreNetworkIndex).init(
+ freeFlowMomentumProblem, poreNetworkProblem,
+ std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<poreNetworkIndex>(gridVarsTuple)),
+ couplingMapper, FFMomPNMSol{ std::get<freeFlowMomentumIndex>(this->curSol()), std::get<poreNetworkIndex>(this->curSol()) }
+ );
+ }
+
};
} // end namespace Dumux
diff --git a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_base.hh b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_base.hh
index 664d051fcc..2c469b314f 100644
--- a/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_base.hh
+++ b/dumux/multidomain/boundary/freeflowporousmedium/couplingmanager_base.hh
@@ -170,6 +170,10 @@ class FreeFlowPorousMediumCouplingManagerBase
template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
using SolutionVector = typename MDTraits::SolutionVector;
+ template<std::size_t id>
+ using SubSolutionVector
+ = std::decay_t<decltype(std::declval<SolutionVector>()[Dune::index_constant<id>()])>;
+
public:
template<std::size_t i, std::size_t j>
@@ -189,30 +193,36 @@ public:
GridVarsTuple&& gridVarsTuple,
const SolutionVector& curSol)
{
- this->updateSolution(curSol); // generic coupling manager stores tuple of shared_ptr
+ // initialize sub coupling manager that are not stationary or transient problem specific
+ this->init_(freeFlowMomentumProblem, freeFlowMassProblem, porousMediumProblem, std::forward<GridVarsTuple>(gridVarsTuple), curSol);
- // initialize the binary sub coupling managers
- typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage ffSolVecTuple;
- std::get<0>(ffSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
- std::get<1>(ffSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
+ // initialize stationary-specific sub coupling manager for free-flow
+ using FFSol = typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage;
this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).init(
freeFlowMomentumProblem, freeFlowMassProblem,
std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<freeFlowMassIndex>(gridVarsTuple)),
- ffSolVecTuple
+ FFSol{ std::get<freeFlowMomentumIndex>(this->curSol()), std::get<freeFlowMassIndex>(this->curSol()) }
);
+ }
- typename SubCouplingManager<freeFlowMassIndex, porousMediumIndex>::SolutionVectorStorage ffMassPmSolVecTuple;
- std::get<0>(ffMassPmSolVecTuple) = std::get<freeFlowMassIndex>(this->curSol());
- std::get<1>(ffMassPmSolVecTuple) = std::get<porousMediumIndex>(this->curSol());
- this->subCouplingManager(freeFlowMassIndex, porousMediumIndex).init(
- freeFlowMassProblem, porousMediumProblem, ffMassPmSolVecTuple
- );
+ template<class GridVarsTuple>
+ void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
+ std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
+ std::shared_ptr<Problem<porousMediumIndex>> porousMediumProblem,
+ GridVarsTuple&& gridVarsTuple,
+ const SolutionVector& curSol,
+ const SolutionVector& prevSol)
+ {
+ // initialize sub coupling manager that are not stationary or transient problem specific
+ this->init_(freeFlowMomentumProblem, freeFlowMassProblem, porousMediumProblem, std::forward<GridVarsTuple>(gridVarsTuple), curSol);
- typename SubCouplingManager<freeFlowMomentumIndex, porousMediumIndex>::SolutionVectorStorage ffMomentumPmSolVecTuple;
- std::get<0>(ffMomentumPmSolVecTuple) = std::get<freeFlowMomentumIndex>(this->curSol());
- std::get<1>(ffMomentumPmSolVecTuple) = std::get<porousMediumIndex>(this->curSol());
- this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).init(
- freeFlowMomentumProblem, porousMediumProblem, ffMomentumPmSolVecTuple
+ using FFSol = typename SubCouplingManager<freeFlowMomentumIndex, freeFlowMassIndex>::SolutionVectorStorage;
+ using FFPrevSol = std::tuple<const SubSolutionVector<freeFlowMomentumIndex>*, const SubSolutionVector<freeFlowMassIndex>*>;
+ this->subCouplingManager(freeFlowMomentumIndex, freeFlowMassIndex).init(
+ freeFlowMomentumProblem, freeFlowMassProblem,
+ std::make_tuple(std::get<freeFlowMomentumIndex>(gridVarsTuple), std::get<freeFlowMassIndex>(gridVarsTuple)),
+ FFSol{ std::get<freeFlowMomentumIndex>(this->curSol()), std::get<freeFlowMassIndex>(this->curSol()) },
+ FFPrevSol{ &prevSol[freeFlowMomentumIndex], &prevSol[freeFlowMassIndex] }
);
}
@@ -271,6 +281,32 @@ public:
return cm.couplingStencil(ii, elementI, scvI, jj);
});
}
+private:
+ /*
+ * \brief Initializes sub-coupling managers for stationary and transient problems
+ */
+ template<class GridVarsTuple>
+ void init_(std::shared_ptr<Problem<freeFlowMomentumIndex>> freeFlowMomentumProblem,
+ std::shared_ptr<Problem<freeFlowMassIndex>> freeFlowMassProblem,
+ std::shared_ptr<Problem<porousMediumIndex>> porousMediumProblem,
+ GridVarsTuple&& gridVarsTuple,
+ const SolutionVector& curSol)
+ {
+ this->updateSolution(curSol); // generic coupling manager stores tuple of shared_ptr
+
+ // initialize the binary sub coupling managers
+ using FFMassPMSol = typename SubCouplingManager<freeFlowMassIndex, porousMediumIndex>::SolutionVectorStorage;
+ this->subCouplingManager(freeFlowMassIndex, porousMediumIndex).init(
+ freeFlowMassProblem, porousMediumProblem,
+ FFMassPMSol{ std::get<freeFlowMassIndex>(this->curSol()), std::get<porousMediumIndex>(this->curSol()) }
+ );
+
+ using FFMomPMSol = typename SubCouplingManager<freeFlowMomentumIndex, porousMediumIndex>::SolutionVectorStorage;
+ this->subCouplingManager(freeFlowMomentumIndex, porousMediumIndex).init(
+ freeFlowMomentumProblem, porousMediumProblem,
+ FFMomPMSol{ std::get<freeFlowMomentumIndex>(this->curSol()), std::get<porousMediumIndex>(this->curSol()) }
+ );
+ }
};
} // end namespace Dumux
diff --git a/dumux/multidomain/freeflow/couplingmanager_staggered.hh b/dumux/multidomain/freeflow/couplingmanager_staggered.hh
index d80e34c868..f2d03395dc 100644
--- a/dumux/multidomain/freeflow/couplingmanager_staggered.hh
+++ b/dumux/multidomain/freeflow/couplingmanager_staggered.hh
@@ -73,6 +73,15 @@ private:
using Scalar = typename Traits::Scalar;
using SolutionVector = typename Traits::SolutionVector;
+ template<std::size_t id>
+ using SubSolutionVector
+ = std::decay_t<decltype(std::declval<SolutionVector>()[Dune::index_constant<id>()])>;
+
+ template<std::size_t id>
+ using ConstSubSolutionVectorPtr = const SubSolutionVector<id>*;
+
+ using PrevSolutionVectorStorage = typename Traits::template Tuple<ConstSubSolutionVectorPtr>;
+
using CouplingStencilType = std::vector<std::size_t>;
using GridVariablesTuple = typename Traits::template TupleOfSharedPtr<GridVariables>;
@@ -131,8 +140,9 @@ public:
const SolutionVector& prevSol)
{
init(momentumProblem, massProblem, std::forward<GridVariablesTuple>(gridVariables), curSol);
- prevSol_ = &prevSol;
- isTransient_ = true;
+
+ Dune::Hybrid::forEach(std::make_index_sequence<Traits::numSubDomains>{}, [&](auto i)
+ { std::get<i>(prevSolutions_) = &prevSol[i]; });
}
//! use as binary coupling manager in multi model context
@@ -148,6 +158,17 @@ public:
computeCouplingStencils_();
}
+ //! use as binary coupling manager in multi model context and for transient problems
+ void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem,
+ std::shared_ptr<Problem<freeFlowMassIndex>> massProblem,
+ GridVariablesTuple&& gridVariables,
+ const typename ParentType::SolutionVectorStorage& curSol,
+ const PrevSolutionVectorStorage& prevSol)
+ {
+ init(momentumProblem, massProblem, std::forward<GridVariablesTuple>(gridVariables), curSol);
+ prevSolutions_ = prevSol;
+ }
+
// \}
using CouplingManager<Traits>::evalCouplingResidual;
@@ -191,7 +212,7 @@ public:
if (!localAssemblerI.assembler().isStationaryProblem())
{
- assert(isTransient_);
+ assert(isTransient_());
localResidual.evalStorage(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, scvI);
}
@@ -234,7 +255,7 @@ public:
const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
const bool considerPreviousTimeStep = false) const
{
- assert(!(considerPreviousTimeStep && !isTransient_));
+ assert(!(considerPreviousTimeStep && !isTransient_()));
bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex());
const auto& insideMomentumScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(insideMomentumScv.elementIndex());
@@ -261,7 +282,7 @@ public:
const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
const bool considerPreviousTimeStep = false) const
{
- assert(!(considerPreviousTimeStep && !isTransient_));
+ assert(!(considerPreviousTimeStep && !isTransient_()));
bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, fvGeometry.elementIndex());
const auto& insideMomentumScv = fvGeometry.scv(scvf.insideScvIdx());
const auto& insideMassScv = momentumCouplingContext_()[0].fvGeometry.scv(insideMomentumScv.elementIndex());
@@ -289,7 +310,7 @@ public:
const SubControlVolume<freeFlowMomentumIndex>& scv,
const bool considerPreviousTimeStep = false) const
{
- assert(!(considerPreviousTimeStep && !isTransient_));
+ assert(!(considerPreviousTimeStep && !isTransient_()));
bindCouplingContext_(Dune::index_constant<freeFlowMomentumIndex>(), element, scv.elementIndex());
const auto& massScv = (*scvs(momentumCouplingContext_()[0].fvGeometry).begin());
@@ -508,11 +529,11 @@ private:
auto curElemVolVars = localView(gridVars_(freeFlowMassIndex).curGridVolVars());
curElemVolVars.bind(elementI, fvGeometry, this->curSol(freeFlowMassIndex));
- auto prevElemVolVars = isTransient_ ? localView(gridVars_(freeFlowMassIndex).prevGridVolVars())
+ auto prevElemVolVars = isTransient_() ? localView(gridVars_(freeFlowMassIndex).prevGridVolVars())
: localView(gridVars_(freeFlowMassIndex).curGridVolVars());
- if (isTransient_)
- prevElemVolVars.bindElement(elementI, fvGeometry, (*prevSol_)[freeFlowMassIndex]);
+ if (isTransient_())
+ prevElemVolVars.bindElement(elementI, fvGeometry, prevSol_(freeFlowMassIndex));
momentumCouplingContext_().emplace_back(MomentumCouplingContext{std::move(fvGeometry), std::move(curElemVolVars), std::move(prevElemVolVars), eIdx});
}
@@ -522,8 +543,8 @@ private:
momentumCouplingContext_()[0].fvGeometry.bind(elementI);
momentumCouplingContext_()[0].curElemVolVars.bind(elementI, momentumCouplingContext_()[0].fvGeometry, this->curSol(freeFlowMassIndex));
- if (isTransient_)
- momentumCouplingContext_()[0].prevElemVolVars.bindElement(elementI, momentumCouplingContext_()[0].fvGeometry, (*prevSol_)[freeFlowMassIndex]);
+ if (isTransient_())
+ momentumCouplingContext_()[0].prevElemVolVars.bindElement(elementI, momentumCouplingContext_()[0].fvGeometry, prevSol_(freeFlowMassIndex));
}
}
@@ -661,8 +682,14 @@ private:
//! A tuple of std::shared_ptrs to the grid variables of the sub problems
GridVariablesTuple gridVariables_;
- const SolutionVector* prevSol_;
- bool isTransient_;
+ bool isTransient_() const
+ { return std::get<0>(prevSolutions_) != nullptr; }
+
+ template<std::size_t i>
+ const SubSolutionVector<i>& prevSol_(Dune::index_constant<i>) const
+ { return *std::get<i>(prevSolutions_); }
+
+ PrevSolutionVectorStorage prevSolutions_;
std::deque<std::vector<ElementSeed<freeFlowMomentumIndex>>> elementSets_;
};
--
GitLab