From 10e896fa6bcc0f6ae71c2b28ad256656c4c96e8f Mon Sep 17 00:00:00 2001 From: kel85uk Date: Mon, 4 Dec 2017 22:26:53 +0100 Subject: [PATCH] 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_));