OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

remove 'flow_solvent' and associated files

Browse Source 
use 'flow'
This commit is contained in:
Arne Morten Kvarving 2017-11-13 14:18:21 +01:00
parent ad23d98726
commit e0907ddc49
13 changed files with 0 additions and 3438 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
CMakeLists_files.cmake
View File

@ -40,7 +40,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/SimulatorIncompTwophaseAd.cpp
opm/autodiff/TransportSolverTwophaseAd.cpp
opm/autodiff/BlackoilPropsAdFromDeck.cpp
opm/autodiff/SolventPropsAdFromDeck.cpp
opm/autodiff/BlackoilModelParameters.cpp
opm/autodiff/WellDensitySegmented.cpp
opm/autodiff/LinearisedBlackoilResidual.cpp
@ -93,7 +92,6 @@ list (APPEND TEST_SOURCE_FILES
tests/test_welldensitysegmented.cpp
tests/test_vfpproperties.cpp
tests/test_singlecellsolves.cpp
tests/test_solventprops_ad.cpp
tests/test_multisegmentwells.cpp
tests/test_multiphaseupwind.cpp
tests/test_wellmodel.cpp
@ -121,7 +119,6 @@ list (APPEND EXAMPLE_SOURCE_FILES
examples/flow_reorder.cpp
examples/flow_sequential.cpp
examples/flow.cpp
examples/flow_solvent.cpp
examples/sim_2p_incomp.cpp
examples/sim_2p_incomp_ad.cpp
examples/sim_2p_comp_reorder.cpp
@ -142,7 +139,6 @@ list (APPEND PROGRAM_SOURCE_FILES
examples/flow_legacy.cpp
examples/flow_reorder.cpp
examples/flow_sequential.cpp
examples/flow_solvent.cpp
examples/opm_init_check.cpp
examples/sim_poly2p_comp_reorder.cpp
examples/sim_poly2p_incomp_reorder.cpp
@ -167,14 +163,11 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/BlackoilModelParameters.hpp
opm/autodiff/BlackoilPressureModel.hpp
opm/autodiff/BlackoilPropsAdFromDeck.hpp
opm/autodiff/SolventPropsAdFromDeck.hpp
opm/autodiff/Compat.hpp
opm/autodiff/CPRPreconditioner.hpp
opm/autodiff/createGlobalCellArray.hpp
opm/autodiff/DefaultBlackoilSolutionState.hpp
opm/autodiff/BlackoilSequentialModel.hpp
opm/autodiff/BlackoilSolventModel.hpp
opm/autodiff/BlackoilSolventModel_impl.hpp
opm/autodiff/BlackoilReorderingTransportModel.hpp
opm/autodiff/BlackoilTransportModel.hpp
opm/autodiff/fastSparseOperations.hpp
@ -185,7 +178,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/FlowMainEbos.hpp
opm/autodiff/FlowMainPolymer.hpp
opm/autodiff/FlowMainSequential.hpp
opm/autodiff/FlowMainSolvent.hpp
opm/autodiff/GeoProps.hpp
opm/autodiff/GridHelpers.hpp
opm/autodiff/GridInit.hpp
@ -212,14 +204,11 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/SimulatorBase_impl.hpp
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp
opm/autodiff/SimulatorIncompTwophaseAd.hpp
opm/autodiff/SimulatorSequentialBlackoil.hpp
opm/autodiff/TransportSolverTwophaseAd.hpp
opm/autodiff/WellDensitySegmented.hpp
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
opm/autodiff/VFPProperties.hpp
opm/autodiff/VFPHelpers.hpp
@ -241,8 +230,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/MSWellHelpers.hpp
opm/autodiff/BlackoilWellModel.hpp
opm/autodiff/BlackoilWellModel_impl.hpp
opm/autodiff/StandardWellsSolvent.hpp
opm/autodiff/StandardWellsSolvent_impl.hpp
opm/autodiff/MissingFeatures.hpp
opm/autodiff/ThreadHandle.hpp
opm/polymer/CompressibleTpfaPolymer.hpp

42
examples/flow_solvent.cpp
View File

@ -1,42 +0,0 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 IRIS 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/>.
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/grid.h>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp>
#include <opm/autodiff/FlowMainSolvent.hpp>
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
typedef UnstructuredGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoilSolvent<Grid> Simulator;
Opm::FlowMainSolvent<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}

269
opm/autodiff/BlackoilSolventModel.hpp
View File

@ -1,269 +0,0 @@
/*
Copyright 2015 IRIS 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/>.
*/
#ifndef OPM_BLACKOILSOLVENTMODEL_HEADER_INCLUDED
#define OPM_BLACKOILSOLVENTMODEL_HEADER_INCLUDED
#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
#include <opm/autodiff/StandardWellsSolvent.hpp>
namespace Opm {
/// A model implementation for three-phase black oil
/// with one extra component.
///
///
/// It uses automatic differentiation via the class AutoDiffBlock
/// to simplify assembly of the jacobian matrix.
template<class Grid>
class BlackoilSolventModel : public BlackoilModelBase<Grid, StandardWellsSolvent, BlackoilSolventModel<Grid> >
{
public:
// --------- Types and enums ---------
typedef BlackoilModelBase<Grid, StandardWellsSolvent, BlackoilSolventModel<Grid> > Base;
typedef typename Base::ReservoirState ReservoirState;
typedef typename Base::WellState WellState;
// The next line requires C++11 support available in g++ 4.7.
// friend Base;
friend class BlackoilModelBase<Grid, StandardWellsSolvent, BlackoilSolventModel<Grid> >;
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
/// \param[in] param parameters
/// \param[in] grid grid data structure
/// \param[in] fluid fluid properties
/// \param[in] geo rock properties
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] solvent_props solvent properties
/// \param[in] wells well structure
/// \param[in] linsolver linear solver
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
/// \param[in] has_solvent turn on solvent feature
/// \param[in] is_miscible turn on miscible feature
BlackoilSolventModel(const typename Base::ModelParameters& param,
const Grid& grid,
const BlackoilPropsAdFromDeck& fluid,
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const SolventPropsAdFromDeck& solvent_props,
const StandardWellsSolvent& well_model,
const NewtonIterationBlackoilInterface& linsolver,
std::shared_ptr< const EclipseState > eclState,
std::shared_ptr< const Schedule> schedule,
std::shared_ptr< const SummaryConfig> summary_config,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output,
const bool has_solvent,
const bool is_miscible);
/// Apply an update to the primary variables, chopped if appropriate.
/// \param[in] dx updates to apply to primary variables
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void updateState(const V& dx,
ReservoirState& reservoir_state,
WellState& well_state);
using Base::wellModel;
std::vector<std::vector<double> >
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum);
protected:
// --------- Types and enums ---------
typedef typename Base::SolutionState SolutionState;
typedef typename Base::DataBlock DataBlock;
enum { Solvent = CanonicalVariablePositions::Next };
// --------- Data members ---------
const bool has_solvent_;
const int solvent_pos_;
const SolventPropsAdFromDeck& solvent_props_;
const bool is_miscible_;
std::vector<ADB> mu_eff_;
std::vector<ADB> b_eff_;
// Need to declare Base members we want to use here.
using Base::grid_;
using Base::fluid_;
using Base::geo_;
using Base::rock_comp_props_;
using Base::linsolver_;
using Base::active_;
using Base::canph_;
using Base::cells_;
using Base::ops_;
using Base::has_disgas_;
using Base::has_vapoil_;
using Base::param_;
using Base::use_threshold_pressure_;
using Base::threshold_pressures_by_connection_;
using Base::sd_;
using Base::phaseCondition_;
using Base::residual_;
using Base::terminal_output_;
using Base::pvdt_;
// --------- Protected methods ---------
// Need to declare Base members we want to use here.
using Base::wells;
using Base::variableState;
using Base::computeGasPressure;
using Base::applyThresholdPressures;
using Base::fluidRsSat;
using Base::fluidRvSat;
using Base::poroMult;
using Base::transMult;
using Base::updatePrimalVariableFromState;
using Base::updatePhaseCondFromPrimalVariable;
using Base::dpMaxRel;
using Base::dsMax;
using Base::drMaxRel;
using Base::maxResidualAllowed;
// using Base::updateWellControls;
// using Base::computeWellConnectionPressures;
// using Base::addWellControlEq;
// using Base::computePropertiesForWellConnectionPressures;
std::vector<ADB>
computeRelPerm(const SolutionState& state) const;
ADB
fluidViscosity(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond) const;
ADB
fluidReciprocFVF(const int phase,
const ADB& p ,
const ADB& temp ,
const ADB& rs ,
const ADB& rv ,
const std::vector<PhasePresence>& cond) const;
ADB
fluidDensity(const int phase,
const ADB& b,
const ADB& rs,
const ADB& rv) const;
void
makeConstantState(SolutionState& state) const;
std::vector<V>
variableStateInitials(const ReservoirState& x,
const WellState& xw) const;
std::vector<int>
variableStateIndices() const;
SolutionState
variableStateExtractVars(const ReservoirState& x,
const std::vector<int>& indices,
std::vector<ADB>& vars) const;
void
computeAccum(const SolutionState& state,
const int aix );
void
assembleMassBalanceEq(const SolutionState& state);
void
addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw);
void updateEquationsScaling();
void
computeMassFlux(const int actph ,
const V& transi,
const ADB& kr ,
const ADB& mu ,
const ADB& rho ,
const ADB& p ,
const SolutionState& state );
const std::vector<PhasePresence>
phaseCondition() const {return this->phaseCondition_;}
// compute effective viscosities (mu_eff_) and effective b factors (b_eff_) using the ToddLongstaff model
void computeEffectiveProperties(const SolutionState& state);
// compute density and viscosity using the ToddLongstaff mixing model
void computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const ADB po, const Opm::PhaseUsage pu);
// compute phase pressures.
std::vector<ADB>
computePressures(const ADB& po,
const ADB& sw,
const ADB& so,
const ADB& sg,
const ADB& ss) const;
};
/// Need to include concentration in our state variables, otherwise all is as
/// the default blackoil model.
struct BlackoilSolventSolutionState : public DefaultBlackoilSolutionState
{
explicit BlackoilSolventSolutionState(const int np)
: DefaultBlackoilSolutionState(np),
solvent_saturation( ADB::null())
{
}
ADB solvent_saturation;
};
/// Providing types by template specialisation of ModelTraits for BlackoilSolventModel.
template <class Grid>
struct ModelTraits< BlackoilSolventModel<Grid> >
{
typedef BlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoilSolvent WellState;
typedef BlackoilModelParameters ModelParameters;
typedef BlackoilSolventSolutionState SolutionState;
};
} // namespace Opm
#include "BlackoilSolventModel_impl.hpp"
#endif // OPM_BLACKOILSOLVENTMODEL_HEADER_INCLUDED

1116
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

File diff suppressed because it is too large Load Diff

127
opm/autodiff/FlowMainSolvent.hpp
View File

@ -1,127 +0,0 @@
/*
Copyright 2015 IRIS AS
Copyright 2015 SINTEF ICT, Applied Mathematics.
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/>.
*/
#ifndef OPM_FLOWMAINSOLVENT_HEADER_INCLUDED
#define OPM_FLOWMAINSOLVENT_HEADER_INCLUDED
#include <opm/autodiff/FlowMain.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
namespace Opm
{
// The FlowMainSolvent class is for a black-oil simulator with solvent.
template <class Grid, class Simulator>
class FlowMainSolvent : public FlowMainBase<FlowMainSolvent<Grid, Simulator>, Grid, Simulator>
{
protected:
using Base = FlowMainBase<FlowMainSolvent<Grid, Simulator>, Grid, Simulator>;
friend Base;
// Set in setupGridAndProps()
std::unique_ptr<SolventPropsAdFromDeck> solvent_props_;
// ------------ Methods ------------
// Print startup message if on output rank.
void printStartupMessage()
{
if (Base::output_cout_) {
const std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow-Solvent (version " << version << ")"
<< std::string(18 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow-Solvent is a simulator for fully implicit three-phase, *\n";
std::cout << "* four-component (black-oil + solvent) flow, and is part of OPM. *\n";
std::cout << "* For more information see http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
}
// Set up grid and property objects, by calling base class
// version and then creating solvent property object.
void setupGridAndProps()
{
Base::setupGridAndProps();
const Grid& grid = Base::grid_init_->grid();
solvent_props_.reset(new SolventPropsAdFromDeck(*Base::deck_,
*Base::eclipse_state_,
UgGridHelpers::numCells(grid),
UgGridHelpers::globalCell(grid)));
}
// Create simulator instance.
// Writes to:
// simulator_
void createSimulator()
{
// Create the simulator instance.
Base::simulator_.reset(new Simulator(Base::param_,
Base::grid_init_->grid(),
*Base::geoprops_,
*Base::fluidprops_,
*solvent_props_,
Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr,
*Base::fis_solver_,
Base::gravity_.data(),
Base::deck_->hasKeyword("DISGAS"),
Base::deck_->hasKeyword("VAPOIL"),
Base::eclipse_state_,
Base::schedule_,
Base::summary_config_,
*Base::output_writer_,
Base::deck_,
Base::threshold_pressures_,
Base::deck_->hasKeyword("SOLVENT")));
}
void setupLinearSolver()
{
// require_full_sparsity_pattern as default for solvent runs
if (Base::deck_->hasKeyword("SOLVENT") && !Base::param_.has("require_full_sparsity_pattern") ) {
Base::param_.insertParameter("require_full_sparsity_pattern","true");
}
Base::setupLinearSolver();
}
};
} // namespace Opm
#endif // OPM_FLOWMAINSOLVENT_HEADER_INCLUDED

147
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp
View File

@ -1,147 +0,0 @@
/*
Copyright 2015 IRIS 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/>.
*/
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_HEADER_INCLUDED
#include <opm/autodiff/SimulatorBase.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/autodiff/BlackoilSolventModel.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/autodiff/NonlinearSolver.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
//#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
//#include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
//#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <algorithm>
#include <cstddef>
#include <cassert>
#include <functional>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm
{
template <class GridT>
class SimulatorFullyImplicitBlackoilSolvent;
class StandardWellsSolvent;
template<class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilSolvent<GridT> >
{
typedef WellStateFullyImplicitBlackoilSolvent WellState;
typedef BlackoilState ReservoirState;
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef BlackoilSolventModel<Grid> Model;
typedef NonlinearSolver<Model> Solver;
typedef StandardWellsSolvent WellModel;
};
/// Class collecting all necessary components for a blackoil simulation with polymer
/// injection.
template <class GridT>
class SimulatorFullyImplicitBlackoilSolvent
: public SimulatorBase<SimulatorFullyImplicitBlackoilSolvent<GridT> >
{
typedef SimulatorFullyImplicitBlackoilSolvent<GridT> ThisType;
typedef SimulatorBase<ThisType> BaseType;
typedef SimulatorTraits<ThisType> Traits;
typedef typename Traits::Solver Solver;
typedef typename Traits::WellModel WellModel;
public:
SimulatorFullyImplicitBlackoilSolvent(const ParameterGroup& param,
const GridT& grid,
DerivedGeology& geo,
BlackoilPropsAdFromDeck& props,
const SolventPropsAdFromDeck& solvent_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool disgas,
const bool vapoil,
std::shared_ptr<EclipseState> eclipse_state,
std::shared_ptr<Schedule> schedule,
std::shared_ptr<SummaryConfig> summary_config,
BlackoilOutputWriter& output_writer,
std::shared_ptr< Deck > deck,
const std::vector<double>& threshold_pressures_by_face,
const bool solvent);
std::unique_ptr<Solver> createSolver(const WellModel& well_model);
void handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& wells_manager,
typename BaseType::WellState& well_state,
const Wells* wells);
private:
bool has_solvent_;
std::shared_ptr< Deck > deck_;
SolventPropsAdFromDeck solvent_props_;
bool is_miscible_;
};
} // namespace Opm
#include "SimulatorFullyImplicitBlackoilSolvent_impl.hpp"
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_HEADER_INCLUDED

147
opm/autodiff/SimulatorFullyImplicitBlackoilSolvent_impl.hpp
View File

@ -1,147 +0,0 @@
/*
Copyright 2015 IRIS 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/>.
*/
#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_IMPL_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_IMPL_HEADER_INCLUDED
namespace Opm
{
template <class GridT>
SimulatorFullyImplicitBlackoilSolvent<GridT>::
SimulatorFullyImplicitBlackoilSolvent(const ParameterGroup& param,
const GridT& grid,
DerivedGeology& geo,
BlackoilPropsAdFromDeck& props,
const SolventPropsAdFromDeck& solvent_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
std::shared_ptr<EclipseState> eclipse_state,
std::shared_ptr<Schedule> schedule,
std::shared_ptr<SummaryConfig> summary_config,
BlackoilOutputWriter& output_writer,
std::shared_ptr< Deck > deck,
const std::vector<double>& threshold_pressures_by_face,
const bool has_solvent)
: BaseType(param,
grid,
geo,
props,
rock_comp_props,
linsolver,
gravity,
has_disgas,
has_vapoil,
eclipse_state,
schedule,
summary_config,
output_writer,
threshold_pressures_by_face,
// names of deactivated wells in parallel run
std::unordered_set<std::string>())
, has_solvent_(has_solvent)
, deck_(deck)
, solvent_props_(solvent_props)
, is_miscible_(false)
{
if(deck->hasKeyword("MISCIBLE")) {
is_miscible_ = true;
}
}
template <class GridT>
auto SimulatorFullyImplicitBlackoilSolvent<GridT>::
createSolver(const WellModel& well_model)
-> std::unique_ptr<Solver>
{
typedef typename Traits::Model Model;
auto model = std::unique_ptr<Model>(new Model(BaseType::model_param_,
BaseType::grid_,
BaseType::props_,
BaseType::geo_,
BaseType::rock_comp_props_,
solvent_props_,
well_model,
BaseType::solver_,
BaseType::eclipse_state_,
BaseType::schedule_,
BaseType::summary_config_,
BaseType::has_disgas_,
BaseType::has_vapoil_,
BaseType::terminal_output_,
has_solvent_,
is_miscible_));
if (!BaseType::threshold_pressures_by_face_.empty()) {
model->setThresholdPressures(BaseType::threshold_pressures_by_face_);
}
return std::unique_ptr<Solver>(new Solver(BaseType::solver_param_, std::move(model)));
}
template <class GridT>
void SimulatorFullyImplicitBlackoilSolvent<GridT>::
handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& /*wells_manager*/,
typename BaseType::WellState& well_state,
const Wells* wells)
{
// compute solvent inflow
const int nw = wells->number_of_wells;
std::vector<double> perfcells_fraction(wells->well_connpos[nw], 0.0);
size_t currentStep = timer.currentStepNum();
for (const auto& well_solvent : BaseType::schedule_->getWells( currentStep )) {
if (well_solvent->getStatus( currentStep ) == WellCommon::SHUT) {
continue;
}
WellInjectionProperties injection = well_solvent->getInjectionProperties(currentStep);
if (injection.injectorType == WellInjector::GAS) {
double solventFraction = well_solvent->getSolventFraction(currentStep);
// Find the solvent well in the well list and add properties to it
int wix = 0;
for (; wix < nw; ++wix) {
if (well_solvent->name() == wells->name[wix]) {
break;
}
}
if (wix == wells->number_of_wells) {
OPM_THROW(std::runtime_error, "Could not find a match for well "
<< well_solvent->name()
<< " from WSOLVENT.");
}
for (int j = wells->well_connpos[wix]; j < wells->well_connpos[wix+1]; ++j) {
perfcells_fraction[j] = solventFraction;
}
}
}
well_state.solventFraction() = perfcells_fraction;
}
} // namespace Opm
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILSOLVENT_IMPL_HEADER_INCLUDED

513
opm/autodiff/SolventPropsAdFromDeck.cpp
View File

@ -1,513 +0,0 @@
/*
Copyright 2015 IRIS
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/>.
*/
#include <config.h>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/utility/extractPvtTableIndex.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdsTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SsfnTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
namespace Opm
{
// Making these typedef to make the code more readable.
typedef SolventPropsAdFromDeck::ADB ADB;
typedef Eigen::SparseMatrix<double> S;
typedef Eigen::DiagonalMatrix<double, Eigen::Dynamic> D;
typedef SolventPropsAdFromDeck::V V;
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Block;
SolventPropsAdFromDeck::SolventPropsAdFromDeck(const Deck& deck,
const EclipseState& eclState,
const int number_of_cells,
const int* global_cell)
{
if (deck.hasKeyword("SOLVENT")) {
// retrieve the cell specific PVT table index from the deck
// and using the grid...
extractPvtTableIndex(cellPvtRegionIdx_, eclState, number_of_cells, global_cell);
extractTableIndex("SATNUM", eclState, number_of_cells, global_cell, cellSatNumRegionIdx_);
// surface densities
if (deck.hasKeyword("SDENSITY")) {
const auto& densityKeyword = deck.getKeyword("SDENSITY");
int numRegions = densityKeyword.size();
solvent_surface_densities_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
solvent_surface_densities_[regionIdx]
= densityKeyword.getRecord(regionIdx).getItem("SOLVENT_DENSITY").getSIDouble(0);
}
} else {
OPM_THROW(std::runtime_error, "SDENSITY must be specified in SOLVENT runs\n");
}
const auto& tables = eclState.getTableManager();
// pvt
const TableContainer& pvdsTables = tables.getPvdsTables();
if (!pvdsTables.empty()) {
int numRegions = pvdsTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
inverseBmu_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdsTable& pvdsTable = pvdsTables.getTable<PvdsTable>(regionIdx);
const auto& press = pvdsTable.getPressureColumn();
const auto& b = pvdsTable.getFormationFactorColumn();
const auto& visc = pvdsTable.getViscosityColumn();
const int sz = b.size();
std::vector<double> inverseBmu(sz);
std::vector<double> inverseB(sz);
for (int i = 0; i < sz; ++i) {
inverseB[i] = 1.0 / b[i];
inverseBmu[i] = 1.0 / (b[i] * visc[i]);
}
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
}
} else {
OPM_THROW(std::runtime_error, "PVDS must be specified in SOLVENT runs\n");
}
const TableContainer& ssfnTables = tables.getSsfnTables();
// relative permeabilty multiplier
if (!ssfnTables.empty()) {
int numRegions = ssfnTables.size();
// resize the attributes of the object
krg_.resize(numRegions);
krs_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::SsfnTable& ssfnTable = ssfnTables.getTable<SsfnTable>(regionIdx);
// Copy data
const auto& solventFraction = ssfnTable.getSolventFractionColumn();
const auto& krg = ssfnTable.getGasRelPermMultiplierColumn();
const auto& krs = ssfnTable.getSolventRelPermMultiplierColumn();
krg_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krg);
krs_[regionIdx] = NonuniformTableLinear<double>(solventFraction, krs);
}
} else {
OPM_THROW(std::runtime_error, "SSFN must be specified in SOLVENT runs\n");
}
if (deck.hasKeyword("MISCIBLE") ) {
// retrieve the cell specific Misc table index from the deck
// and using the grid...
extractTableIndex("MISCNUM", eclState, number_of_cells, global_cell, cellMiscRegionIdx_);
// misicible hydrocabon relative permeability wrt water
const TableContainer& sof2Tables = tables.getSof2Tables();
if (!sof2Tables.empty()) {
int numRegions = sof2Tables.size();
// resize the attributes of the object
krn_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::Sof2Table& sof2Table = sof2Tables.getTable<Sof2Table>(regionIdx);
// Copy data
// Sn = So + Sg + Ss;
const auto& sn = sof2Table.getSoColumn();
const auto& krn = sof2Table.getKroColumn();
krn_[regionIdx] = NonuniformTableLinear<double>(sn, krn);
}
} else {
OPM_THROW(std::runtime_error, "SOF2 must be specified in MISCIBLE (SOLVENT) runs\n");
}
const TableContainer& miscTables = tables.getMiscTables();
if (!miscTables.empty()) {
int numRegions = miscTables.size();
// resize the attributes of the object
misc_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::MiscTable& miscTable = miscTables.getTable<MiscTable>(regionIdx);
// Copy data
// solventFraction = Ss / (Ss + Sg);
const auto& solventFraction = miscTable.getSolventFractionColumn();
const auto& misc = miscTable.getMiscibilityColumn();
misc_[regionIdx] = NonuniformTableLinear<double>(solventFraction, misc);
}
} else {
OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n");
}
const TableContainer& pmiscTables = tables.getPmiscTables();
if (!pmiscTables.empty()) {
int numRegions = pmiscTables.size();
// resize the attributes of the object
pmisc_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PmiscTable& pmiscTable = pmiscTables.getTable<PmiscTable>(regionIdx);
// Copy data
const auto& po = pmiscTable.getOilPhasePressureColumn();
const auto& pmisc = pmiscTable.getMiscibilityColumn();
pmisc_[regionIdx] = NonuniformTableLinear<double>(po, pmisc);
}
}
// miscible relative permeability multipleiers
const TableContainer& msfnTables = tables.getMsfnTables();
if (!msfnTables.empty()) {
int numRegions = msfnTables.size();
// resize the attributes of the object
mkrsg_.resize(numRegions);
mkro_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::MsfnTable& msfnTable = msfnTables.getTable<MsfnTable>(regionIdx);
// Copy data
// Ssg = Ss + Sg;
const auto& Ssg = msfnTable.getGasPhaseFractionColumn();
const auto& krsg = msfnTable.getGasSolventRelpermMultiplierColumn();
const auto& kro = msfnTable.getOilRelpermMultiplierColumn();
mkrsg_[regionIdx] = NonuniformTableLinear<double>(Ssg, krsg);
mkro_[regionIdx] = NonuniformTableLinear<double>(Ssg, kro);
}
}
const TableContainer& sorwmisTables = tables.getSorwmisTables();
if (!sorwmisTables.empty()) {
int numRegions = sorwmisTables.size();
// resize the attributes of the object
sorwmis_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::SorwmisTable& sorwmisTable = sorwmisTables.getTable<SorwmisTable>(regionIdx);
// Copy data
const auto& sw = sorwmisTable.getWaterSaturationColumn();
const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn();
sorwmis_[regionIdx] = NonuniformTableLinear<double>(sw, sorwmis);
}
}
const TableContainer& sgcwmisTables = tables.getSgcwmisTables();
if (!sgcwmisTables.empty()) {
int numRegions = sgcwmisTables.size();
// resize the attributes of the object
sgcwmis_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::SgcwmisTable& sgcwmisTable = sgcwmisTables.getTable<SgcwmisTable>(regionIdx);
// Copy data
const auto& sw = sgcwmisTable.getWaterSaturationColumn();
const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn();
sgcwmis_[regionIdx] = NonuniformTableLinear<double>(sw, sgcwmis);
}
}
if (deck.hasKeyword("TLMIXPAR")) {
const int numRegions = deck.getKeyword("TLMIXPAR").size();
// resize the attributes of the object
mix_param_viscosity_.resize(numRegions);
mix_param_density_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const auto& tlmixparRecord = deck.getKeyword("TLMIXPAR").getRecord(regionIdx);
const auto& mix_params_viscosity = tlmixparRecord.getItem("TL_VISCOSITY_PARAMETER").getSIDoubleData();
mix_param_viscosity_[regionIdx] = mix_params_viscosity[0];
const auto& mix_params_density = tlmixparRecord.getItem("TL_DENSITY_PARAMETER").getSIDoubleData();
const int numDensityItems = mix_params_density.size();
if (numDensityItems == 0) {
mix_param_density_[regionIdx] = mix_param_viscosity_[regionIdx];
} else if (numDensityItems == 1) {
mix_param_density_[regionIdx] = mix_params_density[0];
} else {
OPM_THROW(std::runtime_error, "Only one value can be entered for the TL parameter pr MISC region.");
}
}
}
if (deck.hasKeyword("TLPMIXPA")) {
const TableContainer& tlpmixparTables = tables.getTlpmixpaTables();
if (!tlpmixparTables.empty()) {
int numRegions = tlpmixparTables.size();
// resize the attributes of the object
tlpmix_param_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::TlpmixpaTable& tlpmixparTable = tlpmixparTables.getTable<TlpmixpaTable>(regionIdx);
// Copy data
const auto& po = tlpmixparTable.getOilPhasePressureColumn();
const auto& tlpmixpa = tlpmixparTable.getMiscibilityColumn();
tlpmix_param_[regionIdx] = NonuniformTableLinear<double>(po, tlpmixpa);
}
} else {
// if empty keyword. Try to use the pmisc table as default.
if (pmisc_.size() > 0) {
tlpmix_param_ = pmisc_;
} else {
OPM_THROW(std::invalid_argument, "If the pressure dependent TL values in TLPMIXPA is defaulted (no entries), then the PMISC tables must be specified.");
}
}
}
}
}
}
ADB SolventPropsAdFromDeck::muSolvent(const ADB& pg,
const Cells& cells) const
{
const int n = cells.size();
assert(pg.value().size() == n);
V mu(n);
V dmudp(n);
for (int i = 0; i < n; ++i) {
const double& pg_i = pg.value()[i];
int regionIdx = cellPvtRegionIdx_[cells[i]];
double tempInvB = b_[regionIdx](pg_i);
double tempInvBmu = inverseBmu_[regionIdx](pg_i);
mu[i] = tempInvB / tempInvBmu;
dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(pg_i)
- tempInvB * inverseBmu_[regionIdx].derivative(pg_i)) / (tempInvBmu * tempInvBmu);
}
ADB::M dmudp_diag(dmudp.matrix().asDiagonal());
const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmudp_diag * pg.derivative()[block];
}
return ADB::function(std::move(mu), std::move(jacs));
}
ADB SolventPropsAdFromDeck::bSolvent(const ADB& pg,
const Cells& cells) const
{
return SolventPropsAdFromDeck::makeADBfromTables(pg, cells, cellPvtRegionIdx_, b_);
}
ADB SolventPropsAdFromDeck::gasRelPermMultiplier(const ADB& solventFraction,
const Cells& cells) const
{
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krg_);
}
ADB SolventPropsAdFromDeck::solventRelPermMultiplier(const ADB& solventFraction,
const Cells& cells) const
{
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellSatNumRegionIdx_, krs_);
}
ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn,
const Cells& cells) const
{
return SolventPropsAdFromDeck::makeADBfromTables(Sn, cells, cellSatNumRegionIdx_, krn_);
}
ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg,
const Cells& cells) const
{
if (mkrsg_.size() > 0) {
return SolventPropsAdFromDeck::makeADBfromTables(Ssg, cells, cellSatNumRegionIdx_, mkrsg_);
}
// trivial function if not specified
return Ssg;
}
ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So,
const Cells& cells) const
{
if (mkro_.size() > 0) {
return SolventPropsAdFromDeck::makeADBfromTables(So, cells, cellSatNumRegionIdx_, mkro_);
}
// trivial function if not specified
return So;
}
ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction,
const Cells& cells) const
{
return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellMiscRegionIdx_, misc_);
}
ADB SolventPropsAdFromDeck::pressureMiscibilityFunction(const ADB& po,
const Cells& cells) const
{
if (pmisc_.size() > 0) {
return SolventPropsAdFromDeck::makeADBfromTables(po, cells, cellMiscRegionIdx_, pmisc_);
}
// return ones if not specified i.e. no effect.
return ADB::constant(V::Constant(po.size(), 1.0));
}
ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw,
const Cells& cells) const {
if (sgcwmis_.size()>0) {
return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sgcwmis_);
}
// return zeros if not specified
return ADB::constant(V::Zero(Sw.size()));
}
ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw,
const Cells& cells) const {
if (sorwmis_.size()>0) {
return SolventPropsAdFromDeck::makeADBfromTables(Sw, cells, cellMiscRegionIdx_, sorwmis_);
}
// return zeros if not specified
return ADB::constant(V::Zero(Sw.size()));
}
ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD,
const Cells& cells,
const std::vector<int>& regionIdx,
const std::vector<NonuniformTableLinear<double>>& tables) const {
const int n = cells.size();
assert(X_AD.value().size() == n);
V x(n);
V dx(n);
for (int i = 0; i < n; ++i) {
const double& X_i = X_AD.value()[i];
x[i] = tables[regionIdx[cells[i]]](X_i);
dx[i] = tables[regionIdx[cells[i]]].derivative(X_i);
}
ADB::M dx_diag(dx.matrix().asDiagonal());
const int num_blocks = X_AD.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dx_diag, X_AD.derivative()[block], jacs[block]);
}
return ADB::function(std::move(x), std::move(jacs));
}
V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
const int n = cells.size();
V density(n);
for (int i = 0; i < n; ++i) {
int regionIdx = cellPvtRegionIdx_[cells[i]];
density[i] = solvent_surface_densities_[regionIdx];
}
return density;
}
V SolventPropsAdFromDeck::mixingParameterViscosity(const Cells& cells) const {
const int n = cells.size();
if (mix_param_viscosity_.size() > 0) {
V mix_param(n);
for (int i = 0; i < n; ++i) {
int regionIdx = cellMiscRegionIdx_[cells[i]];
mix_param[i] = mix_param_viscosity_[regionIdx];
}
return mix_param;
}
// return zeros if not specified
return V::Zero(n);
}
V SolventPropsAdFromDeck::mixingParameterDensity(const Cells& cells) const {
const int n = cells.size();
if (mix_param_viscosity_.size() > 0) {
V mix_param(n);
for (int i = 0; i < n; ++i) {
int regionIdx = cellMiscRegionIdx_[cells[i]];
mix_param[i] = mix_param_density_[regionIdx];
}
return mix_param;
}
// return zeros if not specified
return V::Zero(n);
}
ADB SolventPropsAdFromDeck::pressureMixingParameter(const ADB& po,
const Cells& cells) const {
if (tlpmix_param_.size() > 0) {
return SolventPropsAdFromDeck::makeADBfromTables(po, cells, cellMiscRegionIdx_, tlpmix_param_);
}
// return ones if not specified i.e. no pressure effects.
return ADB::constant(V::Constant(po.size(), 1.0));
}
void SolventPropsAdFromDeck::extractTableIndex(const std::string& keyword,
const Opm::EclipseState& eclState,
size_t numCompressed,
const int* compressedToCartesianCellIdx,
std::vector<int>& tableIdx) const {
//Get the Region data
const auto& regionData = eclState.get3DProperties().getIntGridProperty(keyword).getData();
// Convert this into an array of compressed cells
// Eclipse uses Fortran-style indices which start at 1
// instead of 0, we subtract 1.
tableIdx.resize(numCompressed);
for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) {
size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx;
assert(cartesianCellIdx < regionData.size());
tableIdx[cellIdx] = regionData[cartesianCellIdx] - 1;
}
}
} //namespace OPM

202
opm/autodiff/SolventPropsAdFromDeck.hpp
View File

@ -1,202 +0,0 @@
/*
Copyright 2015 IRIS
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/>.
*/
#ifndef SOLVENTPROPSADFROMDECK_HPP
#define SOLVENTPROPSADFROMDECK_HPP
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <cmath>
#include <vector>
#include <opm/common/ErrorMacros.hpp>
namespace Opm
{
class SolventPropsAdFromDeck
{
public:
SolventPropsAdFromDeck(const Deck& deck,
const EclipseState& eclipseState,
const int number_of_cells,
const int* global_cell);
////////////////////////////
// Fluid interface //
////////////////////////////
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef std::vector<int> Cells;
/// Solvent formation volume factor.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n formation volume factor values.
ADB bSolvent(const ADB& pg,
const Cells& cells) const;
/// Solvent viscosity.
/// \param[in] pg Array of n gas pressure values.
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n viscosity values.
ADB muSolvent(const ADB& pg,
const Cells& cells) const;
/// Gas relPerm multipliers
/// \param[in] gasFraction Array of n gas fraction Sg / (sg + Ss) values.
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \return Array of n gas relPerm multiplier values.
ADB gasRelPermMultiplier(const ADB& solventFraction,
const Cells& cells) const;
/// Solvent relPerm multipliers
/// \param[in] solventFraction Array of n solvent fraction Ss / (Sg + Ss) values.
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \return Array of n solvent relPerm multiplier values.
ADB solventRelPermMultiplier(const ADB& solventFraction,
const Cells& cells) const;
/// Miscible hydrocrabon relPerm wrt water
/// \param[in] Sn Array of n total hyrdrocarbon saturation values.
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \return Array of n miscible hyrdrocabon wrt water relPerm values.
ADB misicibleHydrocarbonWaterRelPerm(const ADB& Sn,
const Cells& cells) const;
/// Miscible Solvent + Gas relPerm multiplier
/// \param[in] Ssg Array of n total gas fraction (Sgas + Ssolvent) / Sn values, where
/// Sn = Sgas + Ssolvent + Soil.
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \return Array of n solvent gas relperm multiplier.
ADB miscibleSolventGasRelPermMultiplier(const ADB& Ssg,
const Cells& cells) const;
/// Miscible Oil relPerm multiplier
/// \param[in] So Array of n oil fraction values. Soil / Sn values, where Sn = Sgas + Ssolvent + Soil.
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \return Array of n oil relperm multiplier.
ADB miscibleOilRelPermMultiplier(const ADB& So,
const Cells& cells) const;
/// Miscible function
/// \param[in] solventFraction Array of n solvent fraction Ss / (Sg + Ss) values.
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \return Array of n miscibility values
ADB miscibilityFunction(const ADB& solventFraction,
const Cells& cells) const;
/// Pressure dependent miscibility function
/// \param[in] solventFraction Array of n oil phase pressure .
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n miscibility values
ADB pressureMiscibilityFunction(const ADB& po,
const Cells& cells) const;
/// Miscible critical gas saturation function
/// \param[in] Sw Array of n water saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return Array of n miscible critical gas saturation values
ADB miscibleCriticalGasSaturationFunction(const ADB& Sw,
const Cells& cells) const;
/// Miscible residual oil saturation function
/// \param[in] Sw Array of n water saturation values.
/// \param[in] cells Array of n cell indices to be associated with the saturation values.
/// \return Array of n miscible residual oil saturation values
ADB miscibleResidualOilSaturationFunction(const ADB& Sw,
const Cells& cells) const;
/// Solvent surface density
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \return Array of n solvent density values.
V solventSurfaceDensity(const Cells& cells) const;
/// Todd-Longstaff mixing parameter for viscosity calculation
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// return Array of n mixing paramters for viscosity calculation
V mixingParameterViscosity(const Cells& cells) const;
/// Todd-Longstaff mixing parameter for density calculation
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// return Array of n mixing paramters for density calculation
V mixingParameterDensity(const Cells& cells) const;
/// Todd-Longstaff pressure dependent mixing parameter
/// \param[in] po Array of n oil pressure values
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// return Array of n pressure dependent mixing paramters
ADB pressureMixingParameter(const ADB& po,
const Cells& cells) const;
private:
/// Makes ADB from table values
/// \param[in] X Array of n table lookup values.
/// \param[in] cells Array of n cell indices to be associated with the fraction values.
/// \param[in] tables Vector of tables, one for each PVT region.
/// \return Array of n solvent density values.
ADB makeADBfromTables(const ADB& X,
const Cells& cells,
const std::vector<int>& regionIdx,
const std::vector<NonuniformTableLinear<double>>& tables) const;
/// Helper function to create an array containing the
/// table index of for each compressed cell from an Eclipse deck.
/// \param[in] keyword eclKeyword specifying region (SATNUM etc. )
/// \param[in] eclState eclState from opm-parser
/// \param[in] numCompressed number of compressed cells
/// \param[in] compressedToCartesianCellIdx cartesianCellIdx for each cell in the grid
/// \param[out] tableIdx table index for each compressed cell
void extractTableIndex(const std::string& keyword,
const Opm::EclipseState& eclState,
size_t numCompressed,
const int* compressedToCartesianCellIdx,
std::vector<int>& tableIdx) const;
// The PVT region which is to be used for each cell
std::vector<int> cellPvtRegionIdx_;
std::vector<int> cellMiscRegionIdx_;
std::vector<int> cellSatNumRegionIdx_;
std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_;
std::vector<NonuniformTableLinear<double> > inverseBmu_;
std::vector<double> solvent_surface_densities_;
std::vector<NonuniformTableLinear<double> > krg_;
std::vector<NonuniformTableLinear<double> > krs_;
std::vector<NonuniformTableLinear<double> > krn_;
std::vector<NonuniformTableLinear<double> > mkro_;
std::vector<NonuniformTableLinear<double> > mkrsg_;
std::vector<NonuniformTableLinear<double> > misc_;
std::vector<NonuniformTableLinear<double> > pmisc_;
std::vector<NonuniformTableLinear<double> > sorwmis_;
std::vector<NonuniformTableLinear<double> > sgcwmis_;
std::vector<NonuniformTableLinear<double> > tlpmix_param_;
std::vector<double> mix_param_viscosity_;
std::vector<double> mix_param_density_;
};
} // namespace OPM
#endif // SOLVENTPROPSADFROMDECK_HPP

87
opm/autodiff/StandardWellsSolvent.hpp
View File

@ -1,87 +0,0 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
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/>.
*/
#ifndef OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED
#define OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED
#include <opm/autodiff/StandardWells.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
namespace Opm {
/// Class for handling the standard well model for solvent model
class StandardWellsSolvent : public StandardWells
{
public:
using Base = StandardWells;
using Base::computeWellConnectionDensitesPressures;
// --------- Public methods ---------
StandardWellsSolvent(const Wells* wells_arg, WellCollection* well_collection);
// added the Solvent related
void initSolvent(const SolventPropsAdFromDeck* solvent_props,
const int solvent_pos,
const bool has_solvent);
template <class SolutionState>
void computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
std::vector<ADB>& cq_s) const;
template <class SolutionState, class WellState>
void computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf);
// TODO: fluid and active may be can put in the member list
template <class ReservoirResidualQuant, class SolutionState>
void extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const;
template <class SolutionState, class WellState>
void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw);
protected:
const SolventPropsAdFromDeck* solvent_props_;
int solvent_pos_;
bool has_solvent_;
using Base::phase_condition_;
};
} // namespace Opm
#include "StandardWellsSolvent_impl.hpp"
#endif

458
opm/autodiff/StandardWellsSolvent_impl.hpp
View File

@ -1,458 +0,0 @@
/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
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/>.
*/
#include <opm/autodiff/StandardWellsSolvent.hpp>
namespace Opm
{
StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, WellCollection* well_collection)
: Base(wells_arg, well_collection)
, solvent_props_(nullptr)
, solvent_pos_(-1)
, has_solvent_(false)
{
}
void
StandardWellsSolvent::initSolvent(const SolventPropsAdFromDeck* solvent_props,
const int solvent_pos,
const bool has_solvent)
{
solvent_props_ = solvent_props;
solvent_pos_ = solvent_pos;
has_solvent_ = has_solvent;
}
template<class SolutionState, class WellState>
void
StandardWellsSolvent::
computePropertiesForWellConnectionPressures(const SolutionState& state,
const WellState& xw,
std::vector<double>& b_perf,
std::vector<double>& rsmax_perf,
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf)
{
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects.
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
// Compute the average pressure in each well block
const Vector perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
Vector avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
const std::vector<int>& well_cells = wellOps().well_cells;
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
const ADB avg_press_ad = ADB::constant(avg_press);
std::vector<PhasePresence> perf_cond(nperf);
for (int perf = 0; perf < nperf; ++perf) {
perf_cond[perf] = (*phase_condition_)[well_cells[perf]];
}
const PhaseUsage& pu = fluid_->phaseUsage();
DataBlock b(nperf, pu.num_phases);
const Vector bw = fluid_->bWat(avg_press_ad, perf_temp, well_cells).value();
if (pu.phase_used[BlackoilPhases::Aqua]) {
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
}
assert((*active_)[Oil]);
assert((*active_)[Gas]);
const ADB perf_rv = subset(state.rv, well_cells);
const ADB perf_rs = subset(state.rs, well_cells);
const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) {
const Vector bo = fluid_->bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value();
//const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value();
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
// const Vector rssat = fluidRsSat(avg_press, perf_so, well_cells);
const Vector rssat = fluid_->rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
} else {
rsmax_perf.assign(0.0, nperf);
}
V surf_dens_copy = superset(fluid_->surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) {
continue; // the gas surface density is added after the solvent is accounted for.
}
surf_dens_copy += superset(fluid_->surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
// Unclear wether the effective or the pure values should be used for the wells
// the current usage of unmodified properties values gives best match.
//V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value();
Vector bg = fluid_->bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
Vector rhog = fluid_->surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
// to handle solvent related
if (has_solvent_) {
const Vector bs = solvent_props_->bSolvent(avg_press_ad,well_cells).value();
//const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();
// number of cells
const int nc = state.pressure.size();
const ADB zero = ADB::constant(Vector::Zero(nc));
const ADB& ss = state.solvent_saturation;
const ADB& sg = ((*active_)[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value();
Vector injectedSolventFraction = Eigen::Map<const Vector>(&xw.solventFraction()[0], nperf);
Vector isProducer = Vector::Zero(nperf);
Vector ones = Vector::Constant(nperf,1.0);
for (int w = 0; w < nw; ++w) {
if(wells().type[w] == PRODUCER) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
isProducer[perf] = 1;
}
}
}
F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction;
bg = bg * (ones - F_solvent);
bg = bg + F_solvent * bs;
const Vector& rhos = solvent_props_->solventSurfaceDensity(well_cells);
rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
}
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases);
// const Vector rvsat = fluidRvSat(avg_press, perf_so, well_cells);
const Vector rvsat = fluid_->rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value();
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
} else {
rvmax_perf.assign(0.0, nperf);
}
// b and surf_dens_perf is row major, so can just copy data.
b_perf.assign(b.data(), b.data() + nperf * pu.num_phases);
surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
}
template <class SolutionState>
void
StandardWellsSolvent::
computeWellFlux(const SolutionState& state,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
Vector& aliveWells,
std::vector<ADB>& cq_s) const
{
if( ! localWellsActive() ) return ;
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
Vector Tw = Eigen::Map<const Vector>(wells().WI, nperf);
const std::vector<int>& well_cells = wellOps().well_cells;
// pressure diffs computed already (once per step, not changing per iteration)
const Vector& cdp = wellPerforationPressureDiffs();
// Extract needed quantities for the perforation cells
const ADB& p_perfcells = subset(state.pressure, well_cells);
// Perforation pressure
const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp;
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcells - perfpressure;
// Compute vectors with zero and ones that
// selects the wanted quantities.
// selects injection perforations
Vector selectInjectingPerforations = Vector::Zero(nperf);
// selects producing perforations
Vector selectProducingPerforations = Vector::Zero(nperf);
for (int c = 0; c < nperf; ++c){
if (drawdown.value()[c] < 0)
selectInjectingPerforations[c] = 1;
else
selectProducingPerforations[c] = 1;
}
// Handle cross flow
const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value();
const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value();
for (int w = 0; w < nw; ++w) {
if (!wells().allow_cf[w]) {
for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
// Crossflow is not allowed; reverse flow is prevented.
// At least one of the perforation must be open in order to have a meeningful
// equation to solve. For the special case where all perforations have reverse flow,
// and the target rate is non-zero all of the perforations are keept open.
if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
selectProducingPerforations[perf] = 0.0;
} else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
selectInjectingPerforations[perf] = 0.0;
}
}
}
}
// HANDLE FLOW INTO WELLBORE
// compute phase volumetric rates at standard conditions
std::vector<ADB> cq_p(np, ADB::null());
std::vector<ADB> cq_ps(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
cq_p[phase] = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown);
cq_ps[phase] = b_perfcells[phase] * cq_p[phase];
}
Vector ones = Vector::Constant(nperf,1.0);
ADB F_gas = ADB::constant(ones);
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
if ((*active_)[Oil] && (*active_)[Gas]) {
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
const ADB cq_psOil = cq_ps[oilpos];
ADB cq_psGas = cq_ps[gaspos];
const ADB& rv_perfcells = subset(state.rv, well_cells);
const ADB& rs_perfcells = subset(state.rs, well_cells);
cq_ps[gaspos] += rs_perfcells * cq_psOil;
if(has_solvent_) {
// The solvent gas need to be removed from the gas
// before multiplied with rv.
const ADB& ss = state.solvent_saturation;
const ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
F_gas -= subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
cq_psGas = cq_psGas * F_gas;
}
cq_ps[oilpos] += rv_perfcells * cq_psGas;
}
// HANDLE FLOW OUT FROM WELLBORE
// Using total mobilities
ADB total_mob = mob_perfcells[0];
for (int phase = 1; phase < np; ++phase) {
total_mob += mob_perfcells[phase];
}
// injection perforations total volume rates
const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown);
// Store well perforation total fluxes (reservor volumes) if requested.
if (store_well_perforation_fluxes_) {
// Ugly const-cast, but unappealing alternatives.
Vector& wf = const_cast<Vector&>(well_perforation_fluxes_);
wf = cqt_i.value();
for (int phase = 0; phase < np; ++phase) {
wf += cq_p[phase].value();
}
}
// compute wellbore mixture for injecting perforations
// The wellbore mixture depends on the inflow from the reservoar
// and the well injection rates.
// compute avg. and total wellbore phase volumetric rates at standard conds
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
std::vector<ADB> wbq(np, ADB::null());
ADB wbqt = ADB::constant(Vector::Zero(nw));
for (int phase = 0; phase < np; ++phase) {
const ADB& q_ps = wellOps().p2w * cq_ps[phase];
const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw));
Selector<double> injectingPhase_selector(q_s.value(), Selector<double>::GreaterZero);
const int pos = pu.phase_pos[phase];
wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps;
wbqt += wbq[phase];
}
// compute wellbore mixture at standard conditions.
Selector<double> notDeadWells_selector(wbqt.value(), Selector<double>::Zero);
std::vector<ADB> cmix_s(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase];
cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt);
}
// compute volume ratio between connection at standard conditions
ADB volumeRatio = ADB::constant(Vector::Zero(nperf));
if ((*active_)[Water]) {
const int watpos = pu.phase_pos[Water];
volumeRatio += cmix_s[watpos] / b_perfcells[watpos];
}
if ((*active_)[Oil] && (*active_)[Gas]) {
// Incorporate RS/RV factors if both oil and gas active
const ADB& rv_perfcells = subset(state.rv, well_cells);
const ADB& rs_perfcells = subset(state.rs, well_cells);
const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells;
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
const ADB tmp_oil = (cmix_s[oilpos] - rv_perfcells * F_gas * cmix_s[gaspos]) / d;
volumeRatio += tmp_oil / b_perfcells[oilpos];
const ADB tmp_gas = (cmix_s[gaspos] - rs_perfcells * cmix_s[oilpos]) / d;
volumeRatio += tmp_gas / b_perfcells[gaspos];
}
else {
if ((*active_)[Oil]) {
const int oilpos = pu.phase_pos[Oil];
volumeRatio += cmix_s[oilpos] / b_perfcells[oilpos];
}
if ((*active_)[Gas]) {
const int gaspos = pu.phase_pos[Gas];
volumeRatio += cmix_s[gaspos] / b_perfcells[gaspos];
}
}
// injecting connections total volumerates at standard conditions
ADB cqt_is = cqt_i/volumeRatio;
// connection phase volumerates at standard conditions
cq_s.resize(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is;
}
// check for dead wells (used in the well controll equations)
aliveWells = Vector::Constant(nw, 1.0);
for (int w = 0; w < nw; ++w) {
if (wbqt.value()[w] == 0) {
aliveWells[w] = 0.0;
}
}
}
template <class SolutionState, class WellState>
void
StandardWellsSolvent::
computeWellConnectionPressures(const SolutionState& state,
const WellState& xw)
{
if( ! localWellsActive() ) return ;
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects.
std::vector<double> b_perf;
std::vector<double> rsmax_perf;
std::vector<double> rvmax_perf;
std::vector<double> surf_dens_perf;
computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
const Vector pdepth = perf_cell_depth_;
const int nperf = wells().well_connpos[wells().number_of_wells];
const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
computeWellConnectionDensitesPressures(xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity_);
}
template <class ReservoirResidualQuant, class SolutionState>
void
StandardWellsSolvent::
extractWellPerfProperties(const SolutionState& state,
const std::vector<ReservoirResidualQuant>& rq,
std::vector<ADB>& mob_perfcells,
std::vector<ADB>& b_perfcells) const
{
Base::extractWellPerfProperties(state, rq, mob_perfcells, b_perfcells);
// handle the solvent related
if (has_solvent_) {
const Opm::PhaseUsage& pu = fluid_->phaseUsage();
int gas_pos = pu.phase_pos[Gas];
const std::vector<int>& well_cells = wellOps().well_cells;
const int nperf = well_cells.size();
// Gas and solvent is combinded and solved together
// The input in the well equation is then the
// total gas phase = hydro carbon gas + solvent gas
// The total mobility is the sum of the solvent and gas mobiliy
mob_perfcells[gas_pos] += subset(rq[solvent_pos_].mob, well_cells);
// A weighted sum of the b-factors of gas and solvent are used.
const int nc = rq[solvent_pos_].mob.size();
const ADB zero = ADB::constant(Vector::Zero(nc));
const ADB& ss = state.solvent_saturation;
const ADB& sg = ((*active_)[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ]
: zero);
Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells);
Vector ones = Vector::Constant(nperf,1.0);
b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos];
b_perfcells[gas_pos] += (F_solvent * subset(rq[solvent_pos_].b, well_cells));
}
}
}

73
opm/autodiff/WellStateFullyImplicitBlackoilSolvent.hpp
View File

@ -1,73 +0,0 @@
/*
Copyright 2015 IRIS 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/>.
*/
#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOILSOLVENT_HEADER_INCLUDED
#define OPM_WELLSTATEFULLYIMPLICITBLACKOILSOLVENT_HEADER_INCLUDED
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
namespace Opm
{
class WellStateFullyImplicitBlackoilSolvent : public WellStateFullyImplicitBlackoil
{
typedef WellStateFullyImplicitBlackoil BaseType;
public:
/// One solvent fraction per well connection
std::vector<double>& solventFraction() { return solvent_fraction_; }
const std::vector<double>& solventFraction() const { return solvent_fraction_; }
data::Wells report(const PhaseUsage &pu) const override {
data::Wells res = WellStateFullyImplicitBlackoil::report(pu);
const int nw = WellState::numWells();
// If there are now wells numPhases throws a floating point
// exception.
if (nw == 0) {
return res;
}
const int np = BaseType::numPhases();
assert( np == 3 ); // the solvent model assumes 3 phases in the base model
// completions aren't supported yet
for( auto w = 0; w < nw; ++w ) {
using rt = data::Rates::opt;
double solvent_well_rate = 0.0;
for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
auto solvent_rate_this = BaseType::perfPhaseRates()[np*perf + pu.phase_pos[BlackoilPhases::Vapour]] * solventFraction()[perf];
solvent_well_rate += solvent_rate_this;
}
res.at( wells_->name[ w ]).rates.set( rt::solvent, solvent_well_rate );
}
return res;
}
private:
std::vector<double> solvent_fraction_;
};
} // namespace Opm
#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOILSOLVENT_HEADER_INCLUDED

244
tests/test_solventprops_ad.cpp
View File

@ -1,244 +0,0 @@
/*
Copyright 2016 IRIS 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/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE SolventPropertiesTest
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <fstream>
#include <iostream>
const std::string deckData = "\n\
RUNSPEC \n\
\n\
SOLVENT \n\
\n\
MISCIBLE\n\
1 3 /\n\
\n\
DIMENS \n\
3 1 1 \n\
/\n\
TABDIMS\n\
/\n\
GRID \n\
\n\
DXV \n\
1 1 1\n\
/\n\
DYV \n\
1 \n\
/\n\
DZV \n\
1 \n\
/\n\
TOPS \n\
1 1 1 \n\
/\n";
const std::string solventData = "\n\
SDENSITY \n\
0.1 / \n\
PVDS \n\
1 1 0.1 / \n\
SSFN \n\
0.0 0.0 0.0 \n\
1.0 1.0 1.0 \n\
/ \n\
MISC \n\
0.0 0.0 \n\
1.0 1.0 \n\
/ \n\
SOF2 \n\
0 0 \n\
0.88 1 / \n";
BOOST_AUTO_TEST_CASE(Construction)
{
Opm::ParseContext parseContext;
Opm::Parser parser;
auto deck = parser.parseString(deckData + solventData, parseContext);
Opm::EclipseState eclState(deck , parseContext);
std::vector<int> global_ind = {0 , 1 , 2};
Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind.data());
}
BOOST_AUTO_TEST_CASE(SolventData)
{
Opm::ParseContext parseContext;
Opm::Parser parser;
auto deck = parser.parseString(deckData + solventData, parseContext);
Opm::EclipseState eclState(deck , parseContext);
std::vector<int> global_ind = {0 , 1 , 2};
Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind.data());
const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
typedef Opm::SolventPropsAdFromDeck::V V;
V rho = solventprops.solventSurfaceDensity(cells);
BOOST_REQUIRE_EQUAL(rho.size(), cells.size());
BOOST_CHECK_EQUAL(rho[0], 0.1);
BOOST_CHECK_EQUAL(rho[0], rho[1]);
}
const std::string pmiscData = "\n\
PMISC\n\
100 0.0 \n\
200 0.0 \n\
500 1.0 \n\
1000 1.0 /\n\
\n";
BOOST_AUTO_TEST_CASE(PMISC)
{
Opm::ParseContext parseContext;
Opm::Parser parser;
auto deck = parser.parseString(deckData + solventData + pmiscData, parseContext);
Opm::EclipseState eclState(deck , parseContext);
const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
typedef Opm::SolventPropsAdFromDeck::V V;
std::vector<int> global_ind = {0 , 1 , 2};
Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind.data());
V po(3);
po << 150,250,550;
po = po * Opm::unit::barsa;
BOOST_REQUIRE_EQUAL(po.size(), cells.size());
V pmisc = solventprops.pressureMiscibilityFunction(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
BOOST_REQUIRE_EQUAL(pmisc.size(), cells.size());
BOOST_CHECK_EQUAL(pmisc[0], 0.0);
const double tol = 1e-12;
const double value = (250.0 - 200.0) / (500.0 - 200.0); // linear interpolation
BOOST_CHECK_CLOSE(pmisc[1], value, tol);
BOOST_CHECK_EQUAL(pmisc[2], 1.0);
}
const std::string tlpmixpaData = "\n\
TLPMIXPA\n\
100 0.0 \n\
200 0.0 \n\
500 1.0 \n\
1000 1.0 /\n\
\n";
BOOST_AUTO_TEST_CASE(TLPMIXPA)
{
Opm::ParseContext parseContext;
Opm::Parser parser;
auto deck = parser.parseString(deckData + solventData + tlpmixpaData, parseContext);
Opm::EclipseState eclState(deck , parseContext);
const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
typedef Opm::SolventPropsAdFromDeck::V V;
const int* global_ind = new int[3] {0 , 1 , 2};
Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind);
V po(3);
po << 150,250,550;
po = po * Opm::unit::barsa;
BOOST_REQUIRE_EQUAL(po.size(), cells.size());
V tlpmixpa = solventprops.pressureMixingParameter(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
BOOST_REQUIRE_EQUAL(tlpmixpa.size(), cells.size());
BOOST_CHECK_EQUAL(tlpmixpa[0], 0.0);
const double tol = 1e-12;
const double value = (250.0 - 200.0) / (500.0 - 200.0); // linear interpolation
BOOST_CHECK_CLOSE(tlpmixpa[1], value, tol);
BOOST_CHECK_EQUAL(tlpmixpa[2], 1.0);
}
BOOST_AUTO_TEST_CASE(TLPMIXPA_NOT_SPECIFIED)
{
Opm::ParseContext parseContext;
Opm::Parser parser;
// no pmisc data and default tlpmixdata i.e it should throw
auto deck = parser.parseString(deckData + solventData, parseContext);
Opm::EclipseState eclState(deck , parseContext);
const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
const int* global_ind = new int[3] {0 , 1 , 2};
Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind);
typedef Opm::SolventPropsAdFromDeck::V V;
V po(3);
po << 150,250,550;
po = po * Opm::unit::barsa;
BOOST_REQUIRE_EQUAL(po.size(), cells.size());
V tlpmixpa = solventprops.pressureMixingParameter(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
BOOST_REQUIRE_EQUAL(tlpmixpa.size(), cells.size());
// if not specified tlpmixpa is 1.0 for all cells.
BOOST_CHECK_EQUAL(tlpmixpa[0], 1.0);
BOOST_CHECK_EQUAL(tlpmixpa[1], 1.0);
BOOST_CHECK_EQUAL(tlpmixpa[2], 1.0);
}
const std::string tlpmixpaDataDefault = "\n\
TLPMIXPA\n\
/\n\
\n";
BOOST_AUTO_TEST_CASE(TLPMIXPA_DEFAULT)
{
Opm::ParseContext parseContext;
Opm::Parser parser;
auto deck = parser.parseString(deckData + solventData + pmiscData + tlpmixpaDataDefault, parseContext);
Opm::EclipseState eclState(deck , parseContext);
const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
typedef Opm::SolventPropsAdFromDeck::V V;
const int* global_ind = new int[3] {0 , 1 , 2};
Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind);
V po(3);
po << 150,250,550;
po = po * Opm::unit::barsa;
BOOST_REQUIRE_EQUAL(po.size(), cells.size());
V tlpmixpa = solventprops.pressureMixingParameter(Opm::SolventPropsAdFromDeck::ADB::constant(po), cells).value();
BOOST_REQUIRE_EQUAL(tlpmixpa.size(), cells.size());
BOOST_CHECK_EQUAL(tlpmixpa[0], 0.0);
const double tol = 1e-12;
const double value = (250.0 - 200.0) / (500.0 - 200.0); // linear interpolation
BOOST_CHECK_CLOSE(tlpmixpa[1], value, tol);
BOOST_CHECK_EQUAL(tlpmixpa[2], 1.0);
}
BOOST_AUTO_TEST_CASE(TLPMIXPA_DEFAULT_NOPMISC)
{
Opm::ParseContext parseContext;
Opm::Parser parser;
// no pmisc data and default tlpmixdata i.e it should throw
auto deck = parser.parseString(deckData + solventData + tlpmixpaDataDefault, parseContext);
Opm::EclipseState eclState(deck , parseContext);
const Opm::SolventPropsAdFromDeck::Cells cells(3, 0);
const int* global_ind = new int[3] {0 , 1 , 2};
BOOST_CHECK_THROW(Opm::SolventPropsAdFromDeck solventprops(deck, eclState, 3, global_ind), std::invalid_argument);
}