From 03aff2634f80125f5d5425b377f77e567c0a50a1 Mon Sep 17 00:00:00 2001 From: Jussi Aittoniemi Date: Mon, 25 Mar 2024 13:23:15 +0100 Subject: [PATCH 1/3] Read original global indices from mesh geometry instead of encoding them in positions; use explicit float32 arrays for RESQML properties; some code cleanup --- tests/warmth3d/test_3d.py | 2 +- warmth/build.py | 7 +++---- warmth/forward_modelling.py | 3 ++- warmth/mesh_model.py | 23 +++++------------------ warmth/resqpy_helpers.py | 18 +++++++++--------- warmth3d/fixed_mesh_model.py | 4 ++-- warmth3d/sed_only_mesh_model.py | 6 +++--- 7 files changed, 25 insertions(+), 38 deletions(-) diff --git a/tests/warmth3d/test_3d.py b/tests/warmth3d/test_3d.py index c280115..f214a57 100644 --- a/tests/warmth3d/test_3d.py +++ b/tests/warmth3d/test_3d.py @@ -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: diff --git a/warmth/build.py b/warmth/build.py index 3ad2044..028ebc4 100644 --- a/warmth/build.py +++ b/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: diff --git a/warmth/forward_modelling.py b/warmth/forward_modelling.py index d461b3b..1d72ea2 100644 --- a/warmth/forward_modelling.py +++ b/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] diff --git a/warmth/mesh_model.py b/warmth/mesh_model.py index 36b0d7e..503432c 100644 --- a/warmth/mesh_model.py +++ b/warmth/mesh_model.py @@ -515,13 +515,10 @@ def rhpForLayerID(self, ss:int, node_index:int)->float: return 0 - def buildVertices(self, time_index=0, useFakeEncodedZ=False, optimized=False): + def buildVertices(self, time_index=0, optimized=False): """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 @@ -548,14 +545,10 @@ def buildVertices(self, time_index=0, useFakeEncodedZ=False, optimized=False): zpos = xxT + base_of_current_sediments aa = np.concatenate([aa,np.array([zpos])]) - # zp = base_of_last_sediments+ (base_crust-base_of_last_sediments)*(i/self.numElemInCrust) base_of_last_sediments = (xxT+self.bottom_sed_id_at_nodes[-1][:,tti]) if (len(self.bottom_sed_id_at_nodes)>0) else xxT for i in range(1,self.numElemInCrust+1): zp = base_of_last_sediments+ (bc-base_of_last_sediments)*(i/self.numElemInCrust) - # self.mesh_vertices_0.append( [ node.X, node.Y, base_of_last_sediments+ (base_crust-base_of_last_sediments)*(i/self.numElemInCrust) ] ) - # self.sed_diff_z.append(0.0) - # self.mesh_vertices_age_unsorted.append(1000) aa = np.concatenate([aa,np.array([zp])]) for i in range(1,self.numElemInLith+1): @@ -595,8 +588,6 @@ def buildVertices(self, time_index=0, useFakeEncodedZ=False, optimized=False): self.sed_diff_z.append(-self.minimumCellThick*(self.numberOfSedimentCells - (ss*self.numElemPerSediment+j) )) age_of_previous = node.sediments.baseage[ss-1] if (ss>0) else 0.0 self.mesh_vertices_age_unsorted.append( age_of_previous + ((j+1) / self.numElemPerSediment) * (node.sediments.baseage[ss]-age_of_previous) ) # append interpolatedbase age of current sediment - # if (ind==97) and (tti==0): - # breakpoint() if (ind==0): delta = time.time() - st print("delta 2", delta) @@ -630,8 +621,6 @@ def buildVertices(self, time_index=0, useFakeEncodedZ=False, optimized=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 def updateVertices(self): """Update the mesh vertex positions using the values in self.mesh_vertices, and using the known dolfinx-induded reindexing @@ -651,7 +640,7 @@ def buildMesh(self,tti:int): """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) logger.info("Built vertices") self.constructMesh() logger.info("Built mesh") @@ -665,7 +654,7 @@ def updateMesh(self,tti:int, optimized=False): import time assert self.mesh is not None self.tti = tti - self.buildVertices(time_index=tti, useFakeEncodedZ=False, optimized=optimized) + self.buildVertices(time_index=tti, optimized=optimized) self.updateVertices() self.posarr.append(self.mesh.geometry.x.copy()) self.time_indices.append(self.tti) @@ -789,10 +778,8 @@ def mpi_print(s): mpi_print(f"Dir local verts: {dir(self.mesh.topology.index_map(0))}") # local_to_global ! ? mpi_print(f"Type local verts: {type(self.mesh.topology.index_map(0))}") - # 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) + # store 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() diff --git a/warmth/resqpy_helpers.py b/warmth/resqpy_helpers.py index adea5ff..c72430f 100644 --- a/warmth/resqpy_helpers.py +++ b/warmth/resqpy_helpers.py @@ -382,7 +382,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 +392,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 +413,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 +422,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 +431,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 +443,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 +452,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 +563,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 +574,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', diff --git a/warmth3d/fixed_mesh_model.py b/warmth3d/fixed_mesh_model.py index 88109b0..f8bba47 100644 --- a/warmth3d/fixed_mesh_model.py +++ b/warmth3d/fixed_mesh_model.py @@ -488,7 +488,7 @@ def buildVertices(self, time_index=0, useFakeEncodedZ=False): 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 + self.mesh_vertices[:,2] = np.ceil(self.mesh_vertices[:,2])*1000 + np.array(list(range(self.mesh_vertices.shape[0])))*0.001 # zz = self.mesh_vertices[:,2].copy() # zz2=np.mod(zz,1000) # # mesh_reindex = (1e-4+zz2*10).astype(np.int32) @@ -675,7 +675,7 @@ def constructMesh(self): # 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) + self.mesh_reindex = (1e-4+zz2*1000).astype(np.int32) self.mesh0_geometry_x = self.mesh.geometry.x.copy() # for i,c in enumerate(self.node_index): diff --git a/warmth3d/sed_only_mesh_model.py b/warmth3d/sed_only_mesh_model.py index 5f1b62c..29ad509 100644 --- a/warmth3d/sed_only_mesh_model.py +++ b/warmth3d/sed_only_mesh_model.py @@ -361,7 +361,7 @@ def buildVertices(self, time_index=0, useFakeEncodedZ=False): 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 + self.mesh_vertices[:,2] = np.ceil(self.mesh_vertices[:,2])*10000 + np.array(list(range(self.mesh_vertices.shape[0])))*0.001 # zz = self.mesh_vertices[:,2].copy() # zz2=np.mod(zz,1000) # # mesh_reindex = (1e-4+zz2*10).astype(np.int32) @@ -534,8 +534,8 @@ def constructMesh(self): # # 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) + zz2 = np.mod(zz,10000) + self.mesh_reindex = (1e-4+zz2*1000).astype(np.int32) self.mesh0_geometry_x = self.mesh.geometry.x.copy() def normalizedZ(self, z): From 5bed34756a0768c56612da8080482feee0b2dc66 Mon Sep 17 00:00:00 2001 From: Jussi Aittoniemi Date: Mon, 1 Apr 2024 10:10:33 +0200 Subject: [PATCH 2/3] Read original global indices from mesh geometry also in the legacy mesh classes --- warmth3d/fixed_mesh_model.py | 21 ++++++--------------- warmth3d/sed_only_mesh_model.py | 20 +++++--------------- 2 files changed, 11 insertions(+), 30 deletions(-) diff --git a/warmth3d/fixed_mesh_model.py b/warmth3d/fixed_mesh_model.py index f8bba47..d68c5c4 100644 --- a/warmth3d/fixed_mesh_model.py +++ b/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.001 - # 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,9 +503,8 @@ 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): @@ -519,7 +512,7 @@ def updateMesh(self,tti): """ 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*1000).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): diff --git a/warmth3d/sed_only_mesh_model.py b/warmth3d/sed_only_mesh_model.py index 29ad509..9e69431 100644 --- a/warmth3d/sed_only_mesh_model.py +++ b/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])*10000 + np.array(list(range(self.mesh_vertices.shape[0])))*0.001 - # 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,9 +374,8 @@ 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): @@ -391,7 +383,7 @@ def updateMesh(self,tti): """ 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,10000) - self.mesh_reindex = (1e-4+zz2*1000).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): From c2fb5db12608d975660aab80cb813476478dadca Mon Sep 17 00:00:00 2001 From: Jussi Aittoniemi Date: Tue, 26 Mar 2024 10:41:30 +0100 Subject: [PATCH 3/3] Replace prints with logger in mesh_model --- tests/warmth3d/test_3d.py | 2 +- warmth/mesh_model.py | 76 +++++++++++++++++---------------------- warmth/resqpy_helpers.py | 11 ++---- 3 files changed, 36 insertions(+), 53 deletions(-) diff --git a/tests/warmth3d/test_3d.py b/tests/warmth3d/test_3d.py index f214a57..894573b 100644 --- a/tests/warmth3d/test_3d.py +++ b/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 diff --git a/warmth/mesh_model.py b/warmth/mesh_model.py index 503432c..dfd30a7 100644 --- a/warmth/mesh_model.py +++ b/warmth/mesh_model.py @@ -143,7 +143,7 @@ def boundary(x): lid_to_keep.append(lid0) cell_id_to_keep.append(self.node_index[i]) if abs(self.node_index[i].Y-minY)>1: - print("unusual Y coordinate:", minY, self.node1D[self.node_index[i]].Y, i, self.node_index[i], self.node1D[self.node_index[i]]) + logger.warning( f"unusual Y coordinate:, {minY}, {self.node1D[self.node_index[i]].Y}, {i}, {self.node_index[i]}, {self.node1D[self.node_index[i]]}") for ti in t: p_to_keep.add(ti) poro0_per_cell = np.array( [ self.getSedimentPropForLayerID('phi', lid,cid) for lid,cid in zip(lid_to_keep,cell_id_to_keep) ] ) @@ -213,7 +213,7 @@ def receive_mpi_messages(self): self.age_per_vertex[val] = self.mesh_vertices_age_s[k][ind] delta = time.time() - st - print("receive_mpi_messages delta", delta) + logger.debug( f"receive_mpi_messages delta: {delta}") def get_node_pos_and_temp(self, tti=-1): # @@ -258,14 +258,8 @@ def write_hexa_mesh_resqml( self, out_path, tti): 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]) + logger.warning( f"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() - # if (minY>40000): - # pp = self.getSedimentPropForLayerID('phi', lid0, hex_data_nodeID[i]) - # if (pp<0.7) and lid0>=0: - # print("weird phi: ", pp, 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]) # the actual, Sekiguchi-derived conductivitues @@ -422,17 +416,8 @@ def getSedimentPropForLayerID(self, property, layer_id:int, node_index:int) ->fl """ assert property in ['k_cond', 'rhp', 'phi', 'decay', 'solidus', 'liquidus'], "Unknown property " + property if (layer_id>=0) and (layer_id0.7 and node.Y<40000: - # print("phi", property, phi, node_index, node) - # breakpoint() - # if (property=='phi') and phi<0.7 and node.Y>40000: - # print("phi", property, phi, node_index, node) - # breakpoint() - # phi = self.globalSediments[property][layer_id] - # assert abs(phi-phi0)<1e-6 return prop if (layer_id<=-1) and (layer_id>=-3): lid = -layer_id -1 @@ -590,7 +575,7 @@ def buildVertices(self, time_index=0, optimized=False): self.mesh_vertices_age_unsorted.append( age_of_previous + ((j+1) / self.numElemPerSediment) * (node.sediments.baseage[ss]-age_of_previous) ) # append interpolatedbase age of current sediment if (ind==0): delta = time.time() - st - print("delta 2", delta) + logger.debug(f"delta 2: {delta}") st = time.time() if not self.runSedimentsOnly: base_of_last_sediments = bottom_sed_id(node, self.numberOfSediments-1, tti) if (self.numberOfSediments>0) else top_of_sediments @@ -614,7 +599,7 @@ def buildVertices(self, time_index=0, optimized=False): self.mesh_vertices_age_unsorted.append(1000) if (ind==0): delta = time.time() - st - print("delta 3", delta) + logger.debug(f"delta 3: {delta}") assert len(self.mesh_vertices_0) % self.num_nodes ==0 self.mesh_vertices_0 = np.array(self.mesh_vertices_0) @@ -756,6 +741,7 @@ def constructMesh(self): def mpi_print(s): logger.info(f"Rank {comm.rank}: {s}") + fn = self.modelName+"_mesh.xdmf" if comm.rank==0: mesh.write( fn ) @@ -981,7 +967,7 @@ def updateDirichletBaseTemperature(self): def updateDBC(self): self.averageLABdepth = np.mean(np.array([ top_asth(n, self.tti) for n in self.node1D])) - print("self.averageLABdepth", self.averageLABdepth) + logger.debug(f"self.averageLABdepth: {self.averageLABdepth}") ii = np.where(self.mesh.geometry.x[:,2]>250000) self.bc.value.x.array[ii] = 1330+(self.mesh.geometry.x[ii,2]-self.averageLABdepth)*0.0003 @@ -996,7 +982,7 @@ def buildDirichletBC(self): st = time.time() self.averageLABdepth = np.mean(np.array([ top_asth(n, self.tti) for n in self.node1D])) self.Zmax = self.averageLABdepth + 130000 - print("buildDirichletBC delta 1 ", time.time()-st) + logger.debug(f"buildDirichletBC delta 1: {time.time()-st}") def boundary_D_top_bottom(x): subs0 = self.getSubsidenceAtMultiplePos(x[0,:], x[1,:]) @@ -1014,17 +1000,17 @@ def boundary_D_top(x): if (self.useBaseFlux): st = time.time() dofs_D = dolfinx.fem.locate_dofs_geometrical(self.V, boundary_D_top) - print("buildDirichletBC delta 2 ", time.time()-st) + logger.debug(f"buildDirichletBC delta 2: {time.time()-st}") else: st = time.time() dofs_D = dolfinx.fem.locate_dofs_geometrical(self.V, boundary_D_top_bottom) - dofs_D2 = dolfinx.fem.locate_dofs_geometrical(self.V, boundary_D_bottom) - dofs_D3 = dolfinx.fem.locate_dofs_geometrical(self.V, boundary_D_top) - print("buildDirichletBC delta 3 ", time.time()-st) + # dofs_D2 = dolfinx.fem.locate_dofs_geometrical(self.V, boundary_D_bottom) + # dofs_D3 = dolfinx.fem.locate_dofs_geometrical(self.V, boundary_D_top) + logger.debug(f"buildDirichletBC delta 3: {time.time()-st}") u_bc = dolfinx.fem.Function(self.V) st = time.time() u_bc.interpolate(self.TemperatureGradient) - print("buildDirichletBC delta 4 ", time.time()-st) + logger.debug(f"buildDirichletBC delta 4: {time.time()-st}") bc = dolfinx.fem.dirichletbc(u_bc, dofs_D) return bc @@ -1100,7 +1086,7 @@ def mpi_print(s): nn = self.node1D[i] iix = np.where(self.node_index==i) self.subsidence_at_nodes[iix,:] = nn.subsidence - print("setup delay 1.3", time.time()-st) + logger.debug(f"setup delay 1.3: {time.time()-st}") st = time.time() top_of_sediments = self.subsidence_at_nodes @@ -1119,12 +1105,12 @@ def mpi_print(s): self.top_sed_at_nodes[i,:] = nn.subsidence self.base_crust_at_nodes[i,:] = nn.subsidence + nn.sed_thickness_ls + nn.crust_ls self.base_lith_at_nodes[i,:] = self.base_crust_at_nodes[i,:] + nn.lith_ls - print("setup delay 1.4", time.time()-st) + logger.debug(f"setup delay 1.4: {time.time()-st}") st = time.time() self.thermalCond, self.c_rho, self.layerIDsFcn, self.rhpFcn = self.buildKappaAndLayerIDs() assert not np.any(np.isnan(self.thermalCond.x.array)) - print("setup delay 1", time.time()-st) + logger.debug(f"setup delay 1: {time.time()-st}") st = time.time() # initialise both with initial condition: either a step function, or the solution from another Model instance @@ -1133,8 +1119,7 @@ def mpi_print(s): else: self.u_n.x.array[:] = initial_state_model.uh.x.array[:].copy() self.uh.x.array[:] = self.u_n.x.array[:].copy() - print("setup delay 1.5", time.time()-st) - + logger.debug(f"setup delay 1.5: {time.time()-st}") def setupSolverAndSolve(self, n_steps:int=100, time_step:int=-1, skip_setup:bool = False, initial_state_model = None, update_bc=False): """ Sets up the function spaces, output functions, input function (kappa values), boundary conditions, initial conditions. @@ -1144,7 +1129,7 @@ def setupSolverAndSolve(self, n_steps:int=100, time_step:int=-1, skip_setup:bool """ comm = MPI.COMM_WORLD def mpi_print(s): - print(f"Rank {comm.rank}: {s}") + logger.debug(f"Rank {comm.rank}: {s}") if (not skip_setup): self.setupSolver(initial_state_model) @@ -1155,7 +1140,7 @@ def mpi_print(s): import time st = time.time() self.sedimentsConductivitySekiguchi() - # print("solve delay 2", time.time()-st) + logger.debug(f"solve delay 2: {time.time()-st}") t=0 dt = time_step if (time_step>0) else 3600*24*365 * 5000000 @@ -1184,6 +1169,9 @@ def mpi_print(s): domain_c = dolfinx.fem.Function(self.V) domain_c.x.array[ : ] = 0.0 if (self.CGorder>1): + # + # NOTE: CGorder>1 is under development, not functional + # def marker(x): print(x.shape, x) return x[2,:]>3990 @@ -1216,7 +1204,7 @@ def marker(x): domain_zero = dolfinx.fem.Function(self.V) toppos = self.getSubsidenceAtMultiplePos(self.mesh.geometry.x[:,0], self.mesh.geometry.x[:,1]) domain_zero.x.array[ self.mesh.geometry.x[:,2] < toppos+0.01 ] = 1 - print("Neumann conditions: ", self.tti, np.count_nonzero(domain_c.x.array), np.count_nonzero(domain_zero.x.array)) + logger.debug(f"Neumann conditions: , {self.tti}, {np.count_nonzero(domain_c.x.array)}, {np.count_nonzero(domain_zero.x.array)}") g = (-1.0*baseFlux) * ufl.conditional( domain_c > 0, 1.0, 0.0 ) L = (self.c_rho*self.u_n + dt*f)*v*ufl.dx - dt * g * v * ufl.ds # last term reflects Neumann BC @@ -1289,7 +1277,6 @@ def safeInterpolation(self, interp, pos_x, pos_y, epsilon=1e-2): interp([pos_x+epsilon, pos_y+epsilon])[0]]) res = np.nanmean(manyres) if (np.isnan(res)): - # print(pos_x, pos_y, interp([pos_x, pos_y]), interp([0,0])[0] ) logger.warning(f'NaN encounered in safeInterpolation pos_x {pos_x}: pos_y: {pos_y}; {interp([pos_x, pos_y])} {interp([0,0])[0]} ') # assert not np.isnan(res), "interpolation is nan in safeInterpolation" return res @@ -1529,19 +1516,19 @@ def interpolateResult(self, x): break res.append( val ) else: - print("need to extrapolate cell for point", i, point) + logger.debug(f"need to extrapolate cell for point {i}, {point}") if (point[2]>200000): try: points_cells.append(cell_candidates.links(i)[0]) points_on_proc.append(point) except IndexError: - print("IndexError", point, cell_candidates.links(i) ) + logger.warning(f"IndexError, {point}, {cell_candidates.links(i)}" ) breakpoint() raise else: - print("PING V", point) + logger.debug(f"PING V, {point}") if len(cell_candidates.links(i))==0: - print("PING V V", point) + logger.debug(f"PING V V, {point}") def boundary(x): return np.full(x.shape[1], True) #entities = dolfinx.mesh.locate_entities(self.mesh, 3, boundary ) @@ -1552,7 +1539,7 @@ def boundary(x): aa = np.any(np.isnan(res)) bb = np.any(np.isnan(self.uh.x.array)) if aa or bb: - print(aa,bb) + logger.debug(f"aa {aa}, bb {bb}") breakpoint() if transpose: @@ -1564,6 +1551,7 @@ def boundary(x): def run_3d( builder:Builder, parameters:Parameters, start_time=182, end_time=0, pad_num_nodes=0, out_dir = "out-mapA/",sedimentsOnly=False, writeout=True, base_flux=None): + logger.setLevel(10) # numeric level equals DEBUG comm = MPI.COMM_WORLD builder=interpolate_all_nodes(builder) nums = 4 @@ -1613,7 +1601,7 @@ def run_3d( builder:Builder, parameters:Parameters, start_time=182, end_time=0, mms_tti.append(tti) logger.info(f"Simulated time step {tti}") bar.next() - print("total time solve 3D: " , time_solve) + logger.info(f"total time solve 3D: {time_solve}") if (writeout_final): comm.Barrier() if comm.rank>=1: @@ -1625,8 +1613,8 @@ def run_3d( builder:Builder, parameters:Parameters, start_time=182, end_time=0, if comm.rank==0: mm2.receive_mpi_messages() EPCfilename = mm2.write_hexa_mesh_resqml("temp/", end_time) - print("RESQML model written to: " , EPCfilename) + logger.info(f"RESQML model written to: {EPCfilename}") EPCfilename_ts = mm2.write_hexa_mesh_timeseries("temp/") - print("RESQML partial model with timeseries written to: ", EPCfilename_ts) + logger.info(f"RESQML partial model with timeseries written to: {EPCfilename_ts}") read_mesh_resqml_hexa(EPCfilename) # test reading of the .epc file return mm2 diff --git a/warmth/resqpy_helpers.py b/warmth/resqpy_helpers.py index c72430f..14a2edc 100644 --- a/warmth/resqpy_helpers.py +++ b/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,7 +44,6 @@ 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 @@ -56,13 +51,13 @@ def read_mesh_resqml(epcfilename, meshTitle = 'tetramesh'): # 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,