Get cell volumes from Python

Adds a new method get_cell_volumes() to the opm.simulators Python module
that returns a python list of the cell volumes in the black oil
simulator.

opm/simulators/flow/python/PyBlackOilSimulator.hpp

@@ -44,6 +44,7 @@ public:
         std::shared_ptr<Opm::Schedule> schedule,
         std::shared_ptr<Opm::SummaryConfig> summary_config);
     bool checkSimulationFinished();
+    py::array_t<double> getCellVolumes();
     py::array_t<double> getPorosity();
     int run();
     void setPorosity(
@@ -56,18 +57,18 @@ public:
     const Opm::FlowMainEbos<TypeTag>& getFlowMainEbos() const;
 
 private:
-    const std::string deckFilename_;
-    bool hasRunInit_ = false;
-    bool hasRunCleanup_ = false;
+    const std::string deck_filename_;
+    bool has_run_init_ = false;
+    bool has_run_cleanup_ = false;
     //bool debug_ = false;
     // This *must* be declared before other pointers
     // to simulator objects. This in order to deinitialize
     // MPI at the correct time (ie after the other objects).
     std::unique_ptr<Opm::Main> main_;
 
-    std::unique_ptr<Opm::FlowMainEbos<TypeTag>> mainEbos_;
-    Simulator *ebosSimulator_;
-    std::unique_ptr<PyMaterialState<TypeTag>> materialState_;
+    std::unique_ptr<Opm::FlowMainEbos<TypeTag>> main_ebos_;
+    Simulator *ebos_simulator_;
+    std::unique_ptr<PyMaterialState<TypeTag>> material_state_;
     std::shared_ptr<Opm::Deck> deck_;
     std::shared_ptr<Opm::EclipseState> eclipse_state_;
     std::shared_ptr<Opm::Schedule> schedule_;

opm/simulators/flow/python/PyMaterialState.hpp

@@ -42,14 +42,14 @@ namespace Opm::Pybind
         using GridView = GetPropType<TypeTag, Opm::Properties::GridView>;
 
     public:
-        PyMaterialState(Simulator *ebosSimulator)
-            : ebosSimulator_(ebosSimulator) { }
+        PyMaterialState(Simulator *ebos_simulator)
+            : ebos_simulator_(ebos_simulator) { }
 
         std::unique_ptr<double []> getCellVolumes( std::size_t *size);
         std::unique_ptr<double []> getPorosity( std::size_t *size);
         void setPorosity(const double *poro, std::size_t size);
     private:
-        Simulator *ebosSimulator_;
+        Simulator *ebos_simulator_;
     };
 
 }

opm/simulators/flow/python/PyMaterialState_impl.hpp

@@ -17,6 +17,8 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <fmt/format.h>
+
 namespace Opm::Pybind {
 
 template <class TypeTag>
@@ -24,11 +26,11 @@ std::unique_ptr<double []>
 PyMaterialState<TypeTag>::
 getCellVolumes( std::size_t *size)
 {
-    Model &model = ebosSimulator_->model();
+    Model &model = this->ebos_simulator_->model();
     *size = model.numGridDof();
     auto array = std::make_unique<double []>(*size);
-    for (unsigned dofIdx = 0; dofIdx < *size; ++dofIdx) {
-        array[dofIdx] = model.dofTotalVolume(dofIdx);
+    for (unsigned dof_idx = 0; dof_idx < *size; ++dof_idx) {
+        array[dof_idx] = model.dofTotalVolume(dof_idx);
     }
     return array;
 }
@@ -38,12 +40,12 @@ std::unique_ptr<double []>
 PyMaterialState<TypeTag>::
 getPorosity( std::size_t *size)
 {
-    Problem &problem = ebosSimulator_->problem();
-    Model &model = ebosSimulator_->model();
+    Problem &problem = this->ebos_simulator_->problem();
+    Model &model = this->ebos_simulator_->model();
     *size = model.numGridDof();
     auto array = std::make_unique<double []>(*size);
-    for (unsigned dofIdx = 0; dofIdx < *size; ++dofIdx) {
-        array[dofIdx] = problem.referencePorosity(dofIdx, /*timeIdx*/0);
+    for (unsigned dof_idx = 0; dof_idx < *size; ++dof_idx) {
+        array[dof_idx] = problem.referencePorosity(dof_idx, /*timeIdx*/0);
     }
     return array;
 }
@@ -53,17 +55,17 @@ void
 PyMaterialState<TypeTag>::
 setPorosity(const double *poro, std::size_t size)
 {
-    Problem &problem = ebosSimulator_->problem();
-    Model &model = ebosSimulator_->model();
+    Problem &problem = this->ebos_simulator_->problem();
+    Model &model = this->ebos_simulator_->model();
     auto model_size = model.numGridDof();
     if (model_size != size) {
-        std::ostringstream message;
-        message << "Cannot set porosity. Expected array of size: "
-                << model_size << ", got array of size: " << size;
-        throw std::runtime_error(message.str());
+        const std::string msg = fmt::format(
+            "Cannot set porosity. Expected array of size: {}, got array of size: ",
+            model_size, size);
+        throw std::runtime_error(msg);
     }
-    for (unsigned dofIdx = 0; dofIdx < size; ++dofIdx) {
-        problem.setPorosity(poro[dofIdx], dofIdx);
+    for (unsigned dof_idx = 0; dof_idx < size; ++dof_idx) {
+        problem.setPorosity(poro[dof_idx], dof_idx);
     }
 }
 } //namespace Opm::Pybind

opm/simulators/flow/python/simulators.hpp

@@ -1,64 +0,0 @@
-/*
-  Copyright 2020 Equinor 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 <http://www.gnu.org/licenses/>.
-*/
-
-#ifndef OPM_SIMULATORS_HEADER_INCLUDED
-#define OPM_SIMULATORS_HEADER_INCLUDED
-
-#include <opm/simulators/flow/Main.hpp>
-#include <opm/simulators/flow/FlowMainEbos.hpp>
-#include <opm/models/utils/propertysystem.hh>
-#include <pybind11/pybind11.h>
-#include <pybind11/numpy.h>
-
-namespace py = pybind11;
-
-namespace Opm::Pybind {
-class BlackOilSimulator
-{
-private:
-    using TypeTag = Opm::Properties::TTag::EclFlowProblemTpfa;
-    using Simulator = Opm::GetPropType<TypeTag, Opm::Properties::Simulator>;
-
-public:
-    BlackOilSimulator( const std::string &deckFilename);
-    py::array_t<double> getPorosity();
-    int run();
-    void setPorosity(
-         py::array_t<double, py::array::c_style | py::array::forcecast> array);
-    int step();
-    int stepInit();
-    int stepCleanup();
-
-private:
-    const std::string deckFilename_;
-    bool hasRunInit_ = false;
-    bool hasRunCleanup_ = false;
-
-    // This *must* be declared before other pointers
-    // to simulator objects. This in order to deinitialize
-    // MPI at the correct time (ie after the other objects).
-    std::unique_ptr<Opm::Main> main_;
-
-    std::unique_ptr<Opm::FlowMainEbos<TypeTag>> mainEbos_;
-    Simulator *ebosSimulator_;
-    std::unique_ptr<PyMaterialState<TypeTag>> materialState_;
-};
-
-} // namespace Opm::Pybind
-#endif // OPM_SIMULATORS_HEADER_INCLUDED

python/simulators/PyBlackOilSimulator.cpp

@@ -24,19 +24,21 @@
 #include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 #include <opm/simulators/flow/Main.hpp>
 #include <opm/simulators/flow/FlowMainEbos.hpp>
+//#include <opm/simulators/flow/python/PyFluidState.hpp>
+#include <opm/simulators/flow/python/PyMaterialState.hpp>
+#include <opm/simulators/flow/python/PyBlackOilSimulator.hpp>
 #include <flow/flow_ebos_blackoil.hpp>
 // NOTE: EXIT_SUCCESS, EXIT_FAILURE is defined in cstdlib
 #include <cstdlib>
 #include <iostream>
 #include <stdexcept>
 #include <string>
-#include <opm/simulators/flow/python/PyBlackOilSimulator.hpp>
 
 namespace py = pybind11;
 
 namespace Opm::Pybind {
-PyBlackOilSimulator::PyBlackOilSimulator( const std::string &deckFilename)
-    : deckFilename_{deckFilename}
+PyBlackOilSimulator::PyBlackOilSimulator( const std::string &deck_filename)
+    : deck_filename_{deck_filename}
 {
 }
 
@@ -55,14 +57,14 @@ PyBlackOilSimulator::PyBlackOilSimulator(
 
 bool PyBlackOilSimulator::checkSimulationFinished()
 {
-    return this->mainEbos_->getSimTimer()->done();
+    return this->main_ebos_->getSimTimer()->done();
 }
 
 const Opm::FlowMainEbos<typename Opm::Pybind::PyBlackOilSimulator::TypeTag>&
          PyBlackOilSimulator::getFlowMainEbos() const
 {
-    if (this->mainEbos_) {
-        return *this->mainEbos_;
+    if (this->main_ebos_) {
+        return *this->main_ebos_;
     }
     else {
         throw std::runtime_error("BlackOilSimulator not initialized: "
@@ -73,14 +75,20 @@ const Opm::FlowMainEbos<typename Opm::Pybind::PyBlackOilSimulator::TypeTag>&
 py::array_t<double> PyBlackOilSimulator::getPorosity()
 {
     std::size_t len;
-    auto array = materialState_->getPorosity(&len);
+    auto array = this->material_state_->getPorosity(&len);
+    return py::array(len, array.get());
+}
+
+py::array_t<double> PyBlackOilSimulator::getCellVolumes() {
+    std::size_t len;
+    auto array = this->material_state_->getCellVolumes(&len);
     return py::array(len, array.get());
 }
 
 int PyBlackOilSimulator::run()
 {
-    auto mainObject = Opm::Main( deckFilename_ );
-    return mainObject.runStatic<Opm::Properties::TTag::EclFlowProblemTPFA>();
+    auto main_object = Opm::Main( this->deck_filename_ );
+    return main_object.runStatic<Opm::Properties::TTag::EclFlowProblemTPFA>();
 }
 
 void PyBlackOilSimulator::setPorosity( py::array_t<double,
@@ -88,7 +96,7 @@ void PyBlackOilSimulator::setPorosity( py::array_t<double,
 {
     std::size_t size_ = array.size();
     const double *poro = array.data();
-    materialState_->setPorosity(poro, size_);
+    this->material_state_->setPorosity(poro, size_);
 }
 
 void PyBlackOilSimulator::advance(int report_step)
@@ -100,10 +108,10 @@ void PyBlackOilSimulator::advance(int report_step)
 
 int PyBlackOilSimulator::step()
 {
-    if (!hasRunInit_) {
+    if (!this->has_run_init_) {
         throw std::logic_error("step() called before step_init()");
     }
-    if (hasRunCleanup_) {
+    if (this->has_run_cleanup_) {
         throw std::logic_error("step() called after step_cleanup()");
     }
     if(checkSimulationFinished()) {
@@ -111,7 +119,7 @@ int PyBlackOilSimulator::step()
     }
     //if (this->debug_)
     //    this->mainEbos_->getSimTimer()->report(std::cout);
-    auto result = mainEbos_->executeStep();
+    auto result = this->main_ebos_->executeStep();
     return result;
 }
 
@@ -119,26 +127,26 @@ int PyBlackOilSimulator::step()
 //   is called.
 int PyBlackOilSimulator::currentStep()
 {
-    return this->mainEbos_->getSimTimer()->currentStepNum();
-    // NOTE: this->ebosSimulator_->episodeIndex() would also return the current
+    return this->main_ebos_->getSimTimer()->currentStepNum();
+    // NOTE: this->ebos_simulator_->episodeIndex() would also return the current
     // report step number, but this number is always delayed by 1 step relative
-    // to this->mainEbos_->getSimTimer()->currentStepNum()
+    // to this->main_ebos_->getSimTimer()->currentStepNum()
     // See details in runStep() in file SimulatorFullyImplicitBlackoilEbos.hpp
 }
 
 
 int PyBlackOilSimulator::stepCleanup()
 {
-    hasRunCleanup_ = true;
-    return mainEbos_->executeStepsCleanup();
+    this->has_run_cleanup_ = true;
+    return this->main_ebos_->executeStepsCleanup();
 }
 
 int PyBlackOilSimulator::stepInit()
 {
 
-    if (hasRunInit_) {
+    if (this->has_run_init_) {
         // Running step_init() multiple times is not implemented yet,
-        if (hasRunCleanup_) {
+        if (this->has_run_cleanup_) {
             throw std::logic_error("step_init() called again");
         }
         else {
@@ -146,7 +154,7 @@ int PyBlackOilSimulator::stepInit()
         }
     }
     if (this->deck_) {
-        main_ = std::make_unique<Opm::Main>(
+        this->main_ = std::make_unique<Opm::Main>(
             this->deck_->getDataFile(),
             this->eclipse_state_,
             this->schedule_,
@@ -154,20 +162,19 @@ int PyBlackOilSimulator::stepInit()
         );
     }
     else {
-        main_ = std::make_unique<Opm::Main>( deckFilename_ );
+        this->main_ = std::make_unique<Opm::Main>( this->deck_filename_ );
     }
-    int exitCode = EXIT_SUCCESS;
-    mainEbos_ = main_->initFlowEbosBlackoil(exitCode);
-    if (mainEbos_) {
-        int result = mainEbos_->executeInitStep();
-        hasRunInit_ = true;
-        ebosSimulator_ = mainEbos_->getSimulatorPtr();
-        materialState_ = std::make_unique<PyMaterialState<TypeTag>>(
-            ebosSimulator_);
+    int exit_code = EXIT_SUCCESS;
+    this->main_ebos_ = this->main_->initFlowEbosBlackoil(exit_code);
+    if (this->main_ebos_) {
+        int result = this->main_ebos_->executeInitStep();
+        this->has_run_init_ = true;
+        this->ebos_simulator_ = this->main_ebos_->getSimulatorPtr();
+        this->material_state_ = std::make_unique<PyMaterialState<TypeTag>>(this->ebos_simulator_);
         return result;
     }
     else {
-        return exitCode;
+        return exit_code;
     }
 }
 
@@ -180,6 +187,8 @@ void export_PyBlackOilSimulator(py::module& m)
             std::shared_ptr<Opm::EclipseState>,
             std::shared_ptr<Opm::Schedule>,
             std::shared_ptr<Opm::SummaryConfig> >())
+        .def("get_cell_volumes", &PyBlackOilSimulator::getCellVolumes,
+            py::return_value_policy::copy)
         .def("get_porosity", &PyBlackOilSimulator::getPorosity,
             py::return_value_policy::copy)
         .def("run", &PyBlackOilSimulator::run)

python/simulators/simulators.cpp

@@ -1,125 +0,0 @@
-/*
-  Copyright 2020 Equinor 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 <http://www.gnu.org/licenses/>.
-*/
-
-#include "config.h"
-#include <opm/input/eclipse/Deck/Deck.hpp>
-#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/input/eclipse/Schedule/Schedule.hpp>
-#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
-#include <opm/simulators/flow/Main.hpp>
-#include <opm/simulators/flow/FlowMainEbos.hpp>
-#include <opm/simulators/flow/python/PyMaterialState.hpp>
-#include <pybind11/pybind11.h>
-#include <pybind11/numpy.h>
-#include <pybind11/embed.h>
-// NOTE: EXIT_SUCCESS, EXIT_FAILURE is defined in cstdlib
-#include <cstdlib>
-#include <iostream>
-#include <string>
-#include <opm/simulators/flow/python/simulators.hpp>
-
-namespace py = pybind11;
-
-namespace Opm::Pybind {
-BlackOilSimulator::BlackOilSimulator( const std::string &deckFilename)
-    : deckFilename_{deckFilename}
-{
-}
-
-py::array_t<double> BlackOilSimulator::getPorosity()
-{
-    std::size_t len;
-    auto array = materialState_->getPorosity(&len);
-    return py::array(len, array.get());
-}
-
-int BlackOilSimulator::run()
-{
-    auto mainObject = Opm::Main( deckFilename_ );
-    return mainObject.runDynamic();
-}
-
-void BlackOilSimulator::setPorosity( py::array_t<double,
-    py::array::c_style | py::array::forcecast> array)
-{
-    std::size_t size_ = array.size();
-    const double *poro = array.data();
-    materialState_->setPorosity(poro, size_);
-}
-
-int BlackOilSimulator::step()
-{
-    if (!hasRunInit_) {
-        throw std::logic_error("step() called before step_init()");
-    }
-    if (hasRunCleanup_) {
-        throw std::logic_error("step() called after step_cleanup()");
-    }
-    return mainEbos_->executeStep();
-}
-
-int BlackOilSimulator::stepCleanup()
-{
-    hasRunCleanup_ = true;
-    return mainEbos_->executeStepsCleanup();
-}
-
-int BlackOilSimulator::stepInit()
-{
-
-    if (hasRunInit_) {
-        // Running step_init() multiple times is not implemented yet,
-        if (hasRunCleanup_) {
-            throw std::logic_error("step_init() called again");
-        }
-        else {
-            return EXIT_SUCCESS;
-        }
-    }
-    main_ = std::make_unique<Opm::Main>( deckFilename_ );
-    int exitCode = EXIT_SUCCESS;
-    mainEbos_ = main_->initFlowEbosBlackoil(exitCode);
-    if (mainEbos_) {
-        int result = mainEbos_->executeInitStep();
-        hasRunInit_ = true;
-        ebosSimulator_ = mainEbos_->getSimulatorPtr();
-        materialState_ = std::make_unique<PyMaterialState<TypeTag>>(
-            ebosSimulator_);
-        return result;
-    }
-    else {
-        return exitCode;
-    }
-}
-
-} // namespace Opm::Pybind
-
-PYBIND11_MODULE(simulators, m)
-{
-    using namespace Opm::Pybind;
-    py::class_<BlackOilSimulator>(m, "BlackOilSimulator")
-        .def(py::init< const std::string& >())
-        .def("get_porosity", &BlackOilSimulator::getPorosity,
-            py::return_value_policy::copy)
-        .def("run", &BlackOilSimulator::run)
-        .def("set_porosity", &BlackOilSimulator::setPorosity)
-        .def("step", &BlackOilSimulator::step)
-        .def("step_init", &BlackOilSimulator::stepInit)
-        .def("step_cleanup", &BlackOilSimulator::stepCleanup);
-}

python/test/test_basic.py

@@ -63,7 +63,11 @@ class TestBasic(unittest.TestCase):
             sim = BlackOilSimulator("SPE1CASE1.DATA")
             sim.step_init()
             sim.step()
-
+            vol = sim.get_cell_volumes()
+            self.assertEqual(len(vol), 300, 'length of volume vector')
+            # NOTE: The volume should be 1000 ft x 1000 ft x 20 ft * 0.3 (porosity) 
+            #  = 600000 ft^3 = 566336.93 m^3
+            self.assertAlmostEqual(vol[0], 566336.93, places=2, msg='value of volume')
             poro = sim.get_porosity()
             self.assertEqual(len(poro), 300, 'length of porosity vector')
             self.assertAlmostEqual(poro[0], 0.3, places=7, msg='value of porosity')