{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Tutorial 18: Transient nonlinear equation"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"In this tutorial, we will learn\n",
"* How to write a nonlinear transient weak form in Gridap\n",
"* How to setup a time-marching scheme for a nonlinear ODE"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We assume that the reader is familiar with Gridap's API for linear transient PDEs, introduced in Tutorial 17. We focus here on more advanced features of the ODE module of Gridap, applied to a nonlinear time-dependent PDE."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Problem statement"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We consider the same problem as in Tutorial 17, and use the same notations: find $u(t): \\Omega \\to \\mathbb{R}$ such that\n",
"$$\n",
"\\left\\lbrace\n",
"\\begin{aligned}\n",
"\\rho(t, x) c(t, x) \\partial_{t} u(t, x) - \\nabla \\cdot (k(t, x) \\nabla u(t, x)) &= q(t, x) & \\text{ in } \\Omega, \\\\\n",
"u(t, x) &= g(t, x) & \\text{ on } \\partial \\Omega, \\\\\n",
"u(t_{0}, x) &= u_{0}(x) & \\text{ in } \\Omega \\\\\n",
"\\end{aligned}\n",
"\\right.\n",
"$$\n",
"In this tutorial we consider a nonlinear (quadratic) conductivity coefficient $\\alpha(t, x, u) = \\alpha_{0}(t, x) + \\alpha_{1}(t, x) u(t, x) + \\alpha_{2}(t, x) u(t, x)^{2}$. Here again, we assume that the $\\alpha_{i}$ are continuous in time. The weak form of the problem reads: find $u(t) \\in U_{g}(t)$ such that $b(t, u, v) = \\ell(t, v)$ for all $t \\geq t_{0}$ and $v \\in V_{0}$, where the time-dependent bilinear and linear forms $b(t, \\cdot, \\cdot)$ and $\\ell(t, \\cdot)$ are defined as\n",
"$$\n",
"\\begin{aligned}\n",
"b(t, u, v) &= m(t, u, v) + a(t, u, v), \\\\\n",
"m(t, u, v) &= \\int_{\\Omega} v \\partial_{t} u(t) \\ {\\rm d} \\Omega, \\\\\n",
"a(t, u, v) &= \\int_{\\Omega} \\nabla v \\cdot [(\\alpha_{0}(t) + \\alpha_{1}(t) u(t) + \\alpha_{2}(t) u(t)^{2}) \\nabla u(t)] \\ {\\rm d} \\Omega, \\\\\n",
"\\ell(t, v) &= \\int_{\\Omega} v f(t) \\ {\\rm d} \\Omega,\n",
"\\end{aligned}\n",
"$$\n",
"and the the functional spaces are $U_{g}(t) = \\{u \\in H^{1}_{g(t)}(\\Omega), u \\nabla u \\in \\boldsymbol{L}^{2}(\\Omega), u^{2} \\nabla u \\in \\boldsymbol{L}^{2}(\\Omega)\\}$ and $V_{0} = H^{1}_{0}(\\Omega)$. In addition to the regularity conditions of Tutorial 17 on $f$, $g$ and $u_{0}$, we assume that for all $t \\geq t_{0}$, it holds $\\alpha_{i}(t) \\in L^{\\infty}(\\Omega)$ and $(x, X) \\mapsto \\alpha_{0}(t, x) + \\alpha_{1}(t, x) X + \\alpha_{2}(t, x) X^{2}$ is uniformly positive in $\\Omega \\times \\mathbb{R}$, i.e. $\\alpha_{2}(t)$ and $4 \\alpha_{0}(t) \\alpha_{2}(t) - \\alpha_{1}^{2}(t)$ are uniformly positive."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Discrete model, FE spaces, triangulation and quadrature"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We consider the same mesh, FE spaces, triangulation and quadrature as in Tutorial 17:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Gridap\n",
"domain = (-1, +1, -1, +1)\n",
"partition = (20, 20)\n",
"model = CartesianDiscreteModel(domain, partition)\n",
"\n",
"order = 1\n",
"reffe = ReferenceFE(lagrangian, Float64, order)\n",
"\n",
"V0 = TestFESpace(model, reffe, dirichlet_tags=\"boundary\")\n",
"\n",
"g(t) = x -> exp(-2 * t) * sinpi(t * x[1]) * (x[2]^2 - 1)\n",
"Ug = TransientTrialFESpace(V0, g)\n",
"\n",
"degree = 2\n",
"Ω = Triangulation(model)\n",
"dΩ = Measure(Ω, degree)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## Nonlinear weak form\n",
"We define the diffusion coefficients $\\alpha$ and $\\beta$, the total nonlinear diffusion coefficient $\\kappa$ as well as the forcing term $f$."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"α₀(t) = x -> 1 + sin(t) * (x[1]^2 + x[2]^2) / 4\n",
"α₁(t) = x -> cos(t) * x[1]^2 / 2\n",
"α₂(t) = x -> 1 + t * (x[1]^2 + x[2]^2)\n",
"α(t, u) = α₀(t) + α₁(t) * u + α₂(t) * u * u\n",
"f(t) = x -> sin(t) * sinpi(x[1]) * sinpi(x[2])"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"We now write the nonlinear weak form. Similar to steady nonlinear problems, we provide the residual and its Jacobian, here with respect to $u$ and $\\partial_{t} u$. The mass, stiffness and forcing terms are written as follows."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"m(t, u, v) = ∫(v * ∂t(u))dΩ\n",
"a(t, u, v) = ∫(∇(v) ⋅ (α(t, u) * ∇(u)))dΩ\n",
"l(t, v) = ∫(v * f(t))dΩ"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The Jacobians of the mass and the stiffness are"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"jac_m(t, u, dtu, v) = ∫(v * dtu)dΩ\n",
"jac_α(t, u, du) = α₁(t) * du + α₂(t) * (2 * u * du)\n",
"jac_a(t, u, du, v) = ∫(∇(v) ⋅ (α(t, u) * ∇(du)))dΩ + ∫(∇(v) ⋅ (jac_α(t, u, du) * ∇(u)))dΩ"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"We can now write the residual and its Jacobians with respect to $u$ and $\\partial_{t} u$ as follows"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"res(t, u, v) = m(t, u, v) + a(t, u, v) - l(t, v)\n",
"jac(t, u, du, v) = jac_a(t, u, du, v)\n",
"jac_t(t, u, dtu, v) = jac_m(t, u, dtu, v)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The most general way of constructing a transient FE operator is by using the `TransientFEOperator` constructor, which receives a residual, a Jacobian with respect to the unknown and a Jacobian with respect to the time derivative."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"op = TransientFEOperator(res, (jac, jac_t), Ug, V0)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"In this example, the mass term is linear so this ODE belongs to the class of quasilinear ODEs. We can indicate this additional structure to Gridap as follows"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"mass_ql(t, u, dtu, v) = ∫(dtu * v)dΩ\n",
"res_ql(t, u, v) = a(t, u, v) - l(t, v)\n",
"jac_ql(t, u, du, v) = jac_a(t, u, du, v)\n",
"jac_t_ql(t, u, dtu, v) = jac_m(t, u, dtu, v)\n",
"op_ql = TransientQuasilinearFEOperator(mass_ql, res_ql, (jac_ql, jac_t_ql), Ug, V0)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"In fact, this ODE further classifies as semilinear because its mass term does not involve $u$. In our case, the mass term is also constant in time, so the optimal operator is as follows. Note that the signature of the mass term does not involve `u` anymore, as this is the condition for an ODE to be semilinear."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"mass_sl(t, dtu, v) = ∫(dtu * v)dΩ\n",
"res_sl(t, u, v) = a(t, u, v) - l(t, v)\n",
"jac_sl(t, u, du, v) = jac_a(t, u, du, v)\n",
"jac_t_sl(t, u, dtu, v) = mass_sl(t, dtu, v)\n",
"op_sl = TransientSemilinearFEOperator(\n",
" mass_sl, res_sl, (jac_sl, jac_t_sl),\n",
" Ug, V0, constant_mass=true\n",
")"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"In all cases above, it is also possible to take advantage of automatic differentiation techniques to compute both Jacobians and build the transient FE operator from the residual and the FE spaces only."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Transient solver"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"We proceed to the definition of the ODE solver. If the ODE is described via a general nonlinear FE operator, we will need to provide these schemes with a nonlinear solver for systems of equations. If the operator is quasilinear and the scheme is explicit, one only needs a linear solver. Here we draw from `NLSolvers.jl` and rely on a Newton-Raphson solver based on Gridap's `LUSolver`."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"For example, for the `ThetaMethod`, one would write"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"lin_solver = LUSolver()\n",
"nl_solver = NLSolver(lin_solver, method=:newton, iterations=10, show_trace=false)\n",
"\n",
"Δt = 0.05\n",
"θ = 0.5\n",
"solver = ThetaMethod(nl_solver, Δt, θ)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"For a two-stage singly-diagonally-implicit scheme (of order 2), it would be"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"tableau = :SDIRK_2_2\n",
"solver_rk = RungeKutta(nl_solver, lin_solver, Δt, tableau)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"We define the initial condition and the solution using the `solve` function as in Tutorial 17:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"t0, tF = 0.0, 10.0\n",
"uh0 = interpolate_everywhere(g(t0), Ug(t0))\n",
"uh = solve(solver, op_sl, t0, tF, uh0)"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"## Postprocessing"
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"Here again, we export the solution at each time step as follows"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"if !isdir(\"tmp_nl\")\n",
" mkdir(\"tmp_nl\")\n",
"end\n",
"\n",
"createpvd(\"results_nl\") do pvd\n",
" pvd[0] = createvtk(Ω, \"tmp_nl/results_0\" * \".vtu\", cellfields=[\"u\" => uh0])\n",
" for (tn, uhn) in uh\n",
" pvd[tn] = createvtk(Ω, \"tmp_nl/results_$tn\" * \".vtu\", cellfields=[\"u\" => uhn])\n",
" end\n",
"end"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
""
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"---\n",
"\n",
"*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.10.8"
},
"kernelspec": {
"name": "julia-1.10",
"display_name": "Julia 1.10.8",
"language": "julia"
}
},
"nbformat": 4
}