diff --git a/docs/notebooks/data/0.gri b/docs/notebooks/data/0.gri index 120ca3b..37d352e 100644 Binary files a/docs/notebooks/data/0.gri and b/docs/notebooks/data/0.gri differ diff --git a/docs/notebooks/data/100.gri b/docs/notebooks/data/100.gri index b899d47..1b24b51 100644 Binary files a/docs/notebooks/data/100.gri and b/docs/notebooks/data/100.gri differ diff --git a/docs/notebooks/data/163.gri b/docs/notebooks/data/163.gri index 8bbc9a4..adab1d6 100644 Binary files a/docs/notebooks/data/163.gri and b/docs/notebooks/data/163.gri differ diff --git a/docs/notebooks/data/168.gri b/docs/notebooks/data/168.gri index 680f5c4..64f91b8 100644 Binary files a/docs/notebooks/data/168.gri and b/docs/notebooks/data/168.gri differ diff --git a/docs/notebooks/data/170.gri b/docs/notebooks/data/170.gri index 9e6d4aa..90529c4 100644 Binary files a/docs/notebooks/data/170.gri and b/docs/notebooks/data/170.gri differ diff --git a/docs/notebooks/data/182.gri b/docs/notebooks/data/182.gri index cbd108f..496b294 100644 Binary files a/docs/notebooks/data/182.gri and b/docs/notebooks/data/182.gri differ diff --git a/docs/notebooks/data/66.gri b/docs/notebooks/data/66.gri index a075e84..4f0a8a3 100644 Binary files a/docs/notebooks/data/66.gri and b/docs/notebooks/data/66.gri differ diff --git a/tests/warmth3d/test_3d.py b/tests/warmth3d/test_3d.py new file mode 100644 index 0000000..2d2a3b8 --- /dev/null +++ b/tests/warmth3d/test_3d.py @@ -0,0 +1,97 @@ +if __name__ == '__main__': + import warmth + from pathlib import Path + + maps_dir = Path("./docs/notebooks/data/") + model = warmth.Model() + + inputs = model.builder.input_horizons_template + import os + print(os.getcwd()) + #Add surface grids to the table. You can use other method as well + inputs.loc[0]=[0,"0.gri",None,"Onlap"] + inputs.loc[1]=[66,"66.gri",None,"Onlap"] + inputs.loc[2]=[100,"100.gri",None,"Onlap"] + inputs.loc[3]=[163,"163.gri",None,"Erosive"] + inputs.loc[4]=[168,"168.gri",None,"Erosive"] + inputs.loc[5]=[170,"170.gri",None,"Onlap"] + inputs.loc[6]=[182,"182.gri",None,"Erosive"] + model.builder.input_horizons=inputs + + + inc = 2000 + model.builder.define_geometry(maps_dir/"0.gri",xinc=inc,yinc=inc,fformat="irap_binary") + + model.builder.extract_nodes(4,maps_dir) + + from warmth.data import haq87 + model.builder.set_eustatic_sea_level(haq87) + + for i in model.builder.iter_node(): + i.rift=[[182,175]] + # + # set 1D node parameters to be most similar to those in the (later) 3D simulation, for better comparison + # + i.bflux = False + i.adiab = 0.3e-3 + # i.crustRHP = 0.0 + # i.sediments['rhp']=[0,0,0,0,0,0] + + model.simulator.simulate_every = 1 + + # + # set 1D simulation parameters to be most similar to those in the (later) 3D simulation, for better comparison + # + model.parameters.HPdcr = 1e32 # "infinite" decay of crustal HP + model.parameters.bflux = False + model.parameters.tetha = 0 + model.parameters.alphav = 0 + + model.simulator.run(save=True,purge=True) + + + + + from warmth.mesh_model import run_3d + import os + try: + os.mkdir('mesh') + except FileExistsError: + pass + try: + os.mkdir('temp') + except FileExistsError: + pass + mm2 = run_3d(model.builder,model.parameters,start_time=model.parameters.time_start,end_time=0) + print("done simulation") + import numpy as np + hx = model.builder.grid.num_nodes_x // 2 + hy = model.builder.grid.num_nodes_y // 2 + # hx = model.builder.grid.num_nodes_x - 1 - pad + # hy = model.builder.grid.num_nodes_y - 1 - pad + + nn = model.builder.nodes[hy][hx] + dd = nn._depth_out[:,0] + + node_ind = hy*model.builder.grid.num_nodes_x + hx + v_per_n = int(mm2.mesh_vertices.shape[0]/(model.builder.grid.num_nodes_y*model.builder.grid.num_nodes_x)) + + temp_1d = np.nan_to_num(nn.temperature_out[:,0], nan=5.0) + temp_3d_ind = np.array([ np.where([mm2.mesh_reindex==i])[1][0] for i in range(node_ind*v_per_n, (node_ind+1)*v_per_n) ] ) + dd_mesh = mm2.mesh.geometry.x[temp_3d_ind,2] + temp_3d_mesh = mm2.u_n.x.array[temp_3d_ind] + + temp_1d_at_mesh_pos = np.interp(dd_mesh, dd, temp_1d) + dd_subset = np.where(dd_mesh<5000) + + 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(f'Max. absolute error in temperature at 3D mesh vertex positions: {max_abs_error}') + print(f'Max. absolute error at depths < 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." + + + + diff --git a/warmth/forward_modelling.py b/warmth/forward_modelling.py index baef0a3..d461b3b 100644 --- a/warmth/forward_modelling.py +++ b/warmth/forward_modelling.py @@ -1586,11 +1586,7 @@ def calculate_new_temperature(self, mean_porosity_arr, self.current_node.sediments["solidus"].values[sed_idx_arr]) conductivity_sed = self._sediment_conductivity_sekiguchi( mean_porosity_arr, self.current_node.sediments["k_cond"].values[sed_idx_arr],Tsed) - if (self.current_node.X==12150) and (self.current_node.Y==12000): - print("conductivity_sed", conductivity_sed[0:20], conductivity_sed[-20:] ) - # print("conductivity_sed", np.mean(conductivity_sed), np.nanmin(conductivity_sed),np.nanmax(conductivity_sed) ) - # print("conductivity_sed", np.mean(xsed), np.nanmin(xsed),np.nanmax(xsed) ) - print("conductivity_sed", len(conductivity_sed), len(xsed)) + conductivity_crust_lith = self._build_crust_lithosphere_properties( coord_crust_lith, hc, hLith, self.current_node.kCrust, self.current_node.kLith, self.current_node.kAsth) diff --git a/warmth/mesh_model.py b/warmth/mesh_model.py index 4dfd636..c9c87ca 100644 --- a/warmth/mesh_model.py +++ b/warmth/mesh_model.py @@ -1,15 +1,19 @@ +from typing import Tuple import numpy as np from mpi4py import MPI import meshio +from pyparsing import Literal import dolfinx from petsc4py import PETSc import ufl from scipy.interpolate import LinearNDInterpolator -from warmth.build import single_node +from warmth.build import single_node, Builder +from warmth.forward_modelling import Forward_model +from .parameters import Parameters from .model import Model 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,NodeGrid +from .mesh_utils import top_crust,top_sed,thick_crust, top_lith, top_asth, top_sed_id, bottom_sed_id,interpolate_all_nodes from .resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties,read_mesh_resqml_hexa def tic(): #Homemade version of matlab tic and toc functions @@ -36,12 +40,14 @@ class UniformNodeGridFixedSizeMeshModel: point_domain_edge_map = {} point_top_vertex_map = {} point_bottom_vertex_map = {} - def __init__(self, model:Model, modelName="test", sedimentsOnly = False): - self.node1D = [n for n in model.builder.iter_node()] + def __init__(self, builder:Builder,parameters:Parameters, sedimentsOnly = False): + self._builder = builder + self._parameters=parameters + self.node1D = [n for n in self._builder.iter_node()] self.num_nodes = len(self.node1D) self.mesh = None - self.modelName = modelName + self.modelName = self._parameters.name self.Temp0 = 5 self.TempBase = 1369 self.verbose = True @@ -55,8 +61,8 @@ def __init__(self, model:Model, modelName="test", sedimentsOnly = False): self.numElemInAsth = 0 if self.runSedimentsOnly else 2 # split asth hexahedron into pieces - self.num_nodes_x = model.builder.grid.num_nodes_x - self.num_nodes_y = model.builder.grid.num_nodes_y + self.num_nodes_x = self._builder.grid.num_nodes_x + self.num_nodes_y = self._builder.grid.num_nodes_y self.convexHullEdges = [] for i in range(self.num_nodes_x-1): edge = [i, i+1] @@ -79,7 +85,7 @@ def __init__(self, model:Model, modelName="test", sedimentsOnly = False): self.thermalCond = None self.mean_porosity = None self.c_rho = None - self.numberOfSediments = model.builder.input_horizons.shape[0]-1 #skip basement + self.numberOfSediments = self._builder.input_horizons.shape[0]-1 #skip basement self.numberOfSedimentCells = self.numberOfSediments * self.numElemPerSediment self.interpolators = {} @@ -260,14 +266,14 @@ def getTopOfAsthAtNode(self, tti, node:single_node): return z0 # - def getSedimentPropForLayerID(self, property, layer_id, node_index): + def getSedimentPropForLayerID(self, property, layer_id:int, node_index:int) ->float: """ """ 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() @@ -276,7 +282,7 @@ def getSedimentPropForLayerID(self, property, layer_id, node_index): # breakpoint() # phi = self.globalSediments[property][layer_id] # assert abs(phi-phi0)<1e-6 - return phi + return prop if (layer_id<=-1) and (layer_id>=-3): lid = -layer_id -1 node = self.node1D[node_index] @@ -294,7 +300,7 @@ def getSedimentPropForLayerID(self, property, layer_id, node_index): return [node.crustliquid,node.lithliquid,node.asthliquid][lid] # liquid density for crust, lith, aest return np.nan - def porosity0ForLayerID(self, layer_id, node_index): + def porosity0ForLayerID(self, layer_id:int, node_index:int)->Tuple[float, float]: """Porosity (at surface) conductivity value for the given layer index """ if (layer_id==-1): @@ -311,37 +317,24 @@ def porosity0ForLayerID(self, layer_id, node_index): return phi, decay return 0.0, 0.0 - def cRhoForLayerID(self, ss, node_index): - # - # prefactor 1000 is the heat capacity.. assumed constant - # + def cRhoForLayerID(self, ss:int, node_index:int)->float: node = self.node1D[node_index] if (ss==-1): - return 1000*node.crustsolid + return self._parameters.cp*node.crustsolid if (ss==-2): - return 1000*node.lithsolid + return self._parameters.cp*node.lithsolid if (ss==-3): - return 1000*node.asthsolid + return self._parameters.cp*node.asthsolid if (ss>=0) and (ssfloat: """Thermal conductivity for a layer ID index """ - # ind0 = cell_id - # ind1 = self.cell_index[cell_id] - # ind2 = np.where(self.cell_index==cell_id)[0][0] - # if self.cell_data_layerID[ind2] != ss: - # print("ind", cell_id, ind0, ind1, ind2) - # print(len(self.cell_index), np.amin(self.cell_index), np.amax(self.cell_index)) - # print(len(self.cell_data_layerID), np.amin(self.cell_data_layerID), np.amax(self.cell_data_layerID)) - # print("kForLayerID", ss, ind0, ind1, ind2, self.cell_data_layerID[self.cell_index[cell_id]] ) - # assert self.cell_data_layerID[ind2] == ss if (node_index > len(self.node1D)-1): - print("cell ID", node_index, len(self.node1D)) - breakpoint() + raise Exception(f"Node index {node_index} larger then node length {len(self.node1D)}") node = self.node1D[node_index] if (ss==-1): return node.kCrust @@ -356,17 +349,13 @@ def kForLayerID(self, ss, node_index): return kc - def rhpForLayerID(self, ss, node_index): + def rhpForLayerID(self, ss:int, node_index:int)->float: """Radiogenic heat production for a layer ID index """ if (ss==-1): node = self.node1D[node_index] kc = node.crustRHP * node._upperCrust_ratio return kc - elif (ss==-2): - return 0 - elif (ss==-3): - return 0 elif (ss>=0) and (ss0) else 3600*24*365 * 5000000 - num_steps = no_steps + num_steps = n_steps # # solver setup, see: @@ -932,15 +917,16 @@ def setupSolverAndSolve(self, time_step=-1, no_steps=100, skip_setup = False, in # source = self.globalSediments.rhp[self.numberOfSediments-1] * 1e-6 # conversion from uW/m^3 # f = dolfinx.fem.Constant(self.mesh, PETSc.ScalarType(source)) # source term - f = self.rhpFcn # * 1e-6 # conversion from uW/m^3 - print("mean RHP", np.mean(self.rhpFcn.x.array[:])) + f = self.rhpFcn + logger.info(f"mean RHP {np.mean(self.rhpFcn.x.array[:])}") + if ( self.useBaseFlux ): # baseFlux = 0.03 if (self.tti>50) else 0.03 baseFlux = self.baseFluxMagnitude # define Neumann condition: constant flux at base # expression g defines values of Neumann BC (heat flux at base) - x = ufl.SpatialCoordinate(self.mesh) + #x = ufl.SpatialCoordinate(self.mesh) domain_c = dolfinx.fem.Function(self.V) if (self.CGorder>1): def marker(x): @@ -948,15 +934,17 @@ def marker(x): return x[2,:]>3990 facets = dolfinx.mesh.locate_entities_boundary(self.mesh, dim=(self.mesh.topology.dim - 2), marker=marker ) - print(type(facets), facets.shape) + #print(type(facets), facets.shape) dofs = dolfinx.fem.locate_dofs_topological(V=self.V, entity_dim=1, entities=facets) - print( type(dofs), len(dofs)) - print(facets.shape, dofs.shape) + #print( type(dofs), len(dofs)) + #print(facets.shape, dofs.shape) if (len(facets)>0): - print( np.amax(facets)) + #print( np.amax(facets)) + pass if (len(dofs)>0): - print( np.amax(dofs)) - print(type(domain_c.x.array), len(domain_c.x.array)) + #print( np.amax(dofs)) + pass + #print(type(domain_c.x.array), len(domain_c.x.array)) domain_c.x.array[ dofs ] = 1 else: basepos = self.getBaseAtMultiplePos(self.mesh.geometry.x[:,0], self.mesh.geometry.x[:,1]) @@ -1298,7 +1286,7 @@ def interpolateResult(self, x): print("PING V V", point) def boundary(x): return np.full(x.shape[1], True) - entities = dolfinx.mesh.locate_entities(self.mesh, 3, boundary ) + #entities = dolfinx.mesh.locate_entities(self.mesh, 3, boundary ) breakpoint() points_on_proc.append(point) points_cells.append(cell_candidates.links(i)[0]) @@ -1330,11 +1318,11 @@ def pointIsOnDomainEdge(self, pt, node0, node1, weight): return False -def run( model:Model, run_simulation=True, start_time=182, end_time=0, out_dir = "out-mapA/"): - +def run_3d( builder:Builder, parameters:Parameters, run_simulation=True, start_time=182, end_time=0, out_dir = "out-mapA/"): + builder=interpolate_all_nodes(builder) nums = 4 - dt = 314712e8 / nums + dt = parameters.myr2s / nums # time step is 1/4 of 1Ma mms2 = [] mms_tti = [] @@ -1352,7 +1340,7 @@ def run( model:Model, run_simulation=True, start_time=182, end_time=0, out_dir = rebuild_mesh = (tti==start_time) if rebuild_mesh: print("Rebuild/reload mesh at tti=", tti) - mm2 = UniformNodeGridFixedSizeMeshModel(model, modelName="test"+str(tti)) + mm2 = UniformNodeGridFixedSizeMeshModel(builder, parameters) print("builing") mm2.buildMesh(tti) print("done") @@ -1363,12 +1351,12 @@ def run( model:Model, run_simulation=True, start_time=182, end_time=0, out_dir = print("===",tti,"=========== ") if ( len(mms2) == 0): tic() - mm2.setupSolverAndSolve(no_steps=40, time_step = 314712e8 * 2e2, skip_setup=False) + mm2.setupSolverAndSolve(n_steps=40, time_step = 314712e8 * 2e2, skip_setup=False) time_solve = time_solve + toc(msg="setup solver and solve") else: tic() # mm2.setupSolverAndSolve( initial_state_model=mms2[-1], no_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) - mm2.setupSolverAndSolve( no_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) + mm2.setupSolverAndSolve( n_steps=nums, time_step=dt, skip_setup=(not rebuild_mesh)) time_solve = time_solve + toc(msg="setup solver and solve") # subvolumes.append(mm2.evaluateVolumes()) if (writeout): diff --git a/warmth/mesh_utils.py b/warmth/mesh_utils.py index 53b8159..ded61f8 100644 --- a/warmth/mesh_utils.py +++ b/warmth/mesh_utils.py @@ -1,6 +1,8 @@ from dataclasses import dataclass - -from warmth.build import single_node +from typing import List +import itertools +import numpy as np +from .build import Builder, single_node @dataclass class NodeGrid: @@ -110,3 +112,81 @@ def volumeOfTet(points): cd = points[2]-points[3] bdcd = np.cross(bd,cd) return np.linalg.norm(np.dot(ad,bdcd))/6 + +def interpolateNode(interpolationNodes: List[single_node], interpolationWeights=None) -> single_node: + assert len(interpolationNodes)>0 + if interpolationWeights is None: + interpolationWeights = np.ones([len(interpolationNodes),1]) + assert len(interpolationNodes)==len(interpolationWeights) + wsum = np.sum(np.array(interpolationWeights)) + iWeightNorm = [ w/wsum for w in interpolationWeights] + + node = single_node() + node.__dict__.update(interpolationNodes[0].__dict__) + node.X = np.sum( np.array( [node.X * w for node,w in zip(interpolationNodes,iWeightNorm)] ) ) + node.Y = np.sum( np.array( [node.Y * w for node,w in zip(interpolationNodes,iWeightNorm)] ) ) + + times = range(node.result._depth.shape[1]) + if node.subsidence is None: + node.subsidence = np.sum( np.array( [ [node.result.seabed(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.crust_ls is None: + node.crust_ls = np.sum( np.array( [ [node.result.crust_thickness(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.lith_ls is None: + node.lith_ls = np.sum( np.array( [ [node.result.lithosphere_thickness(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + + if node.beta is None: + node.beta = np.sum( np.array( [node.beta * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.kAsth is None: + node.kAsth = np.sum( np.array( [node.kAsth * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node.kLith is None: + node.kLith = np.sum( np.array( [node.kLith * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) + if node._depth_out is None: + node._depth_out = np.sum([node.result._depth_out*w for n,w in zip(interpolationNodes[0:1], [1] )], axis=0) + if node.temperature_out is None: + node.temperature_out = np.sum([n.result.temperature_out*w for n,w in zip(interpolationNodes[0:1], [1] )], axis=0) + + if node.sed is None: + node.sed = np.sum([n.sed*w for n,w in zip(interpolationNodes,iWeightNorm)], axis=0) + if node.sed_thickness_ls is None: + node.sed_thickness_ls = node.sed[-1,1,:] - node.sed[0,0,:] + return node + + +def interpolate_all_nodes(builder:Builder)->Builder: + for ni in range(len(builder.nodes)): + for nj in range(len(builder.nodes[ni])): + if builder.nodes[ni][nj] is False: + closest_x_up = [] + for j in range(ni,len(builder.nodes[nj])): + matching_x = [ i[0] for i in builder.indexer_full_sim if i[0]==j ] + closest_x_up = closest_x_up + list(set(matching_x)) + if len(matching_x)>0: + break + closest_x_down = [] + for j in range(ni-1,-1,-1): + matching_x = [ i[0] for i in builder.indexer_full_sim if i[0]==j ] + closest_x_down = closest_x_down + list(set(matching_x)) + if len(matching_x)>0: + break + closest_y_up = [] + for j in range(nj,len(builder.nodes[ni])): + matching_y = [ i[1] for i in builder.indexer_full_sim if (i[1]==j and ((i[0] in closest_x_up) or i[0] in closest_x_down)) ] + closest_y_up = closest_y_up + list(set(matching_y)) + if len(matching_y)>0: + break + closest_y_down = [] + for j in range(nj-1,-1,-1): + matching_y = [ i[1] for i in builder.indexer_full_sim if (i[1]==j and (i[0] in closest_x_up or i[0] in closest_x_down) ) ] + closest_y_down = closest_y_down + list(set(matching_y)) + if len(matching_y)>0: + break + + interpolationNodes = [ builder.nodes[i[0]][i[1]] for i in itertools.product(closest_x_up+closest_x_down, closest_y_up+closest_y_down) ] + interpolationNodes = [nn for nn in interpolationNodes if nn is not False] + node = interpolateNode(interpolationNodes) + node.X, node.Y = builder.grid.location_grid[ni,nj,:] + builder.nodes[ni][nj] = node + else: + node = interpolateNode([builder.nodes[ni][nj]]) # "interpolate" the node from itself to make sure the same member variables exist at the end + builder.nodes[ni][nj] = node + return builder \ No newline at end of file diff --git a/warmth/parameters.py b/warmth/parameters.py index b9dd4f7..5813750 100644 --- a/warmth/parameters.py +++ b/warmth/parameters.py @@ -45,7 +45,7 @@ def __init__(self) -> None: self.maxContLith: float = 130000.0 self.starting_beta: float = 1.1 self.positive_down = True - + self.name:str = "model" pass @property