Merge pull request from blattms/parallel-tracers

[feature] Adds support for computing tracers in parallel runs.
This commit is contained in:
Markus Blatt 2021-11-23 12:39:48 +01:00 committed by GitHub
commit 6d8223f233
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 262 additions and 54 deletions

View File

@ -303,6 +303,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/utils/ParallelEclipseState.hpp
opm/simulators/utils/ParallelRestart.hpp
opm/simulators/utils/PropsCentroidsDataHandle.hpp
opm/simulators/utils/VectorVectorDataHandle.hpp
opm/simulators/wells/PerfData.hpp
opm/simulators/wells/PerforationData.hpp
opm/simulators/wells/RateConverter.hpp

View File

@ -24,6 +24,8 @@
#include <config.h>
#include <ebos/eclgenerictracermodel.hh>
#include <opm/simulators/linalg/PropertyTree.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/polyhedralgrid.hh>
@ -31,10 +33,13 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TracerVdTable.hpp>
#include <opm/simulators/linalg/FlexibleSolver.hpp>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/schwarz.hh>
#if HAVE_DUNE_FEM
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@ -45,19 +50,58 @@
#include <iostream>
#include <set>
#include <stdexcept>
#include <functional>
#include <array>
#include <string>
namespace Opm {
#if HAVE_MPI
template<class M, class V>
struct TracerSolverSelector
{
using Comm = Dune::OwnerOverlapCopyCommunication<int, int>;
using TracerOperator = Dune::OverlappingSchwarzOperator<M, V, V, Comm>;
using type = Dune::FlexibleSolver<M, V>;
};
template<class Vector, class Grid, class Matrix>
std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
Dune::OwnerOverlapCopyCommunication<int,int>>>,
std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
createParallelFlexibleSolver(const Grid&, const Matrix&, const PropertyTree&)
{
OPM_THROW(std::logic_error, "Grid not supported for parallel Tracers.");
return {nullptr, nullptr};
}
template<class Vector, class Matrix>
std::tuple<std::unique_ptr<Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
Dune::OwnerOverlapCopyCommunication<int,int>>>,
std::unique_ptr<typename TracerSolverSelector<Matrix,Vector>::type>>
createParallelFlexibleSolver(const Dune::CpGrid& grid, const Matrix& M, const PropertyTree& prm)
{
using TracerOperator = Dune::OverlappingSchwarzOperator<Matrix,Vector,Vector,
Dune::OwnerOverlapCopyCommunication<int,int>>;
using TracerSolver = Dune::FlexibleSolver<Matrix, Vector>;
const auto& cellComm = grid.cellCommunication();
auto op = std::make_unique<TracerOperator>(M, cellComm);
auto dummyWeights = [](){ return Vector();};
return {std::move(op), std::make_unique<TracerSolver>(*op, cellComm, prm, dummyWeights, 0)};
}
#endif
template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
EclGenericTracerModel(const GridView& gridView,
const EclipseState& eclState,
const CartesianIndexMapper& cartMapper,
const DofMapper& dofMapper)
const DofMapper& dofMapper,
const std::function<std::array<double,dimWorld>(int)> centroids)
: gridView_(gridView)
, eclState_(eclState)
, cartMapper_(cartMapper)
, dofMapper_(dofMapper)
, centroids_(centroids)
{
}
@ -87,7 +131,6 @@ doInit(bool enabled, size_t numGridDof,
size_t gasPhaseIdx, size_t oilPhaseIdx, size_t waterPhaseIdx)
{
const auto& tracers = eclState_.tracer();
const auto& comm = gridView_.comm();
if (tracers.size() == 0)
return; // tracer treatment is supposed to be disabled
@ -102,14 +145,6 @@ doInit(bool enabled, size_t numGridDof,
return; // Tracer transport must be enabled by the user
}
if (comm.size() > 1) {
tracerNames_.resize(0);
if (comm.rank() == 0)
std::cout << "Warning: The tracer model currently does not work for parallel runs\n"
<< std::flush;
return;
}
// retrieve the number of tracers from the deck
const size_t numTracers = tracers.size();
tracerNames_.resize(numTracers);
@ -145,12 +180,10 @@ doInit(bool enabled, size_t numGridDof,
}
//TVDPF keyword
else {
const auto& eclGrid = eclState_.getInputGrid();
for (size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx){
int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
const auto& center = eclGrid.getCellCenter(cartDofIdx);
tracerConcentration_[tracerIdx][globalDofIdx] = tracer.free_tvdp.evaluate("TRACER_CONCENTRATION", center[2]);
tracerConcentration_[tracerIdx][globalDofIdx] =
tracer.free_tvdp.evaluate("TRACER_CONCENTRATION",
centroids_(globalDofIdx)[2]);
}
}
++tracerIdx;
@ -223,24 +256,50 @@ linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
int maxIter = 100;
int verbosity = 0;
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
PropertyTree prm;
prm.put("maxiter", maxIter);
prm.put("tol", tolerance);
prm.put("verbosity", verbosity);
prm.put("solver", std::string("bicgstab"));
prm.put("preconditioner.type", std::string("ParOverILU0"));
TracerOperator tracerOperator(M);
TracerScalarProduct tracerScalarProduct;
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
#if HAVE_MPI
if(gridView_.grid().comm().size() > 1)
{
auto [tracerOperator, solver] =
createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
(void) tracerOperator;
TracerSolver solver (tracerOperator, tracerScalarProduct,
tracerPreconditioner, tolerance, maxIter,
verbosity);
Dune::InverseOperatorResult result;
solver->apply(x, b, result);
Dune::InverseOperatorResult result;
solver.apply(x, b, result);
// return the result of the solver
return result.converged;
}
else
{
#endif
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
// return the result of the solver
return result.converged;
TracerOperator tracerOperator(M);
TracerScalarProduct tracerScalarProduct;
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
TracerSolver solver (tracerOperator, tracerScalarProduct,
tracerPreconditioner, tolerance, maxIter,
verbosity);
Dune::InverseOperatorResult result;
solver.apply(x, b, result);
// return the result of the solver
return result.converged;
#if HAVE_MPI
}
#endif
}
template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
@ -255,29 +314,57 @@ linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::
int maxIter = 100;
int verbosity = 0;
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
PropertyTree prm;
prm.put("maxiter", maxIter);
prm.put("tol", tolerance);
prm.put("verbosity", verbosity);
prm.put("solver", std::string("bicgstab"));
prm.put("preconditioner.type", std::string("ParOverILU0"));
TracerOperator tracerOperator(M);
TracerScalarProduct tracerScalarProduct;
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
TracerSolver solver (tracerOperator, tracerScalarProduct,
tracerPreconditioner, tolerance, maxIter,
verbosity);
bool converged = true;
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
x[nrhs] = 0.0;
Dune::InverseOperatorResult result;
solver.apply(x[nrhs], b[nrhs], result);
converged = (converged && result.converged);
#if HAVE_MPI
if(gridView_.grid().comm().size() > 1)
{
auto [tracerOperator, solver] =
createParallelFlexibleSolver<TracerVector>(gridView_.grid(), M, prm);
(void) tracerOperator;
bool converged = true;
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
x[nrhs] = 0.0;
Dune::InverseOperatorResult result;
solver->apply(x[nrhs], b[nrhs], result);
converged = (converged && result.converged);
}
return converged;
}
else
{
#endif
using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
// return the result of the solver
return converged;
TracerOperator tracerOperator(M);
TracerScalarProduct tracerScalarProduct;
TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
TracerSolver solver (tracerOperator, tracerScalarProduct,
tracerPreconditioner, tolerance, maxIter,
verbosity);
bool converged = true;
for (size_t nrhs =0; nrhs < b.size(); ++nrhs) {
x[nrhs] = 0.0;
Dune::InverseOperatorResult result;
solver.apply(x[nrhs], b[nrhs], result);
converged = (converged && result.converged);
}
// return the result of the solver
return converged;
#if HAVE_MPI
}
#endif
}
#if HAVE_DUNE_FEM

View File

@ -51,7 +51,7 @@ public:
using TracerMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
static const int dimWorld = Grid::dimensionworld;
/*!
* \brief Return the number of tracers considered by the tracerModel.
*/
@ -78,7 +78,8 @@ protected:
EclGenericTracerModel(const GridView& gridView,
const EclipseState& eclState,
const CartesianIndexMapper& cartMapper,
const DofMapper& dofMapper);
const DofMapper& dofMapper,
const std::function<std::array<double,dimWorld>(int)> centroids);
/*!
* \brief Initialize all internal data structures needed by the tracer module
@ -109,7 +110,8 @@ protected:
// <wellName, tracerIdx> -> wellRate
std::map<std::pair<std::string, std::string>, double> wellTracerRate_;
/// \brief Function returning the cell centers
std::function<std::array<double,dimWorld>(int)> centroids_;
};
} // namespace Opm

View File

@ -31,6 +31,7 @@
#include <ebos/eclgenerictracermodel.hh>
#include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/utils/VectorVectorDataHandle.hpp>
#include <string>
#include <vector>
@ -90,7 +91,8 @@ public:
: BaseType(simulator.vanguard().gridView(),
simulator.vanguard().eclState(),
simulator.vanguard().cartesianIndexMapper(),
simulator.model().dofMapper())
simulator.model().dofMapper(),
simulator.vanguard().cellCentroids())
, simulator_(simulator)
, wat_(TracerBatch<TracerVector>(waterPhaseIdx))
, oil_(TracerBatch<TracerVector>(oilPhaseIdx))
@ -234,6 +236,15 @@ protected:
for (; elemIt != elemEndIt; ++ elemIt) {
elemCtx.updateAll(*elemIt);
size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/0);
if (elemIt->partitionType() != Dune::InteriorEntity)
{
// Dirichlet boundary conditions needed for the parallel matrix
(*this->tracerMatrix_)[I][I][0][0] = 1.;
continue;
}
Scalar extrusionFactor =
elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
Valgrind::CheckDefined(extrusionFactor);
@ -244,7 +255,6 @@ protected:
* extrusionFactor;
Scalar dt = elemCtx.simulator().timeStepSize();
size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/0);
size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/1);
std::vector<Scalar> storageOfTimeIndex1(tr.numTracer());
@ -322,6 +332,12 @@ protected:
}
}
}
// Communicate overlap using grid Communication
auto handle = VectorVectorDataHandle<GridView, std::vector<TracerVector>>(tr.residual_,
simulator_.gridView());
simulator_.gridView().communicate(handle, Dune::InteriorBorder_All_Interface,
Dune::ForwardCommunication);
}
template <class TrRe>

View File

@ -0,0 +1,102 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=4:
/*
Copyright 2021 Equinor AS.
This file is part of the Open Porous Media project (OPM).
OPM 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 3 of the License, or
(at your option) any later version.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/**
* \file
* \brief A datahandle sending data located in multiple vectors
* \author Markus Blatt, OPM-OP AS
*/
namespace Opm
{
/// \brief A data handle sending multiple data store in vectors attached
/// to cells.
///
/// Each data is assumed to a container with operator[] and the
/// class operates on a vector of these.
/// \tparam GridView the type of the grid view the data associated with
/// \tparam The type of the vector of vectors.
template<class GridView, class Vector>
class VectorVectorDataHandle
: public Dune::CommDataHandleIF<VectorVectorDataHandle<GridView,Vector>,
std::decay_t<decltype(Vector()[0][0])>>
{
public:
/// \brief the data type we send
using DataType = std::decay_t<decltype(Vector()[0][0])>;
/// \brief Constructor
/// \param data The vector of data vectors
/// \param gridView The gridview the data is attached to.
VectorVectorDataHandle(Vector& data, const GridView& gridView)
: data_(data), gridView_(gridView)
{}
bool contains(int /* dim */, int codim) const
{
return codim == 0;
}
bool fixedsize(int /* dim */, int /* codim */) const
{
return true;
}
bool fixedSize(int /* dim */, int /* codim */) const
{
return true;
}
template<class EntityType>
std::size_t size(const EntityType /* entity */) const
{
return data_.size();
}
template<class BufferType, class EntityType>
void gather(BufferType& buffer, const EntityType& e) const
{
for(const auto& vec: data_)
{
buffer.write(vec[gridView_.indexSet().index(e)]);
}
}
template<class BufferType, class EntityType>
void scatter(BufferType& buffer, const EntityType& e, std::size_t n)
{
assert(n == data_.size());
for(auto& vec: data_)
{
buffer.read(vec[gridView_.indexSet().index(e)]);
}
}
private:
Vector& data_;
const GridView& gridView_;
};
} // end namespace Opm