Merge pull request #5 from equinor/subsheat3D/test_from_grids

Subsheat3 d/test from grids

docs/notebooks/3D_simulation.ipynb

@@ -149,7 +149,9 @@
    "outputs": [],
    "source": [
     "for i in model.builder.iter_node():\n",
-    "    i.rift=[[182,175]]"
+    "    i.rift=[[182,175]]\n",
+    "    i.bflux = False\n",
+    "    i.adiab = 0.3e-3    "
    ]
   },
   {
@@ -158,7 +160,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "model.simulator.simulate_every = 1"
+    "model.simulator.simulate_every = 1\n",
+    "model.parameters.HPdcr = 1e11  # set to \"infinite\" decay\n",
+    "model.parameters.bflux = False  # set to \"infinite\" decay\n"
    ]
   },
   {
@@ -230,7 +234,7 @@
     "# 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",
+    "from warmth3d.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",
@@ -269,8 +273,7 @@
     "            model.builder.nodes[ni][nj] = node\n",
     "        # if (model.builder.nodes[ni][nj].Y>12000):\n",
     "        #     model.builder.nodes[ni][nj].crustRHP = (2e-6) * 4\n",
-    "        #     model.builder.nodes[ni][nj].rhp = (2e-6) * 4\n",
-    "        model.builder.nodes[ni][nj].crustRHP = (2e0) * 1\n"
+    "        #     model.builder.nodes[ni][nj].rhp = (2e-6) * 4\n"
    ]
   },
   {
@@ -308,7 +311,7 @@
     "    os.mkdir('temp')\n",
     "except FileExistsError:\n",
     "    pass\n",
-    "run(model,start_time=model.parameters.time_start,end_time=0)\n",
+    "mm2 = run(model,start_time=model.parameters.time_start,end_time=0)\n",
     "\n"
    ]
   },
@@ -317,12 +320,51 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "import numpy as np\n",
+    "hx = model.builder.grid.num_nodes_x // 2\n",
+    "hy = model.builder.grid.num_nodes_y // 2\n",
+    "# hx = model.builder.grid.num_nodes_x - 1 - pad\n",
+    "# hy = model.builder.grid.num_nodes_y - 1 - pad\n",
+    "\n",
+    "nn = model.builder.nodes[hy][hx]\n",
+    "dd = nn._depth_out[:,0]\n",
+    "\n",
+    "node_ind = hy*model.builder.grid.num_nodes_x + hx\n",
+    "v_per_n = int(mm2.mesh_vertices.shape[0]/(model.builder.grid.num_nodes_y*model.builder.grid.num_nodes_x))\n",
+    "\n",
+    "temp_1d = np.nan_to_num(nn.temperature_out[:,0], nan=5.0)\n",
+    "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) ] )\n",
+    "dd_mesh = mm2.mesh.geometry.x[temp_3d_ind,2]\n",
+    "temp_3d_mesh = mm2.u_n.x.array[temp_3d_ind]\n",
+    "\n",
+    "temp_1d_at_mesh_pos = np.interp(dd_mesh, dd, temp_1d)\n",
+    "dd_subset = np.where(dd_mesh<5000)\n",
+    "print(f'Max. absolute error in temperature at 3D mesh vertex positions: {np.amax(np.abs(temp_1d_at_mesh_pos-temp_3d_mesh))}')\n",
+    "print(f'Max. absolute error at depths < 5000m: {np.amax(np.abs(temp_1d_at_mesh_pos[dd_subset]-temp_3d_mesh[dd_subset]))}')\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib\n",
+    "import matplotlib as plt\n",
+    "plt.use('TkAgg')\n",
+    "\n",
+    "plt.pyplot.plot( dd, temp_1d, label=f'1D simulation'); \n",
+    "# plt.pyplot.plot( dd, temp_3d, label=f'3D simulation'); \n",
+    "plt.pyplot.plot( dd_mesh, temp_3d_mesh, 'o', label=f'3D simulation (nodes)'); \n",
+    "plt.pyplot.legend(loc=\"lower right\", prop={'size': 7})\n",
+    "plt.pyplot.show()\n"
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": ".venv",
+   "display_name": "Python 3.9.10 64-bit",
    "language": "python",
    "name": "python3"
   },
@@ -336,7 +378,12 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.9.10"
+  },
+  "vscode": {
+   "interpreter": {
+    "hash": "e51a2b4e99f2ffc4c1b139fefabaf83a9a6f27184f9cac248baef6121453d47a"
+   }
   }
  },
  "nbformat": 4,

tests/benchmark/README.txt

@@ -1,7 +1,7 @@
 To run the benchmark:
 
 python parallel-1Dsed-hypothetical.py
-In Linux with dolfinx installed: python subsheat3D_hypothetical_extended.py
+In Linux with dolfinx installed: python warmth3d_hypothetical_extended.py
 python plot-comparisons.py
 
 

tests/benchmark/parallel-1Dsed-hypothetical.py

@@ -5,9 +5,9 @@
 
 import warmth
 from warmth.data import haq87
-from subsheat3D.SedimentStack import SedimentStack
+from warmth3d.SedimentStack import SedimentStack
 
-from subsheat3D.Helpers import NodeGrid
+from warmth3d.Helpers import NodeGrid
 
 # logger = get_logger(__name__)
 

tests/benchmark/temperer3D_hypothetical_extended.py

@@ -1,9 +1,9 @@
 import numpy as np
 import pickle
 import itertools
-from subsheat3D.fixed_mesh_model import UniformNodeGridFixedSizeMeshModel
-from subsheat3D.Helpers import NodeGrid
-from subsheat3D.resqpy_helpers import read_mesh_resqml
+from warmth3d.fixed_mesh_model import UniformNodeGridFixedSizeMeshModel
+from warmth3d.Helpers import NodeGrid
+from warmth3d.resqpy_helpers import read_mesh_resqml
 
 def tic():
     #Homemade version of matlab tic and toc functions
@@ -279,7 +279,7 @@ def run( nodeGrid, run_simulation=True, start_time=182, end_time=0, out_dir = "o
 
 
 #
-# NOTE: to compute the required 1D node solutions, you must first run  subsheat3D/parallel-1Dsed.py using the same NodeGrid parameters as below!
+# NOTE: to compute the required 1D node solutions, you must first run  warmth3d/parallel-1Dsed.py using the same NodeGrid parameters as below!
 #
 
 # ng = NodeGrid(0, 0, 97, 97, 2, 2, 1, 100, 100)

tests/temperer3D_mapA_example.py

@@ -1,9 +1,9 @@
 import numpy as np
 import pickle
 import itertools
-from subsheat3D.fixed_mesh_model import UniformNodeGridFixedSizeMeshModel
-from subsheat3D.Helpers import NodeGrid
-from subsheat3D.resqpy_helpers import read_mesh_resqml
+from warmth3d.fixed_mesh_model import UniformNodeGridFixedSizeMeshModel
+from warmth3d.Helpers import NodeGrid
+from warmth3d.resqpy_helpers import read_mesh_resqml
 
 #
 # a map for how many meters were eroded between 100-120 My
@@ -192,7 +192,7 @@ def run( nodeGrid, run_simulation=True, start_time=182, end_time=0, out_dir = "o
 
 
 #
-# NOTE: to compute the required 1D node solutions, you must first run  subsheat3D/parallel-1Dsed.py using the same NodeGrid parameters as below!
+# NOTE: to compute the required 1D node solutions, you must first run  warmth3d/parallel-1Dsed.py using the same NodeGrid parameters as below!
 #
 ng = NodeGrid(150, 0, 485, 548, 500, 1000, 5, 1000, 1000)
 ng = NodeGrid(0, 0, 11, 11, 500, 1000, 5, 100, 100)