From d171f043e5d7f7f97006642c5d4116be6d5e5978 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 22 May 2013 22:53:40 +0200 Subject: [PATCH] Added simulator class and program using ImpesTPFAAD. This copies similar code from opm-core and replaces CompressibleTpfa with ImpesTPFAAD. Produces same results for initial simple test cases. --- CMakeLists_files.cmake | 3 + examples/sim_2p_comp_ad.cpp | 285 ++++++++++++ opm/autodiff/SimulatorCompressibleAd.cpp | 545 +++++++++++++++++++++++ opm/autodiff/SimulatorCompressibleAd.hpp | 99 ++++ 4 files changed, 932 insertions(+) create mode 100644 examples/sim_2p_comp_ad.cpp create mode 100644 opm/autodiff/SimulatorCompressibleAd.cpp create mode 100644 opm/autodiff/SimulatorCompressibleAd.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 29f89feb3..c4bc41b60 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -29,6 +29,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/BlackoilPropsAd.cpp opm/autodiff/BlackoilPropsAdInterface.cpp opm/autodiff/ImpesTPFAAD.cpp + opm/autodiff/SimulatorCompressibleAd.cpp opm/autodiff/SimulatorIncompTwophaseAdfi.cpp opm/autodiff/TransportSolverTwophaseAd.cpp ) @@ -49,6 +50,7 @@ list (APPEND TEST_DATA_FILES # find tutorials examples -name '*.c*' -printf '\t%p\n' | sort list (APPEND EXAMPLE_SOURCE_FILES examples/find_zero.cpp + examples/sim_2p_comp_ad.cpp examples/sim_2p_incomp_adfi.cpp examples/sim_simple.cpp examples/test_impestpfa_ad.cpp @@ -70,6 +72,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/BlackoilPropsAdInterface.hpp opm/autodiff/GeoProps.hpp opm/autodiff/ImpesTPFAAD.hpp + opm/autodiff/SimulatorCompressibleAd.hpp opm/autodiff/SimulatorIncompTwophaseAdfi.hpp opm/autodiff/TransportSolverTwophaseAd.hpp ) diff --git a/examples/sim_2p_comp_ad.cpp b/examples/sim_2p_comp_ad.cpp new file mode 100644 index 000000000..03074c10c --- /dev/null +++ b/examples/sim_2p_comp_ad.cpp @@ -0,0 +1,285 @@ +/* + Copyright 2013 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 . +*/ + + +#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 + + +namespace +{ + void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) + { + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + } +} // anon namespace + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + std::cout << "\n================ Test program for weakly compressible two-phase flow ===============\n\n"; + parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr deck; + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr rock_comp; + BlackoilState 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"); + deck.reset(new EclipseGridParser(deck_filename)); + // Grid init + grid.reset(new GridManager(*deck)); + // Rock and fluid init + props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param)); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(*deck)); + // Gravity. + gravity[2] = deck->hasField("NOGRAV") ? 0.0 : 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); + } + initBlackoilSurfvol(*grid->c_grid(), *props, 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 GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Rock compressibility. + rock_comp.reset(new 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); + initBlackoilSurfvol(*grid->c_grid(), *props, state); + } + + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + const double *grav = use_gravity ? &gravity[0] : 0; + + // Initialising src + int num_cells = grid->c_grid()->number_of_cells; + std::vector src(num_cells, 0.0); + if (use_deck) { + // Do nothing, wells will be the driving force, not source terms. + } else { + // Compute pore volumes, in order to enable specifying injection rate + // terms of total pore volume. + 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); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + 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/unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + 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(), FlowBCManager::Side(pside), pside_pressure); + } + + // Linear solver. + LinearSolverFactory linsolver(param); + + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::ofstream epoch_os; + std::string output_dir; + if (output) { + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + std::string filename = output_dir + "/epoch_timing.param"; + epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out); + // open file to clean it. The file is appended to in SimulatorTwophase + filename = output_dir + "/step_timing.param"; + std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out); + step_os.close(); + param.writeParam(output_dir + "/simulation.param"); + } + + + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << " (number of epochs: " + << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + + SimulatorReport rep; + if (!use_deck) { + // Simple simulation without a deck. + WellsManager wells; // no wells. + SimulatorCompressibleAd simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); + SimulatorTimer simtimer; + simtimer.init(param); + warnIfUnusedParams(param); + WellState well_state; + well_state.init(0, state); + rep = simulator.run(simtimer, state, well_state); + } else { + // With a deck, we may have more epochs etc. + WellState well_state; + int step = 0; + SimulatorTimer simtimer; + // Use timer for last epoch to obtain total time. + deck->setCurrentEpoch(deck->numberOfEpochs() - 1); + simtimer.init(*deck); + const double total_time = simtimer.totalTime(); + for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { + // Set epoch index. + deck->setCurrentEpoch(epoch); + + // Update the timer. + if (deck->hasField("TSTEP")) { + simtimer.init(*deck); + } else { + if (epoch != 0) { + THROW("No TSTEP in deck for epoch " << epoch); + } + simtimer.init(param); + } + simtimer.setCurrentStepNum(step); + simtimer.setTotalTime(total_time); + + // Report on start of epoch. + std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" + << "\n (number of steps: " + << simtimer.numSteps() - step << ")\n\n" << std::flush; + + // Create new wells, well_state + WellsManager wells(*deck, *grid->c_grid(), props->permeability()); + // @@@ HACK: we should really make a new well state and + // properly transfer old well state to it every epoch, + // since number of wells may change etc. + if (epoch == 0) { + well_state.init(wells.c_wells(), state); + } + + // Create and run simulator. + SimulatorCompressibleAd simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells, + src, + bcs.c_bcs(), + linsolver, + grav); + if (epoch == 0) { + warnIfUnusedParams(param); + } + SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); + if (output) { + epoch_rep.reportParam(epoch_os); + } + // Update total timing report and remember step number. + rep += epoch_rep; + step = simtimer.currentStepNum(); + } + } + + std::cout << "\n\n================ End of simulation ===============\n\n"; + rep.report(std::cout); + + if (output) { + std::string filename = output_dir + "/walltime.param"; + std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); + rep.reportParam(tot_os); + } + +} diff --git a/opm/autodiff/SimulatorCompressibleAd.cpp b/opm/autodiff/SimulatorCompressibleAd.cpp new file mode 100644 index 000000000..83315f3da --- /dev/null +++ b/opm/autodiff/SimulatorCompressibleAd.cpp @@ -0,0 +1,545 @@ +/* + Copyright 2013 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 . +*/ + +#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 + + +namespace Opm +{ + + class SimulatorCompressibleAd::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity); + + SimulatorReport run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state); + + private: + // Data. + + // Parameters for output. + bool output_; + bool output_vtk_; + std::string output_dir_; + int output_interval_; + // Parameters for well control + bool check_well_controls_; + int max_well_control_iterations_; + // Parameters for transport solver. + int num_transport_substeps_; + bool use_segregation_split_; + // Observed objects. + const UnstructuredGrid& grid_; + const BlackoilPropertiesInterface& props_; + const RockCompressibility* rock_comp_props_; + WellsManager& wells_manager_; + const Wells* wells_; + const std::vector& src_; + const FlowBoundaryConditions* bcs_; + const double* gravity_; + // Solvers + BlackoilPropsAd fluid_; + DerivedGeology geo_; + ImpesTPFAAD psolver_; + TransportSolverCompressibleTwophaseReorder tsolver_; + // Needed by column-based gravity segregation solver. + std::vector< std::vector > columns_; + // Misc. data + std::vector allcells_; + }; + + + + + SimulatorCompressibleAd::SimulatorCompressibleAd(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity)); + } + + + + + + SimulatorReport SimulatorCompressibleAd::run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state) + { + return pimpl_->run(timer, state, well_state); + } + + + + static void outputStateVtk(const UnstructuredGrid& grid, + const Opm::BlackoilState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + vtkfilename << "/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); + } + + + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::BlackoilState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["surfvolume"] = &state.surfacevol(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // 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; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + file.precision(15); + 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); + } + + + + // \TODO: make CompressibleTpfa take src and bcs. + SimulatorCompressibleAd::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + rock_comp_props_(rock_comp_props), + wells_manager_(wells_manager), + wells_(wells_manager.c_wells()), + src_(src), + bcs_(bcs), + gravity_(gravity), + fluid_(props_), + geo_(grid_, fluid_, gravity_), + psolver_(grid_, fluid_, geo_, *wells_manager.c_wells(), linsolver), + /* param.getDefault("nl_pressure_residual_tolerance", 0.0), + param.getDefault("nl_pressure_change_tolerance", 1.0), + param.getDefault("nl_pressure_maxiter", 10), + gravity, */ + tsolver_(grid, props, + param.getDefault("nl_tolerance", 1e-9), + param.getDefault("nl_maxiter", 30)) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); + 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", 1); + } + + // Well control related init. + check_well_controls_ = param.getDefault("check_well_controls", false); + max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); + + // Transport related init. + num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); + use_segregation_split_ = param.getDefault("use_segregation_split", false); + if (gravity != 0 && use_segregation_split_){ + tsolver_.initGravity(gravity); + extractColumn(grid_, columns_); + } + + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + SimulatorReport SimulatorCompressibleAd::Impl::run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state) + { + std::vector transport_src; + + // Initialisation. + std::vector porevol; + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } else { + computePorevolume(grid_, props_.porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + std::vector initial_porevol = porevol; + + // Main simulation loop. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + double init_surfvol[2] = { 0.0 }; + double inplace_surfvol[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol); + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); + Opm::WellReport wellreport; + std::vector fractional_flows; + std::vector well_resflows_phase; + if (wells_) { + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + wellreport.push(props_, *wells_, + state.pressure(), state.surfacevol(), state.saturation(), + 0.0, well_state.bhp(), well_state.perfRates()); + } + std::fstream tstep_os; + if (output_) { + std::string filename = output_dir_ + "/step_timing.param"; + tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); + } + for (; !timer.done(); ++timer) { + // Report timestep and (optionally) write state to disk. + step_timer.start(); + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + } + + SimulatorReport sreport; + + // Solve pressure equation. + if (check_well_controls_) { + computeFractionalFlow(props_, allcells_, + state.pressure(), state.surfacevol(), state.saturation(), + fractional_flows); + wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); + } + bool well_control_passed = !check_well_controls_; + int well_control_iteration = 0; + do { + // Run solver. + pressure_timer.start(); + std::vector initial_pressure = state.pressure(); + psolver_.solve(timer.currentStepLength(), state, well_state); + +#if 0 + // Renormalize pressure if both fluids and rock are + // incompressible, and there are no pressure + // conditions (bcs or wells). It is deemed sufficient + // for now to renormalize using geometric volume + // instead of pore volume. + if (psolver_.singularPressure()) { + // Compute average pressures of previous and last + // step, and total volume. + double av_prev_press = 0.0; + double av_press = 0.0; + double tot_vol = 0.0; + const int num_cells = grid_.number_of_cells; + for (int cell = 0; cell < num_cells; ++cell) { + av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell]; + av_press += state.pressure()[cell]*grid_.cell_volumes[cell]; + tot_vol += 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; + } + const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells; + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += ren_const; + } + } +#endif + + // Stop timer and report. + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + sreport.pressure_time = pt; + + // Optionally, check if well controls are satisfied. + if (check_well_controls_) { + Opm::computePhaseFlowRatesPerWell(*wells_, + well_state.perfRates(), + fractional_flows, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells_manager_.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_props_ && rock_comp_props_->isActive()) { + initial_porevol = porevol; + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } + + // Process transport sources from well flows. + Opm::computeTransportSource(props_, wells_, well_state, transport_src); + + // Solve transport. + transport_timer.start(); + double stepsize = timer.currentStepLength(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], + &initial_porevol[0], &porevol[0], &transport_src[0], stepsize, + state.saturation(), state.surfacevol()); + double substep_injected[2] = { 0.0 }; + double substep_produced[2] = { 0.0 }; + Opm::computeInjectedProduced(props_, state, transport_src, stepsize, + substep_injected, substep_produced); + injected[0] += substep_injected[0]; + injected[1] += substep_injected[1]; + produced[0] += substep_produced[0]; + produced[1] += substep_produced[1]; + if (gravity_ != 0 && use_segregation_split_) { + tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol()); + } + } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + sreport.transport_time = tt; + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol); + 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 << "\nMass balance report.\n"; + std::cout << " Injected surface volumes: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] << std::endl; + std::cout << " Produced surface volumes: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] << std::endl; + std::cout << " Total inj surface volumes: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] << std::endl; + std::cout << " Total prod surface volumes: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] << std::endl; + const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0], + init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] }; + std::cout << " Initial - inplace + inj - prod: " + << std::setw(width) << balance[0] + << std::setw(width) << balance[1] + << std::endl; + std::cout << " Relative mass error: " + << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0]) + << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1]) + << std::endl; + std::cout.precision(8); + + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, + state.pressure(), state.surfacevol(), state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } + sreport.total_time = step_timer.secsSinceStart(); + if (output_) { + sreport.reportParam(tstep_os); + } + } + + if (output_) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputWaterCut(watercut, output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } + tstep_os.close(); + } + + total_timer.stop(); + + SimulatorReport report; + report.pressure_time = ptime; + report.transport_time = ttime; + report.total_time = total_timer.secsSinceStart(); + return report; + } + + +} // namespace Opm diff --git a/opm/autodiff/SimulatorCompressibleAd.hpp b/opm/autodiff/SimulatorCompressibleAd.hpp new file mode 100644 index 000000000..98f5470f9 --- /dev/null +++ b/opm/autodiff/SimulatorCompressibleAd.hpp @@ -0,0 +1,99 @@ +/* + Copyright 2013 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 . +*/ + +#ifndef OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED +#define OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED + +#include +#include + +struct UnstructuredGrid; +struct Wells; +struct FlowBoundaryConditions; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class BlackoilPropertiesInterface; + class RockCompressibility; + class WellsManager; + class LinearSolverInterface; + class SimulatorTimer; + class BlackoilState; + class WellState; + struct SimulatorReport; + + /// Class collecting all necessary components for a two-phase simulation. + class SimulatorCompressibleAd + { + public: + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] well_manager well manager, may manage no (null) wells + /// \param[in] src source terms + /// \param[in] bcs boundary conditions, treat as all noflow if null + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorCompressibleAd(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const BlackoilPropertiesInterface& props, + const RockCompressibility* rock_comp_props, + WellsManager& wells_manager, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity); + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + /// \return simulation report, with timing data + SimulatorReport run(SimulatorTimer& timer, + BlackoilState& state, + WellState& well_state); + + private: + class Impl; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + boost::shared_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORCOMPRESSIBLEAD_HEADER_INCLUDED