From 043c76a5d3c4cebcbced8b20ad4ad28731d39594 Mon Sep 17 00:00:00 2001 From: Rohith Nair Date: Tue, 29 May 2018 16:35:26 +0200 Subject: [PATCH 01/31] Add test for carter tracy aquifer --- compareECLFiles.cmake | 7 +++++++ tests/update_reference_data.sh | 11 ++++++++++- 2 files changed, 17 insertions(+), 1 deletion(-) mode change 100644 => 100755 compareECLFiles.cmake diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake old mode 100644 new mode 100755 index b01e3c3fd..13b230119 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -195,6 +195,13 @@ add_test_compareECLFiles(CASENAME spe1_thermal REL_TOL ${rel_tol} DIR spe1) +add_test_compareECLFiles(CASENAME ctaquifer_2d_oilwater + FILENAME 2D_OW_CTAQUIFER + SIMULATOR flow + ABS_TOL ${abs_tol} + REL_TOL ${rel_tol} + DIR spe1) + foreach(SIM flow flow_legacy) add_test_compareECLFiles(CASENAME spe3 FILENAME SPE3CASE1 diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh index ab88a21c9..50f30707f 100755 --- a/tests/update_reference_data.sh +++ b/tests/update_reference_data.sh @@ -20,7 +20,7 @@ copyToReferenceDir () { } tests=${@:2} -test -z "$tests" && tests="spe11 spe12 spe12p spe1oilgas spe1nowells spe1thermal spe3 spe5 spe9 norne_init msw_2d_h msw_3d_hfa polymer2d spe9group polymer_oilwater" +test -z "$tests" && tests="spe11 spe12 spe12p spe1oilgas spe1nowells spe1thermal ctaquifer_2d_oilwater spe3 spe5 spe9 norne_init msw_2d_h msw_3d_hfa polymer2d spe9group polymer_oilwater" if grep -q -i "norne " <<< $ghprbCommentBody then if test -d $WORKSPACE/deps/opm-tests/norne/flow @@ -98,6 +98,15 @@ for test_name in ${tests}; do EGRID INIT SMSPEC UNRST UNSMRY fi + if grep -q "ctaquifer_2d_oilwater" <<< $test_name + then + copyToReferenceDir \ + $configuration/build-opm-simulators/tests/results/flow+ctaquifer_2d_oilwater/ \ + $OPM_TESTS_ROOT/spe1/opm-simulation-reference/flow \ + 2D_OW_CTAQUIFER \ + EGRID INIT SMSPEC UNRST UNSMRY + fi + if grep -q "msw_2d_h" <<< $test_name then copyToReferenceDir \ From 196b5b29e8b55eeff33062b658de28bc2b129de6 Mon Sep 17 00:00:00 2001 From: Rohith Nair Date: Tue, 29 May 2018 16:40:40 +0200 Subject: [PATCH 02/31] edit --- tests/update_reference_data.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh index 50f30707f..8d384e50d 100755 --- a/tests/update_reference_data.sh +++ b/tests/update_reference_data.sh @@ -102,7 +102,7 @@ for test_name in ${tests}; do then copyToReferenceDir \ $configuration/build-opm-simulators/tests/results/flow+ctaquifer_2d_oilwater/ \ - $OPM_TESTS_ROOT/spe1/opm-simulation-reference/flow \ + $OPM_TESTS_ROOT/aquifer-oilwater/opm-simulation-reference/flow \ 2D_OW_CTAQUIFER \ EGRID INIT SMSPEC UNRST UNSMRY fi From 5c7cc91a39e4d6bb1903b5e6517cfcab9c24e230 Mon Sep 17 00:00:00 2001 From: Rohith Nair Date: Tue, 29 May 2018 16:44:28 +0200 Subject: [PATCH 03/31] edit --- compareECLFiles.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake index 13b230119..c4175c03c 100755 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -200,7 +200,7 @@ add_test_compareECLFiles(CASENAME ctaquifer_2d_oilwater SIMULATOR flow ABS_TOL ${abs_tol} REL_TOL ${rel_tol} - DIR spe1) + DIR aquifer-oilwater) foreach(SIM flow flow_legacy) add_test_compareECLFiles(CASENAME spe3 From 05a02b9d68ef2e35b431936b9433763ddcb740c3 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 2 Nov 2017 12:00:02 +0100 Subject: [PATCH 04/31] Added build directory in the .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 585d691bc..3f98f891f 100644 --- a/.gitignore +++ b/.gitignore @@ -32,3 +32,5 @@ test_vec # emacs directory setting: .dir-locals.el + +build From 10e896fa6bcc0f6ae71c2b28ad256656c4c96e8f Mon Sep 17 00:00:00 2001 From: kel85uk Date: Mon, 4 Dec 2017 22:26:53 +0100 Subject: [PATCH 05/31] First addition of the class BlackoilAquiferModel. --- opm/autodiff/AquiferCarterTracy.hpp | 402 ++++++++++++++++++ opm/autodiff/BlackoilAquiferModel.hpp | 184 ++++++++ opm/autodiff/BlackoilAquiferModel_impl.hpp | 287 +++++++++++++ opm/autodiff/BlackoilModelEbos.hpp | 60 ++- opm/autodiff/BlackoilModelParameters.cpp | 2 + opm/autodiff/BlackoilModelParameters.hpp | 3 + opm/autodiff/BlackoilTransportModel.hpp | 1 + .../SimulatorFullyImplicitBlackoilEbos.hpp | 17 +- 8 files changed, 953 insertions(+), 3 deletions(-) create mode 100644 opm/autodiff/AquiferCarterTracy.hpp create mode 100644 opm/autodiff/BlackoilAquiferModel.hpp create mode 100644 opm/autodiff/BlackoilAquiferModel_impl.hpp diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp new file mode 100644 index 000000000..b6f9cfa17 --- /dev/null +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -0,0 +1,402 @@ +/* + Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface + 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_AQUIFERCT_HEADER_INCLUDED +#define OPM_AQUIFERCT_HEADER_INCLUDED + +#include +#include +#include +#include +#include + + +#include +#include + +#include +#include +#include +#include +#include + +namespace Opm +{ + + template + class AquiferCarterTracy + { + + public: + typedef BlackoilModelParameters ModelParameters; + + static const int Water = BlackoilPhases::Aqua; + static const int Oil = BlackoilPhases::Liquid; + static const int Gas = BlackoilPhases::Vapour; + + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + static const int numEq = BlackoilIndices::numEq; + typedef double Scalar; + + typedef DenseAd::Evaluation Eval; + + typedef Ewoms::BlackOilPolymerModule PolymerModule; + + static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); + static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); + static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; + static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; + + + AquiferCarterTracy(const std::vector& cell_id) + : phi_aq_ (1.0), // + C_t_ (1.0), // + r_o_ (1.0), // + k_a_ (1.0), // + c1_ (1.0), + h_ (1.0), // + theta_ (1.0), // + c2_ (1.0), // + d0_ (1.0), + cell_idx_ (cell_id) + { + mu_w_ = 1e-3; + aqutab_td_.push_back(1.0); + aqutab_pi_.push_back(1.0); + aquiferID_ = 1; + inftableID_ = 1; + pvttableID_ = 1; + init_quantities(); + } + + explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const AquiferCT::AQUANCON_data& aquanconParams, + const int numComponents, const Scalar gravity ) + : phi_aq_ (params.phi_aq), // + C_t_ (params.C_t), // + r_o_ (params.r_o), // + k_a_ (params.k_a), // + c1_ (params.c1), + h_ (params.h), // + theta_ (params.theta), // + c2_ (params.c2), // + d0_ (params.d0), + aqutab_td_ (params.td), + aqutab_pi_ (params.pi), + aquiferID_ (params.aquiferID), + inftableID_ (params.inftableID), + pvttableID_ (params.pvttableID), + cell_idx_ (params.cell_id), + num_components_ (numComponents), + gravity_ (gravity) + { + mu_w_ = 1e-3; + init_quantities(aquanconParams); + } + + inline const PhaseUsage& + phaseUsage() const + { + assert(phase_usage_); + + return *phase_usage_; + } + + inline int + flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage(); + if (pu.phase_pos[Water] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + if (pu.phase_pos[Oil] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (pu.phase_pos[Gas] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + + // for other phases return the index + return phaseIdx; + } + + inline int + flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage(); + if (pu.phase_pos[Water] == phaseIdx) { + return FluidSystem::waterPhaseIdx; + } + if (pu.phase_pos[Oil] == phaseIdx) { + return FluidSystem::oilPhaseIdx; + } + if (pu.phase_pos[Gas] == phaseIdx) { + return FluidSystem::gasPhaseIdx; + } + + assert(phaseIdx < 3); + // for other phases return the index + return phaseIdx; + } + + inline void calculateExplicitQuantities(const Simulator& ebosSimulator) + { + std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; + } + + inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + { + std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; + // resAqui_ = 0.0; + dt_ = timer.currentStepLength(); + auto& ebosJac = ebosSimulator.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + // TODO: it probably can be static member for StandardWell + const double volume = 0.002831684659200; // 0.1 cu ft; + + auto cellID = cell_idx_.begin(); + size_t idx; + for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) + { + Eval qinflow = 0.0; + // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to + // IntensiveQuantities of that particular cell_id + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + // This is the pressure at td + dt + get_current_Pressure_cell(pressure_current_,idx,intQuants); + get_current_density_cell(rhow_,idx,intQuants); + calculate_inflow_rate(idx, timer); + qinflow = Qai_[idx]; + ebosResid[*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)] -= qinflow.value(); + + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) + { + // also need to consider the efficiency factor when manipulating the jacobians. + ebosJac[*cellID][*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)][pvIdx] -= qinflow.derivative(pvIdx); + } + std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ + // << " -> P_wat[t+dt] = " << pressure_current_[idx] << std::endl + << " Qai[t+dt] = " << Qai_[idx] << std::endl; + } + } + + inline void before_time_step(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + { + auto cellID = cell_idx_.begin(); + size_t idx; + for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) + { + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + get_current_Pressure_cell(pressure_previous_ ,idx,intQuants); + } + } + + inline void after_time_step() + { + for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) + { + W_flux_ += (*Qai); + } + std::cout << "Aquifer # " << aquiferID_ << ": My cumulative flux = " << W_flux_ << std::endl; + } + + /* Made into public for testing only!!!!!!. Must be protected */ + inline const Scalar time_constant() const + { + Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); + return Tc; + } + + /* Made into public for testing only!!!!!!. Must be protected */ + inline const Scalar aquifer_influx_constant() const + { + Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; + return beta; + } + + // This is another hack to get the face area only for SPE1. + // Ideally it should be a map which given a cell_id, it returns the area fraction + inline const double area_fraction(const int i) + { + return 1000.0*20.0*0.092903/(1000.0*1000.0*0.092903*2 + 1000.0*20.0*0.092903*4); + } + + inline void print_private_members() const + { + std::cout << "Aquifer CT #" << aquiferID_ << std::endl; + auto ita = aqutab_td_.cbegin(); + auto f_lambda = [&ita] (double i) {std::cout << *ita++ << " " << i << std::endl;}; + std::for_each( aqutab_pi_.cbegin(), aqutab_pi_.cend(), f_lambda ); + + for (auto i = coeff_.begin(); i != coeff_.end(); ++i ) + { + std::cout << "Coeff = " << *i << std::endl; + } + } + + /* Made into public for testing only!!!!!!. Must be protected */ + inline const std::vector cell_id() const + { + return cell_idx_; + } + + inline const int& aquiferID() const + { + return aquiferID_; + } + + + + protected: + const PhaseUsage* phase_usage_; + + + // Aquifer ID, and other IDs + int aquiferID_, inftableID_, pvttableID_; + int num_components_; + + // Grid variables + std::vector cell_idx_; + + // Quantities at each grid id + std::vector cell_depth_; + std::vector pressure_previous_; + std::vector pressure_current_; + std::vector Qai_; + std::vector rhow_; + + // Variables constants + Scalar mu_w_ , //water viscosity + phi_aq_ , //aquifer porosity + d0_, // aquifer datum depth + C_t_ , //total compressibility + r_o_ , //aquifer inner radius + k_a_ , //aquifer permeability + c1_, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) + h_ , //aquifer thickness + theta_ , //angle subtended by the aquifer boundary + c2_ ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). + + // Variables for influence table + std::vector aqutab_td_, aqutab_pi_; + + // Cumulative flux + Scalar W_flux_, dt_, pa0_, gravity_; + + // Also return the polynomial fit + std::vector coeff_; + + + // We fit the tabular data using a polynomial fit + // Modified from Copyright (C) 2014 Clifford Wolf + // http://svn.clifford.at/handicraft/2014/polyfit/polyfit.cc + inline void polynomial_fit( const std::vector &X, const std::vector &y, + std::vector &coeff, int order, bool bias) const + { + int colNum = (bias)? order + 1 : order; + Eigen::MatrixXd A(X.size(), colNum); + Eigen::VectorXd y_mapped = Eigen::VectorXd::Map(&y.front(), y.size()); + Eigen::VectorXd result; + + assert(X.size() == y.size()); + assert(X.size() >= colNum); + + // create matrix + for (size_t i = 0; i < X.size(); i++) + for (size_t j = 0; j < colNum; j++) + A(i, j) = (bias)? pow(X.at(i), j) : pow(X.at(i), j+1); + + // solve for linear least squares fit + result = A.householderQr().solve(y_mapped); + + coeff.resize(colNum); + for (size_t i = 0; i < colNum; i++) + coeff[i] = result[i]; + } + + inline void init_quantities(const AquiferCT::AQUANCON_data& aquanconParams) + { + W_flux_ = 0.; + // pa0_ is the initial aquifer water pressure. Must be calculated from equilibrium if left default, + // or we get the information from the deck Hacked to make it at 45e6 Pa + pa0_ = 45e6; + + pressure_previous_.resize(cell_idx_.size(), 0.); + pressure_current_.resize(cell_idx_.size(), 0.); + // We hack the cell depth values for now. We can actually get it from elementcontext pos + cell_depth_.resize(cell_idx_.size(), d0_); + rhow_.resize(cell_idx_.size(), 998.0); + Qai_.resize(cell_idx_.size(), 0.); + + polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 2, true); + } + + inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + pressure_water[idx] = fs.pressure(FluidSystem::waterPhaseIdx).value(); + } + + inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + rho_water[idx] = fs.density(FluidSystem::waterPhaseIdx).value(); + } + + inline Scalar dpai(int idx) + { + Scalar dp = pa0_ - rhow_[idx]*gravity_*(cell_depth_[idx] - d0_) - pressure_previous_[idx]; + return dp; + } + + inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) + { + // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription + Scalar beta = aquifer_influx_constant(); + Scalar Tc = time_constant(); + Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; + Scalar td = timer.simulationTimeElapsed() / Tc; + Scalar PItdprime = coeff_[1] + 2.0*coeff_[2]*(td_plus_dt); + Scalar PItd = coeff_[0] + coeff_[1]*td_plus_dt + coeff_[2]*td_plus_dt*td_plus_dt; + a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_ * PItdprime) ) / ( PItd - td*PItdprime ); + b = beta / Tc / ( PItd - td*PItdprime); + } + + inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer) + { + Scalar a, b; + calculate_a_b_constants(a,b,idx,timer); + // This function implements Eq 5.7 of the EclipseTechnicalDescription + Qai_[idx] = area_fraction(idx)*( a - b * ( pressure_current_[idx] - pressure_previous_[idx] ) ); + } + + + }; // class AquiferCarterTracy + + +} // namespace Opm + +#endif \ No newline at end of file diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp new file mode 100644 index 000000000..a1ed97e8b --- /dev/null +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -0,0 +1,184 @@ +/* +<<<<<<< HEAD + File adapted from BlackoilWellModel.hpp + + Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface + Copyright 2017 Statoil ASA. +======= + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2017 IRIS AS +>>>>>>> 9ccee28... First addition of the class BlackoilAquiferModel. + + 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_BLACKOILAQUIFERMODEL_HEADER_INCLUDED +#define OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED + +#include + +#include +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + + +#include + + +namespace Opm { + + /// Class for handling the blackoil well model. + template + class BlackoilAquiferModel { + + public: + + + // --------- Types --------- + typedef BlackoilModelParameters ModelParameters; + + + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + + static const int numEq = BlackoilIndices::numEq; + static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; + + typedef Ewoms::BlackOilPolymerModule PolymerModule; + + typedef AquiferCarterTracy Aquifer_object; + + BlackoilAquiferModel(Simulator& ebosSimulator, + const ModelParameters& param, + const bool terminal_output); + + // compute the well fluxes and assemble them in to the reservoir equations as source terms + // and in the well equations. + void assemble( const SimulatorTimerInterface& timer, + const int iterationIdx ); + + // called at the beginning of a time step + void beginTimeStep(); + // called at the end of a time step + void timeStepSucceeded(); + + // called at the beginning of a report step + void beginReportStep(const int time_step); + + // called at the end of a report step + void endReportStep(); + + const SimulatorReport& lastReport() const; + + inline const Simulator& simulator() const + { + return ebosSimulator_; + } + + /// Hack function to get what I need from parser + void init(const Simulator& ebosSimulator, std::vector& aquifers); + + protected: + + Simulator& ebosSimulator_; + + const ModelParameters param_; + bool terminal_output_; + bool has_solvent_; + bool has_polymer_; + std::vector pvt_region_idx_; + PhaseUsage phase_usage_; + std::vector active_; + size_t global_nc_; + // the number of the cells in the local grid + size_t number_of_cells_; + double gravity_; + std::vector depth_; + std::vector aquifers_; + + + SimulatorReport last_report_; + + const Schedule& schedule() const + { return ebosSimulator_.gridManager().schedule(); } + + void updatePrimaryVariables(); + + void initPrimaryVariablesEvaluation() const; + + void updateConnectionIntensiveQuantities() const; + + void calculateExplicitQuantities(); + + // The number of components in the model. + int numComponents() const; + + int numAquifers() const; + + int numPhases() const; + + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; + + void assembleAquiferEq(const SimulatorTimerInterface& timer); + + SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer); + + // some preparation work, mostly related to group control and RESV, + // at the beginning of each time step (Not report step) + void prepareTimeStep(const SimulatorTimerInterface& timer); + + const std::vector& aquifers(); + + }; + + +} // namespace Opm + +#include "BlackoilAquiferModel_impl.hpp" +#endif diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp new file mode 100644 index 000000000..f9b2c4806 --- /dev/null +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -0,0 +1,287 @@ +namespace Opm { + + + template + BlackoilAquiferModel:: + BlackoilAquiferModel(Simulator& ebosSimulator, + const ModelParameters& param, + const bool terminal_output) + : ebosSimulator_(ebosSimulator) + , param_(param) + , terminal_output_(terminal_output) + , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) + , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) + { + const auto& eclState = ebosSimulator_.gridManager().eclState(); + phase_usage_ = phaseUsageFromDeck(eclState); + + active_.resize(phase_usage_.MaxNumPhases, false); + for (int p = 0; p < phase_usage_.MaxNumPhases; ++p) { + active_[ p ] = phase_usage_.phase_used[ p ] != 0; + } + + const auto& gridView = ebosSimulator_.gridView(); + + // calculate the number of elements of the compressed sequential grid. this needs + // to be done in two steps because the dune communicator expects a reference as + // argument for sum() + number_of_cells_ = gridView.size(/*codim=*/0); + global_nc_ = gridView.comm().sum(number_of_cells_); + gravity_ = ebosSimulator_.problem().gravity()[2]; + init(ebosSimulator_, aquifers_); + } + + + + // called at the beginning of a time step + template + void + BlackoilAquiferModel:: beginTimeStep() + { + // Right now it doesn't do shit. + } + + // called at the end of a time step + template + void + BlackoilAquiferModel:: timeStepSucceeded() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + { + aquifer->after_time_step(); + } + } + + // called at the beginning of a report step + template + void + BlackoilAquiferModel:: beginReportStep(const int time_step) + { + // Right now it doesn't do shit. + } + + // called at the end of a report step + template + void + BlackoilAquiferModel:: endReportStep() + { + // Right now it just spits out the constants for each aquifers + // We are using the simple integer indexing for the aquifers + for (int i = 0; i < numAquifers(); ++i) + { + std::cout << "Aquifer[" << i << "]" + << " : Tc = " << aquifers()[i].time_constant() + << ", beta = " << aquifers()[i].aquifer_influx_constant() << std::endl; + } + } + + // Get the last report step + template + const SimulatorReport& + BlackoilAquiferModel:: lastReport() const + { + for (auto i = aquifers_.begin(); i != aquifers_.end(); ++i){ + (*i).print_private_members(); + } + return last_report_; + } + + template + void + BlackoilAquiferModel:: + assemble( const SimulatorTimerInterface& timer, + const int iterationIdx ) + { + last_report_ = SimulatorReport(); + + // We need to update the reservoir pressures connected to the aquifer + updateConnectionIntensiveQuantities(); + + if (iterationIdx == 0) { + // We can do the Table check and coefficients update in this function + // For now, it does nothing! + prepareTimeStep(timer); + } + + if (iterationIdx == 0) { + calculateExplicitQuantities(); + } + + if (param_.solve_aquifereq_initially_ && iterationIdx == 0) { + // solve the aquifer equations as a pre-processing step + last_report_ = solveAquiferEq(timer); + } + + assembleAquiferEq(timer); + + last_report_.converged = true; + } + + // Protected function: Update the primary variables + template + void + BlackoilAquiferModel:: updatePrimaryVariables() + { + // Right now it doesn't do shit. + } + + // Protected function: Init the primary variables + template + void + BlackoilAquiferModel:: initPrimaryVariablesEvaluation() const + { + // Right now it doesn't do shit. + } + + template + void + BlackoilAquiferModel:: updateConnectionIntensiveQuantities() const + { + ElementContext elemCtx(ebosSimulator_); + const auto& gridView = ebosSimulator_.gridView(); + const auto& elemEndIt = gridView.template end(); + for (auto elemIt = gridView.template begin(); + elemIt != elemEndIt; + ++elemIt) + { + elemCtx.updatePrimaryStencil(*elemIt); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + } + } + + template + void + BlackoilAquiferModel:: calculateExplicitQuantities() + { + // for (auto aqui = aquifers_.begin(); aqui!= aquifers_.end(); ++aqui) + // { + // std::cout << "calculateExplicitQuantities: Aquifer id = " << aqui->aquiferID() << std::endl; + // aqui->calculateExplicitQuantities(ebosSimulator_); + // } + } + + + template + SimulatorReport + BlackoilAquiferModel:: solveAquiferEq(const SimulatorTimerInterface& timer) + { + // We need to solve the equilibrium equation first to + // obtain the initial pressure of water in the aquifer + SimulatorReport report = SimulatorReport(); + return report; + } + + // Protected function: Return number of components in the model. + template + int + BlackoilAquiferModel:: numComponents() const + { + if (numPhases() == 2) { + return 2; + } + int numComp = FluidSystem::numComponents; + if (has_solvent_) { + numComp ++; + } + + return numComp; + } + + // Protected function: Return number of aquifers in the model. + template + int + BlackoilAquiferModel:: numAquifers() const + { + return aquifers_.size(); + } + + // Protected function: Return number of phases in the model. + template + int + BlackoilAquiferModel:: numPhases() const + { + // Not implemented yet!!!!!!!!!!!! + const auto& pu = phase_usage_; + return pu.num_phases; + } + + + // Protected function: returns the phase index in ebos + template + int + BlackoilAquiferModel:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const auto& pu = phase_usage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return FluidSystem::waterPhaseIdx; + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return FluidSystem::oilPhaseIdx; + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return FluidSystem::gasPhaseIdx; + + assert(phaseIdx < 3); + // for other phases return the index + return phaseIdx; + } + + // Protected function which calls the individual aquifer models + template + void + BlackoilAquiferModel:: assembleAquiferEq(const SimulatorTimerInterface& timer) + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + { + std::cout << "assembleAquiferEq: Aquifer id = " << aquifer->aquiferID() << std::endl; + aquifer->assembleAquiferEq(ebosSimulator_, timer); + } + } + + // Protected function + // some preparation work, mostly related to group control and RESV, + // at the beginning of each time step (Not report step) + template + void BlackoilAquiferModel:: prepareTimeStep(const SimulatorTimerInterface& timer) + { + // Here we can ask each carter tracy aquifers to get the current previous time step's pressure + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + { + aquifer->before_time_step(ebosSimulator_, timer); + } + } + + // Protected function: Returns a reference to the aquifers members in the model + template + const std::vector< AquiferCarterTracy >& + BlackoilAquiferModel:: aquifers() + { + return aquifers_; + } + + + // Initialize the aquifers in the deck + template + void + BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers)//, std::vector< AquiferCarterTracy >& aquifers) + { + const auto& deck = ebosSimulator.gridManager().deck(); + const auto& eclState = ebosSimulator.gridManager().eclState(); + + // Get all the carter tracy aquifer properties data and put it in aquifers vector + AquiferCT aquiferct = AquiferCT(eclState,deck); + + std::vector aquifersData = aquiferct.getAquifers(); + std::vector aquanconData = aquiferct.getAquancon(); + + // for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) + // { + + // } + + auto ita = aquifersData.cbegin(); + auto f_lambda = [&] (AquiferCT::AQUANCON_data i) { + aquifers.push_back( AquiferCarterTracy (*ita++, i, numComponents(), gravity_ ) ); + }; + std::for_each( aquanconData.cbegin(), aquanconData.cend(), f_lambda ); + } + +} // namespace Opm \ No newline at end of file diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index c16461496..c55f337bd 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -30,6 +30,9 @@ #include #include #include +#include +#include +#include #include #include @@ -143,6 +146,7 @@ namespace Opm { BlackoilModelEbos(Simulator& ebosSimulator, const ModelParameters& param, BlackoilWellModel& well_model, + BlackoilAquiferModel& aquifer_model, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output ) @@ -157,6 +161,7 @@ namespace Opm { , has_energy_(GET_PROP_VALUE(TypeTag, EnableEnergy)) , param_( param ) , well_model_ (well_model) + , aquifer_model_(aquifer_model) , terminal_output_ (terminal_output) , current_relaxation_(1.0) , dx_old_(UgGridHelpers::numCells(grid_)) @@ -204,6 +209,7 @@ namespace Opm { wasSwitched_.resize(numDof); std::fill(wasSwitched_.begin(), wasSwitched_.end(), false); + aquiferModel().beginTimeStep(); wellModel().beginTimeStep(); if (param_.update_equations_scaling_) { @@ -349,6 +355,7 @@ namespace Opm { DUNE_UNUSED_PARAMETER(well_state); wellModel().timeStepSucceeded(); + aquiferModel().timeStepSucceeded(); ebosSimulator_.problem().endTimeStep(); } @@ -366,9 +373,22 @@ namespace Opm { ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - // -------- Well equations ---------- + // -------- Well and aquifer common variables ---------- double dt = timer.currentStepLength(); + // -------- Aquifer models ---------- + try + { + // Modify the Jacobian and residuals according to the aquifer models + aquiferModel().assemble(timer, iterationIdx); + } + catch( const Dune::FMatrixError& e ) + { + OPM_THROW(Opm::NumericalProblem,"Error when assembling aquifer models"); + } + + // -------- Well equations ---------- + try { // assembles the well equations and applies the wells to @@ -1081,6 +1101,9 @@ namespace Opm { // Well Model BlackoilWellModel& well_model_; + // Aquifer Model + BlackoilAquiferModel& aquifer_model_; + /// \brief Whether we print something to std::cout bool terminal_output_; /// \brief The number of cells of the global grid. @@ -1100,6 +1123,41 @@ namespace Opm { const BlackoilWellModel& wellModel() const { return well_model_; } + BlackoilAquiferModel& + aquiferModel() { return aquifer_model_; } + + const BlackoilAquiferModel& + aquiferModel() const { return aquifer_model_; } + + int flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + + // for other phases return the index + return phaseIdx; + } + + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return FluidSystem::waterPhaseIdx; + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return FluidSystem::oilPhaseIdx; + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return FluidSystem::gasPhaseIdx; + + assert(phaseIdx < 3); + // for other phases return the index + return phaseIdx; + } + void beginReportStep() { ebosSimulator_.problem().beginEpisode(); diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index f16fc7685..2bd302c2d 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -63,6 +63,7 @@ namespace Opm param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day ); max_strict_iter_ = param.getDefault("max_strict_iter",8); solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_); + solve_aquifereq_initially_ = param.getDefault("solve_aquifereq_initially",solve_aquifereq_initially_); update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_); use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_); deck_file_name_ = param.template get("deck_filename"); @@ -94,6 +95,7 @@ namespace Opm max_inner_iter_ms_wells_ = 10; maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day ); solve_welleq_initially_ = true; + solve_aquifereq_initially_ = true; update_equations_scaling_ = false; use_update_stabilization_ = true; use_multisegment_well_ = false; diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index b01cfdde9..d57552161 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -77,6 +77,9 @@ namespace Opm /// Solve well equation initially bool solve_welleq_initially_; + /// Solve aquifer equation initially + bool solve_aquifereq_initially_; + /// Update scaling factors for mass balance equations bool update_equations_scaling_; diff --git a/opm/autodiff/BlackoilTransportModel.hpp b/opm/autodiff/BlackoilTransportModel.hpp index ffc0ebc92..d6610f0d5 100644 --- a/opm/autodiff/BlackoilTransportModel.hpp +++ b/opm/autodiff/BlackoilTransportModel.hpp @@ -86,6 +86,7 @@ namespace Opm { { Base::prepareStep(timer, reservoir_state, well_state); Base::param_.solve_welleq_initially_ = false; + Base::param_.solve_aquifereq_initially_ = false; SolutionState state0 = variableState(reservoir_state, well_state); asImpl().makeConstantState(state0); asImpl().computeAccum(state0, 0); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 4445eb390..5e5988b22 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -29,6 +29,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -65,6 +68,7 @@ public: typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; typedef BlackoilWellModel WellModel; + typedef BlackoilAquiferModel AquiferModel; /// Initialise from parameters and objects to observe. @@ -186,6 +190,8 @@ public: auto auxMod = std::make_shared >(schedule(), grid()); ebosSimulator_.model().addAuxiliaryModule(auxMod); } + AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_); + // aquifer_model.hack_init(ebosSimulator_); // Main simulation loop. while (!timer.done()) { @@ -202,7 +208,9 @@ public: well_model.beginReportStep(timer.currentStepNum()); - auto solver = createSolver(well_model); + aquifer_model.beginReportStep(timer.currentStepNum()); + + auto solver = createSolver(well_model, aquifer_model); // write the inital state at the report stage if (timer.initialStep()) { @@ -266,6 +274,7 @@ public: } solver->model().endReportStep(); + aquifer_model.endReportStep(); well_model.endReportStep(); // take time that was used to solve system for this reportStep @@ -308,6 +317,9 @@ public: total_timer.stop(); report.total_time = total_timer.secsSinceStart(); report.converged = true; + + auto reportaquifer = aquifer_model.lastReport(); + return report; } @@ -320,11 +332,12 @@ public: protected: - std::unique_ptr createSolver(WellModel& well_model) + std::unique_ptr createSolver(WellModel& well_model, AquiferModel& aquifer_model) { auto model = std::unique_ptr(new Model(ebosSimulator_, model_param_, well_model, + aquifer_model, solver_, terminal_output_)); From ae22d18e20529a68d587e6173b4a1b2cd0c74109 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Mon, 11 Dec 2017 17:38:45 +0100 Subject: [PATCH 06/31] Working and building skeleton which works for multiple aquifers. Also gets the properties of the aquifers using AQUCT keyword. --- opm/autodiff/AQUCT_params.hpp | 31 ++++++++++ opm/autodiff/BlackoilAquiferModel.hpp | 7 +-- opm/autodiff/BlackoilAquiferModel_impl.hpp | 66 ++++++++++++++++++++++ 3 files changed, 99 insertions(+), 5 deletions(-) create mode 100644 opm/autodiff/AQUCT_params.hpp diff --git a/opm/autodiff/AQUCT_params.hpp b/opm/autodiff/AQUCT_params.hpp new file mode 100644 index 000000000..da6c03ff5 --- /dev/null +++ b/opm/autodiff/AQUCT_params.hpp @@ -0,0 +1,31 @@ +#ifndef OPM_AQUCT_HEADER_INCLUDED +#define OPM_AQUCT_HEADER_INCLUDED + + + struct AQUCT_params{ + + // Aquifer ID + int aquiferID; + // Table IDs + int inftableID, pvttableID; + // Perforation cell id + int cell_id; + // Variables constants + double phi_aq , //aquifer porosity + d0, //aquifer datum depth + C_t , //total compressibility + r_o , //aquifer inner radius + k_a , //aquifer permeability + c1, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) + h , //aquifer thickness + theta , //angle subtended by the aquifer boundary + c2 ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). + std::vector td, pi; + + + }; + + + + +#endif \ No newline at end of file diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index a1ed97e8b..03771dced 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -1,15 +1,13 @@ /* -<<<<<<< HEAD File adapted from BlackoilWellModel.hpp Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface Copyright 2017 Statoil ASA. -======= Copyright 2016 SINTEF ICT, Applied Mathematics. Copyright 2016 - 2017 Statoil ASA. Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services Copyright 2016 - 2017 IRIS AS ->>>>>>> 9ccee28... First addition of the class BlackoilAquiferModel. + This file is part of the Open Porous Media project (OPM). @@ -62,7 +60,6 @@ #include #include - #include @@ -74,7 +71,6 @@ namespace Opm { public: - // --------- Types --------- typedef BlackoilModelParameters ModelParameters; @@ -156,6 +152,7 @@ namespace Opm { void calculateExplicitQuantities(); + // The number of components in the model. int numComponents() const; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index f9b2c4806..b9a861e44 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -284,4 +284,70 @@ namespace Opm { std::for_each( aquanconData.cbegin(), aquanconData.cend(), f_lambda ); } + // Begin the hack to initialize the aquifers in the deck + template + std::vector< AquiferCarterTracy > + BlackoilAquiferModel:: hack_init(const Simulator& ebosSimulator)//, std::vector< AquiferCarterTracy >& aquifers) + { + std::vector< AquiferCarterTracy > aquifers; + /** Begin hack!!!!! */ + const auto& deck = ebosSimulator.gridManager().deck(); + const auto& eclState = ebosSimulator.gridManager().eclState(); + + if (!deck.hasKeyword("AQUCT")){ + std::cout << "Nothing is shown! Where is AQUCT!????" << std::endl; + } + + const auto& aquctKeyword = deck.getKeyword("AQUCT"); + std::vector aquctParams; + // Resize the parameter vector container based on row entries in aquct + // We do the same for aquifers too because number of aquifers is assumed to be for each entry in aquct + aquctParams.resize(aquctKeyword.size()); + // aquifers.resize(aquctKeyword.size()); + + const int tableID = 0; + + std::cout << "Parsing AQUCT stuff" << std::endl; + for (size_t aquctRecordIdx = 0; aquctRecordIdx < aquctKeyword.size(); ++ aquctRecordIdx) + { + const auto& aquctRecord = aquctKeyword.getRecord(aquctRecordIdx); + + aquctParams.at(aquctRecordIdx).aquiferID = aquctRecord.getItem("AQUIFER_ID").template get(0); + aquctParams.at(aquctRecordIdx).h = aquctRecord.getItem("THICKNESS_AQ").template get(0); + aquctParams.at(aquctRecordIdx).phi_aq = aquctRecord.getItem("PORO_AQ").template get(0); + aquctParams.at(aquctRecordIdx).d0 = aquctRecord.getItem("DAT_DEPTH").getSIDouble(0); + aquctParams.at(aquctRecordIdx).C_t = aquctRecord.getItem("C_T").template get(0); + aquctParams.at(aquctRecordIdx).r_o = aquctRecord.getItem("RAD").getSIDouble(0); + aquctParams.at(aquctRecordIdx).k_a = aquctRecord.getItem("PERM_AQ").getSIDouble(0); + aquctParams.at(aquctRecordIdx).theta = aquctRecord.getItem("INFLUENCE_ANGLE").template get(0); + aquctParams.at(aquctRecordIdx).c1 = 0.008527; // We are using SI + aquctParams.at(aquctRecordIdx).c2 = 6.283; + aquctParams.at(aquctRecordIdx).inftableID = aquctRecord.getItem("TABLE_NUM_INFLUENCE_FN").template get(0) - 1; + aquctParams.at(aquctRecordIdx).pvttableID = aquctRecord.getItem("TABLE_NUM_WATER_PRESS").template get(0) - 1; + + std::cout << aquctParams.at(aquctRecordIdx).inftableID << std::endl; + // Get the correct influence table values + const auto& aqutabTable = eclState.getTableManager().getAqutabTables().getTable(aquctParams.at(aquctRecordIdx).inftableID); + const auto& aqutab_tdColumn = aqutabTable.getColumn(0); + const auto& aqutab_piColumn = aqutabTable.getColumn(1); + aquctParams.at(aquctRecordIdx).td = aqutab_tdColumn.vectorCopy(); + aquctParams.at(aquctRecordIdx).pi = aqutab_piColumn.vectorCopy(); + + // We determine the cell perforation here. + int cellID = 10 + aquctRecordIdx; + + aquctParams.at(aquctRecordIdx).cell_id = cellID; + + // We do not have mu_w as it has to be calculated from pvttable + aquifers.push_back(Aquifer_object( aquctParams.at(aquctRecordIdx) )); + } + + // I want to deliberately add another aquifer + aquifers.push_back( Aquifer_object(99) ); + + // aquifers_ = aquifers; + + return aquifers; + } + } // namespace Opm \ No newline at end of file From 5ca8bb676dcc59b615adbd1f8d37ce5699b69654 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Fri, 15 Dec 2017 13:24:15 +0100 Subject: [PATCH 07/31] Working aquifer object with equations implemented. Still need to do initialization solve. --- opm/autodiff/BlackoilAquiferModel.hpp | 1 - opm/autodiff/BlackoilAquiferModel_impl.hpp | 55 ++-------------------- 2 files changed, 5 insertions(+), 51 deletions(-) diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index 03771dced..f8b6e5c9f 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -152,7 +152,6 @@ namespace Opm { void calculateExplicitQuantities(); - // The number of components in the model. int numComponents() const; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index b9a861e44..34d02203e 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -294,60 +294,15 @@ namespace Opm { const auto& deck = ebosSimulator.gridManager().deck(); const auto& eclState = ebosSimulator.gridManager().eclState(); - if (!deck.hasKeyword("AQUCT")){ - std::cout << "Nothing is shown! Where is AQUCT!????" << std::endl; - } + // Get all the carter tracy aquifer properties data and put it in aquifers vector + AquiferCT aquiferct = AquiferCT(eclState,deck); - const auto& aquctKeyword = deck.getKeyword("AQUCT"); - std::vector aquctParams; - // Resize the parameter vector container based on row entries in aquct - // We do the same for aquifers too because number of aquifers is assumed to be for each entry in aquct - aquctParams.resize(aquctKeyword.size()); - // aquifers.resize(aquctKeyword.size()); + std::vector aquifersData = aquiferct.getAquifers(); - const int tableID = 0; - - std::cout << "Parsing AQUCT stuff" << std::endl; - for (size_t aquctRecordIdx = 0; aquctRecordIdx < aquctKeyword.size(); ++ aquctRecordIdx) + for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) { - const auto& aquctRecord = aquctKeyword.getRecord(aquctRecordIdx); - - aquctParams.at(aquctRecordIdx).aquiferID = aquctRecord.getItem("AQUIFER_ID").template get(0); - aquctParams.at(aquctRecordIdx).h = aquctRecord.getItem("THICKNESS_AQ").template get(0); - aquctParams.at(aquctRecordIdx).phi_aq = aquctRecord.getItem("PORO_AQ").template get(0); - aquctParams.at(aquctRecordIdx).d0 = aquctRecord.getItem("DAT_DEPTH").getSIDouble(0); - aquctParams.at(aquctRecordIdx).C_t = aquctRecord.getItem("C_T").template get(0); - aquctParams.at(aquctRecordIdx).r_o = aquctRecord.getItem("RAD").getSIDouble(0); - aquctParams.at(aquctRecordIdx).k_a = aquctRecord.getItem("PERM_AQ").getSIDouble(0); - aquctParams.at(aquctRecordIdx).theta = aquctRecord.getItem("INFLUENCE_ANGLE").template get(0); - aquctParams.at(aquctRecordIdx).c1 = 0.008527; // We are using SI - aquctParams.at(aquctRecordIdx).c2 = 6.283; - aquctParams.at(aquctRecordIdx).inftableID = aquctRecord.getItem("TABLE_NUM_INFLUENCE_FN").template get(0) - 1; - aquctParams.at(aquctRecordIdx).pvttableID = aquctRecord.getItem("TABLE_NUM_WATER_PRESS").template get(0) - 1; - - std::cout << aquctParams.at(aquctRecordIdx).inftableID << std::endl; - // Get the correct influence table values - const auto& aqutabTable = eclState.getTableManager().getAqutabTables().getTable(aquctParams.at(aquctRecordIdx).inftableID); - const auto& aqutab_tdColumn = aqutabTable.getColumn(0); - const auto& aqutab_piColumn = aqutabTable.getColumn(1); - aquctParams.at(aquctRecordIdx).td = aqutab_tdColumn.vectorCopy(); - aquctParams.at(aquctRecordIdx).pi = aqutab_piColumn.vectorCopy(); - - // We determine the cell perforation here. - int cellID = 10 + aquctRecordIdx; - - aquctParams.at(aquctRecordIdx).cell_id = cellID; - - // We do not have mu_w as it has to be calculated from pvttable - aquifers.push_back(Aquifer_object( aquctParams.at(aquctRecordIdx) )); + aquifers.push_back( AquiferCarterTracy (*aquiferData, numComponents(), gravity_ ) ); } - - // I want to deliberately add another aquifer - aquifers.push_back( Aquifer_object(99) ); - - // aquifers_ = aquifers; - - return aquifers; } } // namespace Opm \ No newline at end of file From 7e1bb457cd35cba8cda92703359cf36a1add8c17 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 4 Jan 2018 15:50:04 +0100 Subject: [PATCH 08/31] Implemented initialization procedure. --- opm/autodiff/AQUCT_params.hpp | 31 -- opm/autodiff/AquiferCarterTracy.hpp | 349 ++++++++++-------- opm/autodiff/BlackoilAquiferModel.hpp | 25 +- opm/autodiff/BlackoilAquiferModel_impl.hpp | 118 +----- opm/autodiff/BlackoilModelEbos.hpp | 54 +-- .../SimulatorFullyImplicitBlackoilEbos.hpp | 2 +- 6 files changed, 243 insertions(+), 336 deletions(-) delete mode 100644 opm/autodiff/AQUCT_params.hpp diff --git a/opm/autodiff/AQUCT_params.hpp b/opm/autodiff/AQUCT_params.hpp deleted file mode 100644 index da6c03ff5..000000000 --- a/opm/autodiff/AQUCT_params.hpp +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef OPM_AQUCT_HEADER_INCLUDED -#define OPM_AQUCT_HEADER_INCLUDED - - - struct AQUCT_params{ - - // Aquifer ID - int aquiferID; - // Table IDs - int inftableID, pvttableID; - // Perforation cell id - int cell_id; - // Variables constants - double phi_aq , //aquifer porosity - d0, //aquifer datum depth - C_t , //total compressibility - r_o , //aquifer inner radius - k_a , //aquifer permeability - c1, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) - h , //aquifer thickness - theta , //angle subtended by the aquifer boundary - c2 ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). - std::vector td, pi; - - - }; - - - - -#endif \ No newline at end of file diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index b6f9cfa17..4ba7597c4 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -30,10 +31,12 @@ #include #include +#include #include #include #include +#include #include #include @@ -65,44 +68,30 @@ namespace Opm typedef double Scalar; typedef DenseAd::Evaluation Eval; + typedef Opm::BlackOilFluidState FluidState; + typedef typename FluidSystem::WaterPvt WaterPvt; typedef Ewoms::BlackOilPolymerModule PolymerModule; + + + static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; - AquiferCarterTracy(const std::vector& cell_id) - : phi_aq_ (1.0), // - C_t_ (1.0), // - r_o_ (1.0), // - k_a_ (1.0), // - c1_ (1.0), - h_ (1.0), // - theta_ (1.0), // - c2_ (1.0), // - d0_ (1.0), - cell_idx_ (cell_id) - { - mu_w_ = 1e-3; - aqutab_td_.push_back(1.0); - aqutab_pi_.push_back(1.0); - aquiferID_ = 1; - inftableID_ = 1; - pvttableID_ = 1; - init_quantities(); - } - - explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const AquiferCT::AQUANCON_data& aquanconParams, - const int numComponents, const Scalar gravity ) + explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, + const int numComponents, const Scalar gravity, const Simulator& ebosSimulator ) : phi_aq_ (params.phi_aq), // C_t_ (params.C_t), // r_o_ (params.r_o), // k_a_ (params.k_a), // c1_ (params.c1), h_ (params.h), // + p0_defaulted_ (params.p0_defaulted), + pa0_ (params.p0), theta_ (params.theta), // c2_ (params.c2), // d0_ (params.d0), @@ -111,12 +100,11 @@ namespace Opm aquiferID_ (params.aquiferID), inftableID_ (params.inftableID), pvttableID_ (params.pvttableID), - cell_idx_ (params.cell_id), num_components_ (numComponents), - gravity_ (gravity) + gravity_ (gravity), + ebos_simulator_ (ebosSimulator) { - mu_w_ = 1e-3; - init_quantities(aquanconParams); + init_quantities(connection); } inline const PhaseUsage& @@ -127,79 +115,35 @@ namespace Opm return *phase_usage_; } - inline int - flowPhaseToEbosCompIdx( const int phaseIdx ) const - { - const auto& pu = phaseUsage(); - if (pu.phase_pos[Water] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - if (pu.phase_pos[Oil] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (pu.phase_pos[Gas] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - - // for other phases return the index - return phaseIdx; - } - - inline int - flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const auto& pu = phaseUsage(); - if (pu.phase_pos[Water] == phaseIdx) { - return FluidSystem::waterPhaseIdx; - } - if (pu.phase_pos[Oil] == phaseIdx) { - return FluidSystem::oilPhaseIdx; - } - if (pu.phase_pos[Gas] == phaseIdx) { - return FluidSystem::gasPhaseIdx; - } - - assert(phaseIdx < 3); - // for other phases return the index - return phaseIdx; - } - - inline void calculateExplicitQuantities(const Simulator& ebosSimulator) - { - std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; - } inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) { - std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; - // resAqui_ = 0.0; dt_ = timer.currentStepLength(); auto& ebosJac = ebosSimulator.model().linearizer().matrix(); auto& ebosResid = ebosSimulator.model().linearizer().residual(); - // TODO: it probably can be static member for StandardWell - const double volume = 0.002831684659200; // 0.1 cu ft; auto cellID = cell_idx_.begin(); + size_t idx; for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) { Eval qinflow = 0.0; // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to // IntensiveQuantities of that particular cell_id - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + const IntensiveQuantities intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); // This is the pressure at td + dt get_current_Pressure_cell(pressure_current_,idx,intQuants); get_current_density_cell(rhow_,idx,intQuants); calculate_inflow_rate(idx, timer); - qinflow = Qai_[idx]; - ebosResid[*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)] -= qinflow.value(); + qinflow = Qai_.at(idx); + ebosResid[*cellID][(FluidSystem::waterCompIdx)] -= qinflow.value(); for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[*cellID][*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)][pvIdx] -= qinflow.derivative(pvIdx); + ebosJac[*cellID][*cellID][(FluidSystem::waterCompIdx)][pvIdx] -= qinflow.derivative(pvIdx); } - std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ - // << " -> P_wat[t+dt] = " << pressure_current_[idx] << std::endl - << " Qai[t+dt] = " << Qai_[idx] << std::endl; } } @@ -214,50 +158,19 @@ namespace Opm } } - inline void after_time_step() + inline void after_time_step(const SimulatorTimerInterface& timer) { for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) { - W_flux_ += (*Qai); - } - std::cout << "Aquifer # " << aquiferID_ << ": My cumulative flux = " << W_flux_ << std::endl; - } - - /* Made into public for testing only!!!!!!. Must be protected */ - inline const Scalar time_constant() const - { - Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); - return Tc; - } - - /* Made into public for testing only!!!!!!. Must be protected */ - inline const Scalar aquifer_influx_constant() const - { - Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; - return beta; - } - - // This is another hack to get the face area only for SPE1. - // Ideally it should be a map which given a cell_id, it returns the area fraction - inline const double area_fraction(const int i) - { - return 1000.0*20.0*0.092903/(1000.0*1000.0*0.092903*2 + 1000.0*20.0*0.092903*4); - } - - inline void print_private_members() const - { - std::cout << "Aquifer CT #" << aquiferID_ << std::endl; - auto ita = aqutab_td_.cbegin(); - auto f_lambda = [&ita] (double i) {std::cout << *ita++ << " " << i << std::endl;}; - std::for_each( aqutab_pi_.cbegin(), aqutab_pi_.cend(), f_lambda ); - - for (auto i = coeff_.begin(); i != coeff_.end(); ++i ) - { - std::cout << "Coeff = " << *i << std::endl; + W_flux_ += (*Qai)*timer.currentStepLength(); } } - /* Made into public for testing only!!!!!!. Must be protected */ + inline const double area_fraction(const size_t i) + { + return alphai_.at(i); + } + inline const std::vector cell_id() const { return cell_idx_; @@ -269,9 +182,9 @@ namespace Opm } - - protected: + private: const PhaseUsage* phase_usage_; + const Simulator& ebos_simulator_; // Aquifer ID, and other IDs @@ -279,14 +192,17 @@ namespace Opm int num_components_; // Grid variables - std::vector cell_idx_; + + std::vector cell_idx_; + std::vector faceArea_connected_; // Quantities at each grid id std::vector cell_depth_; - std::vector pressure_previous_; - std::vector pressure_current_; - std::vector Qai_; - std::vector rhow_; + std::vector pressure_previous_; + std::vector pressure_current_; + std::vector Qai_; + std::vector rhow_; + std::vector alphai_; // Variables constants Scalar mu_w_ , //water viscosity @@ -304,8 +220,9 @@ namespace Opm std::vector aqutab_td_, aqutab_pi_; // Cumulative flux - Scalar W_flux_, dt_, pa0_, gravity_; - + Scalar dt_, pa0_, gravity_; + bool p0_defaulted_; + Eval W_flux_; // Also return the polynomial fit std::vector coeff_; @@ -316,7 +233,7 @@ namespace Opm inline void polynomial_fit( const std::vector &X, const std::vector &y, std::vector &coeff, int order, bool bias) const { - int colNum = (bias)? order + 1 : order; + size_t colNum = (bias)? order + 1 : order; Eigen::MatrixXd A(X.size(), colNum); Eigen::VectorXd y_mapped = Eigen::VectorXd::Map(&y.front(), y.size()); Eigen::VectorXd result; @@ -337,62 +254,202 @@ namespace Opm coeff[i] = result[i]; } - inline void init_quantities(const AquiferCT::AQUANCON_data& aquanconParams) + inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { + // http://kluge.in-chemnitz.de/opensource/spline/ + } + + inline void init_quantities(const Aquancon::AquanconOutput& connection) + { + // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 W_flux_ = 0.; - // pa0_ is the initial aquifer water pressure. Must be calculated from equilibrium if left default, - // or we get the information from the deck Hacked to make it at 45e6 Pa - pa0_ = 45e6; + + // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function + initialize_connections(connection); + + calculate_aquifer_condition(); pressure_previous_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.); - // We hack the cell depth values for now. We can actually get it from elementcontext pos - cell_depth_.resize(cell_idx_.size(), d0_); - rhow_.resize(cell_idx_.size(), 998.0); - Qai_.resize(cell_idx_.size(), 0.); + Qai_.resize(cell_idx_.size(), 0.0); - polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 2, true); + polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 1, true); } - inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - pressure_water[idx] = fs.pressure(FluidSystem::waterPhaseIdx).value(); + pressure_water.at(idx) = fs.pressure(FluidSystem::waterPhaseIdx); } - inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) + inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - rho_water[idx] = fs.density(FluidSystem::waterPhaseIdx).value(); + rho_water.at(idx) = fs.density(FluidSystem::waterPhaseIdx); } inline Scalar dpai(int idx) { - Scalar dp = pa0_ - rhow_[idx]*gravity_*(cell_depth_[idx] - d0_) - pressure_previous_[idx]; + Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_) - pressure_previous_.at(idx).value(); return dp; } + // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) { - // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription Scalar beta = aquifer_influx_constant(); Scalar Tc = time_constant(); Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; Scalar td = timer.simulationTimeElapsed() / Tc; - Scalar PItdprime = coeff_[1] + 2.0*coeff_[2]*(td_plus_dt); - Scalar PItd = coeff_[0] + coeff_[1]*td_plus_dt + coeff_[2]*td_plus_dt*td_plus_dt; - a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_ * PItdprime) ) / ( PItd - td*PItdprime ); - b = beta / Tc / ( PItd - td*PItdprime); + + Scalar PItdprime = coeff_.at(1); + Scalar PItd = coeff_.at(0) + coeff_.at(1)*td_plus_dt; + a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); + b = beta / (Tc * ( PItd - td*PItdprime)); } + // This function implements Eq 5.7 of the EclipseTechnicalDescription inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer) { Scalar a, b; calculate_a_b_constants(a,b,idx,timer); - // This function implements Eq 5.7 of the EclipseTechnicalDescription - Qai_[idx] = area_fraction(idx)*( a - b * ( pressure_current_[idx] - pressure_previous_[idx] ) ); + Qai_.at(idx) = area_fraction(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx).value() ) ); } - + + inline const Scalar time_constant() const + { + Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); + return Tc; + } + + inline const Scalar aquifer_influx_constant() const + { + Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; + return beta; + } + + // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer + inline void initialize_connections(const Aquancon::AquanconOutput& connection) + { + const auto& eclState = ebos_simulator_.vanguard().eclState(); + const auto& ugrid = ebos_simulator_.vanguard().grid(); + const auto& grid = eclState.getInputGrid(); + + cell_idx_ = connection.global_index; + + assert( cell_idx_ == connection.global_index); + assert( (cell_idx_.size() == connection.influx_coeff.size()) ); + assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) ); + assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) ); + + // We hack the cell depth values for now. We can actually get it from elementcontext pos + cell_depth_.resize(cell_idx_.size(), d0_); + alphai_.resize(cell_idx_.size(), 1.0); + faceArea_connected_.resize(cell_idx_.size(),0.0); + + auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); + auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); + + // Translate the C face tag into the enum used by opm-parser's TransMult class + Opm::FaceDir::DirEnum faceDirection; + + // denom_face_areas is the sum of the areas connected to an aquifer + Scalar denom_face_areas = 0.; + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + { + auto cellFacesRange = cell2Faces[cell_idx_.at(idx)]; + + for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) + { + // The index of the face in the compressed grid + const int faceIdx = *cellFaceIter; + + // the logically-Cartesian direction of the face + const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); + + + if (faceTag == 0) // left + faceDirection = Opm::FaceDir::XMinus; + else if (faceTag == 1) // right + faceDirection = Opm::FaceDir::XPlus; + else if (faceTag == 2) // back + faceDirection = Opm::FaceDir::YMinus; + else if (faceTag == 3) // front + faceDirection = Opm::FaceDir::YPlus; + else if (faceTag == 4) // bottom + faceDirection = Opm::FaceDir::ZMinus; + else if (faceTag == 5) // top + faceDirection = Opm::FaceDir::ZPlus; + + if (faceDirection == connection.reservoir_face_dir.at(idx)) + { + faceArea_connected_.at(idx) = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); + denom_face_areas += faceArea_connected_.at(idx); + } + } + auto cellCenter = grid.getCellCenter(cell_idx_.at(idx)); + cell_depth_.at(idx) = cellCenter[2]; + } + + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + { + alphai_.at(idx) = faceArea_connected_.at(idx)/denom_face_areas; + } + } + + inline void calculate_aquifer_condition() + { + + int pvttableIdx = pvttableID_ - 1; + + rhow_.resize(cell_idx_.size(),0.); + + if (p0_defaulted_) + { + pa0_ = calculate_reservoir_equilibrium(rhow_); + } + + // Initialize a FluidState object first + FluidState fs_aquifer; + // We use the temperature of the first cell connected to the aquifer + // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs + fs_aquifer.assign( ebos_simulator_.model().cachedIntensiveQuantities(cell_idx_.at(0), /*timeIdx=*/ 0)->fluidState() ); + Eval temperature_aq, pa0_mean; + temperature_aq = fs_aquifer.temperature(0); + pa0_mean = pa0_; + + // rho_mean = FluidSystem::referenceDensity(waterPhaseIdx, pvttableIdx) + // *FluidSystem::waterPvt().inverseFormationVolumeFactor(pvttableIdx, temperature_aq, pa0_mean); + + Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); + std::cout << "Pa0 = " << pa0_mean << ", viscosity = " << mu_w_aquifer.value() << std::endl; + + mu_w_ = mu_w_aquifer.value(); + + } + + // This function is for calculating the aquifer properties from equilibrium state with the reservoir + inline Scalar calculate_reservoir_equilibrium(std::vector& rho_water_reservoir) + { + // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices + std::vector water_pressure_reservoir, pw_aquifer; + + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + { + size_t cellIDx = cell_idx_.at(idx); + const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0)); + const auto& fs = intQuants.fluidState(); + + water_pressure_reservoir.push_back( fs.pressure(FluidSystem::waterPhaseIdx).value() ); + rho_water_reservoir.at(idx) = fs.density(FluidSystem::waterPhaseIdx); + pw_aquifer.push_back( (water_pressure_reservoir.at(idx) - rho_water_reservoir.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_))*area_fraction(idx) ); + } + + // We take the average of the calculated equilibrium pressures. + Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.)/pw_aquifer.size(); + return aquifer_pres_avg; + } + }; // class AquiferCarterTracy diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index f8b6e5c9f..706756ecb 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -3,11 +3,6 @@ Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface Copyright 2017 Statoil ASA. - Copyright 2016 SINTEF ICT, Applied Mathematics. - Copyright 2016 - 2017 Statoil ASA. - Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2016 - 2017 IRIS AS - This file is part of the Open Porous Media project (OPM). @@ -40,6 +35,7 @@ #include #include +#include #include @@ -49,7 +45,6 @@ #include #include #include - #include #include @@ -86,8 +81,6 @@ namespace Opm { static const int numEq = BlackoilIndices::numEq; static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; - typedef Ewoms::BlackOilPolymerModule PolymerModule; - typedef AquiferCarterTracy Aquifer_object; BlackoilAquiferModel(Simulator& ebosSimulator, @@ -102,7 +95,7 @@ namespace Opm { // called at the beginning of a time step void beginTimeStep(); // called at the end of a time step - void timeStepSucceeded(); + void timeStepSucceeded(const SimulatorTimerInterface& timer); // called at the beginning of a report step void beginReportStep(const int time_step); @@ -117,7 +110,7 @@ namespace Opm { return ebosSimulator_; } - /// Hack function to get what I need from parser + // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy void init(const Simulator& ebosSimulator, std::vector& aquifers); protected: @@ -141,17 +134,9 @@ namespace Opm { SimulatorReport last_report_; - const Schedule& schedule() const - { return ebosSimulator_.gridManager().schedule(); } - - void updatePrimaryVariables(); - - void initPrimaryVariablesEvaluation() const; void updateConnectionIntensiveQuantities() const; - void calculateExplicitQuantities(); - // The number of components in the model. int numComponents() const; @@ -159,8 +144,6 @@ namespace Opm { int numPhases() const; - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; - void assembleAquiferEq(const SimulatorTimerInterface& timer); SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer); @@ -177,4 +160,4 @@ namespace Opm { } // namespace Opm #include "BlackoilAquiferModel_impl.hpp" -#endif +#endif \ No newline at end of file diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index 34d02203e..c75b3f53b 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -12,7 +12,7 @@ namespace Opm { , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) { - const auto& eclState = ebosSimulator_.gridManager().eclState(); + const auto& eclState = ebosSimulator_.vanguard().eclState(); phase_usage_ = phaseUsageFromDeck(eclState); active_.resize(phase_usage_.MaxNumPhases, false); @@ -38,17 +38,17 @@ namespace Opm { void BlackoilAquiferModel:: beginTimeStep() { - // Right now it doesn't do shit. + } // called at the end of a time step template void - BlackoilAquiferModel:: timeStepSucceeded() + BlackoilAquiferModel:: timeStepSucceeded(const SimulatorTimerInterface& timer) { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->after_time_step(); + aquifer->after_time_step(timer); } } @@ -57,7 +57,7 @@ namespace Opm { void BlackoilAquiferModel:: beginReportStep(const int time_step) { - // Right now it doesn't do shit. + } // called at the end of a report step @@ -65,14 +65,7 @@ namespace Opm { void BlackoilAquiferModel:: endReportStep() { - // Right now it just spits out the constants for each aquifers - // We are using the simple integer indexing for the aquifers - for (int i = 0; i < numAquifers(); ++i) - { - std::cout << "Aquifer[" << i << "]" - << " : Tc = " << aquifers()[i].time_constant() - << ", beta = " << aquifers()[i].aquifer_influx_constant() << std::endl; - } + } // Get the last report step @@ -80,9 +73,6 @@ namespace Opm { const SimulatorReport& BlackoilAquiferModel:: lastReport() const { - for (auto i = aquifers_.begin(); i != aquifers_.end(); ++i){ - (*i).print_private_members(); - } return last_report_; } @@ -93,7 +83,6 @@ namespace Opm { const int iterationIdx ) { last_report_ = SimulatorReport(); - // We need to update the reservoir pressures connected to the aquifer updateConnectionIntensiveQuantities(); @@ -103,35 +92,14 @@ namespace Opm { prepareTimeStep(timer); } - if (iterationIdx == 0) { - calculateExplicitQuantities(); - } - if (param_.solve_aquifereq_initially_ && iterationIdx == 0) { // solve the aquifer equations as a pre-processing step last_report_ = solveAquiferEq(timer); } - assembleAquiferEq(timer); - last_report_.converged = true; } - // Protected function: Update the primary variables - template - void - BlackoilAquiferModel:: updatePrimaryVariables() - { - // Right now it doesn't do shit. - } - - // Protected function: Init the primary variables - template - void - BlackoilAquiferModel:: initPrimaryVariablesEvaluation() const - { - // Right now it doesn't do shit. - } template void @@ -149,17 +117,6 @@ namespace Opm { } } - template - void - BlackoilAquiferModel:: calculateExplicitQuantities() - { - // for (auto aqui = aquifers_.begin(); aqui!= aquifers_.end(); ++aqui) - // { - // std::cout << "calculateExplicitQuantities: Aquifer id = " << aqui->aquiferID() << std::endl; - // aqui->calculateExplicitQuantities(ebosSimulator_); - // } - } - template SimulatorReport @@ -200,30 +157,10 @@ namespace Opm { int BlackoilAquiferModel:: numPhases() const { - // Not implemented yet!!!!!!!!!!!! const auto& pu = phase_usage_; return pu.num_phases; } - - // Protected function: returns the phase index in ebos - template - int - BlackoilAquiferModel:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const auto& pu = phase_usage_; - if (active_[Water] && pu.phase_pos[Water] == phaseIdx) - return FluidSystem::waterPhaseIdx; - if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) - return FluidSystem::oilPhaseIdx; - if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) - return FluidSystem::gasPhaseIdx; - - assert(phaseIdx < 3); - // for other phases return the index - return phaseIdx; - } - // Protected function which calls the individual aquifer models template void @@ -231,7 +168,6 @@ namespace Opm { { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - std::cout << "assembleAquiferEq: Aquifer id = " << aquifer->aquiferID() << std::endl; aquifer->assembleAquiferEq(ebosSimulator_, timer); } } @@ -261,47 +197,27 @@ namespace Opm { // Initialize the aquifers in the deck template void - BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers)//, std::vector< AquiferCarterTracy >& aquifers) + BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers) { - const auto& deck = ebosSimulator.gridManager().deck(); - const auto& eclState = ebosSimulator.gridManager().eclState(); + updateConnectionIntensiveQuantities(); + const auto& deck = ebosSimulator.vanguard().deck(); + const auto& eclState = ebosSimulator.vanguard().eclState(); // Get all the carter tracy aquifer properties data and put it in aquifers vector AquiferCT aquiferct = AquiferCT(eclState,deck); + Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); std::vector aquifersData = aquiferct.getAquifers(); - std::vector aquanconData = aquiferct.getAquancon(); + std::vector aquifer_connection = aquifer_connect.getAquOutput(); - // for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) - // { - - // } + assert( aquifersData.size() == aquifer_connect.size() ); - auto ita = aquifersData.cbegin(); - auto f_lambda = [&] (AquiferCT::AQUANCON_data i) { - aquifers.push_back( AquiferCarterTracy (*ita++, i, numComponents(), gravity_ ) ); - }; - std::for_each( aquanconData.cbegin(), aquanconData.cend(), f_lambda ); - } - // Begin the hack to initialize the aquifers in the deck - template - std::vector< AquiferCarterTracy > - BlackoilAquiferModel:: hack_init(const Simulator& ebosSimulator)//, std::vector< AquiferCarterTracy >& aquifers) - { - std::vector< AquiferCarterTracy > aquifers; - /** Begin hack!!!!! */ - const auto& deck = ebosSimulator.gridManager().deck(); - const auto& eclState = ebosSimulator.gridManager().eclState(); - - // Get all the carter tracy aquifer properties data and put it in aquifers vector - AquiferCT aquiferct = AquiferCT(eclState,deck); - - std::vector aquifersData = aquiferct.getAquifers(); - - for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) + for (int i = 0; i < aquifersData.size(); ++i) { - aquifers.push_back( AquiferCarterTracy (*aquiferData, numComponents(), gravity_ ) ); + aquifers.push_back( + AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), numComponents(), gravity_, ebosSimulator_) + ); } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index c55f337bd..2408bb3cd 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -29,10 +29,10 @@ #include #include -#include #include #include #include +#include #include #include @@ -350,12 +350,12 @@ namespace Opm { const ReservoirState& reservoir_state, WellState& well_state) { - DUNE_UNUSED_PARAMETER(timer); + // DUNE_UNUSED_PARAMETER(timer); DUNE_UNUSED_PARAMETER(reservoir_state); DUNE_UNUSED_PARAMETER(well_state); wellModel().timeStepSucceeded(); - aquiferModel().timeStepSucceeded(); + aquiferModel().timeStepSucceeded(timer); ebosSimulator_.problem().endTimeStep(); } @@ -429,13 +429,13 @@ namespace Opm { if (elem.partitionType() != Dune::InteriorEntity) continue; - unsigned globalElemIdx = elemMapper.index(elem); + unsigned globalElemIdx = elemMapper.index(elem); const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; Scalar pressureNew; - pressureNew = priVarsNew[Indices::pressureSwitchIdx]; + pressureNew = priVarsNew[Indices::pressureSwitchIdx]; - Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; + Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationNew = 1.0; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx]; @@ -478,7 +478,7 @@ namespace Opm { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; - resultDelta += tmp*tmp; + resultDelta += tmp*tmp; resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; } } @@ -486,9 +486,9 @@ namespace Opm { resultDelta = gridView.comm().sum(resultDelta); resultDenom = gridView.comm().sum(resultDenom); - if (resultDenom > 0.0) - return resultDelta/resultDenom; - return 0.0; + if (resultDenom > 0.0) + return resultDelta/resultDenom; + return 0.0; } @@ -599,7 +599,6 @@ namespace Opm { virtual void apply( const X& x, Y& y ) const { A_.mv( x, y ); - // add well model modification to y wellMod_.apply(x, y ); @@ -613,7 +612,6 @@ namespace Opm { virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const { A_.usmv(alpha,x,y); - // add scaled well model modification to y wellMod_.applyScaleAdd( alpha, x, y ); @@ -1129,29 +1127,14 @@ namespace Opm { const BlackoilAquiferModel& aquiferModel() const { return aquifer_model_; } - int flowPhaseToEbosCompIdx( const int phaseIdx ) const + int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const { - const auto& pu = phaseUsage_; - if (active_[Water] && pu.phase_pos[Water] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - - // for other phases return the index - return phaseIdx; - } - - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const auto& pu = phaseUsage_; - if (active_[Water] && pu.phase_pos[Water] == phaseIdx) - return FluidSystem::waterPhaseIdx; - if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) - return FluidSystem::oilPhaseIdx; - if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) - return FluidSystem::gasPhaseIdx; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx) + return Water; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx) + return Oil; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx) + return Gas; assert(phaseIdx < 3); // for other phases return the index @@ -1170,7 +1153,6 @@ namespace Opm { private: - double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } double drMaxRel() const { return param_.dr_max_rel_; } @@ -1181,4 +1163,4 @@ namespace Opm { }; } // namespace Opm -#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED +#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED \ No newline at end of file diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 5e5988b22..cb0c66399 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -190,8 +190,8 @@ public: auto auxMod = std::make_shared >(schedule(), grid()); ebosSimulator_.model().addAuxiliaryModule(auxMod); } + AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_); - // aquifer_model.hack_init(ebosSimulator_); // Main simulation loop. while (!timer.done()) { From 0d55b7aaced72709ef8ba8001c86e0a444563859 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Wed, 28 Mar 2018 20:43:35 +0200 Subject: [PATCH 09/31] Removes Eigen from interpolation of influence tables --- opm/autodiff/AquiferCarterTracy.hpp | 48 +++++------------------------ 1 file changed, 8 insertions(+), 40 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 4ba7597c4..3af93236d 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -21,11 +21,11 @@ #ifndef OPM_AQUIFERCT_HEADER_INCLUDED #define OPM_AQUIFERCT_HEADER_INCLUDED -#include #include #include #include #include +#include #include @@ -223,40 +223,14 @@ namespace Opm Scalar dt_, pa0_, gravity_; bool p0_defaulted_; Eval W_flux_; - // Also return the polynomial fit - std::vector coeff_; + - // We fit the tabular data using a polynomial fit - // Modified from Copyright (C) 2014 Clifford Wolf - // http://svn.clifford.at/handicraft/2014/polyfit/polyfit.cc - inline void polynomial_fit( const std::vector &X, const std::vector &y, - std::vector &coeff, int order, bool bias) const - { - size_t colNum = (bias)? order + 1 : order; - Eigen::MatrixXd A(X.size(), colNum); - Eigen::VectorXd y_mapped = Eigen::VectorXd::Map(&y.front(), y.size()); - Eigen::VectorXd result; - - assert(X.size() == y.size()); - assert(X.size() >= colNum); - - // create matrix - for (size_t i = 0; i < X.size(); i++) - for (size_t j = 0; j < colNum; j++) - A(i, j) = (bias)? pow(X.at(i), j) : pow(X.at(i), j+1); - - // solve for linear least squares fit - result = A.householderQr().solve(y_mapped); - - coeff.resize(colNum); - for (size_t i = 0; i < colNum; i++) - coeff[i] = result[i]; - } - inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { - // http://kluge.in-chemnitz.de/opensource/spline/ + // We use the opm-common numeric linear interpolator + pitd = Opm::linearInterpolation(aqutab_td_, aqutab_pi_, td); + pitd_prime = Opm::linearInterpolationDerivative(aqutab_td_, aqutab_pi_, td); } inline void init_quantities(const Aquancon::AquanconOutput& connection) @@ -272,8 +246,6 @@ namespace Opm pressure_previous_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.); Qai_.resize(cell_idx_.size(), 0.0); - - polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 1, true); } inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) @@ -301,9 +273,9 @@ namespace Opm Scalar Tc = time_constant(); Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; Scalar td = timer.simulationTimeElapsed() / Tc; - - Scalar PItdprime = coeff_.at(1); - Scalar PItd = coeff_.at(0) + coeff_.at(1)*td_plus_dt; + Scalar PItdprime = 0.; + Scalar PItd = 0.; + get_influence_table_values(PItd, PItdprime, td_plus_dt); a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); b = beta / (Tc * ( PItd - td*PItdprime)); } @@ -418,11 +390,7 @@ namespace Opm temperature_aq = fs_aquifer.temperature(0); pa0_mean = pa0_; - // rho_mean = FluidSystem::referenceDensity(waterPhaseIdx, pvttableIdx) - // *FluidSystem::waterPvt().inverseFormationVolumeFactor(pvttableIdx, temperature_aq, pa0_mean); - Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); - std::cout << "Pa0 = " << pa0_mean << ", viscosity = " << mu_w_aquifer.value() << std::endl; mu_w_ = mu_w_aquifer.value(); From ec680664fd6a13c9b51080ab462a9e0c74bc000c Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 5 Apr 2018 15:10:10 +0200 Subject: [PATCH 10/31] Implements the correct check for aquifer face connections to make sure face is connected to the boundary --- opm/autodiff/AquiferCarterTracy.hpp | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 3af93236d..d1a655df9 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -308,6 +308,14 @@ namespace Opm const auto& grid = eclState.getInputGrid(); cell_idx_ = connection.global_index; + auto globalCellIdx = ugrid.globalCell(); + // for (auto globalCells : globalCellIdx){ + // std::cout << "global id = " << globalCells << std::endl; + // } + + // for (auto cellidx : cell_idx_){ + // std::cout << "aqucell id = " << cellidx << std::endl; + // } assert( cell_idx_ == connection.global_index); assert( (cell_idx_.size() == connection.influx_coeff.size()) ); @@ -318,10 +326,14 @@ namespace Opm cell_depth_.resize(cell_idx_.size(), d0_); alphai_.resize(cell_idx_.size(), 1.0); faceArea_connected_.resize(cell_idx_.size(),0.0); + Scalar faceArea; auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); + + // static_assert(decltype(faceCells)::dummy_error, "DUMP MY TYPE" ); + // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; @@ -355,7 +367,11 @@ namespace Opm if (faceDirection == connection.reservoir_face_dir.at(idx)) { - faceArea_connected_.at(idx) = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); + // Check now if the face is outside of the reservoir, or if it adjoins an inactive cell + // Do not make the connection if the product of the two cellIdx > 0. This is because the + // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) + faceArea = (faceCells(faceIdx,0)*faceCells(faceIdx,1) > 0)? 0. : Opm::UgGridHelpers::faceArea(ugrid, faceIdx); + faceArea_connected_.at(idx) = faceArea; denom_face_areas += faceArea_connected_.at(idx); } } From 8d3d7ccc2ffb48fc551c19f2c5fca003fa1829fd Mon Sep 17 00:00:00 2001 From: kel85uk Date: Tue, 10 Apr 2018 13:05:30 +0200 Subject: [PATCH 11/31] Cleans up the code and removes unused functions --- opm/autodiff/AquiferCarterTracy.hpp | 40 ++------- opm/autodiff/BlackoilAquiferModel.hpp | 27 +----- opm/autodiff/BlackoilAquiferModel_impl.hpp | 86 +------------------ opm/autodiff/BlackoilModelEbos.hpp | 15 ---- opm/autodiff/BlackoilModelParameters.cpp | 2 - opm/autodiff/BlackoilModelParameters.hpp | 3 - opm/autodiff/BlackoilTransportModel.hpp | 1 - .../SimulatorFullyImplicitBlackoilEbos.hpp | 5 -- 8 files changed, 7 insertions(+), 172 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index d1a655df9..2d4c35dcb 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -69,21 +69,12 @@ namespace Opm typedef DenseAd::Evaluation Eval; typedef Opm::BlackOilFluidState FluidState; - typedef typename FluidSystem::WaterPvt WaterPvt; - - typedef Ewoms::BlackOilPolymerModule PolymerModule; - static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); - static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); - static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; - static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; - - explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, - const int numComponents, const Scalar gravity, const Simulator& ebosSimulator ) + const Scalar gravity, const Simulator& ebosSimulator ) : phi_aq_ (params.phi_aq), // C_t_ (params.C_t), // r_o_ (params.r_o), // @@ -100,22 +91,12 @@ namespace Opm aquiferID_ (params.aquiferID), inftableID_ (params.inftableID), pvttableID_ (params.pvttableID), - num_components_ (numComponents), gravity_ (gravity), ebos_simulator_ (ebosSimulator) { init_quantities(connection); } - inline const PhaseUsage& - phaseUsage() const - { - assert(phase_usage_); - - return *phase_usage_; - } - - inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) { dt_ = timer.currentStepLength(); @@ -166,10 +147,6 @@ namespace Opm } } - inline const double area_fraction(const size_t i) - { - return alphai_.at(i); - } inline const std::vector cell_id() const { @@ -183,13 +160,11 @@ namespace Opm private: - const PhaseUsage* phase_usage_; const Simulator& ebos_simulator_; // Aquifer ID, and other IDs int aquiferID_, inftableID_, pvttableID_; - int num_components_; // Grid variables @@ -225,6 +200,10 @@ namespace Opm Eval W_flux_; + inline const double area_fraction(const size_t i) + { + return alphai_.at(i); + } inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { @@ -309,13 +288,6 @@ namespace Opm cell_idx_ = connection.global_index; auto globalCellIdx = ugrid.globalCell(); - // for (auto globalCells : globalCellIdx){ - // std::cout << "global id = " << globalCells << std::endl; - // } - - // for (auto cellidx : cell_idx_){ - // std::cout << "aqucell id = " << cellidx << std::endl; - // } assert( cell_idx_ == connection.global_index); assert( (cell_idx_.size() == connection.influx_coeff.size()) ); @@ -332,8 +304,6 @@ namespace Opm auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); - // static_assert(decltype(faceCells)::dummy_error, "DUMP MY TYPE" ); - // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index 706756ecb..485c3dca6 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -92,19 +92,9 @@ namespace Opm { void assemble( const SimulatorTimerInterface& timer, const int iterationIdx ); - // called at the beginning of a time step - void beginTimeStep(); // called at the end of a time step void timeStepSucceeded(const SimulatorTimerInterface& timer); - // called at the beginning of a report step - void beginReportStep(const int time_step); - - // called at the end of a report step - void endReportStep(); - - const SimulatorReport& lastReport() const; - inline const Simulator& simulator() const { return ebosSimulator_; @@ -119,11 +109,7 @@ namespace Opm { const ModelParameters param_; bool terminal_output_; - bool has_solvent_; - bool has_polymer_; - std::vector pvt_region_idx_; - PhaseUsage phase_usage_; - std::vector active_; + size_t global_nc_; // the number of the cells in the local grid size_t number_of_cells_; @@ -132,23 +118,12 @@ namespace Opm { std::vector aquifers_; - SimulatorReport last_report_; - - void updateConnectionIntensiveQuantities() const; - // The number of components in the model. - int numComponents() const; - int numAquifers() const; - int numPhases() const; - void assembleAquiferEq(const SimulatorTimerInterface& timer); - SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer); - - // some preparation work, mostly related to group control and RESV, // at the beginning of each time step (Not report step) void prepareTimeStep(const SimulatorTimerInterface& timer); diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index c75b3f53b..3087d951e 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -9,22 +9,11 @@ namespace Opm { : ebosSimulator_(ebosSimulator) , param_(param) , terminal_output_(terminal_output) - , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) - , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) { const auto& eclState = ebosSimulator_.vanguard().eclState(); - phase_usage_ = phaseUsageFromDeck(eclState); - - active_.resize(phase_usage_.MaxNumPhases, false); - for (int p = 0; p < phase_usage_.MaxNumPhases; ++p) { - active_[ p ] = phase_usage_.phase_used[ p ] != 0; - } const auto& gridView = ebosSimulator_.gridView(); - // calculate the number of elements of the compressed sequential grid. this needs - // to be done in two steps because the dune communicator expects a reference as - // argument for sum() number_of_cells_ = gridView.size(/*codim=*/0); global_nc_ = gridView.comm().sum(number_of_cells_); gravity_ = ebosSimulator_.problem().gravity()[2]; @@ -32,15 +21,6 @@ namespace Opm { } - - // called at the beginning of a time step - template - void - BlackoilAquiferModel:: beginTimeStep() - { - - } - // called at the end of a time step template void @@ -52,37 +32,12 @@ namespace Opm { } } - // called at the beginning of a report step - template - void - BlackoilAquiferModel:: beginReportStep(const int time_step) - { - - } - - // called at the end of a report step - template - void - BlackoilAquiferModel:: endReportStep() - { - - } - - // Get the last report step - template - const SimulatorReport& - BlackoilAquiferModel:: lastReport() const - { - return last_report_; - } - template void BlackoilAquiferModel:: assemble( const SimulatorTimerInterface& timer, const int iterationIdx ) { - last_report_ = SimulatorReport(); // We need to update the reservoir pressures connected to the aquifer updateConnectionIntensiveQuantities(); @@ -92,12 +47,7 @@ namespace Opm { prepareTimeStep(timer); } - if (param_.solve_aquifereq_initially_ && iterationIdx == 0) { - // solve the aquifer equations as a pre-processing step - last_report_ = solveAquiferEq(timer); - } assembleAquiferEq(timer); - last_report_.converged = true; } @@ -118,32 +68,6 @@ namespace Opm { } - template - SimulatorReport - BlackoilAquiferModel:: solveAquiferEq(const SimulatorTimerInterface& timer) - { - // We need to solve the equilibrium equation first to - // obtain the initial pressure of water in the aquifer - SimulatorReport report = SimulatorReport(); - return report; - } - - // Protected function: Return number of components in the model. - template - int - BlackoilAquiferModel:: numComponents() const - { - if (numPhases() == 2) { - return 2; - } - int numComp = FluidSystem::numComponents; - if (has_solvent_) { - numComp ++; - } - - return numComp; - } - // Protected function: Return number of aquifers in the model. template int @@ -152,14 +76,6 @@ namespace Opm { return aquifers_.size(); } - // Protected function: Return number of phases in the model. - template - int - BlackoilAquiferModel:: numPhases() const - { - const auto& pu = phase_usage_; - return pu.num_phases; - } // Protected function which calls the individual aquifer models template @@ -216,7 +132,7 @@ namespace Opm { for (int i = 0; i < aquifersData.size(); ++i) { aquifers.push_back( - AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), numComponents(), gravity_, ebosSimulator_) + AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), gravity_, ebosSimulator_) ); } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 2408bb3cd..4d91fa59d 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -209,7 +209,6 @@ namespace Opm { wasSwitched_.resize(numDof); std::fill(wasSwitched_.begin(), wasSwitched_.end(), false); - aquiferModel().beginTimeStep(); wellModel().beginTimeStep(); if (param_.update_equations_scaling_) { @@ -1127,20 +1126,6 @@ namespace Opm { const BlackoilAquiferModel& aquiferModel() const { return aquifer_model_; } - int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const - { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx) - return Water; - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx) - return Oil; - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx) - return Gas; - - assert(phaseIdx < 3); - // for other phases return the index - return phaseIdx; - } - void beginReportStep() { ebosSimulator_.problem().beginEpisode(); diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index 2bd302c2d..f16fc7685 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -63,7 +63,6 @@ namespace Opm param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day ); max_strict_iter_ = param.getDefault("max_strict_iter",8); solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_); - solve_aquifereq_initially_ = param.getDefault("solve_aquifereq_initially",solve_aquifereq_initially_); update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_); use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_); deck_file_name_ = param.template get("deck_filename"); @@ -95,7 +94,6 @@ namespace Opm max_inner_iter_ms_wells_ = 10; maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day ); solve_welleq_initially_ = true; - solve_aquifereq_initially_ = true; update_equations_scaling_ = false; use_update_stabilization_ = true; use_multisegment_well_ = false; diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index d57552161..b01cfdde9 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -77,9 +77,6 @@ namespace Opm /// Solve well equation initially bool solve_welleq_initially_; - /// Solve aquifer equation initially - bool solve_aquifereq_initially_; - /// Update scaling factors for mass balance equations bool update_equations_scaling_; diff --git a/opm/autodiff/BlackoilTransportModel.hpp b/opm/autodiff/BlackoilTransportModel.hpp index d6610f0d5..ffc0ebc92 100644 --- a/opm/autodiff/BlackoilTransportModel.hpp +++ b/opm/autodiff/BlackoilTransportModel.hpp @@ -86,7 +86,6 @@ namespace Opm { { Base::prepareStep(timer, reservoir_state, well_state); Base::param_.solve_welleq_initially_ = false; - Base::param_.solve_aquifereq_initially_ = false; SolutionState state0 = variableState(reservoir_state, well_state); asImpl().makeConstantState(state0); asImpl().computeAccum(state0, 0); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index cb0c66399..bfbbc0e6b 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -208,8 +208,6 @@ public: well_model.beginReportStep(timer.currentStepNum()); - aquifer_model.beginReportStep(timer.currentStepNum()); - auto solver = createSolver(well_model, aquifer_model); // write the inital state at the report stage @@ -274,7 +272,6 @@ public: } solver->model().endReportStep(); - aquifer_model.endReportStep(); well_model.endReportStep(); // take time that was used to solve system for this reportStep @@ -318,8 +315,6 @@ public: report.total_time = total_timer.secsSinceStart(); report.converged = true; - auto reportaquifer = aquifer_model.lastReport(); - return report; } From 2170a00b21c11b4f3703dcd84b8168b1452674d3 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Wed, 11 Apr 2018 17:30:23 +0200 Subject: [PATCH 12/31] Adds the influx coefficient multipler to the face area fraction calculation. Also removes more unnecessary variables and reordered initialization list to prevent ordering warnings in gcc compiler --- opm/autodiff/AquiferCarterTracy.hpp | 45 ++++++++++++---------- opm/autodiff/BlackoilAquiferModel.hpp | 7 ---- opm/autodiff/BlackoilAquiferModel_impl.hpp | 10 ++--- 3 files changed, 29 insertions(+), 33 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 2d4c35dcb..9eabb301f 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -50,10 +50,6 @@ namespace Opm public: typedef BlackoilModelParameters ModelParameters; - static const int Water = BlackoilPhases::Aqua; - static const int Oil = BlackoilPhases::Liquid; - static const int Gas = BlackoilPhases::Vapour; - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; @@ -70,29 +66,31 @@ namespace Opm typedef DenseAd::Evaluation Eval; typedef Opm::BlackOilFluidState FluidState; + static const auto waterCompIdx = FluidSystem::waterCompIdx; + static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx; explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, const Scalar gravity, const Simulator& ebosSimulator ) - : phi_aq_ (params.phi_aq), // + : ebos_simulator_ (ebosSimulator), + aquiferID_ (params.aquiferID), + inftableID_ (params.inftableID), + pvttableID_ (params.pvttableID), + phi_aq_ (params.phi_aq), // + d0_ (params.d0), C_t_ (params.C_t), // r_o_ (params.r_o), // k_a_ (params.k_a), // c1_ (params.c1), h_ (params.h), // - p0_defaulted_ (params.p0_defaulted), - pa0_ (params.p0), theta_ (params.theta), // c2_ (params.c2), // - d0_ (params.d0), aqutab_td_ (params.td), aqutab_pi_ (params.pi), - aquiferID_ (params.aquiferID), - inftableID_ (params.inftableID), - pvttableID_ (params.pvttableID), + pa0_ (params.p0), gravity_ (gravity), - ebos_simulator_ (ebosSimulator) + p0_defaulted_ (params.p0_defaulted) { init_quantities(connection); } @@ -118,12 +116,12 @@ namespace Opm get_current_density_cell(rhow_,idx,intQuants); calculate_inflow_rate(idx, timer); qinflow = Qai_.at(idx); - ebosResid[*cellID][(FluidSystem::waterCompIdx)] -= qinflow.value(); + ebosResid[*cellID][waterCompIdx] -= qinflow.value(); for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[*cellID][*cellID][(FluidSystem::waterCompIdx)][pvIdx] -= qinflow.derivative(pvIdx); + ebosJac[*cellID][*cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx); } } } @@ -230,13 +228,13 @@ namespace Opm inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - pressure_water.at(idx) = fs.pressure(FluidSystem::waterPhaseIdx); + pressure_water.at(idx) = fs.pressure(waterPhaseIdx); } inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - rho_water.at(idx) = fs.density(FluidSystem::waterPhaseIdx); + rho_water.at(idx) = fs.density(waterPhaseIdx); } inline Scalar dpai(int idx) @@ -303,6 +301,13 @@ namespace Opm auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); + for (auto influxCoeff: connection.influx_coeff){ + std::cout << "influx_coeff = " << influxCoeff << std::endl; + } + + for (auto influxMult: connection.influx_multiplier){ + std::cout << "influx_multiplier = " << influxMult << std::endl; + } // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; @@ -342,7 +347,7 @@ namespace Opm // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) faceArea = (faceCells(faceIdx,0)*faceCells(faceIdx,1) > 0)? 0. : Opm::UgGridHelpers::faceArea(ugrid, faceIdx); faceArea_connected_.at(idx) = faceArea; - denom_face_areas += faceArea_connected_.at(idx); + denom_face_areas += ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) ); } } auto cellCenter = grid.getCellCenter(cell_idx_.at(idx)); @@ -351,7 +356,7 @@ namespace Opm for (size_t idx = 0; idx < cell_idx_.size(); ++idx) { - alphai_.at(idx) = faceArea_connected_.at(idx)/denom_face_areas; + alphai_.at(idx) = ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas; } } @@ -394,8 +399,8 @@ namespace Opm const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0)); const auto& fs = intQuants.fluidState(); - water_pressure_reservoir.push_back( fs.pressure(FluidSystem::waterPhaseIdx).value() ); - rho_water_reservoir.at(idx) = fs.density(FluidSystem::waterPhaseIdx); + water_pressure_reservoir.push_back( fs.pressure(waterPhaseIdx).value() ); + rho_water_reservoir.at(idx) = fs.density(waterPhaseIdx); pw_aquifer.push_back( (water_pressure_reservoir.at(idx) - rho_water_reservoir.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_))*area_fraction(idx) ); } diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index 485c3dca6..a2d8a5b59 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -78,9 +78,6 @@ namespace Opm { typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - static const int numEq = BlackoilIndices::numEq; - static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; - typedef AquiferCarterTracy Aquifer_object; BlackoilAquiferModel(Simulator& ebosSimulator, @@ -110,11 +107,7 @@ namespace Opm { const ModelParameters param_; bool terminal_output_; - size_t global_nc_; - // the number of the cells in the local grid - size_t number_of_cells_; double gravity_; - std::vector depth_; std::vector aquifers_; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index 3087d951e..be1493550 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -10,12 +10,10 @@ namespace Opm { , param_(param) , terminal_output_(terminal_output) { - const auto& eclState = ebosSimulator_.vanguard().eclState(); + // const auto& gridView = ebosSimulator_.gridView(); - const auto& gridView = ebosSimulator_.gridView(); - - number_of_cells_ = gridView.size(/*codim=*/0); - global_nc_ = gridView.comm().sum(number_of_cells_); + // number_of_cells_ = gridView.size(/*codim=*/0); + // global_nc_ = gridView.comm().sum(number_of_cells_); gravity_ = ebosSimulator_.problem().gravity()[2]; init(ebosSimulator_, aquifers_); } @@ -129,7 +127,7 @@ namespace Opm { assert( aquifersData.size() == aquifer_connect.size() ); - for (int i = 0; i < aquifersData.size(); ++i) + for (size_t i = 0; i < aquifersData.size(); ++i) { aquifers.push_back( AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), gravity_, ebosSimulator_) From 7f2c2444f0a9756663fcbc5bba36c24df8064b02 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 26 Apr 2018 16:22:26 +0200 Subject: [PATCH 13/31] Cleans up the Aquifer model classes. --- opm/autodiff/AquiferCarterTracy.hpp | 254 +++++++----------- opm/autodiff/BlackoilAquiferModel.hpp | 57 +--- opm/autodiff/BlackoilAquiferModel_impl.hpp | 56 +--- opm/autodiff/BlackoilModelEbos.hpp | 14 +- .../SimulatorFullyImplicitBlackoilEbos.hpp | 4 +- 5 files changed, 125 insertions(+), 260 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 9eabb301f..ff9841016 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -24,21 +24,14 @@ #include #include #include -#include #include -#include - #include #include #include -#include -#include #include #include -#include -#include namespace Opm { @@ -48,17 +41,11 @@ namespace Opm { public: - typedef BlackoilModelParameters ModelParameters; - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; - typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; static const int numEq = BlackoilIndices::numEq; typedef double Scalar; @@ -71,73 +58,56 @@ namespace Opm - explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, - const Scalar gravity, const Simulator& ebosSimulator ) - : ebos_simulator_ (ebosSimulator), - aquiferID_ (params.aquiferID), - inftableID_ (params.inftableID), - pvttableID_ (params.pvttableID), - phi_aq_ (params.phi_aq), // - d0_ (params.d0), - C_t_ (params.C_t), // - r_o_ (params.r_o), // - k_a_ (params.k_a), // - c1_ (params.c1), - h_ (params.h), // - theta_ (params.theta), // - c2_ (params.c2), // - aqutab_td_ (params.td), - aqutab_pi_ (params.pi), - pa0_ (params.p0), - gravity_ (gravity), - p0_defaulted_ (params.p0_defaulted) + AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data, + const Aquancon::AquanconOutput& connection, + Simulator& ebosSimulator ) + : ebos_simulator_ (ebosSimulator), + aquct_data_ (aquct_data), + gravity_ (ebos_simulator_.problem().gravity()[2]) { - init_quantities(connection); + initQuantities(connection); } - inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + inline void assembleAquiferEq(const SimulatorTimerInterface& timer) { - dt_ = timer.currentStepLength(); - auto& ebosJac = ebosSimulator.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator.model().linearizer().residual(); + auto& ebosJac = ebos_simulator_.model().linearizer().matrix(); + auto& ebosResid = ebos_simulator_.model().linearizer().residual(); - - auto cellID = cell_idx_.begin(); - - size_t idx; - for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) + size_t cellID; + for ( size_t idx = 0; idx < cell_idx_.size(); ++idx ) { Eval qinflow = 0.0; + cellID = cell_idx_.at(idx); // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to // IntensiveQuantities of that particular cell_id - const IntensiveQuantities intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + const IntensiveQuantities intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellID, /*timeIdx=*/ 0)); // This is the pressure at td + dt - get_current_Pressure_cell(pressure_current_,idx,intQuants); - get_current_density_cell(rhow_,idx,intQuants); - calculate_inflow_rate(idx, timer); + updateCellPressure(pressure_current_,idx,intQuants); + updateCellDensity(idx,intQuants); + calculateInflowRate(idx, timer); qinflow = Qai_.at(idx); - ebosResid[*cellID][waterCompIdx] -= qinflow.value(); + ebosResid[cellID][waterCompIdx] -= qinflow.value(); for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[*cellID][*cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx); + ebosJac[cellID][cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx); } } } - inline void before_time_step(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + inline void beforeTimeStep(const SimulatorTimerInterface& timer) { auto cellID = cell_idx_.begin(); size_t idx; for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) { - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); - get_current_Pressure_cell(pressure_previous_ ,idx,intQuants); + const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + updateCellPressure(pressure_previous_ ,idx,intQuants); } } - inline void after_time_step(const SimulatorTimerInterface& timer) + inline void afterTimeStep(const SimulatorTimerInterface& timer) { for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) { @@ -146,139 +116,117 @@ namespace Opm } - inline const std::vector cell_id() const - { - return cell_idx_; - } - - inline const int& aquiferID() const - { - return aquiferID_; - } - - private: - const Simulator& ebos_simulator_; - - - // Aquifer ID, and other IDs - int aquiferID_, inftableID_, pvttableID_; + Simulator& ebos_simulator_; // Grid variables - std::vector cell_idx_; std::vector faceArea_connected_; // Quantities at each grid id std::vector cell_depth_; - std::vector pressure_previous_; + std::vector pressure_previous_; std::vector pressure_current_; std::vector Qai_; std::vector rhow_; std::vector alphai_; // Variables constants - Scalar mu_w_ , //water viscosity - phi_aq_ , //aquifer porosity - d0_, // aquifer datum depth - C_t_ , //total compressibility - r_o_ , //aquifer inner radius - k_a_ , //aquifer permeability - c1_, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) - h_ , //aquifer thickness - theta_ , //angle subtended by the aquifer boundary - c2_ ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). + const AquiferCT::AQUCT_data aquct_data_; - // Variables for influence table - std::vector aqutab_td_, aqutab_pi_; + Scalar mu_w_ , //water viscosity + beta_ , // Influx constant + Tc_ , // Time constant + pa0_ , // initial aquifer pressure + gravity_ ; // gravitational acceleration - // Cumulative flux - Scalar dt_, pa0_, gravity_; - bool p0_defaulted_; Eval W_flux_; - - inline const double area_fraction(const size_t i) - { - return alphai_.at(i); - } - inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) + inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { // We use the opm-common numeric linear interpolator - pitd = Opm::linearInterpolation(aqutab_td_, aqutab_pi_, td); - pitd_prime = Opm::linearInterpolationDerivative(aqutab_td_, aqutab_pi_, td); + pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td); + pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td); } - inline void init_quantities(const Aquancon::AquanconOutput& connection) + inline void initQuantities(const Aquancon::AquanconOutput& connection) { // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 W_flux_ = 0.; // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function - initialize_connections(connection); + initializeConnections(connection); - calculate_aquifer_condition(); + calculateAquiferCondition(); + + calculateAquiferConstants(); pressure_previous_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.); Qai_.resize(cell_idx_.size(), 0.0); } - inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); pressure_water.at(idx) = fs.pressure(waterPhaseIdx); } - inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) + inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - rho_water.at(idx) = fs.density(waterPhaseIdx); + pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); + } + + inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + rhow_.at(idx) = fs.density(waterPhaseIdx); } inline Scalar dpai(int idx) { - Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_) - pressure_previous_.at(idx).value(); + Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0) - pressure_previous_.at(idx); return dp; } // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription - inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) + inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) { - Scalar beta = aquifer_influx_constant(); - Scalar Tc = time_constant(); - Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; - Scalar td = timer.simulationTimeElapsed() / Tc; + const Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc_; + const Scalar td = timer.simulationTimeElapsed() / Tc_; Scalar PItdprime = 0.; Scalar PItd = 0.; - get_influence_table_values(PItd, PItdprime, td_plus_dt); - a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); - b = beta / (Tc * ( PItd - td*PItdprime)); + getInfluenceTableValues(PItd, PItdprime, td_plus_dt); + a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); + b = beta_ / (Tc_ * ( PItd - td*PItdprime)); } // This function implements Eq 5.7 of the EclipseTechnicalDescription - inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer) + inline void calculateInflowRate(int idx, const SimulatorTimerInterface& timer) { Scalar a, b; - calculate_a_b_constants(a,b,idx,timer); - Qai_.at(idx) = area_fraction(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx).value() ) ); + calculateEqnConstants(a,b,idx,timer); + Qai_.at(idx) = alphai_.at(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx) ) ); } - inline const Scalar time_constant() const + inline void calculateAquiferConstants() { - Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); - return Tc; - } - - inline const Scalar aquifer_influx_constant() const - { - Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; - return beta; + // We calculate the influx constant + beta_ = aquct_data_.c2 * aquct_data_.h + * aquct_data_.theta * aquct_data_.phi_aq + * aquct_data_.C_t + * aquct_data_.r_o * aquct_data_.r_o; + // We calculate the time constant + Tc_ = mu_w_ * aquct_data_.phi_aq + * aquct_data_.C_t + * aquct_data_.r_o * aquct_data_.r_o + / ( aquct_data_.k_a * aquct_data_.c1 ); } // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer - inline void initialize_connections(const Aquancon::AquanconOutput& connection) + inline void initializeConnections(const Aquancon::AquanconOutput& connection) { const auto& eclState = ebos_simulator_.vanguard().eclState(); const auto& ugrid = ebos_simulator_.vanguard().grid(); @@ -293,7 +241,7 @@ namespace Opm assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) ); // We hack the cell depth values for now. We can actually get it from elementcontext pos - cell_depth_.resize(cell_idx_.size(), d0_); + cell_depth_.resize(cell_idx_.size(), aquct_data_.d0); alphai_.resize(cell_idx_.size(), 1.0); faceArea_connected_.resize(cell_idx_.size(),0.0); Scalar faceArea; @@ -301,14 +249,6 @@ namespace Opm auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); - for (auto influxCoeff: connection.influx_coeff){ - std::cout << "influx_coeff = " << influxCoeff << std::endl; - } - - for (auto influxMult: connection.influx_multiplier){ - std::cout << "influx_multiplier = " << influxMult << std::endl; - } - // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; @@ -326,19 +266,22 @@ namespace Opm // the logically-Cartesian direction of the face const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); - - if (faceTag == 0) // left - faceDirection = Opm::FaceDir::XMinus; - else if (faceTag == 1) // right - faceDirection = Opm::FaceDir::XPlus; - else if (faceTag == 2) // back - faceDirection = Opm::FaceDir::YMinus; - else if (faceTag == 3) // front - faceDirection = Opm::FaceDir::YPlus; - else if (faceTag == 4) // bottom - faceDirection = Opm::FaceDir::ZMinus; - else if (faceTag == 5) // top - faceDirection = Opm::FaceDir::ZPlus; + switch(faceTag) + { + case 0: faceDirection = Opm::FaceDir::XMinus; + break; + case 1: faceDirection = Opm::FaceDir::XPlus; + break; + case 2: faceDirection = Opm::FaceDir::YMinus; + break; + case 3: faceDirection = Opm::FaceDir::YPlus; + break; + case 4: faceDirection = Opm::FaceDir::ZMinus; + break; + case 5: faceDirection = Opm::FaceDir::ZPlus; + break; + default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined"); + } if (faceDirection == connection.reservoir_face_dir.at(idx)) { @@ -360,16 +303,20 @@ namespace Opm } } - inline void calculate_aquifer_condition() + inline void calculateAquiferCondition() { - int pvttableIdx = pvttableID_ - 1; + int pvttableIdx = aquct_data_.pvttableID - 1; rhow_.resize(cell_idx_.size(),0.); - if (p0_defaulted_) + if (aquct_data_.p0 < 1.0) { - pa0_ = calculate_reservoir_equilibrium(rhow_); + pa0_ = calculateReservoirEquilibrium(); + } + else + { + pa0_ = aquct_data_.p0; } // Initialize a FluidState object first @@ -388,10 +335,11 @@ namespace Opm } // This function is for calculating the aquifer properties from equilibrium state with the reservoir - inline Scalar calculate_reservoir_equilibrium(std::vector& rho_water_reservoir) + inline Scalar calculateReservoirEquilibrium() { // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices - std::vector water_pressure_reservoir, pw_aquifer; + std::vector pw_aquifer; + Scalar water_pressure_reservoir; for (size_t idx = 0; idx < cell_idx_.size(); ++idx) { @@ -399,9 +347,9 @@ namespace Opm const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0)); const auto& fs = intQuants.fluidState(); - water_pressure_reservoir.push_back( fs.pressure(waterPhaseIdx).value() ); - rho_water_reservoir.at(idx) = fs.density(waterPhaseIdx); - pw_aquifer.push_back( (water_pressure_reservoir.at(idx) - rho_water_reservoir.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_))*area_fraction(idx) ); + water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); + rhow_.at(idx) = fs.density(waterPhaseIdx); + pw_aquifer.push_back( (water_pressure_reservoir - rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0))*alphai_.at(idx) ); } // We take the average of the calculated equilibrium pressures. diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index a2d8a5b59..a10f4500e 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -24,40 +24,12 @@ #ifndef OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED #define OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED -#include - -#include -#include - -#include -#include - -#include - #include #include - -#include - #include - -#include -#include -#include -#include #include - -#include -#include -#include - -#include -#include -#include - #include - namespace Opm { /// Class for handling the blackoil well model. @@ -67,22 +39,13 @@ namespace Opm { public: // --------- Types --------- - typedef BlackoilModelParameters ModelParameters; - - - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef AquiferCarterTracy Aquifer_object; - BlackoilAquiferModel(Simulator& ebosSimulator, - const ModelParameters& param, - const bool terminal_output); + explicit BlackoilAquiferModel(Simulator& ebosSimulator); // compute the well fluxes and assemble them in to the reservoir equations as source terms // and in the well equations. @@ -92,36 +55,22 @@ namespace Opm { // called at the end of a time step void timeStepSucceeded(const SimulatorTimerInterface& timer); - inline const Simulator& simulator() const - { - return ebosSimulator_; - } - - // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy - void init(const Simulator& ebosSimulator, std::vector& aquifers); - protected: Simulator& ebosSimulator_; - const ModelParameters param_; - bool terminal_output_; - - double gravity_; std::vector aquifers_; + // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy + void init(); void updateConnectionIntensiveQuantities() const; - int numAquifers() const; - void assembleAquiferEq(const SimulatorTimerInterface& timer); // at the beginning of each time step (Not report step) void prepareTimeStep(const SimulatorTimerInterface& timer); - const std::vector& aquifers(); - }; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index be1493550..9671c1691 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -3,19 +3,10 @@ namespace Opm { template BlackoilAquiferModel:: - BlackoilAquiferModel(Simulator& ebosSimulator, - const ModelParameters& param, - const bool terminal_output) + BlackoilAquiferModel(Simulator& ebosSimulator) : ebosSimulator_(ebosSimulator) - , param_(param) - , terminal_output_(terminal_output) { - // const auto& gridView = ebosSimulator_.gridView(); - - // number_of_cells_ = gridView.size(/*codim=*/0); - // global_nc_ = gridView.comm().sum(number_of_cells_); - gravity_ = ebosSimulator_.problem().gravity()[2]; - init(ebosSimulator_, aquifers_); + init(); } @@ -26,7 +17,7 @@ namespace Opm { { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->after_time_step(timer); + aquifer->afterTimeStep(timer); } } @@ -65,16 +56,6 @@ namespace Opm { } } - - // Protected function: Return number of aquifers in the model. - template - int - BlackoilAquiferModel:: numAquifers() const - { - return aquifers_.size(); - } - - // Protected function which calls the individual aquifer models template void @@ -82,7 +63,7 @@ namespace Opm { { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->assembleAquiferEq(ebosSimulator_, timer); + aquifer->assembleAquiferEq(timer); } } @@ -95,43 +76,34 @@ namespace Opm { // Here we can ask each carter tracy aquifers to get the current previous time step's pressure for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->before_time_step(ebosSimulator_, timer); + aquifer->beforeTimeStep(timer); } } - // Protected function: Returns a reference to the aquifers members in the model - template - const std::vector< AquiferCarterTracy >& - BlackoilAquiferModel:: aquifers() - { - return aquifers_; - } - - // Initialize the aquifers in the deck template void - BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers) + BlackoilAquiferModel:: init() { updateConnectionIntensiveQuantities(); - const auto& deck = ebosSimulator.vanguard().deck(); - const auto& eclState = ebosSimulator.vanguard().eclState(); + const auto& deck = ebosSimulator_.vanguard().deck(); + const auto& eclState = ebosSimulator_.vanguard().eclState(); // Get all the carter tracy aquifer properties data and put it in aquifers vector - AquiferCT aquiferct = AquiferCT(eclState,deck); - Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); + const AquiferCT aquiferct = AquiferCT(eclState,deck); + const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); std::vector aquifersData = aquiferct.getAquifers(); std::vector aquifer_connection = aquifer_connect.getAquOutput(); - assert( aquifersData.size() == aquifer_connect.size() ); + assert( aquifersData.size() == aquifer_connection.size() ); for (size_t i = 0; i < aquifersData.size(); ++i) { - aquifers.push_back( - AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), gravity_, ebosSimulator_) - ); + aquifers_.push_back( + AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), ebosSimulator_) + ); } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 4d91fa59d..d8a8fcfde 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -30,8 +30,6 @@ #include #include #include -#include -#include #include #include #include @@ -349,7 +347,7 @@ namespace Opm { const ReservoirState& reservoir_state, WellState& well_state) { - // DUNE_UNUSED_PARAMETER(timer); + DUNE_UNUSED_PARAMETER(timer); DUNE_UNUSED_PARAMETER(reservoir_state); DUNE_UNUSED_PARAMETER(well_state); @@ -372,18 +370,18 @@ namespace Opm { ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - // -------- Well and aquifer common variables ---------- - double dt = timer.currentStepLength(); - + // -------- Current time step length ---------- + const double dt = timer.currentStepLength(); + // -------- Aquifer models ---------- try { // Modify the Jacobian and residuals according to the aquifer models aquiferModel().assemble(timer, iterationIdx); } - catch( const Dune::FMatrixError& e ) + catch( ... ) { - OPM_THROW(Opm::NumericalProblem,"Error when assembling aquifer models"); + OPM_THROW(Opm::NumericalIssue,"Error when assembling aquifer models"); } // -------- Well equations ---------- diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index bfbbc0e6b..674f14918 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -30,8 +30,6 @@ #include #include #include -#include -#include #include #include #include @@ -191,7 +189,7 @@ public: ebosSimulator_.model().addAuxiliaryModule(auxMod); } - AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_); + AquiferModel aquifer_model(ebosSimulator_); // Main simulation loop. while (!timer.done()) { From e99fb80ce33b7be22b21e24060146a634fa34797 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Tue, 29 May 2018 14:16:19 +0200 Subject: [PATCH 14/31] Corrects the whitespace again for blackoilmodelebos.hpp --- opm/autodiff/BlackoilModelEbos.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index d8a8fcfde..75b6c3401 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -426,13 +426,13 @@ namespace Opm { if (elem.partitionType() != Dune::InteriorEntity) continue; - unsigned globalElemIdx = elemMapper.index(elem); + unsigned globalElemIdx = elemMapper.index(elem); const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; Scalar pressureNew; - pressureNew = priVarsNew[Indices::pressureSwitchIdx]; + pressureNew = priVarsNew[Indices::pressureSwitchIdx]; - Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; + Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationNew = 1.0; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx]; @@ -475,7 +475,7 @@ namespace Opm { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; - resultDelta += tmp*tmp; + resultDelta += tmp*tmp; resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; } } @@ -483,9 +483,9 @@ namespace Opm { resultDelta = gridView.comm().sum(resultDelta); resultDenom = gridView.comm().sum(resultDenom); - if (resultDenom > 0.0) - return resultDelta/resultDenom; - return 0.0; + if (resultDenom > 0.0) + return resultDelta/resultDenom; + return 0.0; } From 677255055d2e6fa7bcff762b13ba90bbaa584805 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Wed, 30 May 2018 14:26:07 +0200 Subject: [PATCH 15/31] Cleans up the BlackoilModelEbos.hpp file to address the final comments. --- opm/autodiff/BlackoilModelEbos.hpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 75b6c3401..9f37945d0 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -369,9 +369,6 @@ namespace Opm { ebosSimulator_.problem().beginIteration(); ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - - // -------- Current time step length ---------- - const double dt = timer.currentStepLength(); // -------- Aquifer models ---------- try @@ -384,6 +381,9 @@ namespace Opm { OPM_THROW(Opm::NumericalIssue,"Error when assembling aquifer models"); } + // -------- Current time step length ---------- + const double dt = timer.currentStepLength(); + // -------- Well equations ---------- try @@ -596,6 +596,7 @@ namespace Opm { virtual void apply( const X& x, Y& y ) const { A_.mv( x, y ); + // add well model modification to y wellMod_.apply(x, y ); @@ -609,6 +610,7 @@ namespace Opm { virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const { A_.usmv(alpha,x,y); + // add scaled well model modification to y wellMod_.applyScaleAdd( alpha, x, y ); @@ -1121,9 +1123,6 @@ namespace Opm { BlackoilAquiferModel& aquiferModel() { return aquifer_model_; } - const BlackoilAquiferModel& - aquiferModel() const { return aquifer_model_; } - void beginReportStep() { ebosSimulator_.problem().beginEpisode(); From 477740ca38df4c6eeb7928b0ffc4a5ed57c96cea Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 31 May 2018 09:51:54 +0200 Subject: [PATCH 16/31] introducing aquiferActive() to reduce of the overhead for the models that does not include a aquifer model. --- opm/autodiff/BlackoilAquiferModel.hpp | 4 +++- opm/autodiff/BlackoilAquiferModel_impl.hpp | 24 ++++++++++++++++++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index a10f4500e..694223ea2 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -71,10 +71,12 @@ namespace Opm { // at the beginning of each time step (Not report step) void prepareTimeStep(const SimulatorTimerInterface& timer); + bool aquiferActive() const; + }; } // namespace Opm #include "BlackoilAquiferModel_impl.hpp" -#endif \ No newline at end of file +#endif diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index 9671c1691..e92c952cc 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -15,6 +15,10 @@ namespace Opm { void BlackoilAquiferModel:: timeStepSucceeded(const SimulatorTimerInterface& timer) { + if ( !aquiferActive() ) { + return; + } + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { aquifer->afterTimeStep(timer); @@ -27,6 +31,10 @@ namespace Opm { assemble( const SimulatorTimerInterface& timer, const int iterationIdx ) { + if ( !aquiferActive() ) { + return; + } + // We need to update the reservoir pressures connected to the aquifer updateConnectionIntensiveQuantities(); @@ -85,8 +93,13 @@ namespace Opm { void BlackoilAquiferModel:: init() { - updateConnectionIntensiveQuantities(); const auto& deck = ebosSimulator_.vanguard().deck(); + + if ( !deck.hasKeyword("AQUCT") ) { + return ; + } + + updateConnectionIntensiveQuantities(); const auto& eclState = ebosSimulator_.vanguard().eclState(); // Get all the carter tracy aquifer properties data and put it in aquifers vector @@ -107,4 +120,11 @@ namespace Opm { } } -} // namespace Opm \ No newline at end of file + template + bool + BlackoilAquiferModel:: aquiferActive() const + { + return !aquifers_.empty(); + } + +} // namespace Opm From 78caac2c2066dccc1736017822f7a39f6c59f9d8 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 2 Nov 2017 12:00:02 +0100 Subject: [PATCH 17/31] Added build directory in the .gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 585d691bc..3f98f891f 100644 --- a/.gitignore +++ b/.gitignore @@ -32,3 +32,5 @@ test_vec # emacs directory setting: .dir-locals.el + +build From 5a87818bb2522b1e1c68ad73bf399ce197a568d1 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Mon, 4 Dec 2017 22:26:53 +0100 Subject: [PATCH 18/31] First addition of the class BlackoilAquiferModel. --- opm/autodiff/AquiferCarterTracy.hpp | 402 ++++++++++++++++++ opm/autodiff/BlackoilAquiferModel.hpp | 184 ++++++++ opm/autodiff/BlackoilAquiferModel_impl.hpp | 287 +++++++++++++ opm/autodiff/BlackoilModelEbos.hpp | 60 ++- opm/autodiff/BlackoilModelParameters.cpp | 2 + opm/autodiff/BlackoilModelParameters.hpp | 3 + opm/autodiff/BlackoilTransportModel.hpp | 1 + .../SimulatorFullyImplicitBlackoilEbos.hpp | 17 +- 8 files changed, 953 insertions(+), 3 deletions(-) create mode 100644 opm/autodiff/AquiferCarterTracy.hpp create mode 100644 opm/autodiff/BlackoilAquiferModel.hpp create mode 100644 opm/autodiff/BlackoilAquiferModel_impl.hpp diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp new file mode 100644 index 000000000..b6f9cfa17 --- /dev/null +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -0,0 +1,402 @@ +/* + Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface + 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_AQUIFERCT_HEADER_INCLUDED +#define OPM_AQUIFERCT_HEADER_INCLUDED + +#include +#include +#include +#include +#include + + +#include +#include + +#include +#include +#include +#include +#include + +namespace Opm +{ + + template + class AquiferCarterTracy + { + + public: + typedef BlackoilModelParameters ModelParameters; + + static const int Water = BlackoilPhases::Aqua; + static const int Oil = BlackoilPhases::Liquid; + static const int Gas = BlackoilPhases::Vapour; + + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + typedef typename GridView::template Codim<0>::Entity Element; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + static const int numEq = BlackoilIndices::numEq; + typedef double Scalar; + + typedef DenseAd::Evaluation Eval; + + typedef Ewoms::BlackOilPolymerModule PolymerModule; + + static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); + static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); + static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; + static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; + + + AquiferCarterTracy(const std::vector& cell_id) + : phi_aq_ (1.0), // + C_t_ (1.0), // + r_o_ (1.0), // + k_a_ (1.0), // + c1_ (1.0), + h_ (1.0), // + theta_ (1.0), // + c2_ (1.0), // + d0_ (1.0), + cell_idx_ (cell_id) + { + mu_w_ = 1e-3; + aqutab_td_.push_back(1.0); + aqutab_pi_.push_back(1.0); + aquiferID_ = 1; + inftableID_ = 1; + pvttableID_ = 1; + init_quantities(); + } + + explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const AquiferCT::AQUANCON_data& aquanconParams, + const int numComponents, const Scalar gravity ) + : phi_aq_ (params.phi_aq), // + C_t_ (params.C_t), // + r_o_ (params.r_o), // + k_a_ (params.k_a), // + c1_ (params.c1), + h_ (params.h), // + theta_ (params.theta), // + c2_ (params.c2), // + d0_ (params.d0), + aqutab_td_ (params.td), + aqutab_pi_ (params.pi), + aquiferID_ (params.aquiferID), + inftableID_ (params.inftableID), + pvttableID_ (params.pvttableID), + cell_idx_ (params.cell_id), + num_components_ (numComponents), + gravity_ (gravity) + { + mu_w_ = 1e-3; + init_quantities(aquanconParams); + } + + inline const PhaseUsage& + phaseUsage() const + { + assert(phase_usage_); + + return *phase_usage_; + } + + inline int + flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage(); + if (pu.phase_pos[Water] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + if (pu.phase_pos[Oil] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (pu.phase_pos[Gas] == phaseIdx) + return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + + // for other phases return the index + return phaseIdx; + } + + inline int + flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage(); + if (pu.phase_pos[Water] == phaseIdx) { + return FluidSystem::waterPhaseIdx; + } + if (pu.phase_pos[Oil] == phaseIdx) { + return FluidSystem::oilPhaseIdx; + } + if (pu.phase_pos[Gas] == phaseIdx) { + return FluidSystem::gasPhaseIdx; + } + + assert(phaseIdx < 3); + // for other phases return the index + return phaseIdx; + } + + inline void calculateExplicitQuantities(const Simulator& ebosSimulator) + { + std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; + } + + inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + { + std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; + // resAqui_ = 0.0; + dt_ = timer.currentStepLength(); + auto& ebosJac = ebosSimulator.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + // TODO: it probably can be static member for StandardWell + const double volume = 0.002831684659200; // 0.1 cu ft; + + auto cellID = cell_idx_.begin(); + size_t idx; + for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) + { + Eval qinflow = 0.0; + // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to + // IntensiveQuantities of that particular cell_id + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + // This is the pressure at td + dt + get_current_Pressure_cell(pressure_current_,idx,intQuants); + get_current_density_cell(rhow_,idx,intQuants); + calculate_inflow_rate(idx, timer); + qinflow = Qai_[idx]; + ebosResid[*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)] -= qinflow.value(); + + for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) + { + // also need to consider the efficiency factor when manipulating the jacobians. + ebosJac[*cellID][*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)][pvIdx] -= qinflow.derivative(pvIdx); + } + std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ + // << " -> P_wat[t+dt] = " << pressure_current_[idx] << std::endl + << " Qai[t+dt] = " << Qai_[idx] << std::endl; + } + } + + inline void before_time_step(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + { + auto cellID = cell_idx_.begin(); + size_t idx; + for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) + { + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + get_current_Pressure_cell(pressure_previous_ ,idx,intQuants); + } + } + + inline void after_time_step() + { + for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) + { + W_flux_ += (*Qai); + } + std::cout << "Aquifer # " << aquiferID_ << ": My cumulative flux = " << W_flux_ << std::endl; + } + + /* Made into public for testing only!!!!!!. Must be protected */ + inline const Scalar time_constant() const + { + Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); + return Tc; + } + + /* Made into public for testing only!!!!!!. Must be protected */ + inline const Scalar aquifer_influx_constant() const + { + Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; + return beta; + } + + // This is another hack to get the face area only for SPE1. + // Ideally it should be a map which given a cell_id, it returns the area fraction + inline const double area_fraction(const int i) + { + return 1000.0*20.0*0.092903/(1000.0*1000.0*0.092903*2 + 1000.0*20.0*0.092903*4); + } + + inline void print_private_members() const + { + std::cout << "Aquifer CT #" << aquiferID_ << std::endl; + auto ita = aqutab_td_.cbegin(); + auto f_lambda = [&ita] (double i) {std::cout << *ita++ << " " << i << std::endl;}; + std::for_each( aqutab_pi_.cbegin(), aqutab_pi_.cend(), f_lambda ); + + for (auto i = coeff_.begin(); i != coeff_.end(); ++i ) + { + std::cout << "Coeff = " << *i << std::endl; + } + } + + /* Made into public for testing only!!!!!!. Must be protected */ + inline const std::vector cell_id() const + { + return cell_idx_; + } + + inline const int& aquiferID() const + { + return aquiferID_; + } + + + + protected: + const PhaseUsage* phase_usage_; + + + // Aquifer ID, and other IDs + int aquiferID_, inftableID_, pvttableID_; + int num_components_; + + // Grid variables + std::vector cell_idx_; + + // Quantities at each grid id + std::vector cell_depth_; + std::vector pressure_previous_; + std::vector pressure_current_; + std::vector Qai_; + std::vector rhow_; + + // Variables constants + Scalar mu_w_ , //water viscosity + phi_aq_ , //aquifer porosity + d0_, // aquifer datum depth + C_t_ , //total compressibility + r_o_ , //aquifer inner radius + k_a_ , //aquifer permeability + c1_, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) + h_ , //aquifer thickness + theta_ , //angle subtended by the aquifer boundary + c2_ ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). + + // Variables for influence table + std::vector aqutab_td_, aqutab_pi_; + + // Cumulative flux + Scalar W_flux_, dt_, pa0_, gravity_; + + // Also return the polynomial fit + std::vector coeff_; + + + // We fit the tabular data using a polynomial fit + // Modified from Copyright (C) 2014 Clifford Wolf + // http://svn.clifford.at/handicraft/2014/polyfit/polyfit.cc + inline void polynomial_fit( const std::vector &X, const std::vector &y, + std::vector &coeff, int order, bool bias) const + { + int colNum = (bias)? order + 1 : order; + Eigen::MatrixXd A(X.size(), colNum); + Eigen::VectorXd y_mapped = Eigen::VectorXd::Map(&y.front(), y.size()); + Eigen::VectorXd result; + + assert(X.size() == y.size()); + assert(X.size() >= colNum); + + // create matrix + for (size_t i = 0; i < X.size(); i++) + for (size_t j = 0; j < colNum; j++) + A(i, j) = (bias)? pow(X.at(i), j) : pow(X.at(i), j+1); + + // solve for linear least squares fit + result = A.householderQr().solve(y_mapped); + + coeff.resize(colNum); + for (size_t i = 0; i < colNum; i++) + coeff[i] = result[i]; + } + + inline void init_quantities(const AquiferCT::AQUANCON_data& aquanconParams) + { + W_flux_ = 0.; + // pa0_ is the initial aquifer water pressure. Must be calculated from equilibrium if left default, + // or we get the information from the deck Hacked to make it at 45e6 Pa + pa0_ = 45e6; + + pressure_previous_.resize(cell_idx_.size(), 0.); + pressure_current_.resize(cell_idx_.size(), 0.); + // We hack the cell depth values for now. We can actually get it from elementcontext pos + cell_depth_.resize(cell_idx_.size(), d0_); + rhow_.resize(cell_idx_.size(), 998.0); + Qai_.resize(cell_idx_.size(), 0.); + + polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 2, true); + } + + inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + pressure_water[idx] = fs.pressure(FluidSystem::waterPhaseIdx).value(); + } + + inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + rho_water[idx] = fs.density(FluidSystem::waterPhaseIdx).value(); + } + + inline Scalar dpai(int idx) + { + Scalar dp = pa0_ - rhow_[idx]*gravity_*(cell_depth_[idx] - d0_) - pressure_previous_[idx]; + return dp; + } + + inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) + { + // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription + Scalar beta = aquifer_influx_constant(); + Scalar Tc = time_constant(); + Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; + Scalar td = timer.simulationTimeElapsed() / Tc; + Scalar PItdprime = coeff_[1] + 2.0*coeff_[2]*(td_plus_dt); + Scalar PItd = coeff_[0] + coeff_[1]*td_plus_dt + coeff_[2]*td_plus_dt*td_plus_dt; + a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_ * PItdprime) ) / ( PItd - td*PItdprime ); + b = beta / Tc / ( PItd - td*PItdprime); + } + + inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer) + { + Scalar a, b; + calculate_a_b_constants(a,b,idx,timer); + // This function implements Eq 5.7 of the EclipseTechnicalDescription + Qai_[idx] = area_fraction(idx)*( a - b * ( pressure_current_[idx] - pressure_previous_[idx] ) ); + } + + + }; // class AquiferCarterTracy + + +} // namespace Opm + +#endif \ No newline at end of file diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp new file mode 100644 index 000000000..a1ed97e8b --- /dev/null +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -0,0 +1,184 @@ +/* +<<<<<<< HEAD + File adapted from BlackoilWellModel.hpp + + Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface + Copyright 2017 Statoil ASA. +======= + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2017 IRIS AS +>>>>>>> 9ccee28... First addition of the class BlackoilAquiferModel. + + 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_BLACKOILAQUIFERMODEL_HEADER_INCLUDED +#define OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED + +#include + +#include +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + + +#include + + +namespace Opm { + + /// Class for handling the blackoil well model. + template + class BlackoilAquiferModel { + + public: + + + // --------- Types --------- + typedef BlackoilModelParameters ModelParameters; + + + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + + static const int numEq = BlackoilIndices::numEq; + static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; + + typedef Ewoms::BlackOilPolymerModule PolymerModule; + + typedef AquiferCarterTracy Aquifer_object; + + BlackoilAquiferModel(Simulator& ebosSimulator, + const ModelParameters& param, + const bool terminal_output); + + // compute the well fluxes and assemble them in to the reservoir equations as source terms + // and in the well equations. + void assemble( const SimulatorTimerInterface& timer, + const int iterationIdx ); + + // called at the beginning of a time step + void beginTimeStep(); + // called at the end of a time step + void timeStepSucceeded(); + + // called at the beginning of a report step + void beginReportStep(const int time_step); + + // called at the end of a report step + void endReportStep(); + + const SimulatorReport& lastReport() const; + + inline const Simulator& simulator() const + { + return ebosSimulator_; + } + + /// Hack function to get what I need from parser + void init(const Simulator& ebosSimulator, std::vector& aquifers); + + protected: + + Simulator& ebosSimulator_; + + const ModelParameters param_; + bool terminal_output_; + bool has_solvent_; + bool has_polymer_; + std::vector pvt_region_idx_; + PhaseUsage phase_usage_; + std::vector active_; + size_t global_nc_; + // the number of the cells in the local grid + size_t number_of_cells_; + double gravity_; + std::vector depth_; + std::vector aquifers_; + + + SimulatorReport last_report_; + + const Schedule& schedule() const + { return ebosSimulator_.gridManager().schedule(); } + + void updatePrimaryVariables(); + + void initPrimaryVariablesEvaluation() const; + + void updateConnectionIntensiveQuantities() const; + + void calculateExplicitQuantities(); + + // The number of components in the model. + int numComponents() const; + + int numAquifers() const; + + int numPhases() const; + + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; + + void assembleAquiferEq(const SimulatorTimerInterface& timer); + + SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer); + + // some preparation work, mostly related to group control and RESV, + // at the beginning of each time step (Not report step) + void prepareTimeStep(const SimulatorTimerInterface& timer); + + const std::vector& aquifers(); + + }; + + +} // namespace Opm + +#include "BlackoilAquiferModel_impl.hpp" +#endif diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp new file mode 100644 index 000000000..f9b2c4806 --- /dev/null +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -0,0 +1,287 @@ +namespace Opm { + + + template + BlackoilAquiferModel:: + BlackoilAquiferModel(Simulator& ebosSimulator, + const ModelParameters& param, + const bool terminal_output) + : ebosSimulator_(ebosSimulator) + , param_(param) + , terminal_output_(terminal_output) + , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) + , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) + { + const auto& eclState = ebosSimulator_.gridManager().eclState(); + phase_usage_ = phaseUsageFromDeck(eclState); + + active_.resize(phase_usage_.MaxNumPhases, false); + for (int p = 0; p < phase_usage_.MaxNumPhases; ++p) { + active_[ p ] = phase_usage_.phase_used[ p ] != 0; + } + + const auto& gridView = ebosSimulator_.gridView(); + + // calculate the number of elements of the compressed sequential grid. this needs + // to be done in two steps because the dune communicator expects a reference as + // argument for sum() + number_of_cells_ = gridView.size(/*codim=*/0); + global_nc_ = gridView.comm().sum(number_of_cells_); + gravity_ = ebosSimulator_.problem().gravity()[2]; + init(ebosSimulator_, aquifers_); + } + + + + // called at the beginning of a time step + template + void + BlackoilAquiferModel:: beginTimeStep() + { + // Right now it doesn't do shit. + } + + // called at the end of a time step + template + void + BlackoilAquiferModel:: timeStepSucceeded() + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + { + aquifer->after_time_step(); + } + } + + // called at the beginning of a report step + template + void + BlackoilAquiferModel:: beginReportStep(const int time_step) + { + // Right now it doesn't do shit. + } + + // called at the end of a report step + template + void + BlackoilAquiferModel:: endReportStep() + { + // Right now it just spits out the constants for each aquifers + // We are using the simple integer indexing for the aquifers + for (int i = 0; i < numAquifers(); ++i) + { + std::cout << "Aquifer[" << i << "]" + << " : Tc = " << aquifers()[i].time_constant() + << ", beta = " << aquifers()[i].aquifer_influx_constant() << std::endl; + } + } + + // Get the last report step + template + const SimulatorReport& + BlackoilAquiferModel:: lastReport() const + { + for (auto i = aquifers_.begin(); i != aquifers_.end(); ++i){ + (*i).print_private_members(); + } + return last_report_; + } + + template + void + BlackoilAquiferModel:: + assemble( const SimulatorTimerInterface& timer, + const int iterationIdx ) + { + last_report_ = SimulatorReport(); + + // We need to update the reservoir pressures connected to the aquifer + updateConnectionIntensiveQuantities(); + + if (iterationIdx == 0) { + // We can do the Table check and coefficients update in this function + // For now, it does nothing! + prepareTimeStep(timer); + } + + if (iterationIdx == 0) { + calculateExplicitQuantities(); + } + + if (param_.solve_aquifereq_initially_ && iterationIdx == 0) { + // solve the aquifer equations as a pre-processing step + last_report_ = solveAquiferEq(timer); + } + + assembleAquiferEq(timer); + + last_report_.converged = true; + } + + // Protected function: Update the primary variables + template + void + BlackoilAquiferModel:: updatePrimaryVariables() + { + // Right now it doesn't do shit. + } + + // Protected function: Init the primary variables + template + void + BlackoilAquiferModel:: initPrimaryVariablesEvaluation() const + { + // Right now it doesn't do shit. + } + + template + void + BlackoilAquiferModel:: updateConnectionIntensiveQuantities() const + { + ElementContext elemCtx(ebosSimulator_); + const auto& gridView = ebosSimulator_.gridView(); + const auto& elemEndIt = gridView.template end(); + for (auto elemIt = gridView.template begin(); + elemIt != elemEndIt; + ++elemIt) + { + elemCtx.updatePrimaryStencil(*elemIt); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + } + } + + template + void + BlackoilAquiferModel:: calculateExplicitQuantities() + { + // for (auto aqui = aquifers_.begin(); aqui!= aquifers_.end(); ++aqui) + // { + // std::cout << "calculateExplicitQuantities: Aquifer id = " << aqui->aquiferID() << std::endl; + // aqui->calculateExplicitQuantities(ebosSimulator_); + // } + } + + + template + SimulatorReport + BlackoilAquiferModel:: solveAquiferEq(const SimulatorTimerInterface& timer) + { + // We need to solve the equilibrium equation first to + // obtain the initial pressure of water in the aquifer + SimulatorReport report = SimulatorReport(); + return report; + } + + // Protected function: Return number of components in the model. + template + int + BlackoilAquiferModel:: numComponents() const + { + if (numPhases() == 2) { + return 2; + } + int numComp = FluidSystem::numComponents; + if (has_solvent_) { + numComp ++; + } + + return numComp; + } + + // Protected function: Return number of aquifers in the model. + template + int + BlackoilAquiferModel:: numAquifers() const + { + return aquifers_.size(); + } + + // Protected function: Return number of phases in the model. + template + int + BlackoilAquiferModel:: numPhases() const + { + // Not implemented yet!!!!!!!!!!!! + const auto& pu = phase_usage_; + return pu.num_phases; + } + + + // Protected function: returns the phase index in ebos + template + int + BlackoilAquiferModel:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const auto& pu = phase_usage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return FluidSystem::waterPhaseIdx; + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return FluidSystem::oilPhaseIdx; + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return FluidSystem::gasPhaseIdx; + + assert(phaseIdx < 3); + // for other phases return the index + return phaseIdx; + } + + // Protected function which calls the individual aquifer models + template + void + BlackoilAquiferModel:: assembleAquiferEq(const SimulatorTimerInterface& timer) + { + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + { + std::cout << "assembleAquiferEq: Aquifer id = " << aquifer->aquiferID() << std::endl; + aquifer->assembleAquiferEq(ebosSimulator_, timer); + } + } + + // Protected function + // some preparation work, mostly related to group control and RESV, + // at the beginning of each time step (Not report step) + template + void BlackoilAquiferModel:: prepareTimeStep(const SimulatorTimerInterface& timer) + { + // Here we can ask each carter tracy aquifers to get the current previous time step's pressure + for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) + { + aquifer->before_time_step(ebosSimulator_, timer); + } + } + + // Protected function: Returns a reference to the aquifers members in the model + template + const std::vector< AquiferCarterTracy >& + BlackoilAquiferModel:: aquifers() + { + return aquifers_; + } + + + // Initialize the aquifers in the deck + template + void + BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers)//, std::vector< AquiferCarterTracy >& aquifers) + { + const auto& deck = ebosSimulator.gridManager().deck(); + const auto& eclState = ebosSimulator.gridManager().eclState(); + + // Get all the carter tracy aquifer properties data and put it in aquifers vector + AquiferCT aquiferct = AquiferCT(eclState,deck); + + std::vector aquifersData = aquiferct.getAquifers(); + std::vector aquanconData = aquiferct.getAquancon(); + + // for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) + // { + + // } + + auto ita = aquifersData.cbegin(); + auto f_lambda = [&] (AquiferCT::AQUANCON_data i) { + aquifers.push_back( AquiferCarterTracy (*ita++, i, numComponents(), gravity_ ) ); + }; + std::for_each( aquanconData.cbegin(), aquanconData.cend(), f_lambda ); + } + +} // namespace Opm \ No newline at end of file diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index c16461496..c55f337bd 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -30,6 +30,9 @@ #include #include #include +#include +#include +#include #include #include @@ -143,6 +146,7 @@ namespace Opm { BlackoilModelEbos(Simulator& ebosSimulator, const ModelParameters& param, BlackoilWellModel& well_model, + BlackoilAquiferModel& aquifer_model, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output ) @@ -157,6 +161,7 @@ namespace Opm { , has_energy_(GET_PROP_VALUE(TypeTag, EnableEnergy)) , param_( param ) , well_model_ (well_model) + , aquifer_model_(aquifer_model) , terminal_output_ (terminal_output) , current_relaxation_(1.0) , dx_old_(UgGridHelpers::numCells(grid_)) @@ -204,6 +209,7 @@ namespace Opm { wasSwitched_.resize(numDof); std::fill(wasSwitched_.begin(), wasSwitched_.end(), false); + aquiferModel().beginTimeStep(); wellModel().beginTimeStep(); if (param_.update_equations_scaling_) { @@ -349,6 +355,7 @@ namespace Opm { DUNE_UNUSED_PARAMETER(well_state); wellModel().timeStepSucceeded(); + aquiferModel().timeStepSucceeded(); ebosSimulator_.problem().endTimeStep(); } @@ -366,9 +373,22 @@ namespace Opm { ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - // -------- Well equations ---------- + // -------- Well and aquifer common variables ---------- double dt = timer.currentStepLength(); + // -------- Aquifer models ---------- + try + { + // Modify the Jacobian and residuals according to the aquifer models + aquiferModel().assemble(timer, iterationIdx); + } + catch( const Dune::FMatrixError& e ) + { + OPM_THROW(Opm::NumericalProblem,"Error when assembling aquifer models"); + } + + // -------- Well equations ---------- + try { // assembles the well equations and applies the wells to @@ -1081,6 +1101,9 @@ namespace Opm { // Well Model BlackoilWellModel& well_model_; + // Aquifer Model + BlackoilAquiferModel& aquifer_model_; + /// \brief Whether we print something to std::cout bool terminal_output_; /// \brief The number of cells of the global grid. @@ -1100,6 +1123,41 @@ namespace Opm { const BlackoilWellModel& wellModel() const { return well_model_; } + BlackoilAquiferModel& + aquiferModel() { return aquifer_model_; } + + const BlackoilAquiferModel& + aquiferModel() const { return aquifer_model_; } + + int flowPhaseToEbosCompIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + + // for other phases return the index + return phaseIdx; + } + + int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const + { + const auto& pu = phaseUsage_; + if (active_[Water] && pu.phase_pos[Water] == phaseIdx) + return FluidSystem::waterPhaseIdx; + if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) + return FluidSystem::oilPhaseIdx; + if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) + return FluidSystem::gasPhaseIdx; + + assert(phaseIdx < 3); + // for other phases return the index + return phaseIdx; + } + void beginReportStep() { ebosSimulator_.problem().beginEpisode(); diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index f16fc7685..2bd302c2d 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -63,6 +63,7 @@ namespace Opm param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day ); max_strict_iter_ = param.getDefault("max_strict_iter",8); solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_); + solve_aquifereq_initially_ = param.getDefault("solve_aquifereq_initially",solve_aquifereq_initially_); update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_); use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_); deck_file_name_ = param.template get("deck_filename"); @@ -94,6 +95,7 @@ namespace Opm max_inner_iter_ms_wells_ = 10; maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day ); solve_welleq_initially_ = true; + solve_aquifereq_initially_ = true; update_equations_scaling_ = false; use_update_stabilization_ = true; use_multisegment_well_ = false; diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index b01cfdde9..d57552161 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -77,6 +77,9 @@ namespace Opm /// Solve well equation initially bool solve_welleq_initially_; + /// Solve aquifer equation initially + bool solve_aquifereq_initially_; + /// Update scaling factors for mass balance equations bool update_equations_scaling_; diff --git a/opm/autodiff/BlackoilTransportModel.hpp b/opm/autodiff/BlackoilTransportModel.hpp index ffc0ebc92..d6610f0d5 100644 --- a/opm/autodiff/BlackoilTransportModel.hpp +++ b/opm/autodiff/BlackoilTransportModel.hpp @@ -86,6 +86,7 @@ namespace Opm { { Base::prepareStep(timer, reservoir_state, well_state); Base::param_.solve_welleq_initially_ = false; + Base::param_.solve_aquifereq_initially_ = false; SolutionState state0 = variableState(reservoir_state, well_state); asImpl().makeConstantState(state0); asImpl().computeAccum(state0, 0); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 4445eb390..5e5988b22 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -29,6 +29,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -65,6 +68,7 @@ public: typedef BlackoilModelParameters ModelParameters; typedef NonlinearSolver Solver; typedef BlackoilWellModel WellModel; + typedef BlackoilAquiferModel AquiferModel; /// Initialise from parameters and objects to observe. @@ -186,6 +190,8 @@ public: auto auxMod = std::make_shared >(schedule(), grid()); ebosSimulator_.model().addAuxiliaryModule(auxMod); } + AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_); + // aquifer_model.hack_init(ebosSimulator_); // Main simulation loop. while (!timer.done()) { @@ -202,7 +208,9 @@ public: well_model.beginReportStep(timer.currentStepNum()); - auto solver = createSolver(well_model); + aquifer_model.beginReportStep(timer.currentStepNum()); + + auto solver = createSolver(well_model, aquifer_model); // write the inital state at the report stage if (timer.initialStep()) { @@ -266,6 +274,7 @@ public: } solver->model().endReportStep(); + aquifer_model.endReportStep(); well_model.endReportStep(); // take time that was used to solve system for this reportStep @@ -308,6 +317,9 @@ public: total_timer.stop(); report.total_time = total_timer.secsSinceStart(); report.converged = true; + + auto reportaquifer = aquifer_model.lastReport(); + return report; } @@ -320,11 +332,12 @@ public: protected: - std::unique_ptr createSolver(WellModel& well_model) + std::unique_ptr createSolver(WellModel& well_model, AquiferModel& aquifer_model) { auto model = std::unique_ptr(new Model(ebosSimulator_, model_param_, well_model, + aquifer_model, solver_, terminal_output_)); From ac7204fb2bf25be4c3cbfbe24410cd0d3f5f1e86 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Mon, 11 Dec 2017 17:38:45 +0100 Subject: [PATCH 19/31] Working and building skeleton which works for multiple aquifers. Also gets the properties of the aquifers using AQUCT keyword. --- opm/autodiff/AQUCT_params.hpp | 31 ++++++++++ opm/autodiff/BlackoilAquiferModel.hpp | 7 +-- opm/autodiff/BlackoilAquiferModel_impl.hpp | 66 ++++++++++++++++++++++ 3 files changed, 99 insertions(+), 5 deletions(-) create mode 100644 opm/autodiff/AQUCT_params.hpp diff --git a/opm/autodiff/AQUCT_params.hpp b/opm/autodiff/AQUCT_params.hpp new file mode 100644 index 000000000..da6c03ff5 --- /dev/null +++ b/opm/autodiff/AQUCT_params.hpp @@ -0,0 +1,31 @@ +#ifndef OPM_AQUCT_HEADER_INCLUDED +#define OPM_AQUCT_HEADER_INCLUDED + + + struct AQUCT_params{ + + // Aquifer ID + int aquiferID; + // Table IDs + int inftableID, pvttableID; + // Perforation cell id + int cell_id; + // Variables constants + double phi_aq , //aquifer porosity + d0, //aquifer datum depth + C_t , //total compressibility + r_o , //aquifer inner radius + k_a , //aquifer permeability + c1, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) + h , //aquifer thickness + theta , //angle subtended by the aquifer boundary + c2 ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). + std::vector td, pi; + + + }; + + + + +#endif \ No newline at end of file diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index a1ed97e8b..03771dced 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -1,15 +1,13 @@ /* -<<<<<<< HEAD File adapted from BlackoilWellModel.hpp Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface Copyright 2017 Statoil ASA. -======= Copyright 2016 SINTEF ICT, Applied Mathematics. Copyright 2016 - 2017 Statoil ASA. Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services Copyright 2016 - 2017 IRIS AS ->>>>>>> 9ccee28... First addition of the class BlackoilAquiferModel. + This file is part of the Open Porous Media project (OPM). @@ -62,7 +60,6 @@ #include #include - #include @@ -74,7 +71,6 @@ namespace Opm { public: - // --------- Types --------- typedef BlackoilModelParameters ModelParameters; @@ -156,6 +152,7 @@ namespace Opm { void calculateExplicitQuantities(); + // The number of components in the model. int numComponents() const; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index f9b2c4806..b9a861e44 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -284,4 +284,70 @@ namespace Opm { std::for_each( aquanconData.cbegin(), aquanconData.cend(), f_lambda ); } + // Begin the hack to initialize the aquifers in the deck + template + std::vector< AquiferCarterTracy > + BlackoilAquiferModel:: hack_init(const Simulator& ebosSimulator)//, std::vector< AquiferCarterTracy >& aquifers) + { + std::vector< AquiferCarterTracy > aquifers; + /** Begin hack!!!!! */ + const auto& deck = ebosSimulator.gridManager().deck(); + const auto& eclState = ebosSimulator.gridManager().eclState(); + + if (!deck.hasKeyword("AQUCT")){ + std::cout << "Nothing is shown! Where is AQUCT!????" << std::endl; + } + + const auto& aquctKeyword = deck.getKeyword("AQUCT"); + std::vector aquctParams; + // Resize the parameter vector container based on row entries in aquct + // We do the same for aquifers too because number of aquifers is assumed to be for each entry in aquct + aquctParams.resize(aquctKeyword.size()); + // aquifers.resize(aquctKeyword.size()); + + const int tableID = 0; + + std::cout << "Parsing AQUCT stuff" << std::endl; + for (size_t aquctRecordIdx = 0; aquctRecordIdx < aquctKeyword.size(); ++ aquctRecordIdx) + { + const auto& aquctRecord = aquctKeyword.getRecord(aquctRecordIdx); + + aquctParams.at(aquctRecordIdx).aquiferID = aquctRecord.getItem("AQUIFER_ID").template get(0); + aquctParams.at(aquctRecordIdx).h = aquctRecord.getItem("THICKNESS_AQ").template get(0); + aquctParams.at(aquctRecordIdx).phi_aq = aquctRecord.getItem("PORO_AQ").template get(0); + aquctParams.at(aquctRecordIdx).d0 = aquctRecord.getItem("DAT_DEPTH").getSIDouble(0); + aquctParams.at(aquctRecordIdx).C_t = aquctRecord.getItem("C_T").template get(0); + aquctParams.at(aquctRecordIdx).r_o = aquctRecord.getItem("RAD").getSIDouble(0); + aquctParams.at(aquctRecordIdx).k_a = aquctRecord.getItem("PERM_AQ").getSIDouble(0); + aquctParams.at(aquctRecordIdx).theta = aquctRecord.getItem("INFLUENCE_ANGLE").template get(0); + aquctParams.at(aquctRecordIdx).c1 = 0.008527; // We are using SI + aquctParams.at(aquctRecordIdx).c2 = 6.283; + aquctParams.at(aquctRecordIdx).inftableID = aquctRecord.getItem("TABLE_NUM_INFLUENCE_FN").template get(0) - 1; + aquctParams.at(aquctRecordIdx).pvttableID = aquctRecord.getItem("TABLE_NUM_WATER_PRESS").template get(0) - 1; + + std::cout << aquctParams.at(aquctRecordIdx).inftableID << std::endl; + // Get the correct influence table values + const auto& aqutabTable = eclState.getTableManager().getAqutabTables().getTable(aquctParams.at(aquctRecordIdx).inftableID); + const auto& aqutab_tdColumn = aqutabTable.getColumn(0); + const auto& aqutab_piColumn = aqutabTable.getColumn(1); + aquctParams.at(aquctRecordIdx).td = aqutab_tdColumn.vectorCopy(); + aquctParams.at(aquctRecordIdx).pi = aqutab_piColumn.vectorCopy(); + + // We determine the cell perforation here. + int cellID = 10 + aquctRecordIdx; + + aquctParams.at(aquctRecordIdx).cell_id = cellID; + + // We do not have mu_w as it has to be calculated from pvttable + aquifers.push_back(Aquifer_object( aquctParams.at(aquctRecordIdx) )); + } + + // I want to deliberately add another aquifer + aquifers.push_back( Aquifer_object(99) ); + + // aquifers_ = aquifers; + + return aquifers; + } + } // namespace Opm \ No newline at end of file From 89255da464cefd4aab83e706a44f2b5741e553d5 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Fri, 15 Dec 2017 13:24:15 +0100 Subject: [PATCH 20/31] Working aquifer object with equations implemented. Still need to do initialization solve. --- opm/autodiff/BlackoilAquiferModel.hpp | 1 - opm/autodiff/BlackoilAquiferModel_impl.hpp | 55 ++-------------------- 2 files changed, 5 insertions(+), 51 deletions(-) diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index 03771dced..f8b6e5c9f 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -152,7 +152,6 @@ namespace Opm { void calculateExplicitQuantities(); - // The number of components in the model. int numComponents() const; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index b9a861e44..34d02203e 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -294,60 +294,15 @@ namespace Opm { const auto& deck = ebosSimulator.gridManager().deck(); const auto& eclState = ebosSimulator.gridManager().eclState(); - if (!deck.hasKeyword("AQUCT")){ - std::cout << "Nothing is shown! Where is AQUCT!????" << std::endl; - } + // Get all the carter tracy aquifer properties data and put it in aquifers vector + AquiferCT aquiferct = AquiferCT(eclState,deck); - const auto& aquctKeyword = deck.getKeyword("AQUCT"); - std::vector aquctParams; - // Resize the parameter vector container based on row entries in aquct - // We do the same for aquifers too because number of aquifers is assumed to be for each entry in aquct - aquctParams.resize(aquctKeyword.size()); - // aquifers.resize(aquctKeyword.size()); + std::vector aquifersData = aquiferct.getAquifers(); - const int tableID = 0; - - std::cout << "Parsing AQUCT stuff" << std::endl; - for (size_t aquctRecordIdx = 0; aquctRecordIdx < aquctKeyword.size(); ++ aquctRecordIdx) + for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) { - const auto& aquctRecord = aquctKeyword.getRecord(aquctRecordIdx); - - aquctParams.at(aquctRecordIdx).aquiferID = aquctRecord.getItem("AQUIFER_ID").template get(0); - aquctParams.at(aquctRecordIdx).h = aquctRecord.getItem("THICKNESS_AQ").template get(0); - aquctParams.at(aquctRecordIdx).phi_aq = aquctRecord.getItem("PORO_AQ").template get(0); - aquctParams.at(aquctRecordIdx).d0 = aquctRecord.getItem("DAT_DEPTH").getSIDouble(0); - aquctParams.at(aquctRecordIdx).C_t = aquctRecord.getItem("C_T").template get(0); - aquctParams.at(aquctRecordIdx).r_o = aquctRecord.getItem("RAD").getSIDouble(0); - aquctParams.at(aquctRecordIdx).k_a = aquctRecord.getItem("PERM_AQ").getSIDouble(0); - aquctParams.at(aquctRecordIdx).theta = aquctRecord.getItem("INFLUENCE_ANGLE").template get(0); - aquctParams.at(aquctRecordIdx).c1 = 0.008527; // We are using SI - aquctParams.at(aquctRecordIdx).c2 = 6.283; - aquctParams.at(aquctRecordIdx).inftableID = aquctRecord.getItem("TABLE_NUM_INFLUENCE_FN").template get(0) - 1; - aquctParams.at(aquctRecordIdx).pvttableID = aquctRecord.getItem("TABLE_NUM_WATER_PRESS").template get(0) - 1; - - std::cout << aquctParams.at(aquctRecordIdx).inftableID << std::endl; - // Get the correct influence table values - const auto& aqutabTable = eclState.getTableManager().getAqutabTables().getTable(aquctParams.at(aquctRecordIdx).inftableID); - const auto& aqutab_tdColumn = aqutabTable.getColumn(0); - const auto& aqutab_piColumn = aqutabTable.getColumn(1); - aquctParams.at(aquctRecordIdx).td = aqutab_tdColumn.vectorCopy(); - aquctParams.at(aquctRecordIdx).pi = aqutab_piColumn.vectorCopy(); - - // We determine the cell perforation here. - int cellID = 10 + aquctRecordIdx; - - aquctParams.at(aquctRecordIdx).cell_id = cellID; - - // We do not have mu_w as it has to be calculated from pvttable - aquifers.push_back(Aquifer_object( aquctParams.at(aquctRecordIdx) )); + aquifers.push_back( AquiferCarterTracy (*aquiferData, numComponents(), gravity_ ) ); } - - // I want to deliberately add another aquifer - aquifers.push_back( Aquifer_object(99) ); - - // aquifers_ = aquifers; - - return aquifers; } } // namespace Opm \ No newline at end of file From e2513038c4cfc921c3768047c0b704585200ad08 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 4 Jan 2018 15:50:04 +0100 Subject: [PATCH 21/31] Implemented initialization procedure. --- opm/autodiff/AQUCT_params.hpp | 31 -- opm/autodiff/AquiferCarterTracy.hpp | 349 ++++++++++-------- opm/autodiff/BlackoilAquiferModel.hpp | 25 +- opm/autodiff/BlackoilAquiferModel_impl.hpp | 118 +----- opm/autodiff/BlackoilModelEbos.hpp | 54 +-- .../SimulatorFullyImplicitBlackoilEbos.hpp | 2 +- 6 files changed, 243 insertions(+), 336 deletions(-) delete mode 100644 opm/autodiff/AQUCT_params.hpp diff --git a/opm/autodiff/AQUCT_params.hpp b/opm/autodiff/AQUCT_params.hpp deleted file mode 100644 index da6c03ff5..000000000 --- a/opm/autodiff/AQUCT_params.hpp +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef OPM_AQUCT_HEADER_INCLUDED -#define OPM_AQUCT_HEADER_INCLUDED - - - struct AQUCT_params{ - - // Aquifer ID - int aquiferID; - // Table IDs - int inftableID, pvttableID; - // Perforation cell id - int cell_id; - // Variables constants - double phi_aq , //aquifer porosity - d0, //aquifer datum depth - C_t , //total compressibility - r_o , //aquifer inner radius - k_a , //aquifer permeability - c1, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) - h , //aquifer thickness - theta , //angle subtended by the aquifer boundary - c2 ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). - std::vector td, pi; - - - }; - - - - -#endif \ No newline at end of file diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index b6f9cfa17..4ba7597c4 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -30,10 +31,12 @@ #include #include +#include #include #include #include +#include #include #include @@ -65,44 +68,30 @@ namespace Opm typedef double Scalar; typedef DenseAd::Evaluation Eval; + typedef Opm::BlackOilFluidState FluidState; + typedef typename FluidSystem::WaterPvt WaterPvt; typedef Ewoms::BlackOilPolymerModule PolymerModule; + + + static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; - AquiferCarterTracy(const std::vector& cell_id) - : phi_aq_ (1.0), // - C_t_ (1.0), // - r_o_ (1.0), // - k_a_ (1.0), // - c1_ (1.0), - h_ (1.0), // - theta_ (1.0), // - c2_ (1.0), // - d0_ (1.0), - cell_idx_ (cell_id) - { - mu_w_ = 1e-3; - aqutab_td_.push_back(1.0); - aqutab_pi_.push_back(1.0); - aquiferID_ = 1; - inftableID_ = 1; - pvttableID_ = 1; - init_quantities(); - } - - explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const AquiferCT::AQUANCON_data& aquanconParams, - const int numComponents, const Scalar gravity ) + explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, + const int numComponents, const Scalar gravity, const Simulator& ebosSimulator ) : phi_aq_ (params.phi_aq), // C_t_ (params.C_t), // r_o_ (params.r_o), // k_a_ (params.k_a), // c1_ (params.c1), h_ (params.h), // + p0_defaulted_ (params.p0_defaulted), + pa0_ (params.p0), theta_ (params.theta), // c2_ (params.c2), // d0_ (params.d0), @@ -111,12 +100,11 @@ namespace Opm aquiferID_ (params.aquiferID), inftableID_ (params.inftableID), pvttableID_ (params.pvttableID), - cell_idx_ (params.cell_id), num_components_ (numComponents), - gravity_ (gravity) + gravity_ (gravity), + ebos_simulator_ (ebosSimulator) { - mu_w_ = 1e-3; - init_quantities(aquanconParams); + init_quantities(connection); } inline const PhaseUsage& @@ -127,79 +115,35 @@ namespace Opm return *phase_usage_; } - inline int - flowPhaseToEbosCompIdx( const int phaseIdx ) const - { - const auto& pu = phaseUsage(); - if (pu.phase_pos[Water] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - if (pu.phase_pos[Oil] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (pu.phase_pos[Gas] == phaseIdx) - return BlackoilIndices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - - // for other phases return the index - return phaseIdx; - } - - inline int - flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const auto& pu = phaseUsage(); - if (pu.phase_pos[Water] == phaseIdx) { - return FluidSystem::waterPhaseIdx; - } - if (pu.phase_pos[Oil] == phaseIdx) { - return FluidSystem::oilPhaseIdx; - } - if (pu.phase_pos[Gas] == phaseIdx) { - return FluidSystem::gasPhaseIdx; - } - - assert(phaseIdx < 3); - // for other phases return the index - return phaseIdx; - } - - inline void calculateExplicitQuantities(const Simulator& ebosSimulator) - { - std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; - } inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) { - std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ << std::endl; - // resAqui_ = 0.0; dt_ = timer.currentStepLength(); auto& ebosJac = ebosSimulator.model().linearizer().matrix(); auto& ebosResid = ebosSimulator.model().linearizer().residual(); - // TODO: it probably can be static member for StandardWell - const double volume = 0.002831684659200; // 0.1 cu ft; auto cellID = cell_idx_.begin(); + size_t idx; for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) { Eval qinflow = 0.0; // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to // IntensiveQuantities of that particular cell_id - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + const IntensiveQuantities intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); // This is the pressure at td + dt get_current_Pressure_cell(pressure_current_,idx,intQuants); get_current_density_cell(rhow_,idx,intQuants); calculate_inflow_rate(idx, timer); - qinflow = Qai_[idx]; - ebosResid[*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)] -= qinflow.value(); + qinflow = Qai_.at(idx); + ebosResid[*cellID][(FluidSystem::waterCompIdx)] -= qinflow.value(); for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[*cellID][*cellID][flowPhaseToEbosCompIdx(FluidSystem::waterPhaseIdx)][pvIdx] -= qinflow.derivative(pvIdx); + ebosJac[*cellID][*cellID][(FluidSystem::waterCompIdx)][pvIdx] -= qinflow.derivative(pvIdx); } - std::cout << "In CarterTracy: I am aquifer #" << aquiferID_ - // << " -> P_wat[t+dt] = " << pressure_current_[idx] << std::endl - << " Qai[t+dt] = " << Qai_[idx] << std::endl; } } @@ -214,50 +158,19 @@ namespace Opm } } - inline void after_time_step() + inline void after_time_step(const SimulatorTimerInterface& timer) { for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) { - W_flux_ += (*Qai); - } - std::cout << "Aquifer # " << aquiferID_ << ": My cumulative flux = " << W_flux_ << std::endl; - } - - /* Made into public for testing only!!!!!!. Must be protected */ - inline const Scalar time_constant() const - { - Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); - return Tc; - } - - /* Made into public for testing only!!!!!!. Must be protected */ - inline const Scalar aquifer_influx_constant() const - { - Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; - return beta; - } - - // This is another hack to get the face area only for SPE1. - // Ideally it should be a map which given a cell_id, it returns the area fraction - inline const double area_fraction(const int i) - { - return 1000.0*20.0*0.092903/(1000.0*1000.0*0.092903*2 + 1000.0*20.0*0.092903*4); - } - - inline void print_private_members() const - { - std::cout << "Aquifer CT #" << aquiferID_ << std::endl; - auto ita = aqutab_td_.cbegin(); - auto f_lambda = [&ita] (double i) {std::cout << *ita++ << " " << i << std::endl;}; - std::for_each( aqutab_pi_.cbegin(), aqutab_pi_.cend(), f_lambda ); - - for (auto i = coeff_.begin(); i != coeff_.end(); ++i ) - { - std::cout << "Coeff = " << *i << std::endl; + W_flux_ += (*Qai)*timer.currentStepLength(); } } - /* Made into public for testing only!!!!!!. Must be protected */ + inline const double area_fraction(const size_t i) + { + return alphai_.at(i); + } + inline const std::vector cell_id() const { return cell_idx_; @@ -269,9 +182,9 @@ namespace Opm } - - protected: + private: const PhaseUsage* phase_usage_; + const Simulator& ebos_simulator_; // Aquifer ID, and other IDs @@ -279,14 +192,17 @@ namespace Opm int num_components_; // Grid variables - std::vector cell_idx_; + + std::vector cell_idx_; + std::vector faceArea_connected_; // Quantities at each grid id std::vector cell_depth_; - std::vector pressure_previous_; - std::vector pressure_current_; - std::vector Qai_; - std::vector rhow_; + std::vector pressure_previous_; + std::vector pressure_current_; + std::vector Qai_; + std::vector rhow_; + std::vector alphai_; // Variables constants Scalar mu_w_ , //water viscosity @@ -304,8 +220,9 @@ namespace Opm std::vector aqutab_td_, aqutab_pi_; // Cumulative flux - Scalar W_flux_, dt_, pa0_, gravity_; - + Scalar dt_, pa0_, gravity_; + bool p0_defaulted_; + Eval W_flux_; // Also return the polynomial fit std::vector coeff_; @@ -316,7 +233,7 @@ namespace Opm inline void polynomial_fit( const std::vector &X, const std::vector &y, std::vector &coeff, int order, bool bias) const { - int colNum = (bias)? order + 1 : order; + size_t colNum = (bias)? order + 1 : order; Eigen::MatrixXd A(X.size(), colNum); Eigen::VectorXd y_mapped = Eigen::VectorXd::Map(&y.front(), y.size()); Eigen::VectorXd result; @@ -337,62 +254,202 @@ namespace Opm coeff[i] = result[i]; } - inline void init_quantities(const AquiferCT::AQUANCON_data& aquanconParams) + inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { + // http://kluge.in-chemnitz.de/opensource/spline/ + } + + inline void init_quantities(const Aquancon::AquanconOutput& connection) + { + // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 W_flux_ = 0.; - // pa0_ is the initial aquifer water pressure. Must be calculated from equilibrium if left default, - // or we get the information from the deck Hacked to make it at 45e6 Pa - pa0_ = 45e6; + + // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function + initialize_connections(connection); + + calculate_aquifer_condition(); pressure_previous_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.); - // We hack the cell depth values for now. We can actually get it from elementcontext pos - cell_depth_.resize(cell_idx_.size(), d0_); - rhow_.resize(cell_idx_.size(), 998.0); - Qai_.resize(cell_idx_.size(), 0.); + Qai_.resize(cell_idx_.size(), 0.0); - polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 2, true); + polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 1, true); } - inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - pressure_water[idx] = fs.pressure(FluidSystem::waterPhaseIdx).value(); + pressure_water.at(idx) = fs.pressure(FluidSystem::waterPhaseIdx); } - inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) + inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - rho_water[idx] = fs.density(FluidSystem::waterPhaseIdx).value(); + rho_water.at(idx) = fs.density(FluidSystem::waterPhaseIdx); } inline Scalar dpai(int idx) { - Scalar dp = pa0_ - rhow_[idx]*gravity_*(cell_depth_[idx] - d0_) - pressure_previous_[idx]; + Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_) - pressure_previous_.at(idx).value(); return dp; } + // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) { - // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription Scalar beta = aquifer_influx_constant(); Scalar Tc = time_constant(); Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; Scalar td = timer.simulationTimeElapsed() / Tc; - Scalar PItdprime = coeff_[1] + 2.0*coeff_[2]*(td_plus_dt); - Scalar PItd = coeff_[0] + coeff_[1]*td_plus_dt + coeff_[2]*td_plus_dt*td_plus_dt; - a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_ * PItdprime) ) / ( PItd - td*PItdprime ); - b = beta / Tc / ( PItd - td*PItdprime); + + Scalar PItdprime = coeff_.at(1); + Scalar PItd = coeff_.at(0) + coeff_.at(1)*td_plus_dt; + a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); + b = beta / (Tc * ( PItd - td*PItdprime)); } + // This function implements Eq 5.7 of the EclipseTechnicalDescription inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer) { Scalar a, b; calculate_a_b_constants(a,b,idx,timer); - // This function implements Eq 5.7 of the EclipseTechnicalDescription - Qai_[idx] = area_fraction(idx)*( a - b * ( pressure_current_[idx] - pressure_previous_[idx] ) ); + Qai_.at(idx) = area_fraction(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx).value() ) ); } - + + inline const Scalar time_constant() const + { + Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); + return Tc; + } + + inline const Scalar aquifer_influx_constant() const + { + Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; + return beta; + } + + // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer + inline void initialize_connections(const Aquancon::AquanconOutput& connection) + { + const auto& eclState = ebos_simulator_.vanguard().eclState(); + const auto& ugrid = ebos_simulator_.vanguard().grid(); + const auto& grid = eclState.getInputGrid(); + + cell_idx_ = connection.global_index; + + assert( cell_idx_ == connection.global_index); + assert( (cell_idx_.size() == connection.influx_coeff.size()) ); + assert( (connection.influx_coeff.size() == connection.influx_multiplier.size()) ); + assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) ); + + // We hack the cell depth values for now. We can actually get it from elementcontext pos + cell_depth_.resize(cell_idx_.size(), d0_); + alphai_.resize(cell_idx_.size(), 1.0); + faceArea_connected_.resize(cell_idx_.size(),0.0); + + auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); + auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); + + // Translate the C face tag into the enum used by opm-parser's TransMult class + Opm::FaceDir::DirEnum faceDirection; + + // denom_face_areas is the sum of the areas connected to an aquifer + Scalar denom_face_areas = 0.; + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + { + auto cellFacesRange = cell2Faces[cell_idx_.at(idx)]; + + for(auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) + { + // The index of the face in the compressed grid + const int faceIdx = *cellFaceIter; + + // the logically-Cartesian direction of the face + const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); + + + if (faceTag == 0) // left + faceDirection = Opm::FaceDir::XMinus; + else if (faceTag == 1) // right + faceDirection = Opm::FaceDir::XPlus; + else if (faceTag == 2) // back + faceDirection = Opm::FaceDir::YMinus; + else if (faceTag == 3) // front + faceDirection = Opm::FaceDir::YPlus; + else if (faceTag == 4) // bottom + faceDirection = Opm::FaceDir::ZMinus; + else if (faceTag == 5) // top + faceDirection = Opm::FaceDir::ZPlus; + + if (faceDirection == connection.reservoir_face_dir.at(idx)) + { + faceArea_connected_.at(idx) = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); + denom_face_areas += faceArea_connected_.at(idx); + } + } + auto cellCenter = grid.getCellCenter(cell_idx_.at(idx)); + cell_depth_.at(idx) = cellCenter[2]; + } + + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + { + alphai_.at(idx) = faceArea_connected_.at(idx)/denom_face_areas; + } + } + + inline void calculate_aquifer_condition() + { + + int pvttableIdx = pvttableID_ - 1; + + rhow_.resize(cell_idx_.size(),0.); + + if (p0_defaulted_) + { + pa0_ = calculate_reservoir_equilibrium(rhow_); + } + + // Initialize a FluidState object first + FluidState fs_aquifer; + // We use the temperature of the first cell connected to the aquifer + // Here we copy the fluidstate of the first cell, so we do not accidentally mess up the reservoir fs + fs_aquifer.assign( ebos_simulator_.model().cachedIntensiveQuantities(cell_idx_.at(0), /*timeIdx=*/ 0)->fluidState() ); + Eval temperature_aq, pa0_mean; + temperature_aq = fs_aquifer.temperature(0); + pa0_mean = pa0_; + + // rho_mean = FluidSystem::referenceDensity(waterPhaseIdx, pvttableIdx) + // *FluidSystem::waterPvt().inverseFormationVolumeFactor(pvttableIdx, temperature_aq, pa0_mean); + + Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); + std::cout << "Pa0 = " << pa0_mean << ", viscosity = " << mu_w_aquifer.value() << std::endl; + + mu_w_ = mu_w_aquifer.value(); + + } + + // This function is for calculating the aquifer properties from equilibrium state with the reservoir + inline Scalar calculate_reservoir_equilibrium(std::vector& rho_water_reservoir) + { + // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices + std::vector water_pressure_reservoir, pw_aquifer; + + for (size_t idx = 0; idx < cell_idx_.size(); ++idx) + { + size_t cellIDx = cell_idx_.at(idx); + const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0)); + const auto& fs = intQuants.fluidState(); + + water_pressure_reservoir.push_back( fs.pressure(FluidSystem::waterPhaseIdx).value() ); + rho_water_reservoir.at(idx) = fs.density(FluidSystem::waterPhaseIdx); + pw_aquifer.push_back( (water_pressure_reservoir.at(idx) - rho_water_reservoir.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_))*area_fraction(idx) ); + } + + // We take the average of the calculated equilibrium pressures. + Scalar aquifer_pres_avg = std::accumulate(pw_aquifer.begin(), pw_aquifer.end(), 0.)/pw_aquifer.size(); + return aquifer_pres_avg; + } + }; // class AquiferCarterTracy diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index f8b6e5c9f..706756ecb 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -3,11 +3,6 @@ Copyright 2017 TNO - Heat Transfer & Fluid Dynamics, Modelling & Optimization of the Subsurface Copyright 2017 Statoil ASA. - Copyright 2016 SINTEF ICT, Applied Mathematics. - Copyright 2016 - 2017 Statoil ASA. - Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services - Copyright 2016 - 2017 IRIS AS - This file is part of the Open Porous Media project (OPM). @@ -40,6 +35,7 @@ #include #include +#include #include @@ -49,7 +45,6 @@ #include #include #include - #include #include @@ -86,8 +81,6 @@ namespace Opm { static const int numEq = BlackoilIndices::numEq; static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; - typedef Ewoms::BlackOilPolymerModule PolymerModule; - typedef AquiferCarterTracy Aquifer_object; BlackoilAquiferModel(Simulator& ebosSimulator, @@ -102,7 +95,7 @@ namespace Opm { // called at the beginning of a time step void beginTimeStep(); // called at the end of a time step - void timeStepSucceeded(); + void timeStepSucceeded(const SimulatorTimerInterface& timer); // called at the beginning of a report step void beginReportStep(const int time_step); @@ -117,7 +110,7 @@ namespace Opm { return ebosSimulator_; } - /// Hack function to get what I need from parser + // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy void init(const Simulator& ebosSimulator, std::vector& aquifers); protected: @@ -141,17 +134,9 @@ namespace Opm { SimulatorReport last_report_; - const Schedule& schedule() const - { return ebosSimulator_.gridManager().schedule(); } - - void updatePrimaryVariables(); - - void initPrimaryVariablesEvaluation() const; void updateConnectionIntensiveQuantities() const; - void calculateExplicitQuantities(); - // The number of components in the model. int numComponents() const; @@ -159,8 +144,6 @@ namespace Opm { int numPhases() const; - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; - void assembleAquiferEq(const SimulatorTimerInterface& timer); SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer); @@ -177,4 +160,4 @@ namespace Opm { } // namespace Opm #include "BlackoilAquiferModel_impl.hpp" -#endif +#endif \ No newline at end of file diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index 34d02203e..c75b3f53b 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -12,7 +12,7 @@ namespace Opm { , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) { - const auto& eclState = ebosSimulator_.gridManager().eclState(); + const auto& eclState = ebosSimulator_.vanguard().eclState(); phase_usage_ = phaseUsageFromDeck(eclState); active_.resize(phase_usage_.MaxNumPhases, false); @@ -38,17 +38,17 @@ namespace Opm { void BlackoilAquiferModel:: beginTimeStep() { - // Right now it doesn't do shit. + } // called at the end of a time step template void - BlackoilAquiferModel:: timeStepSucceeded() + BlackoilAquiferModel:: timeStepSucceeded(const SimulatorTimerInterface& timer) { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->after_time_step(); + aquifer->after_time_step(timer); } } @@ -57,7 +57,7 @@ namespace Opm { void BlackoilAquiferModel:: beginReportStep(const int time_step) { - // Right now it doesn't do shit. + } // called at the end of a report step @@ -65,14 +65,7 @@ namespace Opm { void BlackoilAquiferModel:: endReportStep() { - // Right now it just spits out the constants for each aquifers - // We are using the simple integer indexing for the aquifers - for (int i = 0; i < numAquifers(); ++i) - { - std::cout << "Aquifer[" << i << "]" - << " : Tc = " << aquifers()[i].time_constant() - << ", beta = " << aquifers()[i].aquifer_influx_constant() << std::endl; - } + } // Get the last report step @@ -80,9 +73,6 @@ namespace Opm { const SimulatorReport& BlackoilAquiferModel:: lastReport() const { - for (auto i = aquifers_.begin(); i != aquifers_.end(); ++i){ - (*i).print_private_members(); - } return last_report_; } @@ -93,7 +83,6 @@ namespace Opm { const int iterationIdx ) { last_report_ = SimulatorReport(); - // We need to update the reservoir pressures connected to the aquifer updateConnectionIntensiveQuantities(); @@ -103,35 +92,14 @@ namespace Opm { prepareTimeStep(timer); } - if (iterationIdx == 0) { - calculateExplicitQuantities(); - } - if (param_.solve_aquifereq_initially_ && iterationIdx == 0) { // solve the aquifer equations as a pre-processing step last_report_ = solveAquiferEq(timer); } - assembleAquiferEq(timer); - last_report_.converged = true; } - // Protected function: Update the primary variables - template - void - BlackoilAquiferModel:: updatePrimaryVariables() - { - // Right now it doesn't do shit. - } - - // Protected function: Init the primary variables - template - void - BlackoilAquiferModel:: initPrimaryVariablesEvaluation() const - { - // Right now it doesn't do shit. - } template void @@ -149,17 +117,6 @@ namespace Opm { } } - template - void - BlackoilAquiferModel:: calculateExplicitQuantities() - { - // for (auto aqui = aquifers_.begin(); aqui!= aquifers_.end(); ++aqui) - // { - // std::cout << "calculateExplicitQuantities: Aquifer id = " << aqui->aquiferID() << std::endl; - // aqui->calculateExplicitQuantities(ebosSimulator_); - // } - } - template SimulatorReport @@ -200,30 +157,10 @@ namespace Opm { int BlackoilAquiferModel:: numPhases() const { - // Not implemented yet!!!!!!!!!!!! const auto& pu = phase_usage_; return pu.num_phases; } - - // Protected function: returns the phase index in ebos - template - int - BlackoilAquiferModel:: flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const auto& pu = phase_usage_; - if (active_[Water] && pu.phase_pos[Water] == phaseIdx) - return FluidSystem::waterPhaseIdx; - if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) - return FluidSystem::oilPhaseIdx; - if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) - return FluidSystem::gasPhaseIdx; - - assert(phaseIdx < 3); - // for other phases return the index - return phaseIdx; - } - // Protected function which calls the individual aquifer models template void @@ -231,7 +168,6 @@ namespace Opm { { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - std::cout << "assembleAquiferEq: Aquifer id = " << aquifer->aquiferID() << std::endl; aquifer->assembleAquiferEq(ebosSimulator_, timer); } } @@ -261,47 +197,27 @@ namespace Opm { // Initialize the aquifers in the deck template void - BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers)//, std::vector< AquiferCarterTracy >& aquifers) + BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers) { - const auto& deck = ebosSimulator.gridManager().deck(); - const auto& eclState = ebosSimulator.gridManager().eclState(); + updateConnectionIntensiveQuantities(); + const auto& deck = ebosSimulator.vanguard().deck(); + const auto& eclState = ebosSimulator.vanguard().eclState(); // Get all the carter tracy aquifer properties data and put it in aquifers vector AquiferCT aquiferct = AquiferCT(eclState,deck); + Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); std::vector aquifersData = aquiferct.getAquifers(); - std::vector aquanconData = aquiferct.getAquancon(); + std::vector aquifer_connection = aquifer_connect.getAquOutput(); - // for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) - // { - - // } + assert( aquifersData.size() == aquifer_connect.size() ); - auto ita = aquifersData.cbegin(); - auto f_lambda = [&] (AquiferCT::AQUANCON_data i) { - aquifers.push_back( AquiferCarterTracy (*ita++, i, numComponents(), gravity_ ) ); - }; - std::for_each( aquanconData.cbegin(), aquanconData.cend(), f_lambda ); - } - // Begin the hack to initialize the aquifers in the deck - template - std::vector< AquiferCarterTracy > - BlackoilAquiferModel:: hack_init(const Simulator& ebosSimulator)//, std::vector< AquiferCarterTracy >& aquifers) - { - std::vector< AquiferCarterTracy > aquifers; - /** Begin hack!!!!! */ - const auto& deck = ebosSimulator.gridManager().deck(); - const auto& eclState = ebosSimulator.gridManager().eclState(); - - // Get all the carter tracy aquifer properties data and put it in aquifers vector - AquiferCT aquiferct = AquiferCT(eclState,deck); - - std::vector aquifersData = aquiferct.getAquifers(); - - for (auto aquiferData = aquifersData.begin(); aquiferData != aquifersData.end(); ++aquiferData) + for (int i = 0; i < aquifersData.size(); ++i) { - aquifers.push_back( AquiferCarterTracy (*aquiferData, numComponents(), gravity_ ) ); + aquifers.push_back( + AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), numComponents(), gravity_, ebosSimulator_) + ); } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index c55f337bd..2408bb3cd 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -29,10 +29,10 @@ #include #include -#include #include #include #include +#include #include #include @@ -350,12 +350,12 @@ namespace Opm { const ReservoirState& reservoir_state, WellState& well_state) { - DUNE_UNUSED_PARAMETER(timer); + // DUNE_UNUSED_PARAMETER(timer); DUNE_UNUSED_PARAMETER(reservoir_state); DUNE_UNUSED_PARAMETER(well_state); wellModel().timeStepSucceeded(); - aquiferModel().timeStepSucceeded(); + aquiferModel().timeStepSucceeded(timer); ebosSimulator_.problem().endTimeStep(); } @@ -429,13 +429,13 @@ namespace Opm { if (elem.partitionType() != Dune::InteriorEntity) continue; - unsigned globalElemIdx = elemMapper.index(elem); + unsigned globalElemIdx = elemMapper.index(elem); const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; Scalar pressureNew; - pressureNew = priVarsNew[Indices::pressureSwitchIdx]; + pressureNew = priVarsNew[Indices::pressureSwitchIdx]; - Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; + Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationNew = 1.0; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx]; @@ -478,7 +478,7 @@ namespace Opm { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; - resultDelta += tmp*tmp; + resultDelta += tmp*tmp; resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; } } @@ -486,9 +486,9 @@ namespace Opm { resultDelta = gridView.comm().sum(resultDelta); resultDenom = gridView.comm().sum(resultDenom); - if (resultDenom > 0.0) - return resultDelta/resultDenom; - return 0.0; + if (resultDenom > 0.0) + return resultDelta/resultDenom; + return 0.0; } @@ -599,7 +599,6 @@ namespace Opm { virtual void apply( const X& x, Y& y ) const { A_.mv( x, y ); - // add well model modification to y wellMod_.apply(x, y ); @@ -613,7 +612,6 @@ namespace Opm { virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const { A_.usmv(alpha,x,y); - // add scaled well model modification to y wellMod_.applyScaleAdd( alpha, x, y ); @@ -1129,29 +1127,14 @@ namespace Opm { const BlackoilAquiferModel& aquiferModel() const { return aquifer_model_; } - int flowPhaseToEbosCompIdx( const int phaseIdx ) const + int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const { - const auto& pu = phaseUsage_; - if (active_[Water] && pu.phase_pos[Water] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) - return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - - // for other phases return the index - return phaseIdx; - } - - int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const - { - const auto& pu = phaseUsage_; - if (active_[Water] && pu.phase_pos[Water] == phaseIdx) - return FluidSystem::waterPhaseIdx; - if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx) - return FluidSystem::oilPhaseIdx; - if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx) - return FluidSystem::gasPhaseIdx; + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx) + return Water; + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx) + return Oil; + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx) + return Gas; assert(phaseIdx < 3); // for other phases return the index @@ -1170,7 +1153,6 @@ namespace Opm { private: - double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } double drMaxRel() const { return param_.dr_max_rel_; } @@ -1181,4 +1163,4 @@ namespace Opm { }; } // namespace Opm -#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED +#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED \ No newline at end of file diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index 5e5988b22..cb0c66399 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -190,8 +190,8 @@ public: auto auxMod = std::make_shared >(schedule(), grid()); ebosSimulator_.model().addAuxiliaryModule(auxMod); } + AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_); - // aquifer_model.hack_init(ebosSimulator_); // Main simulation loop. while (!timer.done()) { From d8f864645223f97bb5bc579fbb3da54014f988cd Mon Sep 17 00:00:00 2001 From: kel85uk Date: Wed, 28 Mar 2018 20:43:35 +0200 Subject: [PATCH 22/31] Removes Eigen from interpolation of influence tables --- opm/autodiff/AquiferCarterTracy.hpp | 48 +++++------------------------ 1 file changed, 8 insertions(+), 40 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 4ba7597c4..3af93236d 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -21,11 +21,11 @@ #ifndef OPM_AQUIFERCT_HEADER_INCLUDED #define OPM_AQUIFERCT_HEADER_INCLUDED -#include #include #include #include #include +#include #include @@ -223,40 +223,14 @@ namespace Opm Scalar dt_, pa0_, gravity_; bool p0_defaulted_; Eval W_flux_; - // Also return the polynomial fit - std::vector coeff_; + - // We fit the tabular data using a polynomial fit - // Modified from Copyright (C) 2014 Clifford Wolf - // http://svn.clifford.at/handicraft/2014/polyfit/polyfit.cc - inline void polynomial_fit( const std::vector &X, const std::vector &y, - std::vector &coeff, int order, bool bias) const - { - size_t colNum = (bias)? order + 1 : order; - Eigen::MatrixXd A(X.size(), colNum); - Eigen::VectorXd y_mapped = Eigen::VectorXd::Map(&y.front(), y.size()); - Eigen::VectorXd result; - - assert(X.size() == y.size()); - assert(X.size() >= colNum); - - // create matrix - for (size_t i = 0; i < X.size(); i++) - for (size_t j = 0; j < colNum; j++) - A(i, j) = (bias)? pow(X.at(i), j) : pow(X.at(i), j+1); - - // solve for linear least squares fit - result = A.householderQr().solve(y_mapped); - - coeff.resize(colNum); - for (size_t i = 0; i < colNum; i++) - coeff[i] = result[i]; - } - inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { - // http://kluge.in-chemnitz.de/opensource/spline/ + // We use the opm-common numeric linear interpolator + pitd = Opm::linearInterpolation(aqutab_td_, aqutab_pi_, td); + pitd_prime = Opm::linearInterpolationDerivative(aqutab_td_, aqutab_pi_, td); } inline void init_quantities(const Aquancon::AquanconOutput& connection) @@ -272,8 +246,6 @@ namespace Opm pressure_previous_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.); Qai_.resize(cell_idx_.size(), 0.0); - - polynomial_fit(aqutab_td_, aqutab_pi_, coeff_, 1, true); } inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) @@ -301,9 +273,9 @@ namespace Opm Scalar Tc = time_constant(); Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; Scalar td = timer.simulationTimeElapsed() / Tc; - - Scalar PItdprime = coeff_.at(1); - Scalar PItd = coeff_.at(0) + coeff_.at(1)*td_plus_dt; + Scalar PItdprime = 0.; + Scalar PItd = 0.; + get_influence_table_values(PItd, PItdprime, td_plus_dt); a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); b = beta / (Tc * ( PItd - td*PItdprime)); } @@ -418,11 +390,7 @@ namespace Opm temperature_aq = fs_aquifer.temperature(0); pa0_mean = pa0_; - // rho_mean = FluidSystem::referenceDensity(waterPhaseIdx, pvttableIdx) - // *FluidSystem::waterPvt().inverseFormationVolumeFactor(pvttableIdx, temperature_aq, pa0_mean); - Eval mu_w_aquifer = FluidSystem::waterPvt().viscosity(pvttableIdx, temperature_aq, pa0_mean); - std::cout << "Pa0 = " << pa0_mean << ", viscosity = " << mu_w_aquifer.value() << std::endl; mu_w_ = mu_w_aquifer.value(); From 0df51a8797366560abd51fe0715ad000c2143a31 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 5 Apr 2018 15:10:10 +0200 Subject: [PATCH 23/31] Implements the correct check for aquifer face connections to make sure face is connected to the boundary --- opm/autodiff/AquiferCarterTracy.hpp | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 3af93236d..d1a655df9 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -308,6 +308,14 @@ namespace Opm const auto& grid = eclState.getInputGrid(); cell_idx_ = connection.global_index; + auto globalCellIdx = ugrid.globalCell(); + // for (auto globalCells : globalCellIdx){ + // std::cout << "global id = " << globalCells << std::endl; + // } + + // for (auto cellidx : cell_idx_){ + // std::cout << "aqucell id = " << cellidx << std::endl; + // } assert( cell_idx_ == connection.global_index); assert( (cell_idx_.size() == connection.influx_coeff.size()) ); @@ -318,10 +326,14 @@ namespace Opm cell_depth_.resize(cell_idx_.size(), d0_); alphai_.resize(cell_idx_.size(), 1.0); faceArea_connected_.resize(cell_idx_.size(),0.0); + Scalar faceArea; auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); + + // static_assert(decltype(faceCells)::dummy_error, "DUMP MY TYPE" ); + // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; @@ -355,7 +367,11 @@ namespace Opm if (faceDirection == connection.reservoir_face_dir.at(idx)) { - faceArea_connected_.at(idx) = Opm::UgGridHelpers::faceArea(ugrid, faceIdx); + // Check now if the face is outside of the reservoir, or if it adjoins an inactive cell + // Do not make the connection if the product of the two cellIdx > 0. This is because the + // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) + faceArea = (faceCells(faceIdx,0)*faceCells(faceIdx,1) > 0)? 0. : Opm::UgGridHelpers::faceArea(ugrid, faceIdx); + faceArea_connected_.at(idx) = faceArea; denom_face_areas += faceArea_connected_.at(idx); } } From 1a33e81d9cd65816f1f86cf1a80c2f7c4950ac69 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Tue, 10 Apr 2018 13:05:30 +0200 Subject: [PATCH 24/31] Cleans up the code and removes unused functions --- opm/autodiff/AquiferCarterTracy.hpp | 40 ++------- opm/autodiff/BlackoilAquiferModel.hpp | 27 +----- opm/autodiff/BlackoilAquiferModel_impl.hpp | 86 +------------------ opm/autodiff/BlackoilModelEbos.hpp | 15 ---- opm/autodiff/BlackoilModelParameters.cpp | 2 - opm/autodiff/BlackoilModelParameters.hpp | 3 - opm/autodiff/BlackoilTransportModel.hpp | 1 - .../SimulatorFullyImplicitBlackoilEbos.hpp | 5 -- 8 files changed, 7 insertions(+), 172 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index d1a655df9..2d4c35dcb 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -69,21 +69,12 @@ namespace Opm typedef DenseAd::Evaluation Eval; typedef Opm::BlackOilFluidState FluidState; - typedef typename FluidSystem::WaterPvt WaterPvt; - - typedef Ewoms::BlackOilPolymerModule PolymerModule; - static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent); - static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer); - static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx; - static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx; - - explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, - const int numComponents, const Scalar gravity, const Simulator& ebosSimulator ) + const Scalar gravity, const Simulator& ebosSimulator ) : phi_aq_ (params.phi_aq), // C_t_ (params.C_t), // r_o_ (params.r_o), // @@ -100,22 +91,12 @@ namespace Opm aquiferID_ (params.aquiferID), inftableID_ (params.inftableID), pvttableID_ (params.pvttableID), - num_components_ (numComponents), gravity_ (gravity), ebos_simulator_ (ebosSimulator) { init_quantities(connection); } - inline const PhaseUsage& - phaseUsage() const - { - assert(phase_usage_); - - return *phase_usage_; - } - - inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) { dt_ = timer.currentStepLength(); @@ -166,10 +147,6 @@ namespace Opm } } - inline const double area_fraction(const size_t i) - { - return alphai_.at(i); - } inline const std::vector cell_id() const { @@ -183,13 +160,11 @@ namespace Opm private: - const PhaseUsage* phase_usage_; const Simulator& ebos_simulator_; // Aquifer ID, and other IDs int aquiferID_, inftableID_, pvttableID_; - int num_components_; // Grid variables @@ -225,6 +200,10 @@ namespace Opm Eval W_flux_; + inline const double area_fraction(const size_t i) + { + return alphai_.at(i); + } inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { @@ -309,13 +288,6 @@ namespace Opm cell_idx_ = connection.global_index; auto globalCellIdx = ugrid.globalCell(); - // for (auto globalCells : globalCellIdx){ - // std::cout << "global id = " << globalCells << std::endl; - // } - - // for (auto cellidx : cell_idx_){ - // std::cout << "aqucell id = " << cellidx << std::endl; - // } assert( cell_idx_ == connection.global_index); assert( (cell_idx_.size() == connection.influx_coeff.size()) ); @@ -332,8 +304,6 @@ namespace Opm auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); - // static_assert(decltype(faceCells)::dummy_error, "DUMP MY TYPE" ); - // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index 706756ecb..485c3dca6 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -92,19 +92,9 @@ namespace Opm { void assemble( const SimulatorTimerInterface& timer, const int iterationIdx ); - // called at the beginning of a time step - void beginTimeStep(); // called at the end of a time step void timeStepSucceeded(const SimulatorTimerInterface& timer); - // called at the beginning of a report step - void beginReportStep(const int time_step); - - // called at the end of a report step - void endReportStep(); - - const SimulatorReport& lastReport() const; - inline const Simulator& simulator() const { return ebosSimulator_; @@ -119,11 +109,7 @@ namespace Opm { const ModelParameters param_; bool terminal_output_; - bool has_solvent_; - bool has_polymer_; - std::vector pvt_region_idx_; - PhaseUsage phase_usage_; - std::vector active_; + size_t global_nc_; // the number of the cells in the local grid size_t number_of_cells_; @@ -132,23 +118,12 @@ namespace Opm { std::vector aquifers_; - SimulatorReport last_report_; - - void updateConnectionIntensiveQuantities() const; - // The number of components in the model. - int numComponents() const; - int numAquifers() const; - int numPhases() const; - void assembleAquiferEq(const SimulatorTimerInterface& timer); - SimulatorReport solveAquiferEq(const SimulatorTimerInterface& timer); - - // some preparation work, mostly related to group control and RESV, // at the beginning of each time step (Not report step) void prepareTimeStep(const SimulatorTimerInterface& timer); diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index c75b3f53b..3087d951e 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -9,22 +9,11 @@ namespace Opm { : ebosSimulator_(ebosSimulator) , param_(param) , terminal_output_(terminal_output) - , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent)) - , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer)) { const auto& eclState = ebosSimulator_.vanguard().eclState(); - phase_usage_ = phaseUsageFromDeck(eclState); - - active_.resize(phase_usage_.MaxNumPhases, false); - for (int p = 0; p < phase_usage_.MaxNumPhases; ++p) { - active_[ p ] = phase_usage_.phase_used[ p ] != 0; - } const auto& gridView = ebosSimulator_.gridView(); - // calculate the number of elements of the compressed sequential grid. this needs - // to be done in two steps because the dune communicator expects a reference as - // argument for sum() number_of_cells_ = gridView.size(/*codim=*/0); global_nc_ = gridView.comm().sum(number_of_cells_); gravity_ = ebosSimulator_.problem().gravity()[2]; @@ -32,15 +21,6 @@ namespace Opm { } - - // called at the beginning of a time step - template - void - BlackoilAquiferModel:: beginTimeStep() - { - - } - // called at the end of a time step template void @@ -52,37 +32,12 @@ namespace Opm { } } - // called at the beginning of a report step - template - void - BlackoilAquiferModel:: beginReportStep(const int time_step) - { - - } - - // called at the end of a report step - template - void - BlackoilAquiferModel:: endReportStep() - { - - } - - // Get the last report step - template - const SimulatorReport& - BlackoilAquiferModel:: lastReport() const - { - return last_report_; - } - template void BlackoilAquiferModel:: assemble( const SimulatorTimerInterface& timer, const int iterationIdx ) { - last_report_ = SimulatorReport(); // We need to update the reservoir pressures connected to the aquifer updateConnectionIntensiveQuantities(); @@ -92,12 +47,7 @@ namespace Opm { prepareTimeStep(timer); } - if (param_.solve_aquifereq_initially_ && iterationIdx == 0) { - // solve the aquifer equations as a pre-processing step - last_report_ = solveAquiferEq(timer); - } assembleAquiferEq(timer); - last_report_.converged = true; } @@ -118,32 +68,6 @@ namespace Opm { } - template - SimulatorReport - BlackoilAquiferModel:: solveAquiferEq(const SimulatorTimerInterface& timer) - { - // We need to solve the equilibrium equation first to - // obtain the initial pressure of water in the aquifer - SimulatorReport report = SimulatorReport(); - return report; - } - - // Protected function: Return number of components in the model. - template - int - BlackoilAquiferModel:: numComponents() const - { - if (numPhases() == 2) { - return 2; - } - int numComp = FluidSystem::numComponents; - if (has_solvent_) { - numComp ++; - } - - return numComp; - } - // Protected function: Return number of aquifers in the model. template int @@ -152,14 +76,6 @@ namespace Opm { return aquifers_.size(); } - // Protected function: Return number of phases in the model. - template - int - BlackoilAquiferModel:: numPhases() const - { - const auto& pu = phase_usage_; - return pu.num_phases; - } // Protected function which calls the individual aquifer models template @@ -216,7 +132,7 @@ namespace Opm { for (int i = 0; i < aquifersData.size(); ++i) { aquifers.push_back( - AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), numComponents(), gravity_, ebosSimulator_) + AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), gravity_, ebosSimulator_) ); } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 2408bb3cd..4d91fa59d 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -209,7 +209,6 @@ namespace Opm { wasSwitched_.resize(numDof); std::fill(wasSwitched_.begin(), wasSwitched_.end(), false); - aquiferModel().beginTimeStep(); wellModel().beginTimeStep(); if (param_.update_equations_scaling_) { @@ -1127,20 +1126,6 @@ namespace Opm { const BlackoilAquiferModel& aquiferModel() const { return aquifer_model_; } - int ebosPhaseToFlowCanonicalPhaseIdx( const int phaseIdx ) const - { - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && FluidSystem::waterPhaseIdx == phaseIdx) - return Water; - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::oilPhaseIdx == phaseIdx) - return Oil; - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && FluidSystem::gasPhaseIdx == phaseIdx) - return Gas; - - assert(phaseIdx < 3); - // for other phases return the index - return phaseIdx; - } - void beginReportStep() { ebosSimulator_.problem().beginEpisode(); diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index 2bd302c2d..f16fc7685 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -63,7 +63,6 @@ namespace Opm param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day ); max_strict_iter_ = param.getDefault("max_strict_iter",8); solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_); - solve_aquifereq_initially_ = param.getDefault("solve_aquifereq_initially",solve_aquifereq_initially_); update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_); use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_); deck_file_name_ = param.template get("deck_filename"); @@ -95,7 +94,6 @@ namespace Opm max_inner_iter_ms_wells_ = 10; maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day ); solve_welleq_initially_ = true; - solve_aquifereq_initially_ = true; update_equations_scaling_ = false; use_update_stabilization_ = true; use_multisegment_well_ = false; diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index d57552161..b01cfdde9 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -77,9 +77,6 @@ namespace Opm /// Solve well equation initially bool solve_welleq_initially_; - /// Solve aquifer equation initially - bool solve_aquifereq_initially_; - /// Update scaling factors for mass balance equations bool update_equations_scaling_; diff --git a/opm/autodiff/BlackoilTransportModel.hpp b/opm/autodiff/BlackoilTransportModel.hpp index d6610f0d5..ffc0ebc92 100644 --- a/opm/autodiff/BlackoilTransportModel.hpp +++ b/opm/autodiff/BlackoilTransportModel.hpp @@ -86,7 +86,6 @@ namespace Opm { { Base::prepareStep(timer, reservoir_state, well_state); Base::param_.solve_welleq_initially_ = false; - Base::param_.solve_aquifereq_initially_ = false; SolutionState state0 = variableState(reservoir_state, well_state); asImpl().makeConstantState(state0); asImpl().computeAccum(state0, 0); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index cb0c66399..bfbbc0e6b 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -208,8 +208,6 @@ public: well_model.beginReportStep(timer.currentStepNum()); - aquifer_model.beginReportStep(timer.currentStepNum()); - auto solver = createSolver(well_model, aquifer_model); // write the inital state at the report stage @@ -274,7 +272,6 @@ public: } solver->model().endReportStep(); - aquifer_model.endReportStep(); well_model.endReportStep(); // take time that was used to solve system for this reportStep @@ -318,8 +315,6 @@ public: report.total_time = total_timer.secsSinceStart(); report.converged = true; - auto reportaquifer = aquifer_model.lastReport(); - return report; } From 2fe24e16cf2027075ba7af2cbdd4786601275231 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Wed, 11 Apr 2018 17:30:23 +0200 Subject: [PATCH 25/31] Adds the influx coefficient multipler to the face area fraction calculation. Also removes more unnecessary variables and reordered initialization list to prevent ordering warnings in gcc compiler --- opm/autodiff/AquiferCarterTracy.hpp | 45 ++++++++++++---------- opm/autodiff/BlackoilAquiferModel.hpp | 7 ---- opm/autodiff/BlackoilAquiferModel_impl.hpp | 10 ++--- 3 files changed, 29 insertions(+), 33 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 2d4c35dcb..9eabb301f 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -50,10 +50,6 @@ namespace Opm public: typedef BlackoilModelParameters ModelParameters; - static const int Water = BlackoilPhases::Aqua; - static const int Oil = BlackoilPhases::Liquid; - static const int Gas = BlackoilPhases::Vapour; - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; @@ -70,29 +66,31 @@ namespace Opm typedef DenseAd::Evaluation Eval; typedef Opm::BlackOilFluidState FluidState; + static const auto waterCompIdx = FluidSystem::waterCompIdx; + static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx; explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, const Scalar gravity, const Simulator& ebosSimulator ) - : phi_aq_ (params.phi_aq), // + : ebos_simulator_ (ebosSimulator), + aquiferID_ (params.aquiferID), + inftableID_ (params.inftableID), + pvttableID_ (params.pvttableID), + phi_aq_ (params.phi_aq), // + d0_ (params.d0), C_t_ (params.C_t), // r_o_ (params.r_o), // k_a_ (params.k_a), // c1_ (params.c1), h_ (params.h), // - p0_defaulted_ (params.p0_defaulted), - pa0_ (params.p0), theta_ (params.theta), // c2_ (params.c2), // - d0_ (params.d0), aqutab_td_ (params.td), aqutab_pi_ (params.pi), - aquiferID_ (params.aquiferID), - inftableID_ (params.inftableID), - pvttableID_ (params.pvttableID), + pa0_ (params.p0), gravity_ (gravity), - ebos_simulator_ (ebosSimulator) + p0_defaulted_ (params.p0_defaulted) { init_quantities(connection); } @@ -118,12 +116,12 @@ namespace Opm get_current_density_cell(rhow_,idx,intQuants); calculate_inflow_rate(idx, timer); qinflow = Qai_.at(idx); - ebosResid[*cellID][(FluidSystem::waterCompIdx)] -= qinflow.value(); + ebosResid[*cellID][waterCompIdx] -= qinflow.value(); for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[*cellID][*cellID][(FluidSystem::waterCompIdx)][pvIdx] -= qinflow.derivative(pvIdx); + ebosJac[*cellID][*cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx); } } } @@ -230,13 +228,13 @@ namespace Opm inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - pressure_water.at(idx) = fs.pressure(FluidSystem::waterPhaseIdx); + pressure_water.at(idx) = fs.pressure(waterPhaseIdx); } inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - rho_water.at(idx) = fs.density(FluidSystem::waterPhaseIdx); + rho_water.at(idx) = fs.density(waterPhaseIdx); } inline Scalar dpai(int idx) @@ -303,6 +301,13 @@ namespace Opm auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); + for (auto influxCoeff: connection.influx_coeff){ + std::cout << "influx_coeff = " << influxCoeff << std::endl; + } + + for (auto influxMult: connection.influx_multiplier){ + std::cout << "influx_multiplier = " << influxMult << std::endl; + } // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; @@ -342,7 +347,7 @@ namespace Opm // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell adjoining) faceArea = (faceCells(faceIdx,0)*faceCells(faceIdx,1) > 0)? 0. : Opm::UgGridHelpers::faceArea(ugrid, faceIdx); faceArea_connected_.at(idx) = faceArea; - denom_face_areas += faceArea_connected_.at(idx); + denom_face_areas += ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) ); } } auto cellCenter = grid.getCellCenter(cell_idx_.at(idx)); @@ -351,7 +356,7 @@ namespace Opm for (size_t idx = 0; idx < cell_idx_.size(); ++idx) { - alphai_.at(idx) = faceArea_connected_.at(idx)/denom_face_areas; + alphai_.at(idx) = ( connection.influx_multiplier.at(idx) * faceArea_connected_.at(idx) )/denom_face_areas; } } @@ -394,8 +399,8 @@ namespace Opm const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0)); const auto& fs = intQuants.fluidState(); - water_pressure_reservoir.push_back( fs.pressure(FluidSystem::waterPhaseIdx).value() ); - rho_water_reservoir.at(idx) = fs.density(FluidSystem::waterPhaseIdx); + water_pressure_reservoir.push_back( fs.pressure(waterPhaseIdx).value() ); + rho_water_reservoir.at(idx) = fs.density(waterPhaseIdx); pw_aquifer.push_back( (water_pressure_reservoir.at(idx) - rho_water_reservoir.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_))*area_fraction(idx) ); } diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index 485c3dca6..a2d8a5b59 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -78,9 +78,6 @@ namespace Opm { typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - static const int numEq = BlackoilIndices::numEq; - static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx; - typedef AquiferCarterTracy Aquifer_object; BlackoilAquiferModel(Simulator& ebosSimulator, @@ -110,11 +107,7 @@ namespace Opm { const ModelParameters param_; bool terminal_output_; - size_t global_nc_; - // the number of the cells in the local grid - size_t number_of_cells_; double gravity_; - std::vector depth_; std::vector aquifers_; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index 3087d951e..be1493550 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -10,12 +10,10 @@ namespace Opm { , param_(param) , terminal_output_(terminal_output) { - const auto& eclState = ebosSimulator_.vanguard().eclState(); + // const auto& gridView = ebosSimulator_.gridView(); - const auto& gridView = ebosSimulator_.gridView(); - - number_of_cells_ = gridView.size(/*codim=*/0); - global_nc_ = gridView.comm().sum(number_of_cells_); + // number_of_cells_ = gridView.size(/*codim=*/0); + // global_nc_ = gridView.comm().sum(number_of_cells_); gravity_ = ebosSimulator_.problem().gravity()[2]; init(ebosSimulator_, aquifers_); } @@ -129,7 +127,7 @@ namespace Opm { assert( aquifersData.size() == aquifer_connect.size() ); - for (int i = 0; i < aquifersData.size(); ++i) + for (size_t i = 0; i < aquifersData.size(); ++i) { aquifers.push_back( AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), gravity_, ebosSimulator_) From b7fbe66c919cf3fac81af08e16cca63a33f09cd5 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Thu, 26 Apr 2018 16:22:26 +0200 Subject: [PATCH 26/31] Cleans up the Aquifer model classes. --- opm/autodiff/AquiferCarterTracy.hpp | 254 +++++++----------- opm/autodiff/BlackoilAquiferModel.hpp | 57 +--- opm/autodiff/BlackoilAquiferModel_impl.hpp | 56 +--- opm/autodiff/BlackoilModelEbos.hpp | 14 +- .../SimulatorFullyImplicitBlackoilEbos.hpp | 4 +- 5 files changed, 125 insertions(+), 260 deletions(-) diff --git a/opm/autodiff/AquiferCarterTracy.hpp b/opm/autodiff/AquiferCarterTracy.hpp index 9eabb301f..ff9841016 100644 --- a/opm/autodiff/AquiferCarterTracy.hpp +++ b/opm/autodiff/AquiferCarterTracy.hpp @@ -24,21 +24,14 @@ #include #include #include -#include #include -#include - #include #include #include -#include -#include #include #include -#include -#include namespace Opm { @@ -48,17 +41,11 @@ namespace Opm { public: - typedef BlackoilModelParameters ModelParameters; - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; - typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; - typedef typename GridView::template Codim<0>::Entity Element; - typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; static const int numEq = BlackoilIndices::numEq; typedef double Scalar; @@ -71,73 +58,56 @@ namespace Opm - explicit AquiferCarterTracy( const AquiferCT::AQUCT_data& params, const Aquancon::AquanconOutput& connection, - const Scalar gravity, const Simulator& ebosSimulator ) - : ebos_simulator_ (ebosSimulator), - aquiferID_ (params.aquiferID), - inftableID_ (params.inftableID), - pvttableID_ (params.pvttableID), - phi_aq_ (params.phi_aq), // - d0_ (params.d0), - C_t_ (params.C_t), // - r_o_ (params.r_o), // - k_a_ (params.k_a), // - c1_ (params.c1), - h_ (params.h), // - theta_ (params.theta), // - c2_ (params.c2), // - aqutab_td_ (params.td), - aqutab_pi_ (params.pi), - pa0_ (params.p0), - gravity_ (gravity), - p0_defaulted_ (params.p0_defaulted) + AquiferCarterTracy( const AquiferCT::AQUCT_data& aquct_data, + const Aquancon::AquanconOutput& connection, + Simulator& ebosSimulator ) + : ebos_simulator_ (ebosSimulator), + aquct_data_ (aquct_data), + gravity_ (ebos_simulator_.problem().gravity()[2]) { - init_quantities(connection); + initQuantities(connection); } - inline void assembleAquiferEq(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + inline void assembleAquiferEq(const SimulatorTimerInterface& timer) { - dt_ = timer.currentStepLength(); - auto& ebosJac = ebosSimulator.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator.model().linearizer().residual(); + auto& ebosJac = ebos_simulator_.model().linearizer().matrix(); + auto& ebosResid = ebos_simulator_.model().linearizer().residual(); - - auto cellID = cell_idx_.begin(); - - size_t idx; - for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) + size_t cellID; + for ( size_t idx = 0; idx < cell_idx_.size(); ++idx ) { Eval qinflow = 0.0; + cellID = cell_idx_.at(idx); // We are dereferencing the value of IntensiveQuantities because cachedIntensiveQuantities return a const pointer to // IntensiveQuantities of that particular cell_id - const IntensiveQuantities intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + const IntensiveQuantities intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellID, /*timeIdx=*/ 0)); // This is the pressure at td + dt - get_current_Pressure_cell(pressure_current_,idx,intQuants); - get_current_density_cell(rhow_,idx,intQuants); - calculate_inflow_rate(idx, timer); + updateCellPressure(pressure_current_,idx,intQuants); + updateCellDensity(idx,intQuants); + calculateInflowRate(idx, timer); qinflow = Qai_.at(idx); - ebosResid[*cellID][waterCompIdx] -= qinflow.value(); + ebosResid[cellID][waterCompIdx] -= qinflow.value(); for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) { // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[*cellID][*cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx); + ebosJac[cellID][cellID][waterCompIdx][pvIdx] -= qinflow.derivative(pvIdx); } } } - inline void before_time_step(Simulator& ebosSimulator, const SimulatorTimerInterface& timer) + inline void beforeTimeStep(const SimulatorTimerInterface& timer) { auto cellID = cell_idx_.begin(); size_t idx; for ( idx = 0; cellID != cell_idx_.end(); ++cellID, ++idx ) { - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); - get_current_Pressure_cell(pressure_previous_ ,idx,intQuants); + const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(*cellID, /*timeIdx=*/ 0)); + updateCellPressure(pressure_previous_ ,idx,intQuants); } } - inline void after_time_step(const SimulatorTimerInterface& timer) + inline void afterTimeStep(const SimulatorTimerInterface& timer) { for (auto Qai = Qai_.begin(); Qai != Qai_.end(); ++Qai) { @@ -146,139 +116,117 @@ namespace Opm } - inline const std::vector cell_id() const - { - return cell_idx_; - } - - inline const int& aquiferID() const - { - return aquiferID_; - } - - private: - const Simulator& ebos_simulator_; - - - // Aquifer ID, and other IDs - int aquiferID_, inftableID_, pvttableID_; + Simulator& ebos_simulator_; // Grid variables - std::vector cell_idx_; std::vector faceArea_connected_; // Quantities at each grid id std::vector cell_depth_; - std::vector pressure_previous_; + std::vector pressure_previous_; std::vector pressure_current_; std::vector Qai_; std::vector rhow_; std::vector alphai_; // Variables constants - Scalar mu_w_ , //water viscosity - phi_aq_ , //aquifer porosity - d0_, // aquifer datum depth - C_t_ , //total compressibility - r_o_ , //aquifer inner radius - k_a_ , //aquifer permeability - c1_, // 0.008527 (METRIC, PVT-M); 0.006328 (FIELD); 3.6 (LAB) - h_ , //aquifer thickness - theta_ , //angle subtended by the aquifer boundary - c2_ ; //6.283 (METRIC, PVT-M); 1.1191 (FIELD); 6.283 (LAB). + const AquiferCT::AQUCT_data aquct_data_; - // Variables for influence table - std::vector aqutab_td_, aqutab_pi_; + Scalar mu_w_ , //water viscosity + beta_ , // Influx constant + Tc_ , // Time constant + pa0_ , // initial aquifer pressure + gravity_ ; // gravitational acceleration - // Cumulative flux - Scalar dt_, pa0_, gravity_; - bool p0_defaulted_; Eval W_flux_; - - inline const double area_fraction(const size_t i) - { - return alphai_.at(i); - } - inline void get_influence_table_values(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) + inline void getInfluenceTableValues(Scalar& pitd, Scalar& pitd_prime, const Scalar& td) { // We use the opm-common numeric linear interpolator - pitd = Opm::linearInterpolation(aqutab_td_, aqutab_pi_, td); - pitd_prime = Opm::linearInterpolationDerivative(aqutab_td_, aqutab_pi_, td); + pitd = Opm::linearInterpolation(aquct_data_.td, aquct_data_.pi, td); + pitd_prime = Opm::linearInterpolationDerivative(aquct_data_.td, aquct_data_.pi, td); } - inline void init_quantities(const Aquancon::AquanconOutput& connection) + inline void initQuantities(const Aquancon::AquanconOutput& connection) { // We reset the cumulative flux at the start of any simulation, so, W_flux = 0 W_flux_ = 0.; // We next get our connections to the aquifer and initialize these quantities using the initialize_connections function - initialize_connections(connection); + initializeConnections(connection); - calculate_aquifer_condition(); + calculateAquiferCondition(); + + calculateAquiferConstants(); pressure_previous_.resize(cell_idx_.size(), 0.); pressure_current_.resize(cell_idx_.size(), 0.); Qai_.resize(cell_idx_.size(), 0.0); } - inline void get_current_Pressure_cell(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) + inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); pressure_water.at(idx) = fs.pressure(waterPhaseIdx); } - inline void get_current_density_cell(std::vector& rho_water, const int idx, const IntensiveQuantities& intQuants) + inline void updateCellPressure(std::vector& pressure_water, const int idx, const IntensiveQuantities& intQuants) { const auto& fs = intQuants.fluidState(); - rho_water.at(idx) = fs.density(waterPhaseIdx); + pressure_water.at(idx) = fs.pressure(waterPhaseIdx).value(); + } + + inline void updateCellDensity(const int idx, const IntensiveQuantities& intQuants) + { + const auto& fs = intQuants.fluidState(); + rhow_.at(idx) = fs.density(waterPhaseIdx); } inline Scalar dpai(int idx) { - Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_) - pressure_previous_.at(idx).value(); + Scalar dp = pa0_ + rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0) - pressure_previous_.at(idx); return dp; } // This function implements Eqs 5.8 and 5.9 of the EclipseTechnicalDescription - inline void calculate_a_b_constants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) + inline void calculateEqnConstants(Scalar& a, Scalar& b, const int idx, const SimulatorTimerInterface& timer) { - Scalar beta = aquifer_influx_constant(); - Scalar Tc = time_constant(); - Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc; - Scalar td = timer.simulationTimeElapsed() / Tc; + const Scalar td_plus_dt = (timer.currentStepLength() + timer.simulationTimeElapsed()) / Tc_; + const Scalar td = timer.simulationTimeElapsed() / Tc_; Scalar PItdprime = 0.; Scalar PItd = 0.; - get_influence_table_values(PItd, PItdprime, td_plus_dt); - a = 1.0/Tc * ( (beta * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); - b = beta / (Tc * ( PItd - td*PItdprime)); + getInfluenceTableValues(PItd, PItdprime, td_plus_dt); + a = 1.0/Tc_ * ( (beta_ * dpai(idx)) - (W_flux_.value() * PItdprime) ) / ( PItd - td*PItdprime ); + b = beta_ / (Tc_ * ( PItd - td*PItdprime)); } // This function implements Eq 5.7 of the EclipseTechnicalDescription - inline void calculate_inflow_rate(int idx, const SimulatorTimerInterface& timer) + inline void calculateInflowRate(int idx, const SimulatorTimerInterface& timer) { Scalar a, b; - calculate_a_b_constants(a,b,idx,timer); - Qai_.at(idx) = area_fraction(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx).value() ) ); + calculateEqnConstants(a,b,idx,timer); + Qai_.at(idx) = alphai_.at(idx)*( a - b * ( pressure_current_.at(idx) - pressure_previous_.at(idx) ) ); } - inline const Scalar time_constant() const + inline void calculateAquiferConstants() { - Scalar Tc = mu_w_*phi_aq_*C_t_*r_o_*r_o_/(k_a_*c1_); - return Tc; - } - - inline const Scalar aquifer_influx_constant() const - { - Scalar beta = c2_*h_*theta_*phi_aq_*C_t_*r_o_*r_o_; - return beta; + // We calculate the influx constant + beta_ = aquct_data_.c2 * aquct_data_.h + * aquct_data_.theta * aquct_data_.phi_aq + * aquct_data_.C_t + * aquct_data_.r_o * aquct_data_.r_o; + // We calculate the time constant + Tc_ = mu_w_ * aquct_data_.phi_aq + * aquct_data_.C_t + * aquct_data_.r_o * aquct_data_.r_o + / ( aquct_data_.k_a * aquct_data_.c1 ); } // This function is used to initialize and calculate the alpha_i for each grid connection to the aquifer - inline void initialize_connections(const Aquancon::AquanconOutput& connection) + inline void initializeConnections(const Aquancon::AquanconOutput& connection) { const auto& eclState = ebos_simulator_.vanguard().eclState(); const auto& ugrid = ebos_simulator_.vanguard().grid(); @@ -293,7 +241,7 @@ namespace Opm assert( (connection.influx_multiplier.size() == connection.reservoir_face_dir.size()) ); // We hack the cell depth values for now. We can actually get it from elementcontext pos - cell_depth_.resize(cell_idx_.size(), d0_); + cell_depth_.resize(cell_idx_.size(), aquct_data_.d0); alphai_.resize(cell_idx_.size(), 1.0); faceArea_connected_.resize(cell_idx_.size(),0.0); Scalar faceArea; @@ -301,14 +249,6 @@ namespace Opm auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid); auto faceCells = Opm::AutoDiffGrid::faceCells(ugrid); - for (auto influxCoeff: connection.influx_coeff){ - std::cout << "influx_coeff = " << influxCoeff << std::endl; - } - - for (auto influxMult: connection.influx_multiplier){ - std::cout << "influx_multiplier = " << influxMult << std::endl; - } - // Translate the C face tag into the enum used by opm-parser's TransMult class Opm::FaceDir::DirEnum faceDirection; @@ -326,19 +266,22 @@ namespace Opm // the logically-Cartesian direction of the face const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter); - - if (faceTag == 0) // left - faceDirection = Opm::FaceDir::XMinus; - else if (faceTag == 1) // right - faceDirection = Opm::FaceDir::XPlus; - else if (faceTag == 2) // back - faceDirection = Opm::FaceDir::YMinus; - else if (faceTag == 3) // front - faceDirection = Opm::FaceDir::YPlus; - else if (faceTag == 4) // bottom - faceDirection = Opm::FaceDir::ZMinus; - else if (faceTag == 5) // top - faceDirection = Opm::FaceDir::ZPlus; + switch(faceTag) + { + case 0: faceDirection = Opm::FaceDir::XMinus; + break; + case 1: faceDirection = Opm::FaceDir::XPlus; + break; + case 2: faceDirection = Opm::FaceDir::YMinus; + break; + case 3: faceDirection = Opm::FaceDir::YPlus; + break; + case 4: faceDirection = Opm::FaceDir::ZMinus; + break; + case 5: faceDirection = Opm::FaceDir::ZPlus; + break; + default: OPM_THROW(Opm::NumericalIssue,"Initialization of Aquifer Carter Tracy problem. Make sure faceTag is correctly defined"); + } if (faceDirection == connection.reservoir_face_dir.at(idx)) { @@ -360,16 +303,20 @@ namespace Opm } } - inline void calculate_aquifer_condition() + inline void calculateAquiferCondition() { - int pvttableIdx = pvttableID_ - 1; + int pvttableIdx = aquct_data_.pvttableID - 1; rhow_.resize(cell_idx_.size(),0.); - if (p0_defaulted_) + if (aquct_data_.p0 < 1.0) { - pa0_ = calculate_reservoir_equilibrium(rhow_); + pa0_ = calculateReservoirEquilibrium(); + } + else + { + pa0_ = aquct_data_.p0; } // Initialize a FluidState object first @@ -388,10 +335,11 @@ namespace Opm } // This function is for calculating the aquifer properties from equilibrium state with the reservoir - inline Scalar calculate_reservoir_equilibrium(std::vector& rho_water_reservoir) + inline Scalar calculateReservoirEquilibrium() { // Since the global_indices are the reservoir index, we just need to extract the fluidstate at those indices - std::vector water_pressure_reservoir, pw_aquifer; + std::vector pw_aquifer; + Scalar water_pressure_reservoir; for (size_t idx = 0; idx < cell_idx_.size(); ++idx) { @@ -399,9 +347,9 @@ namespace Opm const auto& intQuants = *(ebos_simulator_.model().cachedIntensiveQuantities(cellIDx, /*timeIdx=*/ 0)); const auto& fs = intQuants.fluidState(); - water_pressure_reservoir.push_back( fs.pressure(waterPhaseIdx).value() ); - rho_water_reservoir.at(idx) = fs.density(waterPhaseIdx); - pw_aquifer.push_back( (water_pressure_reservoir.at(idx) - rho_water_reservoir.at(idx).value()*gravity_*(cell_depth_.at(idx) - d0_))*area_fraction(idx) ); + water_pressure_reservoir = fs.pressure(waterPhaseIdx).value(); + rhow_.at(idx) = fs.density(waterPhaseIdx); + pw_aquifer.push_back( (water_pressure_reservoir - rhow_.at(idx).value()*gravity_*(cell_depth_.at(idx) - aquct_data_.d0))*alphai_.at(idx) ); } // We take the average of the calculated equilibrium pressures. diff --git a/opm/autodiff/BlackoilAquiferModel.hpp b/opm/autodiff/BlackoilAquiferModel.hpp index a2d8a5b59..a10f4500e 100644 --- a/opm/autodiff/BlackoilAquiferModel.hpp +++ b/opm/autodiff/BlackoilAquiferModel.hpp @@ -24,40 +24,12 @@ #ifndef OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED #define OPM_BLACKOILAQUIFERMODEL_HEADER_INCLUDED -#include - -#include -#include - -#include -#include - -#include - #include #include - -#include - #include - -#include -#include -#include -#include #include - -#include -#include -#include - -#include -#include -#include - #include - namespace Opm { /// Class for handling the blackoil well model. @@ -67,22 +39,13 @@ namespace Opm { public: // --------- Types --------- - typedef BlackoilModelParameters ModelParameters; - - - typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; - typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; - typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; - typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices; typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef AquiferCarterTracy Aquifer_object; - BlackoilAquiferModel(Simulator& ebosSimulator, - const ModelParameters& param, - const bool terminal_output); + explicit BlackoilAquiferModel(Simulator& ebosSimulator); // compute the well fluxes and assemble them in to the reservoir equations as source terms // and in the well equations. @@ -92,36 +55,22 @@ namespace Opm { // called at the end of a time step void timeStepSucceeded(const SimulatorTimerInterface& timer); - inline const Simulator& simulator() const - { - return ebosSimulator_; - } - - // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy - void init(const Simulator& ebosSimulator, std::vector& aquifers); - protected: Simulator& ebosSimulator_; - const ModelParameters param_; - bool terminal_output_; - - double gravity_; std::vector aquifers_; + // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy + void init(); void updateConnectionIntensiveQuantities() const; - int numAquifers() const; - void assembleAquiferEq(const SimulatorTimerInterface& timer); // at the beginning of each time step (Not report step) void prepareTimeStep(const SimulatorTimerInterface& timer); - const std::vector& aquifers(); - }; diff --git a/opm/autodiff/BlackoilAquiferModel_impl.hpp b/opm/autodiff/BlackoilAquiferModel_impl.hpp index be1493550..9671c1691 100644 --- a/opm/autodiff/BlackoilAquiferModel_impl.hpp +++ b/opm/autodiff/BlackoilAquiferModel_impl.hpp @@ -3,19 +3,10 @@ namespace Opm { template BlackoilAquiferModel:: - BlackoilAquiferModel(Simulator& ebosSimulator, - const ModelParameters& param, - const bool terminal_output) + BlackoilAquiferModel(Simulator& ebosSimulator) : ebosSimulator_(ebosSimulator) - , param_(param) - , terminal_output_(terminal_output) { - // const auto& gridView = ebosSimulator_.gridView(); - - // number_of_cells_ = gridView.size(/*codim=*/0); - // global_nc_ = gridView.comm().sum(number_of_cells_); - gravity_ = ebosSimulator_.problem().gravity()[2]; - init(ebosSimulator_, aquifers_); + init(); } @@ -26,7 +17,7 @@ namespace Opm { { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->after_time_step(timer); + aquifer->afterTimeStep(timer); } } @@ -65,16 +56,6 @@ namespace Opm { } } - - // Protected function: Return number of aquifers in the model. - template - int - BlackoilAquiferModel:: numAquifers() const - { - return aquifers_.size(); - } - - // Protected function which calls the individual aquifer models template void @@ -82,7 +63,7 @@ namespace Opm { { for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->assembleAquiferEq(ebosSimulator_, timer); + aquifer->assembleAquiferEq(timer); } } @@ -95,43 +76,34 @@ namespace Opm { // Here we can ask each carter tracy aquifers to get the current previous time step's pressure for (auto aquifer = aquifers_.begin(); aquifer != aquifers_.end(); ++aquifer) { - aquifer->before_time_step(ebosSimulator_, timer); + aquifer->beforeTimeStep(timer); } } - // Protected function: Returns a reference to the aquifers members in the model - template - const std::vector< AquiferCarterTracy >& - BlackoilAquiferModel:: aquifers() - { - return aquifers_; - } - - // Initialize the aquifers in the deck template void - BlackoilAquiferModel:: init(const Simulator& ebosSimulator, std::vector< AquiferCarterTracy >& aquifers) + BlackoilAquiferModel:: init() { updateConnectionIntensiveQuantities(); - const auto& deck = ebosSimulator.vanguard().deck(); - const auto& eclState = ebosSimulator.vanguard().eclState(); + const auto& deck = ebosSimulator_.vanguard().deck(); + const auto& eclState = ebosSimulator_.vanguard().eclState(); // Get all the carter tracy aquifer properties data and put it in aquifers vector - AquiferCT aquiferct = AquiferCT(eclState,deck); - Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); + const AquiferCT aquiferct = AquiferCT(eclState,deck); + const Aquancon aquifer_connect = Aquancon(eclState.getInputGrid(), deck); std::vector aquifersData = aquiferct.getAquifers(); std::vector aquifer_connection = aquifer_connect.getAquOutput(); - assert( aquifersData.size() == aquifer_connect.size() ); + assert( aquifersData.size() == aquifer_connection.size() ); for (size_t i = 0; i < aquifersData.size(); ++i) { - aquifers.push_back( - AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), gravity_, ebosSimulator_) - ); + aquifers_.push_back( + AquiferCarterTracy (aquifersData.at(i), aquifer_connection.at(i), ebosSimulator_) + ); } } diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 4d91fa59d..d8a8fcfde 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -30,8 +30,6 @@ #include #include #include -#include -#include #include #include #include @@ -349,7 +347,7 @@ namespace Opm { const ReservoirState& reservoir_state, WellState& well_state) { - // DUNE_UNUSED_PARAMETER(timer); + DUNE_UNUSED_PARAMETER(timer); DUNE_UNUSED_PARAMETER(reservoir_state); DUNE_UNUSED_PARAMETER(well_state); @@ -372,18 +370,18 @@ namespace Opm { ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - // -------- Well and aquifer common variables ---------- - double dt = timer.currentStepLength(); - + // -------- Current time step length ---------- + const double dt = timer.currentStepLength(); + // -------- Aquifer models ---------- try { // Modify the Jacobian and residuals according to the aquifer models aquiferModel().assemble(timer, iterationIdx); } - catch( const Dune::FMatrixError& e ) + catch( ... ) { - OPM_THROW(Opm::NumericalProblem,"Error when assembling aquifer models"); + OPM_THROW(Opm::NumericalIssue,"Error when assembling aquifer models"); } // -------- Well equations ---------- diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index bfbbc0e6b..674f14918 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -30,8 +30,6 @@ #include #include #include -#include -#include #include #include #include @@ -191,7 +189,7 @@ public: ebosSimulator_.model().addAuxiliaryModule(auxMod); } - AquiferModel aquifer_model(ebosSimulator_, model_param_, terminal_output_); + AquiferModel aquifer_model(ebosSimulator_); // Main simulation loop. while (!timer.done()) { From 963d42aab3b4ea6d9f7f12cf4371088f7395327b Mon Sep 17 00:00:00 2001 From: kel85uk Date: Tue, 29 May 2018 14:16:19 +0200 Subject: [PATCH 27/31] Corrects the whitespace again for blackoilmodelebos.hpp --- opm/autodiff/BlackoilModelEbos.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index d8a8fcfde..75b6c3401 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -426,13 +426,13 @@ namespace Opm { if (elem.partitionType() != Dune::InteriorEntity) continue; - unsigned globalElemIdx = elemMapper.index(elem); + unsigned globalElemIdx = elemMapper.index(elem); const auto& priVarsNew = ebosSimulator_.model().solution(/*timeIdx=*/0)[globalElemIdx]; Scalar pressureNew; - pressureNew = priVarsNew[Indices::pressureSwitchIdx]; + pressureNew = priVarsNew[Indices::pressureSwitchIdx]; - Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; + Scalar saturationsNew[FluidSystem::numPhases] = { 0.0 }; Scalar oilSaturationNew = 1.0; if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx]; @@ -475,7 +475,7 @@ namespace Opm { for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++ phaseIdx) { Scalar tmp = saturationsNew[phaseIdx] - saturationsOld[phaseIdx]; - resultDelta += tmp*tmp; + resultDelta += tmp*tmp; resultDenom += saturationsNew[phaseIdx]*saturationsNew[phaseIdx]; } } @@ -483,9 +483,9 @@ namespace Opm { resultDelta = gridView.comm().sum(resultDelta); resultDenom = gridView.comm().sum(resultDenom); - if (resultDenom > 0.0) - return resultDelta/resultDenom; - return 0.0; + if (resultDenom > 0.0) + return resultDelta/resultDenom; + return 0.0; } From 1237724c9a1051ef45b8bf103a6e8e66d0593ab3 Mon Sep 17 00:00:00 2001 From: kel85uk Date: Wed, 30 May 2018 14:26:07 +0200 Subject: [PATCH 28/31] Cleans up the BlackoilModelEbos.hpp file to address the final comments. --- opm/autodiff/BlackoilModelEbos.hpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 75b6c3401..9f37945d0 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -369,9 +369,6 @@ namespace Opm { ebosSimulator_.problem().beginIteration(); ebosSimulator_.model().linearizer().linearize(); ebosSimulator_.problem().endIteration(); - - // -------- Current time step length ---------- - const double dt = timer.currentStepLength(); // -------- Aquifer models ---------- try @@ -384,6 +381,9 @@ namespace Opm { OPM_THROW(Opm::NumericalIssue,"Error when assembling aquifer models"); } + // -------- Current time step length ---------- + const double dt = timer.currentStepLength(); + // -------- Well equations ---------- try @@ -596,6 +596,7 @@ namespace Opm { virtual void apply( const X& x, Y& y ) const { A_.mv( x, y ); + // add well model modification to y wellMod_.apply(x, y ); @@ -609,6 +610,7 @@ namespace Opm { virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const { A_.usmv(alpha,x,y); + // add scaled well model modification to y wellMod_.applyScaleAdd( alpha, x, y ); @@ -1121,9 +1123,6 @@ namespace Opm { BlackoilAquiferModel& aquiferModel() { return aquifer_model_; } - const BlackoilAquiferModel& - aquiferModel() const { return aquifer_model_; } - void beginReportStep() { ebosSimulator_.problem().beginEpisode(); From a0db7c3229db7ae584ccf710bff24d5ac7e138a1 Mon Sep 17 00:00:00 2001 From: Rohith Nair Date: Tue, 29 May 2018 16:35:26 +0200 Subject: [PATCH 29/31] Add test for carter tracy aquifer --- compareECLFiles.cmake | 7 +++++++ tests/update_reference_data.sh | 11 ++++++++++- 2 files changed, 17 insertions(+), 1 deletion(-) mode change 100644 => 100755 compareECLFiles.cmake diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake old mode 100644 new mode 100755 index b01e3c3fd..13b230119 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -195,6 +195,13 @@ add_test_compareECLFiles(CASENAME spe1_thermal REL_TOL ${rel_tol} DIR spe1) +add_test_compareECLFiles(CASENAME ctaquifer_2d_oilwater + FILENAME 2D_OW_CTAQUIFER + SIMULATOR flow + ABS_TOL ${abs_tol} + REL_TOL ${rel_tol} + DIR spe1) + foreach(SIM flow flow_legacy) add_test_compareECLFiles(CASENAME spe3 FILENAME SPE3CASE1 diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh index ab88a21c9..50f30707f 100755 --- a/tests/update_reference_data.sh +++ b/tests/update_reference_data.sh @@ -20,7 +20,7 @@ copyToReferenceDir () { } tests=${@:2} -test -z "$tests" && tests="spe11 spe12 spe12p spe1oilgas spe1nowells spe1thermal spe3 spe5 spe9 norne_init msw_2d_h msw_3d_hfa polymer2d spe9group polymer_oilwater" +test -z "$tests" && tests="spe11 spe12 spe12p spe1oilgas spe1nowells spe1thermal ctaquifer_2d_oilwater spe3 spe5 spe9 norne_init msw_2d_h msw_3d_hfa polymer2d spe9group polymer_oilwater" if grep -q -i "norne " <<< $ghprbCommentBody then if test -d $WORKSPACE/deps/opm-tests/norne/flow @@ -98,6 +98,15 @@ for test_name in ${tests}; do EGRID INIT SMSPEC UNRST UNSMRY fi + if grep -q "ctaquifer_2d_oilwater" <<< $test_name + then + copyToReferenceDir \ + $configuration/build-opm-simulators/tests/results/flow+ctaquifer_2d_oilwater/ \ + $OPM_TESTS_ROOT/spe1/opm-simulation-reference/flow \ + 2D_OW_CTAQUIFER \ + EGRID INIT SMSPEC UNRST UNSMRY + fi + if grep -q "msw_2d_h" <<< $test_name then copyToReferenceDir \ From 123d3c68c5f08d234df50c088991c101e721fa31 Mon Sep 17 00:00:00 2001 From: Rohith Nair Date: Tue, 29 May 2018 16:40:40 +0200 Subject: [PATCH 30/31] edit --- tests/update_reference_data.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh index 50f30707f..8d384e50d 100755 --- a/tests/update_reference_data.sh +++ b/tests/update_reference_data.sh @@ -102,7 +102,7 @@ for test_name in ${tests}; do then copyToReferenceDir \ $configuration/build-opm-simulators/tests/results/flow+ctaquifer_2d_oilwater/ \ - $OPM_TESTS_ROOT/spe1/opm-simulation-reference/flow \ + $OPM_TESTS_ROOT/aquifer-oilwater/opm-simulation-reference/flow \ 2D_OW_CTAQUIFER \ EGRID INIT SMSPEC UNRST UNSMRY fi From 755b3cc814652fed2047922de09a34c8452d9633 Mon Sep 17 00:00:00 2001 From: Rohith Nair Date: Tue, 29 May 2018 16:44:28 +0200 Subject: [PATCH 31/31] edit --- compareECLFiles.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake index 13b230119..c4175c03c 100755 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -200,7 +200,7 @@ add_test_compareECLFiles(CASENAME ctaquifer_2d_oilwater SIMULATOR flow ABS_TOL ${abs_tol} REL_TOL ${rel_tol} - DIR spe1) + DIR aquifer-oilwater) foreach(SIM flow flow_legacy) add_test_compareECLFiles(CASENAME spe3