From 07f713e8e527b2b6e8ac2b9e4c0cd3be6e91a204 Mon Sep 17 00:00:00 2001 From: lokmanb Date: Wed, 22 Apr 2026 21:56:55 +0200 Subject: [PATCH 1/8] typo correction in heat_equation.jl --- src/heat_equation/heat_equation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index 76143a68..023935db 100644 --- a/src/heat_equation/heat_equation.jl +++ b/src/heat_equation/heat_equation.jl @@ -11,7 +11,7 @@ println("Running heat equation tutorial...") #hide # 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 \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)} # ``` # 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. From 3856788a6ecfbc41be831bd08295ba3cb0accb06 Mon Sep 17 00:00:00 2001 From: lokmanb Date: Wed, 22 Apr 2026 22:30:28 +0200 Subject: [PATCH 2/8] Changing vtk output to use FEFunctions directly in heat_equation.jl --- src/heat_equation/heat_equation.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index 023935db..ef16e3d6 100644 --- a/src/heat_equation/heat_equation.jl +++ b/src/heat_equation/heat_equation.jl @@ -22,10 +22,10 @@ using Bcube using BcubeGmsh using BcubeVTK using LinearAlgebra +using Test const is_tested = get(ENV, "TestMode", "false") == "true" #src if is_tested #src - using Test import ..Tester: test_ref #src end #src @@ -65,27 +65,27 @@ l(v) = ∫(q * v)dΩ # 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) +Tn = Bcube.solve(sys) +Tn_dofs = get_dof_values(Tn) # 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)) +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. To # write the data on mesh centers, we need to wrap them in a `MeshCellData` object. -using BcubeVTK 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) +@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 # The code for the unsteady case if of course very similar to the steady case, at least for the @@ -134,7 +134,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 @@ -158,7 +158,7 @@ while t <= totalTime ## 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, From 1425ae567bbd81a650a44a92738da101b7dc8332 Mon Sep 17 00:00:00 2001 From: lokmanb Date: Wed, 22 Apr 2026 22:51:21 +0200 Subject: [PATCH 3/8] use Cholesky for linsolve! in heat_equation.jl --- src/heat_equation/heat_equation.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index ef16e3d6..9f001c9b 100644 --- a/src/heat_equation/heat_equation.jl +++ b/src/heat_equation/heat_equation.jl @@ -60,27 +60,27 @@ 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) +# 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) = ldiv!(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. Apply the analytical solution -# on mesh centers +# 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. To -# write the data on mesh centers, we need to wrap them in a `MeshCellData` object. +# Write both the obtained FE solution and the analytical solution to a vtk file. mkpath(outputpath) 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 +# 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 From 548225a9694754d74fc6f52259883a3596d683fb Mon Sep 17 00:00:00 2001 From: lokmanb Date: Fri, 24 Apr 2026 21:41:55 +0200 Subject: [PATCH 4/8] changing tolerance for check in helmholtz.jl --- src/helmholtz/helmholtz.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/helmholtz/helmholtz.jl b/src/helmholtz/helmholtz.jl index ed1fceff..43753ff8 100644 --- a/src/helmholtz/helmholtz.jl +++ b/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 From 0805ded06a97634e0da5ab6a9787337c3ad5195d Mon Sep 17 00:00:00 2001 From: lokmanb Date: Sun, 26 Apr 2026 15:27:03 +0200 Subject: [PATCH 5/8] Using functions to solve steady and unsteady problems in heat_equation.jl. Corrections to the documentation of the case --- src/heat_equation/heat_equation.jl | 288 +++++++++++++++-------------- 1 file changed, 147 insertions(+), 141 deletions(-) diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index 9f001c9b..75df0d9a 100644 --- a/src/heat_equation/heat_equation.jl +++ b/src/heat_equation/heat_equation.jl @@ -4,25 +4,24 @@ 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 \rho C_p \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. +# To numerically solve this problem we seek an approximate solution using Lagrange $$P^2$$ elements. -# # Steady case +# # As usual, start by importing the necessary packages. using Bcube using BcubeGmsh using BcubeVTK using LinearAlgebra -using Test +using Test const is_tested = get(ENV, "TestMode", "false") == "true" #src if is_tested #src @@ -30,148 +29,155 @@ if is_tested #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 -# Read 2D mesh -mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri.msh") -mesh = read_mesh(mesh_path) - -# Build function space and associated Trial and Test FE spaces. -# We impose a Dirichlet condition with a temperature of 260K -# on boundary "West" -fs = FunctionSpace(:Lagrange, degree) -U = TrialFESpace(fs, mesh, Dict("West" => 260.0)) -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Ω -l(v) = ∫(q * v)dΩ - -# Create an affine FE system and solve it using the `AffineFESystem` structure. -# 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) = ldiv!(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 (numerical)" => Tn, "Temperature (analytical)" => Ta) -write_file(outputpath * "result_steady_heat_equation.pvd", mesh, dict_vars) - -# 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(Tn_dofs .- Ta_dofs, Inf) / norm(Ta_dofs, Inf) < 2e-14 #src -end #src +# # Steady case +function solve_steady() + println(" ----- Solving steady problem -----") + # Read 2D mesh + mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri.msh") + mesh = read_mesh(mesh_path) + + # Build function space and associated Trial and Test FE spaces. + # We impose a Dirichlet condition with a temperature of 260K + # on boundary "West" + fs = FunctionSpace(:Lagrange, degree) + U = TrialFESpace(fs, mesh, Dict("West" => 260.0)) + 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Ω + l(v) = ∫(q * v)dΩ + + # Create an affine FE system and solve it using the `AffineFESystem` structure. + # 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) = ldiv!(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 (numerical)" => Tn, "Temperature (analytical)" => Ta) + write_file(outputpath * "result_steady_heat_equation.pvd", mesh, dict_vars) + + # 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(Tn_dofs .- Ta_dofs, Inf) / norm(Ta_dofs, Inf) < 2e-14 #src + end #src +end # # Unsteady case # 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: -totalTime = 100.0 -Δt = 0.1 - -# Read a slightly different mesh -mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri_2.msh") -mesh = read_mesh(mesh_path) - -# The rest is similar to the steady case -fs = FunctionSpace(:Lagrange, degree) -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Ω -m(u, v) = ∫(ρCp * u ⋅ v)dΩ -l(v) = ∫(q * v)dΩ - -A = assemble_bilinear(a, U, V) -M = assemble_bilinear(m, U, V) -L = assemble_linear(l, V) - -# Compute a vector of dofs whose values are zeros everywhere -# except on dofs lying on a Dirichlet boundary, where they -# take the Dirichlet value -Ud = assemble_dirichlet_vector(U, V, mesh) - -# Apply lift -L = L - A * Ud - -# Apply homogeneous dirichlet condition -apply_homogeneous_dirichlet_to_vector!(L, U, V, mesh) -apply_dirichlet_to_matrix!((A, M), U, V, mesh) - -# Form time iteration matrix -# (note that this is bad for performance since up to now, -# M and A are sparse matrices) -Miter = factorize(M + Δt * A) - -# Init the solution with a constant temperature of 260K -ϕ = FEFunction(U, 260.0) - -# Write initial solution to a file -mkpath(outputpath) -dict_vars = Dict("Temperature" => ϕ) -write_file(outputpath * "result_unsteady_heat_equation.pvd", mesh, dict_vars, 0, 0.0) - -# Time loop -itime = 0 -t = 0.0 -while t <= totalTime - global t, itime - t += Δt - itime = itime + 1 - #! format: off - if !is_tested #src - @show t, itime - end #src - #! format: on - - ## 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" => ϕ) - write_file( - outputpath * "result_unsteady_heat_equation.pvd", - mesh, - dict_vars, - itime, - t; - collection_append = true, - ) +# beginning. +function solve_unsteady() + println(" ----- Solving unsteady problem -----") + # Start by defining two additional parameters: + totalTime = 100.0 + Δt = 0.1 + + # Read a slightly different mesh + mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri_2.msh") + mesh = read_mesh(mesh_path) + + # The rest is similar to the steady case + fs = FunctionSpace(:Lagrange, degree) + 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Ω + m(u, v) = ∫(ρCp * u ⋅ v)dΩ + l(v) = ∫(q * v)dΩ + + A = assemble_bilinear(a, U, V) + M = assemble_bilinear(m, U, V) + L = assemble_linear(l, V) + + # Compute a vector of dofs whose values are zeros everywhere + # except on dofs lying on a Dirichlet boundary, where they + # take the Dirichlet value + Ud = assemble_dirichlet_vector(U, V, mesh) + + # Apply lift + L = L - A * Ud + + # Apply homogeneous dirichlet condition + apply_homogeneous_dirichlet_to_vector!(L, U, V, mesh) + apply_dirichlet_to_matrix!((A, M), U, V, mesh) + + # Form time iteration matrix + # (note that this is bad for performance since up to now, + # M and A are sparse matrices) + Miter = factorize(M + Δt * A) + + # Init the solution with a constant temperature of 260K + ϕ = FEFunction(U, 260.0) + + # Write initial solution to a file + mkpath(outputpath) + dict_vars = Dict("Temperature" => ϕ) + write_file(outputpath * "result_unsteady_heat_equation.pvd", mesh, dict_vars, 0, 0.0) + + # Time loop + itime = 0 + t = 0.0 + while t <= totalTime + t, itime + t += Δt + itime = itime + 1 + #! format: off + if !is_tested #src + @show t, itime + end #src + #! format: on + + ## 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" => ϕ) + write_file( + outputpath * "result_unsteady_heat_equation.pvd", + mesh, + dict_vars, + itime, + t; + collection_append = true, + ) + end end + + if is_tested #src + test_ref("heat_equation_100s.jld2", get_dof_values(ϕ)) #src + end #src end -if is_tested #src - test_ref("heat_equation_100s.jld2", get_dof_values(ϕ)) #src -end #src +solve_steady() +solve_unsteady() end #hide From 3740e2057b06fd42c854558b38d6d65d84e5bb2f Mon Sep 17 00:00:00 2001 From: lokmanb Date: Sun, 26 Apr 2026 15:38:38 +0200 Subject: [PATCH 6/8] slight change to hide test in doc for heat_equation.jl --- src/heat_equation/heat_equation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index 75df0d9a..c6bd03cb 100644 --- a/src/heat_equation/heat_equation.jl +++ b/src/heat_equation/heat_equation.jl @@ -21,7 +21,7 @@ using Bcube using BcubeGmsh using BcubeVTK using LinearAlgebra -using Test +using Test #src const is_tested = get(ENV, "TestMode", "false") == "true" #src if is_tested #src From 606d29dee150801b7de86dd03a3c6718c6077385 Mon Sep 17 00:00:00 2001 From: lokmanb Date: Tue, 2 Jun 2026 23:16:42 +0200 Subject: [PATCH 7/8] removing two function approach for the heat equation tutorial --- src/heat_equation/heat_equation.jl | 260 ++++++++++++++--------------- 1 file changed, 129 insertions(+), 131 deletions(-) diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index c6bd03cb..de984467 100644 --- a/src/heat_equation/heat_equation.jl +++ b/src/heat_equation/heat_equation.jl @@ -13,9 +13,8 @@ println("Running heat equation tutorial...") #hide # ```math # \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^2$$ elements. - -# +# 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 @@ -37,115 +36,118 @@ const outputpath = joinpath(@__DIR__, "..", "..", "myout", "heat_equation/") mkpath(outputpath) #hide # # Steady case -function solve_steady() - println(" ----- Solving steady problem -----") - # Read 2D mesh - mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri.msh") - mesh = read_mesh(mesh_path) - - # Build function space and associated Trial and Test FE spaces. - # We impose a Dirichlet condition with a temperature of 260K - # on boundary "West" - fs = FunctionSpace(:Lagrange, degree) - U = TrialFESpace(fs, mesh, Dict("West" => 260.0)) - 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Ω - l(v) = ∫(q * v)dΩ - - # Create an affine FE system and solve it using the `AffineFESystem` structure. - # 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) = ldiv!(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 (numerical)" => Tn, "Temperature (analytical)" => Ta) - write_file(outputpath * "result_steady_heat_equation.pvd", mesh, dict_vars) - - # 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(Tn_dofs .- Ta_dofs, Inf) / norm(Ta_dofs, Inf) < 2e-14 #src - end #src -end + +# 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) + +# Build function space and associated Trial and Test FE spaces. +# We impose a Dirichlet condition with a temperature of 260K +# on boundary "West" +fs = FunctionSpace(:Lagrange, degree) +U = TrialFESpace(fs, mesh, Dict("West" => 260.0)) +V = TestFESpace(U) + +# Define measures for cell integration +dΩ = Measure(CellDomain(mesh), 2 * degree + 1) + +# 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. +# 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) = ldiv!(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 (numerical)" => Tn, "Temperature (analytical)" => Ta) +write_file(outputpath * "result_steady_heat_equation.pvd", mesh, dict_vars) + +# 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(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. -function solve_unsteady() - println(" ----- Solving unsteady problem -----") - # Start by defining two additional parameters: - totalTime = 100.0 - Δt = 0.1 - - # Read a slightly different mesh - mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri_2.msh") - mesh = read_mesh(mesh_path) - - # The rest is similar to the steady case - fs = FunctionSpace(:Lagrange, degree) - 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Ω - m(u, v) = ∫(ρCp * u ⋅ v)dΩ - l(v) = ∫(q * v)dΩ - - A = assemble_bilinear(a, U, V) - M = assemble_bilinear(m, U, V) - L = assemble_linear(l, V) - - # Compute a vector of dofs whose values are zeros everywhere - # except on dofs lying on a Dirichlet boundary, where they - # take the Dirichlet value - Ud = assemble_dirichlet_vector(U, V, mesh) - - # Apply lift - L = L - A * Ud - - # Apply homogeneous dirichlet condition - apply_homogeneous_dirichlet_to_vector!(L, U, V, mesh) - apply_dirichlet_to_matrix!((A, M), U, V, mesh) - - # Form time iteration matrix - # (note that this is bad for performance since up to now, - # M and A are sparse matrices) - Miter = factorize(M + Δt * A) - - # Init the solution with a constant temperature of 260K - ϕ = FEFunction(U, 260.0) - - # Write initial solution to a file - mkpath(outputpath) - dict_vars = Dict("Temperature" => ϕ) - write_file(outputpath * "result_unsteady_heat_equation.pvd", mesh, dict_vars, 0, 0.0) - - # Time loop - itime = 0 - t = 0.0 - while t <= totalTime - t, itime - t += Δt - itime = itime + 1 + +println(" ----- Solving unsteady problem -----") +# Start by defining two additional parameters: +totalTime = 100.0 +Δt = 0.1 + +# Read a slightly different mesh +mesh_path = joinpath(@__DIR__, "..", "..", "input", "mesh", "domainSquare_tri_2.msh") +mesh = read_mesh(mesh_path) + +# The rest is similar to the steady case +fs = FunctionSpace(:Lagrange, degree) +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. +# 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Ω + +A = assemble_bilinear(a, U, V) +M = assemble_bilinear(m, U, V) +L = assemble_linear(l, V) + +# Compute a vector of dofs whose values are zeros everywhere +# except on dofs lying on a Dirichlet boundary, where they +# take the Dirichlet value +Ud = assemble_dirichlet_vector(U, V, mesh) + +# Apply lift +L = L - A * Ud + +# Apply homogeneous dirichlet condition +apply_homogeneous_dirichlet_to_vector!(L, U, V, mesh) +apply_dirichlet_to_matrix!((A, M), U, V, mesh) + +# Form time iteration matrix +# (note that this is bad for performance since up to now, +# M and A are sparse matrices) +Miter = factorize(M + Δt * A) + +# Init the solution with a constant temperature of 260K +ϕ = FEFunction(U, 260.0) + +# Write initial solution to a file +mkpath(outputpath) +dict_vars = Dict("Temperature" => ϕ) +write_file(outputpath * "result_unsteady_heat_equation.pvd", mesh, dict_vars, 0, 0.0) + +# Time loop +itime = 0 +t = 0.0 +while t <= totalTime + global t, itime + t += Δt + itime = itime + 1 #! format: off if !is_tested #src @show t, itime @@ -153,31 +155,27 @@ function solve_unsteady() #! format: on ## 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" => ϕ) - write_file( - outputpath * "result_unsteady_heat_equation.pvd", - mesh, - dict_vars, - itime, - t; - collection_append = true, - ) - end + 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" => ϕ) + write_file( + outputpath * "result_unsteady_heat_equation.pvd", + mesh, + dict_vars, + itime, + t; + collection_append = true, + ) end - - if is_tested #src - test_ref("heat_equation_100s.jld2", get_dof_values(ϕ)) #src - end #src end -solve_steady() -solve_unsteady() +if is_tested #src + test_ref("heat_equation_100s.jld2", get_dof_values(ϕ)) #src +end #src end #hide From eaf5e504db1b3af79c1131ab45c742dc9882b929 Mon Sep 17 00:00:00 2001 From: lokmanb Date: Thu, 4 Jun 2026 21:22:49 +0200 Subject: [PATCH 8/8] using backslash instead of ldiv in heat equation tutorial --- src/heat_equation/heat_equation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index de984467..90735683 100644 --- a/src/heat_equation/heat_equation.jl +++ b/src/heat_equation/heat_equation.jl @@ -64,7 +64,7 @@ l(v) = ∫(q * v)dΩ # 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) = ldiv!(y, cholesky(Symmetric(A)), x) +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)