{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 15: Adding parabolic source 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": [
    "Source terms which depend on the gradient of the solution are available by specifying\n",
    "`source_terms_parabolic`. This demo illustrates parabolic source terms for the\n",
    "advection-diffusion equation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEqLowStorageRK\n",
    "using Trixi\n",
    "\n",
    "const a = 0.1\n",
    "const nu = 0.1\n",
    "const beta = 0.3\n",
    "\n",
    "equations = LinearScalarAdvectionEquation1D(a)\n",
    "equations_parabolic = LaplaceDiffusion1D(nu, equations)\n",
    "\n",
    "solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Gradient-dependent source terms"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For a mixed hyperbolic-parabolic system, one can specify source terms which depend\n",
    "on the gradient of the solution. Here, we solve a steady advection-diffusion\n",
    "equation with both solution and gradient-dependent source terms. The exact solution\n",
    "`u(x) = sin(x)` is given by `initial_condition`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "initial_condition = function (x, t, equations::LinearScalarAdvectionEquation1D)\n",
    "    return SVector(sin(x[1]))\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "For standard hyperbolic source terms, we pass in the solution, the spatial coordinate,\n",
    "the current time, and the hyperbolic equations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "source_terms = function (u, x, t, equations::LinearScalarAdvectionEquation1D)\n",
    "    f = a * cos(x[1]) + nu * sin(x[1]) - beta * (cos(x[1])^2)\n",
    "    return SVector(f)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "For gradient-dependent source terms, we also pass the solution gradients and the\n",
    "parabolic equations instead of the hyperbolic equations. Note that all parabolic\n",
    "equations have `equations_hyperbolic` as a solution field.\n",
    "\n",
    "For advection-diffusion, the gradients are gradients of the solution `u`. However,\n",
    "for systems such as `CompressibleNavierStokesDiffusion1D`, different gradient\n",
    "variables can be selected through the `gradient_variables` keyword option. The\n",
    "choice of `gradient_variables` will also determine the variables whose gradients\n",
    "are passed into `source_terms_parabolic`.\n",
    "\n",
    "The `gradients` passed to the `source_terms_parabolic` are a tuple of vectors;\n",
    "`gradients[1]` are the gradients in the first coordinate direction,\n",
    "`gradients[1][1]` is the gradient of the first (and only in this case) variable\n",
    "in the first coordinate direction."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "source_terms_parabolic = function (u, gradients, x, t, equations::LaplaceDiffusion1D)\n",
    "    dudx = gradients[1][1]\n",
    "    return SVector(beta * dudx^2)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The parabolic source terms can then be supplied to `SemidiscretizationHyperbolicParabolic`\n",
    "by setting the optional keyword argument `source_terms_parabolic`.\n",
    "The rest of the code is identical to standard hyperbolic-parabolic cases. We create a\n",
    "system of ODEs through `semidiscretize`, define callbacks, and then passing the system\n",
    "to OrdinaryDiffEq.jl.\n",
    "\n",
    "Note that for this problem, since viscosity `nu` is relatively large, we utilize\n",
    "`ViscousFormulationLocalDG` instead of the default `ViscousFormulationBassiRebay1`\n",
    "parabolic solver, since the Bassi-Rebay 1 formulation is not accurate when the\n",
    "diffusivity is large relative to the mesh size."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "mesh = TreeMesh(-Float64(pi), Float64(pi);\n",
    "                initial_refinement_level = 4,\n",
    "                n_cells_max = 30_000,\n",
    "                periodicity = true)\n",
    "\n",
    "boundary_conditions = boundary_condition_periodic\n",
    "boundary_conditions_parabolic = boundary_condition_periodic\n",
    "\n",
    "semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),\n",
    "                                             initial_condition, solver;\n",
    "                                             solver_parabolic = ViscousFormulationLocalDG(),\n",
    "                                             source_terms = source_terms,\n",
    "                                             source_terms_parabolic = source_terms_parabolic,\n",
    "                                             boundary_conditions = (boundary_conditions,\n",
    "                                                                    boundary_conditions_parabolic))\n",
    "\n",
    "tspan = (0.0, 1.0)\n",
    "ode = semidiscretize(semi, tspan)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Finally, we note that while we solve the ODE system using explicit time-stepping, the maximum\n",
    "stable time-step is $O(h^2)$ due to the dominant parabolic term. We enforce this more stringent\n",
    "parabolic CFL condition using a diffusion-aware `StepsizeCallback`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "cfl_advective = 0.5\n",
    "cfl_diffusive = 0.05\n",
    "stepsize_callback = StepsizeCallback(cfl = cfl_advective,\n",
    "                                     cfl_diffusive = cfl_diffusive)\n",
    "callbacks = CallbackSet(SummaryCallback(), stepsize_callback)\n",
    "sol = solve(ode, RDPK3SpFSAL35(); adaptive = false, dt = stepsize_callback(ode),\n",
    "            ode_default_options()..., callback = callbacks)"
   ],
   "metadata": {},
   "execution_count": null
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.10"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.10",
   "language": "julia"
  }
 },
 "nbformat": 4
}