Merge pull request #18 from equinor/replace_print_with_logger

Replace prints with logger in mesh_model

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

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_id<self.numberOfSediments):
-            # node_index = ind2 // (self.numberOfSediments+6)  # +6 because crust, lith, aest are each cut into two
             node = self.node1D[node_index]
             prop = node.sediments[property][layer_id]
-            # if (property=='phi') and phi>0.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

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,