Add resqml timeseries code

tests/warmth3d/test_3d.py

@@ -71,7 +71,7 @@ def test_3d_compare():
     comm.Barrier()
     # if comm.rank == 0:            
     # mm2 = 
-    run_3d(model.builder,model.parameters,start_time=model.parameters.time_start,end_time=0, pad_num_nodes=2,writeout=False, base_flux=None)
+    mm2, posarr, Tarr = run_3d(model.builder,model.parameters,start_time=model.parameters.time_start,end_time=0, pad_num_nodes=2,writeout=False, base_flux=None)
 
 
     nnx = (model.builder.grid.num_nodes_x+2*mm2.padX) 
@@ -99,6 +99,9 @@ def test_3d_compare():
     max_abs_error = np.amax(np.abs(temp_1d_at_mesh_pos-temp_3d_mesh))
     max_abs_error_shallow = np.amax(np.abs(temp_1d_at_mesh_pos[dd_subset]-temp_3d_mesh[dd_subset]))    
 
+    print("Max error overall:", max_abs_error)
+    print("Max error <5000m:", max_abs_error_shallow)
+
     # assert (max_abs_error<25.0), "Temperature difference between 3D and 1D simulations is >25"
     # assert (max_abs_error_shallow<5.0),  "Temperature difference between 3D and 1D simulations is >5 in the sediments."
 

warmth/mesh_model.py

@@ -11,7 +11,7 @@
 from .parameters import Parameters
 from warmth.logging import logger
 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_properties,read_mesh_resqml_hexa
+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():
     #Homemade version of matlab tic and toc functions
     import time
@@ -233,13 +233,6 @@ def write_hexa_mesh_resqml( self, out_path):
             minY = np.amin(np.array ( [np.array(points_cached)[hi,1] for hi in h] ))
             poro0 = poro0_per_cell[i]
             lid0  = lid_to_keep[i]
-            # if (minY>40000) and poro0 < 0.8 and lid0>=0:
-            #     print("problem A", minY, poro0, i, h)
-            #     breakpoint()
-            # if (minY<40000) and poro0 > 0.8:
-            #     print("problem B", minY, poro0, i, h)
-            #     breakpoint()
-
 
         T_per_vertex = [ self.uh.x.array[reverse_reindex_order[i]] for i in range(self.mesh.geometry.x.shape[0]) if i in p_to_keep  ]
         age_per_vertex = [ self.mesh_vertices_age[reverse_reindex_order[i]] for i in range(self.mesh.geometry.x.shape[0]) if i in p_to_keep  ]
@@ -251,6 +244,86 @@ def write_hexa_mesh_resqml( self, out_path):
             cond_per_cell, rhp_per_cell, lid_per_cell)
         return filename_hex
 
+    def write_hexa_mesh_timeseries( self, out_path, posarr, Tarr):
+        """Prepares arrays and calls the RESQML output helper function for hexa meshes:  the lith and aesth are removed, and the remaining
+           vertices and cells are renumbered;  the sediment properties are prepared for output.
+
+           out_path: string: path to write the resqml model to (.epc and .h5 files)
+
+           returns the filename (of the .epc file) that was written 
+        """            
+        x_original_order = self.mesh.geometry.x[:].copy()
+        reverse_reindex_order = np.arange( self.mesh_vertices.shape[0] )
+        for ind,val in enumerate(self.mesh_reindex):
+            x_original_order[val,:] = self.mesh.geometry.x[ind,:]
+            reverse_reindex_order[val] = ind
+        hexaHedra, hex_data_layerID, hex_data_nodeID = self.buildHexahedra(keep_padding=False)
+
+        hexa_to_keep = []
+        p_to_keep = set()
+        lid_to_keep = []
+        cond_per_cell = []
+        cell_id_to_keep = []
+        for i,h in enumerate(hexaHedra):
+            lid0 = hex_data_layerID[i]
+            # 
+            # discard aesth and lith (layer IDs -2, -3)
+            #
+            if (lid0>=-1) and (lid0<100):
+                hexa_to_keep.append(h)
+                lid_to_keep.append(lid0)
+                # cell_id_to_keep.append(self.node_index[i])
+                cell_id_to_keep.append(hex_data_nodeID[i])
+                # minY = np.amin(np.array ( [x_original_order[hi,1] for hi in h] ))
+                # if abs( self.node1D[hex_data_nodeID[i]].Y - minY)>1:
+                #     print("weird Y:", minY, self.node1D[hex_data_nodeID[i]].Y, abs( self.node1D[hex_data_nodeID[i]].Y - minY), i, hex_data_nodeID[i])
+                #     breakpoint()
+                # k_cond_mean = []
+                for hi in h:
+                    p_to_keep.add(hi)
+                #     k_cond_mean.append(self.thermalCond.x.array[hi])
+                # cond_per_cell.append( np.mean(np.array(k_cond_mean)))
+
+        # poro0_per_cell = np.array( [ self.getSedimentPropForLayerID('phi', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ] )
+        # decay_per_cell = np.array( [ self.getSedimentPropForLayerID('decay', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
+        # density_per_cell = np.array( [ self.getSedimentPropForLayerID('solidus', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
+        # cond_per_cell = np.array( [ self.getSedimentPropForLayerID('k_cond', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
+        # rhp_per_cell = np.array( [ self.getSedimentPropForLayerID('rhp', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ])
+        # lid_per_cell = np.array(lid_to_keep)
+
+        points_cached_series = []
+        for j in range(len(posarr)):
+            x_original_order = posarr[j].copy()
+            # reverse_reindex_order = np.arange( self.mesh_vertices.shape[0] )
+            for ind,val in enumerate(self.mesh_reindex):
+                x_original_order[val,:] = posarr[j][ind,:]            
+            points_cached = []
+            point_original_to_cached = np.ones(self.mesh.geometry.x.shape[0], dtype = np.int32)  * (-1)
+            for i in range(self.mesh.geometry.x.shape[0]):
+                if (i in p_to_keep):
+                    point_original_to_cached[i] = len(points_cached)
+                    points_cached.append(x_original_order[i,:])
+            points_cached_series.append(np.array(points_cached))
+        hexa_renumbered = [ [point_original_to_cached[i] for i in hexa] for hexa in hexa_to_keep ]
+
+        # for i,h in enumerate(hexa_renumbered):
+        #     minY = np.amin(np.array ( [np.array(points_cached)[hi,1] for hi in h] ))
+        #     poro0 = poro0_per_cell[i]
+        #     lid0  = lid_to_keep[i]
+
+        Temp_per_vertex_series = []
+        for j in range(len(Tarr)):            
+            # T_per_vertex = [ self.uh.x.array[reverse_reindex_order[i]] for i in range(self.mesh.geometry.x.shape[0]) if i in p_to_keep  ]
+            T_per_vertex = [ Tarr[j][reverse_reindex_order[i]] for i in range(self.mesh.geometry.x.shape[0]) if i in p_to_keep  ]
+            Temp_per_vertex_series.append(np.array(T_per_vertex))
+        # age_per_vertex = [ self.mesh_vertices_age[reverse_reindex_order[i]] for i in range(self.mesh.geometry.x.shape[0]) if i in p_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 )
+        return filename_hex
+
     def heatflow_at_crust_sed_boundary(self):
         hf_array = np.zeros([self.num_nodes_x-2*self.padX, self.num_nodes_y-2*self.padX])
         for hy in range(self.padX, self.num_nodes_y-self.padX):
@@ -1625,6 +1698,7 @@ def run_3d( builder:Builder, parameters:Parameters,  start_time=182, end_time=0,
     mms_tti = []
     tti = 0
     # base_flux = 0.0033
+    writeout_final = True
     time_solve = 0.0    
     posarr = []
     Tarr = []
@@ -1673,8 +1747,10 @@ def run_3d( builder:Builder, parameters:Parameters,  start_time=182, end_time=0,
             logger.info(f"Simulated time step {tti}")
             bar.next()
     print("total time solve: " , time_solve)
-    if (writeout):
+    if (writeout_final):
         EPCfilename = mm2.write_hexa_mesh_resqml("temp/")
         print("RESQML model written to: " , EPCfilename)
+        EPCfilename_ts = mm2.write_hexa_mesh_timeseries("temp/", posarr, Tarr)
+        print("RESQML partial model with timeseries written to: ", EPCfilename_ts)
         read_mesh_resqml_hexa(EPCfilename)  # test reading of the .epc file
     return mm2,posarr,Tarr

warmth/mesh_utils.py

@@ -113,6 +113,18 @@ def volumeOfTet(points):
     bdcd = np.cross(bd,cd)
     return np.linalg.norm(np.dot(ad,bdcd))/6
 
+
+def volumeOfHex(points):
+    """ Computes the volume of a hexahedron, given as eight 3D-points
+    """ 
+    import numpy as np
+    tetsplit1 = [ [1,2,4,8], [1,2,5,8], [4,8,2,3], [2,3,7,8], [2,5,6,8], [2,6,7,8] ]
+    vol = 0.0
+    for f in tetsplit1:
+        tet = points[[p-1 for p in f],:]
+        vol = vol + volumeOfTet(tet)
+    return vol
+
 def interpolateNode(interpolationNodes: List[single_node], interpolationWeights=None) -> single_node:
     assert len(interpolationNodes)>0
     if interpolationWeights is None:

warmth/resqpy_helpers.py

@@ -463,4 +463,129 @@ def write_hexa_grid_with_properties(filename, nodes, cells, modelTitle = "hexame
     model.store_epc()
 
 
+def write_hexa_grid_with_timeseries(filename, nodes_series, cells, modelTitle = "hexamesh",
+    Temp_per_vertex_series=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.
+ 
+       cells is an array of 8-arrays in which the nodes are ordered:     
+               7------6
+              /      /|
+             /      / |
+            4------5  |
+            |         |
+            |  3------2
+            | /      /
+            |/      /
+            0------1
+
+       NOTE: writing properties that are defines per-node (have 'nodes' as indexable element) requires a patched version of resqpy!
+    """
+
+    nodes = nodes_series[0]
+    node_count = len(nodes)
+    faces_per_cell = []
+    nodes_per_face = []
+    faces_dict = {}
+    faces_repeat = np.zeros(node_count*100, dtype = bool)
+
+    cell_face_is_right_handed = np.zeros( len(cells)*6, dtype = bool)
+    for ih,hexa in enumerate(cells):
+        faces= [[0,3,2,1], [0,1,5,4], [1,2,6,5], [2,3,7,6], [3,0,4,7], [4,5,6,7]]
+        for iq,quad in enumerate(faces):
+            face0 = [hexa[x] for x in quad ]
+            assert -1 not in face0
+            fkey0 = ( x for x in sorted(face0) )
+            #
+            # keep track of which faces are encountered once vs. more than once
+            # faces that are encountered the second time will need to use the reverse handedness
+            #
+            face_is_repeated = False
+            if (fkey0 not in faces_dict):
+                faces_dict[fkey0] = len(nodes_per_face)
+                nodes_per_face.extend(face0)
+                cell_face_is_right_handed[(ih*6 + iq)] = False
+            else:
+                face_is_repeated = True
+                cell_face_is_right_handed[(ih*6 + iq)] = True
+            fidx0 = faces_dict.get(fkey0)            
+            faces_per_cell.append(fidx0/4)
+            faces_repeat[int(fidx0/4)] = face_is_repeated
+
+    set_cell_count = int(len(faces_per_cell)/6)
+    face_count = int(len(nodes_per_face)/4)
+
+
+    model = rq.new_model(filename)
+    crs = rqc.Crs(model)
+    crs.create_xml()
+
+    times_in_years = [ max(int(t*1e6),1) for t in list(range(len(nodes_series)-1,-1,-1))]
+    gts = rts.GeologicTimeSeries.from_year_list(model, times_in_years, title="warmth simulation")
+    gts.create_xml()
+    rts.timeframe_for_time_series_uuid(model, gts.uuid)
+
+    # create an empty HexaGrid
+    hexa = rug.HexaGrid(model, title = modelTitle)
+    assert hexa.cell_shape == 'hexahedral'
+
+    # hand craft all attribute data
+    hexa.crs_uuid = model.uuid(obj_type = 'LocalDepth3dCrs')
+    assert hexa.crs_uuid is not None
+    assert bu.matching_uuids(hexa.crs_uuid, crs.uuid)
+    hexa.set_cell_count(set_cell_count)
+    # faces
+    hexa.face_count = face_count
+    hexa.faces_per_cell_cl = np.arange(6, 6 * set_cell_count + 1, 6, dtype = int)
+    hexa.faces_per_cell = np.array(faces_per_cell)
+
+    # nodes
+    hexa.node_count = node_count
+    hexa.nodes_per_face_cl = np.arange(4, 4 * face_count + 1, 4, dtype = int)
+    hexa.nodes_per_face = np.array(nodes_per_face)
+
+    # face handedness
+    hexa.cell_face_is_right_handed = cell_face_is_right_handed  # False for all faces for external cells
+
+    # points
+    hexa.points_cached = nodes
+
+    # basic validity check
+    hexa.check_hexahedral()
+
+    hexa.create_xml()
+    hexa.write_hdf5()
+
+    if hexa.property_collection is None:
+        hexa.property_collection = rqp.PropertyCollection(support = hexa)
+    pc = hexa.property_collection
+
+    # 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]
+        pc.add_cached_array_to_imported_list(nodes2,
+                                                'dynamic nodes',
+                                                'points',
+                                                uom = 'm',
+                                                property_kind = 'length',
+                                                realization = 0,
+                                                time_index = time_index,
+                                                indexable_element = 'nodes',
+                                                points = True)
+        # active_array = np.ones([2160], dtype = bool)
+        tt = Temp_per_vertex_series[time_index]
+        pc.add_cached_array_to_imported_list(tt,
+                                                'Temperature',
+                                                'Temperature',
+                                                uom = 'degC',
+                                                property_kind = 'thermodynamic temperature',
+                                                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)
+
+    model.store_epc()