From 35b5f9e6ce90f86dbb1e52ee180734940d247923 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 6 Feb 2015 13:50:57 +0000
Subject: [PATCH] [el2p, multidomain] fix residual calculation, generalize
multidomain convergence writer
- The residual calculation for el2p was wrong, probably never used. It now
forwards the call to the GridOperator functionality.
- The multidomain convergence writer assumed that sub-domain solution
vectors can be default-contructed. This is not the case if a PDELab
model is used for the subdomain.
Reviewed by Dennis.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14200 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/geomechanics/el2p/el2passembler.hh | 2 +
dumux/geomechanics/el2p/el2pbasemodel.hh | 62 +++++--------------
.../common/multidomainconvergencewriter.hh | 20 ++----
3 files changed, 25 insertions(+), 59 deletions(-)
diff --git a/dumux/geomechanics/el2p/el2passembler.hh b/dumux/geomechanics/el2p/el2passembler.hh
index dfad36b3fd..b7d11e2998 100644
--- a/dumux/geomechanics/el2p/el2passembler.hh
+++ b/dumux/geomechanics/el2p/el2passembler.hh
@@ -629,6 +629,8 @@ private:
{ return problem_().vertexMapper(); }
const ElementMapper &elementMapper_() const
{ return problem_().elementMapper(); }
+ const GridOperator &gridOperator() const
+ { return *gridOperator_;}
Problem *problemPtr_;
diff --git a/dumux/geomechanics/el2p/el2pbasemodel.hh b/dumux/geomechanics/el2p/el2pbasemodel.hh
index 4bcf558b93..31ca84ce22 100644
--- a/dumux/geomechanics/el2p/el2pbasemodel.hh
+++ b/dumux/geomechanics/el2p/el2pbasemodel.hh
@@ -225,46 +225,7 @@ public:
Scalar globalResidual(SolutionVector &residual,
const SolutionVector &u)
{
- SolutionVector tmp(curSol());
- curSol() = u;
- Scalar res = globalResidual(residual);
- curSol() = tmp;
- return res;
- }
-
- /*!
- * \brief Compute the global residual for the current solution
- * vector.
- *
- * \param residual Stores the result
- */
- Scalar globalResidual(SolutionVector &residual)
- {
- residual = 0;
-
- ElementIterator eIt = gridView_().template begin<0>();
- const ElementIterator eEndIt = gridView_().template end<0>();
- for (; eIt != eEndIt; ++eIt) {
- localResidual().eval(*eIt);
-
- if (isBox)
- {
-#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
- for (int i = 0; i < eIt->subEntities(dim); ++i)
-#else
- for (int i = 0; i < eIt->template count<dim>(); ++i)
-#endif
- {
- int globalI = vertexMapper().map(*eIt, i, dim);
- residual[globalI] += localResidual().residual(i);
- }
- }
- else
- {
- int globalI = elementMapper().map(*eIt);
- residual[globalI] = localResidual().residual(0);
- }
- }
+ jacAsm_.gridOperator().residual(u, residual);
// calculate the square norm of the residual
Scalar result2 = residual.two_norm2();
@@ -283,6 +244,17 @@ public:
return std::sqrt(result2);
}
+ /*!
+ * \brief Compute the global residual for the current solution
+ * vector.
+ *
+ * \param residual Stores the result
+ */
+ Scalar globalResidual(SolutionVector &residual)
+ {
+ return globalResidual(residual, curSol());
+ }
+
/*!
* \brief Compute the integral over the domain of the storage
* terms of all conservation quantities.
@@ -711,13 +683,13 @@ public:
def[eqIdx] = writer.allocateManagedBuffer(numDofs);
}
- for (unsigned int vIdxGlobal = 0; vIdxGlobal < u.size(); vIdxGlobal++)
+ for (unsigned int vIdxGlobal = 0; vIdxGlobal < u.base().size(); vIdxGlobal++)
{
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
- (*x[eqIdx])[vIdxGlobal] = u[vIdxGlobal][eqIdx];
- (*delta[eqIdx])[vIdxGlobal] = - deltaU[vIdxGlobal][eqIdx];
- (*def[eqIdx])[vIdxGlobal] = residual[vIdxGlobal][eqIdx];
+ (*x[eqIdx])[vIdxGlobal] = u.base()[vIdxGlobal][eqIdx];
+ (*delta[eqIdx])[vIdxGlobal] = - deltaU.base()[vIdxGlobal][eqIdx];
+ (*def[eqIdx])[vIdxGlobal] = residual.base()[vIdxGlobal][eqIdx];
}
}
diff --git a/dumux/multidomain/common/multidomainconvergencewriter.hh b/dumux/multidomain/common/multidomainconvergencewriter.hh
index 68816b7630..b49b3f27d0 100644
--- a/dumux/multidomain/common/multidomainconvergencewriter.hh
+++ b/dumux/multidomain/common/multidomainconvergencewriter.hh
@@ -116,21 +116,13 @@ struct MultiDomainConvergenceWriter
void writeFields(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
- SolutionVector1 uLastIter1;
- SolutionVector2 uLastIter2;
- SolutionVector1 deltaU1;
- SolutionVector2 deltaU2;
-
- uLastIter1.resize(ctl_.method().model().sdModel1().numDofs());
- uLastIter2.resize(ctl_.method().model().sdModel2().numDofs());
- deltaU1.resize(ctl_.method().model().sdModel1().numDofs());
- deltaU2.resize(ctl_.method().model().sdModel2().numDofs());
-
- typedef Dumux::SplitAndMerge<TypeTag> Common;
-
- Common::splitSolVector(uLastIter, uLastIter1, uLastIter2);
- Common::splitSolVector(deltaU, deltaU1, deltaU2);
+ SolutionVector1 uLastIter1(ctl_.method().model().sdModel1().curSol());
+ SolutionVector2 uLastIter2(ctl_.method().model().sdModel2().curSol());
+ SolutionVector1 deltaU1(uLastIter1);
+ SolutionVector2 deltaU2(uLastIter2);
+ SplitAndMerge<TypeTag>::splitSolVector(uLastIter, uLastIter1, uLastIter2);
+ SplitAndMerge<TypeTag>::splitSolVector(deltaU, deltaU1, deltaU2);
std::cout << "\n writing convergence file of current Newton iteration \n";
ctl_.method().model().sdModel1().addConvergenceVtkFields(*vtkMultiWriter1_, uLastIter1, deltaU1);
--
GitLab