diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 1f814c52..7290f255 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -161,11 +161,8 @@ list (APPEND EXAMPLE_SOURCE_FILES examples/compute_tof.cpp examples/compute_tof_from_files.cpp examples/import_rewrite.cpp - examples/refine_wells.cpp - examples/scaneclipsedeck.c examples/sim_2p_comp_reorder.cpp examples/sim_2p_incomp.cpp - examples/spu_2p.cpp examples/wells_example.cpp tutorials/tutorial1.cpp tutorials/tutorial2.cpp diff --git a/examples/refine_wells.cpp b/examples/refine_wells.cpp deleted file mode 100644 index 0ddb28fd..00000000 --- a/examples/refine_wells.cpp +++ /dev/null @@ -1,50 +0,0 @@ -/* - Copyright 2012 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include - -// Double I and J coordinates of wells and completions. -// Do not change any well productivity indices or any other data. -int main(int argc, char** argv) -{ - using namespace Opm; - if (argc != 2) { - std::cerr << "Usage: " << argv[0] << " deck\n"; - return 1; - } - EclipseGridParser deck(argv[1], false); - - WELSPECS ws = deck.getWELSPECS(); - const int nw = ws.welspecs.size(); - for (int w = 0; w < nw; ++w) { - ws.welspecs[w].I_ *= 2; - ws.welspecs[w].J_ *= 2; - } - - COMPDAT cd = deck.getCOMPDAT(); - const int nc = cd.compdat.size(); - for (int c = 0; c < nc; ++c) { - cd.compdat[c].grid_ind_[0] *= 2; - cd.compdat[c].grid_ind_[1] *= 2; - } - - ws.write(std::cout); - std::cout << '\n'; - cd.write(std::cout); -} diff --git a/examples/scaneclipsedeck.c b/examples/scaneclipsedeck.c deleted file mode 100644 index e306f9ed..00000000 --- a/examples/scaneclipsedeck.c +++ /dev/null @@ -1,134 +0,0 @@ -/* scaneclipsedec finds technically valid Eclipse keywords in an ascii file. - * Copyright (c) 2010 Jostein R. Natvig - * - * The MIT License - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the "Software"), - * to deal in the Software without restriction, including without limitation - * the rights to use, copy, modify, merge, publish, distribute, sublicense, - * and/or sell copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR - * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, - * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL - * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR - * OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, - * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR - * OTHER DEALINGS IN THE SOFTWARE. - */ -#include -#include -#include -#include - - - -static char* -read_keyword(FILE *fp, char *buf) -{ - int i, j, c; - - /* Clear buf */ - for (i=0; i<9; ++i) { - buf[i] = '\0'; - } - - /* Read first character and check if it is uppercase*/ - buf[0] = fgetc(fp); - if ( !isupper( buf[0] ) ) { - return NULL; /* NOT VALID CHARACTER */ - ungetc(buf[0], fp); - } - - - /* Scan as much as possible possible keyword, 8 characters long */ - i = 1; - while ( (c = fgetc(fp)) && - (c != EOF ) && - (c != '\n' ) && - (c != '/' ) && - (i < 8 )) { - buf[i++] = c; - } - - /* Skip rest of line */ - if (c != '\n'){ - while ( (c = fgetc(fp)) && - (c != EOF ) && - (c != '\n' )) { - ; - } - } - - if(c == '\n') { - ungetc(c, fp); - } - - /* Find first non-uppercase or non-digit character */ - for (i=0; i<8; ++i) { - if ( !(isupper(buf[i]) || isdigit(buf[i])) ) { - break; - } - } - - /* Check if remaining characters are blank */ - for (j = i; j<8; ++j) { - if(!isspace(buf[j]) && buf[j] != '\0') { - return NULL; /* CHARACTER AFTER SPACE OR INVALID CHARACTER */ - } - buf[j] = '\0'; - } - return buf; -} - - -int main(int argc, char *argv[]) -{ - int c, lineno, nkw; - FILE *fp; - char buf[10]; - - if (argc != 2) - { - fprintf(stderr, "Usage: filename.grdecl\n"); - exit(EXIT_FAILURE); - } - fp = fopen(argv[1], "ra"); - if (fp == NULL) - { - fprintf(stderr, "No such file...\n"); - exit(EXIT_FAILURE); - } - - lineno = nkw = 0; - - if (read_keyword(fp, buf) != NULL) { - ++nkw; - fprintf(stderr, "%s\n", buf); - } - - - - while ((c = getc(fp)) != EOF) { /* Eat large chunks */ - if ( c == '\n') { - ++lineno; - - if (read_keyword(fp, buf) != NULL) { - fprintf(stderr, "%s\n", buf); - ++nkw; - } - } - } - - fprintf(stderr, "Scanned %d lines, found %d keywords.\n", lineno, nkw); - - return 0; -} -/* Local Variables: */ -/* c-basic-offset:4 */ -/* End: */ diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp deleted file mode 100644 index d7995626..00000000 --- a/examples/spu_2p.cpp +++ /dev/null @@ -1,716 +0,0 @@ -/*=========================================================================== -// -// File: spu_2p.cpp -// -// Created: 2011-10-05 10:29:01+0200 -// -// Authors: Ingeborg S. Ligaarden -// Jostein R. Natvig -// Halvor M. Nilsen -// Atgeirr F. Rasmussen -// Bård Skaflestad -// -//==========================================================================*/ - - -/* - Copyright 2011, 2012 SINTEF ICT, Applied Mathematics. - Copyright 2011, 2012 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 . -*/ - -#if HAVE_CONFIG_H -#include "config.h" -#endif // HAVE_CONFIG_H - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -//#include -#include -#include -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include - -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -#ifdef HAVE_ERT -#include -#endif - - - -static void outputState(const UnstructuredGrid& grid, - const Opm::TwophaseState& state, - const Opm::SimulatorTimer& simtimer, - const std::string& output_dir) -{ - // Write data in VTK format. - int step = simtimer.currentStepNum(); - std::ostringstream vtkfilename; - vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; - std::ofstream vtkfile(vtkfilename.str().c_str()); - if (!vtkfile) { - THROW("Failed to open " << vtkfilename.str()); - } - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - std::vector cell_velocity; - Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - Opm::writeVtkData(grid, dm, vtkfile); -#ifdef HAVE_ERT - Opm::writeECLData(grid, dm, simtimer.currentStepNum(), simtimer.currentTime(), simtimer.currentDateTime(), output_dir, "OPM" ); -#endif - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; - std::ofstream file(fname.str().c_str()); - if (!file) { - THROW("Failed to open " << fname.str()); - } - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } -} - - -static void outputWaterCut(const Opm::Watercut& watercut, - const std::string& output_dir) -{ - // Write water cut curve. - std::string fname = output_dir + "/watercut.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - THROW("Failed to open " << fname); - } - watercut.write(os); -} - - -static void outputWellReport(const Opm::WellReport& wellreport, - const std::string& output_dir) -{ - // Write well report. - std::string fname = output_dir + "/wellreport.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - THROW("Failed to open " << fname); - } - wellreport.write(os); -} - - - -// --------------- Types needed to define transport solver --------------- - -class SimpleFluid2pWrappingProps -{ -public: - SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props) - : props_(props), - smin_(props.numCells()*props.numPhases()), - smax_(props.numCells()*props.numPhases()) - { - if (props.numPhases() != 2) { - THROW("SimpleFluid2pWrapper requires 2 phases."); - } - const int num_cells = props.numCells(); - std::vector cells(num_cells); - for (int c = 0; c < num_cells; ++c) { - cells[c] = c; - } - props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]); - } - - double density(int phase) const - { - return props_.density()[phase]; - } - - template - void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const - { - props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]); - const double* mu = props_.viscosity(); - mob[0] /= mu[0]; - mob[1] /= mu[1]; - // Recall that we use Fortran ordering for kr derivatives, - // therefore dmob[i*2 + j] is row j and column i of the - // matrix. - // Each row corresponds to a kr function, so which mu to - // divide by also depends on the row, j. - dmob[0*2 + 0] /= mu[0]; - dmob[0*2 + 1] /= mu[1]; - dmob[1*2 + 0] /= mu[0]; - dmob[1*2 + 1] /= mu[1]; - } - - template - void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const - { - double pcow[2]; - double dpcow[4]; - props_.capPress(1, &s[0], &c, pcow, dpcow); - pcap = pcow[0]; - ASSERT(pcow[1] == 0.0); - dpcap = dpcow[0]; - ASSERT(dpcow[1] == 0.0); - ASSERT(dpcow[2] == 0.0); - ASSERT(dpcow[3] == 0.0); - } - - double s_min(int c) const - { - return smin_[2*c + 0]; - } - - double s_max(int c) const - { - return smax_[2*c + 0]; - } - -private: - const Opm::IncompPropertiesInterface& props_; - std::vector smin_; - std::vector smax_; -}; - -typedef SimpleFluid2pWrappingProps TwophaseFluid; -typedef Opm::SinglePointUpwindTwoPhase TransportModel; - -using namespace Opm::ImplicitTransportDefault; - -typedef NewtonVectorCollection< ::std::vector > NVecColl; -typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys; - -template -class MaxNorm { -public: - static double - norm(const Vector& v) { - return AccumulationNorm ::norm(v); - } -}; - -typedef Opm::ImplicitTransport TransportSolver; - - - -// ----------------- Main program ----------------- -int -main(int argc, char** argv) -{ - using namespace Opm; - - std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; - Opm::parameter::ParameterGroup param(argc, argv, false); - std::cout << "--------------- Reading parameters ---------------" << std::endl; - - // Reading various control parameters. - const bool guess_old_solution = param.getDefault("guess_old_solution", false); - const bool use_reorder = param.getDefault("use_reorder", true); - const bool output = param.getDefault("output", true); - std::string output_dir; - int output_interval = 1; - if (output) { - output_dir = param.getDefault("output_dir", std::string("output")); - // Ensure that output dir exists - boost::filesystem::path fpath(output_dir); - try { - create_directories(fpath); - } - catch (...) { - THROW("Creating directories failed: " << fpath); - } - output_interval = param.getDefault("output_interval", output_interval); - } - const int num_transport_substeps = param.getDefault("num_transport_substeps", 1); - - // If we have a "deck_filename", grid and props will be read from that. - bool use_deck = param.has("deck_filename"); - boost::scoped_ptr grid; - boost::scoped_ptr props; - boost::scoped_ptr wells; - boost::scoped_ptr rock_comp; - Opm::SimulatorTimer simtimer; - Opm::TwophaseState state; - bool check_well_controls = false; - int max_well_control_iterations = 0; - double gravity[3] = { 0.0 }; - if (use_deck) { - std::string deck_filename = param.get("deck_filename"); - Opm::EclipseGridParser deck(deck_filename); - // Grid init - grid.reset(new Opm::GridManager(deck)); - // Rock and fluid init - props.reset(new Opm::IncompPropertiesFromDeck(deck, *grid->c_grid())); - // Wells init. - wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); - check_well_controls = param.getDefault("check_well_controls", false); - max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); - // Timer init. - if (deck.hasField("TSTEP")) { - simtimer.init(deck); - } else { - simtimer.init(param); - } - // Rock compressibility. - rock_comp.reset(new Opm::RockCompressibility(deck)); - // Gravity. - gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; - // Init state variables (saturation and pressure). - if (param.has("init_saturation")) { - initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); - } else { - initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); - } - } else { - // Grid init. - const int nx = param.getDefault("nx", 100); - const int ny = param.getDefault("ny", 100); - const int nz = param.getDefault("nz", 1); - const double dx = param.getDefault("dx", 1.0); - const double dy = param.getDefault("dy", 1.0); - const double dz = param.getDefault("dz", 1.0); - grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); - // Rock and fluid init. - props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - // Wells init. - wells.reset(new Opm::WellsManager()); - // Timer init. - simtimer.init(param); - // Rock compressibility. - rock_comp.reset(new Opm::RockCompressibility(param)); - // Gravity. - gravity[2] = param.getDefault("gravity", 0.0); - // Init state variables (saturation and pressure). - initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); - } - - // Extra fluid init for transport solver. - TwophaseFluid fluid(*props); - - // Warn if gravity but no density difference. - bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); - if (use_gravity) { - if (props->density()[0] == props->density()[1]) { - std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; - } - } - bool use_segregation_split = false; - bool use_column_solver = false; - bool use_gauss_seidel_gravity = false; - if (use_gravity && use_reorder) { - use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split); - if (use_segregation_split) { - use_column_solver = param.getDefault("use_column_solver", use_column_solver); - if (use_column_solver) { - use_gauss_seidel_gravity = param.getDefault("use_gauss_seidel_gravity", use_gauss_seidel_gravity); - } - } - } - - // Check that rock compressibility is not used with solvers that do not handle it. - int nl_pressure_maxiter = 0; - double nl_pressure_residual_tolerance = 0.0; - double nl_pressure_change_tolerance = 0.0; - if (rock_comp->isActive()) { - if (!use_reorder) { - THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet."); - } - nl_pressure_residual_tolerance = param.getDefault("nl_pressure_residual_tolerance", 0.0); - nl_pressure_change_tolerance = param.getDefault("nl_pressure_change_tolerance", 1.0); // In Pascal. - nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10); - } - - // Source-related variables init. - int num_cells = grid->c_grid()->number_of_cells; - - // Extra rock init. - std::vector porevol; - if (rock_comp->isActive()) { - computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); - } else { - computePorevolume(*grid->c_grid(), props->porosity(), porevol); - } - double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); - - // Initialising src - std::vector src(num_cells, 0.0); - if (wells->c_wells()) { - // Do nothing, wells will be the driving force, not source terms. - // Opm::wellsToSrc(*wells->c_wells(), num_cells, src); - } else { - const double default_injection = use_gravity ? 0.0 : 0.1; - const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) - *tot_porevol_init/Opm::unit::day; - src[0] = flow_per_sec; - src[num_cells - 1] = -flow_per_sec; - } - - TransportSource* tsrc = create_transport_source(2, 2); - double ssrc[] = { 1.0, 0.0 }; - double ssink[] = { 0.0, 1.0 }; - double zdummy[] = { 0.0, 0.0 }; - for (int cell = 0; cell < num_cells; ++cell) { - if (src[cell] > 0.0) { - append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc); - } else if (src[cell] < 0.0) { - append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc); - } - } - std::vector reorder_src = src; - - // Boundary conditions. - Opm::FlowBCManager bcs; - if (param.getDefault("use_pside", false)) { - int pside = param.get("pside"); - double pside_pressure = param.get("pside_pressure"); - bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure); - } - - // Solvers init. - // Linear solver. - Opm::LinearSolverFactory linsolver(param); - // Pressure solver. - const double *grav = use_gravity ? &gravity[0] : 0; - Opm::IncompTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver, - nl_pressure_residual_tolerance, nl_pressure_change_tolerance, - nl_pressure_maxiter, - grav, wells->c_wells(), src, bcs.c_bcs()); - // Reordering solver. - const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); - const int nl_maxiter = param.getDefault("nl_maxiter", 30); - const double* transport_grav = use_gauss_seidel_gravity ? grav : NULL; - Opm::TransportSolverTwophaseReorder reorder_model(*grid->c_grid(), *props, transport_grav, nl_tolerance, nl_maxiter); - - // Non-reordering solver. - TransportModel model (fluid, *grid->c_grid(), porevol, grav, guess_old_solution); - if (use_gravity) { - model.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans()); - } - TransportSolver tsolver(model); - // Column-based gravity segregation solver. - std::vector > columns; - if (use_column_solver) { - Opm::extractColumn(*grid->c_grid(), columns); - } - Opm::GravityColumnSolver colsolver(model, *grid->c_grid(), nl_tolerance, nl_maxiter); - - // Control init. - Opm::ImplicitTransportDetails::NRReport rpt; - Opm::ImplicitTransportDetails::NRControl ctrl; - if (!use_reorder || use_segregation_split) { - ctrl.max_it = param.getDefault("max_it", 20); - ctrl.verbosity = param.getDefault("verbosity", 0); - ctrl.max_it_ls = param.getDefault("max_it_ls", 5); - } - - // Linear solver init. - using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver; - CSRMatrixUmfpackSolver linsolve; - - // The allcells vector is used in calls to computeTotalMobility() - // and computeTotalMobilityOmega(). - std::vector allcells(num_cells); - for (int cell = 0; cell < num_cells; ++cell) { - allcells[cell] = cell; - } - - // Warn if any parameters are unused. - if (param.anyUnused()) { - std::cout << "-------------------- Unused parameters: --------------------\n"; - param.displayUsage(); - std::cout << "----------------------------------------------------------------" << std::endl; - } - - // Write parameters used for later reference. - if (output) { - param.writeParam(output_dir + "/simulation.param"); - } - - // Main simulation loop. - Opm::time::StopWatch pressure_timer; - double ptime = 0.0; - Opm::time::StopWatch transport_timer; - double ttime = 0.0; - Opm::time::StopWatch total_timer; - total_timer.start(); - std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl; - double init_satvol[2] = { 0.0 }; - double satvol[2] = { 0.0 }; - double injected[2] = { 0.0 }; - double produced[2] = { 0.0 }; - double tot_injected[2] = { 0.0 }; - double tot_produced[2] = { 0.0 }; - Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol); - std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init - << " " << init_satvol[1]/tot_porevol_init << std::endl; - Opm::Watercut watercut; - watercut.push(0.0, 0.0, 0.0); - Opm::WellReport wellreport; - Opm::WellState well_state; - well_state.init(wells->c_wells(), state); - std::vector fractional_flows; - std::vector well_resflows_phase; - int num_wells = 0; - if (wells->c_wells()) { - num_wells = wells->c_wells()->number_of_wells; - well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); - wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); - } - for (; !simtimer.done(); ++simtimer) { - // Report timestep and (optionally) write state to disk. - simtimer.report(std::cout); - if (output && (simtimer.currentStepNum() % output_interval == 0)) { - outputState(*grid->c_grid(), state, simtimer , output_dir); - } - - // Solve pressure. - if (check_well_controls) { - computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows); - } - if (check_well_controls) { - wells->applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); - } - bool well_control_passed = !check_well_controls; - int well_control_iteration = 0; - do { - pressure_timer.start(); - std::vector initial_pressure = state.pressure(); - psolver.solve(simtimer.currentStepLength(), state, well_state); - if (!rock_comp->isActive()) { - // Compute average pressures of previous and last - // step, and total volume. - double av_prev_press = 0.; - double av_press = 0.; - double tot_vol = 0.; - for (int cell = 0; cell < num_cells; ++cell) { - av_prev_press += initial_pressure[cell]*grid->c_grid()->cell_volumes[cell]; - av_press += state.pressure()[cell]*grid->c_grid()->cell_volumes[cell]; - tot_vol += grid->c_grid()->cell_volumes[cell]; - } - // Renormalization constant - const double ren_const = (av_prev_press - av_press)/tot_vol; - for (int cell = 0; cell < num_cells; ++cell) { - state.pressure()[cell] += ren_const; - } - for (int well = 0; well < num_wells; ++well) { - well_state.bhp()[well] += ren_const; - } - } - pressure_timer.stop(); - double pt = pressure_timer.secsSinceStart(); - std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; - ptime += pt; - - - if (check_well_controls) { - Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), - fractional_flows, - well_state.perfRates(), - well_resflows_phase); - std::cout << "Checking well conditions." << std::endl; - // For testing we set surface := reservoir - well_control_passed = wells->conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); - ++well_control_iteration; - if (!well_control_passed && well_control_iteration > max_well_control_iterations) { - THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); - } - if (!well_control_passed) { - std::cout << "Well controls not passed, solving again." << std::endl; - } else { - std::cout << "Well conditions met." << std::endl; - } - } - } while (!well_control_passed); - - // Update pore volumes if rock is compressible. - if (rock_comp->isActive()) { - computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); - } - - // Process transport sources (to include bdy terms and well flows). - Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, - wells->c_wells(), well_state.perfRates(), reorder_src); - if (!use_reorder) { - clear_transport_source(tsrc); - for (int cell = 0; cell < num_cells; ++cell) { - if (reorder_src[cell] > 0.0) { - append_transport_source(cell, 2, 0, reorder_src[cell], ssrc, zdummy, tsrc); - } else if (reorder_src[cell] < 0.0) { - append_transport_source(cell, 2, 0, reorder_src[cell], ssink, zdummy, tsrc); - } - } - } - - // Solve transport. - transport_timer.start(); - double stepsize = simtimer.currentStepLength(); - if (num_transport_substeps != 1) { - stepsize /= double(num_transport_substeps); - std::cout << "Making " << num_transport_substeps << " transport substeps." << std::endl; - } - for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) { - if (use_reorder) { - reorder_model.solve(&porevol[0], &reorder_src[0], stepsize, state); - Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); - if (use_segregation_split) { - if (use_column_solver) { - if (use_gauss_seidel_gravity) { - reorder_model.solveGravity(&porevol[0], stepsize, state); - } else { - colsolver.solve(columns, stepsize, state.saturation()); - } - } else { - std::vector fluxes = state.faceflux(); - std::fill(state.faceflux().begin(), state.faceflux().end(), 0.0); - tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt); - std::cout << rpt; - state.faceflux() = fluxes; - } - } - } else { - tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt); - std::cout << rpt; - Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced); - } - } - transport_timer.stop(); - double tt = transport_timer.secsSinceStart(); - std::cout << "Transport solver took: " << tt << " seconds." << std::endl; - ttime += tt; - - // Report volume balances. - Opm::computeSaturatedVol(porevol, state.saturation(), satvol); - tot_injected[0] += injected[0]; - tot_injected[1] += injected[1]; - tot_produced[0] += produced[0]; - tot_produced[1] += produced[1]; - std::cout.precision(5); - const int width = 18; - std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n"; - std::cout << " Saturated volumes: " - << std::setw(width) << satvol[0]/tot_porevol_init - << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; - std::cout << " Injected volumes: " - << std::setw(width) << injected[0]/tot_porevol_init - << std::setw(width) << injected[1]/tot_porevol_init << std::endl; - std::cout << " Produced volumes: " - << std::setw(width) << produced[0]/tot_porevol_init - << std::setw(width) << produced[1]/tot_porevol_init << std::endl; - std::cout << " Total inj volumes: " - << std::setw(width) << tot_injected[0]/tot_porevol_init - << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; - std::cout << " Total prod volumes: " - << std::setw(width) << tot_produced[0]/tot_porevol_init - << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; - std::cout << " In-place + prod - inj: " - << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init - << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl; - std::cout << " Init - now - pr + inj: " - << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init - << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init - << std::endl; - std::cout.precision(8); - - watercut.push(simtimer.currentTime() + simtimer.currentStepLength(), - produced[0]/(produced[0] + produced[1]), - tot_produced[0]/tot_porevol_init); - if (wells->c_wells()) { - wellreport.push(*props, *wells->c_wells(), state.saturation(), - simtimer.currentTime() + simtimer.currentStepLength(), - well_state.bhp(), well_state.perfRates()); - } - } - total_timer.stop(); - - std::cout << "\n\n================ End of simulation ===============\n" - << "Total time taken: " << total_timer.secsSinceStart() - << "\n Pressure time: " << ptime - << "\n Transport time: " << ttime << std::endl; - - if (output) { - outputState(*grid->c_grid(), state, simtimer, output_dir); - outputWaterCut(watercut, output_dir); - if (wells->c_wells()) { - outputWellReport(wellreport, output_dir); - } - } - - destroy_transport_source(tsrc); -}