{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 13: Adding new 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": [
    "This demo illustrates the steps involved in adding new parabolic terms for the scalar\n",
    "advection equation. In particular, we will add an anisotropic diffusion. We begin by\n",
    "defining the hyperbolic (advection) part of the advection-diffusion equation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq\n",
    "using Trixi\n",
    "\n",
    "\n",
    "advection_velocity = (1.0, 1.0)\n",
    "equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Define a new parabolic equation type\n",
    "\n",
    "Next, we define a 2D parabolic diffusion term type. This is similar to `LaplaceDiffusion2D`\n",
    "except that the `diffusivity` field refers to a spatially constant diffusivity matrix now. Note that\n",
    "`ConstantAnisotropicDiffusion2D` has a field for `equations_hyperbolic`. It is useful to have\n",
    "information about the hyperbolic system available to the parabolic part so that we can reuse\n",
    "functions defined for hyperbolic equations (such as `varnames`).\n",
    "\n",
    "The abstract type `Trixi.AbstractEquationsParabolic` has three parameters: `NDIMS` (the spatial dimension,\n",
    "e.g., 1D, 2D, or 3D), `NVARS` (the number of variables), and `GradientVariable`, which we set as\n",
    "`GradientVariablesConservative`. This indicates that the gradient should be taken with respect to the\n",
    "conservative variables (e.g., the same variables used in `equations_hyperbolic`). Users can also take\n",
    "the gradient with respect to a different set of variables; see, for example, the implementation of\n",
    "`CompressibleNavierStokesDiffusion2D`, which can utilize either \"primitive\" or \"entropy\" variables."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "struct ConstantAnisotropicDiffusion2D{E, T} <: Trixi.AbstractEquationsParabolic{2, 1, GradientVariablesConservative}\n",
    "  diffusivity::T\n",
    "  equations_hyperbolic::E\n",
    "end\n",
    "\n",
    "varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiffusion2D) =\n",
    "  varnames(variable_mapping, equations_parabolic.equations_hyperbolic)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we define the viscous flux function. We assume that the mixed hyperbolic-parabolic system\n",
    "is of the form\n",
    "$$\n",
    "\\partial_t u(t,x) + \\partial_x (f_1(u) - g_1(u, \\nabla u))\n",
    "                  + \\partial_y (f_2(u) - g_2(u, \\nabla u)) = 0\n",
    "$$\n",
    "where $f_1(u)$, $f_2(u)$ are the hyperbolic fluxes and $g_1(u, \\nabla u)$, $g_2(u, \\nabla u)$ denote\n",
    "the viscous fluxes. For anisotropic diffusion, the viscous fluxes are the first and second components\n",
    "of the matrix-vector product involving `diffusivity` and the gradient vector.\n",
    "\n",
    "Here, we specialize the flux to our new parabolic equation type `ConstantAnisotropicDiffusion2D`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function Trixi.flux(u, gradients, orientation::Integer, equations_parabolic::ConstantAnisotropicDiffusion2D)\n",
    "  @unpack diffusivity = equations_parabolic\n",
    "  dudx, dudy = gradients\n",
    "  if orientation == 1\n",
    "    return SVector(diffusivity[1, 1] * dudx + diffusivity[1, 2] * dudy)\n",
    "  else # if orientation == 2\n",
    "    return SVector(diffusivity[2, 1] * dudx + diffusivity[2, 2] * dudy)\n",
    "  end\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Defining boundary conditions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Trixi.jl's implementation of parabolic terms discretizes both the gradient and divergence\n",
    "using weak formulation. In other words, we discretize the system\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\bm{q} &= \\nabla u \\\\\n",
    "\\bm{\\sigma} &= \\begin{pmatrix} g_1(u, \\bm{q}) \\\\ g_2(u, \\bm{q}) \\end{pmatrix} \\\\\n",
    "\\text{viscous contribution } &= \\nabla \\cdot \\bm{\\sigma}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Boundary data must be specified for all spatial derivatives, e.g., for both the gradient\n",
    "equation $\\bm{q} = \\nabla u$ and the divergence of the viscous flux\n",
    "$\\nabla \\cdot \\bm{\\sigma}$. We account for this by introducing internal `Gradient`\n",
    "and `Divergence` types which are used to dispatch on each type of boundary condition.\n",
    "\n",
    "As an example, let us introduce a Dirichlet boundary condition with constant boundary data."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "struct BoundaryConditionConstantDirichlet{T <: Real}\n",
    "  boundary_value::T\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This boundary condition contains only the field `boundary_value`, which we assume to be some\n",
    "real-valued constant which we will impose as the Dirichlet data on the boundary.\n",
    "\n",
    "Boundary conditions have generally been defined as \"callable structs\" (also known as \"functors\").\n",
    "For each boundary condition, we need to specify the appropriate boundary data to return for both\n",
    "the `Gradient` and `Divergence`. Since the gradient is operating on the solution `u`, the boundary\n",
    "data should be the value of `u`, and we can directly impose Dirichlet data."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector,\n",
    "                                                                          x, t, operator_type::Trixi.Gradient,\n",
    "                                                                          equations_parabolic::ConstantAnisotropicDiffusion2D)\n",
    "  return boundary_condition.boundary_value\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "While the gradient acts on the solution `u`, the divergence acts on the viscous flux $\\bm{\\sigma}$.\n",
    "Thus, we have to supply boundary data for the `Divergence` operator that corresponds to $\\bm{\\sigma}$.\n",
    "However, we've already imposed boundary data on `u` for a Dirichlet boundary condition, and imposing\n",
    "boundary data for $\\bm{\\sigma}$ might overconstrain our problem.\n",
    "\n",
    "Thus, for the `Divergence` boundary data under a Dirichlet boundary condition, we simply return\n",
    "`flux_inner`, which is boundary data for $\\bm{\\sigma}$ computed using the \"inner\" or interior solution.\n",
    "This way, we supply boundary data for the divergence operation without imposing any additional conditions."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector,\n",
    "                                                                          x, t, operator_type::Trixi.Divergence,\n",
    "                                                                          equations_parabolic::ConstantAnisotropicDiffusion2D)\n",
    "  return flux_inner\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "### A note on the choice of gradient variables\n",
    "\n",
    "It is often simpler to transform the solution variables (and solution gradients) to another set of\n",
    "variables prior to computing the viscous fluxes (see `CompressibleNavierStokesDiffusion2D`\n",
    "for an example of this). If this is done, then the boundary condition for the `Gradient` operator\n",
    "should be modified accordingly as well.\n",
    "\n",
    "## Putting things together\n",
    "\n",
    "Finally, we can instantiate our new parabolic equation type, define boundary conditions,\n",
    "and run a simulation. The specific anisotropic diffusion matrix we use produces more\n",
    "dissipation in the direction $(1, -1)$ as an isotropic diffusion.\n",
    "\n",
    "For boundary conditions, we impose that $u=1$ on the left wall, $u=2$ on the bottom\n",
    "wall, and $u = 0$ on the outflow walls. The initial condition is taken to be $u = 0$.\n",
    "Note that we use `BoundaryConditionConstantDirichlet` only for the parabolic boundary\n",
    "conditions, since we have not defined its behavior for the hyperbolic part."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi: SMatrix\n",
    "diffusivity = 5.0e-2 * SMatrix{2, 2}([2 -1; -1 2])\n",
    "equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic);\n",
    "\n",
    "boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)),\n",
    "                                    y_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(2.0)),\n",
    "                                    y_pos = boundary_condition_do_nothing,\n",
    "                                    x_pos = boundary_condition_do_nothing)\n",
    "\n",
    "boundary_conditions_parabolic = (; x_neg = BoundaryConditionConstantDirichlet(1.0),\n",
    "                                   y_neg = BoundaryConditionConstantDirichlet(2.0),\n",
    "                                   y_pos = BoundaryConditionConstantDirichlet(0.0),\n",
    "                                   x_pos = BoundaryConditionConstantDirichlet(0.0));\n",
    "\n",
    "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)\n",
    "\n",
    "semi = SemidiscretizationHyperbolicParabolic(mesh,\n",
    "                                             (equations_hyperbolic, equations_parabolic),\n",
    "                                             initial_condition, solver;\n",
    "                                             boundary_conditions=(boundary_conditions_hyperbolic,\n",
    "                                                                  boundary_conditions_parabolic))\n",
    "\n",
    "tspan = (0.0, 2.0)\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);\n",
    "\n",
    "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
}