module Covo #hide
+Euler equations - covo · BcubeTutorials Euler equations - covo
module Covo #hide
println("Running covo example...") #hide
const dir = string(@__DIR__, "/")
@@ -423,4 +423,4 @@
run_covo()
end
-end #hide
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
module EulerNacaSteady #hide
+Euler equations on a NACA0012 · BcubeTutorials Euler equations on a NACA0012
module EulerNacaSteady #hide
println("Running euler_naca_steady example...") #hide
# # Solve Euler equation around a NACA0012 airfoil
@@ -724,4 +724,4 @@
main(stateInit, stateBcFarfield, degreemax)
-end #hide
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
module linear_elasticity #hide
+Linear elasticity · BcubeTutorials Linear elasticity
module linear_elasticity #hide
println("Running linear elasticity API example...") #hide
# # Linear elasticity
@@ -186,4 +186,4 @@
#run_steady()
run_unsteady()
-end #hide
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
module linear_thermo_elasticity #hide
+Linear thermo-elasticity · BcubeTutorials Linear thermo-elasticity
module linear_thermo_elasticity #hide
println("Running linear thermo-elasticity API example...") #hide
# # Thermo-elasticity
@@ -222,4 +222,4 @@
# <img src="../assets/thermo_elasticity.gif" alt="drawing" width="500"/>
# ```
-end #hide
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
To be completed to answer common user questions.
mainLet's say you want to compare the performance of your current branch (named "target" hereafter) with the main branch (named "baseline" hereafter).
Open from Bcube.jl/ a REPL and type:
pkg> activate --temp
+How to... · BcubeTutorials How to
To be completed to answer common user questions.
Comparing manually the benchmarks with main
Let's say you want to compare the performance of your current branch (named "target" hereafter) with the main branch (named "baseline" hereafter).
Open from Bcube.jl/ a REPL and type:
pkg> activate --temp
pkg> add PkgBenchmark BenchmarkTools
pkg> dev .
using PkgBenchmark
@@ -34,4 +34,4 @@
using PkgBenchmark
import Bcube
results = benchmarkpkg(Bcube, BenchmarkConfig(; env = Dict("JULIA_NUM_THREADS" => "1")); resultfile = joinpath(@__DIR__, "result.json"))
-export_markdown("results.md", results)
This will create the markdown file results.md with the results.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
This will create the markdown file results.md with the results.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
Bcube is a Julia library providing tools for the spatial discretization of partial differential equation(s) (PDE). The main objectives are:
This documentation is organised as follow. Checkout the tutorials to see what Bcube is capable of and/or quickly learn how to use it. Then, some more elaborated examples are provided to demonstrate the library capabilities. The "Manual" part explains how the core is organized. Finally, the "API" section is the low level code documentation.
To write documentation for Bcube, Julia's guidelines should be followed : https://docs.julialang.org/en/v1/manual/documentation/. Moreover, this project tries to apply the SciML Style Guide.
This documentation follows the following notation or naming conventions:
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
Bcube is a Julia library providing tools for the spatial discretization of partial differential equation(s) (PDE). The main objectives are:
This documentation is organised as follow. Checkout the tutorials to see what Bcube is capable of and/or quickly learn how to use it. Then, some more elaborated examples are provided to demonstrate the library capabilities. The "Manual" part explains how the core is organized. Finally, the "API" section is the low level code documentation.
To write documentation for Bcube, Julia's guidelines should be followed : https://docs.julialang.org/en/v1/manual/documentation/. Moreover, this project tries to apply the SciML Style Guide.
This documentation follows the following notation or naming conventions:
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
As explained earlier, at least two coordinates systems exist in Bcube : the "reference" coordinates (ReferenceDomain) and the "physical" coordinates (PhysicalDomain). The evaluation of a function on a point in a cell depends on the way this point has been defined. Hence the definition of CellPoints that embed the coordinate system. Given a CellPoint (or eventually a FacePoint), an AbstractCellFunction will be evaluated and the mapping between the ReferenceDomain to the PhysicalDomain (or reciprocally) will be performed internally if necessary : if an AbstractCellFunction defined in terms of reference coordinates is applied on a CellPoint expressed in the reference coordinates system, no mapping is needed.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
AbstractFunctionSpaceIn Bcube, a FunctionSpace is defined by a type (nodal Lagrange polynomials, modal Taylor expansion, etc) and a degree. For each implemented FunctionSpace, a list of shape functions is associated on a given Shape. For instance, one can get the shape functions associated to the Lagrange polynomials or order 3 on a Square. Note that for "tensor" elements such as Line, Square or Cube; the Lagrange polynomials are available at any order; being computed symbolically.
AbstractFESpaceThen, an FESpace (more precisely SingleFESpace) is a function space associated to a numbering of the degrees of freedom. Note that the numbering may depend on the continuous or discontinuous feature of the space. Hence a SingleFESpace takes basically four input to be built : a FunctionSpace, the number of components of this space (scalar or vector), an indicator of the continuous/discontinuous characteristic, and the mesh. The dof numbering is built by combining the mesh numberings (nodes, cells, faces) and the function space. Note that the degree of the FunctionSpace can differ from the "degree" of the mesh elements : it is possible to build a SingleFESpace with P2 polynomials on a mesh only containing straight lines (defined by only two nodes, Bar2_t). Optionaly, a SingleFESpace can also contain the tags of the boundaries where Dirichlet condition(s) applies. A MultiFESpace is simply a set of SingleFESpace, eventually of different natures. Its befenit is that it allows to build a "global" numbering of all the dofs represented by this space. This is especially convenient to solve systems of equations.
AbstractFEFunctionWith a SingleFESpace, one can build the representation of a function discretized on this space: a FEFunction. This structure stores a vector of values, one for each degree of freedom of the finite element space. To set or get the values of a FEFunction, the functions set_dof_values! and get_dof_values are available respectively. A FEFunction can be projected on another FESpace; or evaluated at some specific mesh location (a coordinates, all the nodes, all the mesh centers, etc).
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
A Mesh is a set basically of nodes (Node), a set of entities (the mesh elements) and a list of connectivies that link the entities between themselves and with the nodes.
In Bcube every mesh entity has corresponding reference Shape, a simplified or canonical representation of this element. A 1D line is mapped on the [-1,1] segment, and a rectangle is mapped on a square for instance. On these reference shapes, (almost) everything is known : the vertices location, the area, the quadrature points... Hence in Bcube we always compute things in the reference shape. For "Lagrange" elements (such as Bar*_t, Tri*_t, Quad*_t, Tetra*_t, Hexa*_t, Penta*_t etc), the mapping from the reference shape to a geometrical element is directly obtained from the corresponding Lagrange polynomials and the element node coordinates. Given a geometrical element with n nodes M_i, the mapping reads:
\[F(\xi) = \sum_{i=1}^n \hat{\lambda}_i(\xi)M_i\]
where $(\lambda)_i$ are the Lagrange polynomials whose order matches the element order.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
To compute an integral on a geometrical element, for instance a curved element, a variable substitution is used to compute the integral on the corresponding reference Shape. This variable substitution reads:
\[\int_\Omega g(x) \mathrm{\,d} \Omega = \int_{\hat{\Omega}} |J(x)| \left(g \circ F \right)(\hat{x}) \mathrm{\,d} \hat{\Omega},\]
where we recall that $F$ is the reference to physical mapping and $J$ is the determinant of the jacobian matrix of this mapping. Depending on the shape and element order, this determinant is either hard-coded or computed with ForwardDiff.
Now, to compute the right side, i.e the integral on the reference shape, quadrature rules are applied:
\[\int_{\hat{\Omega}} g(\hat{x}) \mathrm{\,d} \hat{\Omega} = \sum_{i =1}^{N_q} \omega_i g(\hat{x}_i)\]
A specific procedure is applied to compute integrals on a face of a cell (i.e a surfacic integral on a face of a volumic element).
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
WORK IN PROGRESS
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
![\"drawing\"](\"../assets/helmholtz_x21_y21_vp6.png\")
","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"tutorial/phase_field_supercooled/#Phase-field-model-solidification-of-a-liquid-in-supercooled-state","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"","category":"section"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"In this tutorial, a coupled system of two unsteady equations is solved using finite elements and an imex time scheme. This tutorial doesn't introduce MultiFESpace, check the \"euler\" example for this. Warning : this file is currently quite long to run (a few minutes).","category":"page"},{"location":"tutorial/phase_field_supercooled/#Theory","page":"Phase field model - solidification of a liquid in supercooled state","title":"Theory","text":"","category":"section"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"This case is taken from: Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3-4), 410-423. In particular, the variables of the problem are denoted in the same way (p for the phase indicator and T for temperature). Consider a rectangular domain Omega = 0 L_x times 0 L_y on which we wish to solve the following equations:","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":" tau partial_t p = epsilon^2 Delta p + p (1-p)(p - frac12 + m(T))","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":" partial_t T = Delta T + K partial_t p","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"where m(T) = fracalphapi atan left gamma (T_e - T) right. This set of equations represents the solidification of a liquid in a supercooled state. Here T is a dimensionless temperature and p is the solid volume fraction. Lagrange finite elements are used to discretize both equations. Time marching is performed with a forward Euler scheme for the first equation and a backward Euler scheme for the second one.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"To initiate the solidification process, a Dirichlet boundary condition (p=1,T=1) is applied at x=0 (\"West\" boundary).","category":"page"},{"location":"tutorial/phase_field_supercooled/#Code","page":"Phase field model - solidification of a liquid in supercooled state","title":"Code","text":"","category":"section"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Load the necessary packages","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"using Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing Random\n\nRandom.seed!(1234) # to obtain reproductible results","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Define some physical and numerical constants, as well as the g function appearing in the problem definition.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"const dir = string(@__DIR__, \"/../\") # Bcube dir\nconst ε = 0.01\nconst τ = 0.0003\nconst α = 0.9\nconst γ = 10.0\nconst K = 1.6\nconst Te = 1.0\nconst β = 0.0 # noise amplitude, original value : 0.01\nconst Δt = 0.0001 # time step\nconst totalTime = 1.0 # original value : 1\nconst nout = 50 # Number of iterations to skip before writing file\nconst degree = 1 # function space degree\nconst lx = 3.0\nconst ly = 1.0\nconst nx = 100\nconst ny = 20\n\ng(T) = (α / π) * atan(γ * (Te - T))","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Read the mesh using gmsh","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"const mesh_path = dir * \"input/mesh/domainPhaseField_tri.msh\"\nconst mesh = read_msh(mesh_path)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Noise function : random between [-1/2,1/2]","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"const χ = MeshCellData(rand(ncells(mesh)) .- 0.5)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Build the function space and the FE Spaces. The two unknowns will share the same FE spaces for this tutorial. Note the way we specify the Dirichlet condition in the definition of U.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"fs = FunctionSpace(:Lagrange, degree)\nU = TrialFESpace(fs, mesh, Dict(\"West\" => (x, t) -> 1.0))\nV = TestFESpace(U)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Build FE functions","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"ϕ = FEFunction(U)\nT = FEFunction(U)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Define measures for cell integration","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"dΩ = Measure(CellDomain(mesh), 2 * degree + 1)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Define bilinear and linear forms","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"a(u, v) = ∫(∇(u) ⋅ ∇(v))dΩ\nm(u, v) = ∫(u ⋅ v)dΩ\nl(v) = ∫(v * ϕ * (1.0 - ϕ) * (ϕ - 0.5 + g(T) + β * χ))dΩ","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Assemble the two constant matrices","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"A = assemble_bilinear(a, U, V)\nM = assemble_bilinear(m, U, V)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Create iterative matrices","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"C_ϕ = M + Δt / τ * ε^2 * A\nC_T = M + Δt * A","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Apply Dirichlet conditions. For this example, we don't use a lifting method to impose the Dirichlet, but d is used to initialize the solution.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"d = assemble_dirichlet_vector(U, V, mesh)\napply_dirichlet_to_matrix!((C_ϕ, C_T), U, V, mesh)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Init solution and write it to a VTK file","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"set_dof_values!(ϕ, d)\nset_dof_values!(T, d)\n\ndict_vars = Dict(\n \"Temperature\" => (var_on_vertices(T, mesh), VTKPointData()),\n \"Phi\" => (var_on_vertices(ϕ, mesh), VTKPointData()),\n)\nwrite_vtk(dir * \"myout/result_phaseField_imex_1space\", 0, 0.0, mesh, dict_vars)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Factorize and allocate some vectors to increase performance","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"C_ϕ = factorize(C_ϕ)\nC_T = factorize(C_T)\nL = zero(d)\nrhs = zero(d)\nϕ_new = zero(d)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Time loop (imex time integration)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"t = 0.0\nitime = 0\nwhile t <= totalTime\n global t, itime\n t += Δt\n itime += 1\n @show t, totalTime\n\n # Integrate equation on ϕ\n L .= 0.0 # reset L\n assemble_linear!(L, l, V)\n rhs .= M * get_dof_values(ϕ) .+ Δt / τ .* L\n apply_dirichlet_to_vector!(rhs, U, V, mesh)\n ϕ_new .= C_ϕ \\ rhs\n\n # Integrate equation on T\n rhs .= M * (get_dof_values(T) .+ K .* (ϕ_new .- get_dof_values(ϕ)))\n apply_dirichlet_to_vector!(rhs, U, V, mesh)\n\n # Update solution\n set_dof_values!(ϕ, ϕ_new)\n set_dof_values!(T, C_T \\ rhs)\n\n # write solution in vtk format\n if itime % nout == 0\n dict_vars = Dict(\n \"Temperature\" => (var_on_vertices(T, mesh), VTKPointData()),\n \"Phi\" => (var_on_vertices(ϕ, mesh), VTKPointData()),\n )\n write_vtk(\n dir * \"myout/result_phaseField_imex_1space\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n end\nend","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"And here is an animation of the result:","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"![\"drawing\"](\"../assets/phase-field-supercooled-rectangle.gif\")
","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"This page was generated using Literate.jl.","category":"page"},{"location":"manual/cellfunction/#Cell-function","page":"Cell function","title":"Cell function","text":"","category":"section"},{"location":"manual/cellfunction/","page":"Cell function","title":"Cell function","text":"As explained earlier, at least two coordinates systems exist in Bcube : the \"reference\" coordinates (ReferenceDomain) and the \"physical\" coordinates (PhysicalDomain). The evaluation of a function on a point in a cell depends on the way this point has been defined. Hence the definition of CellPoints that embed the coordinate system. Given a CellPoint (or eventually a FacePoint), an AbstractCellFunction will be evaluated and the mapping between the ReferenceDomain to the PhysicalDomain (or reciprocally) will be performed internally if necessary : if an AbstractCellFunction defined in terms of reference coordinates is applied on a CellPoint expressed in the reference coordinates system, no mapping is needed.","category":"page"},{"location":"howto/howto/#How-to","page":"How to...","title":"How to","text":"","category":"section"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"To be completed to answer common user questions.","category":"page"},{"location":"howto/howto/#Comparing-manually-the-benchmarks-with-main","page":"How to...","title":"Comparing manually the benchmarks with main","text":"","category":"section"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Let's say you want to compare the performance of your current branch (named \"target\" hereafter) with the main branch (named \"baseline\" hereafter).","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Open from Bcube.jl/ a REPL and type:","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"pkg> activate --temp\npkg> add PkgBenchmark BenchmarkTools\npkg> dev .\nusing PkgBenchmark\nimport Bcube\nbenchmarkpkg(Bcube, BenchmarkConfig(; env = Dict(\"JULIA_NUM_THREADS\" => \"1\")); resultfile = joinpath(@__DIR__, \"result-target.json\"))","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create a result-target.json in the current directory.","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Then checkout the main branch. Start a fresh REPL and type (almost the same):","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"pkg> activate --temp\npkg> add PkgBenchmark BenchmarkTools\npkg> dev .\nusing PkgBenchmark\nimport Bcube\nbenchmarkpkg(Bcube, BenchmarkConfig(; env = Dict(\"JULIA_NUM_THREADS\" => \"1\")); resultfile = joinpath(@__DIR__, \"result-baseline.json\"))","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create a result-baseline.json in the current directory.","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"You can now \"compare\" the two files by running (watch-out for the order):","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"target = PkgBenchmark.readresults(\"result-target.json\")\nbaseline = PkgBenchmark.readresults(\"result-baseline.json\")\njudgement = judge(target, baseline)\nexport_markdown(\"judgement.md\", judgement)","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create the markdown file judgement.md with the results.","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"For more details, once you've built the judgement object, you can also type the following code from https://github.com/tkf/BenchmarkCI.jl:","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"open(\"detailed-judgement.md\", \"w\") do io\n println(io, \"# Judge result\")\n export_markdown(io, judgement)\n println(io)\n println(io)\n println(io, \"---\")\n println(io, \"# Target result\")\n export_markdown(io, PkgBenchmark.target_result(judgement))\n println(io)\n println(io)\n println(io, \"---\")\n println(io, \"# Baseline result\")\n export_markdown(io, PkgBenchmark.baseline_result(judgement))\n println(io)\n println(io)\n println(io, \"---\")\nend","category":"page"},{"location":"howto/howto/#Run-the-benchmark-manually","page":"How to...","title":"Run the benchmark manually","text":"","category":"section"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Let's say you want to run the benchmarks locally (without comparing with main)","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Open from Bcube.jl/ a REPL and type:","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"pkg> activate --temp\npkg> add PkgBenchmark\npkg> dev .\nusing PkgBenchmark\nimport Bcube\nresults = benchmarkpkg(Bcube, BenchmarkConfig(; env = Dict(\"JULIA_NUM_THREADS\" => \"1\")); resultfile = joinpath(@__DIR__, \"result.json\"))\nexport_markdown(\"results.md\", results)","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create the markdown file results.md with the results.","category":"page"},{"location":"example/linear_thermoelasticity/#Linear-thermo-elasticity","page":"Linear thermo-elasticity","title":"Linear thermo-elasticity","text":"","category":"section"},{"location":"example/linear_thermoelasticity/","page":"Linear thermo-elasticity","title":"Linear thermo-elasticity","text":"module linear_thermo_elasticity #hide\nprintln(\"Running linear thermo-elasticity API example...\") #hide\n\n# # Thermo-elasticity\n\nconst dir = string(@__DIR__, \"/\") # Bcube dir\nusing Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing StaticArrays\n\n# function space (here we shall use Lagrange P1 elements) and quadrature degree.\nconst fspace = :Lagrange\nconst degree = 1 # FunctionSpace degree\nconst degquad = 2 * degree + 1\n\n# Input and output paths\nconst outputpath = joinpath(dir, \"../myout/elasticity/\")\nconst meshpath = joinpath(dir, \"../input/mesh/domainThermoElast_tri.msh\")\n\n# Time stepping scheme params\nconst α = 0.05\nconst γ = 0.5 + α\nconst β = 0.25 * (1.0 + α)^2\n\nconst totalTime = 10.0\nconst Δt = 1.0e-2\n\n# Material parameters (Young's modulus, Poisson coefficient and deduced Lamé coefficients)\nconst E = 200.0e9\nconst ν = 0.3\nconst λ = E * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))\nconst μ = E / (2.0 * (1.0 + ν))\nconst Kₜ = 1.0e-6\nconst ρ = 2500.0\nconst cₚ = 1000.0\nconst k = 250.0\nconst T₀ = 280.0\n\n# Strain tensor and stress tensor (Hooke's law)\nϵ(u) = 0.5 * (∇(u) + transpose(∇(u)))\nσ(u) = λ * tr(ϵ(u)) * I + 2 * μ * (ϵ(u)) # Elastic stress\nσₜ(T) = (3 * λ + 2 * μ) * Kₜ * (T - T₀) * I # Thermal stress\n\nπ(u, v) = σ(u) ⊡ ϵ(v) # with the chosen contraction convention ϵ should be transposed, but as it is symmetric the expression remains correct\nπₜ(T, v) = σₜ(T) ⊡ ϵ(v)\n\n# materialize for identity operator\nBcube.materialize(A::LinearAlgebra.UniformScaling, B) = A\n\n# Function that performs a time step using a Newmark α-HHT scheme\n# The scheme updates the acceleration G, the velocity V and the displacement U using the following formulas:\n# ```math\n# \\begin{cases}\n# M G^{n+1} +(1-\\alpha)A U^{n+1} + \\alpha A U^{n} = (1-\\alpha) L^{n+1} + \\alpha L^n = L \\textrm{(because here $L$ is time independent)} \\\\\n# V^{n+1} = V^{n} + (1-\\gamma) \\Delta t G^n + \\gamma \\Delta t G^{n+1} \\\\\n# U^{n+1} = U^{n} + \\Delta t V^{n} + (\\frac{1}{2} - \\beta)*\\Delta t^2 G^{n} + \\beta \\Delta t^2 G^{n+1}\n# \\end{cases}\n# ```\n# where $$M$$ is the mass matrix, $$A$$ is the stiffness matrix and $$L$$ is the RHS\n# G is then computed by solving the linear system obtained by inserting the expressions for U and V in the equation for G.\nfunction Newmark_α_HHT(dt, L, A, Mat, U0, V0, G0)\n L1 = L - α * A * U0\n L2 = -(1.0 - α) * (A * U0 + dt * A * V0 + (0.5 - β) * dt * dt * A * G0)\n RHS = L1 .+ L2\n\n G = Mat \\ RHS\n U = U0 + dt * V0 + (0.5 - β) * dt * dt * G0 + β * dt * dt * G\n V = V0 + (1.0 - γ) * dt * G0 + γ * dt * G\n\n return U, V, G\nend\n\n# Function that runs the unsteady case:\nfunction run_unsteady()\n mesh = read_msh(meshpath, 2)\n\n fs = FunctionSpace(fspace, degree)\n U_scal = TrialFESpace(fs, mesh, Dict(\"West1\" => 280.0, \"East1\" => 280.0); size = 1)\n V_scal = TestFESpace(U_scal)\n U_vec = TrialFESpace(\n fs,\n mesh,\n Dict(\"West1\" => SA[0.0, 0.0], \"East1\" => SA[0.0, 0.0]);\n size = 2,\n )\n V_vec = TestFESpace(U_vec)\n\n # Initialize solution\n U = FEFunction(U_vec, 0.0)\n U0 = zeros(Bcube.get_ndofs(U_vec))\n V0 = zeros(Bcube.get_ndofs(U_vec))\n G0 = zeros(Bcube.get_ndofs(U_vec))\n\n T = FEFunction(U_scal, T₀)\n\n # Define measures for cell\n dΩ = Measure(CellDomain(mesh), degquad)\n\n # no volume force term\n f = PhysicalFunction(x -> SA[0.0, 0.0])\n\n q = PhysicalFunction(\n x -> x[1] .* (1.0 .- x[1]) .* x[2] .* (0.2 .- x[2]) .* 1500000000.0,\n )\n\n # Definition of bilinear and linear forms for the elasticity problem\n a(u, v) = ∫(π(u, v))dΩ\n m(u, v) = ∫(ρ * u ⋅ v)dΩ\n l(v) = ∫(πₜ(T, v))dΩ\n\n # An alternative way to define this linear form is to use operator composition:\n # l(v) = ∫( πₜ ∘ (T, v, ∇(v)) )dΩ\n # where πₜ(T, v, ∇v) = σₜ(T) ⊡ ϵ(v, ∇v) and ϵ(v, ∇v) = 0.5*( ∇v + transpose(∇v) )\n\n # Definition of bilinear and linear forms for the heat conduction problem\n aₜ(u, v) = ∫(k * ∇(u) ⋅ ∇(v))dΩ\n mₜ(u, v) = ∫(ρ * cₚ * u ⋅ v)dΩ\n lₜ(v) = ∫(q * v)dΩ\n\n # Assemble matrices and vector\n M = assemble_bilinear(m, U_vec, V_vec)\n A = assemble_bilinear(a, U_vec, V_vec)\n L = assemble_linear(l, V_vec)\n AT = assemble_bilinear(aₜ, U_scal, V_scal)\n MT = assemble_bilinear(mₜ, U_scal, V_scal)\n LT = assemble_linear(lₜ, V_scal)\n\n # Apply homogeneous dirichlet on A and b\n Bcube.apply_homogeneous_dirichlet_to_vector!(L, U_vec, V_vec, mesh)\n Bcube.apply_dirichlet_to_matrix!((A, M), U_vec, V_vec, mesh)\n\n # Compute a vector of dofs whose values are zeros everywhere\n # except on dofs lying on a Dirichlet boundary, where they\n # take the Dirichlet value\n Td = Bcube.assemble_dirichlet_vector(U_scal, V_scal, mesh)\n\n # Apply lift\n LT = LT - AT * Td\n\n # Apply homogeneous dirichlet condition\n Bcube.apply_homogeneous_dirichlet_to_vector!(LT, U_scal, V_scal, mesh)\n Bcube.apply_dirichlet_to_matrix!((AT, MT), U_scal, V_scal, mesh)\n\n # Write initial solution\n Un = var_on_vertices(U, mesh)\n Un = transpose(Un)\n Tn = var_on_vertices(T, mesh)\n mkpath(outputpath)\n dict_vars =\n Dict(\"Displacement\" => (Un, VTKPointData()), \"Temperature\" => (Tn, VTKPointData()))\n # Write the obtained FE solution\n write_vtk(\n outputpath * \"result_thermoelasticity\",\n 0,\n 0.0,\n mesh,\n dict_vars;\n append = false,\n )\n\n # Time loop\n itime = 0\n t = 0.0\n\n # Matrix for time stepping\n Mat = factorize(M + (1.0 - α) * (β * Δt * Δt * A))\n Miter = factorize(MT + Δt * AT)\n\n while t <= totalTime\n t += Δt\n itime = itime + 1\n @show t, itime\n\n # solve time step heat equation\n rhs = Δt * LT + MT * (get_dof_values(T) .- Td)\n set_dof_values!(T, Miter \\ rhs .+ Td)\n\n # solve time step elasticity\n U1, V1, G1 = Newmark_α_HHT(Δt, L, A, Mat, U0, V0, G0)\n\n # Update solution\n U0 .= U1\n V0 .= V1\n G0 .= G1\n\n set_dof_values!(U, U1)\n L = assemble_linear(l, V_vec)\n Bcube.apply_homogeneous_dirichlet_to_vector!(L, U_vec, V_vec, mesh)\n\n # Write solution\n if itime % 10 == 0\n Un = var_on_vertices(U, mesh)\n Un = transpose(Un)\n Tn = var_on_vertices(T, mesh)\n mkpath(outputpath)\n dict_vars = Dict(\n \"Displacement\" => (Un, VTKPointData()),\n \"Temperature\" => (Tn, VTKPointData()),\n )\n # Write the obtained FE solution\n write_vtk(\n outputpath * \"result_thermoelasticity\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n # In order to use the warp function in paraview (solid is deformed using the displacement field)\n # the calculator filter has to be used with the following formula to reconstruct a 3D displacement field\n # with 0 z-component: Displacement_X*iHat+Displacement_Y*jHat+0.0*kHat\n end\n end\nend\n\nrun_unsteady()\n\n# Here is an animation of the obtained result:\n# ```@raw html\n# ![\"drawing\"](\"../assets/thermo_elasticity.gif\")
\n# ```\n\nend #hide","category":"page"},{"location":"manual/function_space/#Function-and-FE-spaces","page":"Function and FE spaces","title":"Function and FE spaces","text":"","category":"section"},{"location":"manual/function_space/#AbstractFunctionSpace","page":"Function and FE spaces","title":"AbstractFunctionSpace","text":"","category":"section"},{"location":"manual/function_space/","page":"Function and FE spaces","title":"Function and FE spaces","text":"In Bcube, a FunctionSpace is defined by a type (nodal Lagrange polynomials, modal Taylor expansion, etc) and a degree. For each implemented FunctionSpace, a list of shape functions is associated on a given Shape. For instance, one can get the shape functions associated to the Lagrange polynomials or order 3 on a Square. Note that for \"tensor\" elements such as Line, Square or Cube; the Lagrange polynomials are available at any order; being computed symbolically.","category":"page"},{"location":"manual/function_space/#AbstractFESpace","page":"Function and FE spaces","title":"AbstractFESpace","text":"","category":"section"},{"location":"manual/function_space/","page":"Function and FE spaces","title":"Function and FE spaces","text":"Then, an FESpace (more precisely SingleFESpace) is a function space associated to a numbering of the degrees of freedom. Note that the numbering may depend on the continuous or discontinuous feature of the space. Hence a SingleFESpace takes basically four input to be built : a FunctionSpace, the number of components of this space (scalar or vector), an indicator of the continuous/discontinuous characteristic, and the mesh. The dof numbering is built by combining the mesh numberings (nodes, cells, faces) and the function space. Note that the degree of the FunctionSpace can differ from the \"degree\" of the mesh elements : it is possible to build a SingleFESpace with P2 polynomials on a mesh only containing straight lines (defined by only two nodes, Bar2_t). Optionaly, a SingleFESpace can also contain the tags of the boundaries where Dirichlet condition(s) applies. A MultiFESpace is simply a set of SingleFESpace, eventually of different natures. Its befenit is that it allows to build a \"global\" numbering of all the dofs represented by this space. This is especially convenient to solve systems of equations.","category":"page"},{"location":"manual/function_space/#AbstractFEFunction","page":"Function and FE spaces","title":"AbstractFEFunction","text":"","category":"section"},{"location":"manual/function_space/","page":"Function and FE spaces","title":"Function and FE spaces","text":"With a SingleFESpace, one can build the representation of a function discretized on this space: a FEFunction. This structure stores a vector of values, one for each degree of freedom of the finite element space. To set or get the values of a FEFunction, the functions set_dof_values! and get_dof_values are available respectively. A FEFunction can be projected on another FESpace; or evaluated at some specific mesh location (a coordinates, all the nodes, all the mesh centers, etc).","category":"page"},{"location":"tutorial/heat_equation/#Heat-equation-(FE)","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"In this tutorial, the heat equation (first steady and then unsteady) is solved using finite-elements.","category":"page"},{"location":"tutorial/heat_equation/#Theory","page":"Heat equation (FE)","title":"Theory","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"This example shows how to solve the heat equation with eventually variable physical properties in steady and unsteady formulations:","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":" rho C_p partial_t u - nabla ( lambda u) = f","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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:","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":" forall v in tildeH^1_0(Omega) underbraceint_Omega partial_t u v dx_m(partial_t uv) + underbraceint_Omega nabla u nabla v dx_a(uv) = underbraceint_Omega f v dx_l(v)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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.","category":"page"},{"location":"tutorial/heat_equation/#Steady-case","page":"Heat equation (FE)","title":"Steady case","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"As usual, start by importing the necessary packages.","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"using Bcube\nusing LinearAlgebra\nusing WriteVTK","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"First we define some physical and numerical constants","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"const htc = 100.0 # Heat transfer coefficient (bnd cdt)\nconst Tr = 268.0 # Recovery temperature (bnd cdt)\nconst phi = 100.0\nconst q = 1500.0\nconst λ = 100.0\nconst η = λ\nconst ρCp = 100.0 * 200.0\nconst degree = 2\nconst outputpath = joinpath(@__DIR__, \"../myout/heat_equation/\")","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Read 2D mesh","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mesh_path = joinpath(@__DIR__, \"../input/mesh/domainSquare_tri.msh\")\nmesh = read_msh(mesh_path)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Build function space and associated Trial and Test FE spaces. We impose a Dirichlet condition with a temperature of 260K on boundary \"West\"","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"fs = FunctionSpace(:Lagrange, degree)\nU = TrialFESpace(fs, mesh, Dict(\"West\" => 260.0))\nV = TestFESpace(U)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Define measures for cell integration","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"dΩ = Measure(CellDomain(mesh), 2 * degree + 1)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Define bilinear and linear forms","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"a(u, v) = ∫(η * ∇(u) ⋅ ∇(v))dΩ\nl(v) = ∫(q * v)dΩ","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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.","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"sys = AffineFESystem(a, l, U, V)\nϕ = solve(sys)\nTcn = var_on_centers(ϕ, mesh)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute analytical solution for comparison. Apply the analytical solution on mesh centers","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"T_analytical = x -> 260.0 + (q / λ) * x[1] * (1.0 - 0.5 * x[1])\nTca = map(T_analytical, get_cell_centers(mesh))","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Write both the obtained FE solution and the analytical solution to a vtk file.","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mkpath(outputpath)\ndict_vars =\n Dict(\"Temperature\" => (Tcn, VTKCellData()), \"Temperature_a\" => (Tca, VTKCellData()))\nwrite_vtk(outputpath * \"result_steady_heat_equation\", 0, 0.0, mesh, dict_vars)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute and display the error","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"@show norm(Tcn .- Tca, Inf) / norm(Tca, Inf)","category":"page"},{"location":"tutorial/heat_equation/#Unsteady-case","page":"Heat equation (FE)","title":"Unsteady case","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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:","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"totalTime = 100.0\nΔt = 0.1","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Read a slightly different mesh","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mesh_path = joinpath(@__DIR__, \"../input/mesh/domainSquare_tri_2.msh\")\nmesh = read_msh(mesh_path)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"The rest is similar to the steady case","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"fs = FunctionSpace(:Lagrange, degree)\nU = TrialFESpace(fs, mesh, Dict(\"West\" => 260.0))\nV = TestFESpace(U)\ndΩ = Measure(CellDomain(mesh), 2 * degree + 1)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute matrices associated to bilinear and linear forms, and assemble","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"a(u, v) = ∫(η * ∇(u) ⋅ ∇(v))dΩ\nm(u, v) = ∫(ρCp * u ⋅ v)dΩ\nl(v) = ∫(q * v)dΩ\n\nA = assemble_bilinear(a, U, V)\nM = assemble_bilinear(m, U, V)\nL = assemble_linear(l, V)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute a vector of dofs whose values are zeros everywhere except on dofs lying on a Dirichlet boundary, where they take the Dirichlet value","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Ud = assemble_dirichlet_vector(U, V, mesh)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Apply lift","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"L = L - A * Ud","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Apply homogeneous dirichlet condition","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"apply_homogeneous_dirichlet_to_vector!(L, U, V, mesh)\napply_dirichlet_to_matrix!((A, M), U, V, mesh)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Form time iteration matrix (note that this is bad for performance since up to now, M and A are sparse matrices)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Miter = factorize(M + Δt * A)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Init the solution with a constant temperature of 260K","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"ϕ = FEFunction(U, 260.0)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Write initial solution to a file","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mkpath(outputpath)\ndict_vars = Dict(\"Temperature\" => (var_on_centers(ϕ, mesh), VTKCellData()))\nwrite_vtk(outputpath * \"result_unsteady_heat_equation\", 0, 0.0, mesh, dict_vars)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Time loop","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"itime = 0\nt = 0.0\nwhile t <= totalTime\n global t, itime\n t += Δt\n itime = itime + 1\n @show t, itime\n\n # Compute rhs\n rhs = Δt * L + M * (get_dof_values(ϕ) .- Ud)\n\n # Invert system and apply inverse shift\n set_dof_values!(ϕ, Miter \\ rhs .+ Ud)\n\n # Write solution (every 10 iterations)\n if itime % 10 == 0\n dict_vars = Dict(\"Temperature\" => (var_on_centers(ϕ, mesh), VTKCellData()))\n write_vtk(\n outputpath * \"result_unsteady_heat_equation\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n end\nend\n","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"example/covo/#Euler-equations-covo","page":"Euler equations - covo","title":"Euler equations - covo","text":"","category":"section"},{"location":"example/covo/","page":"Euler equations - covo","title":"Euler equations - covo","text":"module Covo #hide\nprintln(\"Running covo example...\") #hide\n\nconst dir = string(@__DIR__, \"/\")\nusing Bcube\nusing LinearAlgebra\nusing StaticArrays\nusing WriteVTK # only for 'VTKCellData'\nusing Profile\nusing StaticArrays\nusing InteractiveUtils\nusing BenchmarkTools\nusing UnPack\n\nfunction compute_residual(_u, V, params, cache)\n u = get_fe_functions(_u)\n\n # alias on measures\n @unpack dΩ, dΓ, dΓ_perio_x, dΓ_perio_y = params\n\n # face normals for each face domain (lazy, no computation at this step)\n nΓ = get_face_normals(dΓ)\n nΓ_perio_x = get_face_normals(dΓ_perio_x)\n nΓ_perio_y = get_face_normals(dΓ_perio_y)\n\n # flux residuals from faces for all variables\n function l(v)\n ∫(flux_Ω(u, v))dΩ +\n -∫(flux_Γ(u, v, nΓ))dΓ +\n -∫(flux_Γ(u, v, nΓ_perio_x))dΓ_perio_x +\n -∫(flux_Γ(u, v, nΓ_perio_y))dΓ_perio_y\n end\n\n rhs = assemble_linear(l, V)\n\n return cache.mass \\ rhs\nend\n\n\"\"\"\n flux_Ω(u, v)\n\nCompute volume residual using the lazy-operators approach\n\"\"\"\nflux_Ω(u, v) = _flux_Ω ∘ cellvar(u, v)\ncellvar(u, v) = (u, map(∇, v))\nfunction _flux_Ω(u, ∇v)\n ρ, ρu, ρE, ϕ = u\n ∇λ_ρ, ∇λ_ρu, ∇λ_ρE, ∇λ_ϕ = ∇v\n\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = pressure(ρ, ρu, ρE, γ)\n\n flux_ρ = ρu\n flux_ρu = ρuu + p * I\n flux_ρE = (ρE + p) .* vel\n flux_ϕ = ϕ .* vel\n\n return ∇λ_ρ ⋅ flux_ρ + ∇λ_ρu ⊡ flux_ρu + ∇λ_ρE ⋅ flux_ρE + ∇λ_ϕ ⋅ flux_ϕ\nend\n\n\"\"\"\n flux_Γ(u,v,n)\n\nFlux at the interface is defined by a composition of two functions:\n* facevar(u,v,n) defines the input states which are needed for\n the riemann flux using operator notations\n* flux_roe(w) defines the Riemann flux (as usual)\n\"\"\"\nflux_Γ(u, v, n) = flux_roe ∘ (side⁻(u), side⁺(u), jump(v), side⁻(n))\n\n\"\"\"\n flux_roe(w)\n\"\"\"\nfunction flux_roe(ui, uj, δv, nij)\n # destructuring inputs given by `facevar` function\n\n nx, ny = nij\n ρ1, ρu1, ρE1, ϕ1 = ui\n ρ2, ρu2, ρE2, ϕ2 = uj\n δλ_ρ1, δλ_ρu1, δλ_ρE1, δλ_ϕ1 = δv\n ρux1, ρuy1 = ρu1\n ρux2, ρuy2 = ρu2\n\n # Closure\n u1 = ρux1 / ρ1\n v1 = ρuy1 / ρ1\n u2 = ρux2 / ρ2\n v2 = ρuy2 / ρ2\n p1 = pressure(ρ1, ρu1, ρE1, γ)\n p2 = pressure(ρ2, ρu2, ρE2, γ)\n\n H2 = (γ / (γ - 1)) * p2 / ρ2 + (u2 * u2 + v2 * v2) / 2.0\n H1 = (γ / (γ - 1)) * p1 / ρ1 + (u1 * u1 + v1 * v1) / 2.0\n\n R = √(ρ1 / ρ2)\n invR1 = 1.0 / (R + 1)\n uAv = (R * u1 + u2) * invR1\n vAv = (R * v1 + v2) * invR1\n Hav = (R * H1 + H2) * invR1\n cAv = √(abs((γ - 1) * (Hav - (uAv * uAv + vAv * vAv) / 2.0)))\n ecAv = (uAv * uAv + vAv * vAv) / 2.0\n\n λ1 = nx * uAv + ny * vAv\n λ3 = λ1 + cAv\n λ4 = λ1 - cAv\n\n d1 = ρ1 - ρ2\n d2 = ρ1 * u1 - ρ2 * u2\n d3 = ρ1 * v1 - ρ2 * v2\n d4 = ρE1 - ρE2\n\n # computation of the centered part of the flux\n flu11 = nx * ρ2 * u2 + ny * ρ2 * v2\n flu21 = nx * p2 + flu11 * u2\n flu31 = ny * p2 + flu11 * v2\n flu41 = H2 * flu11\n\n # Temp variables\n rc1 = (γ - 1) / cAv\n rc2 = (γ - 1) / cAv / cAv\n uq41 = ecAv / cAv + cAv / (γ - 1)\n uq42 = nx * uAv + ny * vAv\n\n fdc1 = max(λ1, 0.0) * (d1 + rc2 * (-ecAv * d1 + uAv * d2 + vAv * d3 - d4))\n fdc2 = max(λ1, 0.0) * ((nx * vAv - ny * uAv) * d1 + ny * d2 - nx * d3)\n fdc3 =\n max(λ3, 0.0) * (\n (-uq42 * d1 + nx * d2 + ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n fdc4 =\n max(λ4, 0.0) * (\n (uq42 * d1 - nx * d2 - ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n\n duv1 = fdc1 + (fdc3 + fdc4) / cAv\n duv2 = uAv * fdc1 + ny * fdc2 + (uAv / cAv + nx) * fdc3 + (uAv / cAv - nx) * fdc4\n duv3 = vAv * fdc1 - nx * fdc2 + (vAv / cAv + ny) * fdc3 + (vAv / cAv - ny) * fdc4\n duv4 =\n ecAv * fdc1 +\n (ny * uAv - nx * vAv) * fdc2 +\n (uq41 + uq42) * fdc3 +\n (uq41 - uq42) * fdc4\n\n v₁₂ = 0.5 * ((u1 + u2) * nx + (v1 + v2) * ny)\n fluxϕ = max(0.0, v₁₂) * ϕ1 + min(0.0, v₁₂) * ϕ2\n\n return (\n δλ_ρ1 * (flu11 + duv1) +\n δλ_ρu1 ⋅ (SA[flu21 + duv2, flu31 + duv3]) +\n δλ_ρE1 * (flu41 + duv4) +\n δλ_ϕ1 * (fluxϕ)\n )\nend\n\n\"\"\"\nTime integration of `f(q, t)` over a timestep `Δt`.\n\"\"\"\nforward_euler(q, f::Function, t, Δt) = get_dof_values(q) .+ Δt .* f(q, t)\n\n\"\"\"\n rk3_ssp(q, f::Function, t, Δt)\n\n`f(q, t)` is the function to integrate.\n\"\"\"\nfunction rk3_ssp(q, f::Function, t, Δt)\n stepper(q, t) = forward_euler(q, f, t, Δt)\n _q0 = get_dof_values(q)\n\n _q1 = stepper(q, Δt)\n\n set_dof_values!(q, _q1)\n _q2 = (3 / 4) .* _q0 .+ (1 / 4) .* stepper(q, t + Δt)\n\n set_dof_values!(q, _q2)\n _q1 .= (1 / 3) * _q0 .+ (2 / 3) .* stepper(q, t + Δt / 2)\n\n return _q1\nend\n\n\"\"\"\n pressure(ρ, ρu, ρE, γ)\n\nComputes pressure from perfect gaz law\n\"\"\"\nfunction pressure(ρ::Number, ρu::AbstractVector, ρE::Number, γ)\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = (γ - 1) * (ρE - tr(ρuu) / 2)\n return p\nend\n\n\"\"\"\n Init field with a vortex (for the COVO test case)\n\"\"\"\nfunction covo!(q, dΩ)\n\n # Intermediate vars\n Cₚ = γ * r / (γ - 1)\n\n r²(x) = ((x[1] .- xvc) .^ 2 + (x[2] .- yvc) .^ 2) ./ Rc^2\n # Temperature\n T(x) = T₀ .- β^2 * U₀^2 / (2 * Cₚ) .* exp.(-r²(x))\n # Velocity\n ux(x) = U₀ .- β * U₀ / Rc .* (x[2] .- yvc) .* exp.(-r²(x) ./ 2)\n uy(x) = V₀ .+ β * U₀ / Rc .* (x[1] .- xvc) .* exp.(-r²(x) ./ 2)\n # Density\n ρ(x) = ρ₀ .* (T(x) ./ T₀) .^ (1.0 / (γ - 1))\n # momentum\n ρu(x) = SA[ρ(x) * ux(x), ρ(x) * uy(x)]\n # Energy\n ρE(x) = ρ(x) * ((Cₚ / γ) .* T(x) + (ux(x) .^ 2 + uy(x) .^ 2) ./ 2)\n # Passive scalar\n ϕ(x) = Rc^2 * r²(x) < 0.01 ? exp(-r²(x) ./ 2) : 0.0\n\n f = map(PhysicalFunction, (ρ, ρu, ρE, ϕ))\n projection_l2!(q, f, dΩ)\n return nothing\nend\n\n\"\"\"\n Tiny struct to ease the VTK output\n\"\"\"\nmutable struct VtkHandler\n basename::Any\n ite::Any\n VtkHandler(basename) = new(basename, 0)\nend\n\n\"\"\"\n Write solution to vtk\n Wrapper for `write_vtk`\n\"\"\"\nfunction append_vtk(vtk, mesh, vars, t, params)\n ρ, ρu, ρE, ϕ = vars\n\n mesh_degree = 1\n vtk_degree = maximum(x -> get_degree(Bcube.get_function_space(get_fespace(x))), vars)\n vtk_degree = max(1, mesh_degree, vtk_degree)\n\n _ρ = var_on_nodes_discontinuous(ρ, mesh, vtk_degree)\n _ρu = var_on_nodes_discontinuous(ρu, mesh, vtk_degree)\n _ρE = var_on_nodes_discontinuous(ρE, mesh, vtk_degree)\n _ϕ = var_on_nodes_discontinuous(ϕ, mesh, vtk_degree)\n\n _p = pressure.(_ρ, _ρu, _ρE, γ)\n dict_vars_dg = Dict(\n \"rho\" => (_ρ, VTKPointData()),\n \"rhou\" => (_ρu, VTKPointData()),\n \"rhoE\" => (_ρE, VTKPointData()),\n \"phi\" => (_ϕ, VTKPointData()),\n \"p\" => (_p, VTKPointData()),\n )\n Bcube.write_vtk_discontinuous(\n vtk.basename * \"_DG\",\n vtk.ite,\n t,\n mesh,\n dict_vars_dg,\n vtk_degree;\n append = vtk.ite > 0,\n )\n\n # Update counter\n vtk.ite += 1\nend\n\n# Settings\nif get(ENV, \"BenchmarkMode\", \"false\") == \"false\" #hide\n const cellfactor = 1\n const nx = 32 * cellfactor + 1\n const ny = 32 * cellfactor + 1\n const fspace = :Lagrange\n const timeScheme = :ForwardEuler\nelse #hide\n const nx = 128 + 1 #hide\n const ny = 128 + 1 #hide\n const fspace = :Lagrange\n const timeScheme = :ForwardEuler\nend #hide\nconst nperiod = 1 # number of turn\nconst CFL = 0.1\nconst degree = 1 # FunctionSpace degree\nconst degquad = 2 * degree + 1\nconst γ = 1.4\nconst β = 0.2 # vortex intensity\nconst r = 287.15 # Perfect gaz constant\nconst T₀ = 300 # mean-flow temperature\nconst P₀ = 1e5 # mean-flow pressure\nconst M₀ = 0.5 # mean-flow mach number\nconst ρ₀ = 1.0 # mean-flow density\nconst xvc = 0.0 # x-center of vortex\nconst yvc = 0.0 # y-center of vortex\nconst Rc = 0.005 # Charasteristic vortex radius\nconst c₀ = √(γ * r * T₀) # Sound velocity\nconst U₀ = M₀ * c₀ # mean-flow velocity\nconst V₀ = 0.0 # mean-flow velocity\nconst ϕ₀ = 1.0\nconst l = 0.05 # half-width of the domain\nconst Δt = CFL * 2 * l / (nx - 1) / ((1 + β) * U₀ + c₀) / (2 * degree + 1)\n#const Δt = 5.e-7\nconst nout = 100 # Number of time steps to save\nconst outputpath = \"../myout/covo/\"\nconst output = joinpath(@__DIR__, outputpath, \"covo_deg$degree\")\nconst nite = Int(floor(nperiod * 2 * l / (U₀ * Δt))) + 1\n\nfunction run_covo()\n println(\"Starting run_covo...\")\n\n # Build mesh\n meshParam = (nx = nx, ny = ny, lx = 2l, ly = 2l, xc = 0.0, yc = 0.0)\n tmp_path = \"tmp.msh\"\n if get(ENV, \"BenchmarkMode\", \"false\") == \"false\" #hide\n gen_rectangle_mesh(tmp_path, :quad; meshParam...)\n else #hide\n if get(ENV, \"MeshConfig\", \"quad\") == \"triquad\" #hide\n gen_rectangle_mesh_with_tri_and_quad(tmp_path; meshParam...) #hide\n else #hide\n gen_rectangle_mesh(tmp_path, :quad; meshParam...) #hide\n end #hide\n end #hide\n mesh = read_msh(tmp_path)\n rm(tmp_path)\n\n # Define variables and test functions\n fs = FunctionSpace(fspace, degree)\n U_sca = TrialFESpace(fs, mesh, :discontinuous; size = 1) # DG, scalar\n U_vec = TrialFESpace(fs, mesh, :discontinuous; size = 2) # DG, vectoriel\n V_sca = TestFESpace(U_sca)\n V_vec = TestFESpace(U_vec)\n U = MultiFESpace(U_sca, U_vec, U_sca, U_sca)\n V = MultiFESpace(V_sca, V_vec, V_sca, V_sca)\n u = FEFunction(U)\n\n @show Bcube.get_ndofs(U)\n\n # Define measures for cell and interior face integrations\n dΩ = Measure(CellDomain(mesh), degquad)\n dΓ = Measure(InteriorFaceDomain(mesh), degquad)\n\n # Declare periodic boundary conditions and\n # create associated domains and measures\n periodicBCType_x = PeriodicBCType(Translation(SA[-2l, 0.0]), (\"East\",), (\"West\",))\n periodicBCType_y = PeriodicBCType(Translation(SA[0.0, 2l]), (\"South\",), (\"North\",))\n Γ_perio_x = BoundaryFaceDomain(mesh, periodicBCType_x)\n Γ_perio_y = BoundaryFaceDomain(mesh, periodicBCType_y)\n dΓ_perio_x = Measure(Γ_perio_x, degquad)\n dΓ_perio_y = Measure(Γ_perio_y, degquad)\n\n params = (dΩ = dΩ, dΓ = dΓ, dΓ_perio_x = dΓ_perio_x, dΓ_perio_y = dΓ_perio_y)\n\n # Init vtk\n isdir(joinpath(@__DIR__, outputpath)) || mkpath(joinpath(@__DIR__, outputpath))\n vtk = VtkHandler(output)\n\n # Init solution\n t = 0.0\n\n covo!(u, dΩ)\n\n # cache mass matrices\n cache = (mass = factorize(Bcube.build_mass_matrix(U, V, dΩ)),)\n\n if get(ENV, \"BenchmarkMode\", \"false\") == \"true\" #hide\n return u, U, V, params, cache\n end\n\n # Write initial solution\n append_vtk(vtk, mesh, u, t, params)\n\n # Time loop\n for i in 1:nite\n println(\"\")\n println(\"\")\n println(\"Iteration \", i, \" / \", nite)\n\n ## Step forward in time\n rhs(u, t) = compute_residual(u, V, params, cache)\n if timeScheme == :ForwardEuler\n unew = forward_euler(u, rhs, time, Δt)\n elseif timeScheme == :RK3\n unew = rk3_ssp(u, rhs, time, Δt)\n else\n error(\"Unknown time scheme: $timeScheme\")\n end\n\n set_dof_values!(u, unew)\n\n t += Δt\n\n # Write solution to file\n if (i % Int(max(floor(nite / nout), 1)) == 0)\n println(\"--> VTK export\")\n append_vtk(vtk, mesh, u, t, params)\n end\n end\n\n # Summary and benchmark # ndofs total = 20480\n _rhs(u, t) = compute_residual(u, V, params, cache)\n @btime forward_euler($u, $_rhs, $time, $Δt) # 5.639 ms (1574 allocations: 2.08 MiB)\n # stepper = w -> explicit_step(w, params, cache, Δt)\n # RK3_SSP(stepper, (u, v), cache)\n # @btime RK3_SSP($stepper, ($u, $v), $cache)\n println(\"ndofs total = \", Bcube.get_ndofs(U))\n Profile.init(; n = 10^7) # returns the current settings\n Profile.clear()\n Profile.clear_malloc_data()\n @profile begin\n for i in 1:100\n forward_euler(u, _rhs, time, Δt)\n end\n end\n @show Δt, U₀, U₀ * t\n @show boundary_names(mesh)\n return nothing\nend\n\nif get(ENV, \"BenchmarkMode\", \"false\") == \"false\"\n mkpath(outputpath)\n run_covo()\nend\n\nend #hide","category":"page"},{"location":"tutorial/linear_transport/#Linear-transport-(DG)","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"","category":"section"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"In this tutorial, we show how to solve a linear transport equation using a discontinuous-Galerkin framework with Bcube.","category":"page"},{"location":"tutorial/linear_transport/#Theory","page":"Linear transport (DG)","title":"Theory","text":"","category":"section"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"In this example, we solve the following linear transport equation using discontinuous elements:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"fracpartial phipartial t + nabla cdot (c phi) = 0","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"where c is a constant velocity. Using an explicit time scheme, one obtains:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"phi^n+1 = phi^n - Delta t nabla cdot (c phi^n)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"The corresponding weak form of this equation is:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"int_Omega phi^n+1 v mathrmdOmega = int_Omega phi^n v mathrmdOmega + Delta t left\nint_Omega c phi^n cdot nabla v mathrmdOmega - oint_Gamma left( c phi cdot n right) v mathrmdGamma\nright","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"where Gamma = delta Omega. Adopting the discontinuous Galerkin framework, this equation is written in every mesh cell Omega_i. The cell boundary term involves discontinuous quantities and is replaced by a \"numerical flux\", leading to the expression:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"int_Omega_i phi^n+1 v mathrmdOmega_i = int_Omega_i phi^n v mathrmdOmega_i + Delta t left\nint_Omega_i c phi^n cdot nabla v mathrmdOmega_i - oint_Gamma_i F^*(phi) v mathrmd Gamma_i\nright","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"For this example, an upwind flux will be used for F^*. Using a matrix formulation, the above equation can be written as:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"phi^n+1 = phi^n + M^-1(f_Omega - f_Gamma)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"where M^-1 is the inverse of the mass matrix, f_Omega the volumic flux term and f_Gamma the surfacic flux term.","category":"page"},{"location":"tutorial/linear_transport/#Commented-code","page":"Linear transport (DG)","title":"Commented code","text":"","category":"section"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Start by importing the necessary packages: Load the necessary packages","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"using Bcube\nusing LinearAlgebra\nusing WriteVTK","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Before all, to ease to ease the solution VTK output we will write a structure to store the vtk filename and the number of iteration; and a function that exports the solution on demand. Note the use of var_on_nodes_discontinuous to export the solution on the mesh nodes, respecting the discontinuous feature of the solution.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"mutable struct VtkHandler\n basename::Any\n ite::Any\n mesh::Any\n VtkHandler(basename, mesh) = new(basename, 0, mesh)\nend\n\nfunction append_vtk(vtk, u::Bcube.AbstractFEFunction, t)\n # Values on center\n values = var_on_nodes_discontinuous(u, vtk.mesh)\n\n # Write\n Bcube.write_vtk_discontinuous(\n vtk.basename,\n vtk.ite,\n t,\n vtk.mesh,\n Dict(\"u\" => (values, VTKPointData())),\n 1;\n append = vtk.ite > 0,\n )\n\n # Update counter\n vtk.ite += 1\nend","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"First, we define some physical and numerical constant parameters","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"const degree = 0 # Function-space degree (Taylor(0) = first order Finite Volume)\nconst c = [1.0, 0.0] # Convection velocity (must be a vector)\nconst nite = 100 # Number of time iteration(s)\nconst CFL = 1 # 0.1 for degree 1\nconst nx = 41 # Number of nodes in the x-direction\nconst ny = 41 # Number of nodes in the y-direction\nconst lx = 2.0 # Domain width\nconst ly = 2.0 # Domain height\nconst Δt = CFL * min(lx / nx, ly / ny) / norm(c) # Time step","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Then generate the mesh of a rectangle using Gmsh and read it","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"tmp_path = \"tmp.msh\"\ngen_rectangle_mesh(tmp_path, :quad; nx = nx, ny = ny, lx = lx, ly = ly, xc = 0.0, yc = 0.0)\nmesh = read_msh(tmp_path)\nrm(tmp_path)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"We can now init our VtkHandler","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"out_dir = joinpath(@__DIR__, \"../myout\")\nvtk = VtkHandler(joinpath(out_dir, \"linear_transport\"), mesh)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"As seen in the previous tutorial, the definition of trial and test spaces needs a mesh and a function space. Here, we select Taylor space, and build discontinuous FE spaces with it. Then an FEFunction, that will represent our solution, is created.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"fs = FunctionSpace(:Taylor, degree)\nU = TrialFESpace(fs, mesh, :discontinuous)\nV = TestFESpace(U)\nu = FEFunction(U)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Define measures for cell and interior face integrations","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Γ = InteriorFaceDomain(mesh)\nΓ_in = BoundaryFaceDomain(mesh, \"West\")\nΓ_out = BoundaryFaceDomain(mesh, (\"North\", \"East\", \"South\"))\n\ndΩ = Measure(CellDomain(mesh), 2 * degree + 1)\ndΓ = Measure(Γ, 2 * degree + 1)\ndΓ_in = Measure(Γ_in, 2 * degree + 1)\ndΓ_out = Measure(Γ_out, 2 * degree + 1)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"We will also need the face normals associated to the different face domains. Note that this operation is lazy, nΓ is just an abstract representation on face normals of Γ.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"nΓ = get_face_normals(Γ)\nnΓ_in = get_face_normals(Γ_in)\nnΓ_out = get_face_normals(Γ_out)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Let's move on to the bilinear and linear forms. First, the two easiest ones:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"m(u, v) = ∫(u ⋅ v)dΩ # Mass matrix\nl_Ω(v) = ∫((c * u) ⋅ ∇(v))dΩ # Volumic convective term","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"For the flux term, we first need to define a numerical flux. It is convenient to define it separately in a dedicated function. Here is the definition of simple upwind flux.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"function upwind(ui, uj, nij)\n cij = c ⋅ nij\n if cij > zero(cij)\n flux = cij * ui\n else\n flux = cij * uj\n end\n flux\nend","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"We then define the \"flux\" as the composition of the upwind function and the needed entries: namely the solution on the negative side of the face, the solution on the positive face, and the face normal. The orientation negative/positive is arbitrary, the only convention is that the face normals are oriented from the negative side to the positive side.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"flux = upwind ∘ (side⁻(u), side⁺(u), side⁻(nΓ))\nl_Γ(v) = ∫(flux * jump(v))dΓ","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Finally, we define what to perform on the \"two\" boundaries : inlet / oulet. On the inlet, we directly impose the flux with a user defined function that depends on the time (the input is an oscillating wave). On the outlet, we keep our upwind flux but we impose the ghost cell value.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"bc_in = t -> PhysicalFunction(x -> c .* cos(3 * x[2]) * sin(4 * t)) # flux\nl_Γ_in(v, t) = ∫(side⁻(bc_in(t)) ⋅ side⁻(nΓ_in) * side⁻(v))dΓ_in\nflux_out = upwind ∘ (side⁻(u), 0.0, side⁻(nΓ_out))\nl_Γ_out(v) = ∫(flux_out * side⁻(v))dΓ_out","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Assemble the (constant) mass matrix. The returned matrix is a sparse matrix. To simplify the tutorial, we will directly compute the inverse mass matrix. But note that way more performant strategies should be employed to solve such a problem (since we don't need the inverse, only the matrix-vector product).","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"M = assemble_bilinear(m, U, V)\ninvM = inv(Matrix(M)) #WARNING : really expensive !!!","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Let's also create three vectors to avoid allocating them at each time step","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"nd = get_ndofs(V)\nb_vol = zeros(nd)\nb_fac = similar(b_vol)\nrhs = similar(b_vol)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"The time loop is trivial : at each time step we compute the linear forms using the assemble_ methods, we complete the rhs, perform an explicit step and write the solution.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"t = 0.0\nfor i in 1:nite\n global t\n\n # Reset pre-allocated vectors\n b_vol .= 0.0\n b_fac .= 0.0\n\n # Compute linear forms\n assemble_linear!(b_vol, l_Ω, V)\n assemble_linear!(b_fac, l_Γ, V)\n assemble_linear!(b_fac, v -> l_Γ_in(v, t), V)\n assemble_linear!(b_fac, l_Γ_out, V)\n\n # Assemble rhs\n rhs .= Δt .* invM * (b_vol - b_fac)\n\n # Update solution\n u.dofValues .+= rhs\n\n # Update time\n t += Δt\n\n # Write to file\n append_vtk(vtk, u, t)\nend","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"And here is an animation of the result:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"![\"drawing\"](\"../assets/linear_transport.gif\")
","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"manual/operator/#LazyOperators","page":"LazyOperators","title":"LazyOperators","text":"","category":"section"},{"location":"manual/operator/","page":"LazyOperators","title":"LazyOperators","text":"WORK IN PROGRESS","category":"page"},{"location":"manual/integration/#Integration","page":"Integration","title":"Integration","text":"","category":"section"},{"location":"manual/integration/","page":"Integration","title":"Integration","text":"To compute an integral on a geometrical element, for instance a curved element, a variable substitution is used to compute the integral on the corresponding reference Shape. This variable substitution reads:","category":"page"},{"location":"manual/integration/","page":"Integration","title":"Integration","text":"int_Omega g(x) mathrmd Omega = int_hatOmega J(x) left(g circ F right)(hatx) mathrmd hatOmega","category":"page"},{"location":"manual/integration/","page":"Integration","title":"Integration","text":"where we recall that F is the reference to physical mapping and J is the determinant of the jacobian matrix of this mapping. Depending on the shape and element order, this determinant is either hard-coded or computed with ForwardDiff.","category":"page"},{"location":"manual/integration/","page":"Integration","title":"Integration","text":"Now, to compute the right side, i.e the integral on the reference shape, quadrature rules are applied:","category":"page"},{"location":"manual/integration/","page":"Integration","title":"Integration","text":"int_hatOmega g(hatx) mathrmd hatOmega = sum_i =1^N_q omega_i g(hatx_i)","category":"page"},{"location":"manual/integration/","page":"Integration","title":"Integration","text":"A specific procedure is applied to compute integrals on a face of a cell (i.e a surfacic integral on a face of a volumic element).","category":"page"},{"location":"example/linear_elasticity/#Linear-elasticity","page":"Linear elasticity","title":"Linear elasticity","text":"","category":"section"},{"location":"example/linear_elasticity/","page":"Linear elasticity","title":"Linear elasticity","text":"module linear_elasticity #hide\nprintln(\"Running linear elasticity API example...\") #hide\n\n# # Linear elasticity\n\nconst dir = string(@__DIR__, \"/\") # bcube/example dir\nusing Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing StaticArrays\n\n# function space (here we shall use Lagrange P1 elements) and quadrature degree.\nconst fspace = :Lagrange\nconst degree = 1 # FunctionSpace degree\nconst degquad = 2 * degree + 1\n\n# Input and output paths\nconst outputpath = dir * \"../myout/elasticity/\"\nconst meshpath = dir * \"../input/mesh/domainElast_tri.msh\"\n\n# Time stepping scheme params\nconst α = 0.05\nconst γ = 0.5 + α\nconst β = 0.25 * (1.0 + α)^2\n\n# Material parameters (Young's modulus, Poisson coefficient and deduced Lamé coefficients)\nconst ρ = 2500.0\nconst E = 200.0e9\nconst ν = 0.3\nconst λ = E * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))\nconst μ = E / (2.0 * (1.0 + ν))\n\n# Strain tensor and stress tensor (Hooke's law)\nϵ(u) = 0.5 * (∇(u) + transpose(∇(u)))\nσ(u) = λ * tr(ϵ(u)) * I + 2 * μ * ϵ(u)\n\nπ(u, v) = σ(u) ⊡ ϵ(v) # with the chosen contraction convention ϵ should be transposed, but as it is symmetric the expression remains correct\n\n# materialize for identity operator\nBcube.materialize(A::LinearAlgebra.UniformScaling, B) = A\n\n# Function that runs the steady case:\nfunction run_steady()\n # read mesh, the second argument specifies the spatial dimension\n mesh = read_msh(meshpath, 2)\n\n fs = FunctionSpace(fspace, degree)\n U_vec = TrialFESpace(\n fs,\n mesh,\n Dict(\"West\" => SA[0.0, 0.0], \"East\" => SA[1.0, 0.0]);\n size = 2,\n )\n V_vec = TestFESpace(U_vec)\n\n # Define measures for cell\n dΩ = Measure(CellDomain(mesh), degquad)\n\n # no volume force term\n f = PhysicalFunction(x -> SA[0.0, 0.0])\n\n # definition of bilinear and linear forms\n a(u, v) = ∫(π(u, v))dΩ\n l(v) = ∫(f ⋅ v)dΩ\n\n # solve using AffineFESystem\n sys = Bcube.AffineFESystem(a, l, U_vec, V_vec)\n ϕ = Bcube.solve(sys)\n\n Un = var_on_vertices(ϕ, mesh)\n # Write the obtained FE solution\n dict_vars = Dict(\"Displacement\" => (transpose(Un), VTKPointData()))\n mkpath(outputpath)\n write_vtk(outputpath * \"result_elasticity\", itime, t, mesh, dict_vars; append = false)\nend\n\n# Function that performs a time step using a Newmark α-HHT scheme\n# The scheme updates the acceleration G, the velocity V and the displacement U using the following formulas:\n#\n# M G +(1-α)A U + αA U0 = (1-α) L + α L0 = L (because here L is time independent)\n# V = V0 + (1-γ) Δt G0 + γ Δt G\n# U = U0 + Δt V0 + (0.5-β)*Δt^2 G0 + β Δt^2 G\n#\n# G is then computed by solving the linear system obtained by inserting the expressions for U and V in the equation for G.\nfunction Newmark_α_HHT(dt, L, A, Mat, U0, V0, G0)\n L1 = L - α * A * U0\n L2 = -(1.0 - α) * (A * U0 + dt * A * V0 + (0.5 - β) * dt * dt * A * G0)\n RHS = L1 .+ L2\n\n G = Mat \\ RHS\n V = V0 + (1.0 - γ) * dt * G0 + γ * dt * G\n U = U0 + dt * V0 + (0.5 - β) * dt * dt * G0 + β * dt * dt * G\n\n return U, V, G\nend\n\n# Function that runs the unsteady case:\nfunction run_unsteady()\n # read mesh, the second argument specifies the spatial dimension\n mesh = read_msh(meshpath, 2)\n\n fs = FunctionSpace(fspace, degree)\n U_vec = TrialFESpace(fs, mesh, Dict(\"West\" => SA[0.0, 0.0]); size = 2)\n V_vec = TestFESpace(U_vec)\n\n # Define measures for cell\n dΩ = Measure(CellDomain(mesh), degquad)\n Γ = BoundaryFaceDomain(mesh, (\"East\",))\n dΓ = Measure(Γ, degquad)\n\n # surface force to be applied on East boundary\n f = PhysicalFunction(x -> SA[100000.0, 1000.0])\n\n # Definition of bilinear and linear forms\n a(u, v) = ∫(π(u, v))dΩ\n m(u, v) = ∫(ρ * u ⋅ v)dΩ\n l(v) = ∫(side⁻(f) ⋅ side⁻(v))dΓ\n\n # Assemble matrices and vector\n M = assemble_bilinear(m, U_vec, V_vec)\n A = assemble_bilinear(a, U_vec, V_vec)\n L = assemble_linear(l, V_vec)\n\n # Apply homogeneous dirichlet on A and b\n Bcube.apply_homogeneous_dirichlet_to_vector!(L, U_vec, V_vec, mesh)\n Bcube.apply_dirichlet_to_matrix!((A, M), U_vec, V_vec, mesh)\n\n # Initialize solution\n ϕ = FEFunction(U_vec, 0.0)\n U0 = zeros(Bcube.get_ndofs(U_vec))\n V0 = zeros(Bcube.get_ndofs(U_vec))\n G0 = zeros(Bcube.get_ndofs(U_vec))\n\n # Write initial solution\n Un = var_on_vertices(ϕ, mesh)\n # Write the obtained FE solution\n dict_vars = Dict(\"Displacement\" => (transpose(Un), VTKPointData()))\n mkpath(outputpath)\n write_vtk(outputpath * \"result_elasticity\", 0, 0.0, mesh, dict_vars; append = false)\n\n # Time loop\n totalTime = 1.0e-3\n Δt = 1.0e-6\n itime = 0\n t = 0.0\n\n # Matrix for time stepping\n Mat = factorize(M + (1.0 - α) * (β * Δt * Δt * A))\n\n while t <= totalTime\n t += Δt\n itime = itime + 1\n @show t, itime\n\n # solve time step\n U, V, G = Newmark_α_HHT(Δt, L, A, Mat, U0, V0, G0)\n\n # Update solution\n U0 .= U\n V0 .= V\n G0 .= G\n\n set_dof_values!(ϕ, U)\n\n # Write solution\n if itime % 10 == 0\n Un = var_on_vertices(ϕ, mesh)\n # Write the obtained FE solution\n dict_vars = Dict(\"Displacement\" => (transpose(Un), VTKPointData()))\n write_vtk(\n outputpath * \"result_elasticity\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n # In order to use the warp function in paraview (solid is deformed using the displacement field)\n # the calculator filter has to be used with the following formula to reconstruct a 3D displacement field\n # with 0 z-component: Displacement_X*iHat+Displacement_Y*jHat+0.0*kHat\n end\n end\nend\n\n#run_steady()\nrun_unsteady()\n\nend #hide","category":"page"},{"location":"example/euler_naca_steady/#Euler-equations-on-a-NACA0012","page":"Euler equations on a NACA0012","title":"Euler equations on a NACA0012","text":"","category":"section"},{"location":"example/euler_naca_steady/","page":"Euler equations on a NACA0012","title":"Euler equations on a NACA0012","text":"module EulerNacaSteady #hide\nprintln(\"Running euler_naca_steady example...\") #hide\n# # Solve Euler equation around a NACA0012 airfoil\n\nusing Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing StaticArrays\nusing BenchmarkTools\nusing Roots\nusing SparseArrays\nusing Profile\nusing InteractiveUtils\nusing WriteVTK\nusing DifferentialEquations\nusing Symbolics\nusing SparseDiffTools\n\nconst dir = string(@__DIR__, \"/\")\n\nfunction compute_residual(qdof, Q, V, params)\n q = (FEFunction(Q, qdof)...,)\n\n # alias on measures\n dΓ = params.dΓ\n dΩ = params.dΩ\n dΓ_wall = params.dΓ_wall\n dΓ_farfield = params.dΓ_farfield\n\n # Allocate rhs vectors\n b_vol = zero(qdof)\n b_fac = zero(qdof)\n\n # compute volume residuals\n l_vol(v) = ∫(flux_Ω(q, v))dΩ\n assemble_linear!(b_vol, l_vol, V)\n\n # face normals for each face domain (lazy, no computation at this step)\n nΓ = get_face_normals(dΓ)\n nΓ_wall = get_face_normals(dΓ_wall)\n nΓ_farfield = get_face_normals(dΓ_farfield)\n\n # flux residuals from interior faces for all variables\n l_Γ(v) = ∫(flux_Γ(q, v, nΓ))dΓ\n assemble_linear!(b_fac, l_Γ, V)\n\n # flux residuals from bc faces for all variables\n l_Γ_wall(v) = ∫(flux_Γ_wall(q, v, nΓ_wall))dΓ_wall\n l_Γ_farfield(v) = ∫(flux_Γ_farfield(q, v, nΓ_farfield))dΓ_farfield\n assemble_linear!(b_fac, l_Γ_wall, V)\n assemble_linear!(b_fac, l_Γ_farfield, V)\n dQ = b_vol .- b_fac\n\n return dQ\nend\n\n\"\"\"\n flux_Ω(q, v)\n\nCompute volume residual using the lazy-operators approach\n\"\"\"\nflux_Ω(q, v) = _flux_Ω ∘ (q, map(∇, v))\n\nfunction _flux_Ω(q, ∇v)\n ρ, ρu, ρE = q\n ∇λ_ρ, ∇λ_ρu, ∇λ_ρE = ∇v\n γ = stateInit.γ\n\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = pressure(ρ, ρu, ρE, γ)\n\n flux_ρ = ρu\n flux_ρu = ρuu + p * I\n flux_ρE = (ρE + p) .* vel\n\n return return ∇λ_ρ ⋅ flux_ρ + ∇λ_ρu ⊡ flux_ρu + ∇λ_ρE ⋅ flux_ρE\nend\n\n\"\"\"\n flux_Γ(q, v, n)\n\nFlux at the interface is defined by a composition of two functions:\n* the input states at face sides which are needed for the riemann flux\n* `flux_roe` defines the Riemann flux (as usual)\n\"\"\"\nflux_Γ(q, v, n) = flux_roe ∘ (side⁻(q), side⁺(q), jump(v), side⁻(n))\n\n\"\"\"\n flux_roe(q⁻, q⁺, δv, n)\n\"\"\"\nfunction flux_roe(q⁻, q⁺, δv, n)\n γ = stateInit.γ\n nx, ny = n\n ρ1, (ρu1, ρv1), ρE1 = q⁻\n ρ2, (ρu2, ρv2), ρE2 = q⁺\n δλ_ρ1, δλ_ρu1, δλ_ρE1 = δv\n\n ρ1 = max(eps(ρ1), ρ1)\n ρ2 = max(eps(ρ2), ρ2)\n\n # Closure\n u1 = ρu1 / ρ1\n v1 = ρv1 / ρ1\n u2 = ρu2 / ρ2\n v2 = ρv2 / ρ2\n p1 = pressure(ρ1, SA[ρu1, ρv1], ρE1, γ)\n p2 = pressure(ρ2, SA[ρu2, ρv2], ρE2, γ)\n\n H2 = (γ / (γ - 1)) * p2 / ρ2 + (u2 * u2 + v2 * v2) / 2.0\n H1 = (γ / (γ - 1)) * p1 / ρ1 + (u1 * u1 + v1 * v1) / 2.0\n\n R = √(ρ1 / ρ2)\n invR1 = 1.0 / (R + 1)\n uAv = (R * u1 + u2) * invR1\n vAv = (R * v1 + v2) * invR1\n Hav = (R * H1 + H2) * invR1\n cAv = √(abs((γ - 1) * (Hav - (uAv * uAv + vAv * vAv) / 2.0)))\n ecAv = (uAv * uAv + vAv * vAv) / 2.0\n\n λ1 = nx * uAv + ny * vAv\n λ3 = λ1 + cAv\n λ4 = λ1 - cAv\n\n d1 = ρ1 - ρ2\n d2 = ρ1 * u1 - ρ2 * u2\n d3 = ρ1 * v1 - ρ2 * v2\n d4 = ρE1 - ρE2\n\n # computation of the centered part of the flux\n flux_ρ = nx * ρ2 * u2 + ny * ρ2 * v2\n flux_ρu = nx * p2 + flux_ρ * u2\n flux_ρv = ny * p2 + flux_ρ * v2\n flux_ρE = H2 * flux_ρ\n\n # Temp variables\n rc1 = (γ - 1) / cAv\n rc2 = (γ - 1) / cAv / cAv\n uq41 = ecAv / cAv + cAv / (γ - 1)\n uq42 = nx * uAv + ny * vAv\n\n fdc1 = max(λ1, 0.0) * (d1 + rc2 * (-ecAv * d1 + uAv * d2 + vAv * d3 - d4))\n fdc2 = max(λ1, 0.0) * ((nx * vAv - ny * uAv) * d1 + ny * d2 - nx * d3)\n fdc3 =\n max(λ3, 0.0) * (\n (-uq42 * d1 + nx * d2 + ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n fdc4 =\n max(λ4, 0.0) * (\n (uq42 * d1 - nx * d2 - ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n\n duv1 = fdc1 + (fdc3 + fdc4) / cAv\n duv2 = uAv * fdc1 + ny * fdc2 + (uAv / cAv + nx) * fdc3 + (uAv / cAv - nx) * fdc4\n duv3 = vAv * fdc1 - nx * fdc2 + (vAv / cAv + ny) * fdc3 + (vAv / cAv - ny) * fdc4\n duv4 =\n ecAv * fdc1 +\n (ny * uAv - nx * vAv) * fdc2 +\n (uq41 + uq42) * fdc3 +\n (uq41 - uq42) * fdc4\n\n flux_ρ += duv1\n flux_ρu += duv2\n flux_ρv += duv3\n flux_ρE += duv4\n\n return (δλ_ρ1 ⋅ flux_ρ + δλ_ρu1 ⋅ SA[flux_ρu, flux_ρv] + δλ_ρE1 ⋅ flux_ρE)\nend\n\n\"\"\"\n flux_Γ_farfield(q, v, n)\n\nCompute `Roe` flux on boundary face by imposing\n`stateBcFarfield.u_in` on `side_p`\n\"\"\"\nflux_Γ_farfield(q, v, n) = flux_roe ∘ (side⁻(q), stateBcFarfield.u_inf, side⁻(v), side⁻(n))\n\n\"\"\"\n flux_Γ_wall(q, v, n)\n\"\"\"\nflux_Γ_wall(q, v, n) = _flux_Γ_wall ∘ (side⁻(q), side⁻(v), side⁻(n))\n\nfunction _flux_Γ_wall(q⁻, v⁻, n)\n γ = stateInit.γ\n ρ1, ρu1, ρE1 = q⁻\n λ_ρ1, λ_ρu1, λ_ρE1 = v⁻\n\n p1 = pressure(ρ1, ρu1, ρE1, γ)\n\n flux_ρ = zero(ρ1)\n flux_ρu = p1 * n\n flux_ρE = zero(ρE1)\n\n return (λ_ρ1 ⋅ flux_ρ + λ_ρu1 ⋅ flux_ρu + λ_ρE1 ⋅ flux_ρE)\nend\n\nfunction sparse2vtk(\n a::AbstractSparseMatrix,\n name::String = string(@__DIR__, \"/../myout/sparse\"),\n)\n vtk_write_array(name, Array(a), \"my_property_name\")\nend\n\nmutable struct VtkHandler\n basename::String\n basename_residual::String\n ite::Int\n VtkHandler(basename) = new(basename, basename * \"_residual\", 0)\nend\n\n\"\"\"\n Write solution (at cell centers) to vtk\n Wrapper for `write_vtk`\n\"\"\"\nfunction append_vtk(vtk, mesh, vars, t, params; res = nothing)\n ρ, ρu, ρE = vars\n\n # Mean cell values\n # name2val_mean = (;zip(get_name.(vars), mean_values.(vars, degquad))...)\n # p_mean = pressure.(name2val_mean[:ρ], name2val_mean[:ρu], name2val_mean[:ρE], params.stateInit.γ)\n\n vtk_degree = maximum(x -> get_degree(Bcube.get_function_space(get_fespace(x))), vars)\n vtk_degree = max(1, mesh_degree, vtk_degree)\n _ρ = var_on_nodes_discontinuous(ρ, mesh, vtk_degree)\n _ρu = var_on_nodes_discontinuous(ρu, mesh, vtk_degree)\n _ρE = var_on_nodes_discontinuous(ρE, mesh, vtk_degree)\n\n Cp = pressure_coefficient.(_ρ, _ρu, _ρE)\n Ma = mach.(_ρ, _ρu, _ρE)\n dict_vars_dg = Dict(\n \"rho\" => (_ρ, VTKPointData()),\n \"rhou\" => (_ρu, VTKPointData()),\n \"rhoE\" => (_ρE, VTKPointData()),\n \"Cp\" => (Cp, VTKPointData()),\n \"Mach\" => (Ma, VTKPointData()),\n \"rho_mean\" => (get_values(Bcube.cell_mean(ρ, params.dΩ)), VTKCellData()),\n \"rhou_mean\" => (get_values(Bcube.cell_mean(ρu, params.dΩ)), VTKCellData()),\n \"rhoE_mean\" => (get_values(Bcube.cell_mean(ρE, params.dΩ)), VTKCellData()),\n \"lim_rho\" => (get_values(params.limρ), VTKCellData()),\n \"lim_all\" => (get_values(params.limAll), VTKCellData()),\n )\n Bcube.write_vtk_discontinuous(\n vtk.basename * \"_DG\",\n vtk.ite,\n t,\n mesh,\n dict_vars_dg,\n vtk_degree;\n append = vtk.ite > 0,\n )\n\n _ρ_wall = var_on_bnd_nodes_discontinuous(ρ, params.Γ_wall, vtk_degree)\n _ρu_wall = var_on_bnd_nodes_discontinuous(ρu, params.Γ_wall, vtk_degree)\n _ρE_wall = var_on_bnd_nodes_discontinuous(ρE, params.Γ_wall, vtk_degree)\n\n Cp_wall = pressure_coefficient.(_ρ_wall, _ρu_wall, _ρE_wall)\n Ma_wall = pressure_coefficient.(_ρ_wall, _ρu_wall, _ρE_wall)\n\n dict_vars_wall = Dict(\n \"rho\" => (_ρ_wall, VTKPointData()),\n \"rhou\" => (_ρu_wall, VTKPointData()),\n \"rhoE\" => (_ρE_wall, VTKPointData()),\n \"Cp\" => (Cp_wall, VTKPointData()),\n \"Mach\" => (Ma_wall, VTKPointData()),\n )\n Bcube.write_vtk_bnd_discontinuous(\n vtk.basename * \"_bnd_DG\",\n 1,\n 0.0,\n params.Γ_wall,\n dict_vars_wall,\n vtk_degree;\n append = false,\n )\n\n #residual:\n if !isa(res, Nothing)\n vtkfile = vtk_grid(vtk.basename_residual, Float64.(res.iter), [0.0, 1.0])\n for (k, valₖ) in enumerate(res.val)\n vtkfile[\"res_\" * string(k), VTKPointData()] = [valₖ valₖ]\n end\n vtk_save(vtkfile)\n end\n\n # Update counter\n vtk.ite += 1\n\n return nothing\nend\n\nfunction init!(q, dΩ, initstate)\n AoA = initstate.AoA\n Minf = initstate.M_inf\n Pinf = initstate.P_inf\n Tinf = initstate.T_inf\n r = initstate.r_gas\n γ = initstate.γ\n\n ρinf = Pinf / r / Tinf\n ainf = √(γ * r * Tinf)\n Vinf = Minf * ainf\n ρVxinf = ρinf * Vinf * cos(AoA)\n ρVyinf = ρinf * Vinf * sin(AoA)\n ρEinf = Pinf / (γ - 1) + 0.5 * ρinf * Vinf^2\n\n ρ0 = PhysicalFunction(x -> ρinf)\n ρu0 = PhysicalFunction(x -> SA[ρVxinf, ρVyinf])\n ρE0 = PhysicalFunction(x -> ρEinf)\n projection_l2!(q, (ρ0, ρu0, ρE0), dΩ)\n return nothing\nend\n\nfunction main(stateInit, stateBcFarfield, degree)\n @show degree, degquad\n\n mesh = read_msh(dir * \"../input/mesh/naca0012_o\" * string(mesh_degree) * \".msh\", 2)\n scale!(mesh, 1.0 / 0.5334)\n\n dimcar = compute_dimcar(mesh)\n\n DMPrelax = DMPcurv₀ .* dimcar .^ 2\n\n # Then we create a `NamedTuple` to hold the simulation parameters.\n params = (\n degquad = degquad,\n stateInit = stateInit,\n stateBcFarfield = stateBcFarfield,\n DMPrelax = DMPrelax,\n )\n\n # Define measures for cell and interior face integrations\n dΩ = Measure(CellDomain(mesh), degquad)\n dΓ = Measure(InteriorFaceDomain(mesh), degquad)\n\n # Declare boundary conditions and\n # create associated domains and measures\n Γ_wall = BoundaryFaceDomain(mesh, (\"NACA\",))\n Γ_farfield = BoundaryFaceDomain(mesh, (\"FARFIELD\",))\n dΓ_wall = Measure(Γ_wall, degquad)\n dΓ_farfield = Measure(Γ_farfield, degquad)\n\n params = (\n params...,\n Γ_wall = Γ_wall,\n dΓ = dΓ,\n dΩ = dΩ,\n dΓ_wall = dΓ_wall,\n dΓ_farfield = dΓ_farfield,\n )\n\n qLowOrder = nothing\n\n for deg in 0:degree\n params = (params..., degree = deg)\n\n fs = FunctionSpace(fspace, deg)\n Q_sca = TrialFESpace(fs, mesh, :discontinuous; size = 1) # DG, scalar\n Q_vec = TrialFESpace(fs, mesh, :discontinuous; size = 2) # DG, vectoriel\n V_sca = TestFESpace(Q_sca)\n V_vec = TestFESpace(Q_vec)\n Q = MultiFESpace(Q_sca, Q_vec, Q_sca)\n V = MultiFESpace(V_sca, V_vec, V_sca)\n\n q = FEFunction(Q)\n\n # select an initial configurations:\n if deg == 0\n init!(q, mesh, stateInit)\n else\n println(\"Start projection\")\n projection_l2!(q, qLowOrder, dΩ)\n println(\"End projection\")\n end\n\n # create CellData to store limiter values\n limρ = Bcube.MeshCellData(ones(ncells(mesh)))\n limAll = Bcube.MeshCellData(ones(ncells(mesh)))\n params = (params..., limρ = limρ, limAll = limAll)\n\n # Init vtk handler\n mkpath(outputpath)\n vtk = VtkHandler(\n outputpath * \"euler_naca_mdeg\" * string(mesh_degree) * \"_deg\" * string(deg),\n )\n\n # Init time\n time = 0.0\n\n # Save initial solution\n append_vtk(vtk, mesh, q, time, params)\n\n # Build the cache and store everything you want to compute only once (such as the mass matrice inverse...)\n\n cache = ()\n # Allocate buffer for compute_residual\n b_vol = zeros(Bcube.get_ndofs(Q))\n b_fac = zeros(Bcube.get_ndofs(Q))\n cache = (cache..., b_vol = b_vol, b_fac = b_fac)\n\n cache = (\n cache...,\n cacheCellMean = Bcube.build_cell_mean_cache(q, dΩ),\n mass = factorize(Bcube.build_mass_matrix(Q, V, dΩ)),\n mass_sca = factorize(Bcube.build_mass_matrix(Q_sca, V_sca, dΩ)),\n mass_vec = factorize(Bcube.build_mass_matrix(Q_vec, V_vec, dΩ)),\n )\n\n time, q = steady_solve!(Q, V, q, mesh, params, cache, vtk, deg)\n append_vtk(vtk, mesh, q, time, params)\n println(\"end steady_solve for deg=\", deg, \" !\")\n\n deg < degree && (qLowOrder = deepcopy(q))\n end\n return nothing\nend\n\nfunction steady_solve!(Q, V, q, mesh, params, cache, vtk, deg)\n counter = [0]\n q0 = deepcopy(get_dof_values(q))\n ode_params =\n (Q = Q, V = V, params = params, cache = cache, counter = counter, vtk = vtk)\n\n rhs!(dq, q, p, t) = dq .= compute_residual(q, p.Q, p.V, p.params)\n\n # compute sparsity pattern and coloring\n println(\"computing jacobian cache...\")\n if withbench\n _rhs!(dq, q) = rhs!(dq, q, ode_params, 0.0)\n @btime $_rhs!(similar($q0), $q0)\n q_bench = FEFunction(Q, q0)\n @btime $apply_limitation!($q_bench, $ode_params)\n @show length(q0)\n end\n\n #sparsity_pattern = Symbolics.jacobian_sparsity(_rhs!, similar(Q0), Q0)\n #tjac = @elapsed Symbolics.jacobian_sparsity(_rhs!, similar(Q0), Q0)\n #@show tjac\n sparsity_pattern = Bcube.build_jacobian_sparsity_pattern(Q, mesh)\n println(\"sparsity pattern computed !\")\n display(sparsity_pattern)\n colors = matrix_colors(sparsity_pattern)\n println(\"coloring done!\")\n @show maximum(colors)\n\n ode = ODEFunction(\n rhs!;\n mass_matrix = Bcube.build_mass_matrix(Q, V, params.dΩ),\n jac_prototype = sparsity_pattern,\n colorvec = colors,\n )\n\n Tfinal = Inf\n problem = ODEProblem(ode, q0, (0.0, Tfinal), ode_params)\n timestepper = ImplicitEuler(; nlsolve = NLNewton(; max_iter = 20))\n\n cb_cache = DiscreteCallback(always_true, update_cache!; save_positions = (false, false))\n cb_vtk = DiscreteCallback(always_true, output_vtk; save_positions = (false, false))\n cb_steady = TerminateSteadyState(1e-6, 1e-6, condition_steadystate)\n\n error = 1e-1\n\n sol = solve(\n problem,\n timestepper;\n initializealg = NoInit(),\n adaptive = true,\n abstol = error,\n reltol = error,\n progress = false,\n progress_steps = 1000,\n save_everystep = false,\n save_start = false,\n save_end = false,\n isoutofdomain = isoutofdomain,\n callback = CallbackSet(cb_cache, cb_vtk, cb_steady),\n )\n\n set_dof_values!(q, sol.u[end])\n return sol.t[end], q\nend\n\nalways_true(args...) = true\n\nfunction isoutofdomain(dof, p, t)\n any(isnan, dof) && return true\n\n q = FEFunction(p.Q, dof)\n q_mean = map(get_values, Bcube.cell_mean(q, p.cache.cacheCellMean))\n p_mean = pressure.(q_mean..., stateInit.γ)\n\n negative_ρ = any(x -> x < 0, q_mean[1])\n negative_p = any(x -> x < 0, p_mean)\n isout = negative_ρ || negative_p\n isout && @show negative_ρ, negative_p\n return isout\nend\n\nfunction update_cache!(integrator)\n Q = integrator.p.Q\n Q1, = Q\n deg = get_degree(Bcube.get_function_space(Q1))\n println(\n \"deg=\",\n deg,\n \" update_cache! : iter=\",\n integrator.p.counter[1],\n \" dt=\",\n integrator.dt,\n )\n\n q = FEFunction(integrator.p.Q, integrator.u)\n limiter_projection && apply_limitation!(q, integrator.p)\n return nothing\nend\n\nfunction output_vtk(integrator)\n u_modified!(integrator, false)\n mesh = get_mesh(get_domain(integrator.p.params.dΩ))\n q = FEFunction(integrator.p.Q, integrator.u)\n counter = integrator.p.counter\n counter .+= 1\n if (counter[1] % nout == 0)\n println(\"output_vtk \", counter[1])\n append_vtk(integrator.p.vtk, mesh, q, integrator.t, integrator.p.params)\n end\n return nothing\nend\n\nfunction condition_steadystate(integrator, abstol, reltol, min_t)\n u_modified!(integrator, false)\n if DiffEqBase.isinplace(integrator.sol.prob)\n testval = first(get_tmp_cache(integrator))\n @. testval = (integrator.u - integrator.uprev) / (integrator.t - integrator.tprev)\n else\n testval = (integrator.u - integrator.uprev) / (integrator.t - integrator.tprev)\n end\n\n if typeof(integrator.u) <: Array\n any(\n abs(d) > abstol && abs(d) > reltol * abs(u) for (d, abstol, reltol, u) in\n zip(testval, Iterators.cycle(abstol), Iterators.cycle(reltol), integrator.u)\n ) && (return false)\n else\n any((abs.(testval) .> abstol) .& (abs.(testval) .> reltol .* abs.(integrator.u))) &&\n (return false)\n end\n\n if min_t === nothing\n return true\n else\n return integrator.t >= min_t\n end\nend\n\n\"\"\"\nCompute the characteristic dimension of each cell of `mesh`:\ndimcar = (cell volume) / (cell surface)\n\n# TODO :\nto be moved to Bcube\n\"\"\"\nfunction compute_dimcar(mesh)\n fs = FunctionSpace(:Lagrange, 0)\n V = TestFESpace(fs, mesh; size = 1, isContinuous = false)\n\n # Define measures for cell and interior face integrations\n dΩ = Measure(CellDomain(mesh), degquad)\n dΓ = Measure(InteriorFaceDomain(mesh), degquad)\n dΓ_bc = Measure(BoundaryFaceDomain(mesh), degquad)\n\n f1 = PhysicalFunction(x -> 1.0)\n l(v) = ∫(f1 ⋅ v)dΩ\n l_face(v, dω) = ∫(side⁻(f1) ⋅ side⁻(v) + side⁺(f1) ⋅ side⁺(v))dω\n\n vol = assemble_linear(l, V)\n surf = assemble_linear(Base.Fix2(l_face, dΓ), V)\n surf += assemble_linear(Base.Fix2(l_face, dΓ_bc), V)\n return vol ./ surf\nend\n\n\"\"\"\nReferences:\n* Xiangxiong Zhang, Chi-Wang Shu, On positivity-preserving high order discontinuous\n Galerkin schemes for compressible Euler equations on rectangular meshes,\n Journal of Computational Physics, Volume 229, Issue 23, 2010.\n https://doi.org/10.1016/j.jcp.2010.08.016\n* Zhang, X., Xia, Y. & Shu, CW. Maximum-Principle-Satisfying and Positivity-Preserving\n High Order Discontinuous Galerkin Schemes for Conservation Laws on Triangular Meshes.\n J Sci Comput 50, 29–62 (2012). https://doi.org/10.1007/s10915-011-9472-8\n\"\"\"\nfunction apply_limitation!(q::Bcube.AbstractFEFunction, ode_params)\n params = ode_params.params\n cache = ode_params.cache\n mesh = get_mesh(get_domain(params.dΩ))\n ρ, ρu, ρE = q\n\n ρ_mean, ρu_mean, ρE_mean = Bcube.cell_mean(q, cache.cacheCellMean)\n\n _limρ, ρ_proj = linear_scaling_limiter(\n ρ,\n params.dΩ;\n bounds = (ρmin₀, ρmax₀),\n DMPrelax = params.DMPrelax,\n mass = cache.mass_sca,\n )\n\n op_t = limiter_param_p ∘ (ρ_proj, ρu, ρE, ρ_mean, ρu_mean, ρE_mean)\n t = Bcube._minmax_cells(op_t, mesh, Val(params.degquad))\n tmin = Bcube.MeshCellData(getindex.(t, 1))\n\n if eltype(_limρ) == eltype(params.limρ) # skip Dual number case\n set_values!(params.limρ, get_values(_limρ))\n set_values!(params.limAll, get_values(tmin))\n end\n\n limited_var(u, ū, lim_u) = ū + lim_u * (u - ū)\n projection_l2!(ρ, limited_var(ρ_proj, ρ_mean, tmin), params.dΩ; mass = cache.mass_sca)\n projection_l2!(ρu, limited_var(ρu, ρu_mean, tmin), params.dΩ; mass = cache.mass_vec)\n projection_l2!(ρE, limited_var(ρE, ρE_mean, tmin), params.dΩ; mass = cache.mass_sca)\n return nothing\nend\n\nfunction limiter_param_p(ρ̂, ρu, ρE, ρ_mean, ρu_mean, ρE_mean)\n γ = stateInit.γ\n p = pressure(ρ̂, ρu, ρE, γ)\n\n if p ≥ pmin₀\n t = 1.0\n else\n @show p, ρ̂, ρu, ρE\n @show ρ_mean, ρu_mean, ρE_mean\n @show pressure(ρ_mean, ρu_mean, ρE_mean, γ)\n if pressure(ρ_mean, ρu_mean, ρE_mean, γ) > pmin₀\n fₜ =\n t ->\n pressure(\n t * ρ̂ + (1 - t) * ρ_mean,\n t * ρu + (1 - t) * ρu_mean,\n t * ρE + (1 - t) * ρE_mean,\n γ,\n ) - pmin₀\n bounds = (0.0, 1.0)\n t = find_zero(fₜ, bounds, Bisection())\n else\n t = NaN\n println(\"t = NaN\")\n end\n end\n\n return t\nend\n\nfunction pressure(ρ::Number, ρu::AbstractVector, ρE::Number, γ)\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = (γ - 1) * (ρE - tr(ρuu) / 2)\n return p\nend\n\ncompute_Pᵢ(P, γ, M) = P * (1 + 0.5 * (γ - 1) * M^2)^(γ / (γ - 1))\ncompute_Tᵢ(T, γ, M) = T * (1 + 0.5 * (γ - 1) * M^2)\n\nfunction bc_state_farfield(AoA, M, P, T, r, γ)\n a = √(γ * r * T)\n vn = M * a\n ρ = P / r / T\n ρu = SA[ρ * vn * cos(AoA), ρ * vn * sin(AoA)]\n ρE = P / (γ - 1) + 0.5 * ρ * vn^2\n return (ρ, ρu, ρE)\nend\n\nfunction pressure_coefficient(ρ, ρu, ρE)\n (pressure(ρ, ρu, ρE, stateInit.γ) - stateInit.P_inf) /\n (stateBcFarfield.Pᵢ_inf - stateInit.P_inf)\nend\nfunction mach(ρ, ρu, ρE)\n norm(ρu ./ ρ) / √(stateInit.γ * max(0.0, pressure(ρ, ρu, ρE, stateInit.γ) / ρ))\nend\n\nconst degreemax = 2 # Function-space degree\nconst mesh_degree = 2\nconst fspace = :Lagrange\nconst limiter_projection = true\nconst ρmin₀ = 1.0e-8\nconst ρmax₀ = 1.0e+10\nconst pmin₀ = 1.0e-8\nconst pmax₀ = 1.0e+10\nconst DMPcurv₀ = 10.0e3\nconst withbench = false\n\nconst stateInit = (\n AoA = deg2rad(1.25),\n M_inf = 0.8,\n P_inf = 101325.0,\n T_inf = 275.0,\n r_gas = 287.0,\n γ = 1.4,\n)\nconst nite_max = 300 #300000 # Number of time iteration(s)\nconst nout = 1 # number of step between two vtk outputs\nconst mass_matrix_in_solve = true\nconst degquad = 6\nconst outputpath = string(@__DIR__, \"/../myout/euler_naca_steady/\")\n\nconst stateBcFarfield = (\n AoA = stateInit.AoA,\n M_inf = stateInit.M_inf,\n Pᵢ_inf = compute_Pᵢ(stateInit.P_inf, stateInit.γ, stateInit.M_inf),\n Tᵢ_inf = compute_Tᵢ(stateInit.T_inf, stateInit.γ, stateInit.M_inf),\n u_inf = bc_state_farfield(\n stateInit.AoA,\n stateInit.M_inf,\n stateInit.P_inf,\n stateInit.T_inf,\n stateInit.r_gas,\n stateInit.γ,\n ),\n r_gas = stateInit.r_gas,\n γ = stateInit.γ,\n)\n\nmain(stateInit, stateBcFarfield, degreemax)\n\nend #hide","category":"page"},{"location":"manual/geometry/#Geometry-and-mesh","page":"Geometry and mesh","title":"Geometry and mesh","text":"","category":"section"},{"location":"manual/geometry/","page":"Geometry and mesh","title":"Geometry and mesh","text":"A Mesh is a set basically of nodes (Node), a set of entities (the mesh elements) and a list of connectivies that link the entities between themselves and with the nodes.","category":"page"},{"location":"manual/geometry/","page":"Geometry and mesh","title":"Geometry and mesh","text":"In Bcube every mesh entity has corresponding reference Shape, a simplified or canonical representation of this element. A 1D line is mapped on the [-1,1] segment, and a rectangle is mapped on a square for instance. On these reference shapes, (almost) everything is known : the vertices location, the area, the quadrature points... Hence in Bcube we always compute things in the reference shape. For \"Lagrange\" elements (such as Bar*_t, Tri*_t, Quad*_t, Tetra*_t, Hexa*_t, Penta*_t etc), the mapping from the reference shape to a geometrical element is directly obtained from the corresponding Lagrange polynomials and the element node coordinates. Given a geometrical element with n nodes M_i, the mapping reads:","category":"page"},{"location":"manual/geometry/","page":"Geometry and mesh","title":"Geometry and mesh","text":"F(xi) = sum_i=1^n hatlambda_i(xi)M_i","category":"page"},{"location":"manual/geometry/","page":"Geometry and mesh","title":"Geometry and mesh","text":"where (lambda)_i are the Lagrange polynomials whose order matches the element order.","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = BcubeTutorials","category":"page"},{"location":"#Bcube","page":"Home","title":"Bcube","text":"","category":"section"},{"location":"#Purpose-of-Bcube","page":"Home","title":"Purpose of Bcube","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Bcube is a Julia library providing tools for the spatial discretization of partial differential equation(s) (PDE). The main objectives are:","category":"page"},{"location":"","page":"Home","title":"Home","text":"to provide a set of tools to quickly assemble an algorithm solving partial differential equation(s) (so the main objective is to help building prototypes without thinking about the numerical core)\nto be completed : efficient/performant PDE resolution?","category":"page"},{"location":"","page":"Home","title":"Home","text":"This documentation is organised as follow. Checkout the tutorials to see what Bcube is capable of and/or quickly learn how to use it. Then, some more elaborated examples are provided to demonstrate the library capabilities. The \"Manual\" part explains how the core is organized. Finally, the \"API\" section is the low level code documentation.","category":"page"},{"location":"#Writing-documentation","page":"Home","title":"Writing documentation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"To write documentation for Bcube, Julia's guidelines should be followed : https://docs.julialang.org/en/v1/manual/documentation/. Moreover, this project tries to apply the SciML Style Guide.","category":"page"},{"location":"#Conventions","page":"Home","title":"Conventions","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This documentation follows the following notation or naming conventions:","category":"page"},{"location":"","page":"Home","title":"Home","text":"coordinates inside a reference frame are noted hatx haty or xi eta while coordinates in the physical frame are noted xy\nwhen talking about a mapping, F or sometimes F_rp designates the mapping from the reference element to the physical element. On the other side, F^-1 or sometimes F_pr designates the physical element to the reference element mapping.\n\"dof\" means \"degree of freedom\"","category":"page"}]
+[{"location":"tutorial/heat_equation/#Heat-equation-(FE)","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"In this tutorial, the heat equation (first steady and then unsteady) is solved using finite-elements.","category":"page"},{"location":"tutorial/heat_equation/#Theory","page":"Heat equation (FE)","title":"Theory","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"This example shows how to solve the heat equation with eventually variable physical properties in steady and unsteady formulations:","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":" rho C_p partial_t u - nabla ( lambda u) = f","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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:","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":" forall v in tildeH^1_0(Omega) underbraceint_Omega partial_t u v dx_m(partial_t uv) + underbraceint_Omega nabla u nabla v dx_a(uv) = underbraceint_Omega f v dx_l(v)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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.","category":"page"},{"location":"tutorial/heat_equation/#Steady-case","page":"Heat equation (FE)","title":"Steady case","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"As usual, start by importing the necessary packages.","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"using Bcube\nusing LinearAlgebra\nusing WriteVTK","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"First we define some physical and numerical constants","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"const htc = 100.0 # Heat transfer coefficient (bnd cdt)\nconst Tr = 268.0 # Recovery temperature (bnd cdt)\nconst phi = 100.0\nconst q = 1500.0\nconst λ = 100.0\nconst η = λ\nconst ρCp = 100.0 * 200.0\nconst degree = 2\nconst outputpath = joinpath(@__DIR__, \"../myout/heat_equation/\")","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Read 2D mesh","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mesh_path = joinpath(@__DIR__, \"../input/mesh/domainSquare_tri.msh\")\nmesh = read_msh(mesh_path)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Build function space and associated Trial and Test FE spaces. We impose a Dirichlet condition with a temperature of 260K on boundary \"West\"","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"fs = FunctionSpace(:Lagrange, degree)\nU = TrialFESpace(fs, mesh, Dict(\"West\" => 260.0))\nV = TestFESpace(U)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Define measures for cell integration","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"dΩ = Measure(CellDomain(mesh), 2 * degree + 1)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Define bilinear and linear forms","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"a(u, v) = ∫(η * ∇(u) ⋅ ∇(v))dΩ\nl(v) = ∫(q * v)dΩ","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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.","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"sys = AffineFESystem(a, l, U, V)\nϕ = solve(sys)\nTcn = var_on_centers(ϕ, mesh)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute analytical solution for comparison. Apply the analytical solution on mesh centers","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"T_analytical = x -> 260.0 + (q / λ) * x[1] * (1.0 - 0.5 * x[1])\nTca = map(T_analytical, get_cell_centers(mesh))","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Write both the obtained FE solution and the analytical solution to a vtk file.","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mkpath(outputpath)\ndict_vars =\n Dict(\"Temperature\" => (Tcn, VTKCellData()), \"Temperature_a\" => (Tca, VTKCellData()))\nwrite_vtk(outputpath * \"result_steady_heat_equation\", 0, 0.0, mesh, dict_vars)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute and display the error","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"@show norm(Tcn .- Tca, Inf) / norm(Tca, Inf)","category":"page"},{"location":"tutorial/heat_equation/#Unsteady-case","page":"Heat equation (FE)","title":"Unsteady case","text":"","category":"section"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"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:","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"totalTime = 100.0\nΔt = 0.1","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Read a slightly different mesh","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mesh_path = joinpath(@__DIR__, \"../input/mesh/domainSquare_tri_2.msh\")\nmesh = read_msh(mesh_path)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"The rest is similar to the steady case","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"fs = FunctionSpace(:Lagrange, degree)\nU = TrialFESpace(fs, mesh, Dict(\"West\" => 260.0))\nV = TestFESpace(U)\ndΩ = Measure(CellDomain(mesh), 2 * degree + 1)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute matrices associated to bilinear and linear forms, and assemble","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"a(u, v) = ∫(η * ∇(u) ⋅ ∇(v))dΩ\nm(u, v) = ∫(ρCp * u ⋅ v)dΩ\nl(v) = ∫(q * v)dΩ\n\nA = assemble_bilinear(a, U, V)\nM = assemble_bilinear(m, U, V)\nL = assemble_linear(l, V)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Compute a vector of dofs whose values are zeros everywhere except on dofs lying on a Dirichlet boundary, where they take the Dirichlet value","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Ud = assemble_dirichlet_vector(U, V, mesh)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Apply lift","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"L = L - A * Ud","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Apply homogeneous dirichlet condition","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"apply_homogeneous_dirichlet_to_vector!(L, U, V, mesh)\napply_dirichlet_to_matrix!((A, M), U, V, mesh)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Form time iteration matrix (note that this is bad for performance since up to now, M and A are sparse matrices)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Miter = factorize(M + Δt * A)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Init the solution with a constant temperature of 260K","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"ϕ = FEFunction(U, 260.0)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Write initial solution to a file","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"mkpath(outputpath)\ndict_vars = Dict(\"Temperature\" => (var_on_centers(ϕ, mesh), VTKCellData()))\nwrite_vtk(outputpath * \"result_unsteady_heat_equation\", 0, 0.0, mesh, dict_vars)","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"Time loop","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"itime = 0\nt = 0.0\nwhile t <= totalTime\n global t, itime\n t += Δt\n itime = itime + 1\n @show t, itime\n\n # Compute rhs\n rhs = Δt * L + M * (get_dof_values(ϕ) .- Ud)\n\n # Invert system and apply inverse shift\n set_dof_values!(ϕ, Miter \\ rhs .+ Ud)\n\n # Write solution (every 10 iterations)\n if itime % 10 == 0\n dict_vars = Dict(\"Temperature\" => (var_on_centers(ϕ, mesh), VTKCellData()))\n write_vtk(\n outputpath * \"result_unsteady_heat_equation\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n end\nend\n","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"","category":"page"},{"location":"tutorial/heat_equation/","page":"Heat equation (FE)","title":"Heat equation (FE)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"example/euler_naca_steady/#Euler-equations-on-a-NACA0012","page":"Euler equations on a NACA0012","title":"Euler equations on a NACA0012","text":"","category":"section"},{"location":"example/euler_naca_steady/","page":"Euler equations on a NACA0012","title":"Euler equations on a NACA0012","text":"module EulerNacaSteady #hide\nprintln(\"Running euler_naca_steady example...\") #hide\n# # Solve Euler equation around a NACA0012 airfoil\n\nusing Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing StaticArrays\nusing BenchmarkTools\nusing Roots\nusing SparseArrays\nusing Profile\nusing InteractiveUtils\nusing WriteVTK\nusing DifferentialEquations\nusing Symbolics\nusing SparseDiffTools\n\nconst dir = string(@__DIR__, \"/\")\n\nfunction compute_residual(qdof, Q, V, params)\n q = (FEFunction(Q, qdof)...,)\n\n # alias on measures\n dΓ = params.dΓ\n dΩ = params.dΩ\n dΓ_wall = params.dΓ_wall\n dΓ_farfield = params.dΓ_farfield\n\n # Allocate rhs vectors\n b_vol = zero(qdof)\n b_fac = zero(qdof)\n\n # compute volume residuals\n l_vol(v) = ∫(flux_Ω(q, v))dΩ\n assemble_linear!(b_vol, l_vol, V)\n\n # face normals for each face domain (lazy, no computation at this step)\n nΓ = get_face_normals(dΓ)\n nΓ_wall = get_face_normals(dΓ_wall)\n nΓ_farfield = get_face_normals(dΓ_farfield)\n\n # flux residuals from interior faces for all variables\n l_Γ(v) = ∫(flux_Γ(q, v, nΓ))dΓ\n assemble_linear!(b_fac, l_Γ, V)\n\n # flux residuals from bc faces for all variables\n l_Γ_wall(v) = ∫(flux_Γ_wall(q, v, nΓ_wall))dΓ_wall\n l_Γ_farfield(v) = ∫(flux_Γ_farfield(q, v, nΓ_farfield))dΓ_farfield\n assemble_linear!(b_fac, l_Γ_wall, V)\n assemble_linear!(b_fac, l_Γ_farfield, V)\n dQ = b_vol .- b_fac\n\n return dQ\nend\n\n\"\"\"\n flux_Ω(q, v)\n\nCompute volume residual using the lazy-operators approach\n\"\"\"\nflux_Ω(q, v) = _flux_Ω ∘ (q, map(∇, v))\n\nfunction _flux_Ω(q, ∇v)\n ρ, ρu, ρE = q\n ∇λ_ρ, ∇λ_ρu, ∇λ_ρE = ∇v\n γ = stateInit.γ\n\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = pressure(ρ, ρu, ρE, γ)\n\n flux_ρ = ρu\n flux_ρu = ρuu + p * I\n flux_ρE = (ρE + p) .* vel\n\n return return ∇λ_ρ ⋅ flux_ρ + ∇λ_ρu ⊡ flux_ρu + ∇λ_ρE ⋅ flux_ρE\nend\n\n\"\"\"\n flux_Γ(q, v, n)\n\nFlux at the interface is defined by a composition of two functions:\n* the input states at face sides which are needed for the riemann flux\n* `flux_roe` defines the Riemann flux (as usual)\n\"\"\"\nflux_Γ(q, v, n) = flux_roe ∘ (side⁻(q), side⁺(q), jump(v), side⁻(n))\n\n\"\"\"\n flux_roe(q⁻, q⁺, δv, n)\n\"\"\"\nfunction flux_roe(q⁻, q⁺, δv, n)\n γ = stateInit.γ\n nx, ny = n\n ρ1, (ρu1, ρv1), ρE1 = q⁻\n ρ2, (ρu2, ρv2), ρE2 = q⁺\n δλ_ρ1, δλ_ρu1, δλ_ρE1 = δv\n\n ρ1 = max(eps(ρ1), ρ1)\n ρ2 = max(eps(ρ2), ρ2)\n\n # Closure\n u1 = ρu1 / ρ1\n v1 = ρv1 / ρ1\n u2 = ρu2 / ρ2\n v2 = ρv2 / ρ2\n p1 = pressure(ρ1, SA[ρu1, ρv1], ρE1, γ)\n p2 = pressure(ρ2, SA[ρu2, ρv2], ρE2, γ)\n\n H2 = (γ / (γ - 1)) * p2 / ρ2 + (u2 * u2 + v2 * v2) / 2.0\n H1 = (γ / (γ - 1)) * p1 / ρ1 + (u1 * u1 + v1 * v1) / 2.0\n\n R = √(ρ1 / ρ2)\n invR1 = 1.0 / (R + 1)\n uAv = (R * u1 + u2) * invR1\n vAv = (R * v1 + v2) * invR1\n Hav = (R * H1 + H2) * invR1\n cAv = √(abs((γ - 1) * (Hav - (uAv * uAv + vAv * vAv) / 2.0)))\n ecAv = (uAv * uAv + vAv * vAv) / 2.0\n\n λ1 = nx * uAv + ny * vAv\n λ3 = λ1 + cAv\n λ4 = λ1 - cAv\n\n d1 = ρ1 - ρ2\n d2 = ρ1 * u1 - ρ2 * u2\n d3 = ρ1 * v1 - ρ2 * v2\n d4 = ρE1 - ρE2\n\n # computation of the centered part of the flux\n flux_ρ = nx * ρ2 * u2 + ny * ρ2 * v2\n flux_ρu = nx * p2 + flux_ρ * u2\n flux_ρv = ny * p2 + flux_ρ * v2\n flux_ρE = H2 * flux_ρ\n\n # Temp variables\n rc1 = (γ - 1) / cAv\n rc2 = (γ - 1) / cAv / cAv\n uq41 = ecAv / cAv + cAv / (γ - 1)\n uq42 = nx * uAv + ny * vAv\n\n fdc1 = max(λ1, 0.0) * (d1 + rc2 * (-ecAv * d1 + uAv * d2 + vAv * d3 - d4))\n fdc2 = max(λ1, 0.0) * ((nx * vAv - ny * uAv) * d1 + ny * d2 - nx * d3)\n fdc3 =\n max(λ3, 0.0) * (\n (-uq42 * d1 + nx * d2 + ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n fdc4 =\n max(λ4, 0.0) * (\n (uq42 * d1 - nx * d2 - ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n\n duv1 = fdc1 + (fdc3 + fdc4) / cAv\n duv2 = uAv * fdc1 + ny * fdc2 + (uAv / cAv + nx) * fdc3 + (uAv / cAv - nx) * fdc4\n duv3 = vAv * fdc1 - nx * fdc2 + (vAv / cAv + ny) * fdc3 + (vAv / cAv - ny) * fdc4\n duv4 =\n ecAv * fdc1 +\n (ny * uAv - nx * vAv) * fdc2 +\n (uq41 + uq42) * fdc3 +\n (uq41 - uq42) * fdc4\n\n flux_ρ += duv1\n flux_ρu += duv2\n flux_ρv += duv3\n flux_ρE += duv4\n\n return (δλ_ρ1 ⋅ flux_ρ + δλ_ρu1 ⋅ SA[flux_ρu, flux_ρv] + δλ_ρE1 ⋅ flux_ρE)\nend\n\n\"\"\"\n flux_Γ_farfield(q, v, n)\n\nCompute `Roe` flux on boundary face by imposing\n`stateBcFarfield.u_in` on `side_p`\n\"\"\"\nflux_Γ_farfield(q, v, n) = flux_roe ∘ (side⁻(q), stateBcFarfield.u_inf, side⁻(v), side⁻(n))\n\n\"\"\"\n flux_Γ_wall(q, v, n)\n\"\"\"\nflux_Γ_wall(q, v, n) = _flux_Γ_wall ∘ (side⁻(q), side⁻(v), side⁻(n))\n\nfunction _flux_Γ_wall(q⁻, v⁻, n)\n γ = stateInit.γ\n ρ1, ρu1, ρE1 = q⁻\n λ_ρ1, λ_ρu1, λ_ρE1 = v⁻\n\n p1 = pressure(ρ1, ρu1, ρE1, γ)\n\n flux_ρ = zero(ρ1)\n flux_ρu = p1 * n\n flux_ρE = zero(ρE1)\n\n return (λ_ρ1 ⋅ flux_ρ + λ_ρu1 ⋅ flux_ρu + λ_ρE1 ⋅ flux_ρE)\nend\n\nfunction sparse2vtk(\n a::AbstractSparseMatrix,\n name::String = string(@__DIR__, \"/../myout/sparse\"),\n)\n vtk_write_array(name, Array(a), \"my_property_name\")\nend\n\nmutable struct VtkHandler\n basename::String\n basename_residual::String\n ite::Int\n VtkHandler(basename) = new(basename, basename * \"_residual\", 0)\nend\n\n\"\"\"\n Write solution (at cell centers) to vtk\n Wrapper for `write_vtk`\n\"\"\"\nfunction append_vtk(vtk, mesh, vars, t, params; res = nothing)\n ρ, ρu, ρE = vars\n\n # Mean cell values\n # name2val_mean = (;zip(get_name.(vars), mean_values.(vars, degquad))...)\n # p_mean = pressure.(name2val_mean[:ρ], name2val_mean[:ρu], name2val_mean[:ρE], params.stateInit.γ)\n\n vtk_degree = maximum(x -> get_degree(Bcube.get_function_space(get_fespace(x))), vars)\n vtk_degree = max(1, mesh_degree, vtk_degree)\n _ρ = var_on_nodes_discontinuous(ρ, mesh, vtk_degree)\n _ρu = var_on_nodes_discontinuous(ρu, mesh, vtk_degree)\n _ρE = var_on_nodes_discontinuous(ρE, mesh, vtk_degree)\n\n Cp = pressure_coefficient.(_ρ, _ρu, _ρE)\n Ma = mach.(_ρ, _ρu, _ρE)\n dict_vars_dg = Dict(\n \"rho\" => (_ρ, VTKPointData()),\n \"rhou\" => (_ρu, VTKPointData()),\n \"rhoE\" => (_ρE, VTKPointData()),\n \"Cp\" => (Cp, VTKPointData()),\n \"Mach\" => (Ma, VTKPointData()),\n \"rho_mean\" => (get_values(Bcube.cell_mean(ρ, params.dΩ)), VTKCellData()),\n \"rhou_mean\" => (get_values(Bcube.cell_mean(ρu, params.dΩ)), VTKCellData()),\n \"rhoE_mean\" => (get_values(Bcube.cell_mean(ρE, params.dΩ)), VTKCellData()),\n \"lim_rho\" => (get_values(params.limρ), VTKCellData()),\n \"lim_all\" => (get_values(params.limAll), VTKCellData()),\n )\n Bcube.write_vtk_discontinuous(\n vtk.basename * \"_DG\",\n vtk.ite,\n t,\n mesh,\n dict_vars_dg,\n vtk_degree;\n append = vtk.ite > 0,\n )\n\n _ρ_wall = var_on_bnd_nodes_discontinuous(ρ, params.Γ_wall, vtk_degree)\n _ρu_wall = var_on_bnd_nodes_discontinuous(ρu, params.Γ_wall, vtk_degree)\n _ρE_wall = var_on_bnd_nodes_discontinuous(ρE, params.Γ_wall, vtk_degree)\n\n Cp_wall = pressure_coefficient.(_ρ_wall, _ρu_wall, _ρE_wall)\n Ma_wall = pressure_coefficient.(_ρ_wall, _ρu_wall, _ρE_wall)\n\n dict_vars_wall = Dict(\n \"rho\" => (_ρ_wall, VTKPointData()),\n \"rhou\" => (_ρu_wall, VTKPointData()),\n \"rhoE\" => (_ρE_wall, VTKPointData()),\n \"Cp\" => (Cp_wall, VTKPointData()),\n \"Mach\" => (Ma_wall, VTKPointData()),\n )\n Bcube.write_vtk_bnd_discontinuous(\n vtk.basename * \"_bnd_DG\",\n 1,\n 0.0,\n params.Γ_wall,\n dict_vars_wall,\n vtk_degree;\n append = false,\n )\n\n #residual:\n if !isa(res, Nothing)\n vtkfile = vtk_grid(vtk.basename_residual, Float64.(res.iter), [0.0, 1.0])\n for (k, valₖ) in enumerate(res.val)\n vtkfile[\"res_\" * string(k), VTKPointData()] = [valₖ valₖ]\n end\n vtk_save(vtkfile)\n end\n\n # Update counter\n vtk.ite += 1\n\n return nothing\nend\n\nfunction init!(q, dΩ, initstate)\n AoA = initstate.AoA\n Minf = initstate.M_inf\n Pinf = initstate.P_inf\n Tinf = initstate.T_inf\n r = initstate.r_gas\n γ = initstate.γ\n\n ρinf = Pinf / r / Tinf\n ainf = √(γ * r * Tinf)\n Vinf = Minf * ainf\n ρVxinf = ρinf * Vinf * cos(AoA)\n ρVyinf = ρinf * Vinf * sin(AoA)\n ρEinf = Pinf / (γ - 1) + 0.5 * ρinf * Vinf^2\n\n ρ0 = PhysicalFunction(x -> ρinf)\n ρu0 = PhysicalFunction(x -> SA[ρVxinf, ρVyinf])\n ρE0 = PhysicalFunction(x -> ρEinf)\n projection_l2!(q, (ρ0, ρu0, ρE0), dΩ)\n return nothing\nend\n\nfunction main(stateInit, stateBcFarfield, degree)\n @show degree, degquad\n\n mesh = read_msh(dir * \"../input/mesh/naca0012_o\" * string(mesh_degree) * \".msh\", 2)\n scale!(mesh, 1.0 / 0.5334)\n\n dimcar = compute_dimcar(mesh)\n\n DMPrelax = DMPcurv₀ .* dimcar .^ 2\n\n # Then we create a `NamedTuple` to hold the simulation parameters.\n params = (\n degquad = degquad,\n stateInit = stateInit,\n stateBcFarfield = stateBcFarfield,\n DMPrelax = DMPrelax,\n )\n\n # Define measures for cell and interior face integrations\n dΩ = Measure(CellDomain(mesh), degquad)\n dΓ = Measure(InteriorFaceDomain(mesh), degquad)\n\n # Declare boundary conditions and\n # create associated domains and measures\n Γ_wall = BoundaryFaceDomain(mesh, (\"NACA\",))\n Γ_farfield = BoundaryFaceDomain(mesh, (\"FARFIELD\",))\n dΓ_wall = Measure(Γ_wall, degquad)\n dΓ_farfield = Measure(Γ_farfield, degquad)\n\n params = (\n params...,\n Γ_wall = Γ_wall,\n dΓ = dΓ,\n dΩ = dΩ,\n dΓ_wall = dΓ_wall,\n dΓ_farfield = dΓ_farfield,\n )\n\n qLowOrder = nothing\n\n for deg in 0:degree\n params = (params..., degree = deg)\n\n fs = FunctionSpace(fspace, deg)\n Q_sca = TrialFESpace(fs, mesh, :discontinuous; size = 1) # DG, scalar\n Q_vec = TrialFESpace(fs, mesh, :discontinuous; size = 2) # DG, vectoriel\n V_sca = TestFESpace(Q_sca)\n V_vec = TestFESpace(Q_vec)\n Q = MultiFESpace(Q_sca, Q_vec, Q_sca)\n V = MultiFESpace(V_sca, V_vec, V_sca)\n\n q = FEFunction(Q)\n\n # select an initial configurations:\n if deg == 0\n init!(q, mesh, stateInit)\n else\n println(\"Start projection\")\n projection_l2!(q, qLowOrder, dΩ)\n println(\"End projection\")\n end\n\n # create CellData to store limiter values\n limρ = Bcube.MeshCellData(ones(ncells(mesh)))\n limAll = Bcube.MeshCellData(ones(ncells(mesh)))\n params = (params..., limρ = limρ, limAll = limAll)\n\n # Init vtk handler\n mkpath(outputpath)\n vtk = VtkHandler(\n outputpath * \"euler_naca_mdeg\" * string(mesh_degree) * \"_deg\" * string(deg),\n )\n\n # Init time\n time = 0.0\n\n # Save initial solution\n append_vtk(vtk, mesh, q, time, params)\n\n # Build the cache and store everything you want to compute only once (such as the mass matrice inverse...)\n\n cache = ()\n # Allocate buffer for compute_residual\n b_vol = zeros(Bcube.get_ndofs(Q))\n b_fac = zeros(Bcube.get_ndofs(Q))\n cache = (cache..., b_vol = b_vol, b_fac = b_fac)\n\n cache = (\n cache...,\n cacheCellMean = Bcube.build_cell_mean_cache(q, dΩ),\n mass = factorize(Bcube.build_mass_matrix(Q, V, dΩ)),\n mass_sca = factorize(Bcube.build_mass_matrix(Q_sca, V_sca, dΩ)),\n mass_vec = factorize(Bcube.build_mass_matrix(Q_vec, V_vec, dΩ)),\n )\n\n time, q = steady_solve!(Q, V, q, mesh, params, cache, vtk, deg)\n append_vtk(vtk, mesh, q, time, params)\n println(\"end steady_solve for deg=\", deg, \" !\")\n\n deg < degree && (qLowOrder = deepcopy(q))\n end\n return nothing\nend\n\nfunction steady_solve!(Q, V, q, mesh, params, cache, vtk, deg)\n counter = [0]\n q0 = deepcopy(get_dof_values(q))\n ode_params =\n (Q = Q, V = V, params = params, cache = cache, counter = counter, vtk = vtk)\n\n rhs!(dq, q, p, t) = dq .= compute_residual(q, p.Q, p.V, p.params)\n\n # compute sparsity pattern and coloring\n println(\"computing jacobian cache...\")\n if withbench\n _rhs!(dq, q) = rhs!(dq, q, ode_params, 0.0)\n @btime $_rhs!(similar($q0), $q0)\n q_bench = FEFunction(Q, q0)\n @btime $apply_limitation!($q_bench, $ode_params)\n @show length(q0)\n end\n\n #sparsity_pattern = Symbolics.jacobian_sparsity(_rhs!, similar(Q0), Q0)\n #tjac = @elapsed Symbolics.jacobian_sparsity(_rhs!, similar(Q0), Q0)\n #@show tjac\n sparsity_pattern = Bcube.build_jacobian_sparsity_pattern(Q, mesh)\n println(\"sparsity pattern computed !\")\n display(sparsity_pattern)\n colors = matrix_colors(sparsity_pattern)\n println(\"coloring done!\")\n @show maximum(colors)\n\n ode = ODEFunction(\n rhs!;\n mass_matrix = Bcube.build_mass_matrix(Q, V, params.dΩ),\n jac_prototype = sparsity_pattern,\n colorvec = colors,\n )\n\n Tfinal = Inf\n problem = ODEProblem(ode, q0, (0.0, Tfinal), ode_params)\n timestepper = ImplicitEuler(; nlsolve = NLNewton(; max_iter = 20))\n\n cb_cache = DiscreteCallback(always_true, update_cache!; save_positions = (false, false))\n cb_vtk = DiscreteCallback(always_true, output_vtk; save_positions = (false, false))\n cb_steady = TerminateSteadyState(1e-6, 1e-6, condition_steadystate)\n\n error = 1e-1\n\n sol = solve(\n problem,\n timestepper;\n initializealg = NoInit(),\n adaptive = true,\n abstol = error,\n reltol = error,\n progress = false,\n progress_steps = 1000,\n save_everystep = false,\n save_start = false,\n save_end = false,\n isoutofdomain = isoutofdomain,\n callback = CallbackSet(cb_cache, cb_vtk, cb_steady),\n )\n\n set_dof_values!(q, sol.u[end])\n return sol.t[end], q\nend\n\nalways_true(args...) = true\n\nfunction isoutofdomain(dof, p, t)\n any(isnan, dof) && return true\n\n q = FEFunction(p.Q, dof)\n q_mean = map(get_values, Bcube.cell_mean(q, p.cache.cacheCellMean))\n p_mean = pressure.(q_mean..., stateInit.γ)\n\n negative_ρ = any(x -> x < 0, q_mean[1])\n negative_p = any(x -> x < 0, p_mean)\n isout = negative_ρ || negative_p\n isout && @show negative_ρ, negative_p\n return isout\nend\n\nfunction update_cache!(integrator)\n Q = integrator.p.Q\n Q1, = Q\n deg = get_degree(Bcube.get_function_space(Q1))\n println(\n \"deg=\",\n deg,\n \" update_cache! : iter=\",\n integrator.p.counter[1],\n \" dt=\",\n integrator.dt,\n )\n\n q = FEFunction(integrator.p.Q, integrator.u)\n limiter_projection && apply_limitation!(q, integrator.p)\n return nothing\nend\n\nfunction output_vtk(integrator)\n u_modified!(integrator, false)\n mesh = get_mesh(get_domain(integrator.p.params.dΩ))\n q = FEFunction(integrator.p.Q, integrator.u)\n counter = integrator.p.counter\n counter .+= 1\n if (counter[1] % nout == 0)\n println(\"output_vtk \", counter[1])\n append_vtk(integrator.p.vtk, mesh, q, integrator.t, integrator.p.params)\n end\n return nothing\nend\n\nfunction condition_steadystate(integrator, abstol, reltol, min_t)\n u_modified!(integrator, false)\n if DiffEqBase.isinplace(integrator.sol.prob)\n testval = first(get_tmp_cache(integrator))\n @. testval = (integrator.u - integrator.uprev) / (integrator.t - integrator.tprev)\n else\n testval = (integrator.u - integrator.uprev) / (integrator.t - integrator.tprev)\n end\n\n if typeof(integrator.u) <: Array\n any(\n abs(d) > abstol && abs(d) > reltol * abs(u) for (d, abstol, reltol, u) in\n zip(testval, Iterators.cycle(abstol), Iterators.cycle(reltol), integrator.u)\n ) && (return false)\n else\n any((abs.(testval) .> abstol) .& (abs.(testval) .> reltol .* abs.(integrator.u))) &&\n (return false)\n end\n\n if min_t === nothing\n return true\n else\n return integrator.t >= min_t\n end\nend\n\n\"\"\"\nCompute the characteristic dimension of each cell of `mesh`:\ndimcar = (cell volume) / (cell surface)\n\n# TODO :\nto be moved to Bcube\n\"\"\"\nfunction compute_dimcar(mesh)\n fs = FunctionSpace(:Lagrange, 0)\n V = TestFESpace(fs, mesh; size = 1, isContinuous = false)\n\n # Define measures for cell and interior face integrations\n dΩ = Measure(CellDomain(mesh), degquad)\n dΓ = Measure(InteriorFaceDomain(mesh), degquad)\n dΓ_bc = Measure(BoundaryFaceDomain(mesh), degquad)\n\n f1 = PhysicalFunction(x -> 1.0)\n l(v) = ∫(f1 ⋅ v)dΩ\n l_face(v, dω) = ∫(side⁻(f1) ⋅ side⁻(v) + side⁺(f1) ⋅ side⁺(v))dω\n\n vol = assemble_linear(l, V)\n surf = assemble_linear(Base.Fix2(l_face, dΓ), V)\n surf += assemble_linear(Base.Fix2(l_face, dΓ_bc), V)\n return vol ./ surf\nend\n\n\"\"\"\nReferences:\n* Xiangxiong Zhang, Chi-Wang Shu, On positivity-preserving high order discontinuous\n Galerkin schemes for compressible Euler equations on rectangular meshes,\n Journal of Computational Physics, Volume 229, Issue 23, 2010.\n https://doi.org/10.1016/j.jcp.2010.08.016\n* Zhang, X., Xia, Y. & Shu, CW. Maximum-Principle-Satisfying and Positivity-Preserving\n High Order Discontinuous Galerkin Schemes for Conservation Laws on Triangular Meshes.\n J Sci Comput 50, 29–62 (2012). https://doi.org/10.1007/s10915-011-9472-8\n\"\"\"\nfunction apply_limitation!(q::Bcube.AbstractFEFunction, ode_params)\n params = ode_params.params\n cache = ode_params.cache\n mesh = get_mesh(get_domain(params.dΩ))\n ρ, ρu, ρE = q\n\n ρ_mean, ρu_mean, ρE_mean = Bcube.cell_mean(q, cache.cacheCellMean)\n\n _limρ, ρ_proj = linear_scaling_limiter(\n ρ,\n params.dΩ;\n bounds = (ρmin₀, ρmax₀),\n DMPrelax = params.DMPrelax,\n mass = cache.mass_sca,\n )\n\n op_t = limiter_param_p ∘ (ρ_proj, ρu, ρE, ρ_mean, ρu_mean, ρE_mean)\n t = Bcube._minmax_cells(op_t, mesh, Val(params.degquad))\n tmin = Bcube.MeshCellData(getindex.(t, 1))\n\n if eltype(_limρ) == eltype(params.limρ) # skip Dual number case\n set_values!(params.limρ, get_values(_limρ))\n set_values!(params.limAll, get_values(tmin))\n end\n\n limited_var(u, ū, lim_u) = ū + lim_u * (u - ū)\n projection_l2!(ρ, limited_var(ρ_proj, ρ_mean, tmin), params.dΩ; mass = cache.mass_sca)\n projection_l2!(ρu, limited_var(ρu, ρu_mean, tmin), params.dΩ; mass = cache.mass_vec)\n projection_l2!(ρE, limited_var(ρE, ρE_mean, tmin), params.dΩ; mass = cache.mass_sca)\n return nothing\nend\n\nfunction limiter_param_p(ρ̂, ρu, ρE, ρ_mean, ρu_mean, ρE_mean)\n γ = stateInit.γ\n p = pressure(ρ̂, ρu, ρE, γ)\n\n if p ≥ pmin₀\n t = 1.0\n else\n @show p, ρ̂, ρu, ρE\n @show ρ_mean, ρu_mean, ρE_mean\n @show pressure(ρ_mean, ρu_mean, ρE_mean, γ)\n if pressure(ρ_mean, ρu_mean, ρE_mean, γ) > pmin₀\n fₜ =\n t ->\n pressure(\n t * ρ̂ + (1 - t) * ρ_mean,\n t * ρu + (1 - t) * ρu_mean,\n t * ρE + (1 - t) * ρE_mean,\n γ,\n ) - pmin₀\n bounds = (0.0, 1.0)\n t = find_zero(fₜ, bounds, Bisection())\n else\n t = NaN\n println(\"t = NaN\")\n end\n end\n\n return t\nend\n\nfunction pressure(ρ::Number, ρu::AbstractVector, ρE::Number, γ)\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = (γ - 1) * (ρE - tr(ρuu) / 2)\n return p\nend\n\ncompute_Pᵢ(P, γ, M) = P * (1 + 0.5 * (γ - 1) * M^2)^(γ / (γ - 1))\ncompute_Tᵢ(T, γ, M) = T * (1 + 0.5 * (γ - 1) * M^2)\n\nfunction bc_state_farfield(AoA, M, P, T, r, γ)\n a = √(γ * r * T)\n vn = M * a\n ρ = P / r / T\n ρu = SA[ρ * vn * cos(AoA), ρ * vn * sin(AoA)]\n ρE = P / (γ - 1) + 0.5 * ρ * vn^2\n return (ρ, ρu, ρE)\nend\n\nfunction pressure_coefficient(ρ, ρu, ρE)\n (pressure(ρ, ρu, ρE, stateInit.γ) - stateInit.P_inf) /\n (stateBcFarfield.Pᵢ_inf - stateInit.P_inf)\nend\nfunction mach(ρ, ρu, ρE)\n norm(ρu ./ ρ) / √(stateInit.γ * max(0.0, pressure(ρ, ρu, ρE, stateInit.γ) / ρ))\nend\n\nconst degreemax = 2 # Function-space degree\nconst mesh_degree = 2\nconst fspace = :Lagrange\nconst limiter_projection = true\nconst ρmin₀ = 1.0e-8\nconst ρmax₀ = 1.0e+10\nconst pmin₀ = 1.0e-8\nconst pmax₀ = 1.0e+10\nconst DMPcurv₀ = 10.0e3\nconst withbench = false\n\nconst stateInit = (\n AoA = deg2rad(1.25),\n M_inf = 0.8,\n P_inf = 101325.0,\n T_inf = 275.0,\n r_gas = 287.0,\n γ = 1.4,\n)\nconst nite_max = 300 #300000 # Number of time iteration(s)\nconst nout = 1 # number of step between two vtk outputs\nconst mass_matrix_in_solve = true\nconst degquad = 6\nconst outputpath = string(@__DIR__, \"/../myout/euler_naca_steady/\")\n\nconst stateBcFarfield = (\n AoA = stateInit.AoA,\n M_inf = stateInit.M_inf,\n Pᵢ_inf = compute_Pᵢ(stateInit.P_inf, stateInit.γ, stateInit.M_inf),\n Tᵢ_inf = compute_Tᵢ(stateInit.T_inf, stateInit.γ, stateInit.M_inf),\n u_inf = bc_state_farfield(\n stateInit.AoA,\n stateInit.M_inf,\n stateInit.P_inf,\n stateInit.T_inf,\n stateInit.r_gas,\n stateInit.γ,\n ),\n r_gas = stateInit.r_gas,\n γ = stateInit.γ,\n)\n\nmain(stateInit, stateBcFarfield, degreemax)\n\nend #hide","category":"page"},{"location":"example/covo/#Euler-equations-covo","page":"Euler equations - covo","title":"Euler equations - covo","text":"","category":"section"},{"location":"example/covo/","page":"Euler equations - covo","title":"Euler equations - covo","text":"module Covo #hide\nprintln(\"Running covo example...\") #hide\n\nconst dir = string(@__DIR__, \"/\")\nusing Bcube\nusing LinearAlgebra\nusing StaticArrays\nusing WriteVTK # only for 'VTKCellData'\nusing Profile\nusing StaticArrays\nusing InteractiveUtils\nusing BenchmarkTools\nusing UnPack\n\nfunction compute_residual(_u, V, params, cache)\n u = get_fe_functions(_u)\n\n # alias on measures\n @unpack dΩ, dΓ, dΓ_perio_x, dΓ_perio_y = params\n\n # face normals for each face domain (lazy, no computation at this step)\n nΓ = get_face_normals(dΓ)\n nΓ_perio_x = get_face_normals(dΓ_perio_x)\n nΓ_perio_y = get_face_normals(dΓ_perio_y)\n\n # flux residuals from faces for all variables\n function l(v)\n ∫(flux_Ω(u, v))dΩ +\n -∫(flux_Γ(u, v, nΓ))dΓ +\n -∫(flux_Γ(u, v, nΓ_perio_x))dΓ_perio_x +\n -∫(flux_Γ(u, v, nΓ_perio_y))dΓ_perio_y\n end\n\n rhs = assemble_linear(l, V)\n\n return cache.mass \\ rhs\nend\n\n\"\"\"\n flux_Ω(u, v)\n\nCompute volume residual using the lazy-operators approach\n\"\"\"\nflux_Ω(u, v) = _flux_Ω ∘ cellvar(u, v)\ncellvar(u, v) = (u, map(∇, v))\nfunction _flux_Ω(u, ∇v)\n ρ, ρu, ρE, ϕ = u\n ∇λ_ρ, ∇λ_ρu, ∇λ_ρE, ∇λ_ϕ = ∇v\n\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = pressure(ρ, ρu, ρE, γ)\n\n flux_ρ = ρu\n flux_ρu = ρuu + p * I\n flux_ρE = (ρE + p) .* vel\n flux_ϕ = ϕ .* vel\n\n return ∇λ_ρ ⋅ flux_ρ + ∇λ_ρu ⊡ flux_ρu + ∇λ_ρE ⋅ flux_ρE + ∇λ_ϕ ⋅ flux_ϕ\nend\n\n\"\"\"\n flux_Γ(u,v,n)\n\nFlux at the interface is defined by a composition of two functions:\n* facevar(u,v,n) defines the input states which are needed for\n the riemann flux using operator notations\n* flux_roe(w) defines the Riemann flux (as usual)\n\"\"\"\nflux_Γ(u, v, n) = flux_roe ∘ (side⁻(u), side⁺(u), jump(v), side⁻(n))\n\n\"\"\"\n flux_roe(w)\n\"\"\"\nfunction flux_roe(ui, uj, δv, nij)\n # destructuring inputs given by `facevar` function\n\n nx, ny = nij\n ρ1, ρu1, ρE1, ϕ1 = ui\n ρ2, ρu2, ρE2, ϕ2 = uj\n δλ_ρ1, δλ_ρu1, δλ_ρE1, δλ_ϕ1 = δv\n ρux1, ρuy1 = ρu1\n ρux2, ρuy2 = ρu2\n\n # Closure\n u1 = ρux1 / ρ1\n v1 = ρuy1 / ρ1\n u2 = ρux2 / ρ2\n v2 = ρuy2 / ρ2\n p1 = pressure(ρ1, ρu1, ρE1, γ)\n p2 = pressure(ρ2, ρu2, ρE2, γ)\n\n H2 = (γ / (γ - 1)) * p2 / ρ2 + (u2 * u2 + v2 * v2) / 2.0\n H1 = (γ / (γ - 1)) * p1 / ρ1 + (u1 * u1 + v1 * v1) / 2.0\n\n R = √(ρ1 / ρ2)\n invR1 = 1.0 / (R + 1)\n uAv = (R * u1 + u2) * invR1\n vAv = (R * v1 + v2) * invR1\n Hav = (R * H1 + H2) * invR1\n cAv = √(abs((γ - 1) * (Hav - (uAv * uAv + vAv * vAv) / 2.0)))\n ecAv = (uAv * uAv + vAv * vAv) / 2.0\n\n λ1 = nx * uAv + ny * vAv\n λ3 = λ1 + cAv\n λ4 = λ1 - cAv\n\n d1 = ρ1 - ρ2\n d2 = ρ1 * u1 - ρ2 * u2\n d3 = ρ1 * v1 - ρ2 * v2\n d4 = ρE1 - ρE2\n\n # computation of the centered part of the flux\n flu11 = nx * ρ2 * u2 + ny * ρ2 * v2\n flu21 = nx * p2 + flu11 * u2\n flu31 = ny * p2 + flu11 * v2\n flu41 = H2 * flu11\n\n # Temp variables\n rc1 = (γ - 1) / cAv\n rc2 = (γ - 1) / cAv / cAv\n uq41 = ecAv / cAv + cAv / (γ - 1)\n uq42 = nx * uAv + ny * vAv\n\n fdc1 = max(λ1, 0.0) * (d1 + rc2 * (-ecAv * d1 + uAv * d2 + vAv * d3 - d4))\n fdc2 = max(λ1, 0.0) * ((nx * vAv - ny * uAv) * d1 + ny * d2 - nx * d3)\n fdc3 =\n max(λ3, 0.0) * (\n (-uq42 * d1 + nx * d2 + ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n fdc4 =\n max(λ4, 0.0) * (\n (uq42 * d1 - nx * d2 - ny * d3) / 2.0 +\n rc1 * (ecAv * d1 - uAv * d2 - vAv * d3 + d4) / 2.0\n )\n\n duv1 = fdc1 + (fdc3 + fdc4) / cAv\n duv2 = uAv * fdc1 + ny * fdc2 + (uAv / cAv + nx) * fdc3 + (uAv / cAv - nx) * fdc4\n duv3 = vAv * fdc1 - nx * fdc2 + (vAv / cAv + ny) * fdc3 + (vAv / cAv - ny) * fdc4\n duv4 =\n ecAv * fdc1 +\n (ny * uAv - nx * vAv) * fdc2 +\n (uq41 + uq42) * fdc3 +\n (uq41 - uq42) * fdc4\n\n v₁₂ = 0.5 * ((u1 + u2) * nx + (v1 + v2) * ny)\n fluxϕ = max(0.0, v₁₂) * ϕ1 + min(0.0, v₁₂) * ϕ2\n\n return (\n δλ_ρ1 * (flu11 + duv1) +\n δλ_ρu1 ⋅ (SA[flu21 + duv2, flu31 + duv3]) +\n δλ_ρE1 * (flu41 + duv4) +\n δλ_ϕ1 * (fluxϕ)\n )\nend\n\n\"\"\"\nTime integration of `f(q, t)` over a timestep `Δt`.\n\"\"\"\nforward_euler(q, f::Function, t, Δt) = get_dof_values(q) .+ Δt .* f(q, t)\n\n\"\"\"\n rk3_ssp(q, f::Function, t, Δt)\n\n`f(q, t)` is the function to integrate.\n\"\"\"\nfunction rk3_ssp(q, f::Function, t, Δt)\n stepper(q, t) = forward_euler(q, f, t, Δt)\n _q0 = get_dof_values(q)\n\n _q1 = stepper(q, Δt)\n\n set_dof_values!(q, _q1)\n _q2 = (3 / 4) .* _q0 .+ (1 / 4) .* stepper(q, t + Δt)\n\n set_dof_values!(q, _q2)\n _q1 .= (1 / 3) * _q0 .+ (2 / 3) .* stepper(q, t + Δt / 2)\n\n return _q1\nend\n\n\"\"\"\n pressure(ρ, ρu, ρE, γ)\n\nComputes pressure from perfect gaz law\n\"\"\"\nfunction pressure(ρ::Number, ρu::AbstractVector, ρE::Number, γ)\n vel = ρu ./ ρ\n ρuu = ρu * transpose(vel)\n p = (γ - 1) * (ρE - tr(ρuu) / 2)\n return p\nend\n\n\"\"\"\n Init field with a vortex (for the COVO test case)\n\"\"\"\nfunction covo!(q, dΩ)\n\n # Intermediate vars\n Cₚ = γ * r / (γ - 1)\n\n r²(x) = ((x[1] .- xvc) .^ 2 + (x[2] .- yvc) .^ 2) ./ Rc^2\n # Temperature\n T(x) = T₀ .- β^2 * U₀^2 / (2 * Cₚ) .* exp.(-r²(x))\n # Velocity\n ux(x) = U₀ .- β * U₀ / Rc .* (x[2] .- yvc) .* exp.(-r²(x) ./ 2)\n uy(x) = V₀ .+ β * U₀ / Rc .* (x[1] .- xvc) .* exp.(-r²(x) ./ 2)\n # Density\n ρ(x) = ρ₀ .* (T(x) ./ T₀) .^ (1.0 / (γ - 1))\n # momentum\n ρu(x) = SA[ρ(x) * ux(x), ρ(x) * uy(x)]\n # Energy\n ρE(x) = ρ(x) * ((Cₚ / γ) .* T(x) + (ux(x) .^ 2 + uy(x) .^ 2) ./ 2)\n # Passive scalar\n ϕ(x) = Rc^2 * r²(x) < 0.01 ? exp(-r²(x) ./ 2) : 0.0\n\n f = map(PhysicalFunction, (ρ, ρu, ρE, ϕ))\n projection_l2!(q, f, dΩ)\n return nothing\nend\n\n\"\"\"\n Tiny struct to ease the VTK output\n\"\"\"\nmutable struct VtkHandler\n basename::Any\n ite::Any\n VtkHandler(basename) = new(basename, 0)\nend\n\n\"\"\"\n Write solution to vtk\n Wrapper for `write_vtk`\n\"\"\"\nfunction append_vtk(vtk, mesh, vars, t, params)\n ρ, ρu, ρE, ϕ = vars\n\n mesh_degree = 1\n vtk_degree = maximum(x -> get_degree(Bcube.get_function_space(get_fespace(x))), vars)\n vtk_degree = max(1, mesh_degree, vtk_degree)\n\n _ρ = var_on_nodes_discontinuous(ρ, mesh, vtk_degree)\n _ρu = var_on_nodes_discontinuous(ρu, mesh, vtk_degree)\n _ρE = var_on_nodes_discontinuous(ρE, mesh, vtk_degree)\n _ϕ = var_on_nodes_discontinuous(ϕ, mesh, vtk_degree)\n\n _p = pressure.(_ρ, _ρu, _ρE, γ)\n dict_vars_dg = Dict(\n \"rho\" => (_ρ, VTKPointData()),\n \"rhou\" => (_ρu, VTKPointData()),\n \"rhoE\" => (_ρE, VTKPointData()),\n \"phi\" => (_ϕ, VTKPointData()),\n \"p\" => (_p, VTKPointData()),\n )\n Bcube.write_vtk_discontinuous(\n vtk.basename * \"_DG\",\n vtk.ite,\n t,\n mesh,\n dict_vars_dg,\n vtk_degree;\n append = vtk.ite > 0,\n )\n\n # Update counter\n vtk.ite += 1\nend\n\n# Settings\nif get(ENV, \"BenchmarkMode\", \"false\") == \"false\" #hide\n const cellfactor = 1\n const nx = 32 * cellfactor + 1\n const ny = 32 * cellfactor + 1\n const fspace = :Lagrange\n const timeScheme = :ForwardEuler\nelse #hide\n const nx = 128 + 1 #hide\n const ny = 128 + 1 #hide\n const fspace = :Lagrange\n const timeScheme = :ForwardEuler\nend #hide\nconst nperiod = 1 # number of turn\nconst CFL = 0.1\nconst degree = 1 # FunctionSpace degree\nconst degquad = 2 * degree + 1\nconst γ = 1.4\nconst β = 0.2 # vortex intensity\nconst r = 287.15 # Perfect gaz constant\nconst T₀ = 300 # mean-flow temperature\nconst P₀ = 1e5 # mean-flow pressure\nconst M₀ = 0.5 # mean-flow mach number\nconst ρ₀ = 1.0 # mean-flow density\nconst xvc = 0.0 # x-center of vortex\nconst yvc = 0.0 # y-center of vortex\nconst Rc = 0.005 # Charasteristic vortex radius\nconst c₀ = √(γ * r * T₀) # Sound velocity\nconst U₀ = M₀ * c₀ # mean-flow velocity\nconst V₀ = 0.0 # mean-flow velocity\nconst ϕ₀ = 1.0\nconst l = 0.05 # half-width of the domain\nconst Δt = CFL * 2 * l / (nx - 1) / ((1 + β) * U₀ + c₀) / (2 * degree + 1)\n#const Δt = 5.e-7\nconst nout = 100 # Number of time steps to save\nconst outputpath = \"../myout/covo/\"\nconst output = joinpath(@__DIR__, outputpath, \"covo_deg$degree\")\nconst nite = Int(floor(nperiod * 2 * l / (U₀ * Δt))) + 1\n\nfunction run_covo()\n println(\"Starting run_covo...\")\n\n # Build mesh\n meshParam = (nx = nx, ny = ny, lx = 2l, ly = 2l, xc = 0.0, yc = 0.0)\n tmp_path = \"tmp.msh\"\n if get(ENV, \"BenchmarkMode\", \"false\") == \"false\" #hide\n gen_rectangle_mesh(tmp_path, :quad; meshParam...)\n else #hide\n if get(ENV, \"MeshConfig\", \"quad\") == \"triquad\" #hide\n gen_rectangle_mesh_with_tri_and_quad(tmp_path; meshParam...) #hide\n else #hide\n gen_rectangle_mesh(tmp_path, :quad; meshParam...) #hide\n end #hide\n end #hide\n mesh = read_msh(tmp_path)\n rm(tmp_path)\n\n # Define variables and test functions\n fs = FunctionSpace(fspace, degree)\n U_sca = TrialFESpace(fs, mesh, :discontinuous; size = 1) # DG, scalar\n U_vec = TrialFESpace(fs, mesh, :discontinuous; size = 2) # DG, vectoriel\n V_sca = TestFESpace(U_sca)\n V_vec = TestFESpace(U_vec)\n U = MultiFESpace(U_sca, U_vec, U_sca, U_sca)\n V = MultiFESpace(V_sca, V_vec, V_sca, V_sca)\n u = FEFunction(U)\n\n @show Bcube.get_ndofs(U)\n\n # Define measures for cell and interior face integrations\n dΩ = Measure(CellDomain(mesh), degquad)\n dΓ = Measure(InteriorFaceDomain(mesh), degquad)\n\n # Declare periodic boundary conditions and\n # create associated domains and measures\n periodicBCType_x = PeriodicBCType(Translation(SA[-2l, 0.0]), (\"East\",), (\"West\",))\n periodicBCType_y = PeriodicBCType(Translation(SA[0.0, 2l]), (\"South\",), (\"North\",))\n Γ_perio_x = BoundaryFaceDomain(mesh, periodicBCType_x)\n Γ_perio_y = BoundaryFaceDomain(mesh, periodicBCType_y)\n dΓ_perio_x = Measure(Γ_perio_x, degquad)\n dΓ_perio_y = Measure(Γ_perio_y, degquad)\n\n params = (dΩ = dΩ, dΓ = dΓ, dΓ_perio_x = dΓ_perio_x, dΓ_perio_y = dΓ_perio_y)\n\n # Init vtk\n isdir(joinpath(@__DIR__, outputpath)) || mkpath(joinpath(@__DIR__, outputpath))\n vtk = VtkHandler(output)\n\n # Init solution\n t = 0.0\n\n covo!(u, dΩ)\n\n # cache mass matrices\n cache = (mass = factorize(Bcube.build_mass_matrix(U, V, dΩ)),)\n\n if get(ENV, \"BenchmarkMode\", \"false\") == \"true\" #hide\n return u, U, V, params, cache\n end\n\n # Write initial solution\n append_vtk(vtk, mesh, u, t, params)\n\n # Time loop\n for i in 1:nite\n println(\"\")\n println(\"\")\n println(\"Iteration \", i, \" / \", nite)\n\n ## Step forward in time\n rhs(u, t) = compute_residual(u, V, params, cache)\n if timeScheme == :ForwardEuler\n unew = forward_euler(u, rhs, time, Δt)\n elseif timeScheme == :RK3\n unew = rk3_ssp(u, rhs, time, Δt)\n else\n error(\"Unknown time scheme: $timeScheme\")\n end\n\n set_dof_values!(u, unew)\n\n t += Δt\n\n # Write solution to file\n if (i % Int(max(floor(nite / nout), 1)) == 0)\n println(\"--> VTK export\")\n append_vtk(vtk, mesh, u, t, params)\n end\n end\n\n # Summary and benchmark # ndofs total = 20480\n _rhs(u, t) = compute_residual(u, V, params, cache)\n @btime forward_euler($u, $_rhs, $time, $Δt) # 5.639 ms (1574 allocations: 2.08 MiB)\n # stepper = w -> explicit_step(w, params, cache, Δt)\n # RK3_SSP(stepper, (u, v), cache)\n # @btime RK3_SSP($stepper, ($u, $v), $cache)\n println(\"ndofs total = \", Bcube.get_ndofs(U))\n Profile.init(; n = 10^7) # returns the current settings\n Profile.clear()\n Profile.clear_malloc_data()\n @profile begin\n for i in 1:100\n forward_euler(u, _rhs, time, Δt)\n end\n end\n @show Δt, U₀, U₀ * t\n @show boundary_names(mesh)\n return nothing\nend\n\nif get(ENV, \"BenchmarkMode\", \"false\") == \"false\"\n mkpath(outputpath)\n run_covo()\nend\n\nend #hide","category":"page"},{"location":"tutorial/helmholtz/#Helmholtz-equation-(FE)","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"","category":"section"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"This tutorial shows how to solve the Helmholtz eigen problem with a finite-element approach using Bcube.","category":"page"},{"location":"tutorial/helmholtz/#Theory","page":"Helmholtz equation (FE)","title":"Theory","text":"","category":"section"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"We consider the following Helmholtz equation, representing for instance the acoustic wave propagation with Neuman boundary condition(s):","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"begincases\n Delta u + omega^2 u = 0 \n dfracpartial upartial n = 0 textrm on Gamma\nendcases","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"An analytic solution of this equation can be obtained: for a rectangular domain Omega = 0L_x times 0L_y,","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"u(xy) = cos left( frack_x piL_x x right) cos left( frack_y piL_y y right) mathrmwith k_xk_y in mathbbN","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"with omega^2 = pi^2 left( dfrack_x^2L_x^2 + dfrack_y^2L_y^2 right)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Now, both the finite-element method and the discontinuous Galerkin method requires to write the weak form of the problem:","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"int_Omega (Delta u Delta v + omega^2u)v mathrmdOmega = 0","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"- int_Omega nabla u cdot nabla v mathrmdOmega\n+ underbraceleft (nabla u cdot n) v right_Gamma_=0 + omega^2 int_Omega u v mathrmd Omega = 0","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"int_Omega nabla u cdot nabla v mathrmd Omega = omega^2 int_Omega u v mathrmd Omega","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Introducing to bilinear forms a(uv) and b(uv) for respectively the left and right side terms, this equation can be written","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"a(uv) = omega^2 b(uv)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"This is actually a generalized eigenvalue problem which can be written:","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"A u = alpha B u","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"where","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"A u = int_Omega nabla u cdot nabla v mathrmd Omega B u = int_Omega u v mathrmd Omega alpha = omega^2","category":"page"},{"location":"tutorial/helmholtz/#Uncommented-code","page":"Helmholtz equation (FE)","title":"Uncommented code","text":"","category":"section"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"The code below solves the described Helmholtz eigen problem. The code with detailed comments is provided in the next section.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"using Bcube\nusing LinearAlgebra\n\nmesh = rectangle_mesh(21, 21)\n\ndegree = 1\n\nU = TrialFESpace(FunctionSpace(:Lagrange, degree), mesh)\nV = TestFESpace(U)\n\ndΩ = Measure(CellDomain(mesh), 2 * degree + 1)\n\na(u, v) = ∫(∇(u) ⋅ ∇(v))dΩ\nb(u, v) = ∫(u ⋅ v)dΩ\n\nA = assemble_bilinear(a, U, V)\nB = assemble_bilinear(b, U, V)\n\nvp, vecp = eigen(Matrix(A), Matrix(B))\n@show sqrt.(abs.(vp[3:8]))","category":"page"},{"location":"tutorial/helmholtz/#Commented-code","page":"Helmholtz equation (FE)","title":"Commented code","text":"","category":"section"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Load the necessary packages","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"using Bcube\nusing LinearAlgebra","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Mesh a 2D rectangular domain with quads.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"mesh = rectangle_mesh(21, 21)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Next, create the function space that will be used for the trial and test spaces. The Lagrange polynomial space is used here. The degree is set to 1.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"degree = 1\nfs = FunctionSpace(:Lagrange, degree)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"The trial space is created from the function space and the mesh. By default, a scalar \"continuous\" FESpace is created. For \"discontinuous\" (\"DG\") example, check out the linear transport tutorial.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"U = TrialFESpace(fs, mesh)\nV = TestFESpace(U)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Then, define the geometrical dommain on which the operators will apply. For this finite-element example, we only need a CellDomain (no FaceDomain).","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"domain = CellDomain(mesh)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Now, integrating on a domain necessitates a \"measure\", basically a quadrature of given degree.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"dΩ = Measure(domain, Quadrature(2 * degree + 1))","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"The definition of the two bilinear forms is quite natural. Note that these definitions are lazy, no computation is done at this step : the computations will be triggered by the assembly.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"a(u, v) = ∫(∇(u) ⋅ ∇(v))dΩ\nb(u, v) = ∫(u ⋅ v)dΩ","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"To obtain the two sparse matrices corresponding to the discretisation of these bilinear forms, simply call the assemble_bilinear function, providing the trial and test spaces.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"A = assemble_bilinear(a, U, V)\nB = assemble_bilinear(b, U, V)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Compute eigen-values and vectors : we convert to dense matrix to avoid importing additionnal packages, but it is quite easy to solve it in a \"sparse way\".","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"vp, vecp = eigen(Matrix(A), Matrix(B))","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Display the \"first\" five eigenvalues:","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"@show sqrt.(abs.(vp[3:8]))\nref_results = [\n 3.144823462554393,\n 4.447451992013584,\n 6.309054755690625,\n 6.309054755690786,\n 7.049403274103087,\n 7.049403274103147,","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"Now we can export the solution (the eigenvectors) at nodes of the mesh for several eigenvalues. We will restrict to the first 20 eigenvectors. To do so, we will create a FEFunction for each eigenvector. This FEFunction can then be evaluated on the mesh centers, nodes, etc.","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"ϕ = FEFunction(U)\nnvecs = min(20, get_ndofs(U))\nvalues = zeros(nnodes(mesh), nvecs)\nfor ivec in 1:nvecs\n set_dof_values!(ϕ, vecp[:, ivec])\n values[:, ivec] = var_on_vertices(ϕ, mesh)\nend","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"To write a VTK file, we need to build a dictionnary linking the variable name with its values and type","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"using WriteVTK\nout_dir = joinpath(@__DIR__, \"../myout\") # output directory\ndict_vars = Dict(\"u_$i\" => (values[:, i], VTKPointData()) for i in 1:nvecs)\nwrite_vtk(joinpath(out_dir, \"helmholtz_rectangle_mesh\"), 0, 0.0, mesh, dict_vars)","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"And here is the eigenvector corresponding to the 4th eigenvalue:","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"","category":"page"},{"location":"tutorial/helmholtz/","page":"Helmholtz equation (FE)","title":"Helmholtz equation (FE)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"tutorial/phase_field_supercooled/#Phase-field-model-solidification-of-a-liquid-in-supercooled-state","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"","category":"section"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"In this tutorial, a coupled system of two unsteady equations is solved using finite elements and an imex time scheme. This tutorial doesn't introduce MultiFESpace, check the \"euler\" example for this. Warning : this file is currently quite long to run (a few minutes).","category":"page"},{"location":"tutorial/phase_field_supercooled/#Theory","page":"Phase field model - solidification of a liquid in supercooled state","title":"Theory","text":"","category":"section"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"This case is taken from: Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3-4), 410-423. In particular, the variables of the problem are denoted in the same way (p for the phase indicator and T for temperature). Consider a rectangular domain Omega = 0 L_x times 0 L_y on which we wish to solve the following equations:","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":" tau partial_t p = epsilon^2 Delta p + p (1-p)(p - frac12 + m(T))","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":" partial_t T = Delta T + K partial_t p","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"where m(T) = fracalphapi atan left gamma (T_e - T) right. This set of equations represents the solidification of a liquid in a supercooled state. Here T is a dimensionless temperature and p is the solid volume fraction. Lagrange finite elements are used to discretize both equations. Time marching is performed with a forward Euler scheme for the first equation and a backward Euler scheme for the second one.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"To initiate the solidification process, a Dirichlet boundary condition (p=1,T=1) is applied at x=0 (\"West\" boundary).","category":"page"},{"location":"tutorial/phase_field_supercooled/#Code","page":"Phase field model - solidification of a liquid in supercooled state","title":"Code","text":"","category":"section"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Load the necessary packages","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"using Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing Random\n\nRandom.seed!(1234) # to obtain reproductible results","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Define some physical and numerical constants, as well as the g function appearing in the problem definition.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"const dir = string(@__DIR__, \"/../\") # Bcube dir\nconst ε = 0.01\nconst τ = 0.0003\nconst α = 0.9\nconst γ = 10.0\nconst K = 1.6\nconst Te = 1.0\nconst β = 0.0 # noise amplitude, original value : 0.01\nconst Δt = 0.0001 # time step\nconst totalTime = 1.0 # original value : 1\nconst nout = 50 # Number of iterations to skip before writing file\nconst degree = 1 # function space degree\nconst lx = 3.0\nconst ly = 1.0\nconst nx = 100\nconst ny = 20\n\ng(T) = (α / π) * atan(γ * (Te - T))","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Read the mesh using gmsh","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"const mesh_path = dir * \"input/mesh/domainPhaseField_tri.msh\"\nconst mesh = read_msh(mesh_path)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Noise function : random between [-1/2,1/2]","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"const χ = MeshCellData(rand(ncells(mesh)) .- 0.5)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Build the function space and the FE Spaces. The two unknowns will share the same FE spaces for this tutorial. Note the way we specify the Dirichlet condition in the definition of U.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"fs = FunctionSpace(:Lagrange, degree)\nU = TrialFESpace(fs, mesh, Dict(\"West\" => (x, t) -> 1.0))\nV = TestFESpace(U)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Build FE functions","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"ϕ = FEFunction(U)\nT = FEFunction(U)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Define measures for cell integration","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"dΩ = Measure(CellDomain(mesh), 2 * degree + 1)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Define bilinear and linear forms","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"a(u, v) = ∫(∇(u) ⋅ ∇(v))dΩ\nm(u, v) = ∫(u ⋅ v)dΩ\nl(v) = ∫(v * ϕ * (1.0 - ϕ) * (ϕ - 0.5 + g(T) + β * χ))dΩ","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Assemble the two constant matrices","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"A = assemble_bilinear(a, U, V)\nM = assemble_bilinear(m, U, V)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Create iterative matrices","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"C_ϕ = M + Δt / τ * ε^2 * A\nC_T = M + Δt * A","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Apply Dirichlet conditions. For this example, we don't use a lifting method to impose the Dirichlet, but d is used to initialize the solution.","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"d = assemble_dirichlet_vector(U, V, mesh)\napply_dirichlet_to_matrix!((C_ϕ, C_T), U, V, mesh)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Init solution and write it to a VTK file","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"set_dof_values!(ϕ, d)\nset_dof_values!(T, d)\n\ndict_vars = Dict(\n \"Temperature\" => (var_on_vertices(T, mesh), VTKPointData()),\n \"Phi\" => (var_on_vertices(ϕ, mesh), VTKPointData()),\n)\nwrite_vtk(dir * \"myout/result_phaseField_imex_1space\", 0, 0.0, mesh, dict_vars)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Factorize and allocate some vectors to increase performance","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"C_ϕ = factorize(C_ϕ)\nC_T = factorize(C_T)\nL = zero(d)\nrhs = zero(d)\nϕ_new = zero(d)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"Time loop (imex time integration)","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"t = 0.0\nitime = 0\nwhile t <= totalTime\n global t, itime\n t += Δt\n itime += 1\n @show t, totalTime\n\n # Integrate equation on ϕ\n L .= 0.0 # reset L\n assemble_linear!(L, l, V)\n rhs .= M * get_dof_values(ϕ) .+ Δt / τ .* L\n apply_dirichlet_to_vector!(rhs, U, V, mesh)\n ϕ_new .= C_ϕ \\ rhs\n\n # Integrate equation on T\n rhs .= M * (get_dof_values(T) .+ K .* (ϕ_new .- get_dof_values(ϕ)))\n apply_dirichlet_to_vector!(rhs, U, V, mesh)\n\n # Update solution\n set_dof_values!(ϕ, ϕ_new)\n set_dof_values!(T, C_T \\ rhs)\n\n # write solution in vtk format\n if itime % nout == 0\n dict_vars = Dict(\n \"Temperature\" => (var_on_vertices(T, mesh), VTKPointData()),\n \"Phi\" => (var_on_vertices(ϕ, mesh), VTKPointData()),\n )\n write_vtk(\n dir * \"myout/result_phaseField_imex_1space\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n end\nend","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"And here is an animation of the result:","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"","category":"page"},{"location":"tutorial/phase_field_supercooled/","page":"Phase field model - solidification of a liquid in supercooled state","title":"Phase field model - solidification of a liquid in supercooled state","text":"This page was generated using Literate.jl.","category":"page"},{"location":"tutorial/linear_transport/#Linear-transport-(DG)","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"","category":"section"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"In this tutorial, we show how to solve a linear transport equation using a discontinuous-Galerkin framework with Bcube.","category":"page"},{"location":"tutorial/linear_transport/#Theory","page":"Linear transport (DG)","title":"Theory","text":"","category":"section"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"In this example, we solve the following linear transport equation using discontinuous elements:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"fracpartial phipartial t + nabla cdot (c phi) = 0","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"where c is a constant velocity. Using an explicit time scheme, one obtains:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"phi^n+1 = phi^n - Delta t nabla cdot (c phi^n)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"The corresponding weak form of this equation is:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"int_Omega phi^n+1 v mathrmdOmega = int_Omega phi^n v mathrmdOmega + Delta t left\nint_Omega c phi^n cdot nabla v mathrmdOmega - oint_Gamma left( c phi cdot n right) v mathrmdGamma\nright","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"where Gamma = delta Omega. Adopting the discontinuous Galerkin framework, this equation is written in every mesh cell Omega_i. The cell boundary term involves discontinuous quantities and is replaced by a \"numerical flux\", leading to the expression:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"int_Omega_i phi^n+1 v mathrmdOmega_i = int_Omega_i phi^n v mathrmdOmega_i + Delta t left\nint_Omega_i c phi^n cdot nabla v mathrmdOmega_i - oint_Gamma_i F^*(phi) v mathrmd Gamma_i\nright","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"For this example, an upwind flux will be used for F^*. Using a matrix formulation, the above equation can be written as:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"phi^n+1 = phi^n + M^-1(f_Omega - f_Gamma)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"where M^-1 is the inverse of the mass matrix, f_Omega the volumic flux term and f_Gamma the surfacic flux term.","category":"page"},{"location":"tutorial/linear_transport/#Commented-code","page":"Linear transport (DG)","title":"Commented code","text":"","category":"section"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Start by importing the necessary packages: Load the necessary packages","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"using Bcube\nusing LinearAlgebra\nusing WriteVTK","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Before all, to ease to ease the solution VTK output we will write a structure to store the vtk filename and the number of iteration; and a function that exports the solution on demand. Note the use of var_on_nodes_discontinuous to export the solution on the mesh nodes, respecting the discontinuous feature of the solution.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"mutable struct VtkHandler\n basename::Any\n ite::Any\n mesh::Any\n VtkHandler(basename, mesh) = new(basename, 0, mesh)\nend\n\nfunction append_vtk(vtk, u::Bcube.AbstractFEFunction, t)\n # Values on center\n values = var_on_nodes_discontinuous(u, vtk.mesh)\n\n # Write\n Bcube.write_vtk_discontinuous(\n vtk.basename,\n vtk.ite,\n t,\n vtk.mesh,\n Dict(\"u\" => (values, VTKPointData())),\n 1;\n append = vtk.ite > 0,\n )\n\n # Update counter\n vtk.ite += 1\nend","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"First, we define some physical and numerical constant parameters","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"const degree = 0 # Function-space degree (Taylor(0) = first order Finite Volume)\nconst c = [1.0, 0.0] # Convection velocity (must be a vector)\nconst nite = 100 # Number of time iteration(s)\nconst CFL = 1 # 0.1 for degree 1\nconst nx = 41 # Number of nodes in the x-direction\nconst ny = 41 # Number of nodes in the y-direction\nconst lx = 2.0 # Domain width\nconst ly = 2.0 # Domain height\nconst Δt = CFL * min(lx / nx, ly / ny) / norm(c) # Time step","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Then generate the mesh of a rectangle using Gmsh and read it","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"tmp_path = \"tmp.msh\"\ngen_rectangle_mesh(tmp_path, :quad; nx = nx, ny = ny, lx = lx, ly = ly, xc = 0.0, yc = 0.0)\nmesh = read_msh(tmp_path)\nrm(tmp_path)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"We can now init our VtkHandler","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"out_dir = joinpath(@__DIR__, \"../myout\")\nvtk = VtkHandler(joinpath(out_dir, \"linear_transport\"), mesh)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"As seen in the previous tutorial, the definition of trial and test spaces needs a mesh and a function space. Here, we select Taylor space, and build discontinuous FE spaces with it. Then an FEFunction, that will represent our solution, is created.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"fs = FunctionSpace(:Taylor, degree)\nU = TrialFESpace(fs, mesh, :discontinuous)\nV = TestFESpace(U)\nu = FEFunction(U)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Define measures for cell and interior face integrations","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Γ = InteriorFaceDomain(mesh)\nΓ_in = BoundaryFaceDomain(mesh, \"West\")\nΓ_out = BoundaryFaceDomain(mesh, (\"North\", \"East\", \"South\"))\n\ndΩ = Measure(CellDomain(mesh), 2 * degree + 1)\ndΓ = Measure(Γ, 2 * degree + 1)\ndΓ_in = Measure(Γ_in, 2 * degree + 1)\ndΓ_out = Measure(Γ_out, 2 * degree + 1)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"We will also need the face normals associated to the different face domains. Note that this operation is lazy, nΓ is just an abstract representation on face normals of Γ.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"nΓ = get_face_normals(Γ)\nnΓ_in = get_face_normals(Γ_in)\nnΓ_out = get_face_normals(Γ_out)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Let's move on to the bilinear and linear forms. First, the two easiest ones:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"m(u, v) = ∫(u ⋅ v)dΩ # Mass matrix\nl_Ω(v) = ∫((c * u) ⋅ ∇(v))dΩ # Volumic convective term","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"For the flux term, we first need to define a numerical flux. It is convenient to define it separately in a dedicated function. Here is the definition of simple upwind flux.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"function upwind(ui, uj, nij)\n cij = c ⋅ nij\n if cij > zero(cij)\n flux = cij * ui\n else\n flux = cij * uj\n end\n flux\nend","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"We then define the \"flux\" as the composition of the upwind function and the needed entries: namely the solution on the negative side of the face, the solution on the positive face, and the face normal. The orientation negative/positive is arbitrary, the only convention is that the face normals are oriented from the negative side to the positive side.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"flux = upwind ∘ (side⁻(u), side⁺(u), side⁻(nΓ))\nl_Γ(v) = ∫(flux * jump(v))dΓ","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Finally, we define what to perform on the \"two\" boundaries : inlet / oulet. On the inlet, we directly impose the flux with a user defined function that depends on the time (the input is an oscillating wave). On the outlet, we keep our upwind flux but we impose the ghost cell value.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"bc_in = t -> PhysicalFunction(x -> c .* cos(3 * x[2]) * sin(4 * t)) # flux\nl_Γ_in(v, t) = ∫(side⁻(bc_in(t)) ⋅ side⁻(nΓ_in) * side⁻(v))dΓ_in\nflux_out = upwind ∘ (side⁻(u), 0.0, side⁻(nΓ_out))\nl_Γ_out(v) = ∫(flux_out * side⁻(v))dΓ_out","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Assemble the (constant) mass matrix. The returned matrix is a sparse matrix. To simplify the tutorial, we will directly compute the inverse mass matrix. But note that way more performant strategies should be employed to solve such a problem (since we don't need the inverse, only the matrix-vector product).","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"M = assemble_bilinear(m, U, V)\ninvM = inv(Matrix(M)) #WARNING : really expensive !!!","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"Let's also create three vectors to avoid allocating them at each time step","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"nd = get_ndofs(V)\nb_vol = zeros(nd)\nb_fac = similar(b_vol)\nrhs = similar(b_vol)","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"The time loop is trivial : at each time step we compute the linear forms using the assemble_ methods, we complete the rhs, perform an explicit step and write the solution.","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"t = 0.0\nfor i in 1:nite\n global t\n\n # Reset pre-allocated vectors\n b_vol .= 0.0\n b_fac .= 0.0\n\n # Compute linear forms\n assemble_linear!(b_vol, l_Ω, V)\n assemble_linear!(b_fac, l_Γ, V)\n assemble_linear!(b_fac, v -> l_Γ_in(v, t), V)\n assemble_linear!(b_fac, l_Γ_out, V)\n\n # Assemble rhs\n rhs .= Δt .* invM * (b_vol - b_fac)\n\n # Update solution\n u.dofValues .+= rhs\n\n # Update time\n t += Δt\n\n # Write to file\n append_vtk(vtk, u, t)\nend","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"And here is an animation of the result:","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"","category":"page"},{"location":"tutorial/linear_transport/","page":"Linear transport (DG)","title":"Linear transport (DG)","text":"This page was generated using Literate.jl.","category":"page"},{"location":"example/linear_elasticity/#Linear-elasticity","page":"Linear elasticity","title":"Linear elasticity","text":"","category":"section"},{"location":"example/linear_elasticity/","page":"Linear elasticity","title":"Linear elasticity","text":"module linear_elasticity #hide\nprintln(\"Running linear elasticity API example...\") #hide\n\n# # Linear elasticity\n\nconst dir = string(@__DIR__, \"/\") # bcube/example dir\nusing Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing StaticArrays\n\n# function space (here we shall use Lagrange P1 elements) and quadrature degree.\nconst fspace = :Lagrange\nconst degree = 1 # FunctionSpace degree\nconst degquad = 2 * degree + 1\n\n# Input and output paths\nconst outputpath = dir * \"../myout/elasticity/\"\nconst meshpath = dir * \"../input/mesh/domainElast_tri.msh\"\n\n# Time stepping scheme params\nconst α = 0.05\nconst γ = 0.5 + α\nconst β = 0.25 * (1.0 + α)^2\n\n# Material parameters (Young's modulus, Poisson coefficient and deduced Lamé coefficients)\nconst ρ = 2500.0\nconst E = 200.0e9\nconst ν = 0.3\nconst λ = E * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))\nconst μ = E / (2.0 * (1.0 + ν))\n\n# Strain tensor and stress tensor (Hooke's law)\nϵ(u) = 0.5 * (∇(u) + transpose(∇(u)))\nσ(u) = λ * tr(ϵ(u)) * I + 2 * μ * ϵ(u)\n\nπ(u, v) = σ(u) ⊡ ϵ(v) # with the chosen contraction convention ϵ should be transposed, but as it is symmetric the expression remains correct\n\n# materialize for identity operator\nBcube.materialize(A::LinearAlgebra.UniformScaling, B) = A\n\n# Function that runs the steady case:\nfunction run_steady()\n # read mesh, the second argument specifies the spatial dimension\n mesh = read_msh(meshpath, 2)\n\n fs = FunctionSpace(fspace, degree)\n U_vec = TrialFESpace(\n fs,\n mesh,\n Dict(\"West\" => SA[0.0, 0.0], \"East\" => SA[1.0, 0.0]);\n size = 2,\n )\n V_vec = TestFESpace(U_vec)\n\n # Define measures for cell\n dΩ = Measure(CellDomain(mesh), degquad)\n\n # no volume force term\n f = PhysicalFunction(x -> SA[0.0, 0.0])\n\n # definition of bilinear and linear forms\n a(u, v) = ∫(π(u, v))dΩ\n l(v) = ∫(f ⋅ v)dΩ\n\n # solve using AffineFESystem\n sys = Bcube.AffineFESystem(a, l, U_vec, V_vec)\n ϕ = Bcube.solve(sys)\n\n Un = var_on_vertices(ϕ, mesh)\n # Write the obtained FE solution\n dict_vars = Dict(\"Displacement\" => (transpose(Un), VTKPointData()))\n mkpath(outputpath)\n write_vtk(outputpath * \"result_elasticity\", itime, t, mesh, dict_vars; append = false)\nend\n\n# Function that performs a time step using a Newmark α-HHT scheme\n# The scheme updates the acceleration G, the velocity V and the displacement U using the following formulas:\n#\n# M G +(1-α)A U + αA U0 = (1-α) L + α L0 = L (because here L is time independent)\n# V = V0 + (1-γ) Δt G0 + γ Δt G\n# U = U0 + Δt V0 + (0.5-β)*Δt^2 G0 + β Δt^2 G\n#\n# G is then computed by solving the linear system obtained by inserting the expressions for U and V in the equation for G.\nfunction Newmark_α_HHT(dt, L, A, Mat, U0, V0, G0)\n L1 = L - α * A * U0\n L2 = -(1.0 - α) * (A * U0 + dt * A * V0 + (0.5 - β) * dt * dt * A * G0)\n RHS = L1 .+ L2\n\n G = Mat \\ RHS\n V = V0 + (1.0 - γ) * dt * G0 + γ * dt * G\n U = U0 + dt * V0 + (0.5 - β) * dt * dt * G0 + β * dt * dt * G\n\n return U, V, G\nend\n\n# Function that runs the unsteady case:\nfunction run_unsteady()\n # read mesh, the second argument specifies the spatial dimension\n mesh = read_msh(meshpath, 2)\n\n fs = FunctionSpace(fspace, degree)\n U_vec = TrialFESpace(fs, mesh, Dict(\"West\" => SA[0.0, 0.0]); size = 2)\n V_vec = TestFESpace(U_vec)\n\n # Define measures for cell\n dΩ = Measure(CellDomain(mesh), degquad)\n Γ = BoundaryFaceDomain(mesh, (\"East\",))\n dΓ = Measure(Γ, degquad)\n\n # surface force to be applied on East boundary\n f = PhysicalFunction(x -> SA[100000.0, 1000.0])\n\n # Definition of bilinear and linear forms\n a(u, v) = ∫(π(u, v))dΩ\n m(u, v) = ∫(ρ * u ⋅ v)dΩ\n l(v) = ∫(side⁻(f) ⋅ side⁻(v))dΓ\n\n # Assemble matrices and vector\n M = assemble_bilinear(m, U_vec, V_vec)\n A = assemble_bilinear(a, U_vec, V_vec)\n L = assemble_linear(l, V_vec)\n\n # Apply homogeneous dirichlet on A and b\n Bcube.apply_homogeneous_dirichlet_to_vector!(L, U_vec, V_vec, mesh)\n Bcube.apply_dirichlet_to_matrix!((A, M), U_vec, V_vec, mesh)\n\n # Initialize solution\n ϕ = FEFunction(U_vec, 0.0)\n U0 = zeros(Bcube.get_ndofs(U_vec))\n V0 = zeros(Bcube.get_ndofs(U_vec))\n G0 = zeros(Bcube.get_ndofs(U_vec))\n\n # Write initial solution\n Un = var_on_vertices(ϕ, mesh)\n # Write the obtained FE solution\n dict_vars = Dict(\"Displacement\" => (transpose(Un), VTKPointData()))\n mkpath(outputpath)\n write_vtk(outputpath * \"result_elasticity\", 0, 0.0, mesh, dict_vars; append = false)\n\n # Time loop\n totalTime = 1.0e-3\n Δt = 1.0e-6\n itime = 0\n t = 0.0\n\n # Matrix for time stepping\n Mat = factorize(M + (1.0 - α) * (β * Δt * Δt * A))\n\n while t <= totalTime\n t += Δt\n itime = itime + 1\n @show t, itime\n\n # solve time step\n U, V, G = Newmark_α_HHT(Δt, L, A, Mat, U0, V0, G0)\n\n # Update solution\n U0 .= U\n V0 .= V\n G0 .= G\n\n set_dof_values!(ϕ, U)\n\n # Write solution\n if itime % 10 == 0\n Un = var_on_vertices(ϕ, mesh)\n # Write the obtained FE solution\n dict_vars = Dict(\"Displacement\" => (transpose(Un), VTKPointData()))\n write_vtk(\n outputpath * \"result_elasticity\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n # In order to use the warp function in paraview (solid is deformed using the displacement field)\n # the calculator filter has to be used with the following formula to reconstruct a 3D displacement field\n # with 0 z-component: Displacement_X*iHat+Displacement_Y*jHat+0.0*kHat\n end\n end\nend\n\n#run_steady()\nrun_unsteady()\n\nend #hide","category":"page"},{"location":"howto/howto/#How-to","page":"How to...","title":"How to","text":"","category":"section"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"To be completed to answer common user questions.","category":"page"},{"location":"howto/howto/#Comparing-manually-the-benchmarks-with-main","page":"How to...","title":"Comparing manually the benchmarks with main","text":"","category":"section"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Let's say you want to compare the performance of your current branch (named \"target\" hereafter) with the main branch (named \"baseline\" hereafter).","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Open from Bcube.jl/ a REPL and type:","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"pkg> activate --temp\npkg> add PkgBenchmark BenchmarkTools\npkg> dev .\nusing PkgBenchmark\nimport Bcube\nbenchmarkpkg(Bcube, BenchmarkConfig(; env = Dict(\"JULIA_NUM_THREADS\" => \"1\")); resultfile = joinpath(@__DIR__, \"result-target.json\"))","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create a result-target.json in the current directory.","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Then checkout the main branch. Start a fresh REPL and type (almost the same):","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"pkg> activate --temp\npkg> add PkgBenchmark BenchmarkTools\npkg> dev .\nusing PkgBenchmark\nimport Bcube\nbenchmarkpkg(Bcube, BenchmarkConfig(; env = Dict(\"JULIA_NUM_THREADS\" => \"1\")); resultfile = joinpath(@__DIR__, \"result-baseline.json\"))","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create a result-baseline.json in the current directory.","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"You can now \"compare\" the two files by running (watch-out for the order):","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"target = PkgBenchmark.readresults(\"result-target.json\")\nbaseline = PkgBenchmark.readresults(\"result-baseline.json\")\njudgement = judge(target, baseline)\nexport_markdown(\"judgement.md\", judgement)","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create the markdown file judgement.md with the results.","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"For more details, once you've built the judgement object, you can also type the following code from https://github.com/tkf/BenchmarkCI.jl:","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"open(\"detailed-judgement.md\", \"w\") do io\n println(io, \"# Judge result\")\n export_markdown(io, judgement)\n println(io)\n println(io)\n println(io, \"---\")\n println(io, \"# Target result\")\n export_markdown(io, PkgBenchmark.target_result(judgement))\n println(io)\n println(io)\n println(io, \"---\")\n println(io, \"# Baseline result\")\n export_markdown(io, PkgBenchmark.baseline_result(judgement))\n println(io)\n println(io)\n println(io, \"---\")\nend","category":"page"},{"location":"howto/howto/#Run-the-benchmark-manually","page":"How to...","title":"Run the benchmark manually","text":"","category":"section"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Let's say you want to run the benchmarks locally (without comparing with main)","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"Open from Bcube.jl/ a REPL and type:","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"pkg> activate --temp\npkg> add PkgBenchmark\npkg> dev .\nusing PkgBenchmark\nimport Bcube\nresults = benchmarkpkg(Bcube, BenchmarkConfig(; env = Dict(\"JULIA_NUM_THREADS\" => \"1\")); resultfile = joinpath(@__DIR__, \"result.json\"))\nexport_markdown(\"results.md\", results)","category":"page"},{"location":"howto/howto/","page":"How to...","title":"How to...","text":"This will create the markdown file results.md with the results.","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = BcubeTutorials","category":"page"},{"location":"#Bcube","page":"Home","title":"Bcube","text":"","category":"section"},{"location":"#Purpose-of-Bcube","page":"Home","title":"Purpose of Bcube","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Bcube is a Julia library providing tools for the spatial discretization of partial differential equation(s) (PDE). The main objectives are:","category":"page"},{"location":"","page":"Home","title":"Home","text":"to provide a set of tools to quickly assemble an algorithm solving partial differential equation(s) (so the main objective is to help building prototypes without thinking about the numerical core)\nto be completed : efficient/performant PDE resolution?","category":"page"},{"location":"","page":"Home","title":"Home","text":"This documentation is organised as follow. Checkout the tutorials to see what Bcube is capable of and/or quickly learn how to use it. Then, some more elaborated examples are provided to demonstrate the library capabilities. The \"Manual\" part explains how the core is organized. Finally, the \"API\" section is the low level code documentation.","category":"page"},{"location":"#Writing-documentation","page":"Home","title":"Writing documentation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"To write documentation for Bcube, Julia's guidelines should be followed : https://docs.julialang.org/en/v1/manual/documentation/. Moreover, this project tries to apply the SciML Style Guide.","category":"page"},{"location":"#Conventions","page":"Home","title":"Conventions","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"This documentation follows the following notation or naming conventions:","category":"page"},{"location":"","page":"Home","title":"Home","text":"coordinates inside a reference frame are noted hatx haty or xi eta while coordinates in the physical frame are noted xy\nwhen talking about a mapping, F or sometimes F_rp designates the mapping from the reference element to the physical element. On the other side, F^-1 or sometimes F_pr designates the physical element to the reference element mapping.\n\"dof\" means \"degree of freedom\"","category":"page"},{"location":"example/linear_thermoelasticity/#Linear-thermo-elasticity","page":"Linear thermo-elasticity","title":"Linear thermo-elasticity","text":"","category":"section"},{"location":"example/linear_thermoelasticity/","page":"Linear thermo-elasticity","title":"Linear thermo-elasticity","text":"module linear_thermo_elasticity #hide\nprintln(\"Running linear thermo-elasticity API example...\") #hide\n\n# # Thermo-elasticity\n\nconst dir = string(@__DIR__, \"/\") # Bcube dir\nusing Bcube\nusing LinearAlgebra\nusing WriteVTK\nusing StaticArrays\n\n# function space (here we shall use Lagrange P1 elements) and quadrature degree.\nconst fspace = :Lagrange\nconst degree = 1 # FunctionSpace degree\nconst degquad = 2 * degree + 1\n\n# Input and output paths\nconst outputpath = joinpath(dir, \"../myout/elasticity/\")\nconst meshpath = joinpath(dir, \"../input/mesh/domainThermoElast_tri.msh\")\n\n# Time stepping scheme params\nconst α = 0.05\nconst γ = 0.5 + α\nconst β = 0.25 * (1.0 + α)^2\n\nconst totalTime = 10.0\nconst Δt = 1.0e-2\n\n# Material parameters (Young's modulus, Poisson coefficient and deduced Lamé coefficients)\nconst E = 200.0e9\nconst ν = 0.3\nconst λ = E * ν / ((1.0 + ν) * (1.0 - 2.0 * ν))\nconst μ = E / (2.0 * (1.0 + ν))\nconst Kₜ = 1.0e-6\nconst ρ = 2500.0\nconst cₚ = 1000.0\nconst k = 250.0\nconst T₀ = 280.0\n\n# Strain tensor and stress tensor (Hooke's law)\nϵ(u) = 0.5 * (∇(u) + transpose(∇(u)))\nσ(u) = λ * tr(ϵ(u)) * I + 2 * μ * (ϵ(u)) # Elastic stress\nσₜ(T) = (3 * λ + 2 * μ) * Kₜ * (T - T₀) * I # Thermal stress\n\nπ(u, v) = σ(u) ⊡ ϵ(v) # with the chosen contraction convention ϵ should be transposed, but as it is symmetric the expression remains correct\nπₜ(T, v) = σₜ(T) ⊡ ϵ(v)\n\n# materialize for identity operator\nBcube.materialize(A::LinearAlgebra.UniformScaling, B) = A\n\n# Function that performs a time step using a Newmark α-HHT scheme\n# The scheme updates the acceleration G, the velocity V and the displacement U using the following formulas:\n# ```math\n# \\begin{cases}\n# M G^{n+1} +(1-\\alpha)A U^{n+1} + \\alpha A U^{n} = (1-\\alpha) L^{n+1} + \\alpha L^n = L \\textrm{(because here $L$ is time independent)} \\\\\n# V^{n+1} = V^{n} + (1-\\gamma) \\Delta t G^n + \\gamma \\Delta t G^{n+1} \\\\\n# U^{n+1} = U^{n} + \\Delta t V^{n} + (\\frac{1}{2} - \\beta)*\\Delta t^2 G^{n} + \\beta \\Delta t^2 G^{n+1}\n# \\end{cases}\n# ```\n# where $$M$$ is the mass matrix, $$A$$ is the stiffness matrix and $$L$$ is the RHS\n# G is then computed by solving the linear system obtained by inserting the expressions for U and V in the equation for G.\nfunction Newmark_α_HHT(dt, L, A, Mat, U0, V0, G0)\n L1 = L - α * A * U0\n L2 = -(1.0 - α) * (A * U0 + dt * A * V0 + (0.5 - β) * dt * dt * A * G0)\n RHS = L1 .+ L2\n\n G = Mat \\ RHS\n U = U0 + dt * V0 + (0.5 - β) * dt * dt * G0 + β * dt * dt * G\n V = V0 + (1.0 - γ) * dt * G0 + γ * dt * G\n\n return U, V, G\nend\n\n# Function that runs the unsteady case:\nfunction run_unsteady()\n mesh = read_msh(meshpath, 2)\n\n fs = FunctionSpace(fspace, degree)\n U_scal = TrialFESpace(fs, mesh, Dict(\"West1\" => 280.0, \"East1\" => 280.0); size = 1)\n V_scal = TestFESpace(U_scal)\n U_vec = TrialFESpace(\n fs,\n mesh,\n Dict(\"West1\" => SA[0.0, 0.0], \"East1\" => SA[0.0, 0.0]);\n size = 2,\n )\n V_vec = TestFESpace(U_vec)\n\n # Initialize solution\n U = FEFunction(U_vec, 0.0)\n U0 = zeros(Bcube.get_ndofs(U_vec))\n V0 = zeros(Bcube.get_ndofs(U_vec))\n G0 = zeros(Bcube.get_ndofs(U_vec))\n\n T = FEFunction(U_scal, T₀)\n\n # Define measures for cell\n dΩ = Measure(CellDomain(mesh), degquad)\n\n # no volume force term\n f = PhysicalFunction(x -> SA[0.0, 0.0])\n\n q = PhysicalFunction(\n x -> x[1] .* (1.0 .- x[1]) .* x[2] .* (0.2 .- x[2]) .* 1500000000.0,\n )\n\n # Definition of bilinear and linear forms for the elasticity problem\n a(u, v) = ∫(π(u, v))dΩ\n m(u, v) = ∫(ρ * u ⋅ v)dΩ\n l(v) = ∫(πₜ(T, v))dΩ\n\n # An alternative way to define this linear form is to use operator composition:\n # l(v) = ∫( πₜ ∘ (T, v, ∇(v)) )dΩ\n # where πₜ(T, v, ∇v) = σₜ(T) ⊡ ϵ(v, ∇v) and ϵ(v, ∇v) = 0.5*( ∇v + transpose(∇v) )\n\n # Definition of bilinear and linear forms for the heat conduction problem\n aₜ(u, v) = ∫(k * ∇(u) ⋅ ∇(v))dΩ\n mₜ(u, v) = ∫(ρ * cₚ * u ⋅ v)dΩ\n lₜ(v) = ∫(q * v)dΩ\n\n # Assemble matrices and vector\n M = assemble_bilinear(m, U_vec, V_vec)\n A = assemble_bilinear(a, U_vec, V_vec)\n L = assemble_linear(l, V_vec)\n AT = assemble_bilinear(aₜ, U_scal, V_scal)\n MT = assemble_bilinear(mₜ, U_scal, V_scal)\n LT = assemble_linear(lₜ, V_scal)\n\n # Apply homogeneous dirichlet on A and b\n Bcube.apply_homogeneous_dirichlet_to_vector!(L, U_vec, V_vec, mesh)\n Bcube.apply_dirichlet_to_matrix!((A, M), U_vec, V_vec, mesh)\n\n # Compute a vector of dofs whose values are zeros everywhere\n # except on dofs lying on a Dirichlet boundary, where they\n # take the Dirichlet value\n Td = Bcube.assemble_dirichlet_vector(U_scal, V_scal, mesh)\n\n # Apply lift\n LT = LT - AT * Td\n\n # Apply homogeneous dirichlet condition\n Bcube.apply_homogeneous_dirichlet_to_vector!(LT, U_scal, V_scal, mesh)\n Bcube.apply_dirichlet_to_matrix!((AT, MT), U_scal, V_scal, mesh)\n\n # Write initial solution\n Un = var_on_vertices(U, mesh)\n Un = transpose(Un)\n Tn = var_on_vertices(T, mesh)\n mkpath(outputpath)\n dict_vars =\n Dict(\"Displacement\" => (Un, VTKPointData()), \"Temperature\" => (Tn, VTKPointData()))\n # Write the obtained FE solution\n write_vtk(\n outputpath * \"result_thermoelasticity\",\n 0,\n 0.0,\n mesh,\n dict_vars;\n append = false,\n )\n\n # Time loop\n itime = 0\n t = 0.0\n\n # Matrix for time stepping\n Mat = factorize(M + (1.0 - α) * (β * Δt * Δt * A))\n Miter = factorize(MT + Δt * AT)\n\n while t <= totalTime\n t += Δt\n itime = itime + 1\n @show t, itime\n\n # solve time step heat equation\n rhs = Δt * LT + MT * (get_dof_values(T) .- Td)\n set_dof_values!(T, Miter \\ rhs .+ Td)\n\n # solve time step elasticity\n U1, V1, G1 = Newmark_α_HHT(Δt, L, A, Mat, U0, V0, G0)\n\n # Update solution\n U0 .= U1\n V0 .= V1\n G0 .= G1\n\n set_dof_values!(U, U1)\n L = assemble_linear(l, V_vec)\n Bcube.apply_homogeneous_dirichlet_to_vector!(L, U_vec, V_vec, mesh)\n\n # Write solution\n if itime % 10 == 0\n Un = var_on_vertices(U, mesh)\n Un = transpose(Un)\n Tn = var_on_vertices(T, mesh)\n mkpath(outputpath)\n dict_vars = Dict(\n \"Displacement\" => (Un, VTKPointData()),\n \"Temperature\" => (Tn, VTKPointData()),\n )\n # Write the obtained FE solution\n write_vtk(\n outputpath * \"result_thermoelasticity\",\n itime,\n t,\n mesh,\n dict_vars;\n append = true,\n )\n # In order to use the warp function in paraview (solid is deformed using the displacement field)\n # the calculator filter has to be used with the following formula to reconstruct a 3D displacement field\n # with 0 z-component: Displacement_X*iHat+Displacement_Y*jHat+0.0*kHat\n end\n end\nend\n\nrun_unsteady()\n\n# Here is an animation of the obtained result:\n# ```@raw html\n# \n# ```\n\nend #hide","category":"page"}]
}
diff --git a/previews/PR1/tutorial/heat_equation/index.html b/previews/PR1/tutorial/heat_equation/index.html
index 4547dff..b6d44a8 100644
--- a/previews/PR1/tutorial/heat_equation/index.html
+++ b/previews/PR1/tutorial/heat_equation/index.html
@@ -1,5 +1,5 @@
-In this tutorial, the heat equation (first steady and then unsteady) is solved using finite-elements.
This example shows how to solve the heat equation with eventually variable physical properties in steady and unsteady formulations:
\[ \rho C_p \partial_t u - \nabla . ( \lambda 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:
\[ \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)}\]
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.
As usual, start by importing the necessary packages.
using Bcube
+Heat equation (FE) · BcubeTutorials Heat equation (FE)
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:
\[ \rho C_p \partial_t u - \nabla . ( \lambda 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:
\[ \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)}\]
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
As usual, start by importing the necessary packages.
using Bcube
using LinearAlgebra
using WriteVTK
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)
@@ -60,4 +60,4 @@
)
end
end
-
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
This tutorial shows how to solve the Helmholtz eigen problem with a finite-element approach using Bcube.
We consider the following Helmholtz equation, representing for instance the acoustic wave propagation with Neuman boundary condition(s):
\[\begin{cases} +
This tutorial shows how to solve the Helmholtz eigen problem with a finite-element approach using Bcube.
We consider the following Helmholtz equation, representing for instance the acoustic wave propagation with Neuman boundary condition(s):
\[\begin{cases} \Delta u + \omega^2 u = 0 \\ \dfrac{\partial u}{\partial n} = 0 \textrm{ on } \Gamma \end{cases}\]
An analytic solution of this equation can be obtained: for a rectangular domain $\Omega = [0,L_x] \times [0,L_y]$,
\[u(x,y) = \cos \left( \frac{k_x \pi}{L_x} x \right) \cos \left( \frac{k_y \pi}{L_y} y \right) \mathrm{~~with~~} k_x,~k_y \in \mathbb{N}\]
with $\omega^2 = \pi^2 \left( \dfrac{k_x^2}{L_x^2} + \dfrac{k_y^2}{L_y^2} \right)$
Now, both the finite-element method and the discontinuous Galerkin method requires to write the weak form of the problem:
\[\int_\Omega (\Delta u \Delta v + \omega^2u)v \mathrm{\,d}\Omega = 0\]
\[- \int_\Omega \nabla u \cdot \nabla v \mathrm{\,d}\Omega @@ -43,4 +43,4 @@ end
To write a VTK file, we need to build a dictionnary linking the variable name with its values and type
using WriteVTK
out_dir = joinpath(@__DIR__, "../myout") # output directory
dict_vars = Dict("u_$i" => (values[:, i], VTKPointData()) for i in 1:nvecs)
-write_vtk(joinpath(out_dir, "helmholtz_rectangle_mesh"), 0, 0.0, mesh, dict_vars)And here is the eigenvector corresponding to the 4th eigenvalue:
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
And here is the eigenvector corresponding to the 4th eigenvalue:
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
In this tutorial, we show how to solve a linear transport equation using a discontinuous-Galerkin framework with Bcube.
In this example, we solve the following linear transport equation using discontinuous elements:
\[\frac{\partial \phi}{\partial t} + \nabla \cdot (c \phi) = 0\]
where $c$ is a constant velocity. Using an explicit time scheme, one obtains:
\[\phi^{n+1} = \phi^n - \Delta t \nabla \cdot (c \phi^n)\]
The corresponding weak form of this equation is:
\[\int_\Omega \phi^{n+1} v \mathrm{\,d}\Omega = \int_\Omega \phi^n v \mathrm{\,d}\Omega + \Delta t \left[ +
In this tutorial, we show how to solve a linear transport equation using a discontinuous-Galerkin framework with Bcube.
In this example, we solve the following linear transport equation using discontinuous elements:
\[\frac{\partial \phi}{\partial t} + \nabla \cdot (c \phi) = 0\]
where $c$ is a constant velocity. Using an explicit time scheme, one obtains:
\[\phi^{n+1} = \phi^n - \Delta t \nabla \cdot (c \phi^n)\]
The corresponding weak form of this equation is:
\[\int_\Omega \phi^{n+1} v \mathrm{\,d}\Omega = \int_\Omega \phi^n v \mathrm{\,d}\Omega + \Delta t \left[ \int_\Omega c \phi^n \cdot \nabla v \mathrm{\,d}\Omega - \oint_\Gamma \left( c \phi \cdot n \right) v \mathrm{\,d}\Gamma \right]\]
where $\Gamma = \delta \Omega$. Adopting the discontinuous Galerkin framework, this equation is written in every mesh cell $\Omega_i$. The cell boundary term involves discontinuous quantities and is replaced by a "numerical flux", leading to the expression:
\[\int_{\Omega_i} \phi^{n+1} v \mathrm{\,d}\Omega_i = \int_{\Omega_i} \phi^n v \mathrm{\,d}\Omega_i + \Delta t \left[ \int_{\Omega_i} c \phi^n \cdot \nabla v \mathrm{\,d}\Omega_i - \oint_{\Gamma_i} F^*(\phi) v \mathrm{\,d} \Gamma_i @@ -95,4 +95,4 @@ # Write to file append_vtk(vtk, u, t) -end
And here is an animation of the result:
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
And here is an animation of the result:
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
In this tutorial, a coupled system of two unsteady equations is solved using finite elements and an imex time scheme. This tutorial doesn't introduce MultiFESpace, check the "euler" example for this. Warning : this file is currently quite long to run (a few minutes).
This case is taken from: Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3-4), 410-423. In particular, the variables of the problem are denoted in the same way ($p$ for the phase indicator and $T$ for temperature). Consider a rectangular domain $\Omega = [0, L_x] \times [0, L_y]$ on which we wish to solve the following equations:
\[ \tau \partial_t p = \epsilon^2 \Delta p + p (1-p)(p - \frac{1}{2} + m(T))\]
\[ \partial_t T = \Delta T + K \partial_t p\]
where $m(T) = \frac{\alpha}{\pi} atan \left[ \gamma (T_e - T) \right]$. This set of equations represents the solidification of a liquid in a supercooled state. Here $T$ is a dimensionless temperature and $p$ is the solid volume fraction. Lagrange finite elements are used to discretize both equations. Time marching is performed with a forward Euler scheme for the first equation and a backward Euler scheme for the second one.
To initiate the solidification process, a Dirichlet boundary condition ($p=1$,$T=1$) is applied at $x=0$ ("West" boundary).
Load the necessary packages
using Bcube
+Phase field model - solidification of a liquid in supercooled state · BcubeTutorials Phase field model - solidification of a liquid in supercooled state
In this tutorial, a coupled system of two unsteady equations is solved using finite elements and an imex time scheme. This tutorial doesn't introduce MultiFESpace, check the "euler" example for this. Warning : this file is currently quite long to run (a few minutes).
Theory
This case is taken from: Kobayashi, R. (1993). Modeling and numerical simulations of dendritic crystal growth. Physica D: Nonlinear Phenomena, 63(3-4), 410-423. In particular, the variables of the problem are denoted in the same way ($p$ for the phase indicator and $T$ for temperature). Consider a rectangular domain $\Omega = [0, L_x] \times [0, L_y]$ on which we wish to solve the following equations:
\[ \tau \partial_t p = \epsilon^2 \Delta p + p (1-p)(p - \frac{1}{2} + m(T))\]
\[ \partial_t T = \Delta T + K \partial_t p\]
where $m(T) = \frac{\alpha}{\pi} atan \left[ \gamma (T_e - T) \right]$. This set of equations represents the solidification of a liquid in a supercooled state. Here $T$ is a dimensionless temperature and $p$ is the solid volume fraction. Lagrange finite elements are used to discretize both equations. Time marching is performed with a forward Euler scheme for the first equation and a backward Euler scheme for the second one.
To initiate the solidification process, a Dirichlet boundary condition ($p=1$,$T=1$) is applied at $x=0$ ("West" boundary).
Code
Load the necessary packages
using Bcube
using LinearAlgebra
using WriteVTK
using Random
@@ -79,4 +79,4 @@
append = true,
)
end
-end
And here is an animation of the result:
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.
And here is an animation of the result:
This page was generated using Literate.jl.
Settings
This document was generated with Documenter.jl version 1.1.2 on Sunday 19 November 2023. Using Julia version 1.9.4.