Merge pull request #19 from equinor/fix_node_indexing_and_resqml_types

Fix node indexing and resqml types

tests/warmth3d/test_3d.py

@@ -11,7 +11,7 @@
 @pytest.mark.mpi
 def test_3d_compare():
     comm = MPI.COMM_WORLD
-    inc = 1000
+    inc = 2100
     model_pickled = f"model-out-inc_{inc}.p"
     if comm.rank == 0 and not os.path.isfile(model_pickled):
         global runtime_1D_sim
@@ -61,7 +61,7 @@ def test_3d_compare():
         print("Total time 1D simulations:", runtime_1D_sim)
 
         pickle.dump( model, open( model_pickled, "wb" ) )
-        model = pickle.load( open( model_pickled, "rb" ) )
+        # model = pickle.load( open( model_pickled, "rb" ) )
         try:
             os.mkdir('mesh')
         except FileExistsError:

warmth/build.py

@@ -720,11 +720,10 @@ def _create_nodes(self, all_sediments_grid: List[List[List]]):
                                                 'rhp': rhp, 'phi': phi, 'decay': decay, 'solidus': solidus, 'liquidus': liquidus,'strat':strat,'horizonIndex':inputRef})
                     df = df.sort_values(by=["topage"],ignore_index=True)
 
-
+                    # print(df)
                     #df.reset_index(drop=True,inplace=True)
-                    df.at[2, 'top'] = np.nan
-
-                    df.at[3, 'top'] = np.nan
+                    # df.at[2, 'top'] = np.nan            
+                    # df.at[3, 'top'] = np.nan
                     checker = self._check_nan_sed(df)
 
                     if checker is False:

warmth/forward_modelling.py

@@ -282,7 +282,8 @@ def _remesh_crust_lith_asth(coord: np.ndarray[np.float64], key_depths_arr: np.nd
 
                 lower_idx = np.argwhere(coord_new > key_depth)
                 if lower_idx.size == 0:
-                    raise Exception
+                    # raise Exception
+                    lower_idx = len(coord_new)-1
                 else:
                     lower_idx = lower_idx[0][0]
 

warmth/resqpy_helpers.py

@@ -10,8 +10,6 @@
 import resqpy.time_series as rts
 
 #
-#  Our example resqml model can be read using the read_mesh_resqml function below.. 
-#  read_mesh_resqml("/path/mapA-961-nodes-182_0.epc")
 # 
 #
 
@@ -30,12 +28,10 @@ def read_mesh_resqml(epcfilename, meshTitle = 'tetramesh'):
     assert tetra_uuid is not None
     tetra = rug.TetraGrid(model, uuid = tetra_uuid)
     assert tetra is not None
-    print(tetra.title, tetra.node_count, tetra.cell_count, tetra.cell_shape)
+    print(f"Mesh {tetra.title}: {tetra.node_count} nodes, {tetra.cell_count} cells, {tetra.cell_shape} cell shape")
     assert tetra.cell_shape == 'tetrahedral'
 
-    print( tetra.points_ref().shape )   # numpy array of vertex positions
     cells = np.array( [ tetra.distinct_node_indices_for_cell(i) for i in range(tetra.cell_count) ]  ) # cell indices are read using this function(?)
-    print( cells.shape )   # numpy array of vertex positions
 
     tetra.check_tetra()
 
@@ -48,21 +44,20 @@ def read_mesh_resqml(epcfilename, meshTitle = 'tetramesh'):
     temp_prop = rqp.Property(model, uuid = temp_uuid)
     assert temp_prop.uom() == 'degC'
     assert temp_prop.indexable_element() == 'nodes'   # properties are defined either on nodes or on cells
-    print( temp_prop.array_ref().shape, temp_prop.array_ref()[0:10] )  # .array_ref() exposes the values as numpy array
 
     layerID_uuid = model.uuid(title = 'LayerID')
     assert layerID_uuid is not None
     layerID_prop = rqp.Property(model, uuid = layerID_uuid)
     # assert layerID_prop.uom() == 'Euc'
     assert layerID_prop.is_continuous() == False
     assert layerID_prop.indexable_element() == 'cells'
-    print( layerID_prop.array_ref().shape, layerID_prop.array_ref()[0:10] )  # .array_ref() exposes the values as numpy array
+    # print( layerID_prop.array_ref().shape, layerID_prop.array_ref()[0:10] )  # .array_ref() exposes the values as numpy array
 
     titles=['Temperature', 'Age', 'LayerID', 'Porosity_initial', 'Porosity_decay', 'Density_solid', 'insulance_thermal', 'Radiogenic_heat_production']
     for title in titles:
         prop_uuid = model.uuid(title = title)
         prop = rqp.Property(model, uuid = prop_uuid)
-        print(title, prop.indexable_element(), prop.uom(), prop.array_ref()[0:10] )
+        print(f"Property {title}: defined on {prop.indexable_element()}, unit {prop.uom()}, first values: {prop.array_ref()[0:10]}")
 
 def write_tetra_grid_with_properties(filename, nodes, cells, modelTitle = "tetramesh",
     Temp_per_vertex=None, age_per_vertex=None, poro0_per_cell=None, decay_per_cell=None, density_per_cell=None,
@@ -382,7 +377,7 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
 
     if Temp_per_vertex is not None:
         _ = rqp.Property.from_array(model,
-                                    Temp_per_vertex,
+                                    Temp_per_vertex.astype(np.float32),
                                     source_info = 'SubsHeat',
                                     keyword = 'Temperature',
                                     support_uuid = hexa.uuid,
@@ -392,7 +387,7 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
 
     if age_per_vertex is not None:
         _ = rqp.Property.from_array(model,
-                                    age_per_vertex,
+                                    age_per_vertex.astype(np.float32),
                                     source_info = 'SubsHeat',
                                     keyword = 'Age',
                                     support_uuid = hexa.uuid,
@@ -413,7 +408,7 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
 
     if poro0_per_cell is not None:
         _ = rqp.Property.from_array(model,
-                                    poro0_per_cell,
+                                    poro0_per_cell.astype(np.float32),
                                     source_info = 'SubsHeat',
                                     keyword = 'Porosity_initial',
                                     support_uuid = hexa.uuid,
@@ -422,7 +417,7 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
                                     uom = 'm3/m3')
     if decay_per_cell is not None:
         _ = rqp.Property.from_array(model,
-                                    decay_per_cell,
+                                    decay_per_cell.astype(np.float32),
                                     source_info = 'SubsHeat',
                                     keyword = 'Porosity_decay',
                                     support_uuid = hexa.uuid,
@@ -431,7 +426,7 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
                                     uom = 'Euc')
     if density_per_cell is not None:
         _ = rqp.Property.from_array(model,
-                                    density_per_cell,
+                                    density_per_cell.astype(np.float32),
                                     source_info = 'SubsHeat',
                                     keyword = 'Density_solid',
                                     support_uuid = hexa.uuid,
@@ -443,7 +438,7 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
         # we write thermal conductivity as its inverse, the thermal insulance
         #
         _ = rqp.Property.from_array(model,
-                                    np.reciprocal(cond_per_cell),
+                                    np.reciprocal(cond_per_cell).astype(np.float32),
                                     source_info = 'SubsHeat',
                                     keyword = 'insulance_thermal',
                                     support_uuid = hexa.uuid,
@@ -452,7 +447,7 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
                                     uom = 'deltaK.m2/W')
     if rhp_per_cell is not None:
         _ = rqp.Property.from_array(model,
-                                    rhp_per_cell,
+                                    rhp_per_cell.astype(np.float32),
                                     source_info = 'SubsHeat',
                                     keyword = 'Radiogenic_heat_production',
                                     support_uuid = hexa.uuid,
@@ -563,7 +558,7 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
     # nodes0 = nodes.copy()
     for time_index in range(len(nodes_series)-1,-1,-1):
         # nodes2 = nodes0 + [0,0,time_index*10]
-        nodes2 = nodes_series[time_index]
+        nodes2 = nodes_series[time_index].astype(np.float32)
         pc.add_cached_array_to_imported_list(nodes2,
                                                 'dynamic nodes',
                                                 'points',
@@ -574,7 +569,7 @@ def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle =
                                                 indexable_element = 'nodes',
                                                 points = True)
         # active_array = np.ones([2160], dtype = bool)
-        tt = Temp_per_vertex_series[time_index]
+        tt = Temp_per_vertex_series[time_index].astype(np.float32)
         pc.add_cached_array_to_imported_list(tt,
                                                 'Temperature',
                                                 'Temperature',

warmth3d/fixed_mesh_model.py

@@ -421,13 +421,11 @@ def rhpForLayerID(self, ss, node_index):
         return 0
 
 
-    def buildVertices(self, time_index=0, useFakeEncodedZ=False):
+    def buildVertices(self, time_index=0):
         """Determine vertex positions, node-by-node.
            For every node, the same number of vertices is added (one per sediment, one per crust, lith, asth, and one at the bottom)
            Degenerate vertices (e.g. at nodes where sediment is yet to be deposited or has been eroded) are avoided by a small shift, kept in self.sed_diff_z
            
-           When the option useFakeEncodedZ is set, the z-values are repurposed to encode the index of the vertex in the original, deterministic indexing.
-           This is necessary because dolfinx will re-index the vertices upon mesh generation. 
         """           
         tti = time_index
         self.tti = time_index
@@ -487,11 +485,7 @@ def buildVertices(self, time_index=0, useFakeEncodedZ=False):
         self.sed_diff_z = np.array(self.sed_diff_z)
         self.mesh_vertices = self.mesh_vertices_0.copy()
         self.mesh_vertices[:,2] = self.mesh_vertices_0[:,2] + self.sed_diff_z
-        if (useFakeEncodedZ):
-            self.mesh_vertices[:,2] = np.ceil(self.mesh_vertices[:,2])*1000 + np.array(list(range(self.mesh_vertices.shape[0])))*0.01
-            # zz = self.mesh_vertices[:,2].copy()
-            # zz2=np.mod(zz,1000)
-            # # mesh_reindex = (1e-4+zz2*10).astype(np.int32)
+
 
     def updateVertices(self):
         """Update the mesh vertex positions using the values in self.mesh_vertices, and using the known dolfinx-induded reindexing
@@ -509,17 +503,16 @@ def buildMesh(self,tti):
         """Construct a new mesh at the given time index tti, and determine the vertex re-indexing induced by dolfinx
         """        
         self.tti = tti
-        self.buildVertices(time_index=tti, useFakeEncodedZ=True)
+        self.buildVertices(time_index=tti)
         self.constructMesh()
-        self.buildVertices(time_index=tti, useFakeEncodedZ=False)
         self.updateVertices()        
 
     def updateMesh(self,tti):
         """Construct the mesh positions at the given time index tti, and update the existing mesh with the new values
         """   
         assert self.mesh is not None
         self.tti = tti
-        self.buildVertices(time_index=tti, useFakeEncodedZ=False)
+        self.buildVertices(time_index=tti)
         self.updateVertices()        
 
     def buildHexahedra(self):
@@ -672,10 +665,8 @@ def constructMesh(self):
         # breakpoint()
 
         #
-        # obtain original vertex order as encoded in z-pos digits
-        zz  = self.mesh.geometry.x[:,2].copy()
-        zz2 = np.mod(zz,1000)
-        self.mesh_reindex = (1e-4+zz2*100).astype(np.int32)
+        # obtain original vertex order
+        self.mesh_reindex = np.array(self.mesh.geometry.input_global_indices).astype(np.int32)        
         self.mesh0_geometry_x = self.mesh.geometry.x.copy()
 
         # for i,c in enumerate(self.node_index):

warmth3d/sed_only_mesh_model.py

@@ -297,13 +297,11 @@ def kForLayerID(self, ss):
         return self.parameters1D.kCrust
 
 
-    def buildVertices(self, time_index=0, useFakeEncodedZ=False):
+    def buildVertices(self, time_index=0):
         """Determine vertex positions, node-by-node.
            For every node, the same number of vertices is added (one per sediment, one per crust, lith, asth, and one at the bottom)
            Degenerate vertices (e.g. at nodes where sediment is yet to be deposited or has been eroded) are avoided by a small shift, kept in self.sed_diff_z
            
-           When the option useFakeEncodedZ is set, the z-values are repurposed to encode the index of the vertex in the original, deterministic indexing.
-           This is necessary because dolfinx will re-index the vertices upon mesh generation. 
         """           
         tti = time_index
         self.tti = time_index
@@ -360,11 +358,6 @@ def buildVertices(self, time_index=0, useFakeEncodedZ=False):
         self.sed_diff_z = np.array(self.sed_diff_z)
         self.mesh_vertices = self.mesh_vertices_0.copy()
         self.mesh_vertices[:,2] = self.mesh_vertices_0[:,2] + self.sed_diff_z
-        if (useFakeEncodedZ):
-            self.mesh_vertices[:,2] = np.ceil(self.mesh_vertices[:,2])*1000 + np.array(list(range(self.mesh_vertices.shape[0])))*0.01
-            # zz = self.mesh_vertices[:,2].copy()
-            # zz2=np.mod(zz,1000)
-            # # mesh_reindex = (1e-4+zz2*10).astype(np.int32)
 
     def updateVertices(self):
         """Update the mesh vertex positions using the values in self.mesh_vertices, and using the known dolfinx-induded reindexing
@@ -381,17 +374,16 @@ def buildMesh(self,tti):
         """Construct a new mesh at the given time index tti, and determine the vertex re-indexing induced by dolfinx
         """        
         self.tti = tti
-        self.buildVertices(time_index=tti, useFakeEncodedZ=True)
+        self.buildVertices(time_index=tti)
         self.constructMesh()
-        self.buildVertices(time_index=tti, useFakeEncodedZ=False)
         self.updateVertices()        
 
     def updateMesh(self,tti):
         """Construct the mesh positions at the given time index tti, and update the existing mesh with the new values
         """   
         assert self.mesh is not None
         self.tti = tti
-        self.buildVertices(time_index=tti, useFakeEncodedZ=False)
+        self.buildVertices(time_index=tti)
         self.updateVertices()        
 
     def buildHexahedra(self):
@@ -532,10 +524,8 @@ def constructMesh(self):
             self.cell_data_layerID = aa.values.copy()
 
         #
-        # obtain original vertex order as encoded in z-pos digits
-        zz  = self.mesh.geometry.x[:,2].copy()
-        zz2 = np.mod(zz,1000)
-        self.mesh_reindex = (1e-4+zz2*100).astype(np.int32)
+        # obtain original vertex order
+        self.mesh_reindex = np.array(self.mesh.geometry.input_global_indices).astype(np.int32)        
         self.mesh0_geometry_x = self.mesh.geometry.x.copy()
 
     def normalizedZ(self, z):