Chapter20_Blankenbach_Case1a.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8ea95bde",
   "metadata": {},
   "source": [
    "# Chapter 20 - Blankenbach - Case1a\n",
    "\n",
    "This notebook implements the Blankenbacch Case1a isoviscous thermal convection benchmark.\n",
    "\n",
    "This notebook puts all of the major function definitions into separate source files and contains only the high-level program code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adabca4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary packages\n",
    "using SparseArrays\n",
    "using LinearAlgebra\n",
    "using IterativeSolvers\n",
    "using WriteVTK\n",
    "using Printf\n",
    "using Statistics \n",
    "include(\"Grid.jl\")\n",
    "include(\"Markers.jl\")\n",
    "include(\"Stokes.jl\")\n",
    "include(\"Temperature.jl\")\n",
    "include(\"GridOperations.jl\")\n",
    "\n",
    "# note that we import pyplot last to avoid a name conflict with grid.\n",
    "using PyPlot\n",
    "include(\"Visualization.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60ea66dc",
   "metadata": {},
   "source": [
    "## Model setup\n",
    "\n",
    "Boundary conditions are insulating on the sides, fixed temperature top and bottom.\n",
    "\n",
    "Initial condition is:\n",
    "$$\n",
    "T(x,y) = T_{surf} + \\frac{(T_{btm}-T_{surf})}{H} y + dT\\cos(\\frac{\\pi x}{W}) \\sin(\\frac{\\pi y}{H})\n",
    "$$\n",
    "\n",
    "Note that the origin of the coordinate system is in the upper right. $y$ increases downward, $x$ increases rightward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afccc752",
   "metadata": {},
   "outputs": [],
   "source": [
    "function blankenbach_initial_conditions!(markers::Markers,Ttop,Tbtm,W,H)\n",
    "    eta = markers.scalarFields[\"eta\"]\n",
    "    T = markers.scalarFields[\"T\"]\n",
    "    for i in 1:markers.nmark\n",
    "        mx = markers.x[1,i]\n",
    "        my = markers.x[2,i]\n",
    "        # T(x,y) = Tsurf + (Tbtm-Tsurf)/H * y + dT*cos(pi*x/W)\n",
    "        markers.scalars[eta,i] = 1e23\n",
    "        markers.scalars[T,i] = Ttop + (Tbtm-Ttop)/H*my + 10.0*cos(pi*mx/W)*sin(pi*my/H)\n",
    "    end\n",
    "end\n",
    "\n",
    "function compute_nusselt(grid::CartesianGrid,T::Matrix{Float64},kThermal,Tbcval)\n",
    "        dTdy = (T[2,:]-T[1,:])/(grid.yc[2]-grid.yc[1])\n",
    "        kdTdy = kThermal*dTdy\n",
    "        qcond = kThermal*(Tbcval[4]-Tbcval[3])/grid.H\n",
    "        Nu = mean(kdTdy)/qcond;\n",
    "    return Nu\n",
    "end\n",
    "\n",
    "function run_blankenbach_case1a(end_step=1;end_time=nothing)\n",
    "    nx = 97\n",
    "    ny = 97\n",
    "    W = 1e6\n",
    "    H = 1e6\n",
    "    \n",
    "    gx = 0.0\n",
    "    gy = 10.0\n",
    "    alpha = 2.5e-5; # thermal expansivity\n",
    "\n",
    "    bc = BoundaryConditions(0,0,0,0) # currently does nothing but is required argument to stokes solver.\n",
    "    Tbcval = [0.,0.,273.,1273.]   \n",
    "    Tbctype = [-1,-1,1,1] # insulating left and right, isothermal top and bottom\n",
    "    markx = 6\n",
    "    marky = 6\n",
    "    seconds_in_year = 3.15e7\n",
    "    plot_interval = 100e6*seconds_in_year # plot interval in seconds\n",
    "    dtmax = plot_interval\n",
    "    \n",
    "    grid = CartesianGrid(W,H,nx,ny)\n",
    "    println(\"Creating Markers...\")\n",
    "    @time markers = Markers(grid,[\"T\",\"rho\",\"eta\"],[\"material\"] ; nmx=markx,nmy=marky,random=true)\n",
    "    println(\"Initial condition...\")\n",
    "    @time blankenbach_initial_conditions!(markers,Tbcval[3],Tbcval[4],W,H)\n",
    "\n",
    "    # define arrays for k, rho, cp, H at the basic nodes. Fill them with constant values for now.\n",
    "    kThermal = 5.0 .*ones(grid.ny,grid.nx);\n",
    "    Cp = 1250. .* ones(grid.ny,grid.nx);\n",
    "    H = zeros(grid.ny,grid.nx);\n",
    "    rho0 = 4000.0*ones(grid.ny,grid.nx);\n",
    "    \n",
    "    kappa = kThermal[1,1]/rho0[1,1]/Cp[1,1]\n",
    "    diffusion_timescale = (grid.x[2]-grid.x[1])^2/kappa\n",
    "    dtmax = minimum([dtmax,diffusion_timescale])\n",
    "    \n",
    "    local time = 0.0\n",
    "    iout=0\n",
    "    last_plot = 0.0\n",
    "\n",
    "    local rho\n",
    "    local rho_vx\n",
    "    local rho_vy\n",
    "\n",
    "    local eta_s\n",
    "    local eta_n\n",
    "    local T\n",
    "    local dt\n",
    "    local dTmax\n",
    "    local dTemp\n",
    "    local Tnew\n",
    "    local Tlast\n",
    "    \n",
    "    nu_values=zeros(0)\n",
    "    time_values=zeros(0)\n",
    "    \n",
    "    output_dir = \"case1a\"\n",
    "    local terminate = false\n",
    "    local itime=0\n",
    "    while !terminate\n",
    "        if itime >= end_step || (end_time != nothing && time >= end_time)\n",
    "            terminate = true\n",
    "        end\n",
    "        \n",
    "        rho_var = markers.scalarFields[\"rho\"]\n",
    "        T_var = markers.scalarFields[\"T\"]\n",
    "        for i in 1:markers.nmark\n",
    "            markers.scalars[rho_var,i] = 4000.0*(1.0-alpha*(markers.scalars[T_var,i]-Tbcval[3]))\n",
    "        end\n",
    "        \n",
    "        # Get density and viscosity from the markers.\n",
    "        rho_new,eta_s_new = marker_to_stag(markers,grid,[\"rho\",\"eta\"],\"basic\");        \n",
    "        eta_n_new,Tlast_new = marker_to_stag(markers,grid,[\"eta\",\"T\"],\"center\")\n",
    "        rho_vx_new, = marker_to_stag(markers,grid,[\"rho\",],\"vx\")\n",
    "        rho_vy_new, = marker_to_stag(markers,grid,[\"rho\",],\"vy\")\n",
    "        \n",
    "        if time > 0.0\n",
    "            if any(isnan.(eta_s_new))\n",
    "                println(\"found nan values\")\n",
    "            end\n",
    "            replace_nan!(rho,rho_new)\n",
    "            replace_nan!(Tlast,Tlast_new)\n",
    "            replace_nan!(eta_s,eta_s_new)\n",
    "            replace_nan!(eta_n,eta_n_new)\n",
    "            replace_nan!(rho_vx,rho_vx_new)\n",
    "            replace_nan!(rho_vy,rho_vy_new)\n",
    "        end\n",
    "\n",
    "        Tlast = copy(Tlast_new)\n",
    "        rho   = copy(rho_new)\n",
    "        eta_s = copy(eta_s_new) \n",
    "        eta_n = copy(eta_n_new)\n",
    "        rho_vx = copy(rho_vx_new)\n",
    "        rho_vy = copy(rho_vy_new)\n",
    "\n",
    "        # assemble and solve the stokes equations\n",
    "        L,R = form_stokes(grid,eta_s,eta_n,rho_vx,rho_vy,bc,gx,gy;dt=0.0) \n",
    "        solution = L\\R\n",
    "        vx,vy,P = unpack(solution,grid;ghost=true)\n",
    "        vxc,vyc = velocity_to_centers(grid,vx,vy)\n",
    "        dt = compute_timestep(grid,vxc,vyc;cfl=0.25,dtmax=dtmax)\n",
    "        \n",
    "        Tlast = ghost_temperature_center(grid,Tlast,Tbctype,Tbcval);\n",
    "\n",
    "        dTmax = Inf\n",
    "        for titer=1:2\n",
    "            # assemble and solve the energy equation\n",
    "            println(\"Trying with timestep \",dt)\n",
    "            L,R = assemble_energy_equation_center(grid,rho0,Cp,kThermal,H,Tlast,dt,Tbctype,Tbcval);\n",
    "            Tnew = L\\R;\n",
    "            Tnew = reshape(Tnew,grid.ny,grid.nx);\n",
    "            Tnew = ghost_temperature_center(grid,Tnew,Tbctype,Tbcval);\n",
    "\n",
    "            dTemp = Tnew-Tlast\n",
    "            # compute the maximum temperature change\n",
    "            dTmax = maximum(abs.(dTemp[2:end-1,2:end-1]))\n",
    "            println(\"dTmax=\",dTmax,\" dt=\",dt/3.15e7/1e6)\n",
    "            dt = min(dt,dTmax < 10.0 ? dt : dt*10.0/dTmax)\n",
    "            if dTmax < 10.0\n",
    "                break\n",
    "            end\n",
    "        end\n",
    "        # compute Nusselt number\n",
    "        Nu = compute_nusselt(grid,Tnew,kThermal[1,1],Tbcval)\n",
    "        \n",
    "        # store Nusselt number for plotting\n",
    "        append!(nu_values,Nu)\n",
    "        append!(time_values,time)\n",
    "        \n",
    "        dT_subgrid_node = subgrid_temperature_relaxation_center!(markers,grid,Tlast,Cp[1,1],kThermal[1,1],dt)\n",
    "        dT_remaining = dTemp - dT_subgrid_node\n",
    "\n",
    "        cell_center_change_to_markers!(markers,grid,dT_remaining,\"T\")\n",
    "\n",
    "        if time == 0.0 || time - last_plot >= plot_interval || terminate\n",
    "            last_plot = time \n",
    "            name = @sprintf(\"%s/viz.%04d.vtr\",output_dir,iout)\n",
    "            vn = velocity_to_basic_nodes(grid,vxc,vyc)\n",
    "\n",
    "        visualization(grid,rho,eta_s,vn,P,Tnew[2:end-1,2:end-1],time/seconds_in_year;filename=name)\n",
    "            name = @sprintf(\"%s/markers.%04d.vtp\",output_dir,iout)\n",
    "            visualization(markers,time/seconds_in_year;filename=name)\n",
    "            iout += 1\n",
    "        end\n",
    "\n",
    "        move_markers_rk4!(markers,grid,vx,vy,dt,continuity_weight=1/3)\n",
    "        time += dt\n",
    "        itime += 1\n",
    "        println(\"Finished Step \",itime,\" time=\",time/seconds_in_year/1e6,\" Myr Nu=\",Nu,\" \",extrema(Tnew[2:end-1,2:end-1]))\n",
    "    end\n",
    "    return nu_values, time_values, markers, grid, dTemp, Tnew, Tlast\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48d979ff",
   "metadata": {},
   "source": [
    "Run the benchmark code for 10 timesteps and produce a ProfileSVG view:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "739b1032",
   "metadata": {},
   "outputs": [],
   "source": [
    "using ProfileSVG\n",
    "@time nu,t,markers,grid,dTemp,Tnew,Tlast = run_blankenbach_case1a(10,end_time=6e9*3.15e7)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c070c6cd",
   "metadata": {},
   "source": [
    "Run the benchmark code long enough to approach steady state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c65c0fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "@time nu,t,markers,grid,dTemp,Tnew,Tlast = run_blankenbach_case1a(Inf,end_time=6e9*3.15e7);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccdb6e32",
   "metadata": {},
   "source": [
    "Plot the time-evolution of the Nusselt number. Add a red dashed line that indicates the preferred benchmark value from Gerya."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9188e2d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "figure()\n",
    "plot(t/3.15e7/1e6,nu,\".\")\n",
    "plot([0.,6000.],[4.8844,4.8844],\"r--\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35215835",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "pcolor(Tlast)\n",
    "colorbar()\n",
    "show()\n",
    "# print(size(Tnew))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23d8fcd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure()\n",
    "scatter(markers.x[1,1:markers.nmark],markers.x[2,1:markers.nmark],c=markers.scalars[markers.scalarFields[\"T\",],1:markers.nmark],s=0.1)\n",
    "colorbar()\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90118aea",
   "metadata": {},
   "outputs": [],
   "source": [
    "minimum(markers.scalars[markers.scalarFields[\"T\"],1:markers.nmark])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.5",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}