diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 10bcc22786a..8f6d12179eb 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.21.0 + uses: crate-ci/typos@v1.22.9 diff --git a/AUTHORS.md b/AUTHORS.md index 54d63216335..5ab164c0ed0 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -38,6 +38,7 @@ are listed in alphabetical order: * Julia Odenthal * Sigrun Ortleb * Hendrik Ranocha +* Warisa Roongaraya * Andrés M. Rueda-Ramírez * Felipe Santillan * Michael Schlottke-Lakemper diff --git a/NEWS.md b/NEWS.md index ecbd70ce472..bca7d03e954 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,20 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes when updating to v0.8 from v0.7.x + +#### Added + +#### Changed + +- The specification of boundary names on which `AnalysisSurfaceIntegral`s are computed (such as drag and lift coefficients) has changed from `Symbol` and `Vector{Symbol}` to `NTuple{Symbol}`. + Thus, for one boundary the syntax changes from `:boundary` to `(:boundary,)` and for `Vector`s `[:boundary1, :boundary2]` to `(:boundary1, :boundary2)` ([#1959]). +- The names of output files like the one created from the `SaveSolutionCallback` have changed from `%06d` to `%09d` to allow longer-running simulations ([#1996]). + +#### Deprecated + +#### Removed + ## Changes in the v0.7 lifecycle #### Added diff --git a/Project.toml b/Project.toml index a6911d4e2c3..f942de89309 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -version = "0.7.16-pre" +version = "0.8.2-DEV" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" diff --git a/docs/Project.toml b/docs/Project.toml index 9f9bb956274..f3a1e66d788 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -25,5 +25,5 @@ Measurements = "2.5" OrdinaryDiffEq = "6.49.1" Plots = "1.9" Test = "1" -Trixi2Vtk = "0.3" +Trixi2Vtk = "0.3.16" TrixiBase = "0.1.1" diff --git a/docs/literate/src/files/first_steps/create_first_setup.jl b/docs/literate/src/files/first_steps/create_first_setup.jl index ae78a6a1546..227e4cd2cce 100644 --- a/docs/literate/src/files/first_steps/create_first_setup.jl +++ b/docs/literate/src/files/first_steps/create_first_setup.jl @@ -260,28 +260,28 @@ plot!(getmesh(pd)) # import Pkg # Pkg.add(["Trixi2Vtk"]) # ``` -# Now we load the Trixi2Vtk.jl package and convert the file `out/solution_000032.h5` with +# Now we load the Trixi2Vtk.jl package and convert the file `out/solution_000000032.h5` with # the final solution using the [`trixi2vtk`](@ref) function saving the resulting file in the # `out` folder. using Trixi2Vtk -trixi2vtk(joinpath("out", "solution_000032.h5"), output_directory="out") +trixi2vtk(joinpath("out", "solution_000000032.h5"), output_directory="out") -# Now two files `solution_000032.vtu` and `solution_000032_celldata.vtu` have been generated in the +# Now two files `solution_000000032.vtu` and `solution_000000032_celldata.vtu` have been generated in the # `out` folder. The first one contains all the information for visualizing the solution, the # second one contains all the cell-based or discretization-based information. # Now let's visualize the solution from the generated files in ParaView. Follow this short # instruction to get the visualization. # - Download, install and open [ParaView](https://www.paraview.org/download/). -# - Press `Ctrl+O` and select the generated files `solution_000032.vtu` and -# `solution_000032_celldata.vtu` from the `out` folder. +# - Press `Ctrl+O` and select the generated files `solution_000000032.vtu` and +# `solution_000000032_celldata.vtu` from the `out` folder. # - In the upper-left corner in the Pipeline Browser window, left-click on the eye-icon near -# `solution_000032.vtu`. +# `solution_000000032.vtu`. # - In the lower-left corner in the Properties window, change the Coloring from Solid Color to # scalar. This already generates the visualization of the final solution. # - Now let's add the mesh to the visualization. In the upper-left corner in the -# Pipeline Browser window, left-click on the eye-icon near `solution_000032_celldata.vtu`. +# Pipeline Browser window, left-click on the eye-icon near `solution_000000032_celldata.vtu`. # - In the lower-left corner in the Properties window, change the Representation from Surface # to Wireframe. Then a white grid should appear on the visualization. # Now, if you followed the instructions exactly, you should get a similar image as shown in the diff --git a/docs/literate/src/files/hohqmesh_tutorial.jl b/docs/literate/src/files/hohqmesh_tutorial.jl index dd81f47951e..f5d93428f10 100644 --- a/docs/literate/src/files/hohqmesh_tutorial.jl +++ b/docs/literate/src/files/hohqmesh_tutorial.jl @@ -54,7 +54,7 @@ end #hide #md using Trixi2Vtk redirect_stdio(stdout=devnull, stderr=devnull) do # code that prints annoying stuff we don't want to see here #hide #md -trixi2vtk("out/solution_000180.h5", output_directory="out") +trixi2vtk("out/solution_000000180.h5", output_directory="out") end #hide #md # Note this step takes about 15-30 seconds as the package `Trixi2Vtk` must be precompiled and executed for the first time @@ -64,7 +64,7 @@ end #hide #md # visualization nodes. For instance, if we want to use 12 uniformly spaced nodes for visualization we can execute redirect_stdio(stdout=devnull, stderr=devnull) do # code that prints annoying stuff we don't want to see here #hide #md -trixi2vtk("out/solution_000180.h5", output_directory="out", nvisnodes=12) +trixi2vtk("out/solution_000000180.h5", output_directory="out", nvisnodes=12) end #hide #md # By default `trixi2vtk` sets `nvisnodes` to be the same as the number of nodes specified in diff --git a/docs/make.jl b/docs/make.jl index 73ee86abd8d..e0243a8bca0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -102,11 +102,21 @@ tutorials = create_tutorials(files) # Create changelog Changelog.generate( - Changelog.Documenter(), # output type - joinpath(@__DIR__, "..", "NEWS.md"), # input file - joinpath(@__DIR__, "src", "changelog.md"); # output file - repo = "trixi-framework/Trixi.jl", # default repository for links + Changelog.Documenter(), # output type + joinpath(@__DIR__, "..", "NEWS.md"), # input file + joinpath(@__DIR__, "src", "changelog_tmp.md"); # output file + repo = "trixi-framework/Trixi.jl", # default repository for links + branch = "main", # default branch for links ) +# Fix edit URL of changelog +open(joinpath(@__DIR__, "src", "changelog.md"), "w") do io + for line in eachline(joinpath(@__DIR__, "src", "changelog_tmp.md")) + if startswith(line, "EditURL") + line = "EditURL = \"https://github.com/trixi-framework/Trixi.jl/blob/main/NEWS.md\"" + end + println(io, line) + end +end # Make documentation makedocs( diff --git a/docs/src/restart.md b/docs/src/restart.md index c7cbcd11852..fd2044b6a11 100644 --- a/docs/src/restart.md +++ b/docs/src/restart.md @@ -30,7 +30,7 @@ conditionals like `if restart` with a boolean variable `restart` that is user de First we need to define from which file we want to restart, e.g. ```julia -restart_file = "restart_000021.h5" +restart_file = "restart_000000021.h5" restart_filename = joinpath("out", restart_file) ``` diff --git a/docs/src/visualization.md b/docs/src/visualization.md index 36a7e8f5ac8..95f7ed6b684 100644 --- a/docs/src/visualization.md +++ b/docs/src/visualization.md @@ -415,15 +415,15 @@ julia> using Trixi2Vtk ``` To process an HDF5 file generated by Trixi.jl, execute ```julia -julia> trixi2vtk(joinpath("out", "solution_000000.h5"), output_directory="out") +julia> trixi2vtk(joinpath("out", "solution_000000000.h5"), output_directory="out") ``` This will create two unstructured VTK files in the `out` subdirectory that can -be opened with ParaView or VisIt: `solution_000000.vtu` contains the -discontinuous Galerkin solution data while `solution_000000_celldata.vtu` holds +be opened with ParaView or VisIt: `solution_000000000.vtu` contains the +discontinuous Galerkin solution data while `solution_000000000_celldata.vtu` holds any cell-based values such as the current AMR indicator or the cell refinement level. -!["solution_000000_scalar_mesh"](https://github.com/trixi-framework/Trixi2Vtk.jl/raw/main/docs/src/assets/solution_000000_scalar_mesh.png) +!["solution_000000000_scalar_mesh"](https://github.com/trixi-framework/Trixi2Vtk.jl/raw/main/docs/src/assets/solution_000000_scalar_mesh.png) This allows you to generate VTK files for solution, restart and mesh files. By default, Trixi2Vtk generates `.vtu` (unstructured VTK) files for both @@ -440,7 +440,7 @@ publication-quality images, but increases the file size. If you want to convert multiple solution/restart files at once, you can just supply multiple input files as the positional arguments to `trixi2vtk`, e.g., ```julia -julia> trixi2vtk("out/solution_000000.h5", "out/solution_000040.h5") +julia> trixi2vtk("out/solution_000000000.h5", "out/solution_000000040.h5") ``` You may also use file globbing to select a range of files based on filename patterns, e.g., ```julia @@ -450,7 +450,7 @@ to convert all solution files in the `out/` directory or ```julia julia> trixi2vtk("out/restart_00[0-9]000.h5") ``` -to convert every one-thousandth restart file (`out/restart_000000.h5`, +to convert every one-thousandth restart file (`out/restart_000000000.h5`, `out/restart_001000.h5` etc.). When multiple solution/restart files are provided, Trixi2Vtk will also generate a diff --git a/examples/p4est_2d_dgsem/elixir_advection_restart.jl b/examples/p4est_2d_dgsem/elixir_advection_restart.jl index 0f573714c1f..6fb14c60038 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_restart.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_restart.jl @@ -6,7 +6,7 @@ using Trixi # create a restart file elixir_file = "elixir_advection_extended.jl" -restart_file = "restart_000021.h5" +restart_file = "restart_000000021.h5" trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file)) diff --git a/examples/p4est_2d_dgsem/elixir_advection_restart_amr.jl b/examples/p4est_2d_dgsem/elixir_advection_restart_amr.jl index fd3623dd88b..6b68347493a 100644 --- a/examples/p4est_2d_dgsem/elixir_advection_restart_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_advection_restart_amr.jl @@ -6,7 +6,7 @@ using Trixi # create a restart file elixir_file = "elixir_advection_extended.jl" -restart_file = "restart_000021.h5" +restart_file = "restart_000000021.h5" trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file)) diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index fb5f29bd038..8806f9cc2af 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -84,7 +84,7 @@ analysis_interval = 2000 l_inf = 1.0 # Length of airfoil -force_boundary_names = [:AirfoilBottom, :AirfoilTop] +force_boundary_names = (:AirfoilBottom, :AirfoilTop) drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, DragCoefficientPressure(aoa(), rho_inf(), u_inf(equations), l_inf)) diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index dc23e192de8..e522a0dcfd3 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -89,11 +89,11 @@ rho_inf = 1.4 u_inf = 0.38 l_inf = 1.0 # Diameter of circle -drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, +drag_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,), DragCoefficientPressure(aoa, rho_inf, u_inf, l_inf)) -lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg, +lift_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,), LiftCoefficientPressure(aoa, rho_inf, u_inf, l_inf)) diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 1b485913ab2..3e4679f5afe 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -119,7 +119,7 @@ summary_callback = SummaryCallback() analysis_interval = 2000 -force_boundary_names = [:AirfoilBottom, :AirfoilTop] +force_boundary_names = (:AirfoilBottom, :AirfoilTop) drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names, DragCoefficientPressure(aoa(), rho_inf(), u_inf(equations), diff --git a/examples/p4est_3d_dgsem/elixir_advection_restart.jl b/examples/p4est_3d_dgsem/elixir_advection_restart.jl index b3dead42399..9d19d81cf47 100644 --- a/examples/p4est_3d_dgsem/elixir_advection_restart.jl +++ b/examples/p4est_3d_dgsem/elixir_advection_restart.jl @@ -14,7 +14,7 @@ trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_basic.jl"), # Note: If you get a restart file from somewhere else, you need to provide # appropriate setups in the elixir loading a restart file -restart_filename = joinpath("out", "restart_000010.h5") +restart_filename = joinpath("out", "restart_000000010.h5") mesh = load_mesh(restart_filename) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, diff --git a/examples/structured_2d_dgsem/elixir_advection_restart.jl b/examples/structured_2d_dgsem/elixir_advection_restart.jl index 0accbdba702..896eda4f845 100644 --- a/examples/structured_2d_dgsem/elixir_advection_restart.jl +++ b/examples/structured_2d_dgsem/elixir_advection_restart.jl @@ -6,7 +6,7 @@ using Trixi # create a restart file elixir_file = "elixir_advection_extended.jl" -restart_file = "restart_000021.h5" +restart_file = "restart_000000021.h5" trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file)) diff --git a/examples/structured_3d_dgsem/elixir_advection_restart.jl b/examples/structured_3d_dgsem/elixir_advection_restart.jl index e516d794df8..e70a83cee6a 100644 --- a/examples/structured_3d_dgsem/elixir_advection_restart.jl +++ b/examples/structured_3d_dgsem/elixir_advection_restart.jl @@ -14,7 +14,7 @@ trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_basic.jl"), # Note: If you get a restart file from somewhere else, you need to provide # appropriate setups in the elixir loading a restart file -restart_filename = joinpath("out", "restart_000010.h5") +restart_filename = joinpath("out", "restart_000000010.h5") mesh = load_mesh(restart_filename) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, diff --git a/examples/tree_2d_dgsem/elixir_advection_restart.jl b/examples/tree_2d_dgsem/elixir_advection_restart.jl index 6052632ecad..c1a0cb8a8cc 100644 --- a/examples/tree_2d_dgsem/elixir_advection_restart.jl +++ b/examples/tree_2d_dgsem/elixir_advection_restart.jl @@ -15,7 +15,7 @@ trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_extended.jl"), a # Note: If you get a restart file from somewhere else, you need to provide # appropriate setups in the elixir loading a restart file -restart_filename = joinpath("out", "restart_000040.h5") +restart_filename = joinpath("out", "restart_000000040.h5") mesh = load_mesh(restart_filename) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/tree_2d_dgsem/elixir_advection_restart_amr.jl b/examples/tree_2d_dgsem/elixir_advection_restart_amr.jl index f366640ef51..96b62f4f4c8 100644 --- a/examples/tree_2d_dgsem/elixir_advection_restart_amr.jl +++ b/examples/tree_2d_dgsem/elixir_advection_restart_amr.jl @@ -15,7 +15,7 @@ trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_extended.jl"), a # Note: If you get a restart file from somewhere else, you need to provide # appropriate setups in the elixir loading a restart file -restart_filename = joinpath("out", "restart_000040.h5") +restart_filename = joinpath("out", "restart_000000040.h5") mesh = load_mesh(restart_filename) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/tree_3d_dgsem/elixir_advection_restart.jl b/examples/tree_3d_dgsem/elixir_advection_restart.jl index a53f4ebf5b1..0a3e97ece06 100644 --- a/examples/tree_3d_dgsem/elixir_advection_restart.jl +++ b/examples/tree_3d_dgsem/elixir_advection_restart.jl @@ -13,7 +13,7 @@ trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_extended.jl")) # Note: If you get a restart file from somewhere else, you need to provide # appropriate setups in the elixir loading a restart file -restart_filename = joinpath("out", "restart_000019.h5") +restart_filename = joinpath("out", "restart_000000019.h5") mesh = load_mesh(restart_filename) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl index c1908691902..4676718f72f 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_restart.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_restart.jl @@ -13,7 +13,7 @@ trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_euler_basic.jl")) # Note: If you get a restart file from somewhere else, you need to provide # appropriate setups in the elixir loading a restart file -restart_filename = joinpath("out", "restart_000050.h5") +restart_filename = joinpath("out", "restart_000000050.h5") mesh = load_mesh(restart_filename) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index 2c8497dc28d..860e3fa21d3 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -657,6 +657,15 @@ include("analysis_dg2d_parallel.jl") include("analysis_dg3d.jl") include("analysis_dg3d_parallel.jl") +# This version of `analyze` is used for [`AnalysisSurfaceIntegral`](@ref) which requires +# `semi` to be passed along to retrieve the current boundary indices, which are non-static +# in the case of AMR. +function analyze(quantity::AnalysisSurfaceIntegral, du, u, t, + semi::AbstractSemidiscretization) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + analyze(quantity, du, u, t, mesh, equations, solver, cache, semi) +end + # Special analyze for `SemidiscretizationHyperbolicParabolic` such that # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces. # Note that this needs to be included after `analysis_surface_integral_2d.jl` to @@ -669,6 +678,6 @@ function analyze(quantity::AnalysisSurfaceIntegral{Variable}, mesh, equations, solver, cache = mesh_equations_solver_cache(semi) equations_parabolic = semi.equations_parabolic cache_parabolic = semi.cache_parabolic - analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, + analyze(quantity, du, u, t, mesh, equations, equations_parabolic, solver, cache, semi, cache_parabolic) end diff --git a/src/callbacks_step/analysis_surface_integral_2d.jl b/src/callbacks_step/analysis_surface_integral_2d.jl index 7ae259e5285..d92ac5a2076 100644 --- a/src/callbacks_step/analysis_surface_integral_2d.jl +++ b/src/callbacks_step/analysis_surface_integral_2d.jl @@ -21,30 +21,18 @@ drag coefficient [`DragCoefficientPressure`](@ref) of e.g. an airfoil with the b name `:Airfoil` in 2D. - `semi::Semidiscretization`: Passed in to retrieve boundary condition information -- `boundary_symbol_or_boundary_symbols::Symbol|Vector{Symbol}`: Name(s) of the boundary/boundaries +- `boundary_symbols::NTuple{NBoundaries, Symbol}`: Name(s) of the boundary/boundaries where the quantity of interest is computed - `variable::Variable`: Quantity of interest, like lift or drag """ -struct AnalysisSurfaceIntegral{Variable} - indices::Vector{Int} # Indices in `boundary_condition_indices` where quantity of interest is computed +struct AnalysisSurfaceIntegral{Variable, NBoundaries} variable::Variable # Quantity of interest, like lift or drag + boundary_symbols::NTuple{NBoundaries, Symbol} # Name(s) of the boundary/boundaries - function AnalysisSurfaceIntegral(semi, boundary_symbol, variable) - @unpack boundary_symbol_indices = semi.boundary_conditions - indices = boundary_symbol_indices[boundary_symbol] - - return new{typeof(variable)}(indices, variable) - end - - function AnalysisSurfaceIntegral(semi, boundary_symbols::Vector{Symbol}, variable) - @unpack boundary_symbol_indices = semi.boundary_conditions - indices = Vector{Int}() - for name in boundary_symbols - append!(indices, boundary_symbol_indices[name]) - end - sort!(indices) - - return new{typeof(variable)}(indices, variable) + function AnalysisSurfaceIntegral(semi, + boundary_symbols::NTuple{NBoundaries, Symbol}, + variable) where {NBoundaries} + return new{typeof(variable), NBoundaries}(variable, boundary_symbols) end end @@ -255,14 +243,26 @@ function (drag_coefficient::DragCoefficientShearStress)(u, normal_direction, x, (0.5 * rhoinf * uinf^2 * linf) end +function get_boundary_indices(boundary_symbols, boundary_symbol_indices) + indices = Int[] + for name in boundary_symbols + append!(indices, boundary_symbol_indices[name]) + end + sort!(indices) # Try to achieve some data locality by sorting + + return indices +end + function analyze(surface_variable::AnalysisSurfaceIntegral, du, u, t, mesh::P4estMesh{2}, - equations, dg::DGSEM, cache) + equations, dg::DGSEM, cache, semi) @unpack boundaries = cache @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements @unpack weights = dg.basis - @unpack indices, variable = surface_variable + @unpack variable, boundary_symbols = surface_variable + @unpack boundary_symbol_indices = semi.boundary_conditions + indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices) surface_integral = zero(eltype(u)) index_range = eachnode(dg) @@ -308,13 +308,15 @@ end function analyze(surface_variable::AnalysisSurfaceIntegral{Variable}, du, u, t, mesh::P4estMesh{2}, equations, equations_parabolic, - dg::DGSEM, cache, + dg::DGSEM, cache, semi, cache_parabolic) where {Variable <: VariableViscous} @unpack boundaries = cache @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements @unpack weights = dg.basis - @unpack indices, variable = surface_variable + @unpack variable, boundary_symbols = surface_variable + @unpack boundary_symbol_indices = semi.boundary_conditions + indices = get_boundary_indices(boundary_symbols, boundary_symbol_indices) # Additions for parabolic @unpack viscous_container = cache_parabolic diff --git a/src/callbacks_step/save_restart_dg.jl b/src/callbacks_step/save_restart_dg.jl index b83402c5f86..29b3858c135 100644 --- a/src/callbacks_step/save_restart_dg.jl +++ b/src/callbacks_step/save_restart_dg.jl @@ -14,7 +14,7 @@ function save_restart_file(u, time, dt, timestep, @unpack output_directory = restart_callback # Filename based on current time step - filename = joinpath(output_directory, @sprintf("restart_%06d.h5", timestep)) + filename = joinpath(output_directory, @sprintf("restart_%09d.h5", timestep)) # Restart files always store conservative variables data = u @@ -93,7 +93,7 @@ function save_restart_file(u, time, dt, timestep, restart_callback) @unpack output_directory = restart_callback # Filename based on current time step - filename = joinpath(output_directory, @sprintf("restart_%06d.h5", timestep)) + filename = joinpath(output_directory, @sprintf("restart_%09d.h5", timestep)) if HDF5.has_parallel() save_restart_file_parallel(u, time, dt, timestep, mesh, equations, dg, cache, @@ -337,7 +337,7 @@ function save_adaptive_time_integrator(integrator, timestep = integrator.stats.naccept # Filename based on current time step - filename = joinpath(output_directory, @sprintf("restart_%06d.h5", timestep)) + filename = joinpath(output_directory, @sprintf("restart_%09d.h5", timestep)) # Open file (preserve existing content) h5open(filename, "r+") do file diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index deae8f7c930..57cdb3ef8a2 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -19,10 +19,10 @@ function save_solution_file(u, time, dt, timestep, # Filename based on current time step if isempty(system) - filename = joinpath(output_directory, @sprintf("solution_%06d.h5", timestep)) + filename = joinpath(output_directory, @sprintf("solution_%09d.h5", timestep)) else filename = joinpath(output_directory, - @sprintf("solution_%s_%06d.h5", system, timestep)) + @sprintf("solution_%s_%09d.h5", system, timestep)) end # Convert to different set of variables if requested @@ -101,10 +101,10 @@ function save_solution_file(u, time, dt, timestep, # Filename based on current time step if isempty(system) - filename = joinpath(output_directory, @sprintf("solution_%06d.h5", timestep)) + filename = joinpath(output_directory, @sprintf("solution_%09d.h5", timestep)) else filename = joinpath(output_directory, - @sprintf("solution_%s_%06d.h5", system, timestep)) + @sprintf("solution_%s_%09d.h5", system, timestep)) end # Convert to different set of variables if requested diff --git a/src/callbacks_step/visualization.jl b/src/callbacks_step/visualization.jl index 30ac88e9fd7..f91fe27bd33 100644 --- a/src/callbacks_step/visualization.jl +++ b/src/callbacks_step/visualization.jl @@ -255,7 +255,7 @@ function save_plot(plot_data, variable_names; Plots.plot(plots..., layout = layout) # Determine filename and save plot - filename = joinpath("out", @sprintf("solution_%06d.png", timestep)) + filename = joinpath("out", @sprintf("solution_%09d.png", timestep)) Plots.savefig(filename) end end # @muladd diff --git a/src/equations/compressible_navier_stokes_1d.jl b/src/equations/compressible_navier_stokes_1d.jl index 3dbdf11c56b..024ce2d7e87 100644 --- a/src/equations/compressible_navier_stokes_1d.jl +++ b/src/equations/compressible_navier_stokes_1d.jl @@ -171,7 +171,7 @@ function flux(u, gradients, orientation::Integer, mu = dynamic_viscosity(u, equations) # viscous flux components in the x-direction - f1 = zero(rho) + f1 = 0 f2 = tau_11 * mu f3 = (v1 * tau_11 + q1) * mu @@ -252,7 +252,7 @@ end @inline function temperature(u, equations::CompressibleNavierStokesDiffusion1D) rho, rho_v1, rho_e = u - p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1^2 / rho) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1^2 / rho) T = p / rho return T end diff --git a/src/equations/compressible_navier_stokes_2d.jl b/src/equations/compressible_navier_stokes_2d.jl index 3256343703a..5708abb4e00 100644 --- a/src/equations/compressible_navier_stokes_2d.jl +++ b/src/equations/compressible_navier_stokes_2d.jl @@ -157,13 +157,13 @@ function flux(u, gradients, orientation::Integer, # Components of viscous stress tensor - # (4/3 * (v1)_x - 2/3 * (v2)_y) - tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy + # (4 * (v1)_x / 3 - 2 * (v2)_y / 3) + tau_11 = 4 * dv1dx / 3 - 2 * dv2dy / 3 # ((v1)_y + (v2)_x) # stress tensor is symmetric tau_12 = dv1dy + dv2dx # = tau_21 # (4/3 * (v2)_y - 2/3 * (v1)_x) - tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx + tau_22 = 4 * dv2dy / 3 - 2 * dv1dx / 3 # Fick's law q = -kappa * grad(T) = -kappa * grad(p / (R rho)) # with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr) @@ -181,7 +181,7 @@ function flux(u, gradients, orientation::Integer, if orientation == 1 # viscous flux components in the x-direction - f1 = zero(rho) + f1 = 0 f2 = tau_11 * mu f3 = tau_12 * mu f4 = (v1 * tau_11 + v2 * tau_12 + q1) * mu @@ -190,7 +190,7 @@ function flux(u, gradients, orientation::Integer, else # if orientation == 2 # viscous flux components in the y-direction # Note, symmetry is exploited for tau_12 = tau_21 - g1 = zero(rho) + g1 = 0 g2 = tau_12 * mu # tau_21 * mu g3 = tau_22 * mu g4 = (v1 * tau_12 + v2 * tau_22 + q2) * mu @@ -276,7 +276,7 @@ end @inline function temperature(u, equations::CompressibleNavierStokesDiffusion2D) rho, rho_v1, rho_v2, rho_e = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho) T = p / rho return T end @@ -285,7 +285,7 @@ end # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v omega = vorticity(u, gradients, equations) - return 0.5 * u[1] * omega^2 + return 0.5f0 * u[1] * omega^2 end @inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion2D) diff --git a/src/equations/compressible_navier_stokes_3d.jl b/src/equations/compressible_navier_stokes_3d.jl index 9833122eb32..9622882fc1a 100644 --- a/src/equations/compressible_navier_stokes_3d.jl +++ b/src/equations/compressible_navier_stokes_3d.jl @@ -162,12 +162,12 @@ function flux(u, gradients, orientation::Integer, # Components of viscous stress tensor # Diagonal parts - # (4/3 * (v1)_x - 2/3 * ((v2)_y + (v3)_z) - tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * (dv2dy + dv3dz) - # (4/3 * (v2)_y - 2/3 * ((v1)_x + (v3)_z) - tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * (dv1dx + dv3dz) - # (4/3 * (v3)_z - 2/3 * ((v1)_x + (v2)_y) - tau_33 = 4.0 / 3.0 * dv3dz - 2.0 / 3.0 * (dv1dx + dv2dy) + # (4 * (v1)_x / 3 - 2 * ((v2)_y + (v3)_z)) / 3) + tau_11 = 4 * dv1dx / 3 - 2 * (dv2dy + dv3dz) / 3 + # (4 * (v2)_y / 3 - 2 * ((v1)_x + (v3)_z) / 3) + tau_22 = 4 * dv2dy / 3 - 2 * (dv1dx + dv3dz) / 3 + # (4 * (v3)_z / 3 - 2 * ((v1)_x + (v2)_y) / 3) + tau_33 = 4 * dv3dz / 3 - 2 * (dv1dx + dv2dy) / 3 # Off diagonal parts, exploit that stress tensor is symmetric # ((v1)_y + (v2)_x) @@ -194,7 +194,7 @@ function flux(u, gradients, orientation::Integer, if orientation == 1 # viscous flux components in the x-direction - f1 = zero(rho) + f1 = 0 f2 = tau_11 * mu f3 = tau_12 * mu f4 = tau_13 * mu @@ -204,7 +204,7 @@ function flux(u, gradients, orientation::Integer, elseif orientation == 2 # viscous flux components in the y-direction # Note, symmetry is exploited for tau_12 = tau_21 - g1 = zero(rho) + g1 = 0 g2 = tau_12 * mu # tau_21 * mu g3 = tau_22 * mu g4 = tau_23 * mu @@ -214,7 +214,7 @@ function flux(u, gradients, orientation::Integer, else # if orientation == 3 # viscous flux components in the z-direction # Note, symmetry is exploited for tau_13 = tau_31, tau_23 = tau_32 - h1 = zero(rho) + h1 = 0 h2 = tau_13 * mu # tau_31 * mu h3 = tau_23 * mu # tau_32 * mu h4 = tau_33 * mu @@ -304,7 +304,7 @@ end @inline function temperature(u, equations::CompressibleNavierStokesDiffusion3D) rho, rho_v1, rho_v2, rho_v3, rho_e = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho) T = p / rho return T end @@ -313,7 +313,7 @@ end # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v omega = vorticity(u, gradients, equations) - return 0.5 * u[1] * (omega[1]^2 + omega[2]^2 + omega[3]^2) + return 0.5f0 * u[1] * (omega[1]^2 + omega[2]^2 + omega[3]^2) end @inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion3D) diff --git a/src/equations/hyperbolic_diffusion_1d.jl b/src/equations/hyperbolic_diffusion_1d.jl index 39a555e7c72..2b774470fed 100644 --- a/src/equations/hyperbolic_diffusion_1d.jl +++ b/src/equations/hyperbolic_diffusion_1d.jl @@ -53,13 +53,14 @@ function initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations1D) # elliptic equation: -νΔϕ = f # Taken from Section 6.1 of Nishikawa https://doi.org/10.1016/j.jcp.2007.07.029 - if t == 0.0 + RealT = eltype(x) + if t == 0 # initial "guess" of the solution and its derivative phi = x[1]^2 - x[1] q1 = 2 * x[1] - 1 else phi = sinpi(x[1]) # ϕ - q1 = pi * cospi(x[1]) # ϕ_x + q1 = convert(RealT, pi) * cospi(x[1]) # ϕ_x end return SVector(phi, q1) end @@ -76,9 +77,10 @@ diffusion system that is used with [`initial_condition_poisson_nonperiodic`](@re equations::HyperbolicDiffusionEquations1D) # elliptic equation: -νΔϕ = f # analytical solution: ϕ = sin(πx) and f = π^2sin(πx) + RealT = eltype(u) @unpack inv_Tr = equations - dphi = pi^2 * sinpi(x[1]) + dphi = convert(RealT, pi)^2 * sinpi(x[1]) dq1 = -inv_Tr * u[2] return SVector(dphi, dq1) @@ -96,8 +98,9 @@ function boundary_condition_poisson_nonperiodic(u_inner, orientation, direction, surface_flux_function, equations::HyperbolicDiffusionEquations1D) # elliptic equation: -νΔϕ = f + RealT = eltype(u_inner) phi = sinpi(x[1]) # ϕ - q1 = pi * cospi(x[1]) # ϕ_x + q1 = convert(RealT, pi) * cospi(x[1]) # ϕ_x u_boundary = SVector(phi, q1) # Calculate boundary flux @@ -122,7 +125,7 @@ Source term that only includes the forcing from the hyperbolic diffusion system. dq1 = -inv_Tr * u[2] - return SVector(zero(dq1), dq1) + return SVector(0, dq1) end """ @@ -138,13 +141,14 @@ function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations1D) # Determine phi_x - G = 1.0 # gravitational constant - C = -4.0 * G / pi # -4 * G / ndims * pi - A = 0.1 # perturbation coefficient must match Euler setup + RealT = eltype(x) + G = 1 # gravitational constant + C = -4 * G / convert(RealT, pi) # -4 * G / ndims * pi + A = convert(RealT, 0.1) # perturbation coefficient must match Euler setup rho1 = A * sinpi(x[1] - t) # initialize with ansatz of gravity potential phi = C * rho1 - q1 = C * A * pi * cospi(x[1] - t) # = gravity acceleration in x-direction + q1 = C * A * convert(RealT, pi) * cospi(x[1] - t) # = gravity acceleration in x-direction return SVector(phi, q1) end @@ -196,6 +200,6 @@ end @inline function energy_total(u, equations::HyperbolicDiffusionEquations1D) # energy function as found in equations (2.5.12) in the book "I Do Like CFD, Vol. 1" phi, q1 = u - return 0.5 * (phi^2 + equations.Lr^2 * q1^2) + return 0.5f0 * (phi^2 + equations.Lr^2 * q1^2) end end # @muladd diff --git a/src/equations/hyperbolic_diffusion_2d.jl b/src/equations/hyperbolic_diffusion_2d.jl index 511d1b8935d..28f0c2bf47f 100644 --- a/src/equations/hyperbolic_diffusion_2d.jl +++ b/src/equations/hyperbolic_diffusion_2d.jl @@ -39,17 +39,18 @@ end @inline function initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations2D) # elliptic equation: -ν Δϕ = f in Ω, u = g on ∂Ω + RealT = eltype(x) if iszero(t) - T = eltype(x) - phi = one(T) - q1 = one(T) - q2 = one(T) + phi = one(RealT) + q1 = one(RealT) + q2 = one(RealT) else - sinpi_x1, cospi_x1 = sincos(pi * x[1]) - sinpi_2x2, cospi_2x2 = sincos(pi * 2 * x[2]) + # TODO: sincospi + sinpi_x1, cospi_x1 = sincos(convert(RealT, pi) * x[1]) + sinpi_2x2, cospi_2x2 = sincos(convert(RealT, pi) * 2 * x[2]) phi = 2 * cospi_x1 * sinpi_2x2 + 2 # ϕ - q1 = -2 * pi * sinpi_x1 * sinpi_2x2 # ϕ_x - q2 = 4 * pi * cospi_x1 * cospi_2x2 # ϕ_y + q1 = -2 * convert(RealT, pi) * sinpi_x1 * sinpi_2x2 # ϕ_x + q2 = 4 * convert(RealT, pi) * cospi_x1 * cospi_2x2 # ϕ_y end return SVector(phi, q1, q2) end @@ -58,10 +59,11 @@ end equations::HyperbolicDiffusionEquations2D) # elliptic equation: -ν Δϕ = f in Ω, u = g on ∂Ω # analytical solution: ϕ = 2cos(πx)sin(2πy) + 2 and f = 10π^2cos(πx)sin(2πy) + RealT = eltype(u) @unpack inv_Tr = equations x1, x2 = x - du1 = 10 * pi^2 * cospi(x1) * sinpi(2 * x2) + du1 = 10 * convert(RealT, pi)^2 * cospi(x1) * sinpi(2 * x2) du2 = -inv_Tr * u[2] du3 = -inv_Tr * u[3] @@ -115,13 +117,14 @@ function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations2D) # Determine phi_x, phi_y - G = 1.0 # gravitational constant - C = -2.0 * G / pi - A = 0.1 # perturbation coefficient must match Euler setup - rho1 = A * sin(pi * (x[1] + x[2] - t)) + RealT = eltype(x) + G = 1 # gravitational constant + C = -2 * G / convert(RealT, pi) + A = convert(RealT, 0.1) # perturbation coefficient must match Euler setup + rho1 = A * sinpi(x[1] + x[2] - t) # initialize with ansatz of gravity potential phi = C * rho1 - q1 = C * A * pi * cos(pi * (x[1] + x[2] - t)) # = gravity acceleration in x-direction + q1 = C * A * convert(RealT, pi) * cospi(x[1] + x[2] - t) # = gravity acceleration in x-direction q2 = q1 # = gravity acceleration in y-direction return SVector(phi, q1, q2) @@ -133,13 +136,14 @@ end phi, q1, q2 = u @unpack inv_Tr = equations + RealT = eltype(u) if orientation == 1 f1 = -equations.nu * q1 f2 = -phi * inv_Tr - f3 = zero(phi) + f3 = zero(RealT) else f1 = -equations.nu * q2 - f2 = zero(phi) + f2 = zero(RealT) f3 = -phi * inv_Tr end @@ -181,13 +185,13 @@ end # this is an optimized version of the application of the upwind dissipation matrix: # dissipation = 0.5*R_n*|Λ|*inv(R_n)[[u]] λ_max = sqrt(equations.nu * equations.inv_Tr) - f1 = 1 / 2 * (f_ll[1] + f_rr[1]) - 1 / 2 * λ_max * (phi_rr - phi_ll) + f1 = 0.5f0 * (f_ll[1] + f_rr[1]) - 0.5f0 * λ_max * (phi_rr - phi_ll) if orientation == 1 # x-direction - f2 = 1 / 2 * (f_ll[2] + f_rr[2]) - 1 / 2 * λ_max * (q1_rr - q1_ll) - f3 = 1 / 2 * (f_ll[3] + f_rr[3]) + f2 = 0.5f0 * (f_ll[2] + f_rr[2]) - 0.5f0 * λ_max * (q1_rr - q1_ll) + f3 = 0.5f0 * (f_ll[3] + f_rr[3]) else # y-direction - f2 = 1 / 2 * (f_ll[2] + f_rr[2]) - f3 = 1 / 2 * (f_ll[3] + f_rr[3]) - 1 / 2 * λ_max * (q2_rr - q2_ll) + f2 = 0.5f0 * (f_ll[2] + f_rr[2]) + f3 = 0.5f0 * (f_ll[3] + f_rr[3]) - 0.5f0 * λ_max * (q2_rr - q2_ll) end return SVector(f1, f2, f3) @@ -204,13 +208,13 @@ end # this is an optimized version of the application of the upwind dissipation matrix: # dissipation = 0.5*R_n*|Λ|*inv(R_n)[[u]] λ_max = sqrt(equations.nu * equations.inv_Tr) - f1 = 1 / 2 * (f_ll[1] + f_rr[1]) - - 1 / 2 * λ_max * (phi_rr - phi_ll) * + f1 = 0.5f0 * (f_ll[1] + f_rr[1]) - + 0.5f0 * λ_max * (phi_rr - phi_ll) * sqrt(normal_direction[1]^2 + normal_direction[2]^2) - f2 = 1 / 2 * (f_ll[2] + f_rr[2]) - - 1 / 2 * λ_max * (q1_rr - q1_ll) * normal_direction[1] - f3 = 1 / 2 * (f_ll[3] + f_rr[3]) - - 1 / 2 * λ_max * (q2_rr - q2_ll) * normal_direction[2] + f2 = 0.5f0 * (f_ll[2] + f_rr[2]) - + 0.5f0 * λ_max * (q1_rr - q1_ll) * normal_direction[1] + f3 = 0.5f0 * (f_ll[3] + f_rr[3]) - + 0.5f0 * λ_max * (q2_rr - q2_ll) * normal_direction[2] return SVector(f1, f2, f3) end @@ -244,6 +248,6 @@ end @inline function energy_total(u, equations::HyperbolicDiffusionEquations2D) # energy function as found in equations (2.5.12) in the book "I Do Like CFD, Vol. 1" phi, q1, q2 = u - return 0.5 * (phi^2 + equations.Lr^2 * (q1^2 + q2^2)) + return 0.5f0 * (phi^2 + equations.Lr^2 * (q1^2 + q2^2)) end end # @muladd diff --git a/src/equations/hyperbolic_diffusion_3d.jl b/src/equations/hyperbolic_diffusion_3d.jl index ed807511b67..e9017b1544e 100644 --- a/src/equations/hyperbolic_diffusion_3d.jl +++ b/src/equations/hyperbolic_diffusion_3d.jl @@ -49,16 +49,22 @@ end function initial_condition_poisson_nonperiodic(x, t, equations::HyperbolicDiffusionEquations3D) # elliptic equation: -νΔϕ = f - if t == 0.0 - phi = 1.0 - q1 = 1.0 - q2 = 1.0 - q3 = 1.0 + RealT = eltype(x) + if t == 0 + phi = one(RealT) + q1 = one(RealT) + q2 = one(RealT) + q3 = one(RealT) else - phi = 2.0 * cos(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) + 2.0 # ϕ - q1 = -2.0 * pi * sin(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) # ϕ_x - q2 = 4.0 * pi * cos(pi * x[1]) * cos(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) # ϕ_y - q3 = 4.0 * pi * cos(pi * x[1]) * sin(2.0 * pi * x[2]) * cos(2.0 * pi * x[3]) # ϕ_z + phi = 2 * cospi(x[1]) * + sinpi(2 * x[2]) * + sinpi(2 * x[3]) + 2 # ϕ + q1 = -2 * convert(RealT, pi) * sinpi(x[1]) * + sinpi(2 * x[2]) * sinpi(2 * x[3]) # ϕ_x + q2 = 4 * convert(RealT, pi) * cospi(x[1]) * + cospi(2 * x[2]) * sinpi(2 * x[3]) # ϕ_y + q3 = 4 * convert(RealT, pi) * cospi(x[1]) * + sinpi(2 * x[2]) * cospi(2 * x[3]) # ϕ_z end return SVector(phi, q1, q2, q3) end @@ -67,10 +73,11 @@ end equations::HyperbolicDiffusionEquations3D) # elliptic equation: -νΔϕ = f # analytical solution: ϕ = 2 cos(πx)sin(2πy)sin(2πz) + 2 and f = 18 π^2 cos(πx)sin(2πy)sin(2πz) + RealT = eltype(u) @unpack inv_Tr = equations x1, x2, x3 = x - du1 = 18 * pi^2 * cospi(x1) * sinpi(2 * x2) * sinpi(2 * x3) + du1 = 18 * convert(RealT, pi)^2 * cospi(x1) * sinpi(2 * x2) * sinpi(2 * x3) du2 = -inv_Tr * u[2] du3 = -inv_Tr * u[3] du4 = -inv_Tr * u[4] @@ -82,10 +89,15 @@ function boundary_condition_poisson_nonperiodic(u_inner, orientation, direction, surface_flux_function, equations::HyperbolicDiffusionEquations3D) # elliptic equation: -νΔϕ = f - phi = 2.0 * cos(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) + 2.0 # ϕ - q1 = -2.0 * pi * sin(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) # ϕ_x - q2 = 4.0 * pi * cos(pi * x[1]) * cos(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) # ϕ_y - q3 = 4.0 * pi * cos(pi * x[1]) * sin(2.0 * pi * x[2]) * cos(2.0 * pi * x[3]) # ϕ_z + RealT = eltype(u_inner) + phi = 2 * cospi(x[1]) * sinpi(2 * x[2]) * + sinpi(2 * x[3]) + 2 # ϕ + q1 = -2 * convert(RealT, pi) * sinpi(x[1]) * + sinpi(2 * x[2]) * sinpi(2 * x[3]) # ϕ_x + q2 = 4 * convert(RealT, pi) * cospi(x[1]) * + cospi(2 * x[2]) * sinpi(2 * x[3]) # ϕ_y + q3 = 4 * convert(RealT, pi) * cospi(x[1]) * + sinpi(2 * x[2]) * cospi(2 * x[3]) # ϕ_z u_boundary = SVector(phi, q1, q2, q3) # Calculate boundary flux @@ -108,7 +120,7 @@ Source term that only includes the forcing from the hyperbolic diffusion system. # harmonic solution ϕ = (sinh(πx)sin(πy) + sinh(πy)sin(πx))/sinh(π), so f = 0 @unpack inv_Tr = equations - du1 = zero(u[1]) + du1 = 0 du2 = -inv_Tr * u[2] du3 = -inv_Tr * u[3] du4 = -inv_Tr * u[4] @@ -129,13 +141,15 @@ function initial_condition_eoc_test_coupled_euler_gravity(x, t, equations::HyperbolicDiffusionEquations3D) # Determine phi_x, phi_y - G = 1.0 # gravitational constant - C_grav = -4 * G / (3 * pi) # "3" is the number of spatial dimensions # 2D: -2.0*G/pi - A = 0.1 # perturbation coefficient must match Euler setup - rho1 = A * sin(pi * (x[1] + x[2] + x[3] - t)) + RealT = eltype(x) + G = 1 # gravitational constant + C_grav = -4 * G / (3 * convert(RealT, pi)) # "3" is the number of spatial dimensions # 2D: -2.0*G/pi + A = convert(RealT, 0.1) # perturbation coefficient must match Euler setup + rho1 = A * sinpi(x[1] + x[2] + x[3] - t) # initialize with ansatz of gravity potential phi = C_grav * rho1 - q1 = C_grav * A * pi * cos(pi * (x[1] + x[2] + x[3] - t)) # = gravity acceleration in x-direction + q1 = C_grav * A * convert(RealT, pi) * + cospi(x[1] + x[2] + x[3] - t) # = gravity acceleration in x-direction q2 = q1 # = gravity acceleration in y-direction q3 = q1 # = gravity acceleration in z-direction @@ -147,20 +161,21 @@ end equations::HyperbolicDiffusionEquations3D) phi, q1, q2, q3 = u + RealT = eltype(u) if orientation == 1 f1 = -equations.nu * q1 f2 = -phi * equations.inv_Tr - f3 = zero(phi) - f4 = zero(phi) + f3 = zero(RealT) + f4 = zero(RealT) elseif orientation == 2 f1 = -equations.nu * q2 - f2 = zero(phi) + f2 = zero(RealT) f3 = -phi * equations.inv_Tr - f4 = zero(phi) + f4 = zero(RealT) else f1 = -equations.nu * q3 - f2 = zero(phi) - f3 = zero(phi) + f2 = zero(RealT) + f3 = zero(RealT) f4 = -phi * equations.inv_Tr end @@ -184,19 +199,19 @@ end # this is an optimized version of the application of the upwind dissipation matrix: # dissipation = 0.5*R_n*|Λ|*inv(R_n)[[u]] λ_max = sqrt(equations.nu * equations.inv_Tr) - f1 = 1 / 2 * (f_ll[1] + f_rr[1]) - 1 / 2 * λ_max * (phi_rr - phi_ll) + f1 = 0.5f0 * (f_ll[1] + f_rr[1]) - 0.5f0 * λ_max * (phi_rr - phi_ll) if orientation == 1 # x-direction - f2 = 1 / 2 * (f_ll[2] + f_rr[2]) - 1 / 2 * λ_max * (q1_rr - q1_ll) - f3 = 1 / 2 * (f_ll[3] + f_rr[3]) - f4 = 1 / 2 * (f_ll[4] + f_rr[4]) + f2 = 0.5f0 * (f_ll[2] + f_rr[2]) - 0.5f0 * λ_max * (q1_rr - q1_ll) + f3 = 0.5f0 * (f_ll[3] + f_rr[3]) + f4 = 0.5f0 * (f_ll[4] + f_rr[4]) elseif orientation == 2 # y-direction - f2 = 1 / 2 * (f_ll[2] + f_rr[2]) - f3 = 1 / 2 * (f_ll[3] + f_rr[3]) - 1 / 2 * λ_max * (q2_rr - q2_ll) - f4 = 1 / 2 * (f_ll[4] + f_rr[4]) + f2 = 0.5f0 * (f_ll[2] + f_rr[2]) + f3 = 0.5f0 * (f_ll[3] + f_rr[3]) - 0.5f0 * λ_max * (q2_rr - q2_ll) + f4 = 0.5f0 * (f_ll[4] + f_rr[4]) else # y-direction - f2 = 1 / 2 * (f_ll[2] + f_rr[2]) - f3 = 1 / 2 * (f_ll[3] + f_rr[3]) - f4 = 1 / 2 * (f_ll[4] + f_rr[4]) - 1 / 2 * λ_max * (q3_rr - q3_ll) + f2 = 0.5f0 * (f_ll[2] + f_rr[2]) + f3 = 0.5f0 * (f_ll[3] + f_rr[3]) + f4 = 0.5f0 * (f_ll[4] + f_rr[4]) - 0.5f0 * λ_max * (q3_rr - q3_ll) end return SVector(f1, f2, f3, f4) @@ -232,6 +247,6 @@ end @inline function energy_total(u, equations::HyperbolicDiffusionEquations3D) # energy function as found in equation (2.5.12) in the book "I Do Like CFD, Vol. 1" phi, q1, q2, q3 = u - return 0.5 * (phi^2 + equations.Lr^2 * (q1^2 + q2^2 + q3^2)) + return 0.5f0 * (phi^2 + equations.Lr^2 * (q1^2 + q2^2 + q3^2)) end end # @muladd diff --git a/src/equations/ideal_glm_mhd_1d.jl b/src/equations/ideal_glm_mhd_1d.jl index 5a523daf3f6..3253e4410f8 100644 --- a/src/equations/ideal_glm_mhd_1d.jl +++ b/src/equations/ideal_glm_mhd_1d.jl @@ -42,14 +42,15 @@ end A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations1D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 - rho_v3 = -0.5 - rho_e = 50.0 - B1 = 3.0 - B2 = -1.2 - B3 = 0.5 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = -convert(RealT, 0.2) + rho_v3 = -0.5f0 + rho_e = 50 + B1 = 3 + B2 = -convert(RealT, 1.2) + B3 = 0.5f0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) end @@ -61,14 +62,15 @@ An Alfvén wave as smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations1D) # smooth Alfvén wave test from Derigs et al. FLASH (2016) # domain must be set to [0, 1], γ = 5/3 - rho = 1.0 - v1 = 0.0 + RealT = eltype(x) + rho = 1 + v1 = 0 # TODO: sincospi - si, co = sincos(2 * pi * x[1]) - v2 = 0.1 * si - v3 = 0.1 * co - p = 0.1 - B1 = 1.0 + si, co = sincos(2 * convert(RealT, pi) * x[1]) + v2 = convert(RealT, 0.1) * si + v3 = convert(RealT, 0.1) * co + p = convert(RealT, 0.1) + B1 = 1 B2 = v2 B3 = v3 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations) @@ -86,17 +88,18 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates + RealT = eltype(x) inicenter = (0,) x_norm = x[1] - inicenter[1] r = sqrt(x_norm^2) phi = atan(x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - return prim2cons(SVector(rho, v1, 0.0, 0.0, p, 1.0, 1.0, 1.0, 0.0), equations) + return prim2cons(SVector(rho, v1, 0, 0, p, 1, 1, 1, 0), equations) end # Calculate 1D flux in for a single point @@ -105,8 +108,8 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) p_over_gamma_minus_one = (rho_e - kin_en - mag_en) p = (equations.gamma - 1) * p_over_gamma_minus_one @@ -117,7 +120,7 @@ end f4 = rho_v1 * v3 - B1 * B3 f5 = (kin_en + equations.gamma * p_over_gamma_minus_one + 2 * mag_en) * v1 - B1 * (v1 * B1 + v2 * B2 + v3 * B3) - f6 = 0.0 + f6 = 0 f7 = v1 * B2 - v2 * B1 f8 = v1 * B3 - v3 * B1 @@ -150,44 +153,44 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 p_ll = (equations.gamma - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll) + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll) p_rr = (equations.gamma - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr # for convenience store v⋅B vel_dot_mag_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_mean = 0.5 * rho_avg / beta_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) # Ignore orientation since it is always "1" in 1D f1 = rho_mean * v1_avg - f2 = f1 * v1_avg + p_mean + 0.5 * mag_norm_avg - B1_avg * B1_avg + f2 = f1 * v1_avg + p_mean + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f3 = f1 * v2_avg - B1_avg * B2_avg f4 = f1 * v3_avg - B1_avg * B3_avg - f6 = 0.0 + f6 = 0 f7 = v1_avg * B2_avg - v2_avg * B1_avg f8 = v1_avg * B3_avg - v3_avg * B1_avg # total energy flux is complicated and involves the previous eight components - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + - f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5 * v1_mag_avg + + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg) return SVector(f1, f2, f3, f4, f5, f6, f7, f8) @@ -228,27 +231,27 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on orientation with specific direction averages f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below - f6 = 0.0 - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f6 = 0 + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -273,8 +276,8 @@ function flux_hllc(u_ll, u_rr, orientation::Integer, rho_rr, v1_rr, v2_rr, v3_rr, p_rr, B1_rr, B2_rr, B3_rr = cons2prim(u_rr, equations) # Total pressure, i.e., thermal + magnetic pressures (eq. (12)) - p_tot_ll = p_ll + 0.5 * (B1_ll^2 + B2_ll^2 + B3_ll^2) - p_tot_rr = p_rr + 0.5 * (B1_rr^2 + B2_rr^2 + B3_rr^2) + p_tot_ll = p_ll + 0.5f0 * (B1_ll^2 + B2_ll^2 + B3_ll^2) + p_tot_rr = p_rr + 0.5f0 * (B1_rr^2 + B2_rr^2 + B3_rr^2) # Conserved variables rho_v1_ll = u_ll[2] @@ -502,8 +505,8 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * (rho_e - - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 - + B1 * B1 + B2 * B2 + B3 * B3)) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + + B1 * B1 + B2 * B2 + B3 * B3)) return SVector(rho, v1, v2, v3, p, B1, B2, B3) end @@ -517,11 +520,11 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2)) + (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2)) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) / (equations.gamma - 1) - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) / (equations.gamma - 1) - 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = rho_p * v3 @@ -541,8 +544,8 @@ end rho_v2 = rho * v2 rho_v3 = rho * v3 rho_e = p / (equations.gamma - 1) + - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3) end @@ -554,17 +557,17 @@ end @inline function pressure(u, equations::IdealGlmMhdEquations1D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2)) + 0.5f0 * (B1^2 + B2^2 + B3^2)) return p end @inline function density_pressure(u, equations::IdealGlmMhdEquations1D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2)) + 0.5f0 * (B1^2 + B2^2 + B3^2)) return rho * p end @@ -576,7 +579,7 @@ end v3 = rho_v3 / rho v_mag = sqrt(v1^2 + v2^2 + v3^2) p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2)) + (rho_e - 0.5f0 * rho * v_mag^2 - 0.5f0 * (B1^2 + B2^2 + B3^2)) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -584,8 +587,8 @@ end b3 = B3 / sqrt_rho b_square = b1^2 + b2^2 + b3^2 - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) return c_f end @@ -611,7 +614,7 @@ as given by vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 p_ll = (equations.gamma - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll) + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr @@ -619,17 +622,17 @@ as given by vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 p_rr = (equations.gamma - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr) + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -645,19 +648,19 @@ as given by H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equations (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum - a_square_roe = ((2.0 - equations.gamma) * X + - (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity # Ignore orientation since it is always "1" in 1D c_a_roe = B1_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed - a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4 * a_square_roe * c_a_roe) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) return v1_roe, c_f_roe end @@ -666,9 +669,9 @@ end @inline function entropy_thermodynamic(cons, equations::IdealGlmMhdEquations1D) # Pressure p = (equations.gamma - 1) * - (cons[5] - 1 / 2 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + (cons[5] - 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] - - 1 / 2 * (cons[6]^2 + cons[7]^2 + cons[8]^2)) + 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2)) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) @@ -691,13 +694,13 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(cons, equations::IdealGlmMhdEquations1D) - return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + return 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] end # Calculate the magnetic energy for a conservative state `cons'. # OBS! For non-dinmensional form of the ideal MHD magnetic pressure ≡ magnetic energy @inline function energy_magnetic(cons, ::IdealGlmMhdEquations1D) - return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + return 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) end # Calculate internal energy for a conservative state `cons` diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 622410de855..ab2a4b066a1 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -48,15 +48,16 @@ end A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 - rho_v3 = -0.5 - rho_e = 50.0 - B1 = 3.0 - B2 = -1.2 - B3 = 0.5 - psi = 0.0 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = -convert(RealT, 0.2) + rho_v3 = -0.5f0 + rho_e = 50 + B1 = 3 + B2 = -convert(RealT, 1.2) + B3 = 0.5f0 + psi = 0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -68,18 +69,19 @@ An Alfvén wave as smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations2D) # smooth Alfvén wave test from Derigs et al. FLASH (2016) # domain must be set to [0, 1/cos(α)] x [0, 1/sin(α)], γ = 5/3 - alpha = 0.25 * pi + RealT = eltype(x) + alpha = 0.25f0 * convert(RealT, pi) x_perp = x[1] * cos(alpha) + x[2] * sin(alpha) - B_perp = 0.1 * sin(2.0 * pi * x_perp) - rho = 1.0 + B_perp = convert(RealT, 0.1) * sinpi(2 * x_perp) + rho = 1 v1 = -B_perp * sin(alpha) v2 = B_perp * cos(alpha) - v3 = 0.1 * cos(2.0 * pi * x_perp) - p = 0.1 + v3 = convert(RealT, 0.1) * cospi(2 * x_perp) + p = convert(RealT, 0.1) B1 = cos(alpha) + v1 B2 = sin(alpha) + v2 B3 = v3 - psi = 0.0 + psi = 0 return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) end @@ -95,6 +97,7 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates + RealT = eltype(x) inicenter = (0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] @@ -102,12 +105,12 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations phi = atan(y_norm, x_norm) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) - v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi) - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - return prim2cons(SVector(rho, v1, v2, 0.0, p, 1.0, 1.0, 1.0, 0.0), equations) + return prim2cons(SVector(rho, v1, v2, 0, p, 1, 1, 1, 0), equations) end # Pre-defined source terms should be implemented as @@ -119,9 +122,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one if orientation == 1 f1 = rho_v1 @@ -158,9 +161,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] @@ -495,36 +498,38 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 p_ll = (equations.gamma - 1) * - (rho_e_ll - 0.5 * rho_ll * vel_norm_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll - + 0.5f0 * psi_ll^2) p_rr = (equations.gamma - 1) * - (rho_e_rr - 0.5 * rho_rr * vel_norm_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr - + 0.5f0 * psi_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr # for convenience store v⋅B vel_dot_mag_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_mean = 0.5 * rho_avg / beta_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg - f2 = f1 * v1_avg + p_mean + 0.5 * mag_norm_avg - B1_avg * B1_avg + f2 = f1 * v1_avg + p_mean + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f3 = f1 * v2_avg - B1_avg * B2_avg f4 = f1 * v3_avg - B1_avg * B3_avg f6 = equations.c_h * psi_avg @@ -532,29 +537,29 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f8 = v1_avg * B3_avg - v3_avg * B1_avg f9 = equations.c_h * B1_avg # total energy flux is complicated and involves the previous eight components - psi_B1_avg = 0.5 * (B1_ll * psi_ll + B1_rr * psi_rr) - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B1_avg = 0.5f0 * (B1_ll * psi_ll + B1_rr * psi_rr) + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v1_mag_avg + + 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg - equations.c_h * psi_B1_avg) else f1 = rho_mean * v2_avg f2 = f1 * v1_avg - B1_avg * B2_avg - f3 = f1 * v2_avg + p_mean + 0.5 * mag_norm_avg - B2_avg * B2_avg + f3 = f1 * v2_avg + p_mean + 0.5f0 * mag_norm_avg - B2_avg * B2_avg f4 = f1 * v3_avg - B2_avg * B3_avg f6 = v2_avg * B1_avg - v1_avg * B2_avg f7 = equations.c_h * psi_avg f8 = v2_avg * B3_avg - v3_avg * B2_avg f9 = equations.c_h * B2_avg # total energy flux is complicated and involves the previous eight components - psi_B2_avg = 0.5 * (B2_ll * psi_ll + B2_rr * psi_rr) - v2_mag_avg = 0.5 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr) + v2_mag_avg = 0.5f0 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v2_mag_avg + + 0.5f0 * v2_mag_avg + B2_avg * vel_dot_mag_avg - equations.c_h * psi_B2_avg) end @@ -598,31 +603,31 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below f6 = equations.c_h * psi_avg - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) - f9 = equations.c_h * 0.5 * (B1_ll + B1_rr) + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f9 = equations.c_h * 0.5f0 * (B1_ll + B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -633,20 +638,20 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat equations.c_h * (B1_ll * psi_rr + B1_rr * psi_ll))) else # orientation == 2 f1 = rho_mean * v2_avg - f2 = f1 * v1_avg - 0.5 * (B2_ll * B1_rr + B2_rr * B1_ll) + f2 = f1 * v1_avg - 0.5f0 * (B2_ll * B1_rr + B2_rr * B1_ll) f3 = f1 * v2_avg + p_avg + magnetic_square_avg - - 0.5 * (B2_ll * B2_rr + B2_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B2_ll * B3_rr + B2_rr * B3_ll) + 0.5f0 * (B2_ll * B2_rr + B2_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B2_ll * B3_rr + B2_rr * B3_ll) #f5 below - f6 = 0.5 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) + f6 = 0.5f0 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) f7 = equations.c_h * psi_avg - f8 = 0.5 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) - f9 = equations.c_h * 0.5 * (B2_ll + B2_rr) + f8 = 0.5f0 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) + f9 = equations.c_h * 0.5f0 * (B2_ll + B2_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v2_rr + p_rr * v2_ll + 0.5f0 * (+p_ll * v2_rr + p_rr * v2_ll + (v2_ll * B1_ll * B1_rr + v2_rr * B1_rr * B1_ll) + (v2_ll * B3_ll * B3_rr + v2_rr * B3_rr * B3_ll) - @@ -679,41 +684,41 @@ end # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = (f1 * v1_avg + (p_avg + magnetic_square_avg) * normal_direction[1] - - 0.5 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) + 0.5f0 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) f3 = (f1 * v2_avg + (p_avg + magnetic_square_avg) * normal_direction[2] - - 0.5 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) + 0.5f0 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) f4 = (f1 * v3_avg - - 0.5 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) + 0.5f0 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) #f5 below f6 = (equations.c_h * psi_avg * normal_direction[1] + - 0.5 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + v_dot_n_rr * B1_rr - v1_rr * B_dot_n_rr)) f7 = (equations.c_h * psi_avg * normal_direction[2] + - 0.5 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + v_dot_n_rr * B2_rr - v2_rr * B_dot_n_rr)) - f8 = +0.5 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + + f8 = +0.5f0 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + v_dot_n_rr * B3_rr - v3_rr * B_dot_n_rr) - f9 = equations.c_h * 0.5 * (B_dot_n_ll + B_dot_n_rr) + f9 = equations.c_h * 0.5f0 * (B_dot_n_ll + B_dot_n_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + 0.5f0 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + (v_dot_n_ll * B1_ll * B1_rr + v_dot_n_rr * B1_rr * B1_ll) + (v_dot_n_ll * B2_ll * B2_rr + v_dot_n_rr * B2_rr * B2_ll) + (v_dot_n_ll * B3_ll * B3_rr + v_dot_n_rr * B3_rr * B3_ll) @@ -1029,7 +1034,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * (rho_e - - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + B1 * B1 + B2 * B2 + B3 * B3 + psi * psi)) @@ -1045,11 +1050,12 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = rho_p * v3 @@ -1097,8 +1103,8 @@ end rho_v2 = rho * v2 rho_v3 = rho * v3 rho_e = p * equations.inv_gamma_minus_one + - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5 * psi^2 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) + 0.5f0 * psi^2 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -1110,11 +1116,11 @@ end @inline function pressure(u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return p end @@ -1128,17 +1134,17 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 - return (equations.gamma - 1.0) * - SVector(0.5 * v_square, -v1, -v2, -v3, 1.0, -B1, -B2, -B3, -psi) + return (equations.gamma - 1) * + SVector(0.5f0 * v_square, -v1, -v2, -v3, 1, -B1, -B2, -B3, -psi) end @inline function density_pressure(u, equations::IdealGlmMhdEquations2D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return rho * p end @@ -1149,9 +1155,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1159,11 +1165,11 @@ end b3 = B3 / sqrt_rho b_square = b1 * b1 + b2 * b2 + b3 * b3 if orientation == 1 # x-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) else - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b2^2)) end return c_f end @@ -1174,9 +1180,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1188,8 +1194,8 @@ end b_dot_n_squared = (b1 * normal_direction[1] + b2 * normal_direction[2])^2 / norm_squared - c_f = sqrt((0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * + c_f = sqrt((0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * norm_squared) return c_f end @@ -1213,28 +1219,28 @@ as given by v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1250,26 +1256,26 @@ as given by H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum - a_square_roe = ((2.0 - equations.gamma) * X + - (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) if orientation == 1 # x-direction c_a_roe = B1_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + 4 * a_square_roe * c_a_roe) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v1_roe else # y-direction c_a_roe = B2_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + 4 * a_square_roe * c_a_roe) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v2_roe end @@ -1285,28 +1291,28 @@ end v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1322,13 +1328,13 @@ end H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum - a_square_roe = ((2.0 - equations.gamma) * X + - (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) @@ -1339,7 +1345,7 @@ end c_a_roe = B_roe_dot_n_squared * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) vel_out_roe = (v1_roe * normal_direction[1] + v2_roe * normal_direction[2]) @@ -1350,11 +1356,11 @@ end @inline function entropy_thermodynamic(cons, equations::IdealGlmMhdEquations2D) # Pressure p = (equations.gamma - 1) * - (cons[5] - 1 / 2 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + (cons[5] - 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] - - 1 / 2 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) - - 1 / 2 * cons[9]^2) + 0.5f0 * cons[9]^2) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) @@ -1378,13 +1384,13 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(cons, equations::IdealGlmMhdEquations2D) - return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + return 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] end # Calculate the magnetic energy for a conservative state `cons'. # OBS! For non-dinmensional form of the ideal MHD magnetic pressure ≡ magnetic energy @inline function energy_magnetic(cons, ::IdealGlmMhdEquations2D) - return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + return 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) end # Calculate internal energy for a conservative state `cons` @@ -1401,7 +1407,7 @@ end # State validation for Newton-bisection method of subcell IDP limiting @inline function Base.isvalid(u, equations::IdealGlmMhdEquations2D) p = pressure(u, equations) - if u[1] <= 0.0 || p <= 0.0 + if u[1] <= 0 || p <= 0 return false end return true diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 321e501b051..2ffaa575243 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -46,15 +46,16 @@ end A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D) - rho = 1.0 - rho_v1 = 0.1 - rho_v2 = -0.2 - rho_v3 = -0.5 - rho_e = 50.0 - B1 = 3.0 - B2 = -1.2 - B3 = 0.5 - psi = 0.0 + RealT = eltype(x) + rho = 1 + rho_v1 = convert(RealT, 0.1) + rho_v2 = -convert(RealT, 0.2) + rho_v3 = -0.5f0 + rho_e = 50 + B1 = 3 + B2 = -convert(RealT, 1.2) + B3 = 0.5f0 + psi = 0 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -67,19 +68,20 @@ function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquation # Alfvén wave in three space dimensions # Altmann thesis http://dx.doi.org/10.18419/opus-3895 # domain must be set to [-1, 1]^3, γ = 5/3 + RealT = eltype(x) p = 1 - omega = 2 * pi # may be multiplied by frequency + omega = 2 * convert(RealT, pi) # may be multiplied by frequency # r: length-variable = length of computational domain - r = 2 + r = convert(RealT, 2) # e: epsilon = 0.2 - e = 0.2 + e = convert(RealT, 0.2) nx = 1 / sqrt(r^2 + 1) ny = r / sqrt(r^2 + 1) sqr = 1 Va = omega / (ny * sqr) - phi_alv = omega / ny * (nx * (x[1] - 0.5 * r) + ny * (x[2] - 0.5 * r)) - Va * t + phi_alv = omega / ny * (nx * (x[1] - 0.5f0 * r) + ny * (x[2] - 0.5f0 * r)) - Va * t - rho = 1.0 + rho = 1 v1 = -e * ny * cos(phi_alv) / rho v2 = e * nx * cos(phi_alv) / rho v3 = e * sin(phi_alv) / rho @@ -103,22 +105,23 @@ function initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdEquations # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3) # Same discontinuity in the velocities but with magnetic fields # Set up polar coordinates + RealT = eltype(x) inicenter = (0, 0, 0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] z_norm = x[3] - inicenter[3] r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) phi = atan(y_norm, x_norm) - theta = iszero(r) ? 0.0 : acos(z_norm / r) + theta = iszero(r) ? zero(RealT) : acos(z_norm / r) # Calculate primitive variables - rho = r > 0.5 ? 1.0 : 1.1691 - v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi) * sin(theta) - v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi) * sin(theta) - v3 = r > 0.5 ? 0.0 : 0.1882 * cos(theta) - p = r > 0.5 ? 1.0 : 1.245 + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) * sin(theta) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) * sin(theta) + v3 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(theta) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) - return prim2cons(SVector(rho, v1, v2, v3, p, 1.0, 1.0, 1.0, 0.0), equations) + return prim2cons(SVector(rho, v1, v2, v3, p, 1, 1, 1, 0), equations) end # Pre-defined source terms should be implemented as @@ -130,9 +133,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one if orientation == 1 f1 = rho_v1 @@ -180,9 +183,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p_over_gamma_minus_one = (rho_e - kin_en - mag_en - 0.5f0 * psi^2) p = (equations.gamma - 1) * p_over_gamma_minus_one v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + @@ -342,33 +345,33 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 - beta_ll = 0.5 * rho_ll / p_ll - beta_rr = 0.5 * rho_rr / p_rr + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr # for convenience store v⋅B vel_dot_mag_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr # Compute the necessary mean values needed for either direction - rho_avg = 0.5 * (rho_ll + rho_rr) + rho_avg = 0.5f0 * (rho_ll + rho_rr) rho_mean = ln_mean(rho_ll, rho_rr) beta_mean = ln_mean(beta_ll, beta_rr) - beta_avg = 0.5 * (beta_ll + beta_rr) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_mean = 0.5 * rho_avg / beta_avg - B1_avg = 0.5 * (B1_ll + B1_rr) - B2_avg = 0.5 * (B2_ll + B2_rr) - B3_avg = 0.5 * (B3_ll + B3_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr) - mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr) - vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg - f2 = f1 * v1_avg + p_mean + 0.5 * mag_norm_avg - B1_avg * B1_avg + f2 = f1 * v1_avg + p_mean + 0.5f0 * mag_norm_avg - B1_avg * B1_avg f3 = f1 * v2_avg - B1_avg * B2_avg f4 = f1 * v3_avg - B1_avg * B3_avg f6 = equations.c_h * psi_avg @@ -376,46 +379,46 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer, f8 = v1_avg * B3_avg - v3_avg * B1_avg f9 = equations.c_h * B1_avg # total energy flux is complicated and involves the previous eight components - psi_B1_avg = 0.5 * (B1_ll * psi_ll + B1_rr * psi_rr) - v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B1_avg = 0.5f0 * (B1_ll * psi_ll + B1_rr * psi_rr) + v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v1_mag_avg + + 0.5f0 * v1_mag_avg + B1_avg * vel_dot_mag_avg - equations.c_h * psi_B1_avg) elseif orientation == 2 f1 = rho_mean * v2_avg f2 = f1 * v1_avg - B2_avg * B1_avg - f3 = f1 * v2_avg + p_mean + 0.5 * mag_norm_avg - B2_avg * B2_avg + f3 = f1 * v2_avg + p_mean + 0.5f0 * mag_norm_avg - B2_avg * B2_avg f4 = f1 * v3_avg - B2_avg * B3_avg f6 = v2_avg * B1_avg - v1_avg * B2_avg f7 = equations.c_h * psi_avg f8 = v2_avg * B3_avg - v3_avg * B2_avg f9 = equations.c_h * B2_avg # total energy flux is complicated and involves the previous eight components - psi_B2_avg = 0.5 * (B2_ll * psi_ll + B2_rr * psi_rr) - v2_mag_avg = 0.5 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr) + v2_mag_avg = 0.5f0 * (v2_ll * mag_norm_ll + v2_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v2_mag_avg + + 0.5f0 * v2_mag_avg + B2_avg * vel_dot_mag_avg - equations.c_h * psi_B2_avg) else f1 = rho_mean * v3_avg f2 = f1 * v1_avg - B3_avg * B1_avg f3 = f1 * v2_avg - B3_avg * B2_avg - f4 = f1 * v3_avg + p_mean + 0.5 * mag_norm_avg - B3_avg * B3_avg + f4 = f1 * v3_avg + p_mean + 0.5f0 * mag_norm_avg - B3_avg * B3_avg f6 = v3_avg * B1_avg - v1_avg * B3_avg f7 = v3_avg * B2_avg - v2_avg * B3_avg f8 = equations.c_h * psi_avg f9 = equations.c_h * B3_avg # total energy flux is complicated and involves the previous eight components - psi_B3_avg = 0.5 * (B3_ll * psi_ll + B3_rr * psi_rr) - v3_mag_avg = 0.5 * (v3_ll * mag_norm_ll + v3_rr * mag_norm_rr) - f5 = (f1 * 0.5 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + + psi_B3_avg = 0.5f0 * (B3_ll * psi_ll + B3_rr * psi_rr) + v3_mag_avg = 0.5f0 * (v3_ll * mag_norm_ll + v3_rr * mag_norm_rr) + f5 = (f1 * 0.5f0 * (1 / (equations.gamma - 1) / beta_mean - vel_norm_avg) + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + f6 * B1_avg + f7 * B2_avg + f8 * B3_avg + f9 * psi_avg - - 0.5 * v3_mag_avg + + 0.5f0 * v3_mag_avg + B3_avg * vel_dot_mag_avg - equations.c_h * psi_B3_avg) end @@ -459,31 +462,31 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on orientation with specific direction averages if orientation == 1 f1 = rho_mean * v1_avg f2 = f1 * v1_avg + p_avg + magnetic_square_avg - - 0.5 * (B1_ll * B1_rr + B1_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll) + 0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll) #f5 below f6 = equations.c_h * psi_avg - f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) - f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) - f9 = equations.c_h * 0.5 * (B1_ll + B1_rr) + f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr) + f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr) + f9 = equations.c_h * 0.5f0 * (B1_ll + B1_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v1_rr + p_rr * v1_ll + 0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll + (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll) + (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll) - @@ -494,20 +497,20 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat equations.c_h * (B1_ll * psi_rr + B1_rr * psi_ll))) elseif orientation == 2 f1 = rho_mean * v2_avg - f2 = f1 * v1_avg - 0.5 * (B2_ll * B1_rr + B2_rr * B1_ll) + f2 = f1 * v1_avg - 0.5f0 * (B2_ll * B1_rr + B2_rr * B1_ll) f3 = f1 * v2_avg + p_avg + magnetic_square_avg - - 0.5 * (B2_ll * B2_rr + B2_rr * B2_ll) - f4 = f1 * v3_avg - 0.5 * (B2_ll * B3_rr + B2_rr * B3_ll) + 0.5f0 * (B2_ll * B2_rr + B2_rr * B2_ll) + f4 = f1 * v3_avg - 0.5f0 * (B2_ll * B3_rr + B2_rr * B3_ll) #f5 below - f6 = 0.5 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) + f6 = 0.5f0 * (v2_ll * B1_ll - v1_ll * B2_ll + v2_rr * B1_rr - v1_rr * B2_rr) f7 = equations.c_h * psi_avg - f8 = 0.5 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) - f9 = equations.c_h * 0.5 * (B2_ll + B2_rr) + f8 = 0.5f0 * (v2_ll * B3_ll - v3_ll * B2_ll + v2_rr * B3_rr - v3_rr * B2_rr) + f9 = equations.c_h * 0.5f0 * (B2_ll + B2_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v2_rr + p_rr * v2_ll + 0.5f0 * (+p_ll * v2_rr + p_rr * v2_ll + (v2_ll * B1_ll * B1_rr + v2_rr * B1_rr * B1_ll) + (v2_ll * B3_ll * B3_rr + v2_rr * B3_rr * B3_ll) - @@ -518,20 +521,20 @@ Hindenlang and Gassner (2019), extending [`flux_ranocha`](@ref) to the MHD equat equations.c_h * (B2_ll * psi_rr + B2_rr * psi_ll))) else # orientation == 3 f1 = rho_mean * v3_avg - f2 = f1 * v1_avg - 0.5 * (B3_ll * B1_rr + B3_rr * B1_ll) - f3 = f1 * v2_avg - 0.5 * (B3_ll * B2_rr + B3_rr * B2_ll) + f2 = f1 * v1_avg - 0.5f0 * (B3_ll * B1_rr + B3_rr * B1_ll) + f3 = f1 * v2_avg - 0.5f0 * (B3_ll * B2_rr + B3_rr * B2_ll) f4 = f1 * v3_avg + p_avg + magnetic_square_avg - - 0.5 * (B3_ll * B3_rr + B3_rr * B3_ll) + 0.5f0 * (B3_ll * B3_rr + B3_rr * B3_ll) #f5 below - f6 = 0.5 * (v3_ll * B1_ll - v1_ll * B3_ll + v3_rr * B1_rr - v1_rr * B3_rr) - f7 = 0.5 * (v3_ll * B2_ll - v2_ll * B3_ll + v3_rr * B2_rr - v2_rr * B3_rr) + f6 = 0.5f0 * (v3_ll * B1_ll - v1_ll * B3_ll + v3_rr * B1_rr - v1_rr * B3_rr) + f7 = 0.5f0 * (v3_ll * B2_ll - v2_ll * B3_ll + v3_rr * B2_rr - v2_rr * B3_rr) f8 = equations.c_h * psi_avg - f9 = equations.c_h * 0.5 * (B3_ll + B3_rr) + f9 = equations.c_h * 0.5f0 * (B3_ll + B3_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v3_rr + p_rr * v3_ll + 0.5f0 * (+p_ll * v3_rr + p_rr * v3_ll + (v3_ll * B1_ll * B1_rr + v3_rr * B1_rr * B1_ll) + (v3_ll * B2_ll * B2_rr + v3_rr * B2_rr * B2_ll) - @@ -568,43 +571,43 @@ end # log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ) # = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ) inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll) - v1_avg = 0.5 * (v1_ll + v1_rr) - v2_avg = 0.5 * (v2_ll + v2_rr) - v3_avg = 0.5 * (v3_ll + v3_rr) - p_avg = 0.5 * (p_ll + p_rr) - psi_avg = 0.5 * (psi_ll + psi_rr) - velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) - magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + p_avg = 0.5f0 * (p_ll + p_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr) + magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr) # Calculate fluxes depending on normal_direction - f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr) + f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr) f2 = (f1 * v1_avg + (p_avg + magnetic_square_avg) * normal_direction[1] - - 0.5 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) + 0.5f0 * (B_dot_n_ll * B1_rr + B_dot_n_rr * B1_ll)) f3 = (f1 * v2_avg + (p_avg + magnetic_square_avg) * normal_direction[2] - - 0.5 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) + 0.5f0 * (B_dot_n_ll * B2_rr + B_dot_n_rr * B2_ll)) f4 = (f1 * v3_avg + (p_avg + magnetic_square_avg) * normal_direction[3] - - 0.5 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) + 0.5f0 * (B_dot_n_ll * B3_rr + B_dot_n_rr * B3_ll)) #f5 below f6 = (equations.c_h * psi_avg * normal_direction[1] + - 0.5 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B1_ll - v1_ll * B_dot_n_ll + v_dot_n_rr * B1_rr - v1_rr * B_dot_n_rr)) f7 = (equations.c_h * psi_avg * normal_direction[2] + - 0.5 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B2_ll - v2_ll * B_dot_n_ll + v_dot_n_rr * B2_rr - v2_rr * B_dot_n_rr)) f8 = (equations.c_h * psi_avg * normal_direction[3] + - 0.5 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + + 0.5f0 * (v_dot_n_ll * B3_ll - v3_ll * B_dot_n_ll + v_dot_n_rr * B3_rr - v3_rr * B_dot_n_rr)) - f9 = equations.c_h * 0.5 * (B_dot_n_ll + B_dot_n_rr) + f9 = equations.c_h * 0.5f0 * (B_dot_n_ll + B_dot_n_rr) # total energy flux is complicated and involves the previous components f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) + - 0.5 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + 0.5f0 * (+p_ll * v_dot_n_rr + p_rr * v_dot_n_ll + (v_dot_n_ll * B1_ll * B1_rr + v_dot_n_rr * B1_rr * B1_ll) + (v_dot_n_ll * B2_ll * B2_rr + v_dot_n_rr * B2_rr * B2_ll) + (v_dot_n_ll * B3_ll * B3_rr + v_dot_n_rr * B3_rr * B3_ll) @@ -951,7 +954,7 @@ end v2 = rho_v2 / rho v3 = rho_v3 / rho p = (equations.gamma - 1) * (rho_e - - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + B1 * B1 + B2 * B2 + B3 * B3 + psi * psi)) @@ -966,11 +969,12 @@ end v3 = rho_v3 / rho v_square = v1^2 + v2^2 + v3^2 p = (equations.gamma - 1) * - (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2) - 0.5 * psi^2) + (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2) - 0.5f0 * psi^2) s = log(p) - equations.gamma * log(rho) rho_p = rho / p - w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - 0.5 * rho_p * v_square + w1 = (equations.gamma - s) * equations.inv_gamma_minus_one - + 0.5f0 * rho_p * v_square w2 = rho_p * v1 w3 = rho_p * v2 w4 = rho_p * v3 @@ -991,8 +995,8 @@ end rho_v2 = rho * v2 rho_v3 = rho * v3 rho_e = p * equations.inv_gamma_minus_one + - 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5 * psi^2 + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) + 0.5f0 * psi^2 return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi) end @@ -1003,21 +1007,21 @@ end @inline function pressure(u, equations::IdealGlmMhdEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return p end @inline function density_pressure(u, equations::IdealGlmMhdEquations3D) rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u - p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho - - 0.5 * (B1^2 + B2^2 + B3^2) + 0.5f0 * (B1^2 + B2^2 + B3^2) - - 0.5 * psi^2) + 0.5f0 * psi^2) return rho * p end @@ -1028,9 +1032,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1038,14 +1042,14 @@ end b3 = B3 / sqrt_rho b_square = b1 * b1 + b2 * b2 + b3 * b3 if orientation == 1 # x-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)) elseif orientation == 2 # y-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b2^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b2^2)) else # z-direction - c_f = sqrt(0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b3^2)) + c_f = sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b3^2)) end return c_f end @@ -1056,9 +1060,9 @@ end v1 = rho_v1 / rho v2 = rho_v2 / rho v3 = rho_v3 / rho - kin_en = 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) - mag_en = 0.5 * (B1 * B1 + B2 * B2 + B3 * B3) - p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5 * psi^2) + kin_en = 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + mag_en = 0.5f0 * (B1 * B1 + B2 * B2 + B3 * B3) + p = (equations.gamma - 1) * (rho_e - kin_en - mag_en - 0.5f0 * psi^2) a_square = equations.gamma * p / rho sqrt_rho = sqrt(rho) b1 = B1 / sqrt_rho @@ -1072,8 +1076,8 @@ end b2 * normal_direction[2] + b3 * normal_direction[3])^2 / norm_squared - c_f = sqrt((0.5 * (a_square + b_square) + - 0.5 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * + c_f = sqrt((0.5f0 * (a_square + b_square) + + 0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b_dot_n_squared)) * norm_squared) return c_f end @@ -1096,28 +1100,28 @@ Compute the fast magnetoacoustic wave speed using Roe averages as given by v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1133,32 +1137,32 @@ Compute the fast magnetoacoustic wave speed using Roe averages as given by H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum - a_square_roe = ((2.0 - equations.gamma) * X + - (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) if orientation == 1 # x-direction c_a_roe = B1_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + 4 * a_square_roe * c_a_roe) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v1_roe elseif orientation == 2 # y-direction c_a_roe = B2_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + 4 * a_square_roe * c_a_roe) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v2_roe else # z-direction c_a_roe = B3_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - - 4.0 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe)) + 4 * a_square_roe * c_a_roe) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe)) vel_out_roe = v3_roe end @@ -1174,28 +1178,28 @@ end v1_ll = rho_v1_ll / rho_ll v2_ll = rho_v2_ll / rho_ll v3_ll = rho_v3_ll / rho_ll - kin_en_ll = 0.5 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) + kin_en_ll = 0.5f0 * (rho_v1_ll * v1_ll + rho_v2_ll * v2_ll + rho_v3_ll * v3_ll) mag_norm_ll = B1_ll * B1_ll + B2_ll * B2_ll + B3_ll * B3_ll p_ll = (equations.gamma - 1) * - (rho_e_ll - kin_en_ll - 0.5 * mag_norm_ll - 0.5 * psi_ll^2) + (rho_e_ll - kin_en_ll - 0.5f0 * mag_norm_ll - 0.5f0 * psi_ll^2) v1_rr = rho_v1_rr / rho_rr v2_rr = rho_v2_rr / rho_rr v3_rr = rho_v3_rr / rho_rr - kin_en_rr = 0.5 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) + kin_en_rr = 0.5f0 * (rho_v1_rr * v1_rr + rho_v2_rr * v2_rr + rho_v3_rr * v3_rr) mag_norm_rr = B1_rr * B1_rr + B2_rr * B2_rr + B3_rr * B3_rr p_rr = (equations.gamma - 1) * - (rho_e_rr - kin_en_rr - 0.5 * mag_norm_rr - 0.5 * psi_rr^2) + (rho_e_rr - kin_en_rr - 0.5f0 * mag_norm_rr - 0.5f0 * psi_rr^2) # compute total pressure which is thermal + magnetic pressures - p_total_ll = p_ll + 0.5 * mag_norm_ll - p_total_rr = p_rr + 0.5 * mag_norm_rr + p_total_ll = p_ll + 0.5f0 * mag_norm_ll + p_total_rr = p_rr + 0.5f0 * mag_norm_rr # compute the Roe density averages sqrt_rho_ll = sqrt(rho_ll) sqrt_rho_rr = sqrt(rho_rr) - inv_sqrt_rho_add = 1.0 / (sqrt_rho_ll + sqrt_rho_rr) - inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr) + inv_sqrt_rho_add = 1 / (sqrt_rho_ll + sqrt_rho_rr) + inv_sqrt_rho_prod = 1 / (sqrt_rho_ll * sqrt_rho_rr) rho_ll_roe = sqrt_rho_ll * inv_sqrt_rho_add rho_rr_roe = sqrt_rho_rr * inv_sqrt_rho_add # Roe averages @@ -1211,13 +1215,13 @@ end H_rr = (rho_e_rr + p_total_rr) / rho_rr H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe # temporary variable see equation (4.12) in Cargo and Gallice - X = 0.5 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * + X = 0.5f0 * ((B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2) * inv_sqrt_rho_add^2 # averaged components needed to compute c_f, the fast magnetoacoustic wave speed b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum - a_square_roe = ((2.0 - equations.gamma) * X + - (equations.gamma - 1.0) * - (H_roe - 0.5 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - + a_square_roe = ((2 - equations.gamma) * X + + (equations.gamma - 1) * + (H_roe - 0.5f0 * (v1_roe^2 + v2_roe^2 + v3_roe^2) - b_square_roe)) # acoustic speed # finally compute the average wave speed and set the output velocity (depends on orientation) @@ -1230,7 +1234,7 @@ end c_a_roe = B_roe_dot_n_squared * inv_sqrt_rho_prod # (squared) Alfvén wave speed a_star_roe = sqrt((a_square_roe + b_square_roe)^2 - 4 * a_square_roe * c_a_roe) - c_f_roe = sqrt(0.5 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) + c_f_roe = sqrt(0.5f0 * (a_square_roe + b_square_roe + a_star_roe) * norm_squared) vel_out_roe = (v1_roe * normal_direction[1] + v2_roe * normal_direction[2] + v3_roe * normal_direction[3]) @@ -1242,11 +1246,11 @@ end @inline function entropy_thermodynamic(cons, equations::IdealGlmMhdEquations3D) # Pressure p = (equations.gamma - 1) * - (cons[5] - 1 / 2 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + (cons[5] - 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] - - 1 / 2 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) - - 1 / 2 * cons[9]^2) + 0.5f0 * cons[9]^2) # Thermodynamic entropy s = log(p) - equations.gamma * log(cons[1]) @@ -1270,13 +1274,13 @@ end # Calculate kinetic energy for a conservative state `cons` @inline function energy_kinetic(cons, equations::IdealGlmMhdEquations3D) - return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] + return 0.5f0 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1] end # Calculate the magnetic energy for a conservative state `cons'. # OBS! For non-dinmensional form of the ideal MHD magnetic pressure ≡ magnetic energy @inline function energy_magnetic(cons, ::IdealGlmMhdEquations3D) - return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2) + return 0.5f0 * (cons[6]^2 + cons[7]^2 + cons[8]^2) end # Calculate internal energy for a conservative state `cons` diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl index 6c6b9dd3721..743d2df870a 100644 --- a/src/equations/linear_scalar_advection_1d.jl +++ b/src/equations/linear_scalar_advection_1d.jl @@ -34,9 +34,10 @@ A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation1D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - return SVector(2.0) + return SVector(RealT(2)) end """ @@ -49,13 +50,14 @@ in non-periodic domains). function initial_condition_convergence_test(x, t, equation::LinearScalarAdvectionEquation1D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - c = 1.0 - A = 0.5 + c = 1 + A = 0.5f0 L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * sum(x_trans)) return SVector(scalar) end @@ -161,7 +163,7 @@ function flux_engquist_osher(u_ll, u_rr, orientation::Int, u_L = u_ll[1] u_R = u_rr[1] - return SVector(0.5 * (flux(u_L, orientation, equation) + + return SVector(0.5f0 * (flux(u_L, orientation, equation) + flux(u_R, orientation, equation) - abs(equation.advection_velocity[orientation]) * (u_R - u_L))) end @@ -200,14 +202,16 @@ end @inline function splitting_lax_friedrichs(u, ::Val{:plus}, orientation::Integer, equations::LinearScalarAdvectionEquation1D) + RealT = eltype(u) a = equations.advection_velocity[1] - return a > 0 ? flux(u, orientation, equations) : zero(u) + return a > 0 ? flux(u, orientation, equations) : SVector(zero(RealT)) end @inline function splitting_lax_friedrichs(u, ::Val{:minus}, orientation::Integer, equations::LinearScalarAdvectionEquation1D) + RealT = eltype(u) a = equations.advection_velocity[1] - return a < 0 ? flux(u, orientation, equations) : zero(u) + return a < 0 ? flux(u, orientation, equations) : SVector(zero(RealT)) end # Convert conservative variables to primitive @@ -217,11 +221,11 @@ end @inline cons2entropy(u, equation::LinearScalarAdvectionEquation1D) = u # Calculate entropy for a conservative state `cons` -@inline entropy(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5 * u^2 +@inline entropy(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5f0 * u^2 @inline entropy(u, equation::LinearScalarAdvectionEquation1D) = entropy(u[1], equation) # Calculate total energy for a conservative state `cons` -@inline energy_total(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5 * u^2 +@inline energy_total(u::Real, ::LinearScalarAdvectionEquation1D) = 0.5f0 * u^2 @inline function energy_total(u, equation::LinearScalarAdvectionEquation1D) energy_total(u[1], equation) end diff --git a/src/equations/linear_scalar_advection_2d.jl b/src/equations/linear_scalar_advection_2d.jl index d90bf0c8793..5e4f8463f52 100644 --- a/src/equations/linear_scalar_advection_2d.jl +++ b/src/equations/linear_scalar_advection_2d.jl @@ -34,8 +34,8 @@ varnames(::typeof(cons2prim), ::LinearScalarAdvectionEquation2D) = ("scalar",) function x_trans_periodic_2d(x, domain_length = SVector(10, 10), center = SVector(0, 0)) x_normalized = x .- center x_shifted = x_normalized .% domain_length - x_offset = ((x_shifted .< -0.5 * domain_length) - - (x_shifted .> 0.5 * domain_length)) .* domain_length + x_offset = ((x_shifted .< -0.5f0 * domain_length) - + (x_shifted .> 0.5f0 * domain_length)) .* domain_length return center + x_shifted + x_offset end @@ -47,9 +47,10 @@ A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation2D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x_trans_periodic_2d(x - equation.advection_velocity * t) - return SVector(2.0) + return SVector(RealT(2)) end """ @@ -60,13 +61,14 @@ A smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equation::LinearScalarAdvectionEquation2D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - c = 1.0 - A = 0.5 + c = 1 + A = 0.5f0 L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * sum(x_trans)) return SVector(scalar) end @@ -280,11 +282,11 @@ end @inline cons2entropy(u, equation::LinearScalarAdvectionEquation2D) = u # Calculate entropy for a conservative state `cons` -@inline entropy(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5 * u^2 +@inline entropy(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5f0 * u^2 @inline entropy(u, equation::LinearScalarAdvectionEquation2D) = entropy(u[1], equation) # Calculate total energy for a conservative state `cons` -@inline energy_total(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5 * u^2 +@inline energy_total(u::Real, ::LinearScalarAdvectionEquation2D) = 0.5f0 * u^2 @inline function energy_total(u, equation::LinearScalarAdvectionEquation2D) energy_total(u[1], equation) end diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index 7b19974eb49..088f934cc3e 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/src/equations/linear_scalar_advection_3d.jl @@ -38,9 +38,10 @@ A constant initial condition to test free-stream preservation. """ function initial_condition_constant(x, t, equation::LinearScalarAdvectionEquation3D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - return SVector(2.0) + return SVector(RealT(2)) end """ @@ -51,13 +52,14 @@ A smooth initial condition used for convergence tests. function initial_condition_convergence_test(x, t, equation::LinearScalarAdvectionEquation3D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - c = 1.0 - A = 0.5 + c = 1 + A = 0.5f0 L = 2 - f = 1 / L - omega = 2 * pi * f + f = 1.0f0 / L + omega = 2 * convert(RealT, pi) * f scalar = c + A * sin(omega * sum(x_trans)) return SVector(scalar) end @@ -82,10 +84,10 @@ A sine wave in the conserved variable. """ function initial_condition_sin(x, t, equation::LinearScalarAdvectionEquation3D) # Store translated coordinate for easy use of exact solution + RealT = eltype(x) x_trans = x - equation.advection_velocity * t - scalar = sin(2 * pi * x_trans[1]) * sin(2 * pi * x_trans[2]) * - sin(2 * pi * x_trans[3]) + scalar = sinpi(2 * x_trans[1]) * sinpi(2 * x_trans[2]) * sinpi(2 * x_trans[3]) return SVector(scalar) end @@ -199,11 +201,11 @@ end @inline cons2entropy(u, equation::LinearScalarAdvectionEquation3D) = u # Calculate entropy for a conservative state `cons` -@inline entropy(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5 * u^2 +@inline entropy(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5f0 * u^2 @inline entropy(u, equation::LinearScalarAdvectionEquation3D) = entropy(u[1], equation) # Calculate total energy for a conservative state `cons` -@inline energy_total(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5 * u^2 +@inline energy_total(u::Real, ::LinearScalarAdvectionEquation3D) = 0.5f0 * u^2 @inline function energy_total(u, equation::LinearScalarAdvectionEquation3D) energy_total(u[1], equation) end diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index b74a3b4d642..127bc420f65 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -18,7 +18,7 @@ function save_mesh_file(mesh::TreeMesh, output_directory, timestep, # Determine file name based on existence of meaningful time step if timestep > 0 - filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep)) + filename = joinpath(output_directory, @sprintf("mesh_%09d.h5", timestep)) else filename = joinpath(output_directory, "mesh.h5") end @@ -57,7 +57,7 @@ function save_mesh_file(mesh::TreeMesh, output_directory, timestep, # Determine file name based on existence of meaningful time step if timestep >= 0 - filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep)) + filename = joinpath(output_directory, @sprintf("mesh_%09d.h5", timestep)) else filename = joinpath(output_directory, "mesh.h5") end @@ -156,8 +156,8 @@ function save_mesh_file(mesh::P4estMesh, output_directory, timestep, # Determine file name based on existence of meaningful time step if timestep > 0 - filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep)) - p4est_filename = @sprintf("p4est_data_%06d", timestep) + filename = joinpath(output_directory, @sprintf("mesh_%09d.h5", timestep)) + p4est_filename = @sprintf("p4est_data_%09d", timestep) else filename = joinpath(output_directory, "mesh.h5") p4est_filename = "p4est_data" @@ -191,8 +191,8 @@ function save_mesh_file(mesh::P4estMesh, output_directory, timestep, mpi_paralle # Determine file name based on existence of meaningful time step if timestep > 0 - filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep)) - p4est_filename = @sprintf("p4est_data_%06d", timestep) + filename = joinpath(output_directory, @sprintf("mesh_%09d.h5", timestep)) + p4est_filename = @sprintf("p4est_data_%09d", timestep) else filename = joinpath(output_directory, "mesh.h5") p4est_filename = "p4est_data" diff --git a/src/meshes/structured_mesh_view.jl b/src/meshes/structured_mesh_view.jl index bd55115cc90..0b0cccfc7fc 100644 --- a/src/meshes/structured_mesh_view.jl +++ b/src/meshes/structured_mesh_view.jl @@ -116,7 +116,7 @@ function save_mesh_file(mesh::StructuredMeshView, output_directory; system = "", # Create output directory (if it does not exist) mkpath(output_directory) - filename = joinpath(output_directory, @sprintf("mesh_%s_%06d.h5", system, timestep)) + filename = joinpath(output_directory, @sprintf("mesh_%s_%09d.h5", system, timestep)) # Open file (clobber existing content) h5open(filename, "w") do file diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index c5c21584dca..6b009cfad20 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -125,7 +125,7 @@ end """ ndofsglobal(semi::SemidiscretizationCoupled) - + Return the global number of degrees of freedom associated with each scalar variable across all MPI ranks, and summed up over all coupled systems. This is the same as [`ndofs`](@ref) for simulations running in serial or parallelized via threads. It will in general be different for simulations @@ -180,12 +180,10 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationCoupled, t) end # Call rhs! for each semidiscretization - @trixi_timeit timer() "copy to coupled boundaries" begin - foreach_enumerate(semi.semis) do (i, semi_) - u_loc = get_system_u_ode(u_ode, i, semi) - du_loc = get_system_u_ode(du_ode, i, semi) - rhs!(du_loc, u_loc, semi_, t) - end + foreach_enumerate(semi.semis) do (i, semi_) + u_loc = get_system_u_ode(u_ode, i, semi) + du_loc = get_system_u_ode(du_ode, i, semi) + rhs!(du_loc, u_loc, semi_, t) end runtime = time_ns() - time_start @@ -431,12 +429,16 @@ BoundaryConditionCoupled(2, (:j, :i_backwards, :end), Float64, fun) !!! warning "Experimental code" This is an experimental feature and can change any time. """ -mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indices, - CouplingConverter} +mutable struct BoundaryConditionCoupled{NDIMS, + # Store the other semi index as type parameter, + # so that retrieving the other semidiscretization + # is type-stable. + # x-ref: https://github.com/trixi-framework/Trixi.jl/pull/1979 + other_semi_index, NDIMST2M1, + uEltype <: Real, Indices, CouplingConverter} # NDIMST2M1 == NDIMS * 2 - 1 # Buffer for boundary values: [variable, nodes_i, nodes_j, cell_i, cell_j] u_boundary :: Array{uEltype, NDIMST2M1} # NDIMS * 2 - 1 - other_semi_index :: Int other_orientation :: Int indices :: Indices coupling_converter :: CouplingConverter @@ -454,9 +456,9 @@ mutable struct BoundaryConditionCoupled{NDIMS, NDIMST2M1, uEltype <: Real, Indic other_orientation = 3 end - new{NDIMS, NDIMS * 2 - 1, uEltype, typeof(indices), + new{NDIMS, other_semi_index, NDIMS * 2 - 1, uEltype, typeof(indices), typeof(coupling_converter)}(u_boundary, - other_semi_index, other_orientation, + other_orientation, indices, coupling_converter) end end @@ -560,28 +562,17 @@ function copy_to_coupled_boundary!(boundary_conditions::Union{Tuple, NamedTuple} boundary_conditions...) end -function mesh_equations_solver_cache(other_semi_index, i, semi_, semi_tuple...) - if i == other_semi_index - return mesh_equations_solver_cache(semi_) - else - # Walk through semidiscretizations until we find `i` - mesh_equations_solver_cache(other_semi_index, i + 1, semi_tuple...) - end -end - # In 2D -function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2}, - u_ode, - semi_coupled, semi) +function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{2, + other_semi_index}, + u_ode, semi_coupled, semi) where {other_semi_index} @unpack u_indices = semi_coupled - @unpack other_semi_index, other_orientation, indices = boundary_condition + @unpack other_orientation, indices = boundary_condition @unpack coupling_converter, u_boundary = boundary_condition mesh_own, equations_own, solver_own, cache_own = mesh_equations_solver_cache(semi) - - mesh_other, equations_other, solver_other, cache_other = mesh_equations_solver_cache(other_semi_index, - 1, - semi_coupled.semis...) + other_semi = semi_coupled.semis[other_semi_index] + mesh_other, equations_other, solver_other, cache_other = mesh_equations_solver_cache(other_semi) node_coordinates_other = cache_other.elements.node_coordinates u_ode_other = get_system_u_ode(u_ode, other_semi_index, semi_coupled) @@ -605,10 +596,14 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ i_cell_start, i_cell_step = index_to_start_step_2d(indices[1], axes(mesh_other, 1)) j_cell_start, j_cell_step = index_to_start_step_2d(indices[2], axes(mesh_other, 2)) - i_cell = i_cell_start - j_cell = j_cell_start + # We need indices starting at 1 for the handling of `i_cell` etc. + Base.require_one_based_indexing(cells) + + @threaded for i in eachindex(cells) + cell = cells[i] + i_cell = i_cell_start + (i - 1) * i_cell_step + j_cell = j_cell_start + (i - 1) * j_cell_step - for cell in cells i_node = i_node_start j_node = j_node_start element_id = linear_indices[i_cell, j_cell] @@ -630,9 +625,6 @@ function copy_to_coupled_boundary!(boundary_condition::BoundaryConditionCoupled{ i_node += i_node_step j_node += j_node_step end - - i_cell += i_cell_step - j_cell += j_cell_step end end diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl index f3b09b01e97..e5b970c6bda 100644 --- a/src/time_integration/methods_2N.jl +++ b/src/time_integration/methods_2N.jl @@ -104,9 +104,8 @@ function Base.getproperty(integrator::SimpleIntegrator2N, field::Symbol) return getfield(integrator, field) end -# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 -function solve(ode::ODEProblem, alg::T; - dt, callback = nothing, kwargs...) where {T <: SimpleAlgorithm2N} +function init(ode::ODEProblem, alg::SimpleAlgorithm2N; + dt, callback = nothing, kwargs...) u = copy(ode.u0) du = similar(u) u_tmp = similar(u) @@ -129,67 +128,84 @@ function solve(ode::ODEProblem, alg::T; error("unsupported") end + return integrator +end + +# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 +function solve(ode::ODEProblem, alg::SimpleAlgorithm2N; + dt, callback = nothing, kwargs...) + integrator = init(ode, alg, dt = dt, callback = callback; kwargs...) + + # Start actual solve solve!(integrator) end function solve!(integrator::SimpleIntegrator2N) @unpack prob = integrator.sol + + integrator.finalstep = false + + @trixi_timeit timer() "main loop" while !integrator.finalstep + step!(integrator) + end # "main loop" timer + + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), + integrator.sol.prob) +end + +function step!(integrator::SimpleIntegrator2N) + @unpack prob = integrator.sol @unpack alg = integrator t_end = last(prob.tspan) callbacks = integrator.opts.callback - integrator.finalstep = false - @trixi_timeit timer() "main loop" while !integrator.finalstep - if isnan(integrator.dt) - error("time step size `dt` is NaN") - end + @assert !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") + end - # if the next iteration would push the simulation beyond the end time, set dt accordingly - if integrator.t + integrator.dt > t_end || - isapprox(integrator.t + integrator.dt, t_end) - integrator.dt = t_end - integrator.t - terminate!(integrator) - end + # if the next iteration would push the simulation beyond the end time, set dt accordingly + if integrator.t + integrator.dt > t_end || + isapprox(integrator.t + integrator.dt, t_end) + integrator.dt = t_end - integrator.t + terminate!(integrator) + end - # one time step - integrator.u_tmp .= 0 - for stage in eachindex(alg.c) - t_stage = integrator.t + integrator.dt * alg.c[stage] - integrator.f(integrator.du, integrator.u, prob.p, t_stage) - - a_stage = alg.a[stage] - b_stage_dt = alg.b[stage] * integrator.dt - @trixi_timeit timer() "Runge-Kutta step" begin - @threaded for i in eachindex(integrator.u) - integrator.u_tmp[i] = integrator.du[i] - - integrator.u_tmp[i] * a_stage - integrator.u[i] += integrator.u_tmp[i] * b_stage_dt - end + # one time step + integrator.u_tmp .= 0 + for stage in eachindex(alg.c) + t_stage = integrator.t + integrator.dt * alg.c[stage] + integrator.f(integrator.du, integrator.u, prob.p, t_stage) + + a_stage = alg.a[stage] + b_stage_dt = alg.b[stage] * integrator.dt + @trixi_timeit timer() "Runge-Kutta step" begin + @threaded for i in eachindex(integrator.u) + integrator.u_tmp[i] = integrator.du[i] - + integrator.u_tmp[i] * a_stage + integrator.u[i] += integrator.u_tmp[i] * b_stage_dt end end - integrator.iter += 1 - integrator.t += integrator.dt - - # handle callbacks - if callbacks isa CallbackSet - foreach(callbacks.discrete_callbacks) do cb - if cb.condition(integrator.u, integrator.t, integrator) - cb.affect!(integrator) - end - return nothing + end + integrator.iter += 1 + integrator.t += integrator.dt + + # handle callbacks + if callbacks isa CallbackSet + foreach(callbacks.discrete_callbacks) do cb + if cb.condition(integrator.u, integrator.t, integrator) + cb.affect!(integrator) end - end - - # respect maximum number of iterations - if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep - @warn "Interrupted. Larger maxiters is needed." - terminate!(integrator) + return nothing end end - return TimeIntegratorSolution((first(prob.tspan), integrator.t), - (prob.u0, integrator.u), - integrator.sol.prob) + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) + end end # get a cache where the RHS can be stored diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl index 7b70466606c..6128d1551d8 100644 --- a/src/time_integration/methods_3Sstar.jl +++ b/src/time_integration/methods_3Sstar.jl @@ -171,9 +171,8 @@ function Base.getproperty(integrator::SimpleIntegrator3Sstar, field::Symbol) return getfield(integrator, field) end -# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 -function solve(ode::ODEProblem, alg::T; - dt, callback = nothing, kwargs...) where {T <: SimpleAlgorithm3Sstar} +function init(ode::ODEProblem, alg::SimpleAlgorithm3Sstar; + dt, callback = nothing, kwargs...) u = copy(ode.u0) du = similar(u) u_tmp1 = similar(u) @@ -199,73 +198,90 @@ function solve(ode::ODEProblem, alg::T; error("unsupported") end + return integrator +end + +# Fakes `solve`: https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1 +function solve(ode::ODEProblem, alg::SimpleAlgorithm3Sstar; + dt, callback = nothing, kwargs...) + integrator = init(ode, alg, dt = dt, callback = callback; kwargs...) + + # Start actual solve solve!(integrator) end function solve!(integrator::SimpleIntegrator3Sstar) @unpack prob = integrator.sol + + integrator.finalstep = false + + @trixi_timeit timer() "main loop" while !integrator.finalstep + step!(integrator) + end # "main loop" timer + + return TimeIntegratorSolution((first(prob.tspan), integrator.t), + (prob.u0, integrator.u), + integrator.sol.prob) +end + +function step!(integrator::SimpleIntegrator3Sstar) + @unpack prob = integrator.sol @unpack alg = integrator t_end = last(prob.tspan) callbacks = integrator.opts.callback - integrator.finalstep = false - @trixi_timeit timer() "main loop" while !integrator.finalstep - if isnan(integrator.dt) - error("time step size `dt` is NaN") - end + @assert !integrator.finalstep + if isnan(integrator.dt) + error("time step size `dt` is NaN") + end - # if the next iteration would push the simulation beyond the end time, set dt accordingly - if integrator.t + integrator.dt > t_end || - isapprox(integrator.t + integrator.dt, t_end) - integrator.dt = t_end - integrator.t - terminate!(integrator) - end + # if the next iteration would push the simulation beyond the end time, set dt accordingly + if integrator.t + integrator.dt > t_end || + isapprox(integrator.t + integrator.dt, t_end) + integrator.dt = t_end - integrator.t + terminate!(integrator) + end - # one time step - integrator.u_tmp1 .= zero(eltype(integrator.u_tmp1)) - integrator.u_tmp2 .= integrator.u - for stage in eachindex(alg.c) - t_stage = integrator.t + integrator.dt * alg.c[stage] - prob.f(integrator.du, integrator.u, prob.p, t_stage) - - delta_stage = alg.delta[stage] - gamma1_stage = alg.gamma1[stage] - gamma2_stage = alg.gamma2[stage] - gamma3_stage = alg.gamma3[stage] - beta_stage_dt = alg.beta[stage] * integrator.dt - @trixi_timeit timer() "Runge-Kutta step" begin - @threaded for i in eachindex(integrator.u) - integrator.u_tmp1[i] += delta_stage * integrator.u[i] - integrator.u[i] = (gamma1_stage * integrator.u[i] + - gamma2_stage * integrator.u_tmp1[i] + - gamma3_stage * integrator.u_tmp2[i] + - beta_stage_dt * integrator.du[i]) - end + # one time step + integrator.u_tmp1 .= zero(eltype(integrator.u_tmp1)) + integrator.u_tmp2 .= integrator.u + for stage in eachindex(alg.c) + t_stage = integrator.t + integrator.dt * alg.c[stage] + prob.f(integrator.du, integrator.u, prob.p, t_stage) + + delta_stage = alg.delta[stage] + gamma1_stage = alg.gamma1[stage] + gamma2_stage = alg.gamma2[stage] + gamma3_stage = alg.gamma3[stage] + beta_stage_dt = alg.beta[stage] * integrator.dt + @trixi_timeit timer() "Runge-Kutta step" begin + @threaded for i in eachindex(integrator.u) + integrator.u_tmp1[i] += delta_stage * integrator.u[i] + integrator.u[i] = (gamma1_stage * integrator.u[i] + + gamma2_stage * integrator.u_tmp1[i] + + gamma3_stage * integrator.u_tmp2[i] + + beta_stage_dt * integrator.du[i]) end end - integrator.iter += 1 - integrator.t += integrator.dt - - # handle callbacks - if callbacks isa CallbackSet - foreach(callbacks.discrete_callbacks) do cb - if cb.condition(integrator.u, integrator.t, integrator) - cb.affect!(integrator) - end - return nothing + end + integrator.iter += 1 + integrator.t += integrator.dt + + # handle callbacks + if callbacks isa CallbackSet + foreach(callbacks.discrete_callbacks) do cb + if cb.condition(integrator.u, integrator.t, integrator) + cb.affect!(integrator) end - end - - # respect maximum number of iterations - if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep - @warn "Interrupted. Larger maxiters is needed." - terminate!(integrator) + return nothing end end - return TimeIntegratorSolution((first(prob.tspan), integrator.t), - (prob.u0, integrator.u), - integrator.sol.prob) + # respect maximum number of iterations + if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep + @warn "Interrupted. Larger maxiters is needed." + terminate!(integrator) + end end # get a cache where the RHS can be stored diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index 1d9680153e6..8a5c2a24617 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -296,8 +296,6 @@ function step!(integrator::PairedExplicitRK2Integrator) t_end = last(prob.tspan) callbacks = integrator.opts.callback - integrator.finalstep = false - @assert !integrator.finalstep if isnan(integrator.dt) error("time step size `dt` is NaN") diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 28b7057090b..f4924bdcf7d 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -604,9 +604,9 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, - semi.cache) + semi.cache, semi) lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, - semi.cache) + semi.cache, semi) @test isapprox(lift, -6.501138753497174e-15, atol = 1e-13) @test isapprox(drag, 2.588589856781827, atol = 1e-13) @@ -617,13 +617,14 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_NACA0012airfoil_mach085.jl"), l2=[ - 5.371568111383228e-7, 6.4158131303956445e-6, - 1.0324346542348325e-5, 0.0006348064933187732, + 5.634402680811982e-7, 6.748066107517321e-6, + 1.091879472416885e-5, 0.0006686372064029146, ], linf=[ - 0.0016263400091978443, 0.028471072159724428, - 0.02986133204785877, 1.9481060511014872, + 0.0021456247890772823, 0.03957142889488085, + 0.03832024233032798, 2.6628739573358495, ], + amr_interval=1, base_level=0, med_level=1, max_level=1, tspan=(0.0, 0.0001), adapt_initial_condition=false, @@ -647,12 +648,12 @@ end u = Trixi.wrap_array(u_ode, semi) du = Trixi.wrap_array(du_ode, semi) drag = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, - semi.cache) + semi.cache, semi) lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, - semi.cache) + semi.cache, semi) - @test isapprox(lift, 0.0262382560809345, atol = 1e-13) - @test isapprox(drag, 0.10898248971932244, atol = 1e-13) + @test isapprox(lift, 0.029076443678087403, atol = 1e-13) + @test isapprox(drag, 0.13564720009197903, atol = 1e-13) end end diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 7749a0c4780..d038354f88a 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -739,16 +739,16 @@ end du = Trixi.wrap_array(du_ode, semi) drag_p = Trixi.analyze(drag_coefficient, du, u, tspan[2], mesh, equations, solver, - semi.cache) + semi.cache, semi) lift_p = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver, - semi.cache) + semi.cache, semi) drag_f = Trixi.analyze(drag_coefficient_shear_force, du, u, tspan[2], mesh, equations, equations_parabolic, solver, - semi.cache, semi.cache_parabolic) + semi.cache, semi, semi.cache_parabolic) lift_f = Trixi.analyze(lift_coefficient_shear_force, du, u, tspan[2], mesh, equations, equations_parabolic, solver, - semi.cache, semi.cache_parabolic) + semi.cache, semi, semi.cache_parabolic) @test isapprox(drag_p, 0.17963843913309516, atol = 1e-13) @test isapprox(lift_p, 0.26462588007949367, atol = 1e-13) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 4b3aa5c87e4..c22a8e038a7 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -256,7 +256,7 @@ end linf=[0.0015194252169410394], rtol=5.0e-5, # Higher tolerance to make tests pass in CI (in particular with macOS) elixir_file="elixir_advection_waving_flag.jl", - restart_file="restart_000021.h5", + restart_file="restart_000000021.h5", # With the default `maxiters = 1` in coverage tests, # there would be no time steps after the restart. coverage_override=(maxiters = 100_000,)) @@ -275,7 +275,7 @@ end l2=[7.841217436552029e-15], linf=[1.0857981180834031e-13], elixir_file="elixir_advection_free_stream.jl", - restart_file="restart_000036.h5", + restart_file="restart_000000036.h5", # With the default `maxiters = 1` in coverage tests, # there would be no time steps after the restart. coverage_override=(maxiters = 100_000,)) diff --git a/test/test_threaded.jl b/test/test_threaded.jl index 7365dcef21c..7fb64d61cb4 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -209,7 +209,7 @@ end linf=[0.0015194252169410394], rtol=5.0e-5, # Higher tolerance to make tests pass in CI (in particular with macOS) elixir_file="elixir_advection_waving_flag.jl", - restart_file="restart_000021.h5", + restart_file="restart_000000021.h5", # With the default `maxiters = 1` in coverage tests, # there would be no time steps after the restart. coverage_override=(maxiters = 100_000,)) diff --git a/test/test_type.jl b/test/test_type.jl index 7bea104bb86..fc9a9561e6c 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -109,10 +109,9 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT)) t = zero(RealT) - u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT)) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT)) orientation = 1 directions = [1, 2] - cons = SVector(one(RealT), one(RealT), one(RealT)) surface_flux_function = flux_lax_friedrichs @@ -204,12 +203,11 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT), zero(RealT)) t = zero(RealT) - u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), - one(RealT)) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT)) orientations = [1, 2] directions = [1, 2, 3, 4] normal_direction = SVector(one(RealT), zero(RealT)) - cons = SVector(one(RealT), one(RealT), one(RealT), one(RealT)) surface_flux_function = flux_lax_friedrichs flux_lmars = FluxLMARS(340) @@ -380,13 +378,11 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT), zero(RealT), zero(RealT)) t = zero(RealT) - u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), - one(RealT), one(RealT)) + u = u_ll = u_rr = u_inner = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), one(RealT)) orientations = [1, 2, 3] directions = [1, 2, 3, 4, 5, 6] normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) - cons = SVector(one(RealT), one(RealT), one(RealT), - one(RealT), one(RealT)) surface_flux_function = flux_lax_friedrichs flux_lmars = FluxLMARS(340) @@ -685,6 +681,918 @@ isdir(outdir) && rm(outdir, recursive = true) @test typeof(@inferred density_pressure(u, equations)) == RealT end end + + @timed_testset "Compressible Navier Stokes Diffusion 1D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerEquations1D(RealT(1.4)) + prandtl_number = RealT(0.72) + mu = RealT(0.01) + equations_parabolic_primitive = @inferred CompressibleNavierStokesDiffusion1D(equations, + mu = mu, + Prandtl = prandtl_number, + gradient_variables = GradientVariablesPrimitive()) + equations_parabolic_entropy = @inferred CompressibleNavierStokesDiffusion1D(equations, + mu = mu, + Prandtl = prandtl_number, + gradient_variables = GradientVariablesEntropy()) + + u = u_transformed = SVector(one(RealT), zero(RealT), + zero(RealT)) + orientation = 1 + gradients = SVector(RealT(0.1), RealT(0.1), RealT(0.1)) + + for equations_parabolic in (equations_parabolic_primitive, + equations_parabolic_entropy) + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == + RealT + + @test eltype(@inferred cons2prim(u, equations_parabolic)) == RealT + @test eltype(@inferred prim2cons(u, equations_parabolic)) == RealT + @test eltype(@inferred cons2entropy(u, equations_parabolic)) == RealT + @test eltype(@inferred entropy2cons(u, equations_parabolic)) == RealT + @test typeof(@inferred Trixi.temperature(u, equations_parabolic)) == RealT + + @test eltype(@inferred Trixi.convert_transformed_to_primitive(u_transformed, + equations_parabolic)) == + RealT + @test eltype(@inferred Trixi.convert_derivative_to_primitive(u, gradients, + equations_parabolic)) == + RealT + end + + # TODO: BC tests for GradientVariablesPrimitive + # TODO: BC tests for GradientVariablesEntropy + end + end + + @timed_testset "Compressible Navier Stokes Diffusion 2D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerEquations2D(RealT(1.4)) + prandtl_number = RealT(0.72) + mu = RealT(0.01) + equations_parabolic_primitive = @inferred CompressibleNavierStokesDiffusion2D(equations, + mu = mu, + Prandtl = prandtl_number, + gradient_variables = GradientVariablesPrimitive()) + equations_parabolic_entropy = @inferred CompressibleNavierStokesDiffusion2D(equations, + mu = mu, + Prandtl = prandtl_number, + gradient_variables = GradientVariablesEntropy()) + + u = u_transformed = SVector(one(RealT), zero(RealT), zero(RealT), zero(RealT)) + orientations = [1, 2] + gradient = SVector(RealT(0.1), RealT(0.1), RealT(0.1), RealT(0.1)) + gradients = SVector(gradient, gradient) + + for equations_parabolic in (equations_parabolic_primitive, + equations_parabolic_entropy) + for orientation in orientations + @test eltype(@inferred flux(u, gradients, orientation, + equations_parabolic)) == RealT + end + + @test eltype(@inferred cons2prim(u, equations_parabolic)) == RealT + @test eltype(@inferred prim2cons(u, equations_parabolic)) == RealT + @test eltype(@inferred cons2entropy(u, equations_parabolic)) == RealT + @test eltype(@inferred entropy2cons(u, equations_parabolic)) == RealT + @test typeof(@inferred Trixi.temperature(u, equations_parabolic)) == RealT + @test typeof(@inferred Trixi.enstrophy(u, gradients, equations_parabolic)) == + RealT + @test typeof(@inferred Trixi.vorticity(u, gradients, equations_parabolic)) == + RealT + + @test eltype(@inferred Trixi.convert_transformed_to_primitive(u_transformed, + equations_parabolic)) == + RealT + @test eltype(@inferred Trixi.convert_derivative_to_primitive(u, gradient, + equations_parabolic)) == + RealT + end + + # TODO: BC tests for GradientVariablesPrimitive + # TODO: BC tests for GradientVariablesEntropy + end + end + + @timed_testset "Compressible Navier Stokes Diffusion 3D" begin + for RealT in (Float32, Float64) + equations = @inferred CompressibleEulerEquations3D(RealT(1.4)) + prandtl_number = RealT(0.72) + mu = RealT(0.01) + equations_parabolic_primitive = @inferred CompressibleNavierStokesDiffusion3D(equations, + mu = mu, + Prandtl = prandtl_number, + gradient_variables = GradientVariablesPrimitive()) + equations_parabolic_entropy = @inferred CompressibleNavierStokesDiffusion3D(equations, + mu = mu, + Prandtl = prandtl_number, + gradient_variables = GradientVariablesEntropy()) + + u = u_transformed = SVector(one(RealT), zero(RealT), zero(RealT), zero(RealT), + zero(RealT)) + orientations = [1, 2, 3] + gradient = SVector(RealT(0.1), RealT(0.1), RealT(0.1), RealT(0.1), RealT(0.1)) + gradients = SVector(gradient, gradient, gradient) + + for equations_parabolic in (equations_parabolic_primitive, + equations_parabolic_entropy) + for orientation in orientations + @test eltype(@inferred flux(u, gradients, orientation, + equations_parabolic)) == RealT + end + + @test eltype(@inferred cons2prim(u, equations_parabolic)) == RealT + @test eltype(@inferred prim2cons(u, equations_parabolic)) == RealT + @test eltype(@inferred cons2entropy(u, equations_parabolic)) == RealT + @test eltype(@inferred entropy2cons(u, equations_parabolic)) == RealT + @test typeof(@inferred Trixi.temperature(u, equations_parabolic)) == RealT + @test typeof(@inferred Trixi.enstrophy(u, gradients, equations_parabolic)) == + RealT + @test eltype(@inferred Trixi.vorticity(u, gradients, equations_parabolic)) == + RealT + + @test eltype(@inferred Trixi.convert_transformed_to_primitive(u_transformed, + equations_parabolic)) == + RealT + @test eltype(@inferred Trixi.convert_derivative_to_primitive(u, gradient, + equations_parabolic)) == + RealT + end + + # TODO: BC tests for GradientVariablesPrimitive + # TODO: BC tests for GradientVariablesEntropy + end + end + + @timed_testset "Hyperbolic Diffusion 1D" begin + for RealT in (Float32, Float64) + nu = one(RealT) + Lr = RealT(inv(2pi)) + equations = @inferred HyperbolicDiffusionEquations1D(nu = nu, Lr = Lr) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = du = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT)) + orientation = 1 + directions = [1, 2] + + surface_flux_function = flux_lax_friedrichs + + @test typeof(@inferred Trixi.residual_steady_state(du, equations)) == RealT + @test eltype(@inferred initial_condition_poisson_nonperiodic(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_poisson_nonperiodic(u, x, t, equations)) == + RealT + @test eltype(@inferred source_terms_harmonic(u, x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_eoc_test_coupled_euler_gravity(x, + t, + equations)) == + RealT + + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + + @timed_testset "Hyperbolic Diffusion 2D" begin + for RealT in (Float32, Float64) + nu = one(RealT) + Lr = RealT(inv(2pi)) + equations = @inferred HyperbolicDiffusionEquations2D(nu = nu, Lr = Lr) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = du = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = SVector(one(RealT), zero(RealT)) + + surface_flux_function = flux_lax_friedrichs + + @test typeof(@inferred Trixi.residual_steady_state(du, equations)) == RealT + @test eltype(@inferred initial_condition_poisson_nonperiodic(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_poisson_nonperiodic(u, x, t, equations)) == + RealT + @test eltype(@inferred source_terms_harmonic(u, x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_eoc_test_coupled_euler_gravity(x, + t, + equations)) == + RealT + + for orientation in orientations + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end + end + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, + equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + + @timed_testset "Hyperbolic Diffusion 3D" begin + for RealT in (Float32, Float64) + nu = one(RealT) + Lr = RealT(inv(2pi)) + equations = @inferred HyperbolicDiffusionEquations3D(nu = nu, Lr = Lr) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = du = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT), + one(RealT)) + orientations = [1, 2, 3] + directions = [1, 2, 3, 4, 5, 6] + + surface_flux_function = flux_lax_friedrichs + + @test typeof(@inferred Trixi.residual_steady_state(du, equations)) == RealT + @test eltype(@inferred initial_condition_poisson_nonperiodic(x, t, equations)) == + RealT + @test eltype(@inferred source_terms_poisson_nonperiodic(u, x, t, equations)) == + RealT + @test eltype(@inferred source_terms_harmonic(u, x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_eoc_test_coupled_euler_gravity(x, + t, + equations)) == + RealT + + for orientation in orientations + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred boundary_condition_poisson_nonperiodic(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + end + end + end + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, + equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + + @timed_testset "Ideal Glm Mhd 1D" begin + for RealT in (Float32, Float64) + equations = @inferred IdealGlmMhdEquations1D(RealT(1.4)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), zero(RealT), zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT)) + orientation = 1 + directions = [1, 2] + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_hllc(u_ll, u_rr, orientation, equations)) == RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == + RealT + @test typeof(@inferred Trixi.max_abs_speeds(u, equations)) == + RealT + @test eltype(@inferred cons2prim(u, equations)) == + RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + + for direction in directions + @test typeof(Trixi.calc_fast_wavespeed(cons, direction, equations)) == RealT + @test eltype(Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, direction, + equations)) == RealT + end + + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_magnetic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred cross_helicity(cons, equations)) == RealT + end + end + + @timed_testset "Ideal Glm Mhd 2D" begin + for RealT in (Float32, Float64) + equations = @inferred IdealGlmMhdEquations2D(RealT(1.4)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), zero(RealT), zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), + zero(RealT)) + nonconservative_type_local = Trixi.NonConservativeLocal() + nonconservative_type_symmetric = Trixi.NonConservativeSymmetric() + nonconservative_terms = [1, 2] + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + normal_direction, + equations)) == RealT + else + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell_local_symmetric(u_ll, + u_rr, + orientation, + equations)) == + RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == + RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, orientation, + equations)) == RealT + end + for nonconservative_term in nonconservative_terms + @test eltype(@inferred flux_nonconservative_powell_local_symmetric(u_ll, + orientation, + equations, + nonconservative_type_local, + nonconservative_term)) == + RealT + @test eltype(@inferred flux_nonconservative_powell_local_symmetric(u_ll, + u_rr, + orientation, + equations, + nonconservative_type_symmetric, + nonconservative_term)) == + RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred entropy2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test eltype(@inferred Trixi.gradient_conservative(pressure, u, equations)) == + RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, normal_direction, + equations)) == RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + normal_direction, + equations)) == RealT + for orientation in orientations + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, orientation, + equations)) == + RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + orientation, + equations)) == RealT + end + + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_magnetic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred cross_helicity(cons, equations)) == RealT + end + end + + @timed_testset "Ideal Glm Mhd 3D" begin + for RealT in (Float32, Float64) + equations = @inferred IdealGlmMhdEquations3D(RealT(1.4)) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), zero(RealT), zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT), + zero(RealT)) + orientations = [1, 2, 3] + directions = [1, 2, 3, 4, 5, 6] + normal_direction = normal_direction_ll = normal_direction_average = SVector(one(RealT), + zero(RealT), + zero(RealT)) + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations)) == RealT + if RealT == Float32 + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + normal_direction, + equations)) == RealT + else + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_powell(u_ll, u_rr, orientation, + equations)) == RealT + if RealT == Float32 + # check `ln_mean` (test broken) + @test_broken eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == + RealT + # check `ln_mean` and `inv_ln_mean` (test broken) + @test_broken eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, + orientation, + equations)) == + RealT + else + @test eltype(@inferred flux_derigs_etal(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test typeof(@inferred pressure(u, equations)) == RealT + @test typeof(@inferred density_pressure(u, equations)) == RealT + + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, normal_direction, + equations)) == RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + normal_direction, + equations)) == RealT + for orientation in orientations + @test typeof(@inferred Trixi.calc_fast_wavespeed(cons, orientation, + equations)) == + RealT + @test eltype(@inferred Trixi.calc_fast_wavespeed_roe(u_ll, u_rr, + orientation, + equations)) == RealT + end + + @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT + @test typeof(@inferred Trixi.entropy_math(cons, equations)) == RealT + @test typeof(@inferred entropy(cons, equations)) == RealT + @test typeof(@inferred energy_total(cons, equations)) == RealT + @test typeof(@inferred energy_kinetic(cons, equations)) == RealT + @test typeof(@inferred energy_magnetic(cons, equations)) == RealT + @test typeof(@inferred energy_internal(cons, equations)) == RealT + @test typeof(@inferred cross_helicity(cons, equations)) == RealT + end + end + + @timed_testset "Linear Scalar Advection 1D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT)) + orientation = 1 + directions = [1, 2] + + surface_flux_function = flux_lax_friedrichs + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_gauss(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_sin(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_linear_x(x, t, equations)) == + RealT + + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, x, + t, + surface_flux_function, + equations)) == + RealT + end + end + + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) == + RealT + @test eltype(@inferred Trixi.flux_engquist_osher(u_ll, u_rr, orientation, + equations)) == RealT + + @test eltype(eltype(@inferred splitting_lax_friedrichs(u, orientation, + equations))) == + RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + + @timed_testset "Linear Scalar Advection 2D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation2D(RealT(1), RealT(1)) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT)) + orientations = [1, 2] + directions = [1, 2, 3, 4] + normal_direction = SVector(one(RealT), zero(RealT)) + + surface_flux_function = flux_lax_friedrichs + + @test eltype(@inferred Trixi.x_trans_periodic_2d(x)) == RealT + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_gauss(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_sin_sin(x, t, equations)) == + RealT + @test eltype(@inferred Trixi.initial_condition_linear_x_y(x, t, equations)) == + RealT + @test eltype(@inferred Trixi.initial_condition_linear_x(x, t, equations)) == + RealT + @test eltype(@inferred Trixi.initial_condition_linear_y(x, t, equations)) == + RealT + + for orientation in orientations + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test_broken eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test_broken eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_x_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_x(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + @test eltype(@inferred Trixi.boundary_condition_linear_y(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + end + end + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == + RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) == + RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end + + @timed_testset "Linear Scalar Advection 3D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation3D(RealT(1), RealT(1), + RealT(1)) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = u_inner = SVector(one(RealT)) + orientations = [1, 2, 3] + directions = [1, 2, 3, 4, 5, 6] + normal_direction = SVector(one(RealT), zero(RealT), zero(RealT)) + + surface_flux_function = flux_lax_friedrichs + + @test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT + @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) == + RealT + @test eltype(@inferred initial_condition_gauss(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_sin(x, t, equations)) == RealT + @test eltype(@inferred Trixi.initial_condition_linear_z(x, t, equations)) == + RealT + + for orientation in orientations + for direction in directions + if RealT == Float32 + # check `surface_flux_function` (test broken) + @test_broken eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + else + @test eltype(@inferred Trixi.boundary_condition_linear_z(u_inner, + orientation, + direction, + x, + t, + surface_flux_function, + equations)) == + RealT + end + end + end + + @test eltype(@inferred flux(u, normal_direction, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, normal_direction, equations)) == + RealT + + @test eltype(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, + equations)) == RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) == + RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred entropy(u, equations)) == RealT + @test typeof(@inferred energy_total(u, equations)) == RealT + end + end end end # module diff --git a/test/test_visualization.jl b/test/test_visualization.jl index 6444dc91d5d..5c7e5dbbd1f 100644 --- a/test/test_visualization.jl +++ b/test/test_visualization.jl @@ -335,9 +335,9 @@ end tspan = (0.0, 3.0)) @testset "elixir_advection_amr_visualization.jl with save_plot" begin - @test isfile(joinpath(outdir, "solution_000000.png")) - @test isfile(joinpath(outdir, "solution_000020.png")) - @test isfile(joinpath(outdir, "solution_000022.png")) + @test isfile(joinpath(outdir, "solution_000000000.png")) + @test isfile(joinpath(outdir, "solution_000000020.png")) + @test isfile(joinpath(outdir, "solution_000000022.png")) end @testset "show" begin diff --git a/utils/plot.gp b/utils/plot.gp index f58cfcdf608..497b08cd8d8 100644 --- a/utils/plot.gp +++ b/utils/plot.gp @@ -2,7 +2,7 @@ set term pdfcairo color enhanced font ",8" fontscale 1.0 lw 0.5 size 14cm,10cm # Set defaults -if (!exists("infile")) infile="solution_000000.txt" +if (!exists("infile")) infile="solution_000000000.txt" if (!exists("eqn")) eqn="linear_scalar_advection" ext = ".pdf"