Skip to content

GitLab

Commit ffaf6bd9 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

Parallelize also adaptive and multiphysics 2p2c. Reviewed by Benjamin. 

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@9433 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent a9d4f2d3
Hide whitespace changes
Inline Side-by-side
dumux/decoupled/2p2c/fvpressure2p2cadaptive.hh
View file @ ffaf6bd9
......@@ -392,84 +392,104 @@ void FVPressure2P2CAdaptive<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);
Dune::FieldVector<Scalar, 2> entries(0.);
// assemble interior element contributions
if (eIt->partitionType() == Dune::InteriorEntity)
{
// get the cell data
CellData& cellDataI = problem().variables().cellData(globalIdxI);
/***** source term ***********/
problem().pressureModel().getSource(entries, *eIt, cellDataI, first);
this->f_[globalIdxI] += entries[rhs];
Dune::FieldVector<Scalar, 2> entries(0.);
/***** flux term ***********/
// iterate over all faces of the cell
IntersectionIterator isItEnd = problem().gridView().template iend(*eIt);
for (IntersectionIterator isIt = problem().gridView().template ibegin(*eIt); isIt != isItEnd; ++isIt)
{
/************* handle interior face *****************/
if (isIt->neighbor())
/***** source term ***********/
problem().pressureModel().getSource(entries, *eIt, cellDataI, first);
this->f_[globalIdxI] += entries[rhs];
/***** flux term ***********/
// iterate over all faces of the cell
IntersectionIterator isItEnd = problem().gridView().template iend(*eIt);
for (IntersectionIterator isIt = problem().gridView().template 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);
//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
// the last condition is needed to properly assemble in the presence
// of ghost elements
if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce)
&& (globalIdxI > globalIdxJ) && haveSameLevel
&& elementNeighbor->partitionType() == Dune::InteriorEntity)
continue;
entries = 0;
//check for hanging nodes
if(!haveSameLevel && enableMPFA)
{
problem().pressureModel().getMpfaFlux(isIt, cellDataI);
}
else
{
problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);
int globalIdxJ = problem().variables().index(*elementNeighbor);
//set right hand side
this->f_[globalIdxI] -= entries[rhs];
//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;
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
entries = 0;
//check for hanging nodes
if(!haveSameLevel && enableMPFA)
{
problem().pressureModel().getMpfaFlux(isIt, cellDataI);
}
// set off-diagonal entry
this->A_[globalIdxI][globalIdxJ] -= entries[matrix];
// The second condition is needed to not spoil the ghost element entries
if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce)
&& elementNeighbor->partitionType() == Dune::InteriorEntity)
{
this->f_[globalIdxJ] += entries[rhs];
this->A_[globalIdxJ][globalIdxJ] += entries[matrix];
this->A_[globalIdxJ][globalIdxI] -= entries[matrix];
}
}
} // end neighbor
/************* boundary face ************************/
else
{
problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);
entries = 0;
problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);
//set right hand side
this->f_[globalIdxI] -= entries[rhs];
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
// set off-diagonal entry
this->A_[globalIdxI][globalIdxJ] -= entries[matrix];
if (GET_PROP_VALUE(TypeTag, VisitFacesOnlyOnce))
{
this->f_[globalIdxJ] += entries[rhs];
this->A_[globalIdxJ][globalIdxJ] += entries[matrix];
this->A_[globalIdxJ][globalIdxI] -= entries[matrix];
}
}
} //end interfaces loop
// printmatrix(std::cout, this->A_, "global stiffness matrix", "row", 11, 3);
/***** storage term ***********/
entries = 0;
problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
}
// assemble overlap and ghost element contributions
else
{
this->A_[globalIdxI] = 0.0;
this->A_[globalIdxI][globalIdxI] = 1.0;
this->f_[globalIdxI] = this->pressure()[globalIdxI];
}
} // end neighbor
/************* boundary face ************************/
else
{
entries = 0;
problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);
//set right hand side
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
}
} //end interfaces loop
// printmatrix(std::cout, this->A_, "global stiffness matrix", "row", 11, 3);
/***** storage term ***********/
entries = 0;
problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
} // end grid traversal
// printmatrix(std::cout, this->A_, "global stiffness matrix after assempling", "row", 11,3);
// printvector(std::cout, this->f_, "right hand side", "row", 10);
......
dumux/decoupled/2p2c/fvpressure2p2cmultiphysics.hh
View file @ ffaf6bd9
......@@ -229,70 +229,83 @@ void FVPressure2P2CMultiPhysics<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);
Dune::FieldVector<Scalar, 2> entries(0.);
// assemble interior element contributions
if (eIt->partitionType() == Dune::InteriorEntity)
{
// get the cell data
CellData& cellDataI = problem().variables().cellData(globalIdxI);
/***** source term ***********/
if(cellDataI.subdomain() != 2)
problem().pressureModel().get1pSource(entries,*eIt, cellDataI);
else
problem().pressureModel().getSource(entries,*eIt, cellDataI, first);
Dune::FieldVector<Scalar, 2> entries(0.);
this->f_[globalIdxI] = entries[rhs];
/***** source term ***********/
if(cellDataI.subdomain() != 2)
problem().pressureModel().get1pSource(entries,*eIt, cellDataI);
else
problem().pressureModel().getSource(entries,*eIt, cellDataI, first);
/***** 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())
this->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)
{
int globalIdxJ = problem().variables().index(*(isIt->outside()));
/************* handle interior face *****************/
if (isIt->neighbor())
{
int globalIdxJ = problem().variables().index(*(isIt->outside()));
if (cellDataI.subdomain() != 2
or problem().variables().cellData(globalIdxJ).subdomain() != 2) // cell in the 1p domain
get1pFlux(entries, *isIt, cellDataI);
else
problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);
if (cellDataI.subdomain() != 2
or problem().variables().cellData(globalIdxJ).subdomain() != 2) // cell in the 1p domain
get1pFlux(entries, *isIt, cellDataI);
else
problem().pressureModel().getFlux(entries, *isIt, cellDataI, first);
//set right hand side
this->f_[globalIdxI] -= entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
// set off-diagonal entry
this->A_[globalIdxI][globalIdxJ] = -entries[matrix];
} // end neighbor
//set right hand side
this->f_[globalIdxI] -= entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
// set off-diagonal entry
this->A_[globalIdxI][globalIdxJ] = -entries[matrix];
} // end neighbor
/************* boundary face ************************/
else
{
if (cellDataI.subdomain() != 2) //the current cell in the 1p domain
get1pFluxOnBoundary(entries, *isIt, cellDataI);
/************* boundary face ************************/
else
problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);
//set right hand side
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
}
} //end interfaces loop
// printmatrix(std::cout, this->A_, "global stiffness matrix", "row", 11, 3);
{
if (cellDataI.subdomain() != 2) //the current cell in the 1p domain
get1pFluxOnBoundary(entries, *isIt, cellDataI);
else
problem().pressureModel().getFluxOnBoundary(entries, *isIt, cellDataI, first);
//set right hand side
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
}
} //end interfaces loop
// printmatrix(std::cout, this->A_, "global stiffness matrix", "row", 11, 3);
/***** storage term ***********/
if (cellDataI.subdomain() != 2) //the current cell in the 1p domain
get1pStorage(entries, *eIt, cellDataI);
else
problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
/***** storage term ***********/
if (cellDataI.subdomain() != 2) //the current cell in the 1p domain
get1pStorage(entries, *eIt, cellDataI);
else
problem().pressureModel().getStorage(entries, *eIt, cellDataI, first);
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
this->f_[globalIdxI] += entries[rhs];
// set diagonal entry
this->A_[globalIdxI][globalIdxI] += entries[matrix];
}
// assemble overlap and ghost element contributions
else
{
this->A_[globalIdxI] = 0.0;
this->A_[globalIdxI][globalIdxI] = 1.0;
this->f_[globalIdxI] = this->pressure()[globalIdxI];
}
} // end grid traversal
// printmatrix(std::cout, this->A_, "global stiffness matrix after assempling", "row", 11,3);
// printvector(std::cout, this->f_, "right hand side", "row", 10);
......
dumux/decoupled/2p2c/fvtransport2p2cadaptive.hh
View file @ ffaf6bd9
......@@ -25,6 +25,7 @@
#include <dune/grid/common/gridenums.hh>
#include <dumux/decoupled/2p2c/fvtransport2p2c.hh>
#include <dumux/common/math.hh>
#include <dumux/linear/vectorexchange.hh>
/**
* @file
......@@ -251,6 +252,22 @@ void FVTransport2P2CAdaptive<TypeTag>::update(const Scalar t, Scalar& dt, Transp
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 = "
......
dumux/decoupled/2p2c/fvtransport2p2cmultiphysics.hh
View file @ ffaf6bd9
......@@ -20,6 +20,7 @@
#define DUMUX_FVTRANSPORT2P2C_MULTIPHYSICS_HH
#include <dumux/decoupled/2p2c/fvtransport2p2c.hh>
#include <dumux/linear/vectorexchange.hh>
/**
* @file
......@@ -197,6 +198,22 @@ void FVTransport2P2CMultiPhysics<TypeTag>::update(const Scalar t, Scalar& dt, Tr
}
}
} // 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 = "<<dt * GET_PARAM_FROM_GROUP(TypeTag, Scalar, Impet, CFLFactor)<< std::endl;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment