Skip to content

Commit

Permalink
Merge branch 'main' into release
Browse files Browse the repository at this point in the history
  • Loading branch information
@cianwilson
cianwilson committed Aug 8, 2024
2 parents 73a37e2 + 582a95b commit 2e67747
Show file tree
Hide file tree
Showing 233 changed files with 4,142 additions and 656 deletions.
6 changes: 6 additions & 0 deletions notebooks/sz_geometry.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -459,6 +459,8 @@
"source": [
"if __name__ == \"__main__\":\n",
" plotter_mesh = utils.plot_mesh(mesh, tags=cell_tags, gather=True, show_edges=True, line_width=1)\n",
" utils.plot_geometry(geom, plotter=plotter_mesh, color='green', width=2)\n",
" utils.plot_couplingdepth(slab, plotter=plotter_mesh, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_mesh)\n",
" utils.plot_save(plotter_mesh, output_folder / \"sz_geometry_benchmark_mesh.png\")"
]
Expand Down Expand Up @@ -591,6 +593,8 @@
"source": [
"if __name__ == \"__main__\":\n",
" plotter_mesh_ak = utils.plot_mesh(mesh_ak, tags=cell_tags_ak, gather=True, show_edges=True, line_width=1)\n",
" utils.plot_geometry(geom_ak, plotter=plotter_mesh_ak, color='green', width=2)\n",
" utils.plot_couplingdepth(slab_ak, plotter=plotter_mesh_ak, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_mesh_ak)\n",
" utils.plot_save(plotter_mesh_ak, output_folder / \"sz_geometry_ak_mesh.png\")"
]
Expand Down Expand Up @@ -716,6 +720,8 @@
"source": [
"if __name__ == \"__main__\":\n",
" plotter_mesh_ant = utils.plot_mesh(mesh_ant, tags=cell_tags_ant, gather=True, show_edges=True, line_width=1)\n",
" utils.plot_geometry(geom_ant, plotter=plotter_mesh_ant, color='green', width=2)\n",
" utils.plot_couplingdepth(slab_ant, plotter=plotter_mesh_ant, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_mesh_ant)\n",
" utils.plot_save(plotter_mesh_ant, output_folder / \"sz_geometry_ant_mesh.png\")"
]
Expand Down
6 changes: 6 additions & 0 deletions notebooks/sz_geometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -306,6 +306,8 @@ def create_sz_geometry(slab, resscale, sztype, io_depth, extra_width,

if __name__ == "__main__":
plotter_mesh = utils.plot_mesh(mesh, tags=cell_tags, gather=True, show_edges=True, line_width=1)
utils.plot_geometry(geom, plotter=plotter_mesh, color='green', width=2)
utils.plot_couplingdepth(slab, plotter=plotter_mesh, render_points_as_spheres=True, point_size=10.0, color='green')
utils.plot_show(plotter_mesh)
utils.plot_save(plotter_mesh, output_folder / "sz_geometry_benchmark_mesh.png")

Expand Down Expand Up @@ -376,6 +378,8 @@ def create_sz_geometry(slab, resscale, sztype, io_depth, extra_width,

if __name__ == "__main__":
plotter_mesh_ak = utils.plot_mesh(mesh_ak, tags=cell_tags_ak, gather=True, show_edges=True, line_width=1)
utils.plot_geometry(geom_ak, plotter=plotter_mesh_ak, color='green', width=2)
utils.plot_couplingdepth(slab_ak, plotter=plotter_mesh_ak, render_points_as_spheres=True, point_size=10.0, color='green')
utils.plot_show(plotter_mesh_ak)
utils.plot_save(plotter_mesh_ak, output_folder / "sz_geometry_ak_mesh.png")

Expand Down Expand Up @@ -443,6 +447,8 @@ def create_sz_geometry(slab, resscale, sztype, io_depth, extra_width,

if __name__ == "__main__":
plotter_mesh_ant = utils.plot_mesh(mesh_ant, tags=cell_tags_ant, gather=True, show_edges=True, line_width=1)
utils.plot_geometry(geom_ant, plotter=plotter_mesh_ant, color='green', width=2)
utils.plot_couplingdepth(slab_ant, plotter=plotter_mesh_ant, render_points_as_spheres=True, point_size=10.0, color='green')
utils.plot_show(plotter_mesh_ant)
utils.plot_save(plotter_mesh_ant, output_folder / "sz_geometry_ant_mesh.png")

Expand Down
76 changes: 64 additions & 12 deletions notebooks/sz_problem.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -836,6 +836,8 @@
" self.update_T_functions()\n",
"\n",
" # just set the boundary conditions on the boundaries for the velocities\n",
" self.wedge_vpw_i.x.array[:] = 0.0\n",
" self.slab_vps_i.x.array[:] = 0.0\n",
" df.fem.set_bc(self.wedge_vpw_i.x.array, self.bcs_vpw)\n",
" df.fem.set_bc(self.slab_vps_i.x.array, self.bcs_vps)\n",
" # and update the interpolated v functions for consistency\n",
Expand Down Expand Up @@ -1128,8 +1130,10 @@
"outputs": [],
"source": [
"if __name__ == \"__main__\":\n",
" plotter_ic = utils.plot_scalar(sz_case1.T_i, scale=sz_case1.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)'})\n",
" plotter_ic = utils.plot_scalar(sz_case1.T_i, scale=sz_case1.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})\n",
" utils.plot_vector_glyphs(sz_case1.vw_i, plotter=plotter_ic, gather=True, factor=0.1, color='k', scale=utils.mps_to_mmpyr(sz_case1.v0))\n",
" utils.plot_geometry(sz_case1.geom, plotter=plotter_ic, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case1.geom.slab_spline, plotter=plotter_ic, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_ic)\n",
" utils.plot_save(plotter_ic, output_folder / \"sz_problem_case1_ics.png\")"
]
Expand Down Expand Up @@ -1614,9 +1618,11 @@
"outputs": [],
"source": [
"if __name__ == \"__main__\":\n",
" plotter_iso = utils.plot_scalar(sz_case1.T_i, scale=sz_case1.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)'})\n",
" plotter_iso = utils.plot_scalar(sz_case1.T_i, scale=sz_case1.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})\n",
" utils.plot_vector_glyphs(sz_case1.vw_i, plotter=plotter_iso, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case1.v0))\n",
" utils.plot_vector_glyphs(sz_case1.vs_i, plotter=plotter_iso, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case1.v0))\n",
" utils.plot_geometry(sz_case1.geom, plotter=plotter_iso, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case1.geom.slab_spline, plotter=plotter_iso, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_iso)\n",
" utils.plot_save(plotter_iso, output_folder / \"sz_problem_case1_solution.png\")"
]
Expand Down Expand Up @@ -2071,9 +2077,11 @@
"outputs": [],
"source": [
"if __name__ == \"__main__\":\n",
" plotter_dis = utils.plot_scalar(sz_case2.T_i, scale=sz_case2.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)'})\n",
" plotter_dis = utils.plot_scalar(sz_case2.T_i, scale=sz_case2.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})\n",
" utils.plot_vector_glyphs(sz_case2.vw_i, plotter=plotter_dis, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case2.v0))\n",
" utils.plot_vector_glyphs(sz_case2.vs_i, plotter=plotter_dis, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case2.v0))\n",
" utils.plot_geometry(sz_case2.geom, plotter=plotter_dis, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case2.geom.slab_spline, plotter=plotter_dis, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_dis)\n",
" utils.plot_save(plotter_dis, output_folder / \"sz_problem_case2_solution.png\")"
]
Expand Down Expand Up @@ -2122,7 +2130,9 @@
"source": [
"if __name__ == \"__main__\":\n",
" eta_i = sz_case2.project_dislocationcreep_viscosity()\n",
" plotter_eta = utils.plot_scalar(eta_i, scale=sz_case2.eta0, log_scale=True, show_edges=True, scalar_bar_args={'title': 'Viscosity (Pa) [log scale]'})\n",
" plotter_eta = utils.plot_scalar(eta_i, scale=sz_case2.eta0, log_scale=True, show_edges=True, scalar_bar_args={'title': 'Viscosity (Pa) [log scale]', 'bold':True})\n",
" utils.plot_geometry(sz_case2.geom, plotter=plotter_eta, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case2.geom.slab_spline, plotter=plotter_eta, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_eta)\n",
" utils.plot_save(plotter_eta, output_folder / \"sz_problem_case2_eta.png\")"
]
Expand Down Expand Up @@ -2243,7 +2253,7 @@
"outputs": [],
"source": [
"class SubductionProblem(SubductionProblem):\n",
" def solve_timedependent_isoviscous(self, tf, dt, theta=0.5, verbosity=2, petsc_options=None):\n",
" def solve_timedependent_isoviscous(self, tf, dt, theta=0.5, verbosity=2, petsc_options=None, plotter=None):\n",
" \"\"\"\n",
" Solve the coupled temperature-velocity-pressure problem assuming an isoviscous rheology with time dependency\n",
"\n",
Expand All @@ -2268,6 +2278,9 @@
" self.dt = df.fem.Constant(self.mesh, df.default_scalar_type(dt/self.t0_Myr))\n",
" self.theta = df.fem.Constant(self.mesh, df.default_scalar_type(theta))\n",
"\n",
" # reset the initial conditions\n",
" self.setup_boundaryconditions()\n",
" \n",
" # first solve both Stokes systems\n",
" self.solve_stokes_isoviscous(petsc_options=petsc_options)\n",
"\n",
Expand All @@ -2285,6 +2298,11 @@
" while t < tf_nd - 1e-9:\n",
" if self.comm.rank == 0 and verbosity>1:\n",
" print(\"Step: {:>6d}, Times: {:>9g} -> {:>9g} Myr\".format(ti, t*self.t0_Myr, (t+self.dt.value)*self.t0_Myr))\n",
" if plotter is not None:\n",
" for mesh in plotter.meshes:\n",
" if self.T_i.name in mesh.point_data:\n",
" mesh.point_data[self.T_i.name][:] = self.T_i.x.array\n",
" plotter.write_frame()\n",
" self.T_n.x.array[:] = self.T_i.x.array\n",
" self.T_i = problem_T.solve()\n",
" ti+=1\n",
Expand Down Expand Up @@ -2317,7 +2335,17 @@
" geom_case1td = create_sz_geometry(slab, resscale, sztype, io_depth_1, extra_width, \n",
" coast_distance, lc_depth, uc_depth)\n",
" sz_case1td = SubductionProblem(geom_case1td, A=A, Vs=Vs, sztype=sztype, qs=qs)\n",
" sz_case1td.solve_timedependent_isoviscous(25, 0.05, theta=0.5)"
"\n",
" fps = 5\n",
" plotter_gif = pv.Plotter(notebook=False, off_screen=True)\n",
" utils.plot_scalar(sz_case1td.T_i, plotter=plotter_gif, scale=sz_case1td.T0, gather=True, cmap='coolwarm', clim=[0.0, sz_case1td.Tm*sz_case1td.T0], scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})\n",
" utils.plot_geometry(sz_case1td.geom, plotter=plotter_gif, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case1td.geom.slab_spline, plotter=plotter_gif, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" plotter_gif.open_gif( str(output_folder / \"sz_problem_case1td_solution.gif\"), fps=fps)\n",
" \n",
" sz_case1td.solve_timedependent_isoviscous(25, 0.05, theta=0.5, plotter=plotter_gif)\n",
" \n",
" plotter_gif.close()"
]
},
{
Expand All @@ -2328,10 +2356,13 @@
"outputs": [],
"source": [
"if __name__ == \"__main__\":\n",
" plotter_isotd = utils.plot_scalar(sz_case1td.T_i, scale=sz_case1td.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)'})\n",
" plotter_isotd = utils.plot_scalar(sz_case1td.T_i, scale=sz_case1td.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})\n",
" utils.plot_vector_glyphs(sz_case1td.vw_i, plotter=plotter_isotd, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case1td.v0))\n",
" utils.plot_vector_glyphs(sz_case1td.vs_i, plotter=plotter_isotd, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case1td.v0))\n",
" utils.plot_show(plotter_isotd)"
" utils.plot_geometry(sz_case1td.geom, plotter=plotter_isotd, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case1td.geom.slab_spline, plotter=plotter_isotd, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_isotd)\n",
" utils.plot_save(plotter_isotd, output_folder / \"sz_problem_case1td_solution.png\")"
]
},
{
Expand All @@ -2351,7 +2382,7 @@
"source": [
"class SubductionProblem(SubductionProblem):\n",
" def solve_timedependent_dislocationcreep(self, tf, dt, theta=0.5, rtol=5.e-6, atol=5.e-9, maxits=50, verbosity=2, \n",
" petsc_options=None):\n",
" petsc_options=None, plotter=None):\n",
" \"\"\"\n",
" Solve the coupled temperature-velocity-pressure problem assuming a dislocation creep rheology with time dependency\n",
"\n",
Expand Down Expand Up @@ -2380,6 +2411,9 @@
" self.dt = df.fem.Constant(self.mesh, df.default_scalar_type(dt/self.t0_Myr))\n",
" self.theta = df.fem.Constant(self.mesh, df.default_scalar_type(theta))\n",
" \n",
" # reset the initial conditions\n",
" self.setup_boundaryconditions()\n",
" \n",
" # first solve the isoviscous problem\n",
" self.solve_stokes_isoviscous(petsc_options=petsc_options)\n",
"\n",
Expand Down Expand Up @@ -2432,6 +2466,11 @@
" while t < tf_nd - 1e-9:\n",
" if self.comm.rank == 0 and verbosity>1:\n",
" print(\"Step: {:>6d}, Times: {:>9g} -> {:>9g} Myr\".format(ti, t*self.t0_Myr, (t+self.dt.value)*self.t0_Myr,))\n",
" if plotter is not None:\n",
" for mesh in plotter.meshes:\n",
" if self.T_i.name in mesh.point_data:\n",
" mesh.point_data[self.T_i.name][:] = self.T_i.x.array\n",
" plotter.write_frame()\n",
" self.T_n.x.array[:] = self.T_i.x.array\n",
" # calculate the initial residual\n",
" r = calculate_residual()\n",
Expand Down Expand Up @@ -2497,7 +2536,17 @@
" geom_case2td = create_sz_geometry(slab, resscale, sztype, io_depth_2, extra_width, \n",
" coast_distance, lc_depth, uc_depth)\n",
" sz_case2td = SubductionProblem(geom_case2td, A=A, Vs=Vs, sztype=sztype, qs=qs)\n",
" sz_case2td.solve_timedependent_dislocationcreep(10, 0.05, theta=0.5, rtol=1.e-3)"
" \n",
" fps = 5\n",
" plotter_gif2 = pv.Plotter(notebook=False, off_screen=True)\n",
" utils.plot_scalar(sz_case2td.T_i, plotter=plotter_gif2, scale=sz_case2td.T0, gather=True, cmap='coolwarm', clim=[0.0, sz_case2td.Tm*sz_case2td.T0], scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})\n",
" utils.plot_geometry(sz_case2td.geom, plotter=plotter_gif2, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case2td.geom.slab_spline, plotter=plotter_gif2, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" plotter_gif2.open_gif( str(output_folder / \"sz_problem_case2td_solution.gif\"), fps=fps)\n",
" \n",
" sz_case2td.solve_timedependent_dislocationcreep(10, 0.05, theta=0.5, rtol=1.e-3, plotter=plotter_gif2)\n",
"\n",
" plotter_gif2.close()"
]
},
{
Expand All @@ -2508,10 +2557,13 @@
"outputs": [],
"source": [
"if __name__ == \"__main__\":\n",
" plotter_distd = utils.plot_scalar(sz_case2td.T_i, scale=sz_case2td.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)'})\n",
" plotter_distd = utils.plot_scalar(sz_case2td.T_i, scale=sz_case2td.T0, gather=True, cmap='coolwarm', scalar_bar_args={'title': 'Temperature (deg C)', 'bold':True})\n",
" utils.plot_vector_glyphs(sz_case2td.vw_i, plotter=plotter_distd, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case2td.v0))\n",
" utils.plot_vector_glyphs(sz_case2td.vs_i, plotter=plotter_distd, factor=0.1, gather=True, color='k', scale=utils.mps_to_mmpyr(sz_case2td.v0))\n",
" utils.plot_show(plotter_distd)"
" utils.plot_geometry(sz_case1td.geom, plotter=plotter_distd, color='green', width=2)\n",
" utils.plot_couplingdepth(sz_case1td.geom.slab_spline, plotter=plotter_distd, render_points_as_spheres=True, point_size=10.0, color='green')\n",
" utils.plot_show(plotter_distd)\n",
" utils.plot_save(plotter_distd, output_folder / \"sz_problem_case2td_solution.png\")"
]
},
{
Expand Down
Loading

0 comments on commit 2e67747

Please sign in to comment.