From 8233791200e9174375986767176a3a6679ab42ed Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 7 Sep 2021 12:02:41 +0200
Subject: [PATCH] [md] Add periodic constraints to multidomain fvassembler
---
dumux/multidomain/fvassembler.hh | 45 ++++++++++++++++++++++++++++++++
1 file changed, 45 insertions(+)
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index 8039452e79..8f88cc3d37 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -203,6 +203,9 @@ public:
auto& jacRow = (*jacobian_)[domainId];
auto& subRes = (*residual_)[domainId];
this->assembleJacobianAndResidual_(domainId, jacRow, subRes, curSol);
+
+ const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
+ enforcePeriodicConstraints_(domainId, jacRow, subRes, *gridGeometry, curSol[domainId]);
});
}
@@ -551,6 +554,48 @@ private:
domainJ, gridGeometry(domainJ));
}
+ // build periodic contraints into the system matrix
+ template<std::size_t i, class JacRow, class Sol, class GG>
+ void enforcePeriodicConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Sol& res, const GG& gridGeometry, const Sol& curSol)
+ {
+ if constexpr (GG::discMethod == DiscretizationMethod::box || GG::discMethod == DiscretizationMethod::fcstaggered)
+ {
+ for (const auto& m : gridGeometry.periodicVertexMap())
+ {
+ if (m.first < m.second)
+ {
+ auto& jac = jacRow[domainI];
+
+ // add the second row to the first
+ res[m.first] += res[m.second];
+
+ // enforce the solution of the first periodic DOF to the second one
+ res[m.second] = curSol[m.second] - curSol[m.first];
+
+ const auto end = jac[m.second].end();
+ for (auto it = jac[m.second].begin(); it != end; ++it)
+ jac[m.first][it.index()] += (*it);
+
+ // enforce constraint in second row
+ for (auto it = jac[m.second].begin(); it != end; ++it)
+ (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0;
+
+ using namespace Dune::Hybrid;
+ forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
+ {
+ auto& jacCoupling = jacRow[couplingDomainId];
+
+ for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
+ jacCoupling[m.first][it.index()] += (*it);
+
+ for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
+ (*it) = 0.0;
+ });
+ }
+ }
+ }
+ }
+
//! pointer to the problem to be solved
ProblemTuple problemTuple_;
--
GitLab