diff --git a/docs/notebooks/3D_simulation.ipynb b/docs/notebooks/3D_simulation.ipynb index f0ab58f..9f93038 100644 --- a/docs/notebooks/3D_simulation.ipynb +++ b/docs/notebooks/3D_simulation.ipynb @@ -220,6 +220,55 @@ "model.simulator.run(save=True,purge=True)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# interpolate and extrapolate the missing nodes\n", + "# find nearby nodes from the array indexer_full_sim, which is sorted by x-index\n", + "import itertools\n", + "from subsheat3D.fixed_mesh_model import interpolateNode\n", + "for ni in range(len(model.builder.nodes)):\n", + " for nj in range(len(model.builder.nodes[ni])):\n", + " if model.builder.nodes[ni][nj] is False:\n", + " closest_x_up = []\n", + " for j in range(ni,len(model.builder.nodes[nj])):\n", + " matching_x = [ i[0] for i in model.builder.indexer_full_sim if i[0]==j ]\n", + " closest_x_up = closest_x_up + list(set(matching_x))\n", + " if len(matching_x)>0:\n", + " break\n", + " closest_x_down = []\n", + " for j in range(ni-1,-1,-1):\n", + " matching_x = [ i[0] for i in model.builder.indexer_full_sim if i[0]==j ]\n", + " closest_x_down = closest_x_down + list(set(matching_x))\n", + " if len(matching_x)>0:\n", + " break\n", + " closest_y_up = []\n", + " for j in range(nj,len(model.builder.nodes[ni])):\n", + " matching_y = [ i[1] for i in model.builder.indexer_full_sim if (i[1]==j and ((i[0] in closest_x_up) or i[0] in closest_x_down)) ]\n", + " closest_y_up = closest_y_up + list(set(matching_y))\n", + " if len(matching_y)>0:\n", + " break\n", + " closest_y_down = []\n", + " for j in range(nj-1,-1,-1):\n", + " matching_y = [ i[1] for i in model.builder.indexer_full_sim if (i[1]==j and (i[0] in closest_x_up or i[0] in closest_x_down) ) ]\n", + " closest_y_down = closest_y_down + list(set(matching_y))\n", + " if len(matching_y)>0:\n", + " break\n", + "\n", + " interpolationNodes = [ model.builder.nodes[i[0]][i[1]] for i in itertools.product(closest_x_up+closest_x_down, closest_y_up+closest_y_down) ]\n", + " interpolationNodes = [nn for nn in interpolationNodes if nn is not False]\n", + " node = interpolateNode(interpolationNodes)\n", + " node.X, node.Y = model.builder.grid.location_grid[ni,nj,:]\n", + " model.builder.nodes[ni][nj] = node\n", + " else:\n", + " node = interpolateNode([model.builder.nodes[ni][nj]]) # \"interpolate\" the node from itself to make sure the same member variables exist at the end\n", + " model.builder.nodes[ni][nj] = node\n" + ] + }, { "cell_type": "code", "execution_count": 15, @@ -246,7 +295,17 @@ ], "source": [ "from warmth.mesh_model import run\n", - "run(model,start_time=model.parameters.time_start,end_time=0)" + "import os\n", + "try:\n", + " os.mkdir('mesh')\n", + "except FileExistsError:\n", + " pass\n", + "try:\n", + " os.mkdir('temp')\n", + "except FileExistsError:\n", + " pass\n", + "run(model,start_time=model.parameters.time_start,end_time=0)\n", + "\n" ] }, { diff --git a/subsheat3D/fixed_mesh_model.py b/subsheat3D/fixed_mesh_model.py index 387f16b..eebdc84 100644 --- a/subsheat3D/fixed_mesh_model.py +++ b/subsheat3D/fixed_mesh_model.py @@ -1,7 +1,9 @@ import numpy as np from mpi4py import MPI from scipy.interpolate import LinearNDInterpolator +from typing import Iterator, List, Literal +from warmth.build import single_node from warmth.logging import logger from subsheat3D.Helpers import NodeParameters1D, top_crust,top_sed,thick_crust, top_lith, top_asth, top_sed_id, bottom_sed_id from subsheat3D.resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties @@ -15,6 +17,46 @@ # NOTE: the default values of surface and base heat flow imply zero RHP in the crust # + +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 + + # class UniformNodeGridFixedSizeMeshModel: """Manages a 3D heat equation computation using dolfinx diff --git a/warmth/mesh_model.py b/warmth/mesh_model.py index 072d36d..ef524b1 100644 --- a/warmth/mesh_model.py +++ b/warmth/mesh_model.py @@ -10,7 +10,7 @@ from .model import Model from warmth.logging import logger from .mesh_utils import NodeParameters1D, top_crust,top_sed,thick_crust, top_lith, top_asth, top_sed_id, bottom_sed_id,NodeGrid -from .resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties,read_mesh_resqml +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 import time @@ -1374,4 +1374,4 @@ def run( model:Model, run_simulation=True, start_time=182, end_time=0, out_dir = print("total time solve: " , time_solve) EPCfilename = mm2.write_hexa_mesh_resqml("temp/") print("RESQML model written to: " , EPCfilename) - read_mesh_resqml(EPCfilename) # test reading of the .epc file \ No newline at end of file + read_mesh_resqml_hexa(EPCfilename) # test reading of the .epc file \ No newline at end of file