{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 12: Parabolic terms"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Note:** To improve responsiveness via caching, the notebooks are updated only once a week. They are only\n",
    "available for the latest stable release of Trixi.jl at the time of caching."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Experimental support for parabolic diffusion terms is available in Trixi.jl.\n",
    "This demo illustrates parabolic terms for the advection-diffusion equation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq\n",
    "using Trixi"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Splitting a system into hyperbolic and parabolic parts."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For a mixed hyperbolic-parabolic system, we represent the hyperbolic and parabolic\n",
    "parts of the system  separately. We first define the hyperbolic (advection) part of\n",
    "the advection-diffusion equation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "advection_velocity = (1.5, 1.0)\n",
    "equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we define the parabolic diffusion term. The constructor requires knowledge of\n",
    "`equations_hyperbolic` to be passed in because the `LaplaceDiffusion2D` applies\n",
    "diffusion to every variable of the hyperbolic system."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "diffusivity = 5.0e-2\n",
    "equations_parabolic = LaplaceDiffusion2D(diffusivity, equations_hyperbolic);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Boundary conditions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As with the equations, we define boundary conditions separately for the hyperbolic and\n",
    "parabolic part of the system. For this example, we impose inflow BCs for the hyperbolic\n",
    "system (no condition is imposed on the outflow), and we impose Dirichlet boundary conditions\n",
    "for the parabolic equations. Both `BoundaryConditionDirichlet` and `BoundaryConditionNeumann`\n",
    "are defined for `LaplaceDiffusion2D`.\n",
    "\n",
    "The hyperbolic and parabolic boundary conditions are assumed to be consistent with each other."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))\n",
    "\n",
    "boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])),\n",
    "                                    y_neg = boundary_condition_zero_dirichlet,\n",
    "                                    y_pos = boundary_condition_do_nothing,\n",
    "                                    x_pos = boundary_condition_do_nothing)\n",
    "\n",
    "boundary_conditions_parabolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])),\n",
    "                                   y_neg = boundary_condition_zero_dirichlet,\n",
    "                                   y_pos = boundary_condition_zero_dirichlet,\n",
    "                                   x_pos = boundary_condition_zero_dirichlet);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Defining the solver and mesh"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The process of creating the DG solver and mesh is the same as for a purely\n",
    "hyperbolic system of equations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)\n",
    "coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))\n",
    "coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))\n",
    "mesh = TreeMesh(coordinates_min, coordinates_max,\n",
    "                initial_refinement_level=4,\n",
    "                periodicity=false, n_cells_max=30_000) # set maximum capacity of tree data structure\n",
    "\n",
    "initial_condition = (x, t, equations) -> SVector(0.0);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Semidiscretizing and solving"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "To semidiscretize a hyperbolic-parabolic system, we create a `SemidiscretizationHyperbolicParabolic`.\n",
    "This differs from a `SemidiscretizationHyperbolic` in that we pass in a `Tuple` containing both the\n",
    "hyperbolic and parabolic equation, as well as a `Tuple` containing the hyperbolic and parabolic\n",
    "boundary conditions."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "semi = SemidiscretizationHyperbolicParabolic(mesh,\n",
    "                                             (equations_hyperbolic, equations_parabolic),\n",
    "                                             initial_condition, solver;\n",
    "                                             boundary_conditions=(boundary_conditions_hyperbolic,\n",
    "                                                                  boundary_conditions_parabolic))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The rest of the code is identical to the hyperbolic case. We create a system of ODEs through\n",
    "`semidiscretize`, defining callbacks, and then passing the system to OrdinaryDiffEq.jl."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tspan = (0.0, 1.5)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "callbacks = CallbackSet(SummaryCallback())\n",
    "time_int_tol = 1.0e-6\n",
    "sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol,\n",
    "            ode_default_options()..., callback=callbacks);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can now visualize the solution, which develops a boundary layer at the outflow boundaries."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Package versions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "These results were obtained using the following versions."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using InteractiveUtils\n",
    "versioninfo()\n",
    "\n",
    "using Pkg\n",
    "Pkg.status([\"Trixi\", \"OrdinaryDiffEq\", \"Plots\"],\n",
    "           mode=PKGMODE_MANIFEST)"
   ],
   "metadata": {},
   "execution_count": null
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}