diff --git a/src/heat_equation/heat_equation.jl b/src/heat_equation/heat_equation.jl index 76143a68..90735683 100644 --- a/src/heat_equation/heat_equation.jl +++ b/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,13 +148,13 @@ 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 @@ -158,7 +162,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, 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