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,