From 0196fdfaa6ef172e57c01d57bbd429a543717428 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Tue, 23 Oct 2012 15:27:09 +0000
Subject: [PATCH] Enable parallelism for decoupled 2p2c. Requires AMGBackend
and patches for pdelab and istl. Can therefore only be used from devel
(currently). Reviewed by Benjamin.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@9414 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/decoupled/2p2c/fvtransport2p2c.hh | 17 ++++
dumux/decoupled/common/fv/fvpressure.hh | 125 +++++++++++++-----------
dumux/linear/vectorexchange.hh | 90 +++++++++++++++++
3 files changed, 176 insertions(+), 56 deletions(-)
create mode 100644 dumux/linear/vectorexchange.hh
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index 4514d0ecf4..01c30333a7 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -23,6 +23,7 @@
#include <dumux/decoupled/2p2c/2p2cproperties.hh>
#include <dumux/material/constraintsolvers/compositionalflash.hh>
#include <dumux/common/math.hh>
+#include <dumux/linear/vectorexchange.hh>
/**
* @file
@@ -300,6 +301,22 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt,
restrictingCell= globalIdxI;
}
} // end grid traversal
+
+#if HAVE_MPI
+ // communicate updated values
+ typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
+ typedef typename SolutionTypes::ElementMapper ElementMapper;
+ typedef VectorExchange<ElementMapper, Dune::BlockVector<Dune::FieldVector<Scalar, 1> > > DataHandle;
+ for (int i = 0; i < updateVec.size(); i++)
+ {
+ DataHandle dataHandle(problem_.variables().elementMapper(), updateVec[i]);
+ problem_.gridView().template communicate<DataHandle>(dataHandle,
+ Dune::InteriorBorder_All_Interface,
+ Dune::ForwardCommunication);
+ }
+ dt = problem_.gridView().comm().min(dt);
+#endif
+
if(impet)
{
Dune::dinfo << "Timestep restricted by CellIdx " << restrictingCell << " leads to dt = "
diff --git a/dumux/decoupled/common/fv/fvpressure.hh b/dumux/decoupled/common/fv/fvpressure.hh
index 2fec2e5975..1b05b20502 100644
--- a/dumux/decoupled/common/fv/fvpressure.hh
+++ b/dumux/decoupled/common/fv/fvpressure.hh
@@ -366,77 +366,90 @@ void FVPressure<TypeTag>::assemble(bool first)
ElementIterator eItEnd = problem_.gridView().template end<0>();
for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
{
- // cell information
+ // get the global index of the cell
int globalIdxI = problem_.variables().index(*eIt);
- CellData& cellDataI = problem_.variables().cellData(globalIdxI);
- EntryType entries(0.);
+ // assemble interior element contributions
+ if (eIt->partitionType() == Dune::InteriorEntity)
+ {
+ // get the cell data
+ CellData& cellDataI = problem_.variables().cellData(globalIdxI);
- /***** source term ***********/
- asImp_().getSource(entries, *eIt, cellDataI, first);
- f_[globalIdxI] += entries[rhs];
+ EntryType entries(0.);
- /***** flux term ***********/
- // iterate over all faces of the cell
- IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
- {
- /************* handle interior face *****************/
- if (isIt->neighbor())
+ /***** source term ***********/
+ asImp_().getSource(entries, *eIt, cellDataI, first);
+ f_[globalIdxI] += entries[rhs];
+
+ /***** flux term ***********/
+ // iterate over all faces of the cell
+ IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
{
- ElementPointer elementNeighbor = isIt->outside();
+ /************* handle interior face *****************/
+ if (isIt->neighbor())
+ {
+ ElementPointer elementNeighbor = isIt->outside();
- int globalIdxJ = problem_.variables().index(*elementNeighbor);
+ int globalIdxJ = problem_.variables().index(*elementNeighbor);
- //check for hanging nodes
- //take a hanging node never from the element with smaller level!
- bool haveSameLevel = (eIt->level() == elementNeighbor->level());
- //calculate only from one side, but add matrix entries for both sides
- if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce) && (globalIdxI > globalIdxJ) && haveSameLevel)
- continue;
+ //check for hanging nodes
+ //take a hanging node never from the element with smaller level!
+ bool haveSameLevel = (eIt->level() == elementNeighbor->level());
+ //calculate only from one side, but add matrix entries for both sides
+ if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce) && (globalIdxI > globalIdxJ) && haveSameLevel)
+ continue;
- //check for hanging nodes
- entries = 0;
- asImp_().getFlux(entries, *isIt, cellDataI, first);
+ //check for hanging nodes
+ entries = 0;
+ asImp_().getFlux(entries, *isIt, cellDataI, first);
- //set right hand side
- f_[globalIdxI] -= entries[rhs];
+ //set right hand side
+ f_[globalIdxI] -= entries[rhs];
- // set diagonal entry
- A_[globalIdxI][globalIdxI] += entries[matrix];
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entries[matrix];
- // set off-diagonal entry
- A_[globalIdxI][globalIdxJ] -= entries[matrix];
+ // set off-diagonal entry
+ A_[globalIdxI][globalIdxJ] -= entries[matrix];
- if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce))
- {
- f_[globalIdxJ] += entries[rhs];
- A_[globalIdxJ][globalIdxJ] += entries[matrix];
- A_[globalIdxJ][globalIdxI] -= entries[matrix];
- }
+ if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce))
+ {
+ f_[globalIdxJ] += entries[rhs];
+ A_[globalIdxJ][globalIdxJ] += entries[matrix];
+ A_[globalIdxJ][globalIdxI] -= entries[matrix];
+ }
- } // end neighbor
+ } // end neighbor
- /************* boundary face ************************/
- else
- {
- entries = 0;
- asImp_().getFluxOnBoundary(entries, *isIt, cellDataI, first);
+ /************* boundary face ************************/
+ else
+ {
+ entries = 0;
+ asImp_().getFluxOnBoundary(entries, *isIt, cellDataI, first);
- //set right hand side
- f_[globalIdxI] += entries[rhs];
- // set diagonal entry
- A_[globalIdxI][globalIdxI] += entries[matrix];
- }
- } //end interfaces loop
-// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
-
- /***** storage term ***********/
- entries = 0;
- asImp_().getStorage(entries, *eIt, cellDataI, first);
- f_[globalIdxI] += entries[rhs];
-// set diagonal entry
- A_[globalIdxI][globalIdxI] += entries[matrix];
+ //set right hand side
+ f_[globalIdxI] += entries[rhs];
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entries[matrix];
+ }
+ } //end interfaces loop
+ // printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
+
+ /***** storage term ***********/
+ entries = 0;
+ asImp_().getStorage(entries, *eIt, cellDataI, first);
+ f_[globalIdxI] += entries[rhs];
+ // set diagonal entry
+ A_[globalIdxI][globalIdxI] += entries[matrix];
+ }
+ // assemble overlap and ghost element contributions
+ else
+ {
+ A_[globalIdxI] = 0.0;
+ A_[globalIdxI][globalIdxI] = 1.0;
+ f_[globalIdxI] = pressure_[globalIdxI];
+ }
} // end grid traversal
// printmatrix(std::cout, A_, "global stiffness matrix after assempling", "row", 11,3);
// printvector(std::cout, f_, "right hand side", "row", 10);
diff --git a/dumux/linear/vectorexchange.hh b/dumux/linear/vectorexchange.hh
new file mode 100644
index 0000000000..ece6fc1104
--- /dev/null
+++ b/dumux/linear/vectorexchange.hh
@@ -0,0 +1,90 @@
+// -*- 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 2 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/>. *
+ *****************************************************************************/
+#ifndef DUMUX_VECTOR_EXCHANGE_HH
+#define DUMUX_VECTOR_EXCHANGE_HH
+
+#include <dune/grid/common/datahandleif.hh>
+
+namespace Dumux
+{
+
+// A DataHandle class to exchange entries of a vector
+template<class Mapper, class Vector> // mapper type and vector type
+class VectorExchange
+ : public Dune::CommDataHandleIF<VectorExchange<Mapper,Vector>,
+ typename Vector::value_type>
+{
+public:
+ //! export type of data for message buffer
+ typedef typename Vector::value_type DataType;
+
+ //! returns true if data for this codim should be communicated
+ bool contains (int dim, int codim) const
+ {
+ return (codim == 0);
+ }
+
+ //! returns true if size per entity of given dim and codim is a constant
+ bool fixedsize (int dim, int codim) const
+ {
+ return true;
+ }
+
+ /*! how many objects of type DataType have to be sent for a given entity
+
+ Note: Only the sender side needs to know this size.
+ */
+ template<class Entity>
+ size_t size (Entity& entity) const
+ {
+ return 1;
+ }
+
+ //! pack data from user to message buffer
+ template<class MessageBuffer, class Entity>
+ void gather (MessageBuffer& buff, const Entity& entity) const
+ {
+ buff.write(dataVector_[mapper_.map(entity)]);
+ }
+
+ /*! unpack data from message buffer to user
+
+ n is the number of objects sent by the sender
+ */
+ template<class MessageBuffer, class Entity>
+ void scatter (MessageBuffer& buff, const Entity& entity, size_t n)
+ {
+ DataType x;
+ buff.read(x);
+ dataVector_[mapper_.map(entity)] = x;
+ }
+
+ //! constructor
+ VectorExchange (const Mapper& mapper, Vector& dataVector)
+ : mapper_(mapper), dataVector_(dataVector)
+ {}
+
+private:
+ const Mapper& mapper_;
+ Vector& dataVector_;
+};
+
+}
+
+#endif
--
GitLab