Example of FK rheology in UWGeodynamics 2.10.2

FK_convection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thermal Convection\n",
    "\n",
    "This example solves 2D dimensionless isoviscous thermal convection with a Rayleigh number of $10^4$, see Blankenbach *et al.* 1989 for details.\n",
    "\n",
    "**References**\n",
    "\n",
    "B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch, G. Marquart, D. Moore, P. Olson, H. Schmeling and T. Schnaubelt. A benchmark comparison for mantle convection codes. Geophysical Journal International, 98, 1, 23–38, 1989\n",
    "http://onlinelibrary.wiley.com/doi/10.1111/j.1365-246X.1989.tb05511.x/abstract"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import UWGeodynamics as GEO\n",
    "from UWGeodynamics import visualisation as vis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import underworld as uw\n",
    "from underworld import function as fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = GEO.u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boxHeight = 1000. * u.kilometer\n",
    "boxLength = 2000. * u.kilometer\n",
    "\n",
    "tempMin = 273.15 * u.degK\n",
    "tempMax = 1273.15 * u.degK\n",
    "\n",
    "refViscosity = 1e23 * u.pascal * u.second\n",
    "\n",
    "KL = boxHeight\n",
    "KT = (tempMax - tempMin)\n",
    "Kt = 1.0 * u.megayear\n",
    "KM = refViscosity * KL * Kt\n",
    "\n",
    "GEO.scaling_coefficients[\"[length]\"] = KL\n",
    "GEO.scaling_coefficients[\"[temperature]\"]= KT\n",
    "GEO.scaling_coefficients[\"[time]\"] = Kt\n",
    "GEO.scaling_coefficients[\"[mass]\"] = KM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model = GEO.Model(elementRes=(64, 64), \n",
    "                  minCoord=(0. * u.kilometer, 0. * u.kilometer), \n",
    "                  maxCoord=(2000. * u.kilometer, 1000. * u .kilometer),\n",
    "                  gravity=(0., -10. * u.meter / u.second**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model.outputDir = \"1_02_ConvectionExample\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Material property"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model.density = GEO.LinearDensity(4000. * u.kilogram / u.metre**3, thermalExpansivity=2.5e-5 / u.degK)\n",
    "Model.diffusivity = 1e-6 * u.metre**2 / u.second\n",
    "# Model.viscosity = 1e23 * u.pascal * u.second"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Raylegh number is defined as \n",
    "\n",
    "$$ Ra_0 = \\frac{\\alpha g \\Delta T h^3}{\\kappa \\mu_0} $$\n",
    "\n",
    "$\\alpha$ is thermal expansion coefficient,\n",
    "$\\kappa$ is diffusivity,\n",
    "$g$ is gravity,\n",
    "$dT$ is the difference in temperature between top and bottom,\n",
    "$h$ is the Model height,\n",
    "$mu0$ is the kinematic viscosity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = Model.thermalExpansivity\n",
    "kappa = Model.diffusivity\n",
    "g = abs(Model.gravity[-1])\n",
    "dT = KT\n",
    "h = Model.height\n",
    "mu0 = (1e23 * u.pascal * u.second) / (4000 * u.kilogram / u.metre**3)\n",
    "\n",
    "Ra0 = (alpha * g * dT * h**3) / (kappa * mu0)\n",
    "\n",
    "print(\"Rayleigh Number:\", Ra0.to_base_units())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define velocity boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VelocityBCs = Model.set_velocityBCs(left=[0., None], right=[0.,None], top=[None,0.], bottom=[None,0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define thermal boundary conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BoundaryBCs = Model.set_temperatureBCs(top=tempMin, bottom=tempMax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import underworld as uw\n",
    "from underworld import function as fn\n",
    "# set up h variable and prefactor, A\n",
    "h = GEO.nd(3.14e-3 * u.K**-1) # must use nd as this goes inside a uw.fn()\n",
    "A = 1e23 * u.Pa * u.sec\n",
    "\n",
    "# import underworld as uw\n",
    "from underworld import function as fn\n",
    "# setup FK viscosity\n",
    "Model.viscosity = A * fn.math.exp( -h * Model.temperature )\n",
    "\n",
    "# setup isoviscous\n",
    "#Model.viscosity = A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Initial Temperature perturbation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "boxLength = GEO.nd(boxLength)\n",
    "boxHeight = GEO.nd(boxHeight)\n",
    "tempMin = GEO.nd(tempMin)\n",
    "tempMax = GEO.nd(tempMax)\n",
    "\n",
    "Model.temperature.data[:] = 0.\n",
    "pertStrength = GEO.nd(200. * u.kilometer)\n",
    "deltaTemp = tempMax - tempMin\n",
    "\n",
    "for index, coord in enumerate(Model.mesh.data):\n",
    "    pertCoeff = math.cos( math.pi * coord[0] ) * math.sin( math.pi * coord[1] )\n",
    "    Model.temperature.data[index] = tempMin + deltaTemp*(boxHeight - coord[1]) + pertStrength * pertCoeff\n",
    "    Model.temperature.data[index] = max(tempMin, min(tempMax, Model.temperature.data[index]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fig = vis.Figure(figsize=(1200,400))\n",
    "Fig.Surface(Model.mesh, Model.temperature, colours=\"coolwarm\")\n",
    "Fig.save(\"Figure_1.png\")\n",
    "Fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model.solver.set_inner_method(\"mumps\")\n",
    "Model.solver.set_penalty(1e6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model.init_model(temperature=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Model.run_for(20000.0 * u.years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fig = vis.Figure(figsize=(1200,400))\n",
    "Fig.Surface(Model.mesh, Model.velocityField[0])\n",
    "Fig.VectorArrows(Model.mesh, Model.velocityField)\n",
    "Fig.save(\"Figure_1.png\")\n",
    "Fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fig = vis.Figure(figsize=(1200,400))\n",
    "Fig.Surface(Model.mesh, GEO.dimensionalise(Model.temperature, u.degK), colours=\"coolwarm\")\n",
    "Fig.VectorArrows(Model.mesh, Model.velocityField)\n",
    "Fig.save(\"Figure_2.png\")\n",
    "Fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fig = vis.Figure(figsize=(1200,400))\n",
    "Fig.Surface(Model.mesh, GEO.dimensionalise(Model.viscosityField, u.Pa * u.sec), logScale=True)\n",
    "Fig.VectorArrows(Model.mesh, Model.velocityField)\n",
    "Fig.save(\"Figure_3.png\")\n",
    "Fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fig = vis.Figure(figsize=(1200,400))\n",
    "# Fig.Surface(Model.mesh, Model.velocityField[0])\n",
    "# Fig.VectorArrows(Model.mesh, Model.velocityField)\n",
    "# Fig.save(\"Figure_2.png\")\n",
    "# Fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}