Merge pull request #22 from equinor/add_vitrinite_reflectance

WIP: Use EasyRO from postprocessing to compute vitrinite relfectance …

warmth/mesh_model.py

@@ -12,6 +12,7 @@
 from warmth.build import single_node, Builder
 from .parameters import Parameters
 from warmth.logging import logger
+from warmth.postprocessing import VR
 from .mesh_utils import  top_crust,top_sed,thick_crust,  top_lith, top_asth, top_sed_id, bottom_sed_id,interpolateNode, interpolate_all_nodes
 from .resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_timeseries, write_hexa_grid_with_properties,read_mesh_resqml_hexa
 def tic():
@@ -336,6 +337,7 @@ def write_hexa_mesh_timeseries( self, out_path):
         age_per_vertex_keep = np.array([ self.age_per_vertex[i] for i in range(n_vertices) if i in p_to_keep ])
         Temp_per_vertex_series = np.empty([len(self.time_indices), len(p_to_keep)])
         points_cached_series = np.empty([len(self.time_indices), len(p_to_keep),3])
+        Ro_per_vertex_series = np.empty([len(self.time_indices), len(p_to_keep)])
 
         for idx, tti in enumerate(self.time_indices): # oldest first
             if idx > 0:
@@ -350,14 +352,20 @@ def write_hexa_mesh_timeseries( self, out_path):
                     point_original_to_cached[i]= count
                     count += 1
 
-
+        for i in range(Temp_per_vertex_series.shape[1]):
+            ts = Temp_per_vertex_series[:,i]
+            ro = VR.easyRoDL(ts)
+            if (i%100==0):
+                print(f"index {i} RO {ro[0:10]} {ro[-10:]}", flush=True)
+            Ro_per_vertex_series[:,i] = ro.flatten()
+
         hexa_renumbered = [ [point_original_to_cached[i] for i in hexa] for hexa in hexa_to_keep ]
 
         from os import path
 
         filename_hex = path.join(out_path, self.modelName+'_hexa_ts_'+str(self.tti)+'.epc')
         write_hexa_grid_with_timeseries(filename_hex, points_cached_series, hexa_renumbered, "hexamesh",
-            Temp_per_vertex_series,age_per_vertex_keep, poro0_per_cell, decay_per_cell, density_per_cell,
+            Temp_per_vertex_series, Ro_per_vertex_series, age_per_vertex_keep, poro0_per_cell, decay_per_cell, density_per_cell,
             cond_per_cell, rhp_per_cell, lid_per_cell )
         return filename_hex
 

warmth/resqpy_helpers.py

@@ -459,7 +459,8 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
 
 
 def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle = "hexamesh",
-    Temp_per_vertex_series=None, age_per_vertex=None, poro0_per_cell=None, decay_per_cell=None, density_per_cell=None,
+    Temp_per_vertex_series=None, Ro_per_vertex_series= None, 
+    age_per_vertex=None, poro0_per_cell=None, decay_per_cell=None, density_per_cell=None,
     cond_per_cell=None, rhp_per_cell=None, lid_per_cell=None ):
     """Writes the given hexahedral mesh, defined by arrays of nodes and cell indices, into a RESQML .epc file
        Given SubsHeat properties are optionally written.
@@ -583,6 +584,19 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
                                                 time_index = time_index,
                                                 indexable_element = 'nodes')
                                                 # points = True)
+        if (Ro_per_vertex_series is not None):
+            ro = Ro_per_vertex_series[time_index,:].astype(np.float32)
+            pc.add_cached_array_to_imported_list(ro,
+                                                    'Vitrinite reflectance',
+                                                    "%Ro",
+                                                    uom = 'percent',
+                                                    property_kind = 'dimensionless',
+                                                    realization = 0,
+                                                    time_index = time_index,
+                                                    indexable_element = 'nodes')
+                                                    # points = True)
+
+
     pc.write_hdf5_for_imported_list()
     pc.create_xml_for_imported_list_and_add_parts_to_model(time_series_uuid = gts.uuid)