Merge branch 'development' of github.com:underworldcode/underworld2 i…

…nto development

CHANGES.md

@@ -1,5 +1,10 @@
 CHANGES: Underworld2
 =======================
+Release 2.16.0 [2024]
+---------------------------
+New:
+* New 3D free surface implementation (tested). The vertical coordinates of the mesh nodes are rebuilt by being advected and interpolated from the surface velocities, and then solving the steady state heat equation to get a uniform distribution. See the related example in docs/UWGeodynamics/examples
+
 
 Release 2.15.0 [2023-04-19]
 ---------------------------
@@ -12,7 +17,7 @@ Changes:
 
 Fixes:
  * UWGeodynamics - add dynamic heating back into the advection diffusion solver,
-  https://github.com/underworldcode/underworld2/issues/669
+    https://github.com/underworldcode/underworld2/issues/669
  * Using updated Badlands-2.2.3 without license issue. 
 
 
@@ -162,7 +167,7 @@ Docker:
 * Images minimised with unnecessary items removed. 
 * XVFB no longer required for image generation within
   container.
- 
+
 API changes:
 * `glucifer` module moved inside `underworld` and 
    renamed `visualisation`.

docs/UWGeodynamics/examples/1_23_05_FreeSurface_3D_Relaxation.ipynb

@@ -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
+}

underworld/UWGeodynamics/_freesurface.py

@@ -1,5 +1,6 @@
 from __future__ import print_function, absolute_import
 from scipy.interpolate import interp1d
+from scipy.interpolate import CloughTocher2DInterpolator
 import underworld as uw
 from underworld.scaling import non_dimensionalise as nd
 
@@ -29,7 +30,7 @@ def __init__(self, model):
 
         # Create the tools
         self.TField = self._init_mesh.add_variable(nodeDofCount=1)
-        self.TField.data[:, 0] = self._init_mesh.data[:, 1].copy()
+        self.TField.data[:, 0] = self._init_mesh.data[:, -1].copy()
 
         self.top = self.model.top_wall
         self.bottom = self.model.bottom_wall
@@ -53,22 +54,42 @@ def _solve_sle(self):
     def _advect_surface(self, dt):
 
         if self.top:
-            # Extract top surface
-            x = self.model.mesh.data[self.top.data][:, 0]
-            y = self.model.mesh.data[self.top.data][:, 1]
-
-            # Extract velocities from top
-            vx = self.model.velocityField.data[self.top.data][:, 0]
-            vy = self.model.velocityField.data[self.top.data][:, 1]
-
-            # Advect top surface
-            x2 = x + vx * nd(dt)
-            y2 = y + vy * nd(dt)
-
-            # Spline top surface
-            f = interp1d(x2, y2, kind='cubic', fill_value='extrapolate')
-
-            self.TField.data[self.top.data, 0] = f(x)
+            if self.model.mesh.dim == 2:
+                # Extract top surface
+                x = self.model.mesh.data[self.top.data][:, 0]
+                y = self.model.mesh.data[self.top.data][:, 1]
+
+                # Extract velocities from top
+                vx = self.model.velocityField.data[self.top.data][:, 0]
+                vy = self.model.velocityField.data[self.top.data][:, 1]
+
+                # Advect top surface
+                x2 = x + vx * nd(dt)
+                y2 = y + vy * nd(dt)
+
+                # Spline top surface
+                f = interp1d(x2, y2, kind='cubic', fill_value='extrapolate')
+
+                self.TField.data[self.top.data, 0] = f(x)
+            else:
+                # Extract top surface
+                x = self.model.mesh.data[self.top.data][:, 0]
+                y = self.model.mesh.data[self.top.data][:, 1]
+                z = self.model.mesh.data[self.top.data][:, -1]
+
+                # Extract velocities from top
+                vx = self.model.velocityField.data[self.top.data][:, 0]
+                vy = self.model.velocityField.data[self.top.data][:, 1]
+                vz = self.model.velocityField.data[self.top.data][:, -1]
+
+                # Advect top surface
+                x2 = x + vx * nd(dt)
+                y2 = y + vy * nd(dt)
+                z2 = z + vz * nd(dt)
+
+                # Spline top surface
+                f = CloughTocher2DInterpolator((x2, y2), z2)
+                self.TField.data[self.top.data, 0] = f((x,y))
         uw.mpi.barrier()
         self.TField.syncronise()