OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-12 08:54:53 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1344 from atgeirr/remove-legacy-polymer

Browse Source 
Remove legacy fully implicit polymer simulator.
This commit is contained in:
Atgeirr Flø Rasmussen authored and GitHub committed 2017-11-27 07:19:01 +01:00
2 parents 3fd27d9bd6 eabd06389c
commit 369e9f5a06
11 files changed
+7 -2274
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

CMakeLists_files.cmake
+2 -11
View File
@@ -58,8 +58,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/polymer/SimulatorPolymer.cpp opm/polymer/SimulatorPolymer.cpp
opm/polymer/TransportSolverTwophaseCompressiblePolymer.cpp opm/polymer/TransportSolverTwophaseCompressiblePolymer.cpp
opm/polymer/TransportSolverTwophasePolymer.cpp opm/polymer/TransportSolverTwophasePolymer.cpp
opm/polymer/fullyimplicit/PolymerPropsAd.cpp opm/simulators/flow_ebos_blackoil.cpp
opm/simulators/flow_ebos_blackoil.cpp
opm/simulators/flow_ebos_gasoil.cpp opm/simulators/flow_ebos_gasoil.cpp
opm/simulators/flow_ebos_oilwater.cpp opm/simulators/flow_ebos_oilwater.cpp
opm/simulators/flow_ebos_polymer.cpp opm/simulators/flow_ebos_polymer.cpp
@@ -122,7 +121,6 @@ list (APPEND EXAMPLE_SOURCE_FILES
examples/sim_simple.cpp examples/sim_simple.cpp
examples/sim_poly2p_comp_reorder.cpp examples/sim_poly2p_comp_reorder.cpp
examples/sim_poly2p_incomp_reorder.cpp examples/sim_poly2p_incomp_reorder.cpp
examples/flow_polymer.cpp
examples/wells_example.cpp examples/wells_example.cpp
) )
@@ -139,7 +137,6 @@ list (APPEND PROGRAM_SOURCE_FILES
examples/opm_init_check.cpp examples/opm_init_check.cpp
examples/sim_poly2p_comp_reorder.cpp examples/sim_poly2p_comp_reorder.cpp
examples/sim_poly2p_incomp_reorder.cpp examples/sim_poly2p_incomp_reorder.cpp
examples/flow_polymer.cpp
) )
# originally generated with the command: # originally generated with the command:
@@ -240,13 +237,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/polymer/TransportSolverTwophaseCompressiblePolymer.hpp opm/polymer/TransportSolverTwophaseCompressiblePolymer.hpp
opm/polymer/Point2D.hpp opm/polymer/Point2D.hpp
opm/polymer/TransportSolverTwophasePolymer.hpp opm/polymer/TransportSolverTwophasePolymer.hpp
opm/polymer/fullyimplicit/PolymerPropsAd.hpp opm/simulators/flow_ebos_blackoil.hpp
opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp
opm/simulators/flow_ebos_blackoil.hpp
opm/simulators/flow_ebos_gasoil.hpp opm/simulators/flow_ebos_gasoil.hpp
opm/simulators/flow_ebos_oilwater.hpp opm/simulators/flow_ebos_oilwater.hpp
opm/simulators/flow_ebos_polymer.hpp opm/simulators/flow_ebos_polymer.hpp
compareECLFiles.cmake
+5 -7
View File
@@ -179,13 +179,11 @@ add_test_compareECLFiles(CASENAME msw_3d_hfa
REL_TOL ${rel_tol} REL_TOL ${rel_tol}
TEST_ARGS use_multisegment_well=true) TEST_ARGS use_multisegment_well=true)
foreach(SIM flow flow_polymer) add_test_compareECLFiles(CASENAME polymer_simple2D
add_test_compareECLFiles(CASENAME polymer_simple2D FILENAME 2D_THREEPHASE_POLY_HETER
FILENAME 2D_THREEPHASE_POLY_HETER SIMULATOR flow
SIMULATOR ${SIM} ABS_TOL ${abs_tol}
ABS_TOL ${abs_tol} REL_TOL ${coarse_rel_tol})
REL_TOL ${coarse_rel_tol})
endforeach()
add_test_compareECLFiles(CASENAME spe5 add_test_compareECLFiles(CASENAME spe5
FILENAME SPE5CASE1 FILENAME SPE5CASE1
examples/flow_polymer.cpp
-42
View File
@@ -1,42 +0,0 @@
/*
Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014 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/>.
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/grid.h>
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
#include <opm/autodiff/FlowMainPolymer.hpp>
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
typedef UnstructuredGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoilPolymer<Grid> Simulator;
Opm::FlowMainPolymer<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}
opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
-308
View File
@@ -1,308 +0,0 @@
/*
Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014 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_BLACKOILPOLYMERMODEL_HEADER_INCLUDED
#define OPM_BLACKOILPOLYMERMODEL_HEADER_INCLUDED
#include <opm/autodiff/BlackoilModelBase.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
#include <opm/autodiff/StandardWells.hpp>
#include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
namespace Opm {
/// A model implementation for three-phase black oil with polymer.
///
/// The simulator is capable of handling three-phase problems
/// where gas can be dissolved in oil and vice versa, with polymer
/// in the water phase. It uses an industry-standard TPFA
/// discretization with per-phase upwind weighting of mobilities.
///
/// It uses automatic differentiation via the class AutoDiffBlock
/// to simplify assembly of the jacobian matrix.
template<class Grid>
class BlackoilPolymerModel : public BlackoilModelBase<Grid, StandardWells, BlackoilPolymerModel<Grid> >
{
public:
// --------- Types and enums ---------
typedef BlackoilModelBase<Grid, StandardWells, BlackoilPolymerModel<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, StandardWells, BlackoilPolymerModel<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] 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] has_polymer turn on polymer feature
/// \param[in] has_plyshlog true when PLYSHLOG keyword available
/// \param[in] has_shrate true when PLYSHLOG keyword available
/// \param[in] wells_rep_radius representative radius of well perforations during shear effects calculation
/// \param[in] wells_perf_length perforation length for well perforations
/// \param[in] wells_bore_diameter wellbore diameters for well performations
/// \param[in] terminal_output request output to cout/cerr
BlackoilPolymerModel(const typename Base::ModelParameters& param,
const Grid& grid,
const BlackoilPropsAdFromDeck& fluid,
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const StandardWells& well_model,
const NewtonIterationBlackoilInterface& linsolver,
std::shared_ptr< const EclipseState > eclipseState,
std::shared_ptr< const Schedule> schedule,
std::shared_ptr< const SummaryConfig> summary_config,
const bool has_disgas,
const bool has_vapoil,
const bool has_polymer,
const bool has_plyshlog,
const bool has_shrate,
const std::vector<double>& wells_rep_radius,
const std::vector<double>& wells_perf_length,
const std::vector<double>& wells_bore_diameter,
const bool terminal_output);
/// Called once before each time step.
/// \param[in] timer simulation timer
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void prepareStep(const SimulatorTimerInterface& timer,
const ReservoirState& reservoir_state,
const WellState& well_state);
/// Called once after each time step.
/// \param[in] timer simulation timer
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void afterStep(const SimulatorTimerInterface& timer,
ReservoirState& reservoir_state,
WellState& well_state);
/// 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);
/// Assemble the residual and Jacobian of the nonlinear system.
/// \param[in] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
/// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep
SimulatorReport
assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly);
using Base::wellModel;
protected:
// --------- Types and enums ---------
typedef typename Base::SolutionState SolutionState;
typedef typename Base::DataBlock DataBlock;
enum { Concentration = CanonicalVariablePositions::Next };
// --------- Data members ---------
const PolymerPropsAd& polymer_props_ad_;
const bool has_polymer_;
const bool has_plyshlog_;
const bool has_shrate_;
const int poly_pos_;
V cmax_;
// representative radius and perforation length of well perforations
// to be used in shear-thinning computation.
std::vector<double> wells_rep_radius_;
std::vector<double> wells_perf_length_;
// wellbore diameters
std::vector<double> wells_bore_diameter_;
// shear-thinning factor for cell faces
std::vector<double> shear_mult_faces_;
// shear-thinning factor for well perforations
std::vector<double> shear_mult_wells_;
// 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_;
using Base::vfp_properties_;
// --------- Protected methods ---------
// Need to declare Base members we want to use here.
using Base::wells;
using Base::wellsActive;
using Base::variableState;
using Base::computePressures;
using Base::computeGasPressure;
using Base::applyThresholdPressures;
using Base::fluidViscosity;
using Base::fluidReciprocFVF;
using Base::fluidDensity;
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::computeRelPerm;
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
computeInjectionMobility(const SolutionState& state,
std::vector<ADB>& mob_perfcells);
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 );
void
computeCmax(ReservoirState& state);
ADB
computeMc(const SolutionState& state) const;
const std::vector<PhasePresence>
phaseCondition() const {return this->phaseCondition_;}
/// Computing the water velocity without shear-thinning for the cell faces.
/// The water velocity will be used for shear-thinning calculation.
void computeWaterShearVelocityFaces(const V& transi,
const std::vector<ADB>& phasePressure, const SolutionState& state,
std::vector<double>& water_vel, std::vector<double>& visc_mult);
/// Computing the water velocity without shear-thinning for the well perforations based on the water flux rate.
/// The water velocity will be used for shear-thinning calculation.
void computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw,
std::vector<double>& water_vel_wells, std::vector<double>& visc_mult_wells);
};
/// Need to include concentration in our state variables, otherwise all is as
/// the default blackoil model.
struct BlackoilPolymerSolutionState : public DefaultBlackoilSolutionState
{
explicit BlackoilPolymerSolutionState(const int np)
: DefaultBlackoilSolutionState(np),
concentration( ADB::null())
{
}
ADB concentration;
};
/// Providing types by template specialisation of ModelTraits for BlackoilPolymerModel.
template <class Grid>
struct ModelTraits< BlackoilPolymerModel<Grid> >
{
typedef PolymerBlackoilState ReservoirState;
typedef WellStateFullyImplicitBlackoilPolymer WellState;
typedef BlackoilModelParameters ModelParameters;
typedef BlackoilPolymerSolutionState SolutionState;
};
} // namespace Opm
#include "BlackoilPolymerModel_impl.hpp"
#endif // OPM_BLACKOILPOLYMERMODEL_HEADER_INCLUDED
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
-861
View File
@@ -1,861 +0,0 @@
/*
Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2014, 2015 Statoil ASA.
Copyright 2015 NTNU
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_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
#define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/well_controls.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <cassert>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <limits>
namespace Opm {
namespace detail {
template <class PU>
int polymerPos(const PU& pu)
{
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
int pos = 0;
for (int phase = 0; phase < maxnp; ++phase) {
if (pu.phase_used[phase]) {
pos++;
}
}
return pos;
}
} // namespace detail
template <class Grid>
BlackoilPolymerModel<Grid>::BlackoilPolymerModel(const typename Base::ModelParameters& param,
const Grid& grid,
const BlackoilPropsAdFromDeck& fluid,
const DerivedGeology& geo,
const RockCompressibility* rock_comp_props,
const PolymerPropsAd& polymer_props_ad,
const StandardWells& well_model,
const NewtonIterationBlackoilInterface& linsolver,
std::shared_ptr< const EclipseState > eclipse_state,
std::shared_ptr< const Schedule> schedule,
std::shared_ptr< const SummaryConfig> summary_config,
const bool has_disgas,
const bool has_vapoil,
const bool has_polymer,
const bool has_plyshlog,
const bool has_shrate,
const std::vector<double>& wells_rep_radius,
const std::vector<double>& wells_perf_length,
const std::vector<double>& wells_bore_diameter,
const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, well_model, linsolver, eclipse_state,
schedule, summary_config, has_disgas, has_vapoil, terminal_output),
polymer_props_ad_(polymer_props_ad),
has_polymer_(has_polymer),
has_plyshlog_(has_plyshlog),
has_shrate_(has_shrate),
poly_pos_(detail::polymerPos(fluid.phaseUsage())),
wells_rep_radius_(wells_rep_radius),
wells_perf_length_(wells_perf_length),
wells_bore_diameter_(wells_bore_diameter)
{
if (has_polymer_) {
if (!active_[Water]) {
OPM_THROW(std::logic_error, "Polymer must solved in water!\n");
}
residual_.matbalscale.resize(fluid_.numPhases() + 1, 1.1169); // use the same as the water phase
// If deck has polymer, residual_ should contain polymer equation.
sd_.rq.resize(fluid_.numPhases() + 1);
residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
Base::material_name_.push_back("Polymer");
assert(poly_pos_ == fluid_.numPhases());
}
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::
prepareStep(const SimulatorTimerInterface& timer,
const ReservoirState& reservoir_state,
const WellState& well_state)
{
Base::prepareStep(timer, reservoir_state, well_state);
auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX );
// Initial max concentration of this time step from PolymerBlackoilState.
cmax_ = Eigen::Map<const V>(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_));
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::
afterStep(const SimulatorTimerInterface& /* timer */,
ReservoirState& reservoir_state,
WellState& /* well_state */)
{
computeCmax(reservoir_state);
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::makeConstantState(SolutionState& state) const
{
Base::makeConstantState(state);
state.concentration = ADB::constant(state.concentration.value());
}
template <class Grid>
std::vector<V>
BlackoilPolymerModel<Grid>::variableStateInitials(const ReservoirState& x,
const WellState& xw) const
{
std::vector<V> vars0 = Base::variableStateInitials(x, xw);
assert(int(vars0.size()) == fluid_.numPhases() + 2);
// Initial polymer concentration.
if (has_polymer_) {
const auto& concentration = x.getCellData( x.CONCENTRATION );
assert (not concentration.empty());
const int nc = concentration.size();
const V c = Eigen::Map<const V>(concentration.data() , nc);
// Concentration belongs after other reservoir vars but before well vars.
auto concentration_pos = vars0.begin() + fluid_.numPhases();
assert(concentration_pos == vars0.end() - 2);
vars0.insert(concentration_pos, c);
}
return vars0;
}
template <class Grid>
std::vector<int>
BlackoilPolymerModel<Grid>::variableStateIndices() const
{
std::vector<int> ind = Base::variableStateIndices();
assert(ind.size() == 5);
if (has_polymer_) {
ind.resize(6);
// Concentration belongs after other reservoir vars but before well vars.
ind[Concentration] = fluid_.numPhases();
// Concentration is pushing back the well vars.
++ind[Qs];
++ind[Bhp];
}
return ind;
}
template <class Grid>
typename BlackoilPolymerModel<Grid>::SolutionState
BlackoilPolymerModel<Grid>::variableStateExtractVars(const ReservoirState& x,
const std::vector<int>& indices,
std::vector<ADB>& vars) const
{
SolutionState state = Base::variableStateExtractVars(x, indices, vars);
if (has_polymer_) {
state.concentration = std::move(vars[indices[Concentration]]);
}
return state;
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::computeAccum(const SolutionState& state,
const int aix )
{
Base::computeAccum(state, aix);
// Compute accumulation of polymer equation only if needed.
if (has_polymer_) {
const ADB& press = state.pressure;
const std::vector<ADB>& sat = state.saturation;
const ADB& c = state.concentration;
const ADB pv_mult = poroMult(press); // also computed in Base::computeAccum, could be optimized.
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
// compute polymer properties.
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
const double rho_rock = polymer_props_ad_.rockDensity();
const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], AutoDiffGrid::numCells(grid_));
const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
// Compute polymer accumulation term.
sd_.rq[poly_pos_].accum[aix] = pv_mult * sd_.rq[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
+ pv_mult * rho_rock * (1. - phi) / phi * ads;
}
}
template <class Grid>
void BlackoilPolymerModel<Grid>::computeCmax(ReservoirState& state)
{
auto& max_concentration = state.getCellData( state.CMAX );
const auto& concentration = state.getCellData( state.CONCENTRATION );
std::transform( max_concentration.begin() ,
max_concentration.end() ,
concentration.begin() ,
max_concentration.begin() ,
[](double c_max , double c) { return std::max( c_max , c ); });
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::
assembleMassBalanceEq(const SolutionState& state)
{
// Base::assembleMassBalanceEq(state);
// Compute b_p and the accumulation term b_p*s_p for each phase,
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
// These quantities are stored in sd_.rq[phase].accum[1].
// The corresponding accumulation terms from the start of
// the timestep (b^0_p*s^0_p etc.) were already computed
// on the initial call to assemble() and stored in sd_.rq[phase].accum[0].
computeAccum(state, 1);
// Set up the common parts of the mass balance equations
// for each active phase.
const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
{
const std::vector<ADB> kr = computeRelPerm(state);
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
}
}
if (has_plyshlog_) {
std::vector<double> water_vel;
std::vector<double> visc_mult;
computeWaterShearVelocityFaces(transi, state.canonical_phase_pressures, state, water_vel, visc_mult);
if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) {
// std::cerr << " failed in calculating the shear-multiplier " << std::endl;
OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. ");
}
}
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
const std::vector<PhasePresence>& cond = phaseCondition();
sd_.rq[phaseIdx].mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
sd_.rq[phaseIdx].rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
computeMassFlux(phaseIdx, transi, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
residual_.material_balance_eq[ phaseIdx ] =
pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
+ ops_.div*sd_.rq[phaseIdx].mflux;
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
// Add the extra (flux) terms to the mass balance equations
// From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv)
// The extra terms in the accumulation part of the equation are already handled.
if (active_[ Oil ] && active_[ Gas ]) {
const int po = fluid_.phaseUsage().phase_pos[ Oil ];
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
const UpwindSelector<double> upwindOil(grid_, ops_,
sd_.rq[po].dh.value());
const ADB rs_face = upwindOil.select(state.rs);
const UpwindSelector<double> upwindGas(grid_, ops_,
sd_.rq[pg].dh.value());
const ADB rv_face = upwindGas.select(state.rv);
residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
// OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
}
// Add polymer equation.
if (has_polymer_) {
residual_.material_balance_eq[poly_pos_] = pvdt_ * (sd_.rq[poly_pos_].accum[1] - sd_.rq[poly_pos_].accum[0])
+ ops_.div*sd_.rq[poly_pos_].mflux;
}
if (param_.update_equations_scaling_) {
updateEquationsScaling();
}
}
template <class Grid>
void BlackoilPolymerModel<Grid>::updateEquationsScaling()
{
Base::updateEquationsScaling();
if (has_polymer_) {
const int water_pos = fluid_.phaseUsage().phase_pos[Water];
residual_.matbalscale[poly_pos_] = residual_.matbalscale[water_pos];
}
}
template <class Grid>
void BlackoilPolymerModel<Grid>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
const SolutionState& state,
WellState& xw)
{
Base::addWellContributionToMassBalanceEq(cq_s, state, xw);
// Add well contributions to polymer mass balance equation
if (has_polymer_) {
const ADB mc = computeMc(state);
const int nc = xw.polymerInflow().size();
const V polyin = Eigen::Map<const V>(xw.polymerInflow().data(), nc);
const int nperf = wells().well_connpos[wells().number_of_wells];
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
const V poly_in_perf = subset(polyin, well_cells);
const V poly_mc_perf = subset(mc.value(), well_cells);
const ADB& cq_s_water = cq_s[fluid_.phaseUsage().phase_pos[Water]];
Selector<double> injector_selector(cq_s_water.value());
const V poly_perf = injector_selector.select(poly_in_perf, poly_mc_perf);
const ADB cq_s_poly = cq_s_water * poly_perf;
residual_.material_balance_eq[poly_pos_] -= superset(cq_s_poly, well_cells, nc);
}
}
template <class Grid>
void BlackoilPolymerModel<Grid>::updateState(const V& dx,
ReservoirState& reservoir_state,
WellState& well_state)
{
if (has_polymer_) {
// Extract concentration change.
const int np = fluid_.numPhases();
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const V zero = V::Zero(nc);
const int concentration_start = nc * np;
const V dc = subset(dx, Span(nc, 1, concentration_start));
// Create new dx with the dc part deleted.
V modified_dx = V::Zero(dx.size() - nc);
modified_dx.head(concentration_start) = dx.head(concentration_start);
const int tail_len = dx.size() - concentration_start - nc;
modified_dx.tail(tail_len) = dx.tail(tail_len);
// Call base version.
Base::updateState(modified_dx, reservoir_state, well_state);
{
auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION );
// Update concentration.
const V c_old = Eigen::Map<const V>(concentration.data(), nc, 1);
const V c = (c_old - dc).max(zero);
std::copy(&c[0], &c[0] + nc, concentration.begin());
}
} else {
// Just forward call to base version.
Base::updateState(dx, reservoir_state, well_state);
}
}
template <class Grid>
void
BlackoilPolymerModel<Grid>::computeMassFlux(const int actph ,
const V& transi,
const ADB& kr ,
const ADB& mu ,
const ADB& rho ,
const ADB& phasePressure,
const SolutionState& state)
{
Base::computeMassFlux(actph, transi, kr, mu, rho, phasePressure, state);
// Polymer treatment.
const int canonicalPhaseIdx = canph_[ actph ];
if (canonicalPhaseIdx == Water) {
if (has_polymer_) {
const ADB tr_mult = transMult(state.pressure);
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
// Reduce mobility of water phase by relperm reduction and effective viscosity increase.
sd_.rq[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
// Compute polymer mobility.
const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value());
sd_.rq[poly_pos_].mob = tr_mult * krw_eff * state.concentration * inv_poly_eff_visc;
sd_.rq[poly_pos_].b = sd_.rq[actph].b;
sd_.rq[poly_pos_].dh = sd_.rq[actph].dh;
UpwindSelector<double> upwind(grid_, ops_, sd_.rq[poly_pos_].dh.value());
// Compute polymer flux.
sd_.rq[poly_pos_].mflux = upwind.select(sd_.rq[poly_pos_].b * sd_.rq[poly_pos_].mob) * (transi * sd_.rq[poly_pos_].dh);
// Must recompute water flux since we have to use modified mobilities.
sd_.rq[ actph ].mflux = upwind.select(sd_.rq[actph].b * sd_.rq[actph].mob) * (transi * sd_.rq[actph].dh);
// applying the shear-thinning factors
if (has_plyshlog_) {
V shear_mult_faces_v = Eigen::Map<V>(shear_mult_faces_.data(), shear_mult_faces_.size());
ADB shear_mult_faces_adb = ADB::constant(shear_mult_faces_v);
sd_.rq[poly_pos_].mflux = sd_.rq[poly_pos_].mflux / shear_mult_faces_adb;
sd_.rq[actph].mflux = sd_.rq[actph].mflux / shear_mult_faces_adb;
}
}
}
}
template <class Grid>
SimulatorReport
BlackoilPolymerModel<Grid>::assemble(const ReservoirState& reservoir_state,
WellState& well_state,
const bool initial_assembly)
{
using namespace Opm::AutoDiffGrid;
SimulatorReport report;
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
// updateWellControls(well_state);
wellModel().updateWellControls(well_state);
// Create the primary variables.
SolutionState state = variableState(reservoir_state, well_state);
if (initial_assembly) {
// Create the (constant, derivativeless) initial state.
SolutionState state0 = state;
makeConstantState(state0);
// Compute initial accumulation contributions
// and well connection pressures.
computeAccum(state0, 0);
// computeWellConnectionPressures(state0, well_state);
wellModel().computeWellConnectionPressures(state0, well_state);
}
// OPM_AD_DISKVAL(state.pressure);
// OPM_AD_DISKVAL(state.saturation[0]);
// OPM_AD_DISKVAL(state.saturation[1]);
// OPM_AD_DISKVAL(state.saturation[2]);
// OPM_AD_DISKVAL(state.rs);
// OPM_AD_DISKVAL(state.rv);
// OPM_AD_DISKVAL(state.qs);
// OPM_AD_DISKVAL(state.bhp);
// -------- Mass balance equations --------
assembleMassBalanceEq(state);
// -------- Well equations ----------
if ( ! wellsActive() ) {
return report;
}
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
// updating the the injection mobility related to polymer injection when necessary
// only the mobility of water phase is updated
computeInjectionMobility(state, mob_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
Base::solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
}
V aliveWells;
std::vector<ADB> cq_s;
wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
if (has_plyshlog_) {
std::vector<double> water_vel_wells;
std::vector<double> visc_mult_wells;
const int water_pos = fluid_.phaseUsage().phase_pos[Water];
computeWaterShearVelocityWells(state, well_state, cq_s[water_pos], water_vel_wells, visc_mult_wells);
if ( !polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_) ) {
OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells ");
}
// applying the shear-thinning to the water phase
V shear_mult_wells_v = Eigen::Map<V>(shear_mult_wells_.data(), shear_mult_wells_.size());
ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v);
mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb;
wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
}
wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
wellModel().addWellFluxEq(cq_s, state, residual_);
addWellContributionToMassBalanceEq(cq_s, state, well_state);
wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
report.converged = true;
return report;
}
template <class Grid>
ADB
BlackoilPolymerModel<Grid>::computeMc(const SolutionState& state) const
{
return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration);
}
template<class Grid>
void
BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi,
const std::vector<ADB>& phasePressure, const SolutionState& state,
std::vector<double>& water_vel, std::vector<double>& visc_mult)
{
const int phase = fluid_.phaseUsage().phase_pos[Water]; // water position
const int canonicalPhaseIdx = canph_[phase];
const std::vector<PhasePresence> cond = phaseCondition();
const ADB tr_mult = transMult(state.pressure);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond);
sd_.rq[phase].mob = tr_mult * sd_.rq[phase].kr / mu;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rho = fluidDensity(canonicalPhaseIdx, sd_.rq[phase].b, state.rs, state.rv);
const ADB rhoavg = ops_.caver * rho;
sd_.rq[ phase ].dh = ops_.ngrad * phasePressure[ canonicalPhaseIdx ] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
if (use_threshold_pressure_) {
applyThresholdPressures(sd_.rq[ phase ].dh);
}
const ADB& b = sd_.rq[ phase ].b;
const ADB& mob = sd_.rq[ phase ].mob;
const ADB& dh = sd_.rq[ phase ].dh;
UpwindSelector<double> upwind(grid_, ops_, dh.value());
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
cmax,
sd_.rq[phase].kr);
ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
sd_.rq[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
const V& polymer_conc = state.concentration.value();
V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
V visc_mult_faces = upwind.select(visc_mult_cells);
size_t nface = visc_mult_faces.size();
visc_mult.resize(nface);
std::copy(visc_mult_faces.data(), visc_mult_faces.data() + nface, visc_mult.begin());
sd_.rq[ phase ].mflux = (transi * upwind.select(b * mob)) * dh;
const auto& b_faces_adb = upwind.select(b);
std::vector<double> b_faces(b_faces_adb.value().data(), b_faces_adb.value().data() + b_faces_adb.size());
const auto& internal_faces = ops_.internal_faces;
std::vector<double> internal_face_areas;
internal_face_areas.resize(internal_faces.size());
for (int i = 0; i < internal_faces.size(); ++i) {
internal_face_areas[i] = grid_.face_areas[internal_faces[i]];
}
const ADB phi = Opm::AutoDiffBlock<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
const ADB phiavg_adb = ops_.caver * phi;
std::vector<double> phiavg(phiavg_adb.value().data(), phiavg_adb.value().data() + phiavg_adb.size());
water_vel.resize(nface);
std::copy(sd_.rq[0].mflux.value().data(), sd_.rq[0].mflux.value().data() + nface, water_vel.begin());
for (size_t i = 0; i < nface; ++i) {
water_vel[i] = water_vel[i] / (b_faces[i] * phiavg[i] * internal_face_areas[i]);
}
// for SHRATE keyword treatment
if (has_shrate_) {
// get the upwind water saturation
const Opm::PhaseUsage pu = fluid_.phaseUsage();
const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
const ADB& sw_upwind_adb = upwind.select(sw);
std::vector<double> sw_upwind(sw_upwind_adb.value().data(), sw_upwind_adb.value().data() + sw_upwind_adb.size());
// get the absolute permeability for the faces
std::vector<double> perm;
perm.resize(transi.size());
for (int i = 0; i < transi.size(); ++i) {
perm[i] = transi[i] / internal_faces[i];
}
// get the upwind krw_eff
const ADB& krw_adb = upwind.select(krw_eff);
std::vector<double> krw_upwind(krw_adb.value().data(), krw_adb.value().data() + krw_adb.size());
const double& shrate_const = polymer_props_ad_.shrate();
const double epsilon = std::numeric_limits<double>::epsilon();
// std::cout << "espilon is " << epsilon << std::endl;
// std::cin.ignore();
for (size_t i = 0; i < water_vel.size(); ++i) {
// assuming only when upwinding water saturation is not zero
// there will be non-zero water velocity
if (std::abs(water_vel[i]) < epsilon) {
continue;
}
water_vel[i] *= shrate_const * std::sqrt(phiavg[i] / (perm[i] * sw_upwind[i] * krw_upwind[i]));
}
}
}
template<class Grid>
void
BlackoilPolymerModel<Grid>::computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw,
std::vector<double>& water_vel_wells, std::vector<double>& visc_mult_wells)
{
if( ! wellsActive() ) return ;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
water_vel_wells.resize(cq_sw.size());
std::copy(cq_sw.value().data(), cq_sw.value().data() + cq_sw.size(), water_vel_wells.begin());
const V& polymer_conc = state.concentration.value();
V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
V visc_mult_wells_v = subset(visc_mult_cells, well_cells);
visc_mult_wells.resize(visc_mult_wells_v.size());
std::copy(visc_mult_wells_v.data(), visc_mult_wells_v.data() + visc_mult_wells_v.size(), visc_mult_wells.begin());
const int water_pos = fluid_.phaseUsage().phase_pos[Water];
ADB b_perfcells = subset(sd_.rq[water_pos].b, well_cells);
const ADB& p_perfcells = subset(state.pressure, well_cells);
const V& cdp = wellModel().wellPerforationPressureDiffs();
const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcells - perfpressure;
// selects injection perforations
V selectInjectingPerforations = V::Zero(nperf);
for (int c = 0; c < nperf; ++c) {
if (drawdown.value()[c] < 0) {
selectInjectingPerforations[c] = 1;
}
}
// for the injection wells
for (size_t i = 0; i < well_cells.size(); ++i) {
if (xw.polymerInflow()[well_cells[i]] == 0. && selectInjectingPerforations[i] == 1) { // maybe comparison with epsilon threshold
visc_mult_wells[i] = 1.;
}
if (selectInjectingPerforations[i] == 1) { // maybe comparison with epsilon threshold
if (xw.polymerInflow()[well_cells[i]] == 0.) {
visc_mult_wells[i] = 1.;
} else {
// TODO: not tested for this assumption yet
const double c_perf = state.concentration.value()[well_cells[i]];
visc_mult_wells[i] = polymer_props_ad_.viscMult(c_perf);
}
}
}
const ADB phi = Opm::AutoDiffBlock<double>::constant(Eigen::Map<const V>(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1));
const ADB phi_wells_adb = subset(phi, well_cells);
std::vector<double> phi_wells(phi_wells_adb.value().data(), phi_wells_adb.value().data() + phi_wells_adb.size());
std::vector<double> b_wells(b_perfcells.value().data(), b_perfcells.value().data() + b_perfcells.size());
for (size_t i = 0; i < water_vel_wells.size(); ++i) {
water_vel_wells[i] = b_wells[i] * water_vel_wells[i] / (phi_wells[i] * 2. * M_PI * wells_rep_radius_[i] * wells_perf_length_[i]);
// TODO: CHECK to make sure this formulation is corectly used. Why muliplied by bW.
// Although this formulation works perfectly with the tests compared with other formulations
}
// for SHRATE treatment
if (has_shrate_) {
const double& shrate_const = polymer_props_ad_.shrate();
for (size_t i = 0; i < water_vel_wells.size(); ++i) {
water_vel_wells[i] = shrate_const * water_vel_wells[i] / wells_bore_diameter_[i];
}
}
return;
}
template<class Grid>
void
BlackoilPolymerModel<Grid>::
computeInjectionMobility(const SolutionState& state,
std::vector<ADB>& mob_perfcells)
{
// Calculating the mobility for the polymer injection peforations
if (has_polymer_ && wellModel().localWellsActive()) {
const std::vector<int> well_cells = wellModel().wellOps().well_cells;
const int nperf = well_cells.size();
// Calculating the drawdown to decide the injection perforation
const ADB& p_perfcells = subset(state.pressure, well_cells);
const V& cdp = wellModel().wellPerforationPressureDiffs();
const ADB perfpressure = (wellModel().wellOps().w2p * state.bhp) + cdp;
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcells - perfpressure;
// Polymer concentration in the perforations
const ADB c_perfcells = subset(state.concentration, well_cells);
// Distinguishing the injection perforation from other perforation
// The value is the location in the well_cell array, not the global index
std::vector<int> polymer_inj_cells;
std::vector<int> other_well_cells;
polymer_inj_cells.reserve(nperf);
other_well_cells.reserve(nperf);
for (int c = 0; c < nperf; ++c) {
// TODO: more tests need to be done for this criterion
if (drawdown.value()[c] < 0.0 && c_perfcells.value()[c] > 0.0) {
polymer_inj_cells.push_back(c);
} else {
other_well_cells.push_back(c);
}
}
// there is some polymer injection process going
if ( !polymer_inj_cells.empty() ) {
// the mobility need to be recalculated for the polymer injection cells
const int water_pos = fluid_.phaseUsage().phase_pos[Water];
const ADB mu_perfcells = subset(sd_.rq[water_pos].mu, well_cells);
const ADB c_poly_inj_cells = subset(c_perfcells, polymer_inj_cells);
const ADB mu_poly_inj_cells = subset(mu_perfcells, polymer_inj_cells); // water viscosity
const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(c_poly_inj_cells, mu_poly_inj_cells.value());
const ADB fully_mixing_visc = polymer_props_ad_.viscMult(c_poly_inj_cells) * mu_poly_inj_cells;
// the original mobility for the polymer injection well cells
ADB mob_polymer_inj = subset(mob_perfcells[water_pos], polymer_inj_cells);
const ADB mob_others = subset(mob_perfcells[water_pos], other_well_cells);
mob_polymer_inj = mob_polymer_inj / inv_wat_eff_visc / fully_mixing_visc;
mob_perfcells[water_pos] = superset(mob_polymer_inj, polymer_inj_cells, nperf) +
superset(mob_others, other_well_cells, nperf);
}
}
}
} // namespace Opm
#endif // OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED
opm/polymer/fullyimplicit/PolymerPropsAd.cpp
-391
View File
@@ -1,391 +0,0 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 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 "config.h"
#include <cmath>
#include <vector>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
namespace Opm {
typedef PolymerPropsAd::ADB ADB;
typedef PolymerPropsAd::V V;
double
PolymerPropsAd::rockDensity() const
{
return polymer_props_.rockDensity();
}
double
PolymerPropsAd::deadPoreVol() const
{
return polymer_props_.deadPoreVol();
}
double
PolymerPropsAd::cMax() const
{
return polymer_props_.cMax();
}
const std::vector<double>&
PolymerPropsAd::shearWaterVelocity() const
{
return polymer_props_.shearWaterVelocity();
}
const std::vector<double>&
PolymerPropsAd::shearViscosityReductionFactor() const
{
return polymer_props_.shearViscosityReductionFactor();
}
double
PolymerPropsAd::plyshlogRefConc() const
{
return polymer_props_.plyshlogRefConc();
}
bool
PolymerPropsAd::hasPlyshlogRefSalinity() const
{
return polymer_props_.hasPlyshlogRefSalinity();
}
bool
PolymerPropsAd::hasPlyshlogRefTemp() const
{
return polymer_props_.hasPlyshlogRefTemp();
}
double
PolymerPropsAd::plyshlogRefSalinity() const
{
return polymer_props_.plyshlogRefSalinity();
}
double
PolymerPropsAd::plyshlogRefTemp() const
{
return polymer_props_.plyshlogRefTemp();
}
double PolymerPropsAd::shrate() const
{
if (polymer_props_.hasShrate()) {
return polymer_props_.shrate();
} else {
OPM_THROW(std::logic_error, "the SHRATE keyword is not specified while requested \n");
}
}
double
PolymerPropsAd::viscMult(double c) const
{
return polymer_props_.viscMult(c);
}
V
PolymerPropsAd::viscMult(const V& c) const
{
int nc = c.size();
V visc_mult(nc);
for (int i = 0; i < nc; ++i) {
visc_mult[i] = polymer_props_.viscMult(c[i]);
}
return visc_mult;
}
ADB
PolymerPropsAd::viscMult(const ADB& c) const
{
const int nc = c.size();
V visc_mult(nc);
V dvisc_mult(nc);
for (int i = 0; i < nc; ++i) {
double im = 0, dim = 0;
im = polymer_props_.viscMultWithDer(c.value()(i), &dim);
visc_mult(i) = im;
dvisc_mult(i) = dim;
}
ADB::M dim_diag(dvisc_mult.matrix().asDiagonal());
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dim_diag * c.derivative()[block];
}
return ADB::function(std::move(visc_mult), std::move(jacs));
}
PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props)
: polymer_props_ (polymer_props)
{
}
PolymerPropsAd::~PolymerPropsAd()
{
}
V PolymerPropsAd::effectiveInvWaterVisc(const V& c,
const V& mu_w) const
{
assert(c.size() == mu_w.size());
const int nc = c.size();
V inv_mu_w_eff(nc);
for (int i = 0; i < nc; ++i) {
double im = 0;
polymer_props_.effectiveInvVisc(c(i), mu_w(i), im);
inv_mu_w_eff(i) = im;
}
return inv_mu_w_eff;
}
ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
const V& mu_w) const
{
assert(c.size() == mu_w.size());
const int nc = c.size();
V inv_mu_w_eff(nc);
V dinv_mu_w_eff(nc);
for (int i = 0; i < nc; ++i) {
double im = 0, dim = 0;
polymer_props_.effectiveInvViscWithDer(c.value()(i), mu_w(i), im, dim);
inv_mu_w_eff(i) = im;
dinv_mu_w_eff(i) = dim;
}
ADB::M dim_diag(dinv_mu_w_eff.matrix().asDiagonal());
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dim_diag * c.derivative()[block];
}
return ADB::function(std::move(inv_mu_w_eff), std::move(jacs));
}
ADB PolymerPropsAd::effectiveInvPolymerVisc(const ADB& c, const V& mu_w) const
{
assert(c.size() == mu_w.size());
const int nc = c.size();
V inv_mu_p_eff(nc);
V dinv_mu_p_eff(nc);
for (int i = 0; i < nc; ++i) {
double im = 0;
double dim = 0;
polymer_props_.effectiveInvPolyViscWithDer(c.value()(i), mu_w(i), im, dim);
inv_mu_p_eff(i) = im;
dinv_mu_p_eff(i) = dim;
}
ADB::M dim_diag(dinv_mu_p_eff.matrix().asDiagonal());
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dim_diag * c.derivative()[block];
}
return ADB::function(std::move(inv_mu_p_eff), std::move(jacs));
}
V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const
{
const int nc = c.size();
V mc(nc);
for (int i = 0; i < nc; ++i) {
double m = 0;
polymer_props_.computeMc(c(i), m);
mc(i) = m;
}
return mc;
}
ADB PolymerPropsAd::polymerWaterVelocityRatio(const ADB& c) const
{
const int nc = c.size();
V mc(nc);
V dmc(nc);
for (int i = 0; i < nc; ++i) {
double m = 0;
double dm = 0;
polymer_props_.computeMcWithDer(c.value()(i), m, dm);
mc(i) = m;
dmc(i) = dm;
}
ADB::M dmc_diag(dmc.matrix().asDiagonal());
const int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dmc_diag * c.derivative()[block];
}
return ADB::function(std::move(mc), std::move(jacs));
}
V PolymerPropsAd::adsorption(const V& c, const V& cmax_cells) const
{
const int nc = c.size();
V ads(nc);
for (int i = 0; i < nc; ++i) {
double c_ads = 0;
polymer_props_.adsorption(c(i), cmax_cells(i), c_ads);
ads(i) = c_ads;
}
return ads;
}
ADB PolymerPropsAd::adsorption(const ADB& c, const ADB& cmax_cells) const
{
const int nc = c.value().size();
V ads(nc);
V dads(nc);
for (int i = 0; i < nc; ++i) {
double c_ads = 0;
double dc_ads = 0;
polymer_props_.adsorptionWithDer(c.value()(i), cmax_cells.value()(i), c_ads, dc_ads);
ads(i) = c_ads;
dads(i) = dc_ads;
}
ADB::M dads_diag(dads.matrix().asDiagonal());
int num_blocks = c.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dads_diag * c.derivative()[block];
}
return ADB::function(std::move(ads), std::move(jacs));
}
V
PolymerPropsAd::effectiveRelPerm(const V& c,
const V& cmax_cells,
const V& krw) const
{
const int nc = c.size();
V one = V::Ones(nc);
V ads = adsorption(c, cmax_cells);
double max_ads = polymer_props_.cMaxAds();
double res_factor = polymer_props_.resFactor();
double factor = (res_factor -1.) / max_ads;
V rk = one + factor * ads;
return krw / rk;
}
ADB
PolymerPropsAd::effectiveRelPerm(const ADB& c,
const ADB& cmax_cells,
const ADB& krw) const
{
const int nc = c.value().size();
V one = V::Ones(nc);
ADB ads = adsorption(c, cmax_cells);
V krw_eff = effectiveRelPerm(c.value(), cmax_cells.value(), krw.value());
double max_ads = polymer_props_.cMaxAds();
double res_factor = polymer_props_.resFactor();
double factor = (res_factor - 1.) / max_ads;
ADB rk = one + ads * factor;
return krw / rk;
}
bool
PolymerPropsAd::computeShearMultLog(std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult) const
{
return polymer_props_.computeShearMultLog(water_vel, visc_mult, shear_mult);
}
}// namespace Opm
opm/polymer/fullyimplicit/PolymerPropsAd.hpp
-159
View File
@@ -1,159 +0,0 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
Copyright 2014 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_POLYMERPROPSAD_HEADED_INLCUDED
#define OPM_POLYMERPROPSAD_HEADED_INLCUDED
#include <cmath>
#include <vector>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/polymer/PolymerProperties.hpp>
namespace Opm {
class PolymerPropsAd
{
public:
/// \return Reference rock density.
double rockDensity() const;
/// \return The value of dead pore volume.
double deadPoreVol() const;
/// \return The max concentration injected.
double cMax() const;
/// \ return The water velcoity or shear rate in the PLYSHLOG table
const std::vector<double>& shearWaterVelocity() const;
/// \ return The viscosity reducation factor in the PLYSHLOG table
const std::vector<double>& shearViscosityReductionFactor() const;
/// \ return The reference polymer concentration for PLYSHLOG table
double plyshlogRefConc() const;
/// \ return The flag indicating if reference salinity is specified in PLYSHLOG keyword
bool hasPlyshlogRefSalinity() const;
/// \ return The flag indicating if reference temperature is specified in PLYSHLOG keyword
bool hasPlyshlogRefTemp() const;
/// \ return The reference salinity in PLYSHLOG keyword
double plyshlogRefSalinity() const;
/// \ return The reference temperature in PLYSHLOG keyword
double plyshlogRefTemp() const;
/// \ return the value of SHRATE
double shrate() const;
double viscMult(double c) const; // multipler interpolated from PLYVISC table
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
V viscMult(const V& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n viscosity multiplier from PLVISC table.
ADB viscMult(const ADB& c) const;
/// Constructor wrapping a polymer props.
PolymerPropsAd(const PolymerProperties& polymer_props);
/// Destructor.
~PolymerPropsAd();
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] mu_w Array of n water viscosity values.
/// \return Array of inverse effective water viscosity.
V
effectiveInvWaterVisc(const V& c, const V& mu_w) const;
/// \param[in] c ADB of polymer concentraion.
/// \param[in] mu_w Array of water viscosity value.
/// \return ADB of inverse effective water viscosity.
ADB
effectiveInvWaterVisc(const ADB& c,const V& mu_w) const;
/// \param[in] c ADB of polymer concentraion values.
/// \param[in] mu_w Array of water viscosity values
/// \return ADB of inverse effective polymer viscosity.
ADB
effectiveInvPolymerVisc(const ADB& c, const V& mu_w) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
V
polymerWaterVelocityRatio(const V& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \return Array of n mc values, here mc means m(c) * c.
ADB
polymerWaterVelocityRatio(const ADB& c) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
V
adsorption(const V& c, const V& cmax_cells) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \return Array of n adsorption values.
ADB
adsorption(const ADB& c, const ADB& cmax_cells) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
V
effectiveRelPerm(const V& c, const V& cmax_cells, const V& relperm) const;
/// \param[in] c Array of n polymer concentraion values.
/// \param[in] cmax_cells Array of n polymer concentraion values
/// that the cell experienced.
/// \param[in] relperm Array of n relative water relperm values.
/// \return Array of n adsorption values.
ADB
effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw) const;
/// \param[in] water_vel Array of the n values of water velocity or shear rate.
/// \param[in] visc_mult Array of the n values of the viscosity multiplier from PLYVISC table.
/// \parma[out] shear_mult Array of the n values of calculated shear multiplier with PLYSHLOG keyword.
/// \return TRUE if the calculation of shear multiplier is sucessful,
/// FALSE if the calculation of shear multplier is failed.
bool computeShearMultLog(std::vector<double>& water_vel, std::vector<double>& visc_mult, std::vector<double>& shear_mult) const;
private:
const PolymerProperties& polymer_props_;
};
} //namespace Opm
#endif// OPM_POLYMERPROPSAD_HEADED_INLCUDED
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp
-175
View File
@@ -1,175 +0,0 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 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_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
#include <opm/autodiff/SimulatorBase.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
#include <opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp>
#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/polymer/PolymerInflow.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/RateConverter.hpp>
#include <opm/autodiff/NonlinearSolver.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/core/simulator/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/core/simulator/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 SimulatorFullyImplicitBlackoilPolymer;
class StandardWells;
template<class GridT>
struct SimulatorTraits<SimulatorFullyImplicitBlackoilPolymer<GridT> >
{
typedef WellStateFullyImplicitBlackoilPolymer WellState;
typedef PolymerBlackoilState ReservoirState;
typedef BlackoilOutputWriter OutputWriter;
typedef GridT Grid;
typedef BlackoilPolymerModel<Grid> Model;
typedef NonlinearSolver<Model> Solver;
typedef StandardWells WellModel;
};
/// Class collecting all necessary components for a blackoil simulation with polymer
/// injection.
template <class GridT>
class SimulatorFullyImplicitBlackoilPolymer
: public SimulatorBase<SimulatorFullyImplicitBlackoilPolymer<GridT> >
{
typedef SimulatorFullyImplicitBlackoilPolymer<GridT> ThisType;
typedef SimulatorBase<ThisType> BaseType;
typedef SimulatorTraits<ThisType> Traits;
typedef typename Traits::Solver Solver;
typedef typename Traits::WellModel WellModel;
public:
SimulatorFullyImplicitBlackoilPolymer(const ParameterGroup& param,
const GridT& grid,
DerivedGeology& geo,
BlackoilPropsAdFromDeck& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool disgas,
const bool vapoil,
const bool polymer,
const bool plyshlog,
const bool shrate,
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);
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:
const PolymerPropsAd& polymer_props_;
bool has_polymer_;
// flag for PLYSHLOG keyword
bool has_plyshlog_;
// flag for SHRATE keyword
bool has_shrate_;
std::shared_ptr< Deck > deck_;
std::vector<double> wells_rep_radius_;
std::vector<double> wells_perf_length_;
std::vector<double> wells_bore_diameter_;
// generate the mapping from Cartesian grid cells to global compressed cells,
// copied from opm-core, to be used in function computeRepRadiusPerfLength()
static void
setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed);
// calculate the representative radius and length for for well peforations
// and store the wellbore diameters
// it will be used in the shear-thinning calcluation only.
void
computeRepRadiusPerfLength(const Schedule& schedule,
const size_t timeStep,
const GridT& grid,
std::vector<double>& wells_rep_radius,
std::vector<double>& wells_perf_length,
std::vector<double>& wells_bore_diameter);
};
} // namespace Opm
#include "SimulatorFullyImplicitBlackoilPolymer_impl.hpp"
#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
-275
View File
@@ -1,275 +0,0 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2014 IRIS AS
Copyright 2014 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/>.
*/
namespace Opm
{
template <class GridT>
SimulatorFullyImplicitBlackoilPolymer<GridT>::
SimulatorFullyImplicitBlackoilPolymer(const ParameterGroup& param,
const GridT& grid,
DerivedGeology& geo,
BlackoilPropsAdFromDeck& props,
const PolymerPropsAd& polymer_props,
const RockCompressibility* rock_comp_props,
NewtonIterationBlackoilInterface& linsolver,
const double* gravity,
const bool has_disgas,
const bool has_vapoil,
const bool has_polymer,
const bool has_plyshlog,
const bool has_shrate,
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)
: 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>())
, polymer_props_(polymer_props)
, has_polymer_(has_polymer)
, has_plyshlog_(has_plyshlog)
, has_shrate_(has_shrate)
, deck_(deck)
{
}
template <class GridT>
auto SimulatorFullyImplicitBlackoilPolymer<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_,
polymer_props_,
well_model,
BaseType::solver_,
BaseType::eclipse_state_,
BaseType::schedule_,
BaseType::summary_config_,
BaseType::has_disgas_,
BaseType::has_vapoil_,
has_polymer_,
has_plyshlog_,
has_shrate_,
wells_rep_radius_,
wells_perf_length_,
wells_bore_diameter_,
BaseType::terminal_output_));
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 SimulatorFullyImplicitBlackoilPolymer<GridT>::
handleAdditionalWellInflow(SimulatorTimer& timer,
WellsManager& wells_manager,
typename BaseType::WellState& well_state,
const Wells* wells)
{
// compute polymer inflow
std::unique_ptr<PolymerInflowInterface> polymer_inflow_ptr;
if (deck_->hasKeyword("WPOLYMER")) {
if (wells_manager.c_wells() == 0) {
OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(*BaseType::schedule_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
} else {
OPM_MESSAGE("Warning: simulating with no WPOLYMER in deck (no polymer will be injected).");
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));
}
std::vector<double> polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_));
polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
polymer_inflow_c);
well_state.polymerInflow() = polymer_inflow_c;
if (has_plyshlog_) {
computeRepRadiusPerfLength(*BaseType::schedule_, timer.currentStepNum(), BaseType::grid_, wells_rep_radius_, wells_perf_length_, wells_bore_diameter_);
}
}
template <class GridT>
void SimulatorFullyImplicitBlackoilPolymer<GridT>::
setupCompressedToCartesian(const int* global_cell, int number_of_cells,
std::map<int,int>& cartesian_to_compressed )
{
if (global_cell) {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
}
}
else {
for (int i = 0; i < number_of_cells; ++i) {
cartesian_to_compressed.insert(std::make_pair(i, i));
}
}
}
template <class GridT>
void SimulatorFullyImplicitBlackoilPolymer<GridT>::
computeRepRadiusPerfLength(const Opm::Schedule& schedule,
const size_t timeStep,
const GridT& grid,
std::vector<double>& wells_rep_radius,
std::vector<double>& wells_perf_length,
std::vector<double>& wells_bore_diameter)
{
// TODO, the function does not work for parallel running
// to be fixed later.
int number_of_cells = Opm::UgGridHelpers::numCells(grid);
const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
auto cell_to_faces = Opm::UgGridHelpers::cell2Faces(grid);
auto begin_face_centroids = Opm::UgGridHelpers::beginFaceCentroids(grid);
if (schedule.numWells() == 0) {
OPM_MESSAGE("No wells specified in Schedule section, "
"initializing no wells");
return;
}
const size_t n_perf = wells_rep_radius.size();
wells_rep_radius.clear();
wells_perf_length.clear();
wells_bore_diameter.clear();
wells_rep_radius.reserve(n_perf);
wells_perf_length.reserve(n_perf);
wells_bore_diameter.reserve(n_perf);
std::map<int,int> cartesian_to_compressed;
setupCompressedToCartesian(global_cell, number_of_cells,
cartesian_to_compressed);
auto wells = schedule.getWells(timeStep);
int well_index = 0;
for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
const auto* well = (*wellIter);
if (well->getStatus(timeStep) == WellCommon::SHUT) {
continue;
}
{ // COMPDAT handling
const auto& completionSet = well->getCompletions(timeStep);
for (size_t c=0; c<completionSet.size(); c++) {
const auto& completion = completionSet.get(c);
if (completion.getState() == WellCompletion::OPEN) {
int i = completion.getI();
int j = completion.getJ();
int k = completion.getK();
const int* cpgdim = cart_dims;
int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' '
<< k << " not found in grid (well = " << well->name() << ')');
}
int cell = cgit->second;
{
double radius = 0.5*completion.getDiameter();
if (radius <= 0.0) {
radius = 0.5*unit::feet;
OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
}
const std::array<double, 3> cubical =
WellsManagerDetail::getCubeDim<3>(cell_to_faces, begin_face_centroids, cell);
WellCompletion::DirectionEnum direction = completion.getDirection();
double re; // area equivalent radius of the grid block
double perf_length; // the length of the well perforation
switch (direction) {
case Opm::WellCompletion::DirectionEnum::X:
re = std::sqrt(cubical[1] * cubical[2] / M_PI);
perf_length = cubical[0];
break;
case Opm::WellCompletion::DirectionEnum::Y:
re = std::sqrt(cubical[0] * cubical[2] / M_PI);
perf_length = cubical[1];
break;
case Opm::WellCompletion::DirectionEnum::Z:
re = std::sqrt(cubical[0] * cubical[1] / M_PI);
perf_length = cubical[2];
break;
default:
OPM_THROW(std::runtime_error, " Dirtecion of well is not supported ");
}
double repR = std::sqrt(re * radius);
wells_rep_radius.push_back(repR);
wells_perf_length.push_back(perf_length);
wells_bore_diameter.push_back(2. * radius);
}
} else {
if (completion.getState() != WellCompletion::SHUT) {
OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled");
}
}
}
}
well_index++;
}
}
} // namespace Opm
opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp
-39
View File
@@ -1,39 +0,0 @@
/*
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_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
#define OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
namespace Opm
{
class WellStateFullyImplicitBlackoilPolymer : public WellStateFullyImplicitBlackoil
{
public:
std::vector<double>& polymerInflow() { return polymer_inflow_; }
const std::vector<double>& polymerInflow() const { return polymer_inflow_; }
private:
std::vector<double> polymer_inflow_;
};
} // namespace Opm
#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
tests/update_reference_data.sh
-6
View File
@@ -91,12 +91,6 @@ for test_name in ${tests}; do
if grep -q "polymer2d" <<< $test_name if grep -q "polymer2d" <<< $test_name
then then
copyToReferenceDir \
$configuration/build-opm-simulators/tests/results/flow_polymer+polymer_simple2D/ \
$OPM_DATA_ROOT/polymer_simple2D/opm-simulation-reference/flow_polymer \
2D_THREEPHASE_POLY_HETER \
EGRID INIT SMSPEC UNRST UNSMRY
copyToReferenceDir \ copyToReferenceDir \
$configuration/build-opm-simulators/tests/results/flow+polymer_simple2D/ \ $configuration/build-opm-simulators/tests/results/flow+polymer_simple2D/ \
$OPM_DATA_ROOT/polymer_simple2D/opm-simulation-reference/flow \ $OPM_DATA_ROOT/polymer_simple2D/opm-simulation-reference/flow \