diff --git a/.github/workflows/website.yml b/.github/workflows/website.yml index 1765a5c..5d6ae31 100644 --- a/.github/workflows/website.yml +++ b/.github/workflows/website.yml @@ -30,6 +30,7 @@ jobs: run: | mkdir -p public cp introduction_to_julia/Introduction_to_Julia.html public/Introduction_to_Julia.html + cp introduction_to_trixi/Introduction_to_Trixi.html public/Introduction_to_Trixi.html - name: Upload static files as artifact id: deployment uses: actions/upload-pages-artifact@v3 diff --git a/README.md b/README.md index 7a75762..fa2da50 100644 --- a/README.md +++ b/README.md @@ -27,3 +27,11 @@ The material for the brief introduction to Julia is contained in the directory [Pluto.jl](https://github.com/fonsp/Pluto.jl) notebook that can be updated interactively. You can also see a [static preview online](https://trixi-framework.github.io/talk-2025-Julia_and_Trixi_in_Frankfurt/Introduction_to_Julia.html). + +### Introduction to Trixi.jl + +The material for the brief introduction to Trixi.jl is contained in the directory +`introduction_to_trixi` in this repository. It is written as a +[Pluto.jl](https://github.com/fonsp/Pluto.jl) notebook that can be updated +interactively. You can also see a +[static preview online](https://trixi-framework.github.io/talk-2025-Julia_and_Trixi_in_Frankfurt/Introduction_to_Trixi.html). diff --git a/introduction_to_trixi/Introduction_to_Trixi.html b/introduction_to_trixi/Introduction_to_Trixi.html new file mode 100644 index 0000000..06e7672 --- /dev/null +++ b/introduction_to_trixi/Introduction_to_Trixi.html @@ -0,0 +1,17 @@ + + + + + + + +
\ No newline at end of file diff --git a/introduction_to_trixi/Introduction_to_Trixi.jl b/introduction_to_trixi/Introduction_to_Trixi.jl new file mode 100644 index 0000000..4b8e89e --- /dev/null +++ b/introduction_to_trixi/Introduction_to_Trixi.jl @@ -0,0 +1,3371 @@ +### A Pluto.jl notebook ### +# v0.20.3 + +using Markdown +using InteractiveUtils + +# ╔═╡ 1a49968b-3b5d-4eb0-b568-122c758a49cd +begin + using PlutoUI + using PlutoUI: Slider + TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true) +end + +# ╔═╡ bd3eea53-7477-448c-b4c8-e220b49e493c +using LinearAlgebra + +# ╔═╡ 208484e0-07fd-4e13-b046-39bea16d9b9b +using Trixi + +# ╔═╡ b2d22b03-00d0-4c6c-bec8-c747835187eb +begin + using Plots + using Plots: Plots, plot, plot!, savefig +end + +# ╔═╡ 771bb8bd-73b1-423b-aa24-2f335e5f8090 +using LaTeXStrings # for L"..." + +# ╔═╡ 275a3f6e-ce7c-4f33-9cd7-7063baff343e +using OrdinaryDiffEq + +# ╔═╡ e97e43e1-1035-4501-a23b-ade6ef3b0931 +begin + # PassiveTracerEquations setup for Kelvin-Helmholtz instability + using Trixi: AbstractEquations, rotate_to_x + import Trixi: nvariables, varnames, flux, + cons2prim, prim2cons, density, pressure, density_pressure, + max_abs_speeds, max_abs_speed_naive, entropy, cons2entropy + + struct PassiveTracerEquations{NDIMS, NVARS, NTRACERS, + FlowEquations <: AbstractEquations{NDIMS}} <: AbstractEquations{NDIMS, NVARS} + flow_equations::FlowEquations + # TODO: Ensure that the number of variables fits + end + + function PassiveTracerEquations(flow_equations::AbstractEquations, ::Val{NTRACERS}) where {NTRACERS} + NVARS = nvariables(flow_equations) + NTRACERS + return PassiveTracerEquations{ndims(flow_equations), NVARS, NTRACERS, typeof(flow_equations)}(flow_equations) + end + + @inline ntracers(::PassiveTracerEquations{NDIMS, NVARS, NTRACERS}) where {NDIMS, NVARS, NTRACERS} = NTRACERS + + function varnames(variables::typeof(cons2cons), + tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + flow_varnames = varnames(variables, flow_equations) + n_tracers = ntracers(tracer_equations) + return (flow_varnames..., ntuple(i -> "rho_xi_$i", Val(n_tracers))...) + end + + function varnames(variables::typeof(cons2prim), + tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + flow_varnames = varnames(variables, flow_equations) + n_tracers = ntracers(tracer_equations) + return (flow_varnames..., ntuple(i -> "rho_xi_$i", Val(n_tracers))...) + end + + @inline function nvariables_flow(tracer_equations::PassiveTracerEquations) + return nvariables(tracer_equations.flow_equations) + end + + @inline function flow_variables(u, tracer_equations::PassiveTracerEquations) + n_flow_variables = nvariables_flow(tracer_equations) + return SVector(ntuple(i -> u[i], Val(n_flow_variables))...) + end + + @inline function tracers(u, tracer_equations::PassiveTracerEquations) + n_flow_variables = nvariables_flow(tracer_equations) + n_tracers = ntracers(tracer_equations) + + rho = density(u, tracer_equations) + return SVector(ntuple(i -> u[n_flow_variables + i] / rho, Val(n_tracers))...) + end + + @inline function cons2prim(u, tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + return SVector(cons2prim(flow_variables(u, tracer_equations), flow_equations)..., + tracers(u, tracer_equations)...) + end + + @inline function prim2cons(u, tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + u_flow = flow_variables(u, tracer_equations) + + n_flow_variables = nvariables_flow(tracer_equations) + n_tracers = ntracers(tracer_equations) + + rho = density(u, tracer_equations) + return SVector(prim2cons(u_flow, flow_equations)..., + (rho * u[n_flow_variables + i] for i=1:n_tracers)...) + end + + for f in [:density, :pressure, :density_pressure, :max_abs_speeds] + @eval @inline function $f(u, tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + u_flow = flow_variables(u, tracer_equations) + return $f(u_flow, flow_equations) + end + end + + @inline function max_abs_speed_naive(u_ll, u_rr, orientation_or_normal, + tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + u_flow_ll = flow_variables(u_ll, tracer_equations) + u_flow_rr = flow_variables(u_rr, tracer_equations) + return max_abs_speed_naive(u_flow_ll, u_flow_rr, orientation_or_normal, flow_equations) + end + + @inline function entropy(cons, tracer_equations::PassiveTracerEquations) + flow_entropy = entropy(flow_variables(cons, tracer_equations), + tracer_equations.flow_equations) + tracer_entropy = 0.5f0 * density(cons, tracer_equations) * sum(abs2, tracers(cons, tracer_equations)) + return flow_entropy + tracer_entropy + end + + @inline function cons2entropy(u, tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + flow_entropy = cons2entropy(flow_variables(u, tracer_equations), flow_equations) + tracers_ = tracers(u, tracer_equations) + + variables = SVector(flow_entropy[1] - 0.5f0 * sum(abs2, tracers_), + (flow_entropy[i] for i=2:nvariables(flow_equations))..., + tracers_...) + return variables + end + + # Calculate flux for a single point + @inline function flux(u, orientation::Integer, + tracer_equations::PassiveTracerEquations) + n_flow_variables = nvariables_flow(tracer_equations) + + u_flow = flow_variables(u, tracer_equations) + + # Awkward way to fix box type instability when using `if` + # TODO: We should fix this when moving it to Trixi.jl. + # Ideally, we could have an interface `velocity(u, equations)` + # that we just need to use here. + v_normal = u[1 + orientation] / density(u, tracer_equations) + + flux_tracer = SVector(ntuple(@inline(v->u[v+n_flow_variables]*v_normal), + Val(ntracers(tracer_equations)))) + + flux_flow = flux(u_flow, orientation, tracer_equations.flow_equations) + + return SVector(flux_flow..., flux_tracer...) + end + + @inline function flux(u, normal_direction::AbstractVector, + tracer_equations::PassiveTracerEquations) + n_flow_variables = nvariables_flow(tracer_equations) + + u_flow = flow_variables(u, tracer_equations) + + # Awkward way to fix box type instability when using `if` + # TODO: We should fix this when moving it to Trixi.jl. + # Ideally, we could have an interface `velocity(u, equations)` + # that we just need to use here. + v_normal = (normal_direction[1] * u[2] + normal_direction[2] * u[3]) / density(u, tracer_equations) + + flux_tracer = SVector(ntuple(@inline(v->u[v+n_flow_variables]*v_normal), + Val(ntracers(tracer_equations)))) + + flux_flow = flux(u_flow, normal_direction, tracer_equations.flow_equations) + + return SVector(flux_flow..., flux_tracer...) + end + + struct FluxTracerEquationsCentral{FlowFlux} + flow_flux::FlowFlux + end + + @inline function (f::FluxTracerEquationsCentral)(u_ll, u_rr, + orientation_or_normal_direction, + tracer_equations::PassiveTracerEquations) + (; flow_equations) = tracer_equations + u_flow_ll = flow_variables(u_ll, tracer_equations) + u_flow_rr = flow_variables(u_rr, tracer_equations) + n_flow_variables = nvariables_flow(tracer_equations) + + flux_flow = f.flow_flux(u_flow_ll, u_flow_rr, orientation_or_normal_direction, + flow_equations) + flux_rho = density(flux_flow, flow_equations) + rho_ll = density(u_ll, flow_equations) + rho_rr = density(u_rr, flow_equations) + flux_tracer = SVector(ntuple(@inline(v -> u_ll[v + n_flow_variables] / rho_ll + + u_rr[v + n_flow_variables] / rho_rr), + Val(ntracers(tracer_equations)))) + flux_tracer = flux_rho * 0.5f0 * flux_tracer + return SVector(flux_flow..., flux_tracer...) + end + + + @inline function boundary_condition_slip_wall(u_inner, orientation::Int, + direction, x, t, + surface_flux_function, + equations::PassiveTracerEquations) + RealT = eltype(u_inner) + if orientation == 1 + normal_direction = SVector(one(RealT), zero(RealT)) + else # orientation == 2 + normal_direction = SVector(zero(RealT), one(RealT)) + end + + # compute and return the flux using `boundary_condition_slip_wall` routine above + return boundary_condition_slip_wall(u_inner, normal_direction, direction, + x, t, surface_flux_function, equations) + end + + @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + direction, x, t, + surface_flux_function, + equations::PassiveTracerEquations) + # flip sign of normal to make it outward pointing, then flip the sign of the normal flux back + # to be inward pointing on the -x and -y sides due to the orientation convention used by StructuredMesh + if isodd(direction) + boundary_flux = -boundary_condition_slip_wall(u_inner, -normal_direction, + x, t, surface_flux_function, + equations) + else + boundary_flux = boundary_condition_slip_wall(u_inner, normal_direction, + x, t, surface_flux_function, + equations) + end + + return boundary_flux + end + + @inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, + x, t, + surface_flux_function, + tracer_equations::PassiveTracerEquations) + flow_equations = tracer_equations.flow_equations + norm_ = norm(normal_direction) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + normal = normal_direction / norm_ + + # rotate the internal solution state + u_inner_flow = flow_variables(u_inner, tracer_equations) + u_local = rotate_to_x(u_inner_flow, normal, flow_equations) + + # compute the primitive variables + rho_local, v_normal, v_tangent, p_local = cons2prim(u_local, flow_equations) + + # Get the solution of the pressure Riemann problem + # See Section 6.3.3 of + # Eleuterio F. Toro (2009) + # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction + # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761) + if v_normal <= 0.0 + sound_speed = sqrt(flow_equations.gamma * p_local / rho_local) # local sound speed + p_star = p_local * + (1 + 0.5 * (flow_equations.gamma - 1) * v_normal / sound_speed)^(2 * + flow_equations.gamma * + flow_equations.inv_gamma_minus_one) + else # v_normal > 0.0 + A = 2 / ((flow_equations.gamma + 1) * rho_local) + B = p_local * (flow_equations.gamma - 1) / (flow_equations.gamma + 1) + p_star = p_local + + 0.5 * v_normal / A * + (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B))) + end + + # For the slip wall we directly set the flux as the normal velocity is zero + return SVector(zero(eltype(u_inner_flow)), + p_star * normal[1], + p_star * normal[2], + zero(eltype(u_inner_flow)), + v_normal * u_inner[5] / u_inner[1]) * norm_ + end + # End of passive tracer extension + ############################################################################### + + + const figdir = joinpath(@__DIR__, "figures") + if !isdir(figdir) + mkdir(figdir) + end + + fontsizes() = (xtickfontsize = 18, + ytickfontsize = 18, + legendfontsize = 18, + xguidefontsize = 20, + yguidefontsize = 20, + titlefontsize = 20,) + + + function cons2density_and_tracers(u, tracer_equations::PassiveTracerEquations) + return SVector(density(u, tracer_equations), tracers(u, tracer_equations)...) + end + + function varnames(variables::typeof(cons2density_and_tracers), + tracer_equations::PassiveTracerEquations) + n_tracers = ntracers(tracer_equations) + return ("rho", ntuple(i -> "xi_$i", Val(n_tracers))...) + end +end + +# ╔═╡ 2d612d09-5f5e-4947-944a-269dc5bf7116 +html""" +

Introduction to Trixi.jl

+
Arpit Babbar
+
Postdoctoral researcher
+
Numerical Mathematics, JGU Mainz
+""" + +# ╔═╡ ed2d8c29-043f-41ee-bddf-dc0792abfb15 +html""" +

General system of equations

+""" + +# ╔═╡ 7fcc373f-ba04-4df9-a4ed-49313999ae6f +begin +md"""$$\Large \boldsymbol{u}_t + \nabla \cdot \boldsymbol{f}_a(\boldsymbol{u}) + A(\boldsymbol{u}) \nabla \boldsymbol{u} = \nabla \cdot \boldsymbol{f}_v(\boldsymbol{u}, \nabla \boldsymbol{u}) + \boldsymbol{s}(\boldsymbol{u})$$ +""" +end # LaTeX \newcommand does not work when exporting as HTML + +# ╔═╡ b94d039d-1114-406a-87da-63f955e68188 +html""" +

Trixi.jl in action
+Flow-Gravity +

+""" + +# ╔═╡ 1b323e44-0362-4c5b-8a54-f85f0c3d67c9 +html"""