From 542da2e12750922a489a385c971a169fd14809b6 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sat, 7 Sep 2024 09:39:12 +1000 Subject: [PATCH] Add runner script for Convection Notebook because memory leaks make for shorter runs and several restarts, unfortunately --- Notebooks/Convection_Cartesian.py | 159 ++- .../NB_Convection_Visualisation-Cart.ipynb | 1031 ++++++++++++----- Notebooks/mpi_runner.sh | 42 + .../Lectures-1.5-TheLithosphere.reveal.qmd | 2 +- 4 files changed, 902 insertions(+), 332 deletions(-) create mode 100755 Notebooks/mpi_runner.sh diff --git a/Notebooks/Convection_Cartesian.py b/Notebooks/Convection_Cartesian.py index 1ef395f..2a58f50 100644 --- a/Notebooks/Convection_Cartesian.py +++ b/Notebooks/Convection_Cartesian.py @@ -3,30 +3,38 @@ # # + -import mpi4py import underworld3 as uw +import os +import numpy as np +import sympy + from underworld3.systems import Stokes from underworld3 import function +import mpi4py -import numpy as np -import sympy -import os +# - # + -# The problem setup +# The problem setup # mesh parameters -width = 6 +width = 1 + +visc_exp = 4 + # Parameters that define the notebook # These can be set when launching the script as # mpirun python3 scriptname -uw_resolution=0.1 etc -ra_expt = uw.options.getReal("ra_expt", default=4) -resolution = uw.options.getInt("resolution", default=8) +ra_expt = uw.options.getReal("ra_expt", default=6) +visc_expt = uw.options.getReal("visc_expt", default=0) +width = uw.options.getInt("width", default=1) +resolution = uw.options.getInt("resolution", default=15) +resolution_in = uw.options.getInt("resolution_in", default=-1) max_steps = uw.options.getInt("max_steps", default=201) restart_step = uw.options.getInt("restart_step", default=-1) expt_desc = uw.options.getString("expt_description", default="") @@ -34,16 +42,21 @@ if expt_desc != "": expt_desc += "_" -print(f"Restarting from step {restart_step}") +if uw.mpi.rank == 0: + print(f"Restarting from step {restart_step}") + +if resolution_in == -1: + resolution_in = resolution # How that works -rayleigh_number = uw.function.expression(r"\textrm{Ra}", - pow(10,ra_expt), # / (r_o-r_i)**3 , - "Rayleigh number") +rayleigh_number = uw.function.expression( + r"\textrm{Ra}", pow(10, ra_expt), "Rayleigh number" # / (r_o-r_i)**3 , +) -expt_name = f"{expt_desc}Ra1e{ra_expt}_res{resolution}" -output_dir = os.path.join("output","cartesian",f"Ra1e{ra_expt}") +old_expt_name = f"{expt_desc}Ra1e{ra_expt}_visc{visc_exp}_res{resolution_in}" +expt_name = f"{expt_desc}Ra1e{ra_expt}_visc{visc_exp}_res{resolution}" +output_dir = os.path.join("output", f"cartesian_{width}x1", f"Ra1e{ra_expt}") os.makedirs(output_dir, exist_ok=True) @@ -51,19 +64,16 @@ # + ## Set up the mesh geometry / discretisation - meshbox = uw.meshing.UnstructuredSimplexBox( - cellSize=1/resolution, - minCoords=(0.0,0.0), - maxCoords=(width,1.0), - degree=1, + cellSize=1 / resolution, + minCoords=(0.0, 0.0), + maxCoords=(width, 1.0), + degree=1, qdegree=3, + regular=False, ) -# meshbox.return_coords_to_bounds = None -meshbox.dm.view() - -x,y = meshbox.CoordinateSystem.X +x, y = meshbox.CoordinateSystem.X x_vector = meshbox.CoordinateSystem.unit_e_0 y_vector = meshbox.CoordinateSystem.unit_e_1 @@ -81,32 +91,41 @@ # Create solver to solver the momentum equation (Stokes flow) stokes = Stokes( - meshbox, - velocityField=v_soln, - pressureField=p_soln, + meshbox, + velocityField=v_soln, + pressureField=p_soln, solver_name="stokes", ) +C = uw.function.expression("C", sympy.log(sympy.Pow(10, sympy.sympify(visc_exp)))) +visc_fn = uw.function.expression( + r"\eta", + sympy.exp(-C.sym * t_soln.sym[0]) * sympy.Pow(10, sympy.sympify(visc_exp)), + "1", +) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel -stokes.constitutive_model.Parameters.viscosity = 1.0 +stokes.constitutive_model.Parameters.shear_viscosity_0 = visc_fn -stokes.tolerance = 0.00001 +stokes.tolerance = 1e-6 +stokes.penalty = 0.0 -penalty = max(1000000, 10*rayleigh_number.sym) +# penalty = max(1000000, 10*rayleigh_number.sym) # Prevent flow crossing the boundaries -# stokes.add_natural_bc(penalty * v_soln.sym[1] * y_vector, "Top") -# stokes.add_natural_bc(penalty * v_soln.sym[1] * y_vector, "Bottom") -# stokes.add_natural_bc(penalty * v_soln.sym[0] * x_vector, "Left") -# stokes.add_natural_bc(penalty * v_soln.sym[0] * x_vector, "Right") - stokes.add_essential_bc((None, 0.0), "Top") stokes.add_essential_bc((None, 0.0), "Bottom") stokes.add_essential_bc((0.0, None), "Left") stokes.add_essential_bc((0.0, None), "Right") stokes.bodyforce = y_vector * rayleigh_number * t_soln.sym[0] +# - + +stokes.constitutive_model.Parameters.shear_viscosity_0 + + +visc_fn._expr_count # + # Create solver for the energy equation (Advection-Diffusion of temperature) @@ -116,7 +135,8 @@ u_Field=t_soln, V_fn=v_soln, solver_name="adv_diff", - order=2, + order=1, + verbose=False, ) adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel @@ -124,42 +144,54 @@ ## Boundary conditions for this solver -adv_diff.add_dirichlet_bc(1.0, "Bottom") -adv_diff.add_dirichlet_bc(0.0, "Top") - +adv_diff.add_dirichlet_bc(+1.0, "Bottom") +adv_diff.add_dirichlet_bc(-0.0, "Top") # + # The advection / diffusion equation is an initial value problem -# We set this up to have +# We set this up with an approximation to the ultimate boundary +# layer structure (you need to provide delta, the TBL thickness) +# +# Add some perturbation and try to offset this on the different boundary +# layers to avoid too much symmetry -init_t = 0.25 * sympy.cos(8.0 * sympy.pi * x) * sympy.sin(sympy.pi * y) + (1-y) +delta = 0.2 +aveT = 0.5 - 0.5 * (sympy.tanh(2 * y / delta) - sympy.tanh(2 * (1 - y) / delta)) + +init_t = ( + 0.0 * sympy.cos(0.5 * sympy.pi * x) * sympy.sin(2 * sympy.pi * y) + + 0.02 * sympy.cos(10.0 * sympy.pi * x) * sympy.sin(2 * sympy.pi * y) + + aveT +) with meshbox.access(t_soln): t_soln.data[...] = uw.function.evaluate(init_t, t_soln.coords).reshape(-1, 1) if restart_step != -1: - print(f"Reading step {restart_step}") - t_soln.read_timestep(data_filename=f"{expt_name}", - data_name="T", - index=restart_step, - outputPath=output_dir) + print(f"Reading step {restart_step} at resolution {resolution_in}") + t_soln.read_timestep( + data_filename=f"{old_expt_name}", + data_name="T", + index=restart_step, + outputPath=output_dir, + ) # + # Solution strategy: solve for the velocity field, then update the temperature field. -# First we check this works for a tiny timestep, later we will repeat this to +# First we check this works for a tiny timestep, later we will repeat this to # model the passage of time. stokes.solve() -adv_diff.solve(timestep=0.00001 * stokes.estimate_dt()) +# adv_diff.solve(timestep=0.00001 * stokes.estimate_dt()) # + if restart_step == -1: timestep = 0 else: timestep = restart_step - -elapsed_time=0.0 + +elapsed_time = 0.0 # + # Convection model / update in time @@ -167,27 +199,28 @@ output = os.path.join(output_dir, expt_name) for step in range(0, max_steps): - - stokes.solve(zero_init_guess=True) - with meshbox.access(v_soln): - v_soln.data[...] = 0.0 - - delta_t = adv_diff.estimate_dt() + + stokes.solve(zero_init_guess=False) + + delta_t = 2.0 * adv_diff.estimate_dt() + delta_ta = stokes.estimate_dt() + adv_diff.solve(timestep=delta_t) # stats then loop tstats = t_soln.stats() if uw.mpi.rank == 0: - print(f"Timestep {timestep}, dt {delta_t}, t {elapsed_time}") + print( + f"Timestep {timestep}, dt {delta_t:.4e}, dta {delta_ta:.4e}, t {elapsed_time:.4e} " + ) - meshbox.write_timestep(filename=f"{expt_name}", - index=timestep, - outputPath=output_dir, - meshVars=[v_soln, p_soln, t_soln], - ) + meshbox.write_timestep( + filename=f"{expt_name}", + index=timestep, + outputPath=output_dir, + meshVars=[v_soln, p_soln, t_soln], + ) timestep += 1 elapsed_time += delta_t - - diff --git a/Notebooks/NB_Convection_Visualisation-Cart.ipynb b/Notebooks/NB_Convection_Visualisation-Cart.ipynb index c5d7c8d..b772d14 100644 --- a/Notebooks/NB_Convection_Visualisation-Cart.ipynb +++ b/Notebooks/NB_Convection_Visualisation-Cart.ipynb @@ -5,20 +5,20 @@ "id": "4d23436d", "metadata": {}, "source": [ - "## Visualise stokes annulus model (flow etc)" + "## Visualise stokes cartesian model (flow etc)" ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 1, "id": "f56963d3", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T01:46:03.387178Z", - "iopub.status.busy": "2024-08-22T01:46:03.386917Z", - "iopub.status.idle": "2024-08-22T01:46:03.400340Z", - "shell.execute_reply": "2024-08-22T01:46:03.391101Z", - "shell.execute_reply.started": "2024-08-22T01:46:03.387152Z" + "iopub.execute_input": "2024-09-06T20:37:02.738720Z", + "iopub.status.busy": "2024-09-06T20:37:02.738516Z", + "iopub.status.idle": "2024-09-06T20:37:02.749244Z", + "shell.execute_reply": "2024-09-06T20:37:02.748849Z", + "shell.execute_reply.started": "2024-09-06T20:37:02.738692Z" } }, "outputs": [], @@ -31,18 +31,26 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 2, "id": "0b026768", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T01:46:03.565791Z", - "iopub.status.busy": "2024-08-22T01:46:03.565556Z", - "iopub.status.idle": "2024-08-22T01:46:03.572286Z", - "shell.execute_reply": "2024-08-22T01:46:03.571313Z", - "shell.execute_reply.started": "2024-08-22T01:46:03.565776Z" + "iopub.execute_input": "2024-09-06T20:37:02.750255Z", + "iopub.status.busy": "2024-09-06T20:37:02.749971Z", + "iopub.status.idle": "2024-09-06T20:37:03.594261Z", + "shell.execute_reply": "2024-09-06T20:37:03.593900Z", + "shell.execute_reply.started": "2024-09-06T20:37:02.750234Z" } }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[Lyrebird.local:16520] shmem: mmap: an error occurred while determining whether or not /var/folders/tx/95gr762j29z4tt5d1dnqlgth0000gn/T//ompi.Lyrebird.501/jf.0/2560229376/sm_segment.Lyrebird.501.989a0000.0 could be created.\n" + ] + } + ], "source": [ "import petsc4py\n", "import underworld3 as uw\n", @@ -52,15 +60,15 @@ }, { "cell_type": "code", - "execution_count": 44, + "execution_count": 3, "id": "8c3f6ccb", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T02:12:08.005644Z", - "iopub.status.busy": "2024-08-22T02:12:08.000042Z", - "iopub.status.idle": "2024-08-22T02:12:08.164006Z", - "shell.execute_reply": "2024-08-22T02:12:08.163065Z", - "shell.execute_reply.started": "2024-08-22T02:12:08.005172Z" + "iopub.execute_input": "2024-09-06T20:37:03.594882Z", + "iopub.status.busy": "2024-09-06T20:37:03.594764Z", + "iopub.status.idle": "2024-09-06T20:37:03.728107Z", + "shell.execute_reply": "2024-09-06T20:37:03.727634Z", + "shell.execute_reply.started": "2024-09-06T20:37:03.594872Z" } }, "outputs": [ @@ -68,31 +76,42 @@ "name": "stdout", "output_type": "stream", "text": [ + "output/cartesian:\n", "total 0\n", - "drwxr-xr-x 1585 lmoresi staff 50720 Aug 20 07:40 \u001b[34mRa1e7.0\u001b[m\u001b[m/\n", - "drwxr-xr-x@ 1208 lmoresi staff 38656 Aug 20 10:19 \u001b[34mRa1e5.0\u001b[m\u001b[m/\n", - "drwxr-xr-x@ 1815 lmoresi staff 58080 Aug 20 18:03 \u001b[34mRa1e6.0\u001b[m\u001b[m/\n", - "drwxr-xr-x@ 808 lmoresi staff 25856 Aug 20 20:12 \u001b[34mRa1e4.0\u001b[m\u001b[m/\n", - "drwxr-xr-x@ 808 lmoresi staff 25856 Aug 20 21:04 \u001b[34mRa1e3.0\u001b[m\u001b[m/\n", - "drwxr-xr-x@ 404 lmoresi staff 12928 Aug 22 11:54 \u001b[34mRa1e2.0\u001b[m\u001b[m/\n" + "drwxr-xr-x@ 294 lmoresi staff 9408 Aug 26 21:59 \u001b[34mRa1e3.0\u001b[m\u001b[m/\n", + "drwxr-xr-x@ 909 lmoresi staff 29088 Aug 27 16:34 \u001b[34mRa1e4.0\u001b[m\u001b[m/\n", + "drwxr-xr-x@ 2009 lmoresi staff 64288 Aug 27 18:41 \u001b[34mRa1e5.0\u001b[m\u001b[m/\n", + "drwxr-xr-x 7517 lmoresi staff 240544 Aug 28 21:19 \u001b[34mRa1e7.0\u001b[m\u001b[m/\n", + "drwxr-xr-x 6758 lmoresi staff 216256 Aug 29 11:44 \u001b[34mRa1e8.0\u001b[m\u001b[m/\n", + "drwxr-xr-x@ 6818 lmoresi staff 218176 Aug 29 11:56 \u001b[34mRa1e6.0\u001b[m\u001b[m/\n", + "\n", + "output/cartesian_2x1:\n", + "total 0\n", + "drwxr-xr-x 6259 lmoresi staff 200288 Aug 31 10:40 \u001b[34mRa1e6.0\u001b[m\u001b[m/\n", + "\n", + "output/cartesian_1x1:\n", + "total 0\n", + "drwxr-xr-x@ 7022 lmoresi staff 224704 Sep 6 19:25 \u001b[34mRa1e6.0\u001b[m\u001b[m/\n", + "drwxr-xr-x@ 5516 lmoresi staff 176512 Sep 6 23:26 \u001b[34mRa1e8.0\u001b[m\u001b[m/\n", + "drwxr-xr-x@ 16529 lmoresi staff 528928 Sep 7 05:42 \u001b[34mRa1e7.0\u001b[m\u001b[m/\n" ] } ], "source": [ - "ls -trl output/cartesian/ | tail" + "ls -trl output/cartesian* " ] }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 4, "id": "bf972017", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T02:12:08.493267Z", - "iopub.status.busy": "2024-08-22T02:12:08.492857Z", - "iopub.status.idle": "2024-08-22T02:12:08.646581Z", - "shell.execute_reply": "2024-08-22T02:12:08.645702Z", - "shell.execute_reply.started": "2024-08-22T02:12:08.493232Z" + "iopub.execute_input": "2024-09-06T20:37:03.728985Z", + "iopub.status.busy": "2024-09-06T20:37:03.728838Z", + "iopub.status.idle": "2024-09-06T20:37:04.050335Z", + "shell.execute_reply": "2024-09-06T20:37:04.049826Z", + "shell.execute_reply.started": "2024-09-06T20:37:03.728968Z" } }, "outputs": [ @@ -100,34 +119,42 @@ "name": "stdout", "output_type": "stream", "text": [ - "-rw-r--r--@ 1 lmoresi staff 234920 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.T.00097.h5\n", - "-rw-r--r--@ 1 lmoresi staff 3076 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.00097.xdmf\n", - "-rw-r--r--@ 1 lmoresi staff 137272 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.U.00098.h5\n", - "-rw-r--r--@ 1 lmoresi staff 51880 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.P.00098.h5\n", - "-rw-r--r--@ 1 lmoresi staff 234920 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.T.00098.h5\n", - "-rw-r--r--@ 1 lmoresi staff 3076 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.00098.xdmf\n", - "-rw-r--r--@ 1 lmoresi staff 137272 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.U.00099.h5\n", - "-rw-r--r--@ 1 lmoresi staff 51880 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.P.00099.h5\n", - "-rw-r--r--@ 1 lmoresi staff 234920 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.T.00099.h5\n", - "-rw-r--r--@ 1 lmoresi staff 3076 Aug 22 11:54 cartesian_Ra1e2.0_res10.mesh.00099.xdmf\n" + "-rw-r--r--@ 1 lmoresi staff 61432 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.U.00299.h5\n", + "-rw-r--r--@ 1 lmoresi staff 28360 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.P.00299.h5\n", + "-rw-r--r--@ 1 lmoresi staff 101480 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.T.00299.h5\n", + "-rw-r--r--@ 1 lmoresi staff 101480 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.W_14_0.00299.h5\n", + "-rw-r--r--@ 1 lmoresi staff 3751 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.00299.xdmf\n", + "-rw-r--r--@ 1 lmoresi staff 61432 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.U.00300.h5\n", + "-rw-r--r--@ 1 lmoresi staff 28360 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.P.00300.h5\n", + "-rw-r--r--@ 1 lmoresi staff 101480 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.T.00300.h5\n", + "-rw-r--r--@ 1 lmoresi staff 101480 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.W_14_0.00300.h5\n", + "-rw-r--r--@ 1 lmoresi staff 3751 Sep 6 23:26 Ra1e8.0_visc4_res15.mesh.00300.xdmf\n" ] } ], "source": [ - "ls -trl output/cartesian/Ra1e2.0/ | tail" + "ls -trl output/cartesian_1x1/Ra1e8.0/ | tail" ] }, { "cell_type": "code", - "execution_count": 46, + "execution_count": null, + "id": "8da0c742-72d5-426c-96cc-517d6f0d1227", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 20, "id": "9a912364", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T02:12:11.699429Z", - "iopub.status.busy": "2024-08-22T02:12:11.699019Z", - "iopub.status.idle": "2024-08-22T02:12:11.709827Z", - "shell.execute_reply": "2024-08-22T02:12:11.709024Z", - "shell.execute_reply.started": "2024-08-22T02:12:11.699390Z" + "iopub.execute_input": "2024-09-06T21:54:48.583673Z", + "iopub.status.busy": "2024-09-06T21:54:48.583249Z", + "iopub.status.idle": "2024-09-06T21:54:48.586866Z", + "shell.execute_reply": "2024-09-06T21:54:48.586554Z", + "shell.execute_reply.started": "2024-09-06T21:54:48.583650Z" }, "lines_to_next_cell": 2 }, @@ -135,9 +162,9 @@ "source": [ "# This is a bit repetitive, sorry \n", "\n", - "checkpoint_dir = f\"output/cartesian/Ra1e2.0\"\n", - "name = \"cartesian_Ra1e2.0\"\n", - "resolution = 10\n", + "checkpoint_dir = \"output/cartesian_1x1/Ra1e8.0\"\n", + "name = \"Ra1e8.0_visc4\"\n", + "resolution = 25\n", "\n", "checkpoint_base = f\"{name}_res{resolution}\"\n", "meshfile = os.path.join(checkpoint_dir, checkpoint_base) + \".mesh.00000.h5\"\n", @@ -147,15 +174,15 @@ }, { "cell_type": "code", - "execution_count": 47, + "execution_count": 21, "id": "b5d63def", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T02:12:12.100297Z", - "iopub.status.busy": "2024-08-22T02:12:12.099804Z", - "iopub.status.idle": "2024-08-22T02:12:12.194200Z", - "shell.execute_reply": "2024-08-22T02:12:12.193898Z", - "shell.execute_reply.started": "2024-08-22T02:12:12.100266Z" + "iopub.execute_input": "2024-09-06T21:54:49.226309Z", + "iopub.status.busy": "2024-09-06T21:54:49.226157Z", + "iopub.status.idle": "2024-09-06T21:54:49.304204Z", + "shell.execute_reply": "2024-09-06T21:54:49.302779Z", + "shell.execute_reply.started": "2024-09-06T21:54:49.226298Z" } }, "outputs": [], @@ -168,124 +195,523 @@ "\n", "v_soln = uw.discretisation.MeshVariable(\"U\", mesh, mesh.dim, degree=2)\n", "t_soln = uw.discretisation.MeshVariable(\"T\", mesh, 1, degree=3)\n", - "flux = uw.discretisation.MeshVariable(r\"dTdz\", mesh, 1, degree=1)" + "flux = uw.discretisation.MeshVariable(r\"dTdz\", mesh, 1, degree=1, continuous=False)" ] }, { "cell_type": "code", - "execution_count": 48, + "execution_count": 28, "id": "942cc670-9f6c-4469-9e89-c614f66b9135", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T02:12:12.537975Z", - "iopub.status.busy": "2024-08-22T02:12:12.537619Z", - "iopub.status.idle": "2024-08-22T02:12:12.542459Z", - "shell.execute_reply": "2024-08-22T02:12:12.541713Z", - "shell.execute_reply.started": "2024-08-22T02:12:12.537949Z" + "iopub.execute_input": "2024-09-06T22:10:15.267987Z", + "iopub.status.busy": "2024-09-06T22:10:15.267576Z", + "iopub.status.idle": "2024-09-06T22:10:15.282273Z", + "shell.execute_reply": "2024-09-06T22:10:15.280979Z", + "shell.execute_reply.started": "2024-09-06T22:10:15.267958Z" + } + }, + "outputs": [], + "source": [ + "flux_solver = uw.systems.Projection(mesh, flux)\n", + "\n", + "flux_fn = -mesh.vector.gradient(t_soln.sym[0]).dot(mesh.CoordinateSystem.unit_e_1) # + v_soln.sym[1] * t_soln.sym[0]\n", + "flux_solver.uw_function = flux_fn\n", + "flux_solver.smoothing = 0.0\n", + "\n", + "\n", + "sample_points = np.zeros((100,2))\n", + "sample_points[:,0] = np.linspace(0,1,100)\n", + "sample_points[:,1] = 1\n", + "\n", + "sample_points_surf = sample_points.copy()\n", + "\n", + "sample_points[:,1] = 0\n", + "sample_points_base = sample_points.copy()\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "f4ee9fc9-ce74-4e06-89a6-338987d702c9", + "metadata": { + "execution": { + "iopub.execute_input": "2024-09-06T22:10:17.490501Z", + "iopub.status.busy": "2024-09-06T22:10:17.490350Z", + "iopub.status.idle": "2024-09-06T22:10:17.501631Z", + "shell.execute_reply": "2024-09-06T22:10:17.501066Z", + "shell.execute_reply.started": "2024-09-06T22:10:17.490490Z" } }, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "DM Object: uw_output/cartesian/Ra1e2.0/cartesian_Ra1e2.0_res10.mesh.00000.h5 1 MPI process\n", - " type: plex\n", - "uw_output/cartesian/Ra1e2.0/cartesian_Ra1e2.0_res10.mesh.00000.h5 in 2 dimensions:\n", - " Number of 0-cells per rank: 793\n", - " Number of 1-cells per rank: 2236\n", - " Number of 2-cells per rank: 1444\n", - "Labels:\n", - " depth: 3 strata with value/size (0 (793), 1 (2236), 2 (1444))\n", - " All_Boundaries: 1 strata with value/size (1001 (140))\n", - " Bottom: 1 strata with value/size (11 (119))\n", - " Elements: 1 strata with value/size (99999 (2097))\n", - " Left: 1 strata with value/size (14 (19))\n", - " Right: 1 strata with value/size (13 (19))\n", - " Top: 1 strata with value/size (12 (119))\n", - " UW_Boundaries: 5 strata with value/size (1001 (140), 11 (119), 12 (119), 13 (19), 14 (19))\n", - " celltype: 3 strata with value/size (0 (793), 1 (2236), 3 (1444))\n", - "Field U:\n", - " adjacency FEM\n", - "Field T:\n", - " adjacency FEM\n", - "Field dTdz:\n", - " adjacency FEM\n" - ] + "data": { + "text/markdown": [ + "**Class**: " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/markdown": [ + "**Poisson system solver**" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/markdown": [ + "Primary problem: " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\tiny \\quad \\nabla \\cdot \\color{Blue}\\left[\\begin{matrix}0 & 0\\end{matrix}\\right]$ + " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/latex": [ + "$\\tiny \\phantom{ \\quad \\nabla \\cdot} \\color{DarkRed}\\left[\\begin{matrix}{ { \\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\! } {T} }_{,1}(\\mathbf{x}) + { { \\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\!\\,\\! } {dTdz} }(\\mathbf{x})\\end{matrix}\\right]\\color{Black} = 0 $" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/markdown": [ + "This solver is formulated in 2 dimensions" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ - "mesh.dm.view()" + "flux_solver.view()" ] }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 44, "id": "b9768720", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T02:12:13.010028Z", - "iopub.status.busy": "2024-08-22T02:12:13.009384Z", - "iopub.status.idle": "2024-08-22T02:12:13.014551Z", - "shell.execute_reply": "2024-08-22T02:12:13.013669Z", - "shell.execute_reply.started": "2024-08-22T02:12:13.009988Z" + "iopub.execute_input": "2024-09-06T23:27:54.209965Z", + "iopub.status.busy": "2024-09-06T23:27:54.208661Z", + "iopub.status.idle": "2024-09-06T23:27:54.233567Z", + "shell.execute_reply": "2024-09-06T23:27:54.232773Z", + "shell.execute_reply.started": "2024-09-06T23:27:54.209937Z" } }, "outputs": [], "source": [ - "steps = range(0,30,1)" + "steps = range(1450,1800, 1)\n", + "calc_flux = True\n", + "add_streamlines = False\n", + "cmap = \"RdBu_r\"" ] }, { "cell_type": "code", - "execution_count": 50, + "execution_count": 45, "id": "25ffe942", "metadata": { "execution": { - "iopub.execute_input": "2024-08-22T02:12:13.500584Z", - "iopub.status.busy": "2024-08-22T02:12:13.500238Z", - "iopub.status.idle": "2024-08-22T02:12:18.148882Z", - "shell.execute_reply": "2024-08-22T02:12:18.148572Z", - "shell.execute_reply.started": "2024-08-22T02:12:13.500556Z" + "iopub.execute_input": "2024-09-06T23:27:54.692849Z", + "iopub.status.busy": "2024-09-06T23:27:54.692636Z", + "iopub.status.idle": "2024-09-06T23:29:28.947607Z", + "shell.execute_reply": "2024-09-06T23:29:28.946725Z", + "shell.execute_reply.started": "2024-09-06T23:27:54.692836Z" }, - "lines_to_next_cell": 0, - "scrolled": true + "lines_to_next_cell": 0 }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Plotting step 0\n", - "Plotting step 1\n", - "Plotting step 2\n", - "Plotting step 3\n", - "Plotting step 4\n", - "Plotting step 5\n", - "Plotting step 6\n", - "Plotting step 7\n", - "Plotting step 8\n", - "Plotting step 9\n", - "Plotting step 10\n", - "Plotting step 11\n", - "Plotting step 12\n", - "Plotting step 13\n", - "Plotting step 14\n", - "Plotting step 15\n", - "Plotting step 16\n", - "Plotting step 17\n", - "Plotting step 18\n", - "Plotting step 19\n", - "Plotting step 20\n", - "Plotting step 21\n", - "Plotting step 22\n", - "Plotting step 23\n", - "Plotting step 24\n", - "Plotting step 25\n", - "Plotting step 26\n", - "Plotting step 27\n", - "Plotting step 28\n", - "Plotting step 29\n" + "01450, 7.3276, 13.6851, 10.5063\n", + "01451, 7.3353, 13.6430, 10.4891\n", + "01452, 7.3432, 13.6123, 10.4778\n", + "01453, 7.3486, 13.5973, 10.4729\n", + "01454, 7.3421, 13.6011, 10.4716\n", + "01455, 7.3297, 13.6033, 10.4665\n", + "01456, 7.3159, 13.5921, 10.4540\n", + "01457, 7.3034, 13.5625, 10.4330\n", + "01458, 7.2932, 13.5192, 10.4062\n", + "01459, 7.2835, 13.4734, 10.3785\n", + "01460, 7.2782, 13.4263, 10.3523\n", + "01461, 7.2785, 13.3781, 10.3283\n", + "01462, 7.2853, 13.3278, 10.3066\n", + "01463, 7.2984, 13.2783, 10.2883\n", + "01464, 7.3059, 13.2471, 10.2765\n", + "01465, 7.3093, 13.2254, 10.2673\n", + "01466, 7.3100, 13.2122, 10.2611\n", + "01467, 7.3044, 13.2121, 10.2583\n", + "01468, 7.2951, 13.2182, 10.2567\n", + "01469, 7.2819, 13.2262, 10.2540\n", + "01470, 7.2697, 13.2288, 10.2492\n", + "01471, 7.2605, 13.2212, 10.2408\n", + "01472, 7.2553, 13.1988, 10.2270\n", + "01473, 7.2518, 13.1712, 10.2115\n", + "01474, 7.2499, 13.1478, 10.1988\n", + "01475, 7.2468, 13.1307, 10.1888\n", + "01476, 7.2419, 13.1211, 10.1815\n", + "01477, 7.2378, 13.1219, 10.1799\n", + "01478, 7.2372, 13.1368, 10.1870\n", + "01479, 7.2377, 13.1772, 10.2075\n", + "01480, 7.2349, 13.2484, 10.2416\n", + "01481, 7.2270, 13.3396, 10.2833\n", + "01482, 7.2087, 13.4357, 10.3222\n", + "01483, 7.1867, 13.5207, 10.3537\n", + "01484, 7.1670, 13.5896, 10.3783\n", + "01485, 7.1524, 13.6421, 10.3973\n", + "01486, 7.1455, 13.6807, 10.4131\n", + "01487, 7.1419, 13.7072, 10.4245\n", + "01488, 7.1396, 13.7228, 10.4312\n", + "01489, 7.1357, 13.7321, 10.4339\n", + "01490, 7.1331, 13.7407, 10.4369\n", + "01491, 7.1295, 13.7577, 10.4436\n", + "01492, 7.1227, 13.7756, 10.4491\n", + "01493, 7.1099, 13.7645, 10.4372\n", + "01494, 7.0997, 13.7100, 10.4048\n", + "01495, 7.0923, 13.6237, 10.3580\n", + "01496, 7.0883, 13.5211, 10.3047\n", + "01497, 7.0876, 13.4335, 10.2605\n", + "01498, 7.0896, 13.3725, 10.2310\n", + "01499, 7.0884, 13.3477, 10.2181\n", + "01500, 7.0863, 13.3044, 10.1954\n", + "01501, 7.0856, 13.2305, 10.1581\n", + "01502, 7.0814, 13.0857, 10.0836\n", + "01503, 7.0781, 13.0412, 10.0597\n", + "01504, 7.0743, 12.9900, 10.0321\n", + "01505, 7.0699, 12.9026, 9.9863\n", + "01506, 7.0653, 12.8068, 9.9361\n", + "01507, 7.0606, 12.7215, 9.8910\n", + "01508, 7.0561, 12.6638, 9.8599\n", + "01509, 7.0517, 12.6432, 9.8475\n", + "01510, 7.0478, 12.5984, 9.8231\n", + "01511, 7.0443, 12.5407, 9.7925\n", + "01512, 7.0411, 12.5173, 9.7792\n", + "01513, 7.0384, 12.5530, 9.7957\n", + "01514, 7.0361, 12.5817, 9.8089\n", + "01515, 7.0342, 12.7380, 9.8861\n", + "01516, 7.0244, 13.0747, 10.0496\n", + "01517, 7.0107, 13.4399, 10.2253\n", + "01518, 6.9952, 13.8069, 10.4011\n", + "01519, 6.9532, 14.1255, 10.5393\n", + "01520, 6.9033, 14.3812, 10.6422\n", + "01521, 6.8553, 14.5738, 10.7146\n", + "01522, 6.8213, 14.7211, 10.7712\n", + "01523, 6.7983, 14.7952, 10.7967\n", + "01524, 6.7897, 14.8864, 10.8380\n", + "01525, 6.7765, 14.9690, 10.8728\n", + "01526, 6.7580, 15.0762, 10.9171\n", + "01527, 6.7307, 15.1853, 10.9580\n", + "01528, 6.6986, 15.2353, 10.9669\n", + "01529, 6.6663, 15.2153, 10.9408\n", + "01530, 6.6344, 15.1640, 10.8992\n", + "01531, 6.6040, 15.0860, 10.8450\n", + "01532, 6.5748, 14.9865, 10.7806\n", + "01533, 6.5509, 14.8833, 10.7171\n", + "01534, 6.5343, 14.7468, 10.6405\n", + "01535, 6.5245, 14.5893, 10.5569\n", + "01536, 6.5207, 14.4410, 10.4809\n", + "01537, 6.5258, 14.2827, 10.4042\n", + "01538, 6.5378, 14.0892, 10.3135\n", + "01539, 6.5555, 13.8836, 10.2196\n", + "01540, 6.5814, 13.7005, 10.1409\n", + "01541, 6.6204, 13.5625, 10.0914\n", + "01542, 6.6563, 13.4796, 10.0680\n", + "01543, 6.6851, 13.4542, 10.0696\n", + "01544, 6.7078, 13.4611, 10.0844\n", + "01545, 6.7254, 13.5200, 10.1227\n", + "01546, 6.7389, 13.6669, 10.2029\n", + "01547, 6.7490, 13.8630, 10.3060\n", + "01548, 6.7565, 14.0646, 10.4105\n", + "01549, 6.7615, 14.2125, 10.4870\n", + "01550, 6.7648, 14.2930, 10.5289\n", + "01551, 6.7665, 14.2707, 10.5186\n", + "01552, 6.7499, 14.2525, 10.5012\n", + "01553, 6.7277, 14.2693, 10.4985\n", + "01554, 6.7041, 14.1092, 10.4066\n", + "01555, 6.6791, 13.9228, 10.3009\n", + "01556, 6.6515, 13.6948, 10.1731\n", + "01557, 6.6218, 13.5708, 10.0963\n", + "01558, 6.5929, 13.4543, 10.0236\n", + "01559, 6.5672, 13.2886, 9.9279\n", + "01560, 6.5410, 13.1359, 9.8385\n", + "01561, 6.5150, 13.0516, 9.7833\n", + "01562, 6.4909, 13.0436, 9.7672\n", + "01563, 6.4694, 13.1042, 9.7868\n", + "01564, 6.4470, 13.2270, 9.8370\n", + "01565, 6.4272, 13.3353, 9.8812\n", + "01566, 6.4100, 13.4542, 9.9321\n", + "01567, 6.3953, 13.5814, 9.9883\n", + "01568, 6.3866, 13.6897, 10.0381\n", + "01569, 6.3801, 13.7783, 10.0792\n", + "01570, 6.3703, 13.8482, 10.1093\n", + "01571, 6.3586, 13.8941, 10.1263\n", + "01572, 6.3457, 13.8587, 10.1022\n", + "01573, 6.3324, 13.7857, 10.0591\n", + "01574, 6.3184, 13.7256, 10.0220\n", + "01575, 6.3041, 13.7071, 10.0056\n", + "01576, 6.2903, 13.6550, 9.9726\n", + "01577, 6.2762, 13.6253, 9.9508\n", + "01578, 6.2612, 13.6469, 9.9540\n", + "01579, 6.2465, 13.6556, 9.9511\n", + "01580, 6.2334, 13.6880, 9.9607\n", + "01581, 6.2194, 13.7989, 10.0091\n", + "01582, 6.2017, 13.9559, 10.0788\n", + "01583, 6.1810, 14.0949, 10.1379\n", + "01584, 6.1473, 14.2213, 10.1843\n", + "01585, 6.0988, 14.3144, 10.2066\n", + "01586, 6.0444, 14.3632, 10.2038\n", + "01587, 5.9865, 14.3743, 10.1804\n", + "01588, 5.9307, 14.3558, 10.1432\n", + "01589, 5.8756, 14.3200, 10.0978\n", + "01590, 5.8223, 14.2704, 10.0463\n", + "01591, 5.7737, 14.2109, 9.9923\n", + "01592, 5.7279, 14.1425, 9.9352\n", + "01593, 5.6888, 14.0667, 9.8777\n", + "01594, 5.6559, 13.9839, 9.8199\n", + "01595, 5.6311, 13.8928, 9.7620\n", + "01596, 5.6171, 13.7927, 9.7049\n", + "01597, 5.6148, 13.6830, 9.6489\n", + "01598, 5.6275, 13.5651, 9.5963\n", + "01599, 5.6534, 13.4379, 9.5457\n", + "01600, 5.6890, 13.3088, 9.4989\n", + "01601, 5.7202, 13.1930, 9.4566\n", + "01602, 5.7434, 13.0926, 9.4180\n", + "01603, 5.7591, 13.0039, 9.3815\n", + "01604, 5.7683, 12.9231, 9.3457\n", + "01605, 5.7719, 12.8464, 9.3091\n", + "01606, 5.7714, 12.7710, 9.2712\n", + "01607, 5.7690, 12.6947, 9.2319\n", + "01608, 5.7655, 12.6170, 9.1913\n", + "01609, 5.7624, 12.5370, 9.1497\n", + "01610, 5.7604, 12.4543, 9.1074\n", + "01611, 5.7585, 12.3692, 9.0638\n", + "01612, 5.7572, 12.2811, 9.0192\n", + "01613, 5.7582, 12.1892, 8.9737\n", + "01614, 5.7617, 12.0937, 8.9277\n", + "01615, 5.7689, 11.9940, 8.8814\n", + "01616, 5.7813, 11.8892, 8.8353\n", + "01617, 5.8003, 11.7789, 8.7896\n", + "01618, 5.8225, 11.6656, 8.7441\n", + "01619, 5.8491, 11.5492, 8.6992\n", + "01620, 5.8744, 11.4355, 8.6550\n", + "01621, 5.8981, 11.3260, 8.6120\n", + "01622, 5.9201, 11.2229, 8.5715\n", + "01623, 5.9405, 11.1286, 8.5346\n", + "01624, 5.9596, 11.0451, 8.5023\n", + "01625, 5.9778, 10.9743, 8.4760\n", + "01626, 5.9957, 10.9176, 8.4566\n", + "01627, 6.0139, 10.8754, 8.4447\n", + "01628, 6.0329, 10.8480, 8.4405\n", + "01629, 6.0530, 10.8368, 8.4449\n", + "01630, 6.0741, 10.8430, 8.4585\n", + "01631, 6.0964, 10.8652, 8.4808\n", + "01632, 6.1191, 10.9043, 8.5117\n", + "01633, 6.1415, 10.9632, 8.5523\n", + "01634, 6.1629, 11.0331, 8.5980\n", + "01635, 6.1829, 11.1116, 8.6473\n", + "01636, 6.2013, 11.2139, 8.7076\n", + "01637, 6.2184, 11.3270, 8.7727\n", + "01638, 6.2345, 11.4218, 8.8281\n", + "01639, 6.2500, 11.5262, 8.8881\n", + "01640, 6.2651, 11.6332, 8.9492\n", + "01641, 6.2801, 11.7078, 8.9939\n", + "01642, 6.2948, 11.7954, 9.0451\n", + "01643, 6.3093, 11.9154, 9.1123\n", + "01644, 6.3236, 11.9737, 9.1486\n", + "01645, 6.3373, 12.0644, 9.2008\n", + "01646, 6.3467, 12.1791, 9.2629\n", + "01647, 6.3503, 12.3836, 9.3670\n", + "01648, 6.3458, 12.7228, 9.5343\n", + "01649, 6.3325, 13.0461, 9.6893\n", + "01650, 6.3023, 13.2997, 9.8010\n", + "01651, 6.2609, 13.5312, 9.8961\n", + "01652, 6.2141, 13.7045, 9.9593\n", + "01653, 6.1655, 13.8618, 10.0137\n", + "01654, 6.1147, 14.0017, 10.0582\n", + "01655, 6.0613, 14.1611, 10.1112\n", + "01656, 6.0085, 14.3277, 10.1681\n", + "01657, 5.9552, 14.4549, 10.2050\n", + "01658, 5.9041, 14.5830, 10.2436\n", + "01659, 5.8556, 14.7579, 10.3068\n", + "01660, 5.8094, 14.9270, 10.3682\n", + "01661, 5.7679, 15.1466, 10.4572\n", + "01662, 5.7291, 15.4172, 10.5731\n", + "01663, 5.6933, 15.6040, 10.6486\n", + "01664, 5.6597, 15.7469, 10.7033\n", + "01665, 5.6280, 15.8980, 10.7630\n", + "01666, 5.5941, 16.0307, 10.8124\n", + "01667, 5.5591, 16.1334, 10.8462\n", + "01668, 5.5229, 16.1974, 10.8601\n", + "01669, 5.4851, 16.2281, 10.8566\n", + "01670, 5.4482, 16.2452, 10.8467\n", + "01671, 5.4133, 16.2505, 10.8319\n", + "01672, 5.3802, 16.2310, 10.8056\n", + "01673, 5.3494, 16.1771, 10.7633\n", + "01674, 5.3242, 16.1015, 10.7129\n", + "01675, 5.3048, 16.0177, 10.6613\n", + "01676, 5.2907, 15.9102, 10.6005\n", + "01677, 5.2856, 15.7714, 10.5285\n", + "01678, 5.2938, 15.5962, 10.4450\n", + "01679, 5.3156, 15.4009, 10.3583\n", + "01680, 5.3582, 15.2005, 10.2793\n", + "01681, 5.4290, 14.9714, 10.2002\n", + "01682, 5.5290, 14.6593, 10.0942\n", + "01683, 5.6237, 14.3481, 9.9859\n", + "01684, 5.7133, 14.0882, 9.9008\n", + "01685, 5.7957, 13.8706, 9.8331\n", + "01686, 5.8756, 13.6722, 9.7739\n", + "01687, 5.9451, 13.5047, 9.7249\n", + "01688, 6.0129, 13.3614, 9.6871\n", + "01689, 6.0717, 13.2757, 9.6737\n", + "01690, 6.1230, 13.2532, 9.6881\n", + "01691, 6.1679, 13.2790, 9.7235\n", + "01692, 6.2078, 13.3393, 9.7736\n", + "01693, 6.2436, 13.4100, 9.8268\n", + "01694, 6.2650, 13.4809, 9.8729\n", + "01695, 6.2788, 13.5367, 9.9077\n", + "01696, 6.2939, 13.5782, 9.9361\n", + "01697, 6.3028, 13.6297, 9.9662\n", + "01698, 6.3062, 13.6959, 10.0010\n", + "01699, 6.3081, 13.7631, 10.0356\n", + "01700, 6.3072, 13.8259, 10.0666\n", + "01701, 6.3053, 13.8777, 10.0915\n", + "01702, 6.3039, 13.9151, 10.1095\n", + "01703, 6.3061, 13.9937, 10.1499\n", + "01704, 6.3055, 14.0607, 10.1831\n", + "01705, 6.3044, 14.1473, 10.2258\n", + "01706, 6.3062, 14.2328, 10.2695\n", + "01707, 6.3104, 14.2968, 10.3036\n", + "01708, 6.3168, 14.3336, 10.3252\n", + "01709, 6.3212, 14.3319, 10.3265\n", + "01710, 6.3268, 14.3001, 10.3134\n", + "01711, 6.3307, 14.2719, 10.3013\n", + "01712, 6.3333, 14.2581, 10.2957\n", + "01713, 6.3356, 14.2528, 10.2942\n", + "01714, 6.3389, 14.2313, 10.2851\n", + "01715, 6.3424, 14.1732, 10.2578\n", + "01716, 6.3455, 14.1055, 10.2255\n", + "01717, 6.3489, 14.0502, 10.1995\n", + "01718, 6.3525, 13.9990, 10.1757\n", + "01719, 6.3564, 13.9043, 10.1303\n", + "01720, 6.3599, 13.7855, 10.0727\n", + "01721, 6.3618, 13.6801, 10.0209\n", + "01722, 6.3633, 13.5948, 9.9790\n", + "01723, 6.3643, 13.4952, 9.9298\n", + "01724, 6.3643, 13.3597, 9.8620\n", + "01725, 6.3622, 13.2464, 9.8043\n", + "01726, 6.3588, 13.1675, 9.7632\n", + "01727, 6.3567, 13.0665, 9.7116\n", + "01728, 6.3543, 13.0176, 9.6859\n", + "01729, 6.3486, 13.0683, 9.7084\n", + "01730, 6.3384, 13.1390, 9.7387\n", + "01731, 6.3230, 13.1970, 9.7600\n", + "01732, 6.2858, 13.2520, 9.7689\n", + "01733, 6.2379, 13.2805, 9.7592\n", + "01734, 6.1877, 13.2798, 9.7338\n", + "01735, 6.1358, 13.2577, 9.6968\n", + "01736, 6.0875, 13.2221, 9.6548\n", + "01737, 6.0413, 13.1793, 9.6103\n", + "01738, 5.9987, 13.1299, 9.5643\n", + "01739, 5.9609, 13.0745, 9.5177\n", + "01740, 5.9286, 13.0134, 9.4710\n", + "01741, 5.9026, 12.9468, 9.4247\n", + "01742, 5.8825, 12.8756, 9.3790\n", + "01743, 5.8694, 12.7990, 9.3342\n", + "01744, 5.8540, 12.7233, 9.2886\n", + "01745, 5.8327, 12.6495, 9.2411\n", + "01746, 5.8072, 12.5763, 9.1918\n", + "01747, 5.7807, 12.5025, 9.1416\n", + "01748, 5.7504, 12.4292, 9.0898\n", + "01749, 5.7213, 12.3547, 9.0380\n", + "01750, 5.6933, 12.2798, 8.9865\n", + "01751, 5.6676, 12.2038, 8.9357\n", + "01752, 5.6499, 12.1245, 8.8872\n", + "01753, 5.6403, 12.0435, 8.8419\n", + "01754, 5.6428, 11.9588, 8.8008\n", + "01755, 5.6637, 11.8717, 8.7677\n", + "01756, 5.6967, 11.7820, 8.7393\n", + "01757, 5.7296, 11.6957, 8.7127\n", + "01758, 5.7625, 11.6125, 8.6875\n", + "01759, 5.7950, 11.5321, 8.6636\n", + "01760, 5.8272, 11.4538, 8.6405\n", + "01761, 5.8561, 11.3775, 8.6168\n", + "01762, 5.8794, 11.3024, 8.5909\n", + "01763, 5.8985, 11.2273, 8.5629\n", + "01764, 5.9140, 11.1520, 8.5330\n", + "01765, 5.9276, 11.0769, 8.5022\n", + "01766, 5.9398, 11.0029, 8.4714\n", + "01767, 5.9507, 10.9316, 8.4411\n", + "01768, 5.9612, 10.8637, 8.4124\n", + "01769, 5.9723, 10.8003, 8.3863\n", + "01770, 5.9860, 10.7422, 8.3641\n", + "01771, 6.0028, 10.6910, 8.3469\n", + "01772, 6.0240, 10.6479, 8.3359\n", + "01773, 6.0478, 10.6152, 8.3315\n", + "01774, 6.0746, 10.5942, 8.3344\n", + "01775, 6.1061, 10.5870, 8.3466\n", + "01776, 6.1348, 10.5964, 8.3656\n", + "01777, 6.1597, 10.6222, 8.3910\n", + "01778, 6.1804, 10.6658, 8.4231\n", + "01779, 6.1973, 10.7289, 8.4631\n", + "01780, 6.2108, 10.8133, 8.5120\n", + "01781, 6.2217, 10.9211, 8.5714\n", + "01782, 6.2306, 11.0556, 8.6431\n", + "01783, 6.2381, 11.2193, 8.7287\n", + "01784, 6.2450, 11.4126, 8.8288\n", + "01785, 6.2519, 11.6315, 8.9417\n", + "01786, 6.2592, 11.8717, 9.0655\n", + "01787, 6.2671, 12.1285, 9.1978\n", + "01788, 6.2750, 12.3879, 9.3314\n", + "01789, 6.2831, 12.6337, 9.4584\n", + "01790, 6.2913, 12.8691, 9.5802\n", + "01791, 6.2995, 13.1040, 9.7017\n", + "01792, 6.3077, 13.3364, 9.8220\n", + "01793, 6.3160, 13.5607, 9.9384\n", + "01794, 6.3246, 13.7717, 10.0481\n", + "01795, 6.3340, 13.9659, 10.1499\n", + "01796, 6.3445, 14.1454, 10.2449\n", + "01797, 6.3565, 14.3126, 10.3345\n", + "01798, 6.3712, 14.4747, 10.4230\n", + "01799, 6.3887, 14.6333, 10.5110\n" ] } ], @@ -294,14 +720,15 @@ "import pyvista as pv\n", "import underworld3.visualisation as vis\n", "\n", - "pl = pv.Plotter(window_size=[1000, 500])\n", + "pl = pv.Plotter(window_size=[750, 750])\n", "\n", "for step in steps:\n", "\n", " try:\n", " v_soln.read_timestep(checkpoint_base, \"U\", step, outputPath=checkpoint_dir)\n", " t_soln.read_timestep(checkpoint_base, \"T\", step, outputPath=checkpoint_dir)\n", - " print(f\"Plotting step {step}\")\n", + " if not calc_flux:\n", + " print(f\"Plotting {step:05d}\")\n", "\n", " except:\n", " print(f\"Failed to read step {step} correctly\")\n", @@ -320,64 +747,86 @@ " pvmesh_t.point_data[\"T\"] = vis.scalar_fn_to_pv_points(pvmesh_t, t_soln.sym[0])\n", " \n", " # point sources at cell centres\n", - " skip = 10\n", + " skip = 19\n", " points = np.zeros((mesh._centroids[::skip].shape[0], 3))\n", " points[:, 0] = mesh._centroids[::skip, 0]\n", " points[:, 1] = mesh._centroids[::skip, 1]\n", " point_cloud = pv.PolyData(points)\n", "\n", - " pvstream = pvmesh.streamlines_from_source(\n", - " point_cloud,\n", - " vectors=\"V\",\n", - " integration_direction=\"both\",\n", - " max_time=1,\n", - " surface_streamlines=True,\n", - " )\n", + "\n", " \n", " pl.add_mesh(\n", " pvmesh_t,\n", - " cmap=\"RdBu_r\",\n", + " cmap=cmap,\n", " edge_color=\"Black\",\n", " edge_opacity=0.1,\n", - " show_edges=True,\n", + " show_edges=False,\n", " scalars=\"T\",\n", " use_transparency=False,\n", " opacity=1.0,\n", " # clim=[0.0,1.0],\n", - " # show_scalar_bar=False\n", + " show_scalar_bar=False\n", " )\n", "\n", "\n", "\n", - " pl.add_mesh(pvmesh, cmap=\"RdBu_r\", scalars=\"T\", opacity=0.5, show_edges=True, edge_color=\"Black\",\n", - " edge_opacity=0.5, show_scalar_bar=True) \n", + " # pl.add_mesh(pvmesh, cmap=\"RdBu_r\", scalars=\"T\", opacity=0.5, show_edges=True, edge_color=\"Black\",\n", + " # edge_opacity=0.5, show_scalar_bar=True) \n", "\n", - " pl.add_mesh(pvstream)\n", + " if add_streamlines:\n", + " pvstream = pvmesh.streamlines_from_source(\n", + " point_cloud,\n", + " vectors=\"V\",\n", + " integration_direction=\"both\",\n", + " max_time=1,\n", + " surface_streamlines=True,\n", + " )\n", + " \n", + " pl.add_mesh(pvstream,\n", + " opacity=0.5,\n", + " show_scalar_bar=False)\n", "\n", - " pl.camera_position = [(3.0, 0.5, 6.750994846396983),\n", - " (3.0, 0.5, 0.0),\n", - " (0.0, 1.0, 0.0)]\n", + " # pl.camera_position = [(3.0, 0.5, 6.750994846396983),\n", + " # (3.0, 0.5, 0.0),\n", + " # (0.0, 1.0, 0.0)]\n", "\n", " imagefile = os.path.join(image_dir, checkpoint_base) + f\"_V_T_{step}.png\"\n", " \n", " pl.screenshot(filename=imagefile,\n", - " window_size=(1000, 500),\n", + " window_size=(1000, 1000),\n", " return_img=False)\n", + "\n", + " if calc_flux:\n", + " flux_solver.solve()\n", + " surf_flux = uw.function.evaluate(flux.sym, sample_points_surf)\n", + " basal_flux = uw.function.evaluate(flux.sym, sample_points_base)\n", + " \n", + " # Heat flux: Upper, Lower, Mean\n", + " print(f\"{step:05d}, {surf_flux.mean():8.4f}, {basal_flux.mean():8.4f}, {((surf_flux.mean() + basal_flux.mean())/2):8.4f}\")\n", + "\n", "\n" ] }, { "cell_type": "code", - "execution_count": 43, + "execution_count": null, + "id": "7e8d7158-a608-41e2-a172-dd73b0383ee4", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 47, "id": "ef5a9c93", "metadata": { "editable": true, "execution": { - "iopub.execute_input": "2024-08-22T01:51:17.880307Z", - "iopub.status.busy": "2024-08-22T01:51:17.880079Z", - "iopub.status.idle": "2024-08-22T01:51:17.970414Z", - "shell.execute_reply": "2024-08-22T01:51:17.961560Z", - "shell.execute_reply.started": "2024-08-22T01:51:17.880290Z" + "iopub.execute_input": "2024-09-06T23:31:45.307179Z", + "iopub.status.busy": "2024-09-06T23:31:45.306732Z", + "iopub.status.idle": "2024-09-06T23:31:45.464296Z", + "shell.execute_reply": "2024-09-06T23:31:45.448349Z", + "shell.execute_reply.started": "2024-09-06T23:31:45.307149Z" }, "slideshow": { "slide_type": "" @@ -385,15 +834,23 @@ "tags": [] }, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "A view with name (P_0x316480f50_8) is already registered\n", + " => returning previous one\n" + ] + }, { "data": { "application/vnd.jupyter.widget-view+json": { - "model_id": "f0832e971e4946f9bf650b0add46fbd7", + "model_id": "9960671e118f465898bbc6b921200138", "version_major": 2, "version_minor": 0 }, "text/plain": [ - "Widget(value='