Merge pull request #5687 from GitPaean/tyring_ecl_compositional

ecl style output for compositional simulation
This commit is contained in:
Bård Skaflestad 2024-10-30 15:49:35 +01:00 committed by GitHub
commit 9e7c3b9254
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
18 changed files with 806 additions and 180 deletions

View File

@ -829,6 +829,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/flow/NewTranFluxModule.hpp
opm/simulators/flow/NonlinearSolver.hpp
opm/simulators/flow/OutputBlackoilModule.hpp
opm/simulators/flow/OutputCompositionalModule.hpp
opm/simulators/flow/partitionCells.hpp
opm/simulators/flow/PolyhedralGridVanguard.hpp
opm/simulators/flow/priVarsPacking.hpp

View File

@ -498,6 +498,26 @@ private:
fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
}
// determine the component total fractions
// TODO: duplicated code here, while should be refactored out when we swithing
// to starting from total mole fractions
Dune::FieldVector<Scalar, numComponents> z(0.0);
Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto saturation = getValue(fs.saturation(phaseIdx));
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar tmp = getValue(fs.molarity(phaseIdx, compIdx)) * saturation;
tmp = max(tmp, 1.e-8);
z[compIdx] += tmp;
sumMoles += tmp;
}
}
z /= sumMoles;
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
fs.setMoleFraction(compIdx, z[compIdx]);
}
// Set initial K and L
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation Ktmp = fs.wilsonK_(compIdx);

View File

@ -0,0 +1,97 @@
/*
Copyright 2024, SINTEF Digital
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/>.
*/
// this is an empty model that having a lot of empty interfaces.
// it is use for the development when some facility class are not ready
#include <opm/output/data/Aquifer.hpp>
#include <opm/output/data/Wells.hpp>
#include <opm/models/discretization/common/baseauxiliarymodule.hh>
#include <opm/input/eclipse/Schedule/Well/WellTestState.hpp>
namespace Opm {
template<typename TypeTag>
class EmptyModel : public BaseAuxiliaryModule<TypeTag>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
public:
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
EmptyModel(Simulator& /*simulator*/)
{
}
void init(){}
template<class Something>
void init(Something /*A*/){}
void prepareTracerBatches(){};
using NeighborSet = std::set<unsigned>;
void linearize(SparseMatrixAdapter& /*matrix*/, GlobalEqVector& /*residual*/){};
unsigned numDofs() const{return 0;};
void addNeighbors(std::vector<NeighborSet>& /*neighbors*/) const{};
//void applyInitial(){};
void initialSolutionApplied(){};
//void initFromRestart(const data::Aquifers& aquiferSoln);
template <class Restarter>
void serialize(Restarter& /*res*/){};
template <class Restarter>
void deserialize(Restarter& /*res*/){};
void beginEpisode(){};
void beginTimeStep(){};
void beginIteration(){};
// add the water rate due to aquifers to the source term.
template<class RateVector, class Context>
void addToSource(RateVector& /*rates*/, const Context& /*context*/,
unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const {}
template<class RateVector>
void addToSource(RateVector& /*rates*/, unsigned /*globalSpaceIdx*/,
unsigned /*timeIdx*/) const {}
void endIteration()const{};
void endTimeStep(){};
void endEpisode(){};
void applyInitial(){};
template<class RateType>
void computeTotalRatesForDof(RateType& /*rate*/, unsigned /*globalIdx*/) const{};
auto wellData() const {
return data::Wells{};
}
auto wellBlockAveragePressures() const {
return data::WellBlockAveragePressures{};
}
auto groupAndNetworkData(const int&) const {
return data::GroupAndNetworkValues{};
}
auto wellTestState() const {
return WellTestState{};
}
auto aquiferData() const {
return data::Aquifers{};
}
};
} // end of namespace Opm

View File

@ -30,6 +30,8 @@
#include <opm/simulators/linalg/parallelbicgstabbackend.hh>
#include <flowexperimental/comp/EmptyModel.hpp>
// // the current code use eclnewtonmethod adding other conditions to proceed_ should do the trick for KA
// // adding linearshe sould be chaning the update_ function in the same class with condition that the error is reduced.
// the trick is to be able to recalculate the residual from here.
@ -37,55 +39,6 @@
// suggestTimeStep is taken from newton solver in problem.limitTimestep
namespace Opm {
template<typename TypeTag>
class EmptyModel : public BaseAuxiliaryModule<TypeTag>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
public:
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
EmptyModel(Simulator& /*simulator*/)
{
}
void init(){}
template<class Something>
void init(Something /*A*/){}
void prepareTracerBatches(){};
using NeighborSet = std::set<unsigned>;
void linearize(SparseMatrixAdapter& /*matrix*/, GlobalEqVector& /*residual*/){};
unsigned numDofs() const{return 0;};
void addNeighbors(std::vector<NeighborSet>& /*neighbors*/) const{};
//void applyInitial(){};
void initialSolutionApplied(){};
//void initFromRestart(const data::Aquifers& aquiferSoln);
template <class Restarter>
void serialize(Restarter& /*res*/){};
template <class Restarter>
void deserialize(Restarter& /*res*/){};
void beginEpisode(){};
void beginTimeStep(){};
void beginIteration(){};
// add the water rate due to aquifers to the source term.
template<class RateVector, class Context>
void addToSource(RateVector& /*rates*/, const Context& /*context*/,
unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const {}
template<class RateVector>
void addToSource(RateVector& /*rates*/, unsigned /*globalSpaceIdx*/,
unsigned /*timeIdx*/) const {}
void endIteration()const{};
void endTimeStep(){};
void endEpisode(){};
void applyInitial(){};
template<class RateType>
void computeTotalRatesForDof(RateType& /*rate*/, unsigned /*globalIdx*/) const{};
};
template<int numComp>
int dispatchFlowExpComp(int argc, char** argv);

View File

@ -113,6 +113,8 @@ public:
const Scalar flashTolerance = Parameters::Get<Parameters::FlashTolerance<Scalar>>();
const int flashVerbosity = Parameters::Get<Parameters::FlashVerbosity>();
const std::string flashTwoPhaseMethod = Parameters::Get<Parameters::FlashTwoPhaseMethod>();
// TODO: the formulation here is still to begin with XMF and YMF values to derive ZMF value
// TODO: we should check how we update ZMF in the newton update, since it is the primary variables.
// extract the total molar densities of the components
ComponentVector z(0.);
@ -132,6 +134,10 @@ public:
z /= sumz;
}
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
fluidState_.setMoleFraction(compIdx, z[compIdx]);
}
Evaluation p = priVars.makeEvaluation(pressure0Idx, timeIdx);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fluidState_.setPressure(phaseIdx, p);

View File

@ -120,22 +120,10 @@ public:
// the energy module
EnergyModule::setPriVarTemperatures(*this, fluidState);
// determine the component fractions
Dune::FieldVector<Scalar, numComponents> z(0.0);
Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar tmp = Opm::getValue(fluidState.molarity(phaseIdx, compIdx) * fluidState.saturation(phaseIdx));
z[compIdx] += Opm::max(tmp, 1e-8);
sumMoles += tmp;
}
}
z /= sumMoles;
for (int i = 0; i < numComponents - 1; ++i)
(*this)[z0Idx + i] = z[i];
(*this)[z0Idx + i] = getValue(fluidState.moleFraction(i));
(*this)[pressure0Idx] = Opm::getValue(fluidState.pressure(0));
(*this)[pressure0Idx] = getValue(fluidState.pressure(0));
}
/*!

View File

@ -669,7 +669,7 @@ public:
// write initial condition
if (problem_->shouldWriteOutput())
EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->writeOutput());
EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->writeOutput(true));
timeStepSize_ = oldTimeStepSize;
timeStepIdx_ = oldTimeStepIdx;
@ -748,7 +748,7 @@ public:
// write the result to disk
writeTimer_.start();
if (problem_->shouldWriteOutput())
EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->writeOutput());
EWOMS_CATCH_PARALLEL_EXCEPTIONS_FATAL(problem_->writeOutput(true));
writeTimer_.stop();
// do the next time integration

View File

@ -30,23 +30,30 @@
#include <dune/grid/common/partitionset.hh>
#include <opm/common/TimingMacros.hpp> // OPM_TIMEBLOCK
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/input/eclipse/Schedule/RPTConfig.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/output/eclipse/Inplace.hpp>
#include <opm/output/eclipse/RestartValue.hpp>
#include <opm/models/blackoil/blackoilproperties.hh> // Properties::EnableMech, EnableTemperature, EnableSolvent
#include <opm/models/common/multiphasebaseproperties.hh> // Properties::FluidSystem
#include <opm/simulators/flow/CollectDataOnIORank.hpp>
#include <opm/simulators/flow/countGlobalCells.hpp>
#include <opm/simulators/flow/EclGenericWriter.hpp>
#include <opm/simulators/flow/FlowBaseVanguard.hpp>
#include <opm/simulators/flow/OutputBlackoilModule.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/utils/ParallelRestart.hpp>
#include <opm/simulators/utils/ParallelSerialization.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <limits>
#include <map>
#include <memory>
@ -101,7 +108,7 @@ namespace Opm {
* - This class requires to use the black oil model with the element
* centered finite volume discretization.
*/
template <class TypeTag>
template <class TypeTag, class OutputModule>
class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>,
GetPropType<TypeTag, Properties::EquilGrid>,
GetPropType<TypeTag, Properties::GridView>,
@ -120,7 +127,7 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
using ElementIterator = typename GridView::template Codim<0>::Iterator;
using BaseType = EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>;
typedef Dune::MultipleCodimMultipleGeomTypeMapper< GridView > VertexMapper;
enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
@ -132,7 +139,7 @@ public:
static void registerParameters()
{
OutputBlackOilModule<TypeTag>::registerParameters();
OutputModule::registerParameters();
Parameters::Register<Parameters::EnableAsyncEclOutput>
("Write the ECL-formated results in a non-blocking way "
@ -169,13 +176,13 @@ public:
eclBroadcast(this->simulator_.vanguard().grid().comm(), smryCfg);
this->outputModule_ = std::make_unique<OutputBlackOilModule<TypeTag>>
this->outputModule_ = std::make_unique<OutputModule>
(simulator, smryCfg, this->collectOnIORank_);
}
else
#endif
{
this->outputModule_ = std::make_unique<OutputBlackOilModule<TypeTag>>
this->outputModule_ = std::make_unique<OutputModule>
(simulator, this->eclIO_->finalSummaryConfig(), this->collectOnIORank_);
}
@ -428,7 +435,6 @@ public:
const bool isFloresn = this->outputModule_->hasFloresn();
auto floresn = this->outputModule_->getFloresn();
// data::Solution localCellData = {};
if (! isSubStep || Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
if (localCellData.empty()) {
@ -661,10 +667,10 @@ public:
}
}
const OutputBlackOilModule<TypeTag>& outputModule() const
const OutputModule& outputModule() const
{ return *outputModule_; }
OutputBlackOilModule<TypeTag>& mutableOutputModule() const
OutputModule& mutableOutputModule() const
{ return *outputModule_; }
Scalar restartTimeStepSize() const
@ -826,7 +832,7 @@ private:
}
Simulator& simulator_;
std::unique_ptr<OutputBlackOilModule<TypeTag> > outputModule_;
std::unique_ptr<OutputModule> outputModule_;
Scalar restartTimeStepSize_;
int rank_ ;
Inplace inplace_;

View File

@ -493,7 +493,7 @@ public:
* \brief Write the requested quantities of the current solution into the output
* files.
*/
void writeOutput(bool verbose = true)
virtual void writeOutput(bool verbose)
{
OPM_TIMEBLOCK(problemWriteOutput);
// use the generic code to prepare the output fields and to

View File

@ -50,6 +50,7 @@
#include <opm/simulators/flow/FlowProblemBlackoilProperties.hpp>
#include <opm/simulators/flow/FlowThresholdPressure.hpp>
#include <opm/simulators/flow/MixingRateControls.hpp>
#include <opm/simulators/flow/OutputBlackoilModule.hpp>
#include <opm/simulators/flow/VtkTracerModule.hpp>
#include <opm/simulators/utils/satfunc/SatfuncConsistencyCheckManager.hpp>
@ -141,7 +142,7 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
using ModuleParams = typename BlackOilLocalResidualTPFA<TypeTag>::ModuleParams;
using InitialFluidState = typename EquilInitializer<TypeTag>::ScalarFluidState;
using EclWriterType = EclWriter<TypeTag>;
using EclWriterType = EclWriter<TypeTag, OutputBlackOilModule<TypeTag> >;
#if HAVE_DAMARIS
using DamarisWriterType = DamarisWriter<TypeTag>;
#endif
@ -432,7 +433,7 @@ public:
/*!
* \brief Called by the simulator after each time integration.
*/
void endTimeStep () override
void endTimeStep() override
{
FlowProblemType::endTimeStep();
@ -517,7 +518,7 @@ public:
* \brief Write the requested quantities of the current solution into the output
* files.
*/
void writeOutput(bool verbose = true)
void writeOutput(bool verbose) override
{
FlowProblemType::writeOutput(verbose);

View File

@ -32,6 +32,9 @@
#include <opm/simulators/flow/FlowProblem.hpp>
#include <opm/simulators/flow/FlowThresholdPressure.hpp>
#include <opm/simulators/flow/OutputCompositionalModule.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/thermal/EclThermalLawManager.hpp>
@ -82,6 +85,7 @@ class FlowProblemComp : public FlowProblem<TypeTag>
using typename FlowProblemType::RateVector;
using InitialFluidState = CompositionalFluidState<Scalar, FluidSystem>;
using EclWriterType = EclWriter<TypeTag, OutputCompositionalModule<TypeTag> >;
public:
using FlowProblemType::porosity;
@ -94,6 +98,8 @@ public:
{
FlowProblemType::registerParameters();
EclWriterType::registerParameters();
// tighter tolerance is needed for compositional modeling here
Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1e-7);
}
@ -104,7 +110,10 @@ public:
*/
explicit FlowProblemComp(Simulator& simulator)
: FlowProblemType(simulator)
, thresholdPressures_(simulator)
{
eclWriter_ = std::make_unique<EclWriterType>(simulator);
enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
}
/*!
@ -127,9 +136,18 @@ public:
[&vg = this->simulator().vanguard()](const unsigned int it) { return vg.gridIdxToEquilGridIdx(it); });
updated = true;
};
// TODO: we might need to do the same with FlowProblemBlackoil for parallel
finishTransmissibilities();
if (enableEclOutput_) {
eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
std::function<unsigned int(unsigned int)> equilGridToGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridEquilIdxToGridIdx(i);
};
eclWriter_->extractOutputTransAndNNC(equilGridToGrid);
}
const auto& eclState = simulator.vanguard().eclState();
const auto& schedule = simulator.vanguard().schedule();
@ -179,6 +197,11 @@ public:
FlowProblemType::readMaterialParameters_();
FlowProblemType::readThermalParameters_();
// write the static output files (EGRID, INIT)
if (enableEclOutput_) {
eclWriter_->writeInit();
}
const auto& initconfig = eclState.getInitConfig();
if (initconfig.restartRequested())
readEclRestartSolution_();
@ -220,6 +243,53 @@ public:
}
}
/*!
* \brief Called by the simulator after each time integration.
*/
void endTimeStep() override
{
FlowProblemType::endTimeStep();
const bool isSubStep = !this->simulator().episodeWillBeOver();
// after the solution is updated, the values in output module also needs to be updated
this->eclWriter_->mutableOutputModule().invalidateLocalData();
// For CpGrid with LGRs, ecl/vtk output is not supported yet.
const auto& grid = this->simulator().vanguard().gridView().grid();
using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
if (!isCpGrid || (grid.maxLevel() == 0)) {
this->eclWriter_->evalSummaryState(isSubStep);
}
}
void writeReports(const SimulatorTimer& timer) {
if (enableEclOutput_){
eclWriter_->writeReports(timer);
}
}
/*!
* \brief Write the requested quantities of the current solution into the output
* files.
*/
void writeOutput(bool verbose) override
{
FlowProblemType::writeOutput(verbose);
const bool isSubStep = !this->simulator().episodeWillBeOver();
data::Solution localCellData = {};
if (enableEclOutput_) {
if (Parameters::Get<Parameters::EnableWriteAllSolutions>() || !isSubStep) {
eclWriter_->writeOutput(std::move(localCellData), isSubStep);
}
}
}
/*!
* \copydoc FvBaseProblem::boundary
*
@ -253,25 +323,26 @@ public:
{
const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
const auto& initial_fs = initialFluidStates_[globalDofIdx];
Opm::CompositionalFluidState<Evaluation, FluidSystem> fs;
using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
// TODO: the current approach is assuming we begin with XMF and YMF.
// TODO: maybe we should make it begin with ZMF
using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
for (unsigned p = 0; p < numPhases; ++p) { // TODO: assuming the phaseidx continuous
ComponentVector evals;
auto& last_eval = evals[numComponents - 1];
ComponentVector vals;
auto& last_eval = vals[numComponents - 1];
last_eval = 1.;
for (unsigned c = 0; c < numComponents - 1; ++c) {
const auto val = initial_fs.moleFraction(p, c);
const Evaluation eval = Evaluation::createVariable(val, c + 1);
evals[c] = eval;
last_eval -= eval;
vals[c] = val;
last_eval -= val;
}
for (unsigned c = 0; c < numComponents; ++c) {
fs.setMoleFraction(p, c, evals[c]);
fs.setMoleFraction(p, c, vals[c]);
}
// pressure
const auto p_val = initial_fs.pressure(p);
fs.setPressure(p, Evaluation::createVariable(p_val, 0));
fs.setPressure(p, p_val);
const auto sat_val = initial_fs.saturation(p);
fs.setSaturation(p, sat_val);
@ -281,7 +352,7 @@ public:
}
{
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
typename FluidSystem::template ParameterCache<Scalar> paramCache;
paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
@ -290,13 +361,31 @@ public:
fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
}
// determine the component fractions
Dune::FieldVector<Scalar, numComponents> z(0.0);
Scalar sumMoles = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto saturation = getValue(fs.saturation(phaseIdx));
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
Scalar tmp = getValue(fs.molarity(phaseIdx, compIdx)) * saturation;
tmp = max(tmp, 1e-8);
z[compIdx] += tmp;
sumMoles += tmp;
}
}
z /= sumMoles;
// Set initial K and L
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const Evaluation Ktmp = fs.wilsonK_(compIdx);
const auto& Ktmp = fs.wilsonK_(compIdx);
fs.setKvalue(compIdx, Ktmp);
}
const Evaluation& Ltmp = -1.0;
for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
fs.setMoleFraction(compIdx, z[compIdx]);
}
const Scalar& Ltmp = -1.0;
fs.setLvalue(Ltmp);
values.assignNaive(fs);
@ -316,11 +405,19 @@ public:
const std::vector<InitialFluidState>& initialFluidStates() const
{ return initialFluidStates_; }
const FlowThresholdPressure<TypeTag>& thresholdPressure() const
{
assert( !thresholdPressures_.enableThresholdPressure() &&
" Threshold Pressures are not supported by compostional simulation ");
return thresholdPressures_;
}
// TODO: do we need this one?
template<class Serializer>
void serializeOp(Serializer& serializer)
{
serializer(static_cast<FlowProblemType&>(*this));
serializer(*eclWriter_);
}
protected:
@ -467,6 +564,8 @@ protected:
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
const std::size_t data_idx = compIdx * numDof + dofIdx;
const Scalar zmf = zmfData[data_idx];
dofFluidState.setMoleFraction(compIdx, zmf);
if (gas_active) {
const auto ymf = (dofFluidState.saturation(FluidSystem::gasPhaseIdx) > 0.) ? zmf : Scalar{0};
dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
@ -492,9 +591,14 @@ private:
throw std::logic_error("polymer is disabled for compositional modeling and you're trying to add polymer to BC");
}
FlowThresholdPressure<TypeTag> thresholdPressures_;
std::vector<InitialFluidState> initialFluidStates_;
bool enableVtkOutput_;
bool enableEclOutput_{false};
std::unique_ptr<EclWriterType> eclWriter_;
bool enableVtkOutput_{false};
};
} // namespace Opm

View File

@ -28,6 +28,8 @@
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
#include <opm/grid/common/CommunicationUtils.hpp>
#include <opm/output/data/Solution.hpp>
@ -201,7 +203,8 @@ GenericOutputBlackoilModule(const EclipseState& eclState,
bool enableBrine,
bool enableSaltPrecipitation,
bool enableExtbo,
bool enableMICP)
bool enableMICP,
bool isCompositional)
: eclState_(eclState)
, schedule_(schedule)
, summaryState_(summaryState)
@ -220,6 +223,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState,
, enableSaltPrecipitation_(enableSaltPrecipitation)
, enableExtbo_(enableExtbo)
, enableMICP_(enableMICP)
, isCompositional_(isCompositional)
, local_data_valid_(false)
{
const auto& fp = eclState_.fieldProps();
@ -497,83 +501,90 @@ assignToSolution(data::Solution& sol)
target);
};
const auto baseSolutionArrays = std::array {
DataEntry{"1OVERBG", UnitSystem::measure::gas_inverse_formation_volume_factor, invB_[gasPhaseIdx]},
DataEntry{"1OVERBO", UnitSystem::measure::oil_inverse_formation_volume_factor, invB_[oilPhaseIdx]},
DataEntry{"1OVERBW", UnitSystem::measure::water_inverse_formation_volume_factor, invB_[waterPhaseIdx]},
DataEntry{"FLRGASI+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XPlus)][gasCompIdx]},
DataEntry{"FLRGASJ+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YPlus)][gasCompIdx]},
DataEntry{"FLRGASK+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][gasCompIdx]},
DataEntry{"FLROILI+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XPlus)][oilCompIdx]},
DataEntry{"FLROILJ+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YPlus)][oilCompIdx]},
DataEntry{"FLROILK+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][oilCompIdx]},
DataEntry{"FLRWATI+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XPlus)][waterCompIdx]},
DataEntry{"FLRWATJ+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YPlus)][waterCompIdx]},
DataEntry{"FLRWATK+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][waterCompIdx]},
DataEntry{"FOAM", UnitSystem::measure::identity, cFoam_},
DataEntry{"GASKR", UnitSystem::measure::identity, relativePermeability_[gasPhaseIdx]},
DataEntry{"GAS_DEN", UnitSystem::measure::density, density_[gasPhaseIdx]},
DataEntry{"GAS_VISC", UnitSystem::measure::viscosity, viscosity_[gasPhaseIdx]},
DataEntry{"OILKR", UnitSystem::measure::identity, relativePermeability_[oilPhaseIdx]},
DataEntry{"OIL_DEN", UnitSystem::measure::density, density_[oilPhaseIdx]},
DataEntry{"OIL_VISC", UnitSystem::measure::viscosity, viscosity_[oilPhaseIdx]},
DataEntry{"PBUB", UnitSystem::measure::pressure, bubblePointPressure_},
DataEntry{"PCGW", UnitSystem::measure::pressure, pcgw_},
DataEntry{"PCOG", UnitSystem::measure::pressure, pcog_},
DataEntry{"PCOW", UnitSystem::measure::pressure, pcow_},
DataEntry{"PDEW", UnitSystem::measure::pressure, dewPointPressure_},
DataEntry{"POLYMER", UnitSystem::measure::identity, cPolymer_},
DataEntry{"PPCW", UnitSystem::measure::pressure, ppcw_},
DataEntry{"PRESROCC", UnitSystem::measure::pressure, minimumOilPressure_},
DataEntry{"PRESSURE", UnitSystem::measure::pressure, fluidPressure_},
DataEntry{"RPORV", UnitSystem::measure::volume, rPorV_},
DataEntry{"RS", UnitSystem::measure::gas_oil_ratio, rs_},
DataEntry{"RSSAT", UnitSystem::measure::gas_oil_ratio, gasDissolutionFactor_},
DataEntry{"RV", UnitSystem::measure::oil_gas_ratio, rv_},
DataEntry{"RVSAT", UnitSystem::measure::oil_gas_ratio, oilVaporizationFactor_},
DataEntry{"SALT", UnitSystem::measure::salinity, cSalt_},
DataEntry{"SGMAX", UnitSystem::measure::identity, sgmax_},
DataEntry{"SHMAX", UnitSystem::measure::identity, shmax_},
DataEntry{"SOMAX", UnitSystem::measure::identity, soMax_},
DataEntry{"SOMIN", UnitSystem::measure::identity, somin_},
DataEntry{"SSOLVENT", UnitSystem::measure::identity, sSol_},
DataEntry{"SWHY1", UnitSystem::measure::identity, swmin_},
DataEntry{"SWMAX", UnitSystem::measure::identity, swMax_},
DataEntry{"WATKR", UnitSystem::measure::identity, relativePermeability_[waterPhaseIdx]},
DataEntry{"WAT_DEN", UnitSystem::measure::density, density_[waterPhaseIdx]},
DataEntry{"WAT_VISC", UnitSystem::measure::viscosity, viscosity_[waterPhaseIdx]},
// if index not specified, we treat it as valid (>= 0)
auto addEntry = [](std::vector<DataEntry>& container, const std::string& name, UnitSystem::measure measure, const auto& flowArray, int index = 1) {
if (index >= 0) { // Only add if index is valid
container.emplace_back(name, measure, flowArray);
}
};
std::vector<DataEntry> baseSolutionVector;
addEntry(baseSolutionVector, "1OVERBG", UnitSystem::measure::gas_inverse_formation_volume_factor, invB_[gasPhaseIdx], gasPhaseIdx);
addEntry(baseSolutionVector, "1OVERBO", UnitSystem::measure::oil_inverse_formation_volume_factor, invB_[oilPhaseIdx], oilPhaseIdx);
addEntry(baseSolutionVector, "1OVERBW", UnitSystem::measure::water_inverse_formation_volume_factor, invB_[waterPhaseIdx], waterPhaseIdx);
addEntry(baseSolutionVector, "FLRGASI+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XPlus)][gasCompIdx], gasCompIdx);
addEntry(baseSolutionVector, "FLRGASJ+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YPlus)][gasCompIdx], gasCompIdx);
addEntry(baseSolutionVector, "FLRGASK+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][gasCompIdx], gasCompIdx);
addEntry(baseSolutionVector, "FLROILI+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XPlus)][oilCompIdx], oilCompIdx);
addEntry(baseSolutionVector, "FLROILJ+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YPlus)][oilCompIdx], oilCompIdx);
addEntry(baseSolutionVector, "FLROILK+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][oilCompIdx], oilCompIdx);
addEntry(baseSolutionVector, "FLRWATI+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XPlus)][waterCompIdx], waterCompIdx);
addEntry(baseSolutionVector, "FLRWATJ+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YPlus)][waterCompIdx], waterCompIdx);
addEntry(baseSolutionVector, "FLRWATK+", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][waterCompIdx], waterCompIdx);
addEntry(baseSolutionVector, "FOAM", UnitSystem::measure::identity, cFoam_);
addEntry(baseSolutionVector, "GASKR", UnitSystem::measure::identity, relativePermeability_[gasPhaseIdx], gasPhaseIdx);
addEntry(baseSolutionVector, "GAS_DEN", UnitSystem::measure::density, density_[gasPhaseIdx], gasPhaseIdx);
addEntry(baseSolutionVector, "GAS_VISC", UnitSystem::measure::viscosity, viscosity_[gasPhaseIdx], gasPhaseIdx);
addEntry(baseSolutionVector, "OILKR", UnitSystem::measure::identity, relativePermeability_[oilPhaseIdx], oilPhaseIdx);
addEntry(baseSolutionVector, "OIL_DEN", UnitSystem::measure::density, density_[oilPhaseIdx], oilPhaseIdx);
addEntry(baseSolutionVector, "OIL_VISC", UnitSystem::measure::viscosity, viscosity_[oilPhaseIdx], oilPhaseIdx);
addEntry(baseSolutionVector, "PBUB", UnitSystem::measure::pressure, bubblePointPressure_);
addEntry(baseSolutionVector, "PCGW", UnitSystem::measure::pressure, pcgw_);
addEntry(baseSolutionVector, "PCOG", UnitSystem::measure::pressure, pcog_);
addEntry(baseSolutionVector, "PCOW", UnitSystem::measure::pressure, pcow_);
addEntry(baseSolutionVector, "PDEW", UnitSystem::measure::pressure, dewPointPressure_);
addEntry(baseSolutionVector, "POLYMER", UnitSystem::measure::identity, cPolymer_);
addEntry(baseSolutionVector, "PPCW", UnitSystem::measure::pressure, ppcw_);
addEntry(baseSolutionVector, "PRESROCC", UnitSystem::measure::pressure, minimumOilPressure_);
addEntry(baseSolutionVector, "PRESSURE", UnitSystem::measure::pressure, fluidPressure_);
addEntry(baseSolutionVector, "RPORV", UnitSystem::measure::volume, rPorV_);
addEntry(baseSolutionVector, "RS", UnitSystem::measure::gas_oil_ratio, rs_);
addEntry(baseSolutionVector, "RSSAT", UnitSystem::measure::gas_oil_ratio, gasDissolutionFactor_);
addEntry(baseSolutionVector, "RV", UnitSystem::measure::oil_gas_ratio, rv_);
addEntry(baseSolutionVector, "RVSAT", UnitSystem::measure::oil_gas_ratio, oilVaporizationFactor_);
addEntry(baseSolutionVector, "SALT", UnitSystem::measure::salinity, cSalt_);
addEntry(baseSolutionVector, "SGMAX", UnitSystem::measure::identity, sgmax_);
addEntry(baseSolutionVector, "SHMAX", UnitSystem::measure::identity, shmax_);
addEntry(baseSolutionVector, "SOMAX", UnitSystem::measure::identity, soMax_);
addEntry(baseSolutionVector, "SOMIN", UnitSystem::measure::identity, somin_);
addEntry(baseSolutionVector, "SSOLVENT", UnitSystem::measure::identity, sSol_);
addEntry(baseSolutionVector, "SWHY1", UnitSystem::measure::identity, swmin_);
addEntry(baseSolutionVector, "SWMAX", UnitSystem::measure::identity, swMax_);
addEntry(baseSolutionVector, "WATKR", UnitSystem::measure::identity, relativePermeability_[waterPhaseIdx], waterPhaseIdx);
addEntry(baseSolutionVector, "WAT_DEN", UnitSystem::measure::density, density_[waterPhaseIdx], waterPhaseIdx);
addEntry(baseSolutionVector, "WAT_VISC", UnitSystem::measure::viscosity, viscosity_[waterPhaseIdx], waterPhaseIdx);
// Separate these as flows*_ may be defined due to BFLOW[I|J|K] even without FLOWS in RPTRST
const auto flowsSolutionArrays = std::array {
DataEntry{"FLOGASI+", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XPlus)][gasCompIdx]},
DataEntry{"FLOGASJ+", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YPlus)][gasCompIdx]},
DataEntry{"FLOGASK+", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][gasCompIdx]},
DataEntry{"FLOOILI+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XPlus)][oilCompIdx]},
DataEntry{"FLOOILJ+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YPlus)][oilCompIdx]},
DataEntry{"FLOOILK+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][oilCompIdx]},
DataEntry{"FLOWATI+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XPlus)][waterCompIdx]},
DataEntry{"FLOWATJ+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YPlus)][waterCompIdx]},
DataEntry{"FLOWATK+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][waterCompIdx]},
DataEntry{"FLOGASI-", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XMinus)][gasCompIdx]},
DataEntry{"FLOGASJ-", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YMinus)][gasCompIdx]},
DataEntry{"FLOGASK-", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][gasCompIdx]},
DataEntry{"FLOOILI-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XMinus)][oilCompIdx]},
DataEntry{"FLOOILJ-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YMinus)][oilCompIdx]},
DataEntry{"FLOOILK-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][oilCompIdx]},
DataEntry{"FLOWATI-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XMinus)][waterCompIdx]},
DataEntry{"FLOWATJ-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YMinus)][waterCompIdx]},
DataEntry{"FLOWATK-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][waterCompIdx]},
DataEntry{"FLRGASI-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XMinus)][gasCompIdx]},
DataEntry{"FLRGASJ-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YMinus)][gasCompIdx]},
DataEntry{"FLRGASK-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][gasCompIdx]},
DataEntry{"FLROILI-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XMinus)][oilCompIdx]},
DataEntry{"FLROILJ-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YMinus)][oilCompIdx]},
DataEntry{"FLROILK-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][oilCompIdx]},
DataEntry{"FLRWATI-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XMinus)][waterCompIdx]},
DataEntry{"FLRWATJ-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YMinus)][waterCompIdx]},
DataEntry{"FLRWATK-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][waterCompIdx]},
};
std::vector<DataEntry> flowsSolutionVector;
addEntry(flowsSolutionVector, "FLOGASI+", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XPlus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLOGASJ+", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YPlus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLOGASK+", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLOOILI+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XPlus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLOOILJ+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YPlus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLOOILK+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLOWATI+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XPlus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLOWATJ+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YPlus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLOWATK+", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZPlus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLOGASI-", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XMinus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLOGASJ-", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YMinus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLOGASK-", UnitSystem::measure::gas_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLOOILI-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XMinus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLOOILJ-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YMinus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLOOILK-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLOWATI-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::XMinus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLOWATJ-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::YMinus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLOWATK-", UnitSystem::measure::liquid_surface_rate, flows_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLRGASI-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XMinus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLRGASJ-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YMinus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLRGASK-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][gasCompIdx], gasCompIdx);
addEntry(flowsSolutionVector, "FLROILI-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XMinus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLROILJ-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YMinus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLROILK-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][oilCompIdx], oilCompIdx);
addEntry(flowsSolutionVector, "FLRWATI-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::XMinus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLRWATJ-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::YMinus)][waterCompIdx], waterCompIdx);
addEntry(flowsSolutionVector, "FLRWATK-", UnitSystem::measure::rate, flores_[FaceDir::ToIntersectionIndex(Dir::ZMinus)][waterCompIdx], waterCompIdx);
const auto extendedSolutionArrays = std::array {
DataEntry{"BIOFILM", UnitSystem::measure::identity, cBiofilm_},
@ -624,12 +635,44 @@ assignToSolution(data::Solution& sol)
DataEntry{"UREA", UnitSystem::measure::density, cUrea_},
};
for (const auto& array : baseSolutionArrays) {
// basically, for compositional, we can not use std::array for this. We need to generate the ZMF1, ZMF2, and so on
// and also, we need to map these values.
// TODO: the following should go to a function
if (this->isCompositional_) {
auto compositionalEntries = std::vector<DataEntry>{};
{
// ZMF
for (int i = 0; i < numComponents; ++i) {
const auto name = fmt::format("ZMF{}", i + 1); // Generate ZMF1, ZMF2, ...
compositionalEntries.emplace_back(name, UnitSystem::measure::identity, moleFractions_[i]);
}
// XMF
for (int i = 0; i < numComponents; ++i) {
const auto name = fmt::format("XMF{}", i + 1); // Generate XMF1, XMF2, ...
compositionalEntries.emplace_back(name, UnitSystem::measure::identity,
phaseMoleFractions_[oilPhaseIdx][i]);
}
// YMF
for (int i = 0; i < numComponents; ++i) {
const auto name = fmt::format("YMF{}", i + 1); // Generate YMF1, YMF2, ...
compositionalEntries.emplace_back(name, UnitSystem::measure::identity,
phaseMoleFractions_[gasPhaseIdx][i]);
}
}
for (const auto& array: compositionalEntries) {
doInsert(array, data::TargetType::RESTART_SOLUTION);
}
}
for (const auto& array : baseSolutionVector) {
doInsert(array, data::TargetType::RESTART_SOLUTION);
}
if (this->enableFlows_) {
for (const auto& array : flowsSolutionArrays) {
for (const auto& array : flowsSolutionVector) {
doInsert(array, data::TargetType::RESTART_SOLUTION);
}
}
@ -660,6 +703,15 @@ assignToSolution(data::Solution& sol)
data::TargetType::RESTART_SOLUTION);
}
if (this->isCompositional_ && FluidSystem::phaseIsActive(oilPhaseIdx) &&
! this->saturation_[oilPhaseIdx].empty())
{
sol.insert("SOIL", UnitSystem::measure::identity,
std::move(this->saturation_[oilPhaseIdx]),
data::TargetType::RESTART_SOLUTION);
}
if ((eclState_.runspec().co2Storage() || eclState_.runspec().h2Storage()) && !rsw_.empty()) {
auto mfrac = std::vector<double>(this->rsw_.size(), 0.0);
@ -1481,6 +1533,30 @@ doAllocBuffers(const unsigned bufferSize,
overburdenPressure_.resize(bufferSize, 0.0);
}
if (this->isCompositional_) {
if (rstKeywords["ZMF"] > 0) {
rstKeywords["ZMF"] = 0;
for (int i = 0; i < numComponents; ++i) {
moleFractions_[i].resize(bufferSize, 0.0);
}
}
if (rstKeywords["XMF"] > 0 && FluidSystem::phaseIsActive(oilPhaseIdx)) {
rstKeywords["XMF"] = 0;
for (int i = 0; i < numComponents; ++i) {
phaseMoleFractions_[oilPhaseIdx][i].resize(bufferSize, 0.0);
}
}
if (rstKeywords["YMF"] > 0 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
rstKeywords["YMF"] = 0;
for (int i = 0; i < numComponents; ++i) {
phaseMoleFractions_[gasPhaseIdx][i].resize(bufferSize, 0.0);
}
}
}
//Warn for any unhandled keyword
if (log) {
for (auto& keyValue: rstKeywords) {
@ -1752,4 +1828,15 @@ INSTANTIATE_TYPE(double)
INSTANTIATE_TYPE(float)
#endif
#define INSTANTIATE_COMP(NUM) \
template<class T> using FS##NUM = GenericOilGasFluidSystem<T, NUM>; \
template class GenericOutputBlackoilModule<FS##NUM<double>>;
INSTANTIATE_COMP(2)
INSTANTIATE_COMP(3)
INSTANTIATE_COMP(4)
INSTANTIATE_COMP(5)
INSTANTIATE_COMP(6)
INSTANTIATE_COMP(7)
} // namespace Opm

View File

@ -335,20 +335,21 @@ protected:
bool enableBrine,
bool enableSaltPrecipitation,
bool enableExtbo,
bool enableMICP);
bool enableMICP,
bool isCompositional = false);
void doAllocBuffers(unsigned bufferSize,
unsigned reportStepNum,
const bool substep,
const bool log,
const bool isRestart,
const bool vapparsActive,
const bool enablePCHysteresis,
const bool enableNonWettingHysteresis,
const bool enableWettingHysteresis,
unsigned numTracers,
const std::vector<bool>& enableSolTracers,
unsigned numOutputNnc);
const bool vapparsActive = false,
const bool enablePCHysteresis = false,
const bool enableNonWettingHysteresis =false,
const bool enableWettingHysteresis = false,
unsigned numTracers = 0,
const std::vector<bool>& enableSolTracers = {},
unsigned numOutputNnc = 0);
void makeRegionSum(Inplace& inplace,
const std::string& region_name,
@ -405,6 +406,7 @@ protected:
bool enableSaltPrecipitation_{false};
bool enableExtbo_{false};
bool enableMICP_{false};
bool isCompositional_{false};
bool forceDisableFipOutput_{false};
bool forceDisableFipresvOutput_{false};
@ -539,6 +541,10 @@ protected:
std::array<ScalarBuffer, numPhases> viscosity_;
std::array<ScalarBuffer, numPhases> relativePermeability_;
// total mole fractions for each component
std::array<ScalarBuffer, numComponents> moleFractions_;
// mole fractions for each component in each phase
std::array<std::array<ScalarBuffer, numComponents>, numPhases> phaseMoleFractions_;
std::vector<ScalarBuffer> freeTracerConcentrations_;
std::vector<ScalarBuffer> solTracerConcentrations_;

View File

@ -81,6 +81,8 @@ public:
//! \details Returns the union of explicitly configured entries and defaulted values.
std::vector<Scalar> getRestartVector() const;
bool enableThresholdPressure() const;
protected:
/*!
* \brief Actually compute the threshold pressures over a face as a pre-compute step.

View File

@ -246,6 +246,15 @@ getRestartVector() const
return result;
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
bool
GenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
enableThresholdPressure() const
{
return this->enableThresholdPressure_;
}
template<class Grid, class GridView, class ElementMapper, class Scalar>
void
GenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::

View File

@ -0,0 +1,346 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
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 2 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
* \copydoc Opm::OutputCompositionalModule
*/
#ifndef OPM_OUTPUT_COMPOSITIONAL_MODULE_HPP
#define OPM_OUTPUT_COMPOSITIONAL_MODULE_HPP
#include <dune/grid/common/gridenums.hh>
#include <opm/simulators/utils/moduleVersion.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/TimingMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/models/utils/parametersystem.hpp>
#include <opm/models/utils/propertysystem.hh>
#include <opm/simulators/flow/FlowBaseVanguard.hpp>
#include <opm/simulators/flow/GenericOutputBlackoilModule.hpp>
#include <algorithm>
#include <cstddef>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>
namespace Opm
{
// forward declaration
template <class TypeTag>
class EcfvDiscretization;
/*!
* \ingroup BlackOilSimulator
*
* \brief Output module for the results black oil model writing in
* ECL binary format.
*/
template <class TypeTag>
class OutputCompositionalModule : public GenericOutputBlackoilModule<GetPropType<TypeTag, Properties::FluidSystem>>
{
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Discretization = GetPropType<TypeTag, Properties::Discretization>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using BaseType = GenericOutputBlackoilModule<FluidSystem>;
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
public:
template <class CollectDataToIORankType>
OutputCompositionalModule(const Simulator& simulator,
const SummaryConfig& smryCfg,
const CollectDataToIORankType& collectToIORank)
: BaseType(simulator.vanguard().eclState(),
simulator.vanguard().schedule(),
smryCfg,
simulator.vanguard().summaryState(),
moduleVersionName(),
getPropValue<TypeTag, Properties::EnableEnergy>(),
getPropValue<TypeTag, Properties::EnableTemperature>(),
getPropValue<TypeTag, Properties::EnableMech>(),
getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
getPropValue<TypeTag, Properties::EnableSaltPrecipitation>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnableMICP>(),
true)
, simulator_(simulator)
{
for (auto& region_pair : this->regions_) {
this->createLocalRegion_(region_pair.second);
}
auto isCartIdxOnThisRank = [&collectToIORank](const int idx) {
return collectToIORank.isCartIdxOnThisRank(idx);
};
this->setupBlockData(isCartIdxOnThisRank);
if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
const std::string msg = "The output code does not support --owner-cells-first=false.";
if (collectToIORank.isIORank()) {
OpmLog::error(msg);
}
OPM_THROW_NOLOG(std::runtime_error, msg);
}
if (smryCfg.match("[FB]PP[OGW]") || smryCfg.match("RPP[OGW]*")) {
auto rset = this->eclState_.fieldProps().fip_regions();
rset.push_back("PVTNUM");
// Note: We explicitly use decltype(auto) here because the
// default scheme (-> auto) will deduce an undesirable type. We
// need the "reference to vector" semantics in this instance.
this->regionAvgDensity_
.emplace(this->simulator_.gridView().comm(),
FluidSystem::numPhases, rset,
[fp = std::cref(this->eclState_.fieldProps())]
(const std::string& rsetName) -> decltype(auto)
{ return fp.get().get_int(rsetName); });
}
}
/*!
* \brief Register all run-time parameters for the Vtk output module.
*/
static void registerParameters()
{
}
/*!
* \brief Allocate memory for the scalar fields we would like to
* write to ECL output files
*/
void
allocBuffers(const unsigned bufferSize,
const unsigned reportStepNum,
const bool substep,
const bool log,
const bool isRestart)
{
if (! std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
return;
}
this->doAllocBuffers(bufferSize,
reportStepNum,
substep,
log,
isRestart);
}
/*!
* \brief Modify the internal buffers according to the intensive
* quanties relevant for an element
*/
void processElement(const ElementContext& elemCtx)
{
OPM_TIMEBLOCK_LOCAL(processElement);
if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
return;
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
// const unsigned pvtRegionIdx = 0; // elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (this->saturation_[phaseIdx].empty())
continue;
this->saturation_[phaseIdx][globalDofIdx] = getValue(fs.saturation(phaseIdx));
Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]);
}
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
if (this->moleFractions_[compIdx].empty()) continue;
this->moleFractions_[compIdx][globalDofIdx] = getValue(fs.moleFraction(compIdx));
}
// XMF and YMF
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
if (this->phaseMoleFractions_[oilPhaseIdx][compIdx].empty()) continue;
this->phaseMoleFractions_[oilPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(oilPhaseIdx, compIdx));
}
if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
if (this->phaseMoleFractions_[gasPhaseIdx][compIdx].empty()) continue;
this->phaseMoleFractions_[gasPhaseIdx][compIdx][globalDofIdx] = getValue(fs.moleFraction(gasPhaseIdx, compIdx));
}
}
if (!this->fluidPressure_.empty()) {
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
// Output oil pressure as default
this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx));
} else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
// Output gas if oil is not present
this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx));
} else {
// Output water if neither oil nor gas is present
this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx));
}
Valgrind::CheckDefined(this->fluidPressure_[globalDofIdx]);
}
if (!this->temperature_.empty()) {
this->temperature_[globalDofIdx] = getValue(fs.temperature(oilPhaseIdx));
Valgrind::CheckDefined(this->temperature_[globalDofIdx]);
}
}
}
void processElementFlows(const ElementContext& /* elemCtx */)
{
OPM_TIMEBLOCK_LOCAL(processElementBlockData);
if (!std::is_same_v<Discretization, EcfvDiscretization<TypeTag>>)
return;
}
void processElementBlockData(const ElementContext& /* elemCtx */)
{
OPM_TIMEBLOCK_LOCAL(processElementBlockData);
if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
return;
}
/*!
* \brief Capture connection fluxes, particularly to account for inter-region flows.
*
* \tparam ActiveIndex Callable type, typically a lambda, that enables
* retrieving the active index, on the local MPI rank, of a
* particular cell/element. Must support a function call operator of
* the form
\code
int operator()(const Element& elem) const
\endcode
*
* \tparam CartesianIndex Callable type, typically a lambda, that
* enables retrieving the globally unique Cartesian index of a
* particular cell/element given its active index on the local MPI
* rank. Must support a function call operator of the form
\code
int operator()(const int activeIndex) const
\endcode
*
* \param[in] elemCtx Primary lookup structure for per-cell/element
* dynamic information.
*
* \param[in] activeIndex Mapping from cell/elements to linear indices
* on local MPI rank.
*
* \param[in] cartesianIndex Mapping from active index on local MPI rank
* to globally unique Cartesian cell/element index.
*/
template <class ActiveIndex, class CartesianIndex>
void processFluxes(const ElementContext& /* elemCtx */,
ActiveIndex&& /* activeIndex*/,
CartesianIndex&& /* cartesianIndex */)
{
}
/*!
* \brief Prepare for capturing connection fluxes, particularly to
* account for inter-region flows.
*/
void initializeFluxData()
{
// Inter-region flow rates. Note: ".clear()" prepares to accumulate
// contributions per bulk connection between FIP regions.
this->interRegionFlows_.clear();
}
/*!
* \brief Finalize capturing connection fluxes.
*/
void finalizeFluxData()
{
this->interRegionFlows_.compress();
}
/*!
* \brief Get read-only access to collection of inter-region flows.
*/
const InterRegFlowMap& getInterRegFlows() const
{
return this->interRegionFlows_;
}
void updateFluidInPlace(const unsigned /* globalDofIdx */,
const IntensiveQuantities& /* intQuants */,
const double /* totVolume */)
{
// this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
}
private:
bool isDefunctParallelWell(std::string wname) const override
{
if (simulator_.gridView().comm().size() == 1)
return false;
const auto& parallelWells = simulator_.vanguard().parallelWells();
std::pair<std::string, bool> value {wname, true};
auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
return candidate == parallelWells.end() || *candidate != value;
}
void createLocalRegion_(std::vector<int>& region)
{
std::size_t elemIdx = 0;
for (const auto& elem : elements(simulator_.gridView())) {
if (elem.partitionType() != Dune::InteriorEntity) {
region[elemIdx] = 0;
}
++elemIdx;
}
}
const Simulator& simulator_;
};
} // namespace Opm
#endif // OPM_OUTPUT_COMPOSITIONAL_MODULE_HPP

View File

@ -283,7 +283,7 @@ public:
simulator_.setEpisodeLength(0.0);
simulator_.setTimeStepSize(0.0);
wellModel_().beginReportStep(timer.currentStepNum());
simulator_.problem().writeOutput();
simulator_.problem().writeOutput(true);
report_.success.output_write_time += perfTimer.stop();
}
@ -370,7 +370,7 @@ public:
perfTimer.start();
const double nextstep = adaptiveTimeStepping_ ? adaptiveTimeStepping_->suggestedNextStep() : -1.0;
simulator_.problem().setNextTimeStepSize(nextstep);
simulator_.problem().writeOutput();
simulator_.problem().writeOutput(true);
report_.success.output_write_time += perfTimer.stop();
solver_->model().endReportStep();

View File

@ -353,7 +353,7 @@ void registerAdaptiveParameters();
time::StopWatch perfTimer;
perfTimer.start();
problem.writeOutput();
problem.writeOutput(true);
report.success.output_write_time += perfTimer.secsSinceStart();
}