{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The rotating cone problem\n", "\n", "Over a unit square a blob of heat is rotated one revolution subject to a Péclet number of 1e6. \n", "The evolution of the heat is governed by the advection diffusion equation. \n", "\n", "In Underworld this equation can be numerically approximated with 2 methods: the Streamline Upwind/Petrov-Galerkin method (**SUPG**); and the Semi-Lagrangian Crank-Nicholson methd (**SLCN**). Various numerical parameters can be investigate in the model below.\n", "\n", "This classical test problem is taken from J. Donea and A. Huerta, _Finite Element Methods for Flow Problems (2003)_ with an initial field (temperature in this case) given:\n", "\n", " \n", "$$\n", "T(x,y,0) = \n", "\\begin{array}{cc}\n", " \\biggl\\{ & \n", " \\begin{array}{cc}\n", " \\frac{1}{4}(1+\\cos(\\pi X))(1+\\cos(\\pi Y)) & X^2 + Y^2 \\leq 1, \\\\\n", " 0 & X^2 + Y^2 \\gt 1, \\\\\n", " \\end{array}\n", "\\end{array}\n", "$$\n", "\n", "where $$\n", "(X, Y) = (\\space (x-x_0) / \\sigma ,\\space (y-y_0) / \\sigma)\\space)\n", "$$ \n", "\n", "with $(x_0, y_0)$ as the center of the blob and the boundary condition is $T=0$ on all walls.\n", "\n", "The domain is a unit square $[-0.5,0.5]x[-0.5,0.5]$. The advection field is of pure rotation with velocity as $V(x,y) = (-y,x)$. A static timestep $dt = \\frac{2\\pi}{200} $.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import underworld as uw\n", "import underworld.function as fn\n", "import underworld.visualisation as vis\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# numerical parameters\n", "nEls = 30 # number of elements in x & y\n", "dt = 2*np.pi/200 # timestep size\n", "viz = (1 if uw.mpi.size==1 else 0)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def runModel(method, nEls = nEls, dt = dt, elementType='Q1', viz = 1, rotations = 1., diffuseOnly=False):\n", " \"\"\"\n", " Defines and executes \"rotating cone problem\" numerical model.\n", " \n", " Args\n", " ----\n", " method : str, \"SUPG\" or \"SLCN\"\n", " The numerical scheme used to evolve the temperature in time.\n", " \n", " nEls : int\n", " The number of elements of that define the numerical domain, \n", " equal in both axis.\n", " \n", " dt : float\n", " The dimensionless timestep size used to integrate the temperature \n", " field one time step.\n", " \n", " Q1 : str\n", " The underworld finite element type used to approximate the temperature field \n", " \n", " viz : bool\n", " Show the initial and final timesteps of the model\n", " \n", " rotations: float\n", " Number of rotations to perform.\n", " \n", " diffuseOnly : bool\n", " Don't advect the model (zero velocity field). Only diffuse process allowed.\n", " \n", " Returns\n", " -------\n", " maxT : float\n", " The maximum temperature after one revolution.\n", " \n", " avgT : float\n", " The average temperature after one revolution.\n", " \n", " num_t : float\n", " Wall time (fractional seconds) taken to complete.\n", " \n", " steps : int\n", " The number of discrete time steps to taken to complete.\n", " \"\"\"\n", " \n", " # default model parameters\n", " sigma = 0.2 # width of blob\n", " x_0 = (1./6, 1./6) # position of blob\n", " kappa = 1e-6 # thermal diffusion (entire domain)\n", " # unit 2D mesh max/min coords\n", " minCoord = [-0.5,-0.5]\n", " maxCoord = [0.5,0.5]\n", " \n", " # build FE mesh, velocity and temperature fields\n", " mesh = uw.mesh.FeMesh_Cartesian(elementType = elementType, \n", " elementRes = (nEls,nEls),\n", " minCoord = minCoord,\n", " maxCoord = maxCoord )\n", "\n", " tField = mesh.add_variable(1)\n", " vField = mesh.add_variable(2)\n", "\n", " all_walls = mesh.specialSets[\"AllWalls_VertexSet\"]\n", " \n", " # initialise fields\n", " x_vel = -1.*uw.function.coord()[1]\n", " y_vel = 1.*uw.function.coord()[0]\n", " vField.data[:,0] = x_vel.evaluate(mesh).reshape(-1)\n", " vField.data[:,1] = y_vel.evaluate(mesh).reshape(-1)\n", " \n", " if diffuseOnly:\n", " vField.data[...] = 0.\n", "\n", " fn_X = (fn.input()[0]-x_0[0])/sigma\n", " fn_Y = (fn.input()[1]-x_0[1])/sigma\n", " \n", " fn_inside = ( fn.math.pow(fn_X,2.) + fn.math.pow(fn_Y,2.) ) <= 1.\n", " fn_hill = 0.25 * ( 1.+fn.math.cos(np.pi*fn_X) ) * ( 1.+fn.math.cos(np.pi*fn_Y) )\n", "\n", " fn_ic = fn.branching.conditional([(fn_inside, fn_hill ), \n", " (True , 0.)])\n", " \n", " tField.data[:] = fn_ic.evaluate(mesh)\n", "\n", " # setup visualisation\n", " fig = vis.Figure()\n", " fig.VectorArrows(mesh,vField)\n", "# fig.Mesh(mesh)\n", " fig.Surface(mesh,tField, valueRange=[0.,1.]);\n", " \n", " # define advection diffusion system\n", " dbc = uw.conditions.DirichletCondition(variable=tField,\n", " indexSetsPerDof=all_walls)\n", "\n", " # extra variable needed with supg\n", " tDotField = None\n", " if method == \"SUPG\":\n", " tDotField = mesh.add_variable( nodeDofCount=1, dataType='double')\n", " \n", " system = uw.systems.AdvectionDiffusion( phiField = tField,\n", " method = method,\n", " phiDotField = tDotField,\n", " velocityField = vField,\n", " fn_diffusivity = kappa,\n", " conditions = [dbc], kwargs=None )\n", " \n", " from timeit import default_timer as timer\n", "\n", " t0 = timer()\n", " num_t = 0\n", " time = 0\n", " steps = 0\n", " end_t = 2.*np.pi*rotations\n", "\n", " if dt < 0:\n", " dt = system.get_max_dt()\n", " \n", " # time loop\n", " while (time+dt)" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "SLCN: t_max=1.000000000000000, t_avg=0.039831195790397, wall time=2.7290e-06, its=0, dt=3.1416e-02\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "SLCN (diffusion only): t_max=0.998415344757853, t_avg=0.039831185188472, wall time=1.5593e+01, its=200, dt=3.1416e-02\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "SLCN: t_max=0.852184812161098, t_avg=0.040071167501389, wall time=1.6513e+01, its=200, dt=3.1416e-02\n" ] }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "SUPG: t_max=0.494862381241538, t_avg=0.039927886622725, wall time=6.6796e+00, its=200, dt=3.1416e-02\n" ] } ], "source": [ "results_slcn = runModel(\"SLCN\", nEls = nEls, dt = dt, viz = viz, rotations=0.) # initial conditions only\n", "results_slcn = runModel(\"SLCN\", nEls = nEls, dt = dt, viz = viz, diffuseOnly=True) # diffuse only\n", "results_slcn = runModel(\"SLCN\", nEls = nEls, dt = dt, viz = viz)\n", "results_supg = runModel(\"SUPG\", nEls = nEls, dt = dt, viz = viz)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# test the max temperature and avg. temperature iff nEls == 30\n", "# against expected results for 1 proc\n", "if uw.mpi.size == 0:\n", " if [results_supg[3],results_supg[4]] == [200,30]:\n", " expected = [0.494862, 0.039927]\n", " if not np.allclose([results_supg[0],results_supg[1]], \n", " expected, rtol=1e-2):\n", " raise RuntimeError(\"Failed SUPG test\")\n", " \n", " if [results_slcn[3],results_slcn[4]] == [200,30]:\n", " expected = [0.852184,0.040071]\n", " if not np.allclose([results_slcn[0], results_slcn[1]],\n", " expected, rtol=1e-2):\n", " raise RuntimeError(\"Failed SLCN test\") " ] } ], "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.6.7" } }, "nbformat": 4, "nbformat_minor": 2 }