-
Notifications
You must be signed in to change notification settings - Fork 61
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge branch 'v2.16.x' of github.com:underworldcode/underworld2 into …
…v2.16.x
- Loading branch information
Showing
12 changed files
with
355 additions
and
79 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
253 changes: 253 additions & 0 deletions
253
docs/UWGeodynamics/examples/1_23_05_FreeSurface_3D_Relaxation.ipynb
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,253 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "attachments": {}, | ||
| "cell_type": "markdown", | ||
| "id": "25381a9f-72ce-4ed9-9d31-5eaf2daa2ae1", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "3D_Relaxation\n", | ||
| "======\n", | ||
| "\n", | ||
| "This notebook reproduce the a cosine perturbation of the surface relaxation with true free surface. Loading of Earth's surface can be described with an initial periodic surface displacement of a viscous fluid within an infinite half space, the solution of which is outlined in Turcotte and Schubert (1982), 6.10 Postglacial Rebound. The surface decreases exponentially with time and is dependent on the magnitude, $w_m$, and wavelength $\\lambda$ of the perturbation, and the viscosity, $\\eta$ and density, $\\rho$ of the fluid,\n", | ||
| "\n", | ||
| "$$ w = w_m exp\\Big(\\frac{-\\lambda \\rho g t}{4\\pi\\eta}\\Big) $$\n", | ||
| "\n", | ||
| "where $w$ is displacement, $w_m$ the initial load magnitude, $g$ gravity, $t$ time. This solution can be charaterised by the relaxation time, $t_{relax} = 4\\pi\\eta / \\rho g \\lambda $, the time taken for the initial load to decrease by $e^{-1}$. The solution for an elastic material with the equivalent load produces the same magnitude of displacement instantaneously.\n", | ||
| "\n", | ||
| "\n", | ||
| "**Keywords:** Free surface\n", | ||
| "\n", | ||
| "<img src=\"./images/3D_FreeSurface.gif\" width = \"400\" height = \"200\" align=center /> \n", | ||
| "<img src=\"./images/3D_FreeSurface_Topography.png\" width = \"400\" height = \"300\" align=center />" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "5f1419bd-076d-4aa8-8093-3f9265f52713", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "from underworld import UWGeodynamics as GEO\n", | ||
| "from underworld import visualisation as vis\n", | ||
| "import underworld as uw\n", | ||
| "import numpy as np" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "bbcb3ae5-34ac-4d7f-879e-b6b4263aaa63", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "u = GEO.UnitRegistry\n", | ||
| "ndim = GEO.non_dimensionalise\n", | ||
| "dim = GEO.dimensionalise\n", | ||
| "\n", | ||
| "# scaling 3: vel\n", | ||
| "half_rate = 1.0 * u.centimeter / u.year\n", | ||
| "model_length = 100. * u.kilometer\n", | ||
| "bodyforce = 3300 * u.kilogram / u.metre**3 * 9.81 * u.meter / u.second**2\n", | ||
| "\n", | ||
| "KL = model_length\n", | ||
| "Kt = KL / half_rate\n", | ||
| "KM = bodyforce * KL**2 * Kt**2\n", | ||
| "\n", | ||
| "GEO.scaling_coefficients[\"[length]\"] = KL\n", | ||
| "GEO.scaling_coefficients[\"[time]\"] = Kt\n", | ||
| "GEO.scaling_coefficients[\"[mass]\"]= KM" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "d00b5e0a-95c3-4540-9024-79f05bc62aad", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "zres = 8\n", | ||
| "xres = zres*4\n", | ||
| "yres = zres*4\n", | ||
| "figsize = (800,400)\n", | ||
| "\n", | ||
| "xmin, xmax = -200 * u.kilometer, 200 * u.kilometer\n", | ||
| "ymin, ymax = -200 * u.kilometer, 200 * u.kilometer\n", | ||
| "zmin, zmax = -100 * u.kilometer, 0 * u.kilometer\n", | ||
| "\n", | ||
| "eta = ndim(1e21 * u.pascal * u.second)\n", | ||
| "density = ndim(3300 * u.kilogram / u.metre**3)\n", | ||
| "gravity = ndim(9.81 * u.meter / u.second**2)\n", | ||
| "w_m = ndim(5.0 * u.kilometer)\n", | ||
| "Lambda = ndim(100.0 * u.kilometer) \n", | ||
| "\n", | ||
| "densityM = density\n", | ||
| "viscM = eta\n", | ||
| "ND_gravity = gravity\n", | ||
| "\n", | ||
| "def perturbation(x):\n", | ||
| " return w_m * np.cos(2.*np.pi*(x)/Lambda)\n", | ||
| "\n", | ||
| "# analytic solution\n", | ||
| "xMax = ndim(xmax - xmin)\n", | ||
| "x = np.linspace(0, xMax, 200+1)\n", | ||
| "w_0 = perturbation(x)\n", | ||
| "t_relax = 4 * np.pi * eta / (Lambda * density * gravity)\n", | ||
| "tMax = t_relax * 5 \n", | ||
| "t = np.linspace(0, tMax, 100 * 10 + 1)\n", | ||
| "w_t = w_m * np.exp(-1.*t/t_relax)\n", | ||
| "\n", | ||
| "max_time = dim(tMax,u.kiloyear)\n", | ||
| "dt_set = dim(t_relax*1e-2,u.kiloyear)\n", | ||
| "save_every = 5\n", | ||
| "checkpoint_interval = dt_set*save_every\n", | ||
| "\n", | ||
| "Model = GEO.Model(elementRes=(xres,yres,zres),\n", | ||
| " minCoord=(xmin,ymin,zmin), \n", | ||
| " maxCoord=(xmax,ymax,zmax),\n", | ||
| " gravity=(0.0, 0.0,-9.81 * u.meter / u.second**2))\n", | ||
| "\n", | ||
| "fdir = \"1_23_05_FreeSurface_3D_Relaxation_zres\"+str(zres)+\"/\"\n", | ||
| "Model.outputDir = fdir" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "60244a41-ea7d-44c4-b5b1-f87ee6ec7690", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "MShape = GEO.shapes.Layer3D(top=Model.top,bottom=Model.bottom)\n", | ||
| "ma = Model.add_material(name=\"material\",shape=MShape)" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "291c22b7-7a3e-41a7-bfb9-d36775cb61d8", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "def perturbation3D(x,y):\n", | ||
| " return w_m * np.cos(2.*np.pi*(x)/Lambda)\n", | ||
| "\n", | ||
| "minCoord = tuple([GEO.nd(val) for val in Model.minCoord])\n", | ||
| "maxCoord = tuple([GEO.nd(val) for val in Model.maxCoord])\n", | ||
| "\n", | ||
| "init_mesh = uw.mesh.FeMesh_Cartesian(elementType=Model.elementType,\n", | ||
| " elementRes=Model.elementRes,\n", | ||
| " minCoord=minCoord,\n", | ||
| " maxCoord=maxCoord,\n", | ||
| " periodic=Model.periodic)\n", | ||
| "\n", | ||
| "TField = init_mesh.add_variable(nodeDofCount=1)\n", | ||
| "TField.data[:, 0] = init_mesh.data[:, -1].copy()\n", | ||
| "\n", | ||
| "top = Model.top_wall\n", | ||
| "bottom = Model.bottom_wall\n", | ||
| "conditions = uw.conditions.DirichletCondition(variable=TField,indexSetsPerDof=(top + bottom,))\n", | ||
| "system = uw.systems.SteadyStateHeat(\n", | ||
| " temperatureField=TField,\n", | ||
| " fn_diffusivity=1.0,\n", | ||
| " conditions=conditions)\n", | ||
| "solver = uw.systems.Solver(system)\n", | ||
| "\n", | ||
| "x = init_mesh.data[top,0]\n", | ||
| "y = init_mesh.data[top,1]\n", | ||
| "TField.data[top, 0] = perturbation3D(x,y)\n", | ||
| "\n", | ||
| "solver.solve()\n", | ||
| "with Model.mesh.deform_mesh():\n", | ||
| " Model.mesh.data[:, -1] = TField.data[:, 0].copy()\n", | ||
| "\n", | ||
| "Model.population_control.repopulate()" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "2b37ca13-01a2-49cf-9f64-32a9d668a04e", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Fig = vis.Figure(resolution=figsize,rulers=True,margin = 80,axis=True)\n", | ||
| "# Fig.Mesh(Model.mesh)\n", | ||
| "# lv = Fig.window()\n", | ||
| "# lv.rotate('x',-45)\n", | ||
| "# lv.redisplay()" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "8b366503-56c0-4b8b-af42-a093d69d5b7a", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "ma.density = 3300. * u.kilogram / u.metre**3\n", | ||
| "ma.viscosity = 1e21 * u.pascal * u.second\n", | ||
| "Model.set_velocityBCs(left=[0.0, None, None],right=[0.0,None, None],\n", | ||
| " front=[None, 0.0, None], back=[None, 0.0, None],\n", | ||
| " bottom=[0.0, 0.0, 0.0],)\n", | ||
| "Model.freeSurface = True" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "1c25abe0-36c7-425a-8db6-ccc6eb1e7bbb", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "# Fig = vis.Figure(resolution=figsize,rulers=True,margin = 80,axis=True)\n", | ||
| "# Fig.Points(Model.swarm, Model.materialField,fn_size=2.0,discrete=True,colourBar=False)\n", | ||
| "# lv = Fig.window()\n", | ||
| "# lv.rotate('x',-45)\n", | ||
| "# lv.redisplay()" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "aa60d111-0ca5-4422-a774-5df68b9c5ab1", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "Model.solver.set_inner_method(\"mg\")" | ||
| ] | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "execution_count": null, | ||
| "id": "dfcbafd2-ea27-494f-aaaf-9b7cd105b94f", | ||
| "metadata": {}, | ||
| "outputs": [], | ||
| "source": [ | ||
| "Model.run_for(max_time, checkpoint_interval=checkpoint_interval,dt= dt_set)" | ||
| ] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "display_name": "env_uw2", | ||
| "language": "python", | ||
| "name": "env_uw2" | ||
| }, | ||
| "language_info": { | ||
| "codemirror_mode": { | ||
| "name": "ipython", | ||
| "version": 3 | ||
| }, | ||
| "file_extension": ".py", | ||
| "mimetype": "text/x-python", | ||
| "name": "python", | ||
| "nbconvert_exporter": "python", | ||
| "pygments_lexer": "ipython3", | ||
| "version": "3.11.6" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 5 | ||
| } |
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.