diff --git a/Makefile.am b/Makefile.am index 1ffe69a8..f5be231e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -94,7 +94,7 @@ opm/core/pressure/tpfa/trans_tpfa.c \ opm/core/pressure/well.c \ opm/core/simulator/SimulatorReport.cpp \ opm/core/simulator/SimulatorTimer.cpp \ -opm/core/simulator/SimulatorTwophase.cpp \ +opm/core/simulator/SimulatorIncompTwophase.cpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \ opm/core/transport/reorder/TransportModelInterface.cpp \ opm/core/transport/reorder/TransportModelTwophase.cpp \ @@ -197,7 +197,7 @@ opm/core/pressure/tpfa/trans_tpfa.h \ opm/core/simulator/BlackoilState.hpp \ opm/core/simulator/SimulatorReport.hpp \ opm/core/simulator/SimulatorTimer.hpp \ -opm/core/simulator/SimulatorTwophase.hpp \ +opm/core/simulator/SimulatorIncompTwophase.hpp \ opm/core/simulator/TwophaseState.hpp \ opm/core/simulator/WellState.hpp \ opm/core/transport/CSRMatrixBlockAssembler.hpp \ diff --git a/examples/sim_2p_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp index eba5dbb7..109e5a94 100644 --- a/examples/sim_2p_incomp_reorder.cpp +++ b/examples/sim_2p_incomp_reorder.cpp @@ -43,7 +43,7 @@ #include #include -#include +#include #include #include @@ -180,7 +180,7 @@ main(int argc, char** argv) catch (...) { THROW("Creating directories failed: " << fpath); } - param.writeParam(output_dir + "/spu_2p.param"); + param.writeParam(output_dir + "/simulation.param"); } @@ -191,15 +191,15 @@ main(int argc, char** argv) SimulatorReport rep; if (!use_deck) { // Simple simulation without a deck. - SimulatorTwophase simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - 0, // wells - src, - bcs.c_bcs(), - linsolver, - grav); + SimulatorIncompTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + 0, // wells + src, + bcs.c_bcs(), + linsolver, + grav); SimulatorTimer simtimer; simtimer.init(param); warnIfUnusedParams(param); @@ -246,15 +246,15 @@ main(int argc, char** argv) } // Create and run simulator. - SimulatorTwophase simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - wells.c_wells(), - src, - bcs.c_bcs(), - linsolver, - grav); + SimulatorIncompTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells.c_wells(), + src, + bcs.c_bcs(), + linsolver, + grav); if (epoch == 0) { warnIfUnusedParams(param); } diff --git a/opm/core/simulator/SimulatorTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp similarity index 82% rename from opm/core/simulator/SimulatorTwophase.cpp rename to opm/core/simulator/SimulatorIncompTwophase.cpp index 489401f0..7a369bc7 100644 --- a/opm/core/simulator/SimulatorTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -21,7 +21,7 @@ #include "config.h" #endif // HAVE_CONFIG_H -#include +#include #include #include @@ -56,7 +56,7 @@ namespace Opm { - class SimulatorTwophase::Impl + class SimulatorIncompTwophase::Impl { public: Impl(const parameter::ParameterGroup& param, @@ -78,6 +78,7 @@ namespace Opm // Parameters for output. bool output_; + bool output_vtk_; std::string output_dir_; int output_interval_; // Parameters for transport solver. @@ -104,15 +105,15 @@ namespace Opm - SimulatorTwophase::SimulatorTwophase(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const RockCompressibility* rock_comp, - const Wells* wells, - const std::vector& src, - const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, - const double* gravity) + SimulatorIncompTwophase::SimulatorIncompTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) { pimpl_.reset(new Impl(param, grid, props, rock_comp, wells, src, bcs, linsolver, gravity)); } @@ -121,19 +122,19 @@ namespace Opm - SimulatorReport SimulatorTwophase::run(SimulatorTimer& timer, - TwophaseState& state, - WellState& well_state) + SimulatorReport SimulatorIncompTwophase::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) { return pimpl_->run(timer, state, well_state); } - static void outputState(const UnstructuredGrid& grid, - const Opm::TwophaseState& state, - const int step, - const std::string& output_dir) + static void outputStateVtk(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) { // Write data in VTK format. std::ostringstream vtkfilename; @@ -157,12 +158,26 @@ namespace Opm Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); dm["velocity"] = &cell_velocity; Opm::writeVtkData(grid, dm, vtkfile); + } + + + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) + { + 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; // 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; - fpath = fname.str(); + boost::filesystem::path fpath = fname.str(); try { create_directories(fpath); } @@ -209,15 +224,15 @@ namespace Opm - SimulatorTwophase::Impl::Impl(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const RockCompressibility* rock_comp, - const Wells* wells, - const std::vector& src, - const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, - const double* gravity) + SimulatorIncompTwophase::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + LinearSolverInterface& linsolver, + const double* gravity) : grid_(grid), props_(props), rock_comp_(rock_comp), @@ -238,6 +253,7 @@ namespace Opm // 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_); @@ -269,9 +285,9 @@ namespace Opm - SimulatorReport SimulatorTwophase::Impl::run(SimulatorTimer& timer, - TwophaseState& state, - WellState& well_state) + SimulatorReport SimulatorIncompTwophase::Impl::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) { std::vector transport_src; @@ -314,7 +330,10 @@ namespace Opm // Report timestep and (optionally) write state to disk. timer.report(std::cout); if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - outputState(grid_, state, timer.currentStepNum(), output_dir_); + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); } // Solve pressure. @@ -400,7 +419,10 @@ namespace Opm } if (output_) { - outputState(grid_, state, timer.currentStepNum(), output_dir_); + 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_); diff --git a/opm/core/simulator/SimulatorTwophase.hpp b/opm/core/simulator/SimulatorIncompTwophase.hpp similarity index 84% rename from opm/core/simulator/SimulatorTwophase.hpp rename to opm/core/simulator/SimulatorIncompTwophase.hpp index 818d40db..e5e0a966 100644 --- a/opm/core/simulator/SimulatorTwophase.hpp +++ b/opm/core/simulator/SimulatorIncompTwophase.hpp @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef OPM_SIMULATORTWOPHASE_HEADER_INCLUDED -#define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED +#ifndef OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED +#define OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED #include #include @@ -39,7 +39,7 @@ namespace Opm struct SimulatorReport; /// Class collecting all necessary components for a two-phase simulation. - class SimulatorTwophase + class SimulatorIncompTwophase { public: /// Initialise from parameters and objects to observe. @@ -66,15 +66,15 @@ namespace Opm /// \param[in] bcs boundary conditions, treat as all noflow if null /// \param[in] linsolver linear solver /// \param[in] gravity if non-null, gravity vector - SimulatorTwophase(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const IncompPropertiesInterface& props, - const RockCompressibility* rock_comp, - const Wells* wells, - const std::vector& src, - const FlowBoundaryConditions* bcs, - LinearSolverInterface& linsolver, - const double* gravity); + SimulatorIncompTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + 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