From 0cf6699591862cb3c5b4a2a20f127a02e7e17d77 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 15 Jun 2017 11:41:03 +0200 Subject: [PATCH] adding StandardWell class copied from the old implementation, which is the starting point for the new refactoring --- CMakeLists_files.cmake | 2 + opm/autodiff/StandardWell.cpp | 292 ++++++++++++++++++++++++++++++++++ opm/autodiff/StandardWell.hpp | 126 +++++++++++++++ 3 files changed, 420 insertions(+) create mode 100644 opm/autodiff/StandardWell.cpp create mode 100644 opm/autodiff/StandardWell.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index c31f4fdd9..85bb58a4d 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -50,6 +50,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/WellMultiSegment.cpp opm/autodiff/MultisegmentWells.cpp opm/autodiff/WellInterface.cpp + opm/autodiff/StandardWell.hpp opm/autodiff/MissingFeatures.cpp opm/polymer/PolymerState.cpp opm/polymer/PolymerBlackoilState.cpp @@ -241,6 +242,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/StandardWells.hpp opm/autodiff/StandardWells_impl.hpp opm/autodiff/WellInterface.hpp + opm/autodiff/StandardWell.cpp opm/autodiff/StandardWellsDense.hpp opm/autodiff/StandardWellsSolvent.hpp opm/autodiff/StandardWellsSolvent_impl.hpp diff --git a/opm/autodiff/StandardWell.cpp b/opm/autodiff/StandardWell.cpp new file mode 100644 index 000000000..5a8568b2b --- /dev/null +++ b/opm/autodiff/StandardWell.cpp @@ -0,0 +1,292 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 2017 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 . +*/ + +#include "config.h" + +#include + + + +namespace Opm +{ + StandardWell:: + StandardWell(const Well* well, const size_t time_step, const Wells* wells) + : WellInterface(well, time_step, wells) + , perf_densities_(number_of_perforations_) + , perf_pressure_diffs_(number_of_perforations_) + , well_variables_(blocksize) // the number of the primary variables + { + dune_B_.setBuildMode( Mat::row_wise ); + dune_C_.setBuildMode( Mat::row_wise ); + inv_dune_D_.setBuildMode( Mat::row_wise ); + } + + + + + + const std::vector& + StandardWell:: + perfDensities() const + { + return perf_densities_; + } + + + + + + std::vector& + StandardWell:: + perfDensities() + { + return perf_densities_; + } + + + + + + const std::vector& + StandardWell:: + perfPressureDiffs() const + { + return perf_pressure_diffs_; + } + + + + + + std::vector& + StandardWell:: + perfPressureDiffs() + { + return perf_pressure_diffs_; + } + + + + + + void + StandardWell:: + assembleWellEq(Simulator& ebos_simulator, + const double dt, + WellState& well_state, + bool only_wells) + { + } + + + + + + void StandardWell:: + setWellVariables(const WellState& well_state) + { + const int np = number_of_phases_; + const int nw = well_state.bhp().size(); + // TODO: it should be the number of primary variables + // TODO: this is from the old version of StandardWellsDense, it is a coincidence, 3 phases and 3 primary variables + // TODO: it needs to be careful. + // TODO: the following code has to be rewritten later for correctness purpose. + for (int phase = 0; phase < np; ++phase) { + well_variables_[phase] = 0.0; + well_variables_[phase].setValue(well_state.wellSolutions()[index_of_well_ + nw * phase]); + well_variables_[phase].setDerivative(blocksize + phase, 1.0); + } + } + + + + + + StandardWell::EvalWell + StandardWell:: + getBhp() const + { + const WellControls* wc = well_controls_; + if (well_controls_get_current_type(wc) == BHP) { + EvalWell bhp = 0.0; + const double target_rate = well_controls_get_current_target(wc); + bhp.setValue(target_rate); + return bhp; + } else if (well_controls_get_current_type(wc) == THP) { + const int control = well_controls_get_current(wc); + const double thp = well_controls_get_current_target(wc); + const double alq = well_controls_iget_alq(wc, control); + const int table_id = well_controls_iget_vfp(wc, control); + EvalWell aqua = 0.0; + EvalWell liquid = 0.0; + EvalWell vapour = 0.0; + EvalWell bhp = 0.0; + double vfp_ref_depth = 0.0; + + const Opm::PhaseUsage& pu = phaseUsage(); + + if (active()[ Water ]) { + aqua = getQs(pu.phase_pos[ Water]); + } + if (active()[ Oil ]) { + liquid = getQs(pu.phase_pos[ Oil ]); + } + if (active()[ Gas ]) { + vapour = getQs(pu.phase_pos[ Gas ]); + } + if (wellType() == INJECTOR) { + bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp); + vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + } else { + bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq); + vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + } + + // pick the density in the top layer + const double rho = perf_densities_[0]; + // TODO: not sure whether it is always correct + const double well_ref_depth = perf_depth_[0]; + // const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_); + const double dp = wellhelpers::computeHydrostaticCorrection(well_ref_depth, vfp_ref_depth, rho, gravity_); + bhp -= dp; + return bhp; + } + + return well_variables_[XvarWell]; + } + + + + + + StandardWell::EvalWell + StandardWell:: + getQs(const int phase) const + { + EvalWell qs = 0.0; + + const WellControls* wc = well_controls_; + const int np = number_of_phases_; + + // the target from the well controls + const double target = well_controls_get_current_target(wc); + + // TODO: the formulation for the injectors decides it only work with single phase + // surface rate injection control. Improvement will be required. + // Injectors + if (wellType() == INJECTOR) { + // TODO: we should not rely on comp_frac anymore, it should depend on the values in distr + const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx]; + if (comp_frac == 0.0) { + return qs; + } + + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) { + return wellVariables_[XvarWell]; + } + + // rate control + // TODO: if it is reservoir volume rate, it should be wrong here + qs.setValue(target); + return qs; + } + + + // Producers + if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) { + return wellVariables_[XvarWell] * wellVolumeFractionScaled(phase); + } + + } + + + + + + + StandardWell::EvalWell + StandardWell:: + wellVolumeFractionScaled(const int phase) const + { + // TODO: we should be able to set the g for the well based on the control type + // instead of using explicit code for g all the times + const WellControls* wc = well_controls_; + if (well_controls_get_current_type(wc) == RESERVOIR_RATE) { + const double* distr = well_controls_get_current_distr(wc); + if (distr[phase] > 0.) { + return wellVolumeFraction(phase) / distr[phase]; + } else { + // TODO: not sure why return EvalWell(0.) causing problem here + // Probably due to the wrong Jacobians. + return wellVolumeFraction(phase); + } + } + std::vector g = {1,1,0.01}; + return (wellVolumeFraction(phase) / g[phase]); + } + + + + + + StandardWell::EvalWell + StandardWell:: + wellVolumeFraction(const int phase) const + { + if (phase == Water) { + return wellVariables_[WFrac]; + } + + if (phase == Gas) { + return wellVariables_[GFrac]; + } + + // Oil fraction + EvalWell well_fraction = 1.0; + if (active_[Water]) { + well_fraction -= wellVariables_[WFrac]; + } + + if (active_[Gas]) { + well_fraction -= wellVariables_[GFrac]; + } + return well_fraction; + } + + + + + + StandardWell::EvalWell + StandardWell:: + wellSurfaceVolumeFraction(const int phase) const + { + EvalWell sum_volume_fraction_scaled = 0.; + const int np = number_of_phases_; + for (int p = 0; p < np; ++p) { + sum_volume_fraction_scaled += wellVolumeFractionScaled(p); + } + + assert(sum_volume_fraction_scaled.value() != 0.); + + return wellVolumeFractionScaled(phase) / sum_volume_fraction_scaled; + } + +} diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp new file mode 100644 index 000000000..f2da9043a --- /dev/null +++ b/opm/autodiff/StandardWell.hpp @@ -0,0 +1,126 @@ +/* + Copyright 2017 SINTEF ICT, Applied Mathematics. + Copyright 2017 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 . +*/ + + +#ifndef OPM_STANDARDWELL_HEADER_INCLUDED +#define OPM_STANDARDWELL_HEADER_INCLUDED + + +#include "config.h" + +#include + +#include +#include +#include + +#include +#include + +namespace Opm +{ + + class StandardWell: public WellInterface + { + + public: + using WellInterface::Simulator; + using WellInterface::WellState; + + // the positions of the primary variables for StandardWell + // there are three primary variables, the second and the third ones are F_w and F_g + // the first one can be total rate (G_t) or bhp, based on the control + enum WellVariablePositions { + XvarWell = 0, + WFrac = 1, + GFrac = 2 + }; + + // for now, using the matrix and block version in StandardWellsDense. + // TODO: for bettern generality, it should contain blocksize_field and blocksize_well. + // They are allowed to be different and it will create four types of matrix blocks and two types of + // vector blocks. + + /* const static int blocksize = 3; + typedef double Scalar; + typedef Dune::FieldVector VectorBlockType; + typedef Dune::FieldMatrix MatrixBlockType; + typedef Dune::BCRSMatrix Mat; + typedef Dune::BlockVector BVector; + typedef DenseAd::Evaluation EvalWell; */ + /* using WellInterface::EvalWell; + using WellInterface::BVector; + using WellInterface::Mat; + using WellInterface::MatrixBlockType; + using WellInterface::VectorBlockType; */ + + StandardWell(const Well* well, const size_t time_step, const Wells* wells); + + /// the densities of the fluid in each perforation + virtual const std::vector& perfDensities() const; + virtual std::vector& perfDensities(); + + /// the pressure difference between different perforations + virtual const std::vector& perfPressureDiffs() const; + virtual std::vector& perfPressureDiffs(); + + virtual void assembleWellEq(Simulator& ebos_simulator, + const double dt, + WellState& well_state, + bool only_wells); + + virtual void setWellVariables(const WellState& well_state); + + EvalWell wellVolumeFractionScaled(const int phase) const; + + EvalWell wellVolumeFraction(const int phase) const; + + EvalWell wellSurfaceVolumeFraction(const int phase) const; + + protected: + // densities of the fluid in each perforation + std::vector perf_densities_; + // pressure drop between different perforations + std::vector perf_pressure_diffs_; + + // TODO: probably, they should be moved to the WellInterface, when + // we decide the template paramters. + // two off-diagonal matrices + Mat dune_B_; + Mat dune_C_; + // diagonal matrix for the well + Mat inv_dune_D_; + + BVector res_well_; + + std::vector well_variables_; + + // TODO: this function should be moved to the base class. + // while it faces chanllenges for MSWell later, since the calculation of bhp + // based on THP is never implemented for MSWell yet. + EvalWell getBhp() const; + + // TODO: it is also possible to be moved to the base class. + EvalWell getQs(const int phase) const; + }; + +} + +#endif // OPM_STANDARDWELL_HEADER_INCLUDED