#12568 Python: add API for sampling given positions in a grid case

ApplicationLibCode/ProjectDataModelCommands/RimcEclipseCase.cpp

@@ -20,23 +20,25 @@
 
 #include "RiaApplication.h"
 #include "RiaGuiApplication.h"
+#include "RiaKeyValueStoreUtil.h"
 
-#include "RicImportSummaryCasesFeature.h"
+#include "RigCaseCellResultsData.h"
+#include "RigEclipseCaseData.h"
+#include "RigEclipseResultAddress.h"
+#include "RigMainGrid.h"
+#include "RigResultAccessor.h"
+#include "RigResultAccessorFactory.h"
 
 #include "RimEclipseCase.h"
-#include "RimFileSummaryCase.h"
-#include "RimOilField.h"
 #include "RimProject.h"
-#include "RimSummaryCase.h"
-#include "RimSurfaceCollection.h"
-#include "RiuPlotMainWindow.h"
 
 #include "cafPdmFieldScriptingCapability.h"
 
 #include <QDir>
 #include <QFileInfo>
 
-#include <memory>
+#include <expected>
+#include <limits>
 
 CAF_PDM_OBJECT_METHOD_SOURCE_INIT( RimEclipseCase, RimcEclipseCase_importProperties, "import_properties" );
 
@@ -75,3 +77,107 @@ std::expected<caf::PdmObjectHandle*, QString> RimcEclipseCase_importProperties::
 
     return nullptr;
 }
+
+CAF_PDM_OBJECT_METHOD_SOURCE_INIT( RimEclipseCase, RimcEclipseCase_exportValuesInternal, "export_values_internal" );
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+RimcEclipseCase_exportValuesInternal::RimcEclipseCase_exportValuesInternal( caf::PdmObjectHandle* self )
+    : caf::PdmVoidObjectMethod( self )
+{
+    CAF_PDM_InitObject( "Export Values Internal", "", "", "Export Values Internal" );
+
+    CAF_PDM_InitScriptableFieldNoDefault( &m_coordinateX, "CoordinateX", "" );
+    CAF_PDM_InitScriptableFieldNoDefault( &m_coordinateY, "CoordinateY", "" );
+    CAF_PDM_InitScriptableFieldNoDefault( &m_coordinateZ, "CoordinateZ", "" );
+    CAF_PDM_InitScriptableFieldNoDefault( &m_resultKey, "ResultKey", "" );
+
+    CAF_PDM_InitScriptableFieldNoDefault( &m_propertyType, "PropertyType", "" );
+    CAF_PDM_InitScriptableFieldNoDefault( &m_propertyName, "PropertyName", "" );
+    CAF_PDM_InitScriptableFieldNoDefault( &m_timeStep, "TimeStep", "" );
+    CAF_PDM_InitScriptableFieldNoDefault( &m_porosityModel, "PorosityModel", "" );
+}
+
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::expected<caf::PdmObjectHandle*, QString> RimcEclipseCase_exportValuesInternal::execute()
+{
+    auto eclipseCase = self<RimEclipseCase>();
+
+    RigEclipseCaseData* eclipseCaseData = eclipseCase->eclipseCaseData();
+    if ( !eclipseCaseData )
+    {
+        return std::unexpected( "Unable to load eclipse case data." );
+    }
+
+    const RigMainGrid* mainGrid = eclipseCase->mainGrid();
+    if ( !mainGrid )
+    {
+        return std::unexpected( "Unable to load main grid for eclipse case." );
+    }
+
+    caf::AppEnum<RiaDefines::ResultCatType> resultCatTypeEnum;
+    if ( !resultCatTypeEnum.setFromText( m_propertyType ) )
+    {
+        return std::unexpected( "Invalid property type." );
+    }
+
+    caf::AppEnum<RiaDefines::PorosityModelType> porosityModel;
+    if ( !porosityModel.setFromText( m_porosityModel ) )
+    {
+        return std::unexpected( "Invalid porosity model." );
+    }
+
+    RigEclipseResultAddress resVarAddr( resultCatTypeEnum, m_propertyName );
+    eclipseCase->results( porosityModel )->ensureKnownResultLoaded( resVarAddr );
+
+    cvf::ref<RigResultAccessor> resultAccessor =
+        RigResultAccessorFactory::createFromResultAddress( eclipseCaseData, 0, porosityModel, m_timeStep, resVarAddr );
+    if ( resultAccessor.isNull() )
+    {
+        return std::unexpected( "Result property not found." );
+    }
+
+    if ( m_timeStep < 0 || m_timeStep >= static_cast<int>( eclipseCase->timeStepDates().size() ) )
+    {
+        return std::unexpected( "Invalid time step." );
+    }
+
+    auto keyValueStore = RiaApplication::instance()->keyValueStore();
+
+    std::vector<float> xs = RiaKeyValueStoreUtil::convertToFloatVector( keyValueStore->get( m_coordinateX().toStdString() ) );
+    std::vector<float> ys = RiaKeyValueStoreUtil::convertToFloatVector( keyValueStore->get( m_coordinateY().toStdString() ) );
+    std::vector<float> zs = RiaKeyValueStoreUtil::convertToFloatVector( keyValueStore->get( m_coordinateZ().toStdString() ) );
+    if ( xs.empty() || xs.size() != ys.size() || xs.size() != zs.size() )
+    {
+        return std::unexpected( "Invalid positions specified." );
+    }
+
+    std::vector<cvf::Vec3d> positions;
+    for ( size_t i = 0; i < xs.size(); i++ )
+    {
+        positions.push_back( cvf::Vec3d( xs[i], ys[i], -zs[i] ) );
+    }
+
+    std::vector<float> values;
+    for ( const cvf::Vec3d& position : positions )
+    {
+        auto cellIdx = mainGrid->findReservoirCellIndexFromPoint( position );
+        if ( cellIdx != cvf::UNDEFINED_SIZE_T )
+        {
+            float valueFromEclipse = static_cast<float>( resultAccessor->cellScalar( cellIdx ) );
+            values.push_back( valueFromEclipse );
+        }
+        else
+        {
+            // Use inf as signal that nothing was found for the position
+            values.push_back( std::numeric_limits<float>::infinity() );
+        }
+    }
+
+    keyValueStore->set( m_resultKey().toStdString(), RiaKeyValueStoreUtil::convertToByteVector( values ) );
+
+    return nullptr;
+}

ApplicationLibCode/ProjectDataModelCommands/RimcEclipseCase.h

@@ -41,3 +41,26 @@ public:
 private:
     caf::PdmField<std::vector<QString>> m_fileNames;
 };
+
+//==================================================================================================
+///
+//==================================================================================================
+class RimcEclipseCase_exportValuesInternal : public caf::PdmVoidObjectMethod
+{
+    CAF_PDM_HEADER_INIT;
+
+public:
+    RimcEclipseCase_exportValuesInternal( caf::PdmObjectHandle* self );
+
+    std::expected<caf::PdmObjectHandle*, QString> execute() override;
+
+private:
+    caf::PdmField<QString> m_coordinateX;
+    caf::PdmField<QString> m_coordinateY;
+    caf::PdmField<QString> m_coordinateZ;
+    caf::PdmField<QString> m_propertyType;
+    caf::PdmField<QString> m_propertyName;
+    caf::PdmField<int>     m_timeStep;
+    caf::PdmField<QString> m_porosityModel;
+    caf::PdmField<QString> m_resultKey;
+};

GrpcInterface/Python/rips/case.py

@@ -36,6 +36,7 @@ result
 """
 
 import grpc
+import uuid
 from typing import List, Tuple
 
 import Case_pb2
@@ -56,6 +57,7 @@ from .resinsight_classes import (
 )
 
 from .grid import Grid as Grid
+from .project import Project as Project
 from .pdmobject import add_method
 from .view import View as View
 from .simulation_well import SimulationWell
@@ -1332,3 +1334,55 @@ def set_nnc_connections_values(
     reply = self.__nnc_properties_stub.SetNNCValues(request_iterator)
     if reply.accepted_value_count < len(values):
         raise IndexError
+
+
+@add_method(EclipseCase)
+def grid_property_for_positions(
+    self,
+    positions: List[List[float]],
+    property_type: str,
+    property_name: str,
+    time_step: int,
+    porosity_model: str = "MATRIX_MODEL",
+) -> List[float]:
+    shared_uuid = uuid.uuid4()
+    coordinate_x = "{}_{}".format(shared_uuid, "coordinate_x")
+    coordinate_y = "{}_{}".format(shared_uuid, "coordinate_y")
+    coordinate_z = "{}_{}".format(shared_uuid, "coordinate_z")
+    result_key = "{}_{}".format(shared_uuid, "result")
+
+    coord_x = []
+    coord_y = []
+    coord_z = []
+    for pos in positions:
+        coord_x.append(pos[0])
+        coord_y.append(pos[1])
+        coord_z.append(pos[2])
+
+        # Get the results from the key-value store.
+    project = self.ancestor(Project)
+    if project:
+        project.set_key_values(coordinate_x, coord_x)
+        project.set_key_values(coordinate_y, coord_y)
+        project.set_key_values(coordinate_z, coord_z)
+
+        self.export_values_internal(
+            coordinate_x=coordinate_x,
+            coordinate_y=coordinate_y,
+            coordinate_z=coordinate_z,
+            result_key=result_key,
+            property_type=property_type,
+            property_name=property_name,
+            time_step=time_step,
+            porosity_model=porosity_model,
+        )
+        result = project.key_values(result_key)
+
+        project.remove_key_values(coordinate_x)
+        project.remove_key_values(coordinate_y)
+        project.remove_key_values(coordinate_z)
+        project.remove_key_values(result_key)
+
+        return result
+    else:
+        return []

GrpcInterface/Python/rips/tests/test_cases.py

@@ -216,3 +216,133 @@ def test_multiple_load_of_same_case(rips_instance, initialize_test):
     for i in range(case_count):
         c = rips_instance.project.load_case(path_name)
         assert c
+
+
+def test_10k_property_for_positions(rips_instance, initialize_test):
+    case_path = dataroot.PATH + "/TEST10K_FLT_LGR_NNC/TEST10K_FLT_LGR_NNC.EGRID"
+    case = rips_instance.project.load_case(path=case_path)
+    assert len(case.grids()) == 2
+    time_steps = case.time_steps()
+    assert len(time_steps) == 9
+
+    # Expected values from ResInsight UI (cell ijk=[23, 44, 19])
+    positions = [
+        [3655.67, 5145.34, 4176.63],
+        [3690.07, 5240.69, 4180.02],
+        [3599.87, 5275.16, 4179.32],
+        [3654.78, 5144.79, 4179.23],
+        [3688.99, 5239.88, 4182.7],
+        [3598.62, 5274.48, 4181.96],
+    ]
+
+    property_type = "STATIC_NATIVE"
+    property_name = "DX"
+    porosity_model = "MATRIX_MODEL"
+    time_step = 0
+    result = case.grid_property_for_positions(
+        positions, property_type, property_name, time_step, porosity_model
+    )
+    assert result
+    assert len(result) == 6
+
+    for dx in result:
+        assert dx > 96 and dx < 99
+
+
+def test_10k_property_for_positions_outside_of_model(rips_instance, initialize_test):
+    case_path = dataroot.PATH + "/TEST10K_FLT_LGR_NNC/TEST10K_FLT_LGR_NNC.EGRID"
+    case = rips_instance.project.load_case(path=case_path)
+
+    positions = [
+        [3655.67, 5145.34, -4176.63],
+        [3690.07, 52400.69, 4180.02],
+    ]
+
+    property_type = "STATIC_NATIVE"
+    property_name = "DX"
+    porosity_model = "MATRIX_MODEL"
+    time_step = 0
+    result = case.grid_property_for_positions(
+        positions, property_type, property_name, time_step, porosity_model
+    )
+
+    assert result
+    assert len(result) == 2
+    for dx in result:
+        assert math.isinf(dx)
+
+
+def test_10k_property_for_positions_no_positions(rips_instance, initialize_test):
+    case_path = dataroot.PATH + "/TEST10K_FLT_LGR_NNC/TEST10K_FLT_LGR_NNC.EGRID"
+    case = rips_instance.project.load_case(path=case_path)
+
+    positions = []
+
+    property_type = "STATIC_NATIVE"
+    property_name = "DX"
+    porosity_model = "MATRIX_MODEL"
+    time_step = 0
+    with pytest.raises(rips.RipsError, match="Invalid positions specified"):
+        case.grid_property_for_positions(
+            positions, property_type, property_name, time_step, porosity_model
+        )
+
+
+def test_10k_property_for_positions_invalid_property_name(
+    rips_instance, initialize_test
+):
+    case_path = dataroot.PATH + "/TEST10K_FLT_LGR_NNC/TEST10K_FLT_LGR_NNC.EGRID"
+    case = rips_instance.project.load_case(path=case_path)
+
+    positions = [
+        [3655.67, 5145.34, 4176.63],
+    ]
+
+    property_type = "STATIC_NATIVE"
+    invalid_property_name = "NON_EXISTING_RESULT"
+    porosity_model = "MATRIX_MODEL"
+    time_step = 0
+    with pytest.raises(rips.RipsError, match="Result property not found."):
+        case.grid_property_for_positions(
+            positions, property_type, invalid_property_name, time_step, porosity_model
+        )
+
+
+def test_10k_property_for_positions_invalid_property_type(
+    rips_instance, initialize_test
+):
+    case_path = dataroot.PATH + "/TEST10K_FLT_LGR_NNC/TEST10K_FLT_LGR_NNC.EGRID"
+    case = rips_instance.project.load_case(path=case_path)
+
+    positions = [
+        [3655.67, 5145.34, 4176.63],
+    ]
+
+    property_type = "NON_EXISTING_PROPERTY_TYPE"
+    invalid_property_name = "PRESSURE"
+    porosity_model = "MATRIX_MODEL"
+    time_step = 0
+    with pytest.raises(rips.RipsError, match="Invalid property type."):
+        case.grid_property_for_positions(
+            positions, property_type, invalid_property_name, time_step, porosity_model
+        )
+
+
+def test_10k_property_for_positions_invalid_time_step(rips_instance, initialize_test):
+    case_path = dataroot.PATH + "/TEST10K_FLT_LGR_NNC/TEST10K_FLT_LGR_NNC.EGRID"
+    case = rips_instance.project.load_case(path=case_path)
+    time_steps = case.time_steps()
+    assert len(time_steps) == 9
+
+    positions = [
+        [3655.67, 5145.34, 4176.63],
+    ]
+
+    property_type = "DYNAMIC_NATIVE"
+    property_name = "PRESSURE"
+    porosity_model = "MATRIX_MODEL"
+    invalid_time_step = 99
+    with pytest.raises(rips.RipsError, match="Invalid time step."):
+        case.grid_property_for_positions(
+            positions, property_type, property_name, invalid_time_step, porosity_model
+        )

GrpcInterface/Python/rips/well_path.py

@@ -48,12 +48,15 @@ def trajectory_properties(
 
     # Get the results from the key-value store.
     project = self.ancestor(Project)
-    result = {
-        prop: project.key_values(temp_key) for prop, temp_key in temp_keys.items()
-    }
+    if project:
+        result = {
+            prop: project.key_values(temp_key) for prop, temp_key in temp_keys.items()
+        }
 
-    # Delete results form key-value store.
-    for temp_key in temp_keys.values():
-        project.remove_key_values(temp_key)
+        # Delete results form key-value store.
+        for temp_key in temp_keys.values():
+            project.remove_key_values(temp_key)
 
-    return result
+        return result
+    else:
+        return {}