Corrections for heat equation tutorial

src/heat_equation/heat_equation.jl

@@ -4,43 +4,42 @@ println("Running heat equation tutorial...") #hide
 # In this tutorial, the heat equation (first steady and then unsteady) is solved using finite-elements.
 #
 # # Theory
-# This example shows how to solve the heat equation with eventually variable physical properties in steady and unsteady formulations:
+# This example shows how to solve the heat equation in steady and unsteady formulations. The unsteady heat equation is given by:
 # ```math
-#   \rho C_p \partial_t u - \nabla . ( \lambda u) = f
+#   \rho C_p \frac{\partial u}{\partial t} - \nabla \cdot ( \lambda \nabla u) = f
 # ```
 # We shall assume that $f, \, \rho, \, C_p, \, \lambda \, \in L^2(\Omega)$. The weak form of the problem is given by: find $ u \in \tilde{H}^1_0(\Omega)$
 # (there will be at least one Dirichlet boundary condition) such that:
 # ```math
-#   \forall v \in  \tilde{H}^1_0(\Omega), \, \, \, \underbrace{\int_\Omega \partial_t u . v dx}_{m(\partial_t u,v)} + \underbrace{\int_\Omega \nabla u . \nabla v dx}_{a(u,v)} = \underbrace{\int_\Omega f v dx}_{l(v)}
+#   \forall v \in  \tilde{H}^1_0(\Omega), \, \, \, \underbrace{\int_\Omega \rho C_p \frac{\partial u}{\partial t}  v dx}_{m(\partial_t u,v)} + \underbrace{\int_\Omega \lambda \nabla u \cdot \nabla v dx}_{a(u,v)} = \underbrace{\int_\Omega f v dx}_{l(v)}
 # ```
-# To numerically solve this problem we seek an approximate solution using Lagrange $$P^1$$ or $$P^2$$ elements.
-# Here we assume that the domain can be split into two domains having different material properties.
-
-# # Steady case
+# To numerically solve this problem we seek an approximate solution using Lagrange $$P^2$$ elements. 
+# 
 # As usual, start by importing the necessary packages.
 using Bcube
 using BcubeGmsh
 using BcubeVTK
 using LinearAlgebra
+using Test #src
 
 const is_tested = get(ENV, "TestMode", "false") == "true" #src
 if is_tested                                              #src
-    using Test
     import ..Tester: test_ref                             #src
 end                                                       #src
 
 # First we define some physical and numerical constants
-const htc = 100.0 # Heat transfer coefficient (bnd cdt)
-const Tr = 268.0 # Recovery temperature (bnd cdt)
-const phi = 100.0
-const q = 1500.0
-const λ = 100.0
-const η = λ
-const ρCp = 100.0 * 200.0
-const degree = 2
+const q = 1500.0 # heat source
+const λ = 100.0 # thermal conductivity
+const ρCp = 100.0 * 200.0 # density times specific heat capacity
+const degree = 2 # Degree of the Lagrange polynomials (for $$P^2$$ elements)
 const outputpath = joinpath(@__DIR__, "..", "..", "myout", "heat_equation/")
 mkpath(outputpath) #hide
 
+# # Steady case
+
+# We will first start by solving the steady case: find $u \in \tilde{H}^1_0(\Omega)$ such that $\forall v \in  \tilde{H}^1_0(\Omega), \, \, \, a(u,v) = l(v)$.
+
+println(" ----- Solving steady problem -----")
 # Read 2D mesh
 mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri.msh")
 mesh = read_mesh(mesh_path)
@@ -55,41 +54,46 @@ V = TestFESpace(U)
 # Define measures for cell integration
 dΩ = Measure(CellDomain(mesh), 2 * degree + 1)
 
-# Define bilinear and linear forms
-a(u, v) = ∫(η * ∇(u) ⋅ ∇(v))dΩ
+# Define bilinear and linear forms. The steady problem is defined by the bilinear form a
+# and the linear form l.
+a(u, v) = ∫(λ * ∇(u) ⋅ ∇(v))dΩ
 l(v) = ∫(q * v)dΩ
 
 # Create an affine FE system and solve it using the `AffineFESystem` structure.
-# The package `LinearSolve` is used behind the scenes, so different solver may
-# be used to invert the system (ex: `solve(...; alg = IterativeSolversJL_GMRES())`)
-# The result is a FEFunction (`ϕ`).
-# We can interpolate it on mesh centers : the result is named `Tcn`.
-sys = AffineFESystem(a, l, U, V)
-ϕ = Bcube.solve(sys)
-Tcn = var_on_centers(ϕ, mesh)
-
-# Compute analytical solution for comparison. Apply the analytical solution
-# on mesh centers
-T_analytical = x -> 260.0 + (q / λ) * x[1] * (1.0 - 0.5 * x[1])
-Tca = map(T_analytical, get_cell_centers(mesh))
-
-# Write both the obtained FE solution and the analytical solution to a vtk file. To
-# write the data on mesh centers, we need to wrap them in a `MeshCellData` object.
-using BcubeVTK
+# By default, an LU decomposition is used to solve the system. In the present case
+# we know that the matrix A is symmetric. Hence a Cholesky decomposition can be used.
+# The result is a FEFunction (`Tn`).
+# We can extract its dof values: the result is named `Tn_dofs`.
+Cholesky_linsolve!(y, A, x) = y .= cholesky(Symmetric(A)) \ x
+sys = AffineFESystem(a, l, U, V, Cholesky_linsolve!)
+Tn = Bcube.solve(sys)
+Tn_dofs = get_dof_values(Tn)
+
+# Compute analytical solution for comparison. An FEFunction is built using
+# the analytical solution and the dof values are extracted. 
+T_analytical = PhysicalFunction(x -> 260.0 + (q / λ) * x[1] * (1.0 - 0.5 * x[1]))
+Ta = FEFunction(U, mesh, T_analytical)
+Ta_dofs = get_dof_values(Ta)
+
+# Write both the obtained FE solution and the analytical solution to a vtk file. 
 mkpath(outputpath)
-dict_vars = Dict("Temperature" => MeshCellData(Tcn), "Temperature_a" => MeshCellData(Tca))
+dict_vars = Dict("Temperature (numerical)" => Tn, "Temperature (analytical)" => Ta)
 write_file(outputpath * "result_steady_heat_equation.pvd", mesh, dict_vars)
 
-# Compute and display the error
-@show norm(Tcn .- Tca, Inf) / norm(Tca, Inf)
+# Compute and display the error, which is computed using the previously extracted dof values.
+@show norm(Tn_dofs .- Ta_dofs, Inf) / norm(Ta_dofs, Inf)
 
-if is_tested                                             #src
-    @test norm(Tcn .- Tca, Inf) / norm(Tca, Inf) < 2e-14 #src
-end                                                      #src
+if is_tested                                                         #src
+    @test norm(Tn_dofs .- Ta_dofs, Inf) / norm(Ta_dofs, Inf) < 2e-14 #src
+end                                                                  #src
 
 # # Unsteady case
+# We now solve the unsteady problem: find $ u $ such that for all $v , \, \, \, m(\frac{\partial u}{\partial t}, v) + a(u,v) = l(v)$.
 # The code for the unsteady case if of course very similar to the steady case, at least for the
-# beginning. Start by defining two additional parameters:
+# beginning. 
+
+println(" ----- Solving unsteady problem -----")
+# Start by defining two additional parameters:
 totalTime = 100.0
 Δt = 0.1
 
@@ -103,10 +107,10 @@ U = TrialFESpace(fs, mesh, Dict("West" => 260.0))
 V = TestFESpace(U)
 dΩ = Measure(CellDomain(mesh), 2 * degree + 1)
 
-# Compute matrices associated to bilinear and linear forms, and assemble
-a(u, v) = ∫(η * ∇(u) ⋅ ∇(v))dΩ
+# Compute matrices associated to bilinear and linear forms, and assemble.
+# The bilinear form a and linear form l have already been defined in the
+# steady problem. We now define the bilinear form m.
 m(u, v) = ∫(ρCp * u ⋅ v)dΩ
-l(v) = ∫(q * v)dΩ
 
 A = assemble_bilinear(a, U, V)
 M = assemble_bilinear(m, U, V)
@@ -134,7 +138,7 @@ Miter = factorize(M + Δt * A)
 
 # Write initial solution to a file
 mkpath(outputpath)
-dict_vars = Dict("Temperature" => MeshCellData(var_on_centers(ϕ, mesh)))
+dict_vars = Dict("Temperature" => ϕ)
 write_file(outputpath * "result_unsteady_heat_equation.pvd", mesh, dict_vars, 0, 0.0)
 
 # Time loop
@@ -144,21 +148,21 @@ while t <= totalTime
     global t, itime
     t += Δt
     itime = itime + 1
-    #! format: off
-    if !is_tested #src
-    @show t, itime
-    end #src
-    #! format: on
+        #! format: off
+        if !is_tested #src
+        @show t, itime
+        end #src
+        #! format: on
 
-    ## Compute rhs
+        ## Compute rhs
     rhs = Δt * L + M * (get_dof_values(ϕ) .- Ud)
 
     ## Invert system and apply inverse shift
     set_dof_values!(ϕ, Miter \ rhs .+ Ud)
 
     ## Write solution (every 10 iterations)
     if itime % 10 == 0
-        dict_vars = Dict("Temperature" => MeshCellData(var_on_centers(ϕ, mesh)))
+        dict_vars = Dict("Temperature" => ϕ)
         write_file(
             outputpath * "result_unsteady_heat_equation.pvd",
             mesh,

src/helmholtz/helmholtz.jl

@@ -109,7 +109,7 @@ write_file(joinpath(@__DIR__, outputvtk), mesh, dict_vars)                  #md
 # ![](../assets/helmholtz_x21_y21_vp6.png)
 
 if get(ENV, "TestMode", "false") == "true"                      #src
-    import ..Tester: test_ref                                   #src
+    import ..Tester: test_ref, compare                          #src
     results = sqrt.(abs.(vp[3:8]))                              #src
     ref_results = [                                             #src
         3.144823462554393,                                      #src
@@ -122,6 +122,10 @@ if get(ENV, "TestMode", "false") == "true"                      #src
     @test all(results .≈ ref_results)                           #src
     test_ref("helmholtz_A.jld2", A)                             #src
     test_ref("helmholtz_B.jld2", B)                             #src
-    test_ref("helmholtz_vp.jld2", vp)                           #src
+    test_ref(                                                   #src
+        "helmholtz_vp.jld2",                                    #src
+        vp,                                                     #src
+        (a, b) -> compare(a, b; atol = 1.0e-11, rtol = 1.0e-11),#src
+    )                                                           #src
 end                                                             #src
 end #hide