diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 343a611db..74d266709 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -31,6 +31,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/FullyImplicitBlackoilSolver.cpp
opm/autodiff/ImpesTPFAAD.cpp
opm/autodiff/SimulatorCompressibleAd.cpp
+ opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
opm/autodiff/SimulatorIncompTwophaseAdfi.cpp
opm/autodiff/TransportSolverTwophaseAd.cpp
)
@@ -76,6 +77,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/ImpesTPFAAD.hpp
opm/autodiff/FullyImplicitBlackoilSolver.hpp
opm/autodiff/SimulatorCompressibleAd.hpp
+ opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
opm/autodiff/SimulatorIncompTwophaseAdfi.hpp
opm/autodiff/TransportSolverTwophaseAd.hpp
)
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
new file mode 100644
index 000000000..01db9f573
--- /dev/null
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp
@@ -0,0 +1,480 @@
+/*
+ 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 SimulatorFullyImplicitBlackoil::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_;
+ FullyImplicitBlackoilSolver solver_;
+ // Misc. data
+ std::vector allcells_;
+ };
+
+
+
+
+ SimulatorFullyImplicitBlackoil::SimulatorFullyImplicitBlackoil(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 SimulatorFullyImplicitBlackoil::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"));
+ }
+ }
+
+
+#if 0
+ 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);
+ }
+#endif
+
+
+ // \TODO: make CompressibleTpfa take src and bcs.
+ SimulatorFullyImplicitBlackoil::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_),
+ solver_(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, */
+ {
+ // 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);
+
+ // 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 SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
+ BlackoilState& state,
+ WellState& well_state)
+ {
+ // 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 solver_timer;
+ double stime = 0.0;
+ Opm::time::StopWatch step_timer;
+ Opm::time::StopWatch total_timer;
+ total_timer.start();
+#if 0
+ // These must be changed for three-phase.
+ 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;
+#endif
+ 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);
+#if 0
+ wellreport.push(props_, *wells_,
+ state.pressure(), state.surfacevol(), state.saturation(),
+ 0.0, well_state.bhp(), well_state.perfRates());
+#endif
+ }
+ 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.
+ solver_timer.start();
+ std::vector initial_pressure = state.pressure();
+ solver_.step(timer.currentStepLength(), state/*, well_state*/);
+
+ // Stop timer and report.
+ solver_timer.stop();
+ const double st = solver_timer.secsSinceStart();
+ std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
+ stime += st;
+ sreport.pressure_time = st;
+
+ // 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);
+ }
+
+ // The reports below are geared towards two phases only.
+#if 0
+ // Report mass balances.
+ double injected[2] = { 0.0 };
+ double produced[2] = { 0.0 };
+ Opm::computeInjectedProduced(props_, state, transport_src, stepsize,
+ injected, produced);
+ 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);
+
+ // Make well reports.
+ 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());
+ }
+#endif
+ 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_);
+#if 0
+ outputWaterCut(watercut, output_dir_);
+ if (wells_) {
+ outputWellReport(wellreport, output_dir_);
+ }
+#endif
+ tstep_os.close();
+ }
+
+ total_timer.stop();
+
+ SimulatorReport report;
+ report.pressure_time = stime;
+ report.transport_time = 0.0;
+ report.total_time = total_timer.secsSinceStart();
+ return report;
+ }
+
+
+} // namespace Opm
diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp
new file mode 100644
index 000000000..5ef8cd7b7
--- /dev/null
+++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.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_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_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 SimulatorFullyImplicitBlackoil
+ {
+ 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
+ SimulatorFullyImplicitBlackoil(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_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED