{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Transient heat equation\n",
"\n",
"\n",
"\n",
"\n",
"*Figure 1*: Visualization of the temperature time evolution on a unit\n",
"square where the prescribed temperature on the upper and lower parts\n",
"of the boundary increase with time."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Introduction\n",
"\n",
"In this example we extend the heat equation by a time dependent term, i.e.\n",
"$$\n",
" \\frac{\\partial u}{\\partial t}-\\nabla \\cdot (k \\nabla u) = f \\quad x \\in \\Omega,\n",
"$$\n",
"\n",
"where $u$ is the unknown temperature field, $k$ the heat conductivity,\n",
"$f$ the heat source and $\\Omega$ the domain. For simplicity, we hard code $f = 0.1$\n",
"and $k = 10^{-3}$. We define homogeneous Dirichlet boundary conditions along the left and right edge of the domain.\n",
"$$\n",
"u(x,t) = 0 \\quad x \\in \\partial \\Omega_1,\n",
"$$\n",
"where $\\partial \\Omega_1$ denotes the left and right boundary of $\\Omega$.\n",
"\n",
"Further, we define heterogeneous Dirichlet boundary conditions at the top and bottom edge $\\partial \\Omega_2$.\n",
"We choose a linearly increasing function $a(t)$ that describes the temperature at this boundary\n",
"$$\n",
"u(x,t) = a(t) \\quad x \\in \\partial \\Omega_2.\n",
"$$\n",
"The semidiscrete weak form is given by\n",
"$$\n",
"\\int_{\\Omega}v \\frac{\\partial u}{\\partial t} \\ \\mathrm{d}\\Omega + \\int_{\\Omega} k \\nabla v \\cdot \\nabla u \\ \\mathrm{d}\\Omega = \\int_{\\Omega} f v \\ \\mathrm{d}\\Omega,\n",
"$$\n",
"where $v$ is a suitable test function. Now, we still need to discretize the time derivative. An implicit Euler scheme is applied,\n",
"which yields:\n",
"$$\n",
"\\int_{\\Omega} v\\, u_{n+1}\\ \\mathrm{d}\\Omega + \\Delta t\\int_{\\Omega} k \\nabla v \\cdot \\nabla u_{n+1} \\ \\mathrm{d}\\Omega = \\Delta t\\int_{\\Omega} f v \\ \\mathrm{d}\\Omega + \\int_{\\Omega} v \\, u_{n} \\ \\mathrm{d}\\Omega.\n",
"$$\n",
"If we assemble the discrete operators, we get the following algebraic system:\n",
"$$\n",
"\\mathbf{M} \\mathbf{u}_{n+1} + Δt \\mathbf{K} \\mathbf{u}_{n+1} = Δt \\mathbf{f} + \\mathbf{M} \\mathbf{u}_{n}\n",
"$$\n",
"In this example we apply the boundary conditions to the assembled discrete operators (mass matrix $\\mathbf{M}$ and stiffnes matrix $\\mathbf{K}$)\n",
"only once. We utilize the fact that in finite element computations Dirichlet conditions can be applied by\n",
"zero out rows and columns that correspond\n",
"to a prescribed dof in the system matrix ($\\mathbf{A} = Δt \\mathbf{K} + \\mathbf{M}$) and setting the value of the right-hand side vector to the value\n",
"of the Dirichlet condition. Thus, we only need to apply in every time step the Dirichlet condition to the right-hand side of the problem."
],
"metadata": {}
},
{
"cell_type": "markdown",
"source": [
"## Commented Program\n",
"\n",
"Now we solve the problem in Ferrite. What follows is a program spliced with comments.\n",
"\n",
"First we load Ferrite, and some other packages we need."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"using Ferrite, SparseArrays, WriteVTK"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"We create the same grid as in the heat equation example."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"grid = generate_grid(Quadrilateral, (100, 100));"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"### Trial and test functions\n",
"Again, we define the structs that are responsible for the `shape_value` and `shape_gradient` evaluation."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"ip = Lagrange{RefQuadrilateral, 1}()\n",
"qr = QuadratureRule{RefQuadrilateral}(2)\n",
"cellvalues = CellValues(qr, ip);"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"### Degrees of freedom\n",
"After this, we can define the `DofHandler` and distribute the DOFs of the problem."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"dh = DofHandler(grid)\n",
"add!(dh, :u, ip)\n",
"close!(dh);"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"By means of the `DofHandler` we can allocate the needed `SparseMatrixCSC`.\n",
"`M` refers here to the so called mass matrix, which always occurs in time related terms, i.e.\n",
"$$\n",
"M_{ij} = \\int_{\\Omega} v_i \\, u_j \\ \\mathrm{d}\\Omega,\n",
"$$\n",
"where $u_i$ and $v_j$ are trial and test functions, respectively."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"K = allocate_matrix(dh);\n",
"M = allocate_matrix(dh);"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"We also preallocate the right hand side"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"f = zeros(ndofs(dh));"
],
"metadata": {},
"execution_count": 6
},
{
"cell_type": "markdown",
"source": [
"### Boundary conditions\n",
"In order to define the time dependent problem, we need some end time `T` and something that describes\n",
"the linearly increasing Dirichlet boundary condition on $\\partial \\Omega_2$."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"max_temp = 100\n",
"Δt = 1\n",
"T = 200\n",
"t_rise = 100\n",
"ch = ConstraintHandler(dh);"
],
"metadata": {},
"execution_count": 7
},
{
"cell_type": "markdown",
"source": [
"Here, we define the boundary condition related to $\\partial \\Omega_1$."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"∂Ω₁ = union(getfacetset.((grid,), [\"left\", \"right\"])...)\n",
"dbc = Dirichlet(:u, ∂Ω₁, (x, t) -> 0)\n",
"add!(ch, dbc);"
],
"metadata": {},
"execution_count": 8
},
{
"cell_type": "markdown",
"source": [
"While the next code block corresponds to the linearly increasing temperature description on $\\partial \\Omega_2$\n",
"until `t=t_rise`, and then keep constant"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"∂Ω₂ = union(getfacetset.((grid,), [\"top\", \"bottom\"])...)\n",
"dbc = Dirichlet(:u, ∂Ω₂, (x, t) -> max_temp * clamp(t / t_rise, 0, 1))\n",
"add!(ch, dbc)\n",
"close!(ch)\n",
"update!(ch, 0.0);"
],
"metadata": {},
"execution_count": 9
},
{
"cell_type": "markdown",
"source": [
"### Assembling the linear system\n",
"As in the heat equation example we define a `doassemble!` function that assembles the diffusion parts of the equation:"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "doassemble_K! (generic function with 1 method)"
},
"metadata": {},
"execution_count": 10
}
],
"cell_type": "code",
"source": [
"function doassemble_K!(K::SparseMatrixCSC, f::Vector, cellvalues::CellValues, dh::DofHandler)\n",
"\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
" Ke = zeros(n_basefuncs, n_basefuncs)\n",
" fe = zeros(n_basefuncs)\n",
"\n",
" assembler = start_assemble(K, f)\n",
"\n",
" for cell in CellIterator(dh)\n",
"\n",
" fill!(Ke, 0)\n",
" fill!(fe, 0)\n",
"\n",
" reinit!(cellvalues, cell)\n",
"\n",
" for q_point in 1:getnquadpoints(cellvalues)\n",
" dΩ = getdetJdV(cellvalues, q_point)\n",
"\n",
" for i in 1:n_basefuncs\n",
" v = shape_value(cellvalues, q_point, i)\n",
" ∇v = shape_gradient(cellvalues, q_point, i)\n",
" fe[i] += 0.1 * v * dΩ\n",
" for j in 1:n_basefuncs\n",
" ∇u = shape_gradient(cellvalues, q_point, j)\n",
" Ke[i, j] += 1.0e-3 * (∇v ⋅ ∇u) * dΩ\n",
" end\n",
" end\n",
" end\n",
"\n",
" assemble!(assembler, celldofs(cell), Ke, fe)\n",
" end\n",
" return K, f\n",
"end"
],
"metadata": {},
"execution_count": 10
},
{
"cell_type": "markdown",
"source": [
"In addition to the diffusive part, we also need a function that assembles the mass matrix `M`."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "doassemble_M! (generic function with 1 method)"
},
"metadata": {},
"execution_count": 11
}
],
"cell_type": "code",
"source": [
"function doassemble_M!(M::SparseMatrixCSC, cellvalues::CellValues, dh::DofHandler)\n",
"\n",
" n_basefuncs = getnbasefunctions(cellvalues)\n",
" Me = zeros(n_basefuncs, n_basefuncs)\n",
"\n",
" assembler = start_assemble(M)\n",
"\n",
" for cell in CellIterator(dh)\n",
"\n",
" fill!(Me, 0)\n",
"\n",
" reinit!(cellvalues, cell)\n",
"\n",
" for q_point in 1:getnquadpoints(cellvalues)\n",
" dΩ = getdetJdV(cellvalues, q_point)\n",
"\n",
" for i in 1:n_basefuncs\n",
" v = shape_value(cellvalues, q_point, i)\n",
" for j in 1:n_basefuncs\n",
" u = shape_value(cellvalues, q_point, j)\n",
" Me[i, j] += (v * u) * dΩ\n",
" end\n",
" end\n",
" end\n",
"\n",
" assemble!(assembler, celldofs(cell), Me)\n",
" end\n",
" return M\n",
"end"
],
"metadata": {},
"execution_count": 11
},
{
"cell_type": "markdown",
"source": [
"### Solution of the system\n",
"We first assemble all parts in the prior allocated `SparseMatrixCSC`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"K, f = doassemble_K!(K, f, cellvalues, dh)\n",
"M = doassemble_M!(M, cellvalues, dh)\n",
"A = (Δt .* K) + M;"
],
"metadata": {},
"execution_count": 12
},
{
"cell_type": "markdown",
"source": [
"Now, we need to save all boundary condition related values of the unaltered system matrix `A`, which is done\n",
"by `get_rhs_data`. The function returns a `RHSData` struct, which contains all needed information to apply\n",
"the boundary conditions solely on the right-hand-side vector of the problem."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"rhsdata = get_rhs_data(ch, A);"
],
"metadata": {},
"execution_count": 13
},
{
"cell_type": "markdown",
"source": [
"We set the values at initial time step, denoted by uₙ, to a bubble-shape described by\n",
"$(x_1^2-1)(x_2^2-1)$, such that it is zero at the boundaries and the maximum temperature in the center."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"uₙ = zeros(length(f));\n",
"apply_analytical!(uₙ, dh, :u, x -> (x[1]^2 - 1) * (x[2]^2 - 1) * max_temp);"
],
"metadata": {},
"execution_count": 14
},
{
"cell_type": "markdown",
"source": [
"Here, we apply **once** the boundary conditions to the system matrix `A`."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"apply!(A, ch);"
],
"metadata": {},
"execution_count": 15
},
{
"cell_type": "markdown",
"source": [
"To store the solution, we initialize a paraview collection (.pvd) file,"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "VTKGridFile for the closed file \"transient-heat-0.vtu\"."
},
"metadata": {},
"execution_count": 16
}
],
"cell_type": "code",
"source": [
"pvd = paraview_collection(\"transient-heat\")\n",
"VTKGridFile(\"transient-heat-0\", dh) do vtk\n",
" write_solution(vtk, dh, uₙ)\n",
" pvd[0.0] = vtk\n",
"end"
],
"metadata": {},
"execution_count": 16
},
{
"cell_type": "markdown",
"source": [
"At this point everything is set up and we can finally approach the time loop."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"for (step, t) in enumerate(Δt:Δt:T)\n",
" #First of all, we need to update the Dirichlet boundary condition values.\n",
" update!(ch, t)\n",
"\n",
" #Secondly, we compute the right-hand-side of the problem.\n",
" b = Δt .* f .+ M * uₙ\n",
" #Then, we can apply the boundary conditions of the current time step.\n",
" apply_rhs!(rhsdata, b, ch)\n",
"\n",
" #Finally, we can solve the time step and save the solution afterwards.\n",
" u = A \\ b\n",
"\n",
" VTKGridFile(\"transient-heat-$step\", dh) do vtk\n",
" write_solution(vtk, dh, u)\n",
" pvd[t] = vtk\n",
" end\n",
" #At the end of the time loop, we set the previous solution to the current one and go to the next time step.\n",
" uₙ .= u\n",
"end"
],
"metadata": {},
"execution_count": 17
},
{
"cell_type": "markdown",
"source": [
"In order to use the .pvd file we need to store it to the disk, which is done by:"
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"vtk_save(pvd);"
],
"metadata": {},
"execution_count": 18
},
{
"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.11.3"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.3",
"language": "julia"
}
},
"nbformat": 4
}