From 1c1d5ea0f7fc763dd3954121fe665c6001e21df1 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Mon, 19 Aug 2024 20:18:01 -1000 Subject: [PATCH 1/6] Add numerical support of other real types (`laplace`) (#2025) * start * reformat * Update test/test_type.jl Co-authored-by: Hendrik Ranocha * Update test/test_type.jl Co-authored-by: Hendrik Ranocha * Update test/test_type.jl Co-authored-by: Hendrik Ranocha --------- Co-authored-by: Hendrik Ranocha --- test/test_type.jl | 63 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index c15ac5f78c..d13b626b06 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1459,6 +1459,69 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Laplace Diffusion 1D" begin + for RealT in (Float32, Float64) + equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) + equations_parabolic = @inferred LaplaceDiffusion1D(RealT(0.1), equations) + + x = SVector(zero(RealT)) + t = zero(RealT) + u = gradients = SVector(one(RealT)) + orientation = 1 + + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == + RealT + + # TODO: BC tests for BoundaryConditionDirichlet + # TODO: BC tests for BoundaryConditionNeumann + end + end + + @timed_testset "Laplace Diffusion 2D" begin + for RealT in (Float32, Float64) + equations = LinearScalarAdvectionEquation2D(RealT(1), RealT(1)) + equations_parabolic = LaplaceDiffusion2D(RealT(0.1), equations) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_inner = u_outer = inv_h = gradients = SVector(one(RealT), one(RealT)) + orientations = [1, 2] + + for orientation in orientations + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == + RealT + end + + parabolic_solver = ViscousFormulationLocalDG(RealT(0.1)) + @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h, + equations_parabolic, parabolic_solver)) == + RealT + end + end + + @timed_testset "Laplace Diffusion 3D" begin + for RealT in (Float32, Float64) + equations = LinearScalarAdvectionEquation3D(RealT(1), RealT(1), RealT(1)) + equations_parabolic = LaplaceDiffusion3D(RealT(0.1), equations) + + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_inner = u_outer = inv_h = gradients = SVector(one(RealT), one(RealT), + one(RealT)) + orientations = [1, 2, 3] + + for orientation in orientations + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == + RealT + end + + parabolic_solver = ViscousFormulationLocalDG(RealT(0.1)) + @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h, + equations_parabolic, parabolic_solver)) == + RealT + end + end + @timed_testset "Lattice Boltzmann 2D" begin for RealT in (Float32, Float64) equations = @inferred LatticeBoltzmannEquations2D(Ma = RealT(0.1), Re = 1000) From e47afa1048f43ac0bad93710afaabcedfb3989d0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 20 Aug 2024 13:41:33 +0200 Subject: [PATCH 2/6] set version to v0.8.8 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 923cf5248a..3283a41c1a 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.8.8-DEV" +version = "0.8.8" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From afd80605ef7e7f4d36e3fc12ddb0280a8db38bfd Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 20 Aug 2024 13:41:47 +0200 Subject: [PATCH 3/6] set development version to v0.8.9-DEV --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3283a41c1a..f8be84e0ce 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.8.8" +version = "0.8.9-DEV" [deps] CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" From 20ab97bea82fbd4c286f8fc4cebad6ee12a1e465 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 21 Aug 2024 09:51:23 +0200 Subject: [PATCH 4/6] computations using only `Float32` and no `Float64` (#2042) * add example using only Float32 * make RHS work for examples/unstructured_2d_dgsem/elixir_shallowwater_ec_float32.jl * more fixes * format * more fixes * examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries_float32.jl * add tests * upwind FDSBP * format * format * fix reference values * adapt SWE EC test value after PR #2038 * more mesh fixes * more fixes * Apply suggestions from code review Co-authored-by: Andrew Winters * remove redundant using Downloads in examples * more comments * format --------- Co-authored-by: Andrew Winters --- .../elixir_euler_NACA0012airfoil_mach085.jl | 2 +- .../elixir_euler_NACA6412airfoil_mach2.jl | 9 +- .../elixir_euler_subsonic_cylinder.jl | 2 +- ...xir_navierstokes_NACA0012airfoil_mach08.jl | 4 +- .../elixir_euler_free_stream_boundaries.jl | 5 +- ...ir_euler_free_stream_boundaries_float32.jl | 62 +++++++++ .../elixir_advection_float32.jl | 61 +++++++++ .../elixir_shallowwater_ec_float32.jl | 124 ++++++++++++++++++ ...elixir_euler_free_stream_upwind_float32.jl | 88 +++++++++++++ src/auxiliary/auxiliary.jl | 2 +- src/meshes/p4est_mesh.jl | 14 +- src/meshes/structured_mesh.jl | 8 +- src/meshes/surface_interpolant.jl | 8 +- src/meshes/t8code_mesh.jl | 4 +- src/meshes/transfinite_mappings_3d.jl | 6 +- .../semidiscretization_coupled.jl | 4 +- src/solvers/dgsem_p4est/dg_2d.jl | 14 +- src/solvers/dgsem_p4est/dg_2d_parabolic.jl | 58 ++++---- src/solvers/dgsem_p4est/dg_3d.jl | 14 +- src/solvers/dgsem_p4est/dg_3d_parabolic.jl | 34 ++--- src/solvers/dgsem_structured/containers_3d.jl | 12 +- src/solvers/dgsem_structured/dg_2d.jl | 26 ++-- .../dg_2d_subcell_limiters.jl | 4 +- src/solvers/dgsem_structured/dg_3d.jl | 34 ++--- src/solvers/dgsem_structured/indicators_1d.jl | 5 +- src/solvers/dgsem_structured/indicators_2d.jl | 10 +- src/solvers/dgsem_structured/indicators_3d.jl | 15 ++- src/solvers/dgsem_t8code/containers_3d.jl | 6 +- src/solvers/dgsem_tree/dg_1d.jl | 16 +-- src/solvers/dgsem_tree/dg_1d_parabolic.jl | 4 +- src/solvers/dgsem_tree/dg_2d.jl | 28 ++-- src/solvers/dgsem_tree/dg_2d_parabolic.jl | 10 +- .../dgsem_tree/dg_2d_subcell_limiters.jl | 4 +- src/solvers/dgsem_tree/dg_3d.jl | 42 +++--- src/solvers/dgsem_tree/dg_3d_parabolic.jl | 10 +- src/solvers/dgsem_tree/indicators.jl | 7 +- src/solvers/dgsem_tree/indicators_1d.jl | 4 +- src/solvers/dgsem_tree/indicators_2d.jl | 12 +- src/solvers/dgsem_tree/indicators_3d.jl | 24 ++-- src/solvers/dgsem_tree/subcell_limiters_2d.jl | 2 +- src/solvers/dgsem_unstructured/dg_2d.jl | 9 +- .../dgsem_unstructured/indicators_2d.jl | 4 +- .../mappings_geometry_curved_2d.jl | 62 +++++---- .../mappings_geometry_straight_2d.jl | 40 +++--- src/visualization/utilities.jl | 2 +- test/test_p4est_3d.jl | 29 ++++ test/test_structured_2d.jl | 17 +++ test/test_trixi.jl | 19 ++- test/test_unstructured_2d.jl | 44 +++++++ 49 files changed, 737 insertions(+), 287 deletions(-) create mode 100644 examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries_float32.jl create mode 100644 examples/structured_2d_dgsem/elixir_advection_float32.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_ec_float32.jl create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind_float32.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl index 8806f9cc2a..45b750728c 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl @@ -1,4 +1,4 @@ -using Downloads: download + using OrdinaryDiffEq using Trixi diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl index 7e55a25959..573d182e15 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl @@ -1,7 +1,6 @@ using Trixi using OrdinaryDiffEq -using Downloads: download ############################################################################### # semidiscretization of the compressible Euler equations @@ -62,7 +61,7 @@ volume_integral = VolumeIntegralShockCapturingHG(shock_indicator; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) -# DG Solver +# DG Solver solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, volume_integral = volume_integral) @@ -70,8 +69,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, # https://gist.githubusercontent.com/DanielDoehring/5ade6d93629f0d8c23a598812dbee2a9/raw/d2bc904fe92146eae1a36156e7f5c535dc1a80f1/NACA6412.geo mesh_file = joinpath(@__DIR__, "mesh_NACA6412.inp") isfile(mesh_file) || - download("https://gist.githubusercontent.com/DanielDoehring/e2a389f04f1e37b33819b9637e8ee4c3/raw/4bf7607a2ce4432fdb5cb87d5e264949b11bd5d7/mesh_NACA6412.inp", - mesh_file) + Trixi.download("https://gist.githubusercontent.com/DanielDoehring/e2a389f04f1e37b33819b9637e8ee4c3/raw/4bf7607a2ce4432fdb5cb87d5e264949b11bd5d7/mesh_NACA6412.inp", + mesh_file) boundary_symbols = [:PhysicalLine1, :PhysicalLine2, :PhysicalLine3, :PhysicalLine4] @@ -79,7 +78,7 @@ mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_sy boundary_conditions = Dict(:PhysicalLine1 => boundary_condition_supersonic_inflow, # Left boundary :PhysicalLine2 => boundary_condition_supersonic_outflow, # Right boundary - :PhysicalLine3 => boundary_condition_supersonic_outflow, # Top and bottom boundary + :PhysicalLine3 => boundary_condition_supersonic_outflow, # Top and bottom boundary :PhysicalLine4 => boundary_condition_slip_wall) # Airfoil semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl index e522a0dcfd..21df724da4 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl @@ -1,4 +1,4 @@ -using Downloads: download + using OrdinaryDiffEq using Trixi diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl index 3e4679f5af..40baef6ef9 100644 --- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl +++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl @@ -1,4 +1,4 @@ -using Downloads: download + using OrdinaryDiffEq using Trixi @@ -15,7 +15,7 @@ using Trixi # Structured and Unstructured Grid Methods (2016) # [https://ntrs.nasa.gov/citations/20160003623] (https://ntrs.nasa.gov/citations/20160003623) # - Deep Ray, Praveen Chandrashekar (2017) -# An entropy stable finite volume scheme for the +# An entropy stable finite volume scheme for the # two dimensional Navier–Stokes equations on triangular grids # [DOI:10.1016/j.amc.2017.07.020](https://doi.org/10.1016/j.amc.2017.07.020) diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl index bdc4da26c1..8c008fac2e 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -18,8 +17,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) default_mesh_file = joinpath(@__DIR__, "mesh_cube_with_boundaries.inp") isfile(default_mesh_file) || - download("https://gist.githubusercontent.com/DanielDoehring/710eab379fe3042dc08af6f2d1076e49/raw/38e9803bc0dab9b32a61d9542feac5343c3e6f4b/mesh_cube_with_boundaries.inp", - default_mesh_file) + Trixi.download("https://gist.githubusercontent.com/DanielDoehring/710eab379fe3042dc08af6f2d1076e49/raw/38e9803bc0dab9b32a61d9542feac5343c3e6f4b/mesh_cube_with_boundaries.inp", + default_mesh_file) mesh_file = default_mesh_file boundary_symbols = [:PhysicalSurface1, :PhysicalSurface2] diff --git a/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries_float32.jl b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries_float32.jl new file mode 100644 index 0000000000..c41ad9a422 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_free_stream_boundaries_float32.jl @@ -0,0 +1,62 @@ +# Similar to p4est_3d_dgsem/elixir_euler_free_stream_boundaries.jl +# but using Float32 instead of the default Float64 + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4f0) + +initial_condition = initial_condition_constant + +polydeg = 3 +solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs, RealT = Float32) + +############################################################################### +# Get the uncurved mesh from a file (downloads the file if not available locally) + +default_mesh_file = joinpath(@__DIR__, "mesh_cube_with_boundaries.inp") +isfile(default_mesh_file) || + Trixi.download("https://gist.githubusercontent.com/DanielDoehring/710eab379fe3042dc08af6f2d1076e49/raw/38e9803bc0dab9b32a61d9542feac5343c3e6f4b/mesh_cube_with_boundaries.inp", + default_mesh_file) +mesh_file = default_mesh_file + +boundary_symbols = [:PhysicalSurface1, :PhysicalSurface2] + +mesh = P4estMesh{3}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols, + RealT = Float32) + +boundary_conditions = Dict(:PhysicalSurface1 => BoundaryConditionDirichlet(initial_condition), + :PhysicalSurface2 => BoundaryConditionDirichlet(initial_condition)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0f0, 1.0f0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.5f0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0f0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_advection_float32.jl b/examples/structured_2d_dgsem/elixir_advection_float32.jl new file mode 100644 index 0000000000..20c30d4aeb --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_advection_float32.jl @@ -0,0 +1,61 @@ +# Similar to structured_2d_dgsem/elixir_advection_basic.jl +# but using Float32 instead of the default Float64 + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2f0, -0.7f0) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, RealT = Float32) + +coordinates_min = (-1.0f0, -1.0f0) # minimum coordinates (min(x), min(y)) +coordinates_max = (1.0f0, 1.0f0) # maximum coordinates (max(x), max(y)) + +cells_per_dimension = (16, 16) + +# Create curved mesh with 16 x 16 elements +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, (0.0f0, 1.0f0)) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 100, + solution_variables = cons2prim) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 1.6f0) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0f0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_float32.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_float32.jl new file mode 100644 index 0000000000..cc527fc00e --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_float32.jl @@ -0,0 +1,124 @@ +# Similar to unstructured_2d_dgsem/elixir_shallowwater_ec_float32.jl +# but using Float32 instead of the default Float64 + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function + +equations = ShallowWaterEquations2D(gravity_constant = 9.81f0) + +# Note, this initial condition is used to compute errors in the analysis callback but the initialization is +# overwritten by `initial_condition_ec_discontinuous_bottom` below. +initial_condition = initial_condition_weak_blast_wave + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 6, + surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux), + RealT = Float32) + +############################################################################### +# This setup is for the curved, split form entropy conservation testing (needs periodic BCs) + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true, RealT = Float32) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0f0, 2.0f0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function and initial condition for this academic testcase of entropy conservation. +# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! +function initial_condition_ec_discontinuous_bottom(x, t, element_id, + equations::ShallowWaterEquations2D) + RealT = eltype(x) + + # Set up polar coordinates + inicenter = SVector{2, RealT}(0.7, 0.7) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Set the background values + H = 3.25f0 + v1 = zero(RealT) + v2 = zero(RealT) + b = zero(RealT) + + # setup the discontinuous water height and velocities + if element_id == 10 + H = 4.0f0 + v1 = convert(RealT, 0.1882) * cos_phi + v2 = convert(RealT, 0.1882) * sin_phi + end + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2 + 0.5f0 * sinpi(2 * x[1]) + 0.5f0 * cospi(2 * x[2]) + end + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for j in eachnode(semi.solver), i in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, + semi.solver, i, j, element) + u_node = initial_condition_ec_discontinuous_bottom(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0f0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0f0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind_float32.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind_float32.jl new file mode 100644 index 0000000000..9cfd45a4d3 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind_float32.jl @@ -0,0 +1,88 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. +# Similar to unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl +# but using Float32 instead of the default Float64 + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4f0) + +initial_condition = initial_condition_constant + +# Boundary conditions for free-stream preservation test +boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition) + +boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream, + :cone1 => boundary_condition_free_stream, + :cone2 => boundary_condition_free_stream, + :iceCream => boundary_condition_free_stream) + +############################################################################### +# Get the Upwind FDSBP approximation space + +# TODO: FDSBP +# Note, one must set `xmin=-1` and `xmax=1` due to the reuse +# of interpolation routines from `calc_node_coordinates!` to create +# the physical coordinates in the mappings. +D_upw = upwind_operators(SummationByPartsOperators.Mattsson2017, + derivative_order = 1, + accuracy_order = 8, + xmin = -1.0f0, xmax = 1.0f0, + N = 17) + +flux_splitting = splitting_vanleer_haenel +solver = FDSBP(D_upw, + surface_integral = SurfaceIntegralStrongForm(FluxUpwind(flux_splitting)), + volume_integral = VolumeIntegralUpwind(flux_splitting)) + +############################################################################### +# Get the curved quad mesh from a file (downloads the file if not available locally) + +# Mesh with second-order boundary polynomials requires an upwind SBP operator +# with (at least) 4th order boundary closure to guarantee the approximation is +# free-stream preserving +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/ec9a345f09199ebe471d35d5c1e4e08f/raw/15975943d8642e42f8292235314b6f1b30aa860d/mesh_inner_outer_boundaries.mesh", + joinpath(@__DIR__, "mesh_inner_outer_boundaries.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, RealT = Float32) + +############################################################################### +# create the semi discretization object + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0f0, 5.0f0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + save_solution, + alive_callback) + +############################################################################### +# run the simulation + +# set small tolerances for the free-stream preservation test +sol = solve(ode, SSPRK43(), abstol = 1.0f-6, reltol = 1.0f-6, + save_everystep = false, callback = callbacks) + +summary_callback() # print the timer summary diff --git a/src/auxiliary/auxiliary.jl b/src/auxiliary/auxiliary.jl index 97263405d2..8ab798dd33 100644 --- a/src/auxiliary/auxiliary.jl +++ b/src/auxiliary/auxiliary.jl @@ -207,7 +207,7 @@ macro threaded(expr) # !!! danger "Heisenbug" # Look at the comments for `wrap_array` when considering to change this macro. expr = if _PREFERENCE_POLYESTER - # Currently using `@batch` from Polyester.jl is more efficient, + # Currently using `@batch` from Polyester.jl is more efficient, # bypasses the Julia task scheduler and provides parallelization with less overhead. quote $Trixi.@batch $(expr) diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl index 2046220aec..526f5d9f23 100644 --- a/src/meshes/p4est_mesh.jl +++ b/src/meshes/p4est_mesh.jl @@ -349,7 +349,7 @@ For example, if a two-dimensional base mesh contains 25 elements then setting independent of domain partitioning. Should be `false` for static meshes to permit more fine-grained partitioning. - `boundary_symbols::Vector{Symbol}`: A vector of symbols that correspond to the boundary names in the `meshfile`. - If `nothing` is passed then all boundaries are named `:all`. + If `nothing` is passed then all boundaries are named `:all`. """ function P4estMesh{NDIMS}(meshfile::String; mapping = nothing, polydeg = 1, RealT = Float64, @@ -534,7 +534,7 @@ function parse_elements(meshfile, n_trees, n_dims) element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", "*ELEMENT, type=S4"] : ["*ELEMENT, type=C3D8"] - # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index + # 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index expected_content_length = n_dims == 2 ? 5 : 9 element_node_matrix = Matrix{Int64}(undef, n_trees, expected_content_length - 1) @@ -617,7 +617,7 @@ function assign_boundaries_standard_abaqus!(boundary_names, n_trees, ::Val{2}) # 2D version for tree in 1:n_trees tree_nodes = element_node_matrix[tree, :] - # For node labeling, see + # For node labeling, see # https://docs.software.vt.edu/abaqusv2022/English/SIMACAEELMRefMap/simaelm-r-2delem.htm#simaelm-r-2delem-t-nodedef1 # and search for "Node ordering and face numbering on elements" for boundary in keys(node_set_dict) # Loop over specified boundaries @@ -658,7 +658,7 @@ function assign_boundaries_standard_abaqus!(boundary_names, n_trees, ::Val{3}) # 3D version for tree in 1:n_trees tree_nodes = element_node_matrix[tree, :] - # For node labeling, see + # For node labeling, see # https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html for boundary in keys(node_set_dict) # Loop over specified boundaries # Check "front face" (y_min) @@ -1460,7 +1460,7 @@ function cubed_sphere_mapping(xi, eta, zeta, inner_radius, thickness, direction) r = sqrt(1 + x^2 + y^2) # Radius of the sphere - R = inner_radius + thickness * (0.5 * (zeta + 1)) + R = inner_radius + thickness * (0.5f0 * (zeta + 1)) # Projection onto the sphere return R / r * cube_coordinates[direction] @@ -1650,7 +1650,7 @@ function calc_tree_node_coordinates!(node_coordinates::AbstractArray{<:Any, 5}, nodes[q]) end else # curved_check[face] == 1 - # Curved face boundary information is supplied by + # Curved face boundary information is supplied by # the mesh file. Just read it into a work array for q in 1:nnodes, p in 1:nnodes file_idx += 1 @@ -1723,7 +1723,7 @@ end # and return the 3D coordinate point (x, y, z) function bilinear_interpolation!(coordinate, face_vertices, u, v) for j in 1:3 - coordinate[j] = 0.25 * (face_vertices[j, 1] * (1 - u) * (1 - v) + coordinate[j] = 0.25f0 * (face_vertices[j, 1] * (1 - u) * (1 - v) + face_vertices[j, 2] * (1 + u) * (1 - v) + face_vertices[j, 3] * (1 + u) * (1 + v) + face_vertices[j, 4] * (1 - u) * (1 + v)) diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index 9e5c55197d..155bad8464 100644 --- a/src/meshes/structured_mesh.jl +++ b/src/meshes/structured_mesh.jl @@ -160,7 +160,7 @@ end # Interpolate linearly between left and right value where s should be between -1 and 1 function linear_interpolate(s, left_value, right_value) - 0.5 * ((1 - s) * left_value + (1 + s) * right_value) + 0.5f0 * ((1 - s) * left_value + (1 + s) * right_value) end # Convert min and max coordinates of a rectangle to the corresponding transformation mapping @@ -197,7 +197,7 @@ function bilinear_mapping(x, y, faces) x3 = faces[1](1) # Top left x4 = faces[2](1) # Top right - return 0.25 * (x1 * (1 - x) * (1 - y) + + return 0.25f0 * (x1 * (1 - x) * (1 - y) + x2 * (1 + x) * (1 - y) + x3 * (1 - x) * (1 + y) + x4 * (1 + x) * (1 + y)) @@ -215,7 +215,7 @@ function trilinear_mapping(x, y, z, faces) x7 = faces[1](1, 1) # mapped from (-1, 1, 1) x8 = faces[2](1, 1) # mapped from ( 1, 1, 1) - return 0.125 * (x1 * (1 - x) * (1 - y) * (1 - z) + + return 0.125f0 * (x1 * (1 - x) * (1 - y) * (1 - z) + x2 * (1 + x) * (1 - y) * (1 - z) + x3 * (1 - x) * (1 + y) * (1 - z) + x4 * (1 + x) * (1 + y) * (1 - z) + @@ -264,7 +264,7 @@ function correction_term_3d(x, y, z, faces) # Each of the 12 edges are counted twice above # so we divide the correction term by two - return 0.5 * (c_x + c_y + c_z) + return 0.5f0 * (c_x + c_y + c_z) end # In 3D diff --git a/src/meshes/surface_interpolant.jl b/src/meshes/surface_interpolant.jl index 6a7c5c7834..20a6611668 100644 --- a/src/meshes/surface_interpolant.jl +++ b/src/meshes/surface_interpolant.jl @@ -66,8 +66,8 @@ function chebyshev_gauss_lobatto_nodes_weights(n_nodes::Integer) nodes[j] = -cospi((j - 1) / N) weights[j] = pi / N end - weights[1] = 0.5 * weights[1] - weights[end] = 0.5 * weights[end] + weights[1] = 0.5f0 * weights[1] + weights[end] = 0.5f0 * weights[end] return nodes, weights end @@ -80,7 +80,9 @@ function lagrange_interpolation(x, nodes, fvals, wbary) denominator = zero(eltype(fvals)) for j in eachindex(nodes) - if isapprox(x, nodes[j], rtol = eps(x)) + # using eps(nodes[j]) instead of eps(x) allows us to use integer + # coordinates for the target location x + if isapprox(x, nodes[j], rtol = eps(nodes[j])) return fvals[j] end t = wbary[j] / (x - nodes[j]) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 10a042be3b..aa49e641dc 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -121,7 +121,7 @@ function T8codeMesh{NDIMS, RealT}(forest::Ptr{t8_forest}, boundary_names; polyde mapping = nothing) where {NDIMS, RealT} # In t8code reference space is [0,1]. basis = LobattoLegendreBasis(RealT, polydeg) - nodes = 0.5 .* (basis.nodes .+ 1.0) + nodes = 0.5f0 .* (basis.nodes .+ 1) cmesh = t8_forest_get_cmesh(forest) number_of_trees = t8_forest_get_num_global_trees(forest) @@ -572,7 +572,7 @@ For example, if a two-dimensional base mesh contains 25 elements then setting - `RealT::Type`: The type that should be used for coordinates. - `initial_refinement_level::Integer`: Refine the mesh uniformly to this level before the simulation starts. - `boundary_symbols::Vector{Symbol}`: A vector of symbols that correspond to the boundary names in the `meshfile`. - If `nothing` is passed then all boundaries are named `:all`. + If `nothing` is passed then all boundaries are named `:all`. """ function T8codeMesh(meshfile::AbaqusFile{NDIMS}; mapping = nothing, polydeg = 1, RealT = Float64, diff --git a/src/meshes/transfinite_mappings_3d.jl b/src/meshes/transfinite_mappings_3d.jl index 59a02f33e1..a82526bc1d 100644 --- a/src/meshes/transfinite_mappings_3d.jl +++ b/src/meshes/transfinite_mappings_3d.jl @@ -52,7 +52,7 @@ function straight_side_hex_map(xi, eta, zeta, corner_points) coordinate = zeros(eltype(xi), 3) for j in 1:3 - coordinate[j] += (0.125 * + coordinate[j] += (0.125f0 * (corner_points[j, 1] * (1 - xi) * (1 - eta) * (1 - zeta) + corner_points[j, 2] * (1 + xi) * (1 - eta) * (1 - zeta) + corner_points[j, 3] * (1 + xi) * (1 + eta) * (1 - zeta) @@ -131,7 +131,7 @@ function transfinite_hex_map(xi, eta, zeta, face_curves::AbstractVector{<:Curved # Compute the transfinite mapping for j in 1:3 # Linear interpolation between opposite faces - coordinate[j] = (0.5 * + coordinate[j] = (0.5f0 * (face_values[j, 6] * (1 - xi) + face_values[j, 4] * (1 + xi) + face_values[j, 1] * (1 - eta) + face_values[j, 2] * (1 + eta) @@ -139,7 +139,7 @@ function transfinite_hex_map(xi, eta, zeta, face_curves::AbstractVector{<:Curved face_values[j, 5] * (1 + zeta))) # Edge corrections to ensure faces match - coordinate[j] -= (0.25 * (edge_values[j, 1] * (1 - eta) * (1 - zeta) + coordinate[j] -= (0.25f0 * (edge_values[j, 1] * (1 - eta) * (1 - zeta) + edge_values[j, 2] * (1 + xi) * (1 - eta) + edge_values[j, 3] * (1 - eta) * (1 + zeta) + edge_values[j, 4] * (1 - xi) * (1 - eta) diff --git a/src/semidiscretization/semidiscretization_coupled.jl b/src/semidiscretization/semidiscretization_coupled.jl index 6b009cfad2..7c1fbef972 100644 --- a/src/semidiscretization/semidiscretization_coupled.jl +++ b/src/semidiscretization/semidiscretization_coupled.jl @@ -486,13 +486,13 @@ function (boundary_condition::BoundaryConditionCoupled)(u_inner, orientation, di if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary flux = (surface_flux_function[1](u_inner, u_boundary, orientation, equations) + - 0.5 * + 0.5f0 * surface_flux_function[2](u_inner, u_boundary, orientation, equations)) else # u_boundary is "left" of boundary, u_inner is "right" of boundary flux = (surface_flux_function[1](u_boundary, u_inner, orientation, equations) + - 0.5 * + 0.5f0 * surface_flux_function[2](u_boundary, u_inner, orientation, equations)) end diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 36624f2ce8..29e94f4370 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -238,10 +238,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, primary_node_index, primary_direction_index, primary_element_index] = (flux_[v] + - 0.5 * + 0.5f0 * noncons_primary[v]) surface_flux_values[v, secondary_node_index, secondary_direction_index, secondary_element_index] = -(flux_[v] + - 0.5 * + 0.5f0 * noncons_secondary[v]) end end @@ -501,10 +501,10 @@ function calc_mortar_flux!(surface_flux_values, # copying in the correct orientation u_buffer = cache.u_threaded[Threads.threadid()] - # in calc_interface_flux!, the interface flux is computed once over each - # interface using the normal from the "primary" element. The result is then - # passed back to the "secondary" element, flipping the sign to account for the - # change in the normal direction. For mortars, this sign flip occurs in + # in calc_interface_flux!, the interface flux is computed once over each + # interface using the normal from the "primary" element. The result is then + # passed back to the "secondary" element, flipping the sign to account for the + # change in the normal direction. For mortars, this sign flip occurs in # "mortar_fluxes_to_elements!" instead. mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, @@ -557,7 +557,7 @@ end noncons = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations) - flux_plus_noncons = flux + 0.5 * noncons + flux_plus_noncons = flux + 0.5f0 * noncons # Copy to buffer set_node_vars!(fstar[position_index], flux_plus_noncons, equations, dg, node_index) diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index ed21f37144..4bf4c79ebe 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -28,8 +28,8 @@ function create_cache_parabolic(mesh::P4estMesh{2}, end #= -Reusing `rhs_parabolic!` for `TreeMesh`es is not easily possible as -for `P4estMesh`es we call +Reusing `rhs_parabolic!` for `TreeMesh`es is not easily possible as +for `P4estMesh`es we call ``` prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic, @@ -38,8 +38,8 @@ for `P4estMesh`es we call calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values, mesh, equations_parabolic, dg.mortar, dg.surface_integral, dg, cache) - ``` -instead of + ``` +instead of ``` prolong2mortars!(cache, flux_viscous, mesh, equations_parabolic, dg.mortar, dg.surface_integral, dg) @@ -221,14 +221,14 @@ function calc_gradient!(gradients, u_transformed, t, end end - # Prolong solution to interfaces. + # Prolong solution to interfaces. # This reuses `prolong2interfaces` for the purely hyperbolic case. @trixi_timeit timer() "prolong2interfaces" begin prolong2interfaces!(cache_parabolic, u_transformed, mesh, equations_parabolic, dg.surface_integral, dg) end - # Calculate interface fluxes for the gradient. + # Calculate interface fluxes for the gradient. # This reuses `calc_interface_flux!` for the purely hyperbolic case. @trixi_timeit timer() "interface flux" begin calc_interface_flux!(cache_parabolic.elements.surface_flux_values, @@ -257,8 +257,8 @@ function calc_gradient!(gradients, u_transformed, t, end # Calculate mortar fluxes. This reuses the hyperbolic version of `calc_mortar_flux`, - # along with a specialization on `calc_mortar_flux!(fstar, ...)` and `mortar_fluxes_to_elements!` for - # AbstractEquationsParabolic. + # along with a specialization on `calc_mortar_flux!(fstar, ...)` and `mortar_fluxes_to_elements!` for + # AbstractEquationsParabolic. @trixi_timeit timer() "mortar flux" begin calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, mesh, False(), # False() = no nonconservative terms @@ -388,11 +388,11 @@ function calc_gradient!(gradients, u_transformed, t, return nothing end -# This version is called during `calc_gradients!` and must be specialized because the -# flux for the gradient is {u}, which doesn't depend on the outward normal. Thus, -# you don't need to scale by 2 (e.g., the scaling factor in the normals (and in the -# contravariant vectors) along large/small elements across a non-conforming -# interface in 2D) and flip the sign when storing the mortar fluxes back +# This version is called during `calc_gradients!` and must be specialized because the +# flux for the gradient is {u}, which doesn't depend on the outward normal. Thus, +# you don't need to scale by 2 (e.g., the scaling factor in the normals (and in the +# contravariant vectors) along large/small elements across a non-conforming +# interface in 2D) and flip the sign when storing the mortar fluxes back # into `surface_flux_values`. @inline function mortar_fluxes_to_elements!(surface_flux_values, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, @@ -462,7 +462,7 @@ end u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_node_index, interface_index) - flux_ = 0.5 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux + flux_ = 0.5f0 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux # Note that we don't flip the sign on the secondary flux. This is because for parabolic terms, # the normals are not embedded in `flux_` for the parabolic gradient computations. @@ -604,7 +604,7 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, return nothing end -# This version is used for divergence flux computations +# This version is used for divergence flux computations function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, equations_parabolic, dg::DG, cache_parabolic) @@ -645,7 +645,7 @@ function calc_interface_flux!(surface_flux_values, end for node in eachnode(dg) - # We prolong the viscous flux dotted with respect the outward normal on the + # We prolong the viscous flux dotted with respect the outward normal on the # primary element. We assume a BR-1 type of flux. viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, @@ -653,7 +653,7 @@ function calc_interface_flux!(surface_flux_values, node, interface) - flux = 0.5 * (viscous_flux_normal_ll + viscous_flux_normal_rr) + flux = 0.5f0 * (viscous_flux_normal_ll + viscous_flux_normal_rr) for v in eachvariable(equations_parabolic) surface_flux_values[v, node, primary_direction_index, primary_element] = flux[v] @@ -741,11 +741,11 @@ function prolong2mortars_divergence!(cache, flux_viscous::Vector{Array{uEltype, flux_viscous = SVector(flux_viscous_x[v, i_large, j_large, element], flux_viscous_y[v, i_large, j_large, element]) - # We prolong the viscous flux dotted with respect the outward normal - # on the small element. We scale by -1/2 here because the normal - # direction on the large element is negative 2x that of the small + # We prolong the viscous flux dotted with respect the outward normal + # on the small element. We scale by -1/2 here because the normal + # direction on the large element is negative 2x that of the small # element (these normal directions are "scaled" by the surface Jacobian) - u_buffer[v, i] = -0.5 * dot(flux_viscous, normal_direction) + u_buffer[v, i] = -0.5f0 * dot(flux_viscous, normal_direction) end i_large += i_large_step j_large += j_large_step @@ -763,8 +763,8 @@ function prolong2mortars_divergence!(cache, flux_viscous::Vector{Array{uEltype, return nothing end -# We specialize `calc_mortar_flux!` for the divergence part of -# the parabolic terms. +# We specialize `calc_mortar_flux!` for the divergence part of +# the parabolic terms. function calc_mortar_flux_divergence!(surface_flux_values, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations::AbstractEquationsParabolic, @@ -789,7 +789,7 @@ function calc_mortar_flux_divergence!(surface_flux_values, mortar] # TODO: parabolic; only BR1 at the moment - fstar[position][v, node] = 0.5 * (viscous_flux_normal_ll + + fstar[position][v, node] = 0.5f0 * (viscous_flux_normal_ll + viscous_flux_normal_rr) end end @@ -808,11 +808,11 @@ function calc_mortar_flux_divergence!(surface_flux_values, return nothing end -# We structure `calc_interface_flux!` similarly to "calc_mortar_flux!" for -# hyperbolic equations with no nonconservative terms. -# The reasoning is that parabolic fluxes are treated like conservative +# We structure `calc_interface_flux!` similarly to "calc_mortar_flux!" for +# hyperbolic equations with no nonconservative terms. +# The reasoning is that parabolic fluxes are treated like conservative # terms (e.g., we compute a viscous conservative "flux") and thus no -# non-conservative terms are present. +# non-conservative terms are present. @inline function calc_mortar_flux!(fstar, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, @@ -827,7 +827,7 @@ end mortar_index) # TODO: parabolic; only BR1 at the moment - flux_ = 0.5 * (u_ll + u_rr) + flux_ = 0.5f0 * (u_ll + u_rr) # Copy flux to buffer set_node_vars!(fstar[position_index], flux_, equations, dg, node_index) diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index 5b3c5ae5ca..ece4840b74 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/src/solvers/dgsem_p4est/dg_3d.jl @@ -305,10 +305,10 @@ end # central fluxes/SATs surface_flux_values[v, primary_i_node_index, primary_j_node_index, primary_direction_index, primary_element_index] = flux_[v] + - 0.5 * noncons_primary[v] + 0.5f0 * noncons_primary[v] surface_flux_values[v, secondary_i_node_index, secondary_j_node_index, secondary_direction_index, secondary_element_index] = -(flux_[v] + - 0.5 * + 0.5f0 * noncons_secondary[v]) end end @@ -580,10 +580,10 @@ function calc_mortar_flux!(surface_flux_values, # copying in the correct orientation u_buffer = cache.u_threaded[Threads.threadid()] - # in calc_interface_flux!, the interface flux is computed once over each - # interface using the normal from the "primary" element. The result is then - # passed back to the "secondary" element, flipping the sign to account for the - # change in the normal direction. For mortars, this sign flip occurs in + # in calc_interface_flux!, the interface flux is computed once over each + # interface using the normal from the "primary" element. The result is then + # passed back to the "secondary" element, flipping the sign to account for the + # change in the normal direction. For mortars, this sign flip occurs in # "mortar_fluxes_to_elements!" instead. mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, @@ -635,7 +635,7 @@ end # central fluxes/SATs noncons = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations) - flux_plus_noncons = flux + 0.5 * noncons + flux_plus_noncons = flux + 0.5f0 * noncons # Copy to buffer set_node_vars!(fstar, flux_plus_noncons, equations, dg, i_node_index, j_node_index, diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index 63d431d35d..3f286ca01f 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -144,7 +144,7 @@ function calc_gradient!(gradients, u_transformed, t, end # Prolong solution to mortars. These should reuse the hyperbolic version of `prolong2mortars` - # !!! NOTE: we reuse the hyperbolic cache here, since it contains both `mortars` and `u_threaded`. + # !!! NOTE: we reuse the hyperbolic cache here, since it contains both `mortars` and `u_threaded`. # !!! should we have a separate mortars/u_threaded in cache_parabolic? @trixi_timeit timer() "prolong2mortars" begin prolong2mortars!(cache, u_transformed, mesh, equations_parabolic, @@ -152,8 +152,8 @@ function calc_gradient!(gradients, u_transformed, t, end # Calculate mortar fluxes. These should reuse the hyperbolic version of `calc_mortar_flux`, - # along with a specialization on `calc_mortar_flux!` and `mortar_fluxes_to_elements!` for - # AbstractEquationsParabolic. + # along with a specialization on `calc_mortar_flux!` and `mortar_fluxes_to_elements!` for + # AbstractEquationsParabolic. @trixi_timeit timer() "mortar flux" begin calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, mesh, False(), # False() = no nonconservative terms @@ -265,7 +265,7 @@ function calc_gradient!(gradients, u_transformed, t, end # This version is called during `calc_gradients!` and must be specialized because the flux -# in the gradient is {u} which doesn't depend on normals. Thus, you don't need to scale by +# in the gradient is {u} which doesn't depend on normals. Thus, you don't need to scale by # 2 and flip the sign when storing the mortar fluxes back into surface_flux_values @inline function mortar_fluxes_to_elements!(surface_flux_values, mesh::P4estMesh{3}, @@ -369,7 +369,7 @@ end primary_j_node_index, interface_index) - flux_ = 0.5 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux + flux_ = 0.5f0 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux # Note that we don't flip the sign on the secondondary flux. This is because for parabolic terms, # the normals are not embedded in `flux_` for the parabolic gradient computations. @@ -558,7 +558,7 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, return nothing end -# This version is used for divergence flux computations +# This version is used for divergence flux computations function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{3}, equations_parabolic, dg::DG, cache_parabolic) @@ -603,7 +603,7 @@ function calc_interface_flux!(surface_flux_values, for j in eachnode(dg) for i in eachnode(dg) - # We prolong the viscous flux dotted with respect the outward normal on the + # We prolong the viscous flux dotted with respect the outward normal on the # primary element. We assume a BR-1 type of flux. viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, @@ -612,7 +612,7 @@ function calc_interface_flux!(surface_flux_values, j, interface) - flux = 0.5 * (viscous_flux_normal_ll + viscous_flux_normal_rr) + flux = 0.5f0 * (viscous_flux_normal_ll + viscous_flux_normal_rr) for v in eachvariable(equations_parabolic) surface_flux_values[v, i, j, primary_direction_index, primary_element] = flux[v] @@ -738,11 +738,11 @@ function prolong2mortars_divergence!(cache, flux_viscous, flux_viscous_z[v, i_large, j_large, k_large, element]) - # We prolong the viscous flux dotted with respect the outward normal - # on the small element. We scale by -1/2 here because the normal - # direction on the large element is negative 2x that of the small + # We prolong the viscous flux dotted with respect the outward normal + # on the small element. We scale by -1/2 here because the normal + # direction on the large element is negative 2x that of the small # element (these normal directions are "scaled" by the surface Jacobian) - u_buffer[v, i, j] = -0.5 * dot(flux_viscous, normal_direction) + u_buffer[v, i, j] = -0.5f0 * dot(flux_viscous, normal_direction) end i_large += i_large_step_i j_large += j_large_step_i @@ -779,8 +779,8 @@ function prolong2mortars_divergence!(cache, flux_viscous, return nothing end -# We specialize `calc_mortar_flux!` for the divergence part of -# the parabolic terms. +# We specialize `calc_mortar_flux!` for the divergence part of +# the parabolic terms. function calc_mortar_flux_divergence!(surface_flux_values, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations::AbstractEquationsParabolic, @@ -821,7 +821,7 @@ function calc_mortar_flux_divergence!(surface_flux_values, mortar] # TODO: parabolic; only BR1 at the moment - fstar[v, i, j, position] = 0.5 * (viscous_flux_normal_ll + + fstar[v, i, j, position] = 0.5f0 * (viscous_flux_normal_ll + viscous_flux_normal_rr) end @@ -849,7 +849,7 @@ function calc_mortar_flux_divergence!(surface_flux_values, end # NOTE: Use analogy to "calc_mortar_flux!" for hyperbolic eqs with no nonconservative terms. -# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for +# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for # hyperbolic terms with conserved terms only, i.e., no nonconservative terms. @inline function calc_mortar_flux!(fstar, mesh::P4estMesh{3}, @@ -865,7 +865,7 @@ end j_node_index, mortar_index) # TODO: parabolic; only BR1 at the moment - flux_ = 0.5 * (u_ll + u_rr) + flux_ = 0.5f0 * (u_ll + u_rr) # Copy flux to buffer set_node_vars!(fstar, flux_, equations, dg, i_node_index, j_node_index, position_index) diff --git a/src/solvers/dgsem_structured/containers_3d.jl b/src/solvers/dgsem_structured/containers_3d.jl index e843e869bf..75cc98bf2b 100644 --- a/src/solvers/dgsem_structured/containers_3d.jl +++ b/src/solvers/dgsem_structured/containers_3d.jl @@ -144,7 +144,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, for ii in eachnode(basis) # Multiply derivative_matrix to j-dimension to differentiate wrt η - result += 0.5 * derivative_matrix[j, ii] * + result += 0.5f0 * derivative_matrix[j, ii] * (node_coordinates[m, i, ii, k, element] * jacobian_matrix[l, 3, i, ii, k, element] - node_coordinates[l, i, ii, k, element] * @@ -160,7 +160,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, for ii in eachnode(basis) # Multiply derivative_matrix to k-dimension to differentiate wrt ζ - result += 0.5 * derivative_matrix[k, ii] * + result += 0.5f0 * derivative_matrix[k, ii] * (node_coordinates[m, i, j, ii, element] * jacobian_matrix[l, 2, i, j, ii, element] - node_coordinates[l, i, j, ii, element] * @@ -178,7 +178,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, for ii in eachnode(basis) # Multiply derivative_matrix to k-dimension to differentiate wrt ζ - result += 0.5 * derivative_matrix[k, ii] * + result += 0.5f0 * derivative_matrix[k, ii] * (node_coordinates[m, i, j, ii, element] * jacobian_matrix[l, 1, i, j, ii, element] - node_coordinates[l, i, j, ii, element] * @@ -194,7 +194,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, for ii in eachnode(basis) # Multiply derivative_matrix to i-dimension to differentiate wrt ξ - result += 0.5 * derivative_matrix[i, ii] * + result += 0.5f0 * derivative_matrix[i, ii] * (node_coordinates[m, ii, j, k, element] * jacobian_matrix[l, 3, ii, j, k, element] - node_coordinates[l, ii, j, k, element] * @@ -212,7 +212,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, for ii in eachnode(basis) # Multiply derivative_matrix to i-dimension to differentiate wrt ξ - result += 0.5 * derivative_matrix[i, ii] * + result += 0.5f0 * derivative_matrix[i, ii] * (node_coordinates[m, ii, j, k, element] * jacobian_matrix[l, 2, ii, j, k, element] - node_coordinates[l, ii, j, k, element] * @@ -228,7 +228,7 @@ function calc_contravariant_vectors!(contravariant_vectors::AbstractArray{<:Any, for ii in eachnode(basis) # Multiply derivative_matrix to j-dimension to differentiate wrt η - result += 0.5 * derivative_matrix[j, ii] * + result += 0.5f0 * derivative_matrix[j, ii] * (node_coordinates[m, i, ii, k, element] * jacobian_matrix[l, 1, i, ii, k, element] - node_coordinates[l, i, ii, k, element] * diff --git a/src/solvers/dgsem_structured/dg_2d.jl b/src/solvers/dgsem_structured/dg_2d.jl index 51b91eac6b..6023c9e619 100644 --- a/src/solvers/dgsem_structured/dg_2d.jl +++ b/src/solvers/dgsem_structured/dg_2d.jl @@ -50,8 +50,8 @@ function rhs!(du, u, t, end #= -`weak_form_kernel!` is only implemented for conserved terms as -non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, see `flux_differencing_kernel!`. This treatment is required to achieve, e.g., entropy-stability or well-balancedness. See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 @@ -128,7 +128,7 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, element) - Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii) + Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations) @@ -144,7 +144,7 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, element) - Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj) + Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations) @@ -193,7 +193,7 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, element) - Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii) + Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # Compute the contravariant nonconservative flux. fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_node, Ja1_avg, equations) @@ -207,7 +207,7 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, element) - Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj) + Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_node, Ja2_avg, @@ -217,7 +217,7 @@ end end # The factor 0.5 cancels the factor 2 in the flux differencing form - multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, + multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution, equations, dg, i, j, element) end end @@ -337,10 +337,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde1_L = ftilde1 + - 0.5 * nonconservative_flux(u_ll, u_rr, normal_direction, + 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations) ftilde1_R = ftilde1 + - 0.5 * nonconservative_flux(u_rr, u_ll, normal_direction, + 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, normal_direction, equations) set_node_vars!(fstar1_L, ftilde1_L, equations, dg, i, j) @@ -377,10 +377,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde2_L = ftilde2 + - 0.5 * nonconservative_flux(u_ll, u_rr, normal_direction, + 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations) ftilde2_R = ftilde2 + - 0.5 * nonconservative_flux(u_rr, u_ll, normal_direction, + 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, normal_direction, equations) set_node_vars!(fstar2_L, ftilde2_L, equations, dg, i, j) @@ -545,10 +545,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, i, right_direction, left_element] = flux[v] + - 0.5 * + 0.5f0 * noncons_left[v] surface_flux_values[v, i, left_direction, right_element] = flux[v] + - 0.5 * + 0.5f0 * noncons_right[v] end end diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index 4da402425e..e10a529545 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -48,7 +48,7 @@ # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, element) - Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii) + Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant sharp flux in the direction of the averaged contravariant vector fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations) @@ -85,7 +85,7 @@ # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, element) - Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj) + Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant sharp flux in the direction of the averaged contravariant vector fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations) multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2, diff --git a/src/solvers/dgsem_structured/dg_3d.jl b/src/solvers/dgsem_structured/dg_3d.jl index 1df9f40889..e8430eaa49 100644 --- a/src/solvers/dgsem_structured/dg_3d.jl +++ b/src/solvers/dgsem_structured/dg_3d.jl @@ -50,8 +50,8 @@ function rhs!(du, u, t, end #= -`weak_form_kernel!` is only implemented for conserved terms as -non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, see `flux_differencing_kernel!`. This treatment is required to achieve, e.g., entropy-stability or well-balancedness. See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 @@ -145,7 +145,7 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, k, element) - Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii) + Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations) @@ -161,7 +161,7 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, k, element) - Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj) + Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations) @@ -177,7 +177,7 @@ end # pull the contravariant vectors and compute the average Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors, i, j, kk, element) - Ja3_avg = 0.5 * (Ja3_node + Ja3_node_kk) + Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk) # compute the contravariant sharp flux in the direction of the # averaged contravariant vector fluxtilde3 = volume_flux(u_node, u_node_kk, Ja3_avg, equations) @@ -225,7 +225,7 @@ end # pull the contravariant vectors and compute the average Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, k, element) - Ja1_avg = 0.5 * (Ja1_node + Ja1_node_ii) + Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector fluxtilde1 = nonconservative_flux(u_node, u_node_ii, Ja1_node, Ja1_avg, @@ -240,7 +240,7 @@ end # pull the contravariant vectors and compute the average Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, k, element) - Ja2_avg = 0.5 * (Ja2_node + Ja2_node_jj) + Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector fluxtilde2 = nonconservative_flux(u_node, u_node_jj, Ja2_node, Ja2_avg, @@ -255,7 +255,7 @@ end # pull the contravariant vectors and compute the average Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors, i, j, kk, element) - Ja3_avg = 0.5 * (Ja3_node + Ja3_node_kk) + Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk) # compute the contravariant nonconservative flux in the direction of the # averaged contravariant vector fluxtilde3 = nonconservative_flux(u_node, u_node_kk, Ja3_node, Ja3_avg, @@ -265,7 +265,7 @@ end end # The factor 0.5 cancels the factor 2 in the flux differencing form - multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, + multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution, equations, dg, i, j, k, element) end end @@ -411,10 +411,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde_L = ftilde + - 0.5 * nonconservative_flux(u_ll, u_rr, normal_direction, + 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations) ftilde_R = ftilde + - 0.5 * nonconservative_flux(u_rr, u_ll, normal_direction, + 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, normal_direction, equations) set_node_vars!(fstar1_L, ftilde_L, equations, dg, i, j, k) @@ -449,10 +449,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde_L = ftilde + - 0.5 * nonconservative_flux(u_ll, u_rr, normal_direction, + 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations) ftilde_R = ftilde + - 0.5 * nonconservative_flux(u_rr, u_ll, normal_direction, + 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, normal_direction, equations) set_node_vars!(fstar2_L, ftilde_L, equations, dg, i, j, k) @@ -487,10 +487,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs ftilde_L = ftilde + - 0.5 * nonconservative_flux(u_ll, u_rr, normal_direction, + 0.5f0 * nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations) ftilde_R = ftilde + - 0.5 * nonconservative_flux(u_rr, u_ll, normal_direction, + 0.5f0 * nonconservative_flux(u_rr, u_ll, normal_direction, normal_direction, equations) set_node_vars!(fstar3_L, ftilde_L, equations, dg, i, j, k) @@ -681,10 +681,10 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, i, j, right_direction, left_element] = flux[v] + - 0.5 * + 0.5f0 * noncons_left[v] surface_flux_values[v, i, j, left_direction, right_element] = flux[v] + - 0.5 * + 0.5f0 * noncons_right[v] end end diff --git a/src/solvers/dgsem_structured/indicators_1d.jl b/src/solvers/dgsem_structured/indicators_1d.jl index 4299ec603d..a6d518699d 100644 --- a/src/solvers/dgsem_structured/indicators_1d.jl +++ b/src/solvers/dgsem_structured/indicators_1d.jl @@ -19,8 +19,9 @@ function apply_smoothing!(mesh::StructuredMesh{1}, alpha, alpha_tmp, dg, cache) left = cache.elements.left_neighbors[1, element] # Apply smoothing - alpha[left] = max(alpha_tmp[left], 0.5 * alpha_tmp[element], alpha[left]) - alpha[element] = max(alpha_tmp[element], 0.5 * alpha_tmp[left], alpha[element]) + alpha[left] = max(alpha_tmp[left], 0.5f0 * alpha_tmp[element], alpha[left]) + alpha[element] = max(alpha_tmp[element], 0.5f0 * alpha_tmp[left], + alpha[element]) end end end # @muladd diff --git a/src/solvers/dgsem_structured/indicators_2d.jl b/src/solvers/dgsem_structured/indicators_2d.jl index f4b07b70cb..52d6ac2a95 100644 --- a/src/solvers/dgsem_structured/indicators_2d.jl +++ b/src/solvers/dgsem_structured/indicators_2d.jl @@ -20,11 +20,13 @@ function apply_smoothing!(mesh::StructuredMesh{2}, alpha, alpha_tmp, dg, cache) lower = cache.elements.left_neighbors[2, element] # Apply smoothing - alpha[left] = max(alpha_tmp[left], 0.5 * alpha_tmp[element], alpha[left]) - alpha[element] = max(alpha_tmp[element], 0.5 * alpha_tmp[left], alpha[element]) + alpha[left] = max(alpha_tmp[left], 0.5f0 * alpha_tmp[element], alpha[left]) + alpha[element] = max(alpha_tmp[element], 0.5f0 * alpha_tmp[left], + alpha[element]) - alpha[lower] = max(alpha_tmp[lower], 0.5 * alpha_tmp[element], alpha[lower]) - alpha[element] = max(alpha_tmp[element], 0.5 * alpha_tmp[lower], alpha[element]) + alpha[lower] = max(alpha_tmp[lower], 0.5f0 * alpha_tmp[element], alpha[lower]) + alpha[element] = max(alpha_tmp[element], 0.5f0 * alpha_tmp[lower], + alpha[element]) end end end # @muladd diff --git a/src/solvers/dgsem_structured/indicators_3d.jl b/src/solvers/dgsem_structured/indicators_3d.jl index 155bf50dc6..8b477da2e8 100644 --- a/src/solvers/dgsem_structured/indicators_3d.jl +++ b/src/solvers/dgsem_structured/indicators_3d.jl @@ -21,14 +21,17 @@ function apply_smoothing!(mesh::StructuredMesh{3}, alpha, alpha_tmp, dg, cache) front = cache.elements.left_neighbors[3, element] # Apply smoothing - alpha[left] = max(alpha_tmp[left], 0.5 * alpha_tmp[element], alpha[left]) - alpha[element] = max(alpha_tmp[element], 0.5 * alpha_tmp[left], alpha[element]) + alpha[left] = max(alpha_tmp[left], 0.5f0 * alpha_tmp[element], alpha[left]) + alpha[element] = max(alpha_tmp[element], 0.5f0 * alpha_tmp[left], + alpha[element]) - alpha[lower] = max(alpha_tmp[lower], 0.5 * alpha_tmp[element], alpha[lower]) - alpha[element] = max(alpha_tmp[element], 0.5 * alpha_tmp[lower], alpha[element]) + alpha[lower] = max(alpha_tmp[lower], 0.5f0 * alpha_tmp[element], alpha[lower]) + alpha[element] = max(alpha_tmp[element], 0.5f0 * alpha_tmp[lower], + alpha[element]) - alpha[front] = max(alpha_tmp[front], 0.5 * alpha_tmp[element], alpha[front]) - alpha[element] = max(alpha_tmp[element], 0.5 * alpha_tmp[front], alpha[element]) + alpha[front] = max(alpha_tmp[front], 0.5f0 * alpha_tmp[element], alpha[front]) + alpha[element] = max(alpha_tmp[element], 0.5f0 * alpha_tmp[front], + alpha[element]) end end end # @muladd diff --git a/src/solvers/dgsem_t8code/containers_3d.jl b/src/solvers/dgsem_t8code/containers_3d.jl index 4d56bc734a..1375782631 100644 --- a/src/solvers/dgsem_t8code/containers_3d.jl +++ b/src/solvers/dgsem_t8code/containers_3d.jl @@ -44,13 +44,13 @@ function calc_node_coordinates!(node_coordinates, pointer(element_coords)) nodes_out_x = (2 * - (element_length * 0.5 * (nodes .+ 1) .+ element_coords[1]) .- + (element_length * 0.5f0 * (nodes .+ 1) .+ element_coords[1]) .- 1) nodes_out_y = (2 * - (element_length * 0.5 * (nodes .+ 1) .+ element_coords[2]) .- + (element_length * 0.5f0 * (nodes .+ 1) .+ element_coords[2]) .- 1) nodes_out_z = (2 * - (element_length * 0.5 * (nodes .+ 1) .+ element_coords[3]) .- + (element_length * 0.5f0 * (nodes .+ 1) .+ element_coords[3]) .- 1) polynomial_interpolation_matrix!(matrix1, mesh.nodes, nodes_out_x, diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 225b85a059..641c944d5f 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -145,8 +145,8 @@ function calc_volume_integral!(du, u, end #= -`weak_form_kernel!` is only implemented for conserved terms as -non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, see `flux_differencing_kernel!`. This treatment is required to achieve, e.g., entropy-stability or well-balancedness. See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 @@ -245,7 +245,7 @@ end end # The factor 0.5 cancels the factor 2 in the flux differencing form - multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, + multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution, equations, dg, i, element) end end @@ -377,8 +377,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - f1_L = f1 + 0.5 * nonconservative_flux(u_ll, u_rr, 1, equations) - f1_R = f1 + 0.5 * nonconservative_flux(u_rr, u_ll, 1, equations) + f1_L = f1 + 0.5f0 * nonconservative_flux(u_ll, u_rr, 1, equations) + f1_R = f1 + 0.5f0 * nonconservative_flux(u_rr, u_ll, 1, equations) # Copy to temporary storage set_node_vars!(fstar1_L, f1_L, equations, dg, i) @@ -471,9 +471,9 @@ function calc_interface_flux!(surface_flux_values, # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, left_direction, left_id] = flux[v] + - 0.5 * noncons_left[v] + 0.5f0 * noncons_left[v] surface_flux_values[v, right_direction, right_id] = flux[v] + - 0.5 * noncons_right[v] + 0.5f0 * noncons_right[v] end end @@ -595,7 +595,7 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A # Copy flux to left and right element storage for v in eachvariable(equations) surface_flux_values[v, direction, neighbor] = flux[v] + - 0.5 * noncons_flux[v] + 0.5f0 * noncons_flux[v] end end diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index 0017f9ca88..fb96e0535a 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -190,7 +190,7 @@ function calc_interface_flux!(surface_flux_values, # Compute interface flux as mean of left and right viscous fluxes # TODO: parabolic; only BR1 at the moment - flux = 0.5 * (flux_ll + flux_rr) + flux = 0.5f0 * (flux_ll + flux_rr) # Copy flux to left and right element storage for v in eachvariable(equations_parabolic) @@ -456,7 +456,7 @@ function calc_gradient!(gradients, u_transformed, t, # Call pointwise Riemann solver u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, interface) - flux = 0.5 * (u_ll + u_rr) + flux = 0.5f0 * (u_ll + u_rr) # Copy flux to left and right element storage for v in eachvariable(equations_parabolic) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 4f2c72cd67..734a31d547 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -196,8 +196,8 @@ function calc_volume_integral!(du, u, end #= -`weak_form_kernel!` is only implemented for conserved terms as -non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, see `flux_differencing_kernel!`. This treatment is required to achieve, e.g., entropy-stability or well-balancedness. See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 @@ -325,7 +325,7 @@ end end # The factor 0.5 cancels the factor 2 in the flux differencing form - multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, + multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution, equations, dg, i, j, element) end end @@ -504,8 +504,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - f1_L = f1 + 0.5 * nonconservative_flux(u_ll, u_rr, 1, equations) - f1_R = f1 + 0.5 * nonconservative_flux(u_rr, u_ll, 1, equations) + f1_L = f1 + 0.5f0 * nonconservative_flux(u_ll, u_rr, 1, equations) + f1_R = f1 + 0.5f0 * nonconservative_flux(u_rr, u_ll, 1, equations) # Copy to temporary storage set_node_vars!(fstar1_L, f1_L, equations, dg, i, j) @@ -530,8 +530,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - f2_L = f2 + 0.5 * nonconservative_flux(u_ll, u_rr, 2, equations) - f2_R = f2 + 0.5 * nonconservative_flux(u_rr, u_ll, 2, equations) + f2_L = f2 + 0.5f0 * nonconservative_flux(u_ll, u_rr, 2, equations) + f2_R = f2 + 0.5f0 * nonconservative_flux(u_rr, u_ll, 2, equations) # Copy to temporary storage set_node_vars!(fstar2_L, f2_L, equations, dg, i, j) @@ -637,10 +637,10 @@ function calc_interface_flux!(surface_flux_values, # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, i, left_direction, left_id] = flux[v] + - 0.5 * + 0.5f0 * noncons_left[v] surface_flux_values[v, i, right_direction, right_id] = flux[v] + - 0.5 * + 0.5f0 * noncons_right[v] end end @@ -790,7 +790,7 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A # Copy flux to left and right element storage for v in eachvariable(equations) surface_flux_values[v, i, direction, neighbor] = flux[v] + - 0.5 * noncons_flux[v] + 0.5f0 * noncons_flux[v] end end end @@ -968,9 +968,9 @@ function calc_mortar_flux!(surface_flux_values, noncons_lower = nonconservative_flux(u_lower_ll, u_lower_rr, orientation, equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper, 0.5, noncons_upper, equations, + multiply_add_to_node_vars!(fstar_upper, 0.5f0, noncons_upper, equations, dg, i) - multiply_add_to_node_vars!(fstar_lower, 0.5, noncons_lower, equations, + multiply_add_to_node_vars!(fstar_lower, 0.5f0, noncons_lower, equations, dg, i) end else # large_sides[mortar] == 2 -> small elements on the left @@ -986,9 +986,9 @@ function calc_mortar_flux!(surface_flux_values, noncons_lower = nonconservative_flux(u_lower_rr, u_lower_ll, orientation, equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper, 0.5, noncons_upper, equations, + multiply_add_to_node_vars!(fstar_upper, 0.5f0, noncons_upper, equations, dg, i) - multiply_add_to_node_vars!(fstar_lower, 0.5, noncons_lower, equations, + multiply_add_to_node_vars!(fstar_lower, 0.5f0, noncons_lower, equations, dg, i) end end diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index a6c962e03c..f84a08cb69 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -228,7 +228,7 @@ function calc_interface_flux!(surface_flux_values, # Compute interface flux as mean of left and right viscous fluxes # TODO: parabolic; only BR1 at the moment - flux = 0.5 * (flux_ll + flux_rr) + flux = 0.5f0 * (flux_ll + flux_rr) # Copy flux to left and right element storage for v in eachvariable(equations_parabolic) @@ -513,8 +513,8 @@ end # `cache` is the hyperbolic cache, i.e., in particular not `cache_parabolic`. # This is because mortar handling is done in the (hyperbolic) `cache`. -# Specialization `flux_viscous::Vector{Array{uEltype, 4}}` needed since -#`prolong2mortars!` in dg_2d.jl is used for both purely hyperbolic and +# Specialization `flux_viscous::Vector{Array{uEltype, 4}}` needed since +#`prolong2mortars!` in dg_2d.jl is used for both purely hyperbolic and # hyperbolic-parabolic systems. function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 4}}, mesh::TreeMesh{2}, @@ -654,7 +654,7 @@ end u_ll, u_rr = get_surface_node_vars(u_interfaces, equations_parabolic, dg, i, interface) # TODO: parabolic; only BR1 at the moment - flux = 0.5 * (u_ll + u_rr) + flux = 0.5f0 * (u_ll + u_rr) # Copy flux to left and right element storage set_node_vars!(destination, flux, equations_parabolic, dg, i) @@ -801,7 +801,7 @@ function calc_gradient!(gradients, u_transformed, t, u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, i, interface) - flux = 0.5 * (u_ll + u_rr) + flux = 0.5f0 * (u_ll + u_rr) # Copy flux to left and right element storage for v in eachvariable(equations_parabolic) diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 2cc5e2cae8..998790ef8e 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -262,11 +262,11 @@ end flux1_noncons = volume_flux_noncons(u_node, u_node_ii, 1, equations, NonConservativeSymmetric(), noncons) multiply_add_to_node_vars!(flux_noncons_temp, - 0.5 * derivative_split[i, ii], + 0.5f0 * derivative_split[i, ii], flux1_noncons, equations, dg, noncons, i, j) multiply_add_to_node_vars!(flux_noncons_temp, - 0.5 * derivative_split[ii, i], + 0.5f0 * derivative_split[ii, i], flux1_noncons, equations, dg, noncons, ii, j) end diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 2e7e882f5e..acd90cee09 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -227,8 +227,8 @@ function calc_volume_integral!(du, u, end #= -`weak_form_kernel!` is only implemented for conserved terms as -non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, +`weak_form_kernel!` is only implemented for conserved terms as +non-conservative terms should always be discretized in conjunction with a flux-splitting scheme, see `flux_differencing_kernel!`. This treatment is required to achieve, e.g., entropy-stability or well-balancedness. See also https://github.com/trixi-framework/Trixi.jl/issues/1671#issuecomment-1765644064 @@ -374,7 +374,7 @@ end end # The factor 0.5 cancels the factor 2 in the flux differencing form - multiply_add_to_node_vars!(du, alpha * 0.5, integral_contribution, equations, + multiply_add_to_node_vars!(du, alpha * 0.5f0, integral_contribution, equations, dg, i, j, k, element) end end @@ -550,8 +550,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - flux_L = flux + 0.5 * nonconservative_flux(u_ll, u_rr, 1, equations) - flux_R = flux + 0.5 * nonconservative_flux(u_rr, u_ll, 1, equations) + flux_L = flux + 0.5f0 * nonconservative_flux(u_ll, u_rr, 1, equations) + flux_R = flux + 0.5f0 * nonconservative_flux(u_rr, u_ll, 1, equations) set_node_vars!(fstar1_L, flux_L, equations, dg, i, j, k) set_node_vars!(fstar1_R, flux_R, equations, dg, i, j, k) @@ -573,8 +573,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - flux_L = flux + 0.5 * nonconservative_flux(u_ll, u_rr, 2, equations) - flux_R = flux + 0.5 * nonconservative_flux(u_rr, u_ll, 2, equations) + flux_L = flux + 0.5f0 * nonconservative_flux(u_ll, u_rr, 2, equations) + flux_R = flux + 0.5f0 * nonconservative_flux(u_rr, u_ll, 2, equations) set_node_vars!(fstar2_L, flux_L, equations, dg, i, j, k) set_node_vars!(fstar2_R, flux_R, equations, dg, i, j, k) @@ -596,8 +596,8 @@ end # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - flux_L = flux + 0.5 * nonconservative_flux(u_ll, u_rr, 3, equations) - flux_R = flux + 0.5 * nonconservative_flux(u_rr, u_ll, 3, equations) + flux_L = flux + 0.5f0 * nonconservative_flux(u_ll, u_rr, 3, equations) + flux_R = flux + 0.5f0 * nonconservative_flux(u_rr, u_ll, 3, equations) set_node_vars!(fstar3_L, flux_L, equations, dg, i, j, k) set_node_vars!(fstar3_R, flux_R, equations, dg, i, j, k) @@ -712,10 +712,10 @@ function calc_interface_flux!(surface_flux_values, # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, i, j, left_direction, left_id] = flux[v] + - 0.5 * + 0.5f0 * noncons_left[v] surface_flux_values[v, i, j, right_direction, right_id] = flux[v] + - 0.5 * + 0.5f0 * noncons_right[v] end end @@ -1145,13 +1145,15 @@ function calc_mortar_flux!(surface_flux_values, u_lower_right_rr, orientation, equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper_left, 0.5, noncons_upper_left, + multiply_add_to_node_vars!(fstar_upper_left, 0.5f0, noncons_upper_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_upper_right, 0.5, noncons_upper_right, + multiply_add_to_node_vars!(fstar_upper_right, 0.5f0, + noncons_upper_right, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_left, 0.5, noncons_lower_left, + multiply_add_to_node_vars!(fstar_lower_left, 0.5f0, noncons_lower_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_right, 0.5, noncons_lower_right, + multiply_add_to_node_vars!(fstar_lower_right, 0.5f0, + noncons_lower_right, equations, dg, i, j) end else # large_sides[mortar] == 2 -> small elements on the left @@ -1185,13 +1187,15 @@ function calc_mortar_flux!(surface_flux_values, u_lower_right_ll, orientation, equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper_left, 0.5, noncons_upper_left, + multiply_add_to_node_vars!(fstar_upper_left, 0.5f0, noncons_upper_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_upper_right, 0.5, noncons_upper_right, + multiply_add_to_node_vars!(fstar_upper_right, 0.5f0, + noncons_upper_right, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_left, 0.5, noncons_lower_left, + multiply_add_to_node_vars!(fstar_lower_left, 0.5f0, noncons_lower_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_right, 0.5, noncons_lower_right, + multiply_add_to_node_vars!(fstar_lower_right, 0.5f0, + noncons_lower_right, equations, dg, i, j) end end diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index d550474474..69956d5834 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -141,7 +141,7 @@ function calc_interface_flux!(surface_flux_values, # Compute interface flux as mean of left and right viscous fluxes # TODO: parabolic; only BR1 at the moment - flux = 0.5 * (flux_ll + flux_rr) + flux = 0.5f0 * (flux_ll + flux_rr) # Copy flux to left and right element storage for v in eachvariable(equations_parabolic) @@ -488,8 +488,8 @@ end # `cache` is the hyperbolic cache, i.e., in particular not `cache_parabolic`. # This is because mortar handling is done in the (hyperbolic) `cache`. -# Specialization `flux_viscous::Vector{Array{uEltype, 4}}` needed since -#`prolong2mortars!` in dg_2d.jl is used for both purely hyperbolic and +# Specialization `flux_viscous::Vector{Array{uEltype, 4}}` needed since +#`prolong2mortars!` in dg_2d.jl is used for both purely hyperbolic and # hyperbolic-parabolic systems. function prolong2mortars!(cache, flux_viscous::Vector{Array{uEltype, 5}}, @@ -773,7 +773,7 @@ end u_ll, u_rr = get_surface_node_vars(u_interfaces, equations_parabolic, dg, i, j, interface) # TODO: parabolic; only BR1 at the moment - flux = 0.5 * (u_ll + u_rr) + flux = 0.5f0 * (u_ll + u_rr) # Copy flux to left and right element storage set_node_vars!(destination, flux, equations_parabolic, dg, i, j) @@ -854,7 +854,7 @@ function calc_gradient!(gradients, u_transformed, t, u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u, equations_parabolic, dg, i, j, interface) - flux = 0.5 * (u_ll + u_rr) + flux = 0.5f0 * (u_ll + u_rr) # Copy flux to left and right element storage for v in eachvariable(equations_parabolic) diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index 9f25a6d2db..323c1236c2 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -114,8 +114,11 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh, equations, dg::DGSEM end # magic parameters - threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25) - parameter_s = log((1 - 0.0001) / 0.0001) + # TODO: Are there better values for Float32? + RealT = real(dg) + threshold = 0.5f0 * 10^(convert(RealT, -1.8) * nnodes(dg)^convert(RealT, 0.25)) + o_0001 = convert(RealT, 0.0001) + parameter_s = log((1 - o_0001) / o_0001) @threaded for element in eachelement(dg, cache) # This is dispatched by mesh dimension. diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 4796ddcc60..9d11d05edc 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -104,8 +104,8 @@ function apply_smoothing!(mesh::Union{TreeMesh{1}, P4estMesh{1}}, alpha, alpha_t right = cache.interfaces.neighbor_ids[2, interface] # Apply smoothing - alpha[left] = max(alpha_tmp[left], 0.5 * alpha_tmp[right], alpha[left]) - alpha[right] = max(alpha_tmp[right], 0.5 * alpha_tmp[left], alpha[right]) + alpha[left] = max(alpha_tmp[left], 0.5f0 * alpha_tmp[right], alpha[left]) + alpha[right] = max(alpha_tmp[right], 0.5f0 * alpha_tmp[left], alpha[right]) end end diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 665d2254e5..6ef65c4e9d 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -111,8 +111,8 @@ function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, right = cache.interfaces.neighbor_ids[2, interface] # Apply smoothing - alpha[left] = max(alpha_tmp[left], 0.5 * alpha_tmp[right], alpha[left]) - alpha[right] = max(alpha_tmp[right], 0.5 * alpha_tmp[left], alpha[right]) + alpha[left] = max(alpha_tmp[left], 0.5f0 * alpha_tmp[right], alpha[left]) + alpha[right] = max(alpha_tmp[right], 0.5f0 * alpha_tmp[left], alpha[right]) end # Loop over L2 mortars @@ -123,10 +123,10 @@ function apply_smoothing!(mesh::Union{TreeMesh{2}, P4estMesh{2}, T8codeMesh{2}}, large = cache.mortars.neighbor_ids[3, mortar] # Apply smoothing - alpha[lower] = max(alpha_tmp[lower], 0.5 * alpha_tmp[large], alpha[lower]) - alpha[upper] = max(alpha_tmp[upper], 0.5 * alpha_tmp[large], alpha[upper]) - alpha[large] = max(alpha_tmp[large], 0.5 * alpha_tmp[lower], alpha[large]) - alpha[large] = max(alpha_tmp[large], 0.5 * alpha_tmp[upper], alpha[large]) + alpha[lower] = max(alpha_tmp[lower], 0.5f0 * alpha_tmp[large], alpha[lower]) + alpha[upper] = max(alpha_tmp[upper], 0.5f0 * alpha_tmp[large], alpha[upper]) + alpha[large] = max(alpha_tmp[large], 0.5f0 * alpha_tmp[lower], alpha[large]) + alpha[large] = max(alpha_tmp[large], 0.5f0 * alpha_tmp[upper], alpha[large]) end return alpha diff --git a/src/solvers/dgsem_tree/indicators_3d.jl b/src/solvers/dgsem_tree/indicators_3d.jl index a11a8e06e4..6d4ee63618 100644 --- a/src/solvers/dgsem_tree/indicators_3d.jl +++ b/src/solvers/dgsem_tree/indicators_3d.jl @@ -116,8 +116,8 @@ function apply_smoothing!(mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, right = cache.interfaces.neighbor_ids[2, interface] # Apply smoothing - alpha[left] = max(alpha_tmp[left], 0.5 * alpha_tmp[right], alpha[left]) - alpha[right] = max(alpha_tmp[right], 0.5 * alpha_tmp[left], alpha[right]) + alpha[left] = max(alpha_tmp[left], 0.5f0 * alpha_tmp[right], alpha[left]) + alpha[right] = max(alpha_tmp[right], 0.5f0 * alpha_tmp[left], alpha[right]) end # Loop over L2 mortars @@ -130,19 +130,23 @@ function apply_smoothing!(mesh::Union{TreeMesh{3}, P4estMesh{3}, T8codeMesh{3}}, large = cache.mortars.neighbor_ids[5, mortar] # Apply smoothing - alpha[lower_left] = max(alpha_tmp[lower_left], 0.5 * alpha_tmp[large], + alpha[lower_left] = max(alpha_tmp[lower_left], 0.5f0 * alpha_tmp[large], alpha[lower_left]) - alpha[lower_right] = max(alpha_tmp[lower_right], 0.5 * alpha_tmp[large], + alpha[lower_right] = max(alpha_tmp[lower_right], 0.5f0 * alpha_tmp[large], alpha[lower_right]) - alpha[upper_left] = max(alpha_tmp[upper_left], 0.5 * alpha_tmp[large], + alpha[upper_left] = max(alpha_tmp[upper_left], 0.5f0 * alpha_tmp[large], alpha[upper_left]) - alpha[upper_right] = max(alpha_tmp[upper_right], 0.5 * alpha_tmp[large], + alpha[upper_right] = max(alpha_tmp[upper_right], 0.5f0 * alpha_tmp[large], alpha[upper_right]) - alpha[large] = max(alpha_tmp[large], 0.5 * alpha_tmp[lower_left], alpha[large]) - alpha[large] = max(alpha_tmp[large], 0.5 * alpha_tmp[lower_right], alpha[large]) - alpha[large] = max(alpha_tmp[large], 0.5 * alpha_tmp[upper_left], alpha[large]) - alpha[large] = max(alpha_tmp[large], 0.5 * alpha_tmp[upper_right], alpha[large]) + alpha[large] = max(alpha_tmp[large], 0.5f0 * alpha_tmp[lower_left], + alpha[large]) + alpha[large] = max(alpha_tmp[large], 0.5f0 * alpha_tmp[lower_right], + alpha[large]) + alpha[large] = max(alpha_tmp[large], 0.5f0 * alpha_tmp[upper_left], + alpha[large]) + alpha[large] = max(alpha_tmp[large], 0.5f0 * alpha_tmp[upper_right], + alpha[large]) end end diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl index 4095f0853f..b185cd4c54 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl @@ -588,7 +588,7 @@ end # Check bounds if (beta < beta_L) || (beta > beta_R) || (dgoal_dbeta == 0) || isnan(beta) # Out of bounds, do a bisection step - beta = 0.5 * (beta_L + beta_R) + beta = 0.5f0 * (beta_L + beta_R) # Get new u u_curr = u + beta * dt * antidiffusive_flux diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index ce602e178d..645bf05f07 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -20,7 +20,8 @@ function create_cache(mesh::UnstructuredMesh2D, equations, # perform a check on the sufficient metric identities condition for free-stream preservation # and halt computation if it fails - if !isapprox(max_discrete_metric_identities(dg, cache), 0, atol = 1e-12) + atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0)) + if !isapprox(max_discrete_metric_identities(dg, cache), 0, atol = atol) error("metric terms fail free-stream preservation check with maximum error $(max_discrete_metric_identities(dg, cache))") end @@ -260,10 +261,10 @@ function calc_interface_flux!(surface_flux_values, # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, primary_index, primary_side, primary_element] = (flux[v] + - 0.5 * + 0.5f0 * noncons_primary[v]) surface_flux_values[v, secondary_index, secondary_side, secondary_element] = -(flux[v] + - 0.5 * + 0.5f0 * noncons_secondary[v]) end @@ -459,7 +460,7 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, node_index, side_index, element_index] = flux[v] + - 0.5 * + 0.5f0 * noncons_flux[v] end end diff --git a/src/solvers/dgsem_unstructured/indicators_2d.jl b/src/solvers/dgsem_unstructured/indicators_2d.jl index 8052534ad4..e331cb5ee7 100644 --- a/src/solvers/dgsem_unstructured/indicators_2d.jl +++ b/src/solvers/dgsem_unstructured/indicators_2d.jl @@ -17,8 +17,8 @@ function apply_smoothing!(mesh::UnstructuredMesh2D, alpha, alpha_tmp, dg, cache) right = cache.interfaces.element_ids[2, interface] # Apply smoothing - alpha[left] = max(alpha_tmp[left], 0.5 * alpha_tmp[right], alpha[left]) - alpha[right] = max(alpha_tmp[right], 0.5 * alpha_tmp[left], alpha[right]) + alpha[left] = max(alpha_tmp[left], 0.5f0 * alpha_tmp[right], alpha[left]) + alpha[right] = max(alpha_tmp[right], 0.5f0 * alpha_tmp[left], alpha[right]) end end end # @muladd diff --git a/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl b/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl index 75b9a1f4da..5949e1f97e 100644 --- a/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl +++ b/src/solvers/dgsem_unstructured/mappings_geometry_curved_2d.jl @@ -11,10 +11,10 @@ function transfinite_quad_map(xi, eta, surface_curves::AbstractVector{<:CurvedSurface}) # evaluate the gamma curves to get the four corner points of the element - x_corner1, y_corner1 = evaluate_at(-1.0, surface_curves[1]) - x_corner2, y_corner2 = evaluate_at(1.0, surface_curves[1]) - x_corner3, y_corner3 = evaluate_at(1.0, surface_curves[3]) - x_corner4, y_corner4 = evaluate_at(-1.0, surface_curves[3]) + x_corner1, y_corner1 = evaluate_at(-1, surface_curves[1]) + x_corner2, y_corner2 = evaluate_at(1, surface_curves[1]) + x_corner3, y_corner3 = evaluate_at(1, surface_curves[3]) + x_corner4, y_corner4 = evaluate_at(-1, surface_curves[3]) # evaluate along the gamma curves at a particular point (ξ, η) in computational space to get # the value (x,y) in physical space @@ -23,17 +23,15 @@ function transfinite_quad_map(xi, eta, surface_curves::AbstractVector{<:CurvedSu x3, y3 = evaluate_at(xi, surface_curves[3]) x4, y4 = evaluate_at(eta, surface_curves[4]) - x = (0.5 * - ((1.0 - xi) * x4 + (1.0 + xi) * x2 + (1.0 - eta) * x1 + (1.0 + eta) * x3) + x = (0.5f0 * ((1 - xi) * x4 + (1 + xi) * x2 + (1 - eta) * x1 + (1 + eta) * x3) - - 0.25 * ((1.0 - xi) * ((1.0 - eta) * x_corner1 + (1.0 + eta) * x_corner4) + - (1.0 + xi) * ((1.0 - eta) * x_corner2 + (1.0 + eta) * x_corner3))) + 0.25f0 * ((1 - xi) * ((1 - eta) * x_corner1 + (1 + eta) * x_corner4) + + (1 + xi) * ((1 - eta) * x_corner2 + (1 + eta) * x_corner3))) - y = (0.5 * - ((1.0 - xi) * y4 + (1.0 + xi) * y2 + (1.0 - eta) * y1 + (1.0 + eta) * y3) + y = (0.5f0 * ((1 - xi) * y4 + (1 + xi) * y2 + (1 - eta) * y1 + (1 + eta) * y3) - - 0.25 * ((1.0 - xi) * ((1.0 - eta) * y_corner1 + (1.0 + eta) * y_corner4) + - (1.0 + xi) * ((1.0 - eta) * y_corner2 + (1.0 + eta) * y_corner3))) + 0.25f0 * ((1 - xi) * ((1 - eta) * y_corner1 + (1 + eta) * y_corner4) + + (1 + xi) * ((1 - eta) * y_corner2 + (1 + eta) * y_corner3))) return x, y end @@ -44,10 +42,10 @@ function transfinite_quad_map_metrics(xi, eta, surface_curves::AbstractVector{<:CurvedSurface}) # evaluate the gamma curves to get the four corner points of the element - x_corner1, y_corner1 = evaluate_at(-1.0, surface_curves[1]) - x_corner2, y_corner2 = evaluate_at(1.0, surface_curves[1]) - x_corner3, y_corner3 = evaluate_at(1.0, surface_curves[3]) - x_corner4, y_corner4 = evaluate_at(-1.0, surface_curves[3]) + x_corner1, y_corner1 = evaluate_at(-1, surface_curves[1]) + x_corner2, y_corner2 = evaluate_at(1, surface_curves[1]) + x_corner3, y_corner3 = evaluate_at(1, surface_curves[3]) + x_corner4, y_corner4 = evaluate_at(-1, surface_curves[3]) # evaluate along the gamma curves at a particular point (ξ, η) in computational space to get # the value (x,y) in physical space @@ -63,25 +61,25 @@ function transfinite_quad_map_metrics(xi, eta, x3_prime, y3_prime = derivative_at(xi, surface_curves[3]) x4_prime, y4_prime = derivative_at(eta, surface_curves[4]) - X_xi = (0.5 * (x2 - x4 + (1.0 - eta) * x1_prime + (1.0 + eta) * x3_prime) + X_xi = (0.5f0 * (x2 - x4 + (1 - eta) * x1_prime + (1 + eta) * x3_prime) - - 0.25 * ((1.0 - eta) * (x_corner2 - x_corner1) + - (1.0 + eta) * (x_corner3 - x_corner4))) + 0.25f0 * ((1 - eta) * (x_corner2 - x_corner1) + + (1 + eta) * (x_corner3 - x_corner4))) - X_eta = (0.5 * ((1.0 - xi) * x4_prime + (1.0 + xi) * x2_prime + x3 - x1) + X_eta = (0.5f0 * ((1 - xi) * x4_prime + (1 + xi) * x2_prime + x3 - x1) - - 0.25 * ((1.0 - xi) * (x_corner4 - x_corner1) + - (1.0 + xi) * (x_corner3 - x_corner2))) + 0.25f0 * ((1 - xi) * (x_corner4 - x_corner1) + + (1 + xi) * (x_corner3 - x_corner2))) - Y_xi = (0.5 * (y2 - y4 + (1.0 - eta) * y1_prime + (1.0 + eta) * y3_prime) + Y_xi = (0.5f0 * (y2 - y4 + (1 - eta) * y1_prime + (1 + eta) * y3_prime) - - 0.25 * ((1.0 - eta) * (y_corner2 - y_corner1) + - (1.0 + eta) * (y_corner3 - y_corner4))) + 0.25f0 * ((1 - eta) * (y_corner2 - y_corner1) + + (1 + eta) * (y_corner3 - y_corner4))) - Y_eta = (0.5 * ((1.0 - xi) * y4_prime + (1.0 + xi) * y2_prime + y3 - y1) + Y_eta = (0.5f0 * ((1 - xi) * y4_prime + (1 + xi) * y2_prime + y3 - y1) - - 0.25 * ((1.0 - xi) * (y_corner4 - y_corner1) + - (1.0 + xi) * (y_corner3 - y_corner2))) + 0.25f0 * ((1 - xi) * (y_corner4 - y_corner1) + + (1 + xi) * (y_corner3 - y_corner2))) return X_xi, X_eta, Y_xi, Y_eta end @@ -127,14 +125,14 @@ function calc_normal_directions!(normal_directions, element, nodes, # normal directions on the boundary for the left (local side 4) and right (local side 2) for j in eachindex(nodes) # side 2 - X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(1.0, nodes[j], + X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(1, nodes[j], surface_curves) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, j, 2, element] = sign(Jtemp) * (Y_eta) normal_directions[2, j, 2, element] = sign(Jtemp) * (-X_eta) # side 4 - X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(-1.0, nodes[j], + X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(-1, nodes[j], surface_curves) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, j, 4, element] = -sign(Jtemp) * (Y_eta) @@ -144,14 +142,14 @@ function calc_normal_directions!(normal_directions, element, nodes, # normal directions on the boundary for the top (local side 3) and bottom (local side 1) for i in eachindex(nodes) # side 1 - X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(nodes[i], -1.0, + X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(nodes[i], -1, surface_curves) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, i, 1, element] = -sign(Jtemp) * (-Y_xi) normal_directions[2, i, 1, element] = -sign(Jtemp) * (X_xi) # side 3 - X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(nodes[i], 1.0, + X_xi, X_eta, Y_xi, Y_eta = transfinite_quad_map_metrics(nodes[i], 1, surface_curves) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, i, 3, element] = sign(Jtemp) * (-Y_xi) diff --git a/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl b/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl index 7ceba93188..4f2f4e6f4b 100644 --- a/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl +++ b/src/solvers/dgsem_unstructured/mappings_geometry_straight_2d.jl @@ -9,15 +9,15 @@ # in physical coordinate space for a quadrilateral element with straight sides # Alg. 95 from the blue book of Kopriva function straight_side_quad_map(xi, eta, corner_points) - x = 0.25 * (corner_points[1, 1] * (1.0 - xi) * (1.0 - eta) - + corner_points[2, 1] * (1.0 + xi) * (1.0 - eta) - + corner_points[3, 1] * (1.0 + xi) * (1.0 + eta) - + corner_points[4, 1] * (1.0 - xi) * (1.0 + eta)) + x = 0.25f0 * (corner_points[1, 1] * (1 - xi) * (1 - eta) + + corner_points[2, 1] * (1 + xi) * (1 - eta) + + corner_points[3, 1] * (1 + xi) * (1 + eta) + + corner_points[4, 1] * (1 - xi) * (1 + eta)) - y = 0.25 * (corner_points[1, 2] * (1.0 - xi) * (1.0 - eta) - + corner_points[2, 2] * (1.0 + xi) * (1.0 - eta) - + corner_points[3, 2] * (1.0 + xi) * (1.0 + eta) - + corner_points[4, 2] * (1.0 - xi) * (1.0 + eta)) + y = 0.25f0 * (corner_points[1, 2] * (1 - xi) * (1 - eta) + + corner_points[2, 2] * (1 + xi) * (1 - eta) + + corner_points[3, 2] * (1 + xi) * (1 + eta) + + corner_points[4, 2] * (1 - xi) * (1 + eta)) return x, y end @@ -25,17 +25,17 @@ end # Compute the metric terms for the straight sided quadrilateral mapping # Alg. 100 from the blue book of Kopriva function straight_side_quad_map_metrics(xi, eta, corner_points) - X_xi = 0.25 * ((1.0 - eta) * (corner_points[2, 1] - corner_points[1, 1]) + - (1.0 + eta) * (corner_points[3, 1] - corner_points[4, 1])) + X_xi = 0.25f0 * ((1 - eta) * (corner_points[2, 1] - corner_points[1, 1]) + + (1 + eta) * (corner_points[3, 1] - corner_points[4, 1])) - X_eta = 0.25 * ((1.0 - xi) * (corner_points[4, 1] - corner_points[1, 1]) + - (1.0 + xi) * (corner_points[3, 1] - corner_points[2, 1])) + X_eta = 0.25f0 * ((1 - xi) * (corner_points[4, 1] - corner_points[1, 1]) + + (1 + xi) * (corner_points[3, 1] - corner_points[2, 1])) - Y_xi = 0.25 * ((1.0 - eta) * (corner_points[2, 2] - corner_points[1, 2]) + - (1.0 + eta) * (corner_points[3, 2] - corner_points[4, 2])) + Y_xi = 0.25f0 * ((1 - eta) * (corner_points[2, 2] - corner_points[1, 2]) + + (1 + eta) * (corner_points[3, 2] - corner_points[4, 2])) - Y_eta = 0.25 * ((1.0 - xi) * (corner_points[4, 2] - corner_points[1, 2]) + - (1.0 + xi) * (corner_points[3, 2] - corner_points[2, 2])) + Y_eta = 0.25f0 * ((1 - xi) * (corner_points[4, 2] - corner_points[1, 2]) + + (1 + xi) * (corner_points[3, 2] - corner_points[2, 2])) return X_xi, X_eta, Y_xi, Y_eta end @@ -78,14 +78,14 @@ function calc_normal_directions!(normal_directions, element, nodes, corners) # normal directions on the boundary for the left (local side 4) and right (local side 2) for j in eachindex(nodes) # side 2 - X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(1.0, nodes[j], + X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(1, nodes[j], corners) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, j, 2, element] = sign(Jtemp) * (Y_eta) normal_directions[2, j, 2, element] = sign(Jtemp) * (-X_eta) # side 4 - X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(-1.0, nodes[j], + X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(-1, nodes[j], corners) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, j, 4, element] = -sign(Jtemp) * (Y_eta) @@ -95,14 +95,14 @@ function calc_normal_directions!(normal_directions, element, nodes, corners) # normal directions on the boundary for the top (local side 3) and bottom (local side 1) for i in eachindex(nodes) # side 1 - X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], -1.0, + X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], -1, corners) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, i, 1, element] = -sign(Jtemp) * (-Y_xi) normal_directions[2, i, 1, element] = -sign(Jtemp) * (X_xi) # side 3 - X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], 1.0, + X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], 1, corners) Jtemp = X_xi * Y_eta - X_eta * Y_xi normal_directions[1, i, 3, element] = sign(Jtemp) * (-Y_xi) diff --git a/src/visualization/utilities.jl b/src/visualization/utilities.jl index c1108128c9..1f843c6a9d 100644 --- a/src/visualization/utilities.jl +++ b/src/visualization/utilities.jl @@ -14,7 +14,7 @@ # using the [Shoelace_formula](https://en.wikipedia.org/wiki/Shoelace_formula). function compute_triangle_area(tri) A, B, C = tri - return 0.5 * (A[1] * (B[2] - C[2]) + B[1] * (C[2] - A[2]) + C[1] * (A[2] - B[2])) + return 0.5f0 * (A[1] * (B[2] - C[2]) + B[1] * (C[2] - A[2]) + C[1] * (A[2] - B[2])) end # reference_plotting_triangulation(reference_plotting_coordinates) diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 7b69d1c0cf..3432bd69b2 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -257,6 +257,35 @@ end end end +@trixi_testset "elixir_euler_free_stream_boundaries_float32.jl" begin + # Expected errors are taken from elixir_euler_free_stream_boundaries.jl + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_free_stream_boundaries_float32.jl"), + l2=[ + Float32(6.530157034651212e-16), + Float32(1.6057829680004379e-15), + Float32(3.31107455378537e-15), + Float32(3.908829498281281e-15), + Float32(5.048390610424672e-15), + ], + linf=[ + Float32(4.884981308350689e-15), + Float32(1.1921019726912618e-14), + Float32(1.5432100042289676e-14), + Float32(2.298161660974074e-14), + Float32(6.039613253960852e-14), + ], + RealT=Float32) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_free_stream_extruded.jl with HLLC FLux" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"), l2=[ diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index c22a8e038a..82c76cc5d2 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -29,6 +29,23 @@ isdir(outdir) && rm(outdir, recursive = true) end end +@trixi_testset "elixir_advection_float32.jl" begin + # Expected errors are taken from elixir_advection_basic.jl + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_float32.jl"), + # Expected errors are taken from elixir_advection_basic.jl + l2=[Float32(8.311947673061856e-6)], + linf=[Float32(6.627000273229378e-5)], + RealT=Float32) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_advection_coupled.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_coupled.jl"), l2=[ diff --git a/test/test_trixi.jl b/test/test_trixi.jl index 9bfd73ea28..be3521ed16 100644 --- a/test/test_trixi.jl +++ b/test/test_trixi.jl @@ -4,8 +4,8 @@ import Trixi # Use a macro to avoid world age issues when defining new initial conditions etc. # inside an elixir. """ - @test_trixi_include(elixir; l2=nothing, linf=nothing, - atol=500*eps(), rtol=sqrt(eps()), + @test_trixi_include(elixir; l2=nothing, linf=nothing, RealT=Float64, + atol=500*eps(RealT), rtol=sqrt(eps(RealT)), parameters...) Test Trixi by calling `trixi_include(elixir; parameters...)`. @@ -17,8 +17,16 @@ as absolute/relative tolerance. macro test_trixi_include(elixir, args...) local l2 = get_kwarg(args, :l2, nothing) local linf = get_kwarg(args, :linf, nothing) - local atol = get_kwarg(args, :atol, 500 * eps()) - local rtol = get_kwarg(args, :rtol, sqrt(eps())) + local RealT = get_kwarg(args, :RealT, :Float64) + if RealT === :Float64 + atol_default = 500 * eps(Float64) + rtol_default = sqrt(eps(Float64)) + elseif RealT === :Float32 + atol_default = 500 * eps(Float32) + rtol_default = sqrt(eps(Float32)) + end + local atol = get_kwarg(args, :atol, atol_default) + local rtol = get_kwarg(args, :rtol, rtol_default) local skip_coverage = get_kwarg(args, :skip_coverage, false) local coverage_override = expr_to_named_tuple(get_kwarg(args, :coverage_override, :())) if !(:maxiters in keys(coverage_override)) @@ -33,7 +41,8 @@ macro test_trixi_include(elixir, args...) local kwargs = Pair{Symbol, Any}[] for arg in args if (arg.head == :(=) && - !(arg.args[1] in (:l2, :linf, :atol, :rtol, :coverage_override, :skip_coverage)) + !(arg.args[1] in (:l2, :linf, :RealT, :atol, :rtol, :coverage_override, + :skip_coverage)) && !(coverage && arg.args[1] in keys(coverage_override))) push!(kwargs, Pair(arg.args...)) end diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 563f99792d..43b911cc66 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -323,6 +323,33 @@ end end end +@trixi_testset "elixir_shallowwater_ec_float32.jl" begin + # Expected errors are nearly all taken from elixir_shallowwater_ec.jl + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec_float32.jl"), + l2=[ + Float32(0.6107326269462766), + Float32(0.48666631722018877), + Float32(0.48309775159067053), + Float32(0.29467422718511704), + ], + linf=[ + Float32(2.776782342826098), + 3.2162943f0, # this needs to be adapted + 3.6683278f0, # this needed to be adapted + Float32(2.052861364219655), + ], + tspan=(0.0f0, 0.25f0), + RealT=Float32) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_shallowwater_well_balanced.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ @@ -713,6 +740,23 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "FDSBP (upwind): elixir_euler_free_stream_upwind_float32.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_free_stream_upwind_float32.jl"), + l2=[0, 0, 0, 0], + linf=[0, 0, 0, 0], + tspan=(0.0f0, 0.05f0), + atol=9.0f-4) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory From 8d7b20aae40fc19428b07f616aee2e755ba6beb8 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 22 Aug 2024 14:12:30 +0200 Subject: [PATCH 5/6] adapt auxiliary math functions for `Float32` (#2048) * adapt auxiliary math functions * adapt tests * more tests and fix shock capturing volume integral * fix typo in tests * p4est_2d_dgsem/elixir_euler_shockcapturing_ec_float32.jl * format * fix stolarsky_mean test * use Float32 for test tolerances * format * increase tolerance for CI * format --- .../elixir_euler_shockcapturing_ec_float32.jl | 66 ++++ src/auxiliary/math.jl | 53 ++- src/solvers/dgmulti/shock_capturing.jl | 6 +- src/solvers/dgsem_p4est/dg_2d.jl | 2 +- src/solvers/dgsem_tree/dg.jl | 6 +- src/solvers/dgsem_unstructured/dg_2d.jl | 2 + test/test_p4est_2d.jl | 28 ++ test/test_type.jl | 349 +++++------------- 8 files changed, 224 insertions(+), 288 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec_float32.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec_float32.jl b/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec_float32.jl new file mode 100644 index 0000000000..4f390f4fc2 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec_float32.jl @@ -0,0 +1,66 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations2D(1.4f0) + +initial_condition = initial_condition_weak_blast_wave + +surface_flux = flux_ranocha +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(Float32, polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 1.0f0, + alpha_min = 0.001f0, + alpha_smooth = true, + variable = density_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral, RealT = Float32) + +############################################################################### + +coordinates_min = (-1.0f0, -1.0f0) +coordinates_max = (+1.0f0, +1.0f0) + +trees_per_dimension = (4, 4) +mesh = P4estMesh(trees_per_dimension, + polydeg = 4, initial_refinement_level = 2, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + periodicity = true, RealT = Float32) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0f0, 2.0f0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.0f0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + stepsize_callback) +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0f0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/auxiliary/math.jl b/src/auxiliary/math.jl index 171c631a84..fa816da9a1 100644 --- a/src/auxiliary/math.jl +++ b/src/auxiliary/math.jl @@ -124,7 +124,7 @@ end end """ - ln_mean(x, y) + Trixi.ln_mean(x::Real, y::Real) Compute the logarithmic mean @@ -171,18 +171,24 @@ Given ε = 1.0e-4, we use the following algorithm. for Intel, AMD, and VIA CPUs. [https://www.agner.org/optimize/instruction_tables.pdf](https://www.agner.org/optimize/instruction_tables.pdf) """ -@inline function ln_mean(x, y) - epsilon_f2 = 1.0e-4 +@inline ln_mean(x::Real, y::Real) = ln_mean(promote(x, y)...) + +@inline function ln_mean(x::RealT, y::RealT) where {RealT <: Real} + epsilon_f2 = convert(RealT, 1.0e-4) f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2 if f2 < epsilon_f2 - return (x + y) / @evalpoly(f2, 2, 2/3, 2/5, 2/7) + return (x + y) / @evalpoly(f2, + 2, + convert(RealT, 2 / 3), + convert(RealT, 2 / 5), + convert(RealT, 2 / 7)) else return (y - x) / log(y / x) end end """ - inv_ln_mean(x, y) + Trixi.inv_ln_mean(x::Real, y::Real) Compute the inverse `1 / ln_mean(x, y)` of the logarithmic mean [`ln_mean`](@ref). @@ -191,11 +197,17 @@ This function may be used to increase performance where the inverse of the logarithmic mean is needed, by replacing a (slow) division by a (fast) multiplication. """ -@inline function inv_ln_mean(x, y) - epsilon_f2 = 1.0e-4 +@inline inv_ln_mean(x::Real, y::Real) = inv_ln_mean(promote(x, y)...) + +@inline function inv_ln_mean(x::RealT, y::RealT) where {RealT <: Real} + epsilon_f2 = convert(RealT, 1.0e-4) f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2 if f2 < epsilon_f2 - return @evalpoly(f2, 2, 2/3, 2/5, 2/7) / (x + y) + return @evalpoly(f2, + 2, + convert(RealT, 2 / 3), + convert(RealT, 2 / 5), + convert(RealT, 2 / 7)) / (x + y) else return log(y / x) / (y - x) end @@ -273,7 +285,7 @@ end # [Fortran](https://godbolt.org/z/Yrsa1js7P) # or [C++](https://godbolt.org/z/674G7Pccv). """ - max(x, y, ...) + Trixi.max(x, y, ...) Return the maximum of the arguments. See also the `maximum` function to take the maximum element from a collection. @@ -292,7 +304,7 @@ julia> max(2, 5, 1) @inline max(args...) = @fastmath max(args...) """ - min(x, y, ...) + Trixi.min(x, y, ...) Return the minimum of the arguments. See also the `minimum` function to take the minimum element from a collection. @@ -311,7 +323,7 @@ julia> min(2, 5, 1) @inline min(args...) = @fastmath min(args...) """ - positive_part(x) + Trixi.positive_part(x) Return `x` if `x` is positive, else zero. In other words, return `(x + abs(x)) / 2` for real numbers `x`. @@ -321,7 +333,7 @@ Return `x` if `x` is positive, else zero. In other words, return end """ - negative_part(x) + Trixi.negative_part(x) Return `x` if `x` is negative, else zero. In other words, return `(x - abs(x)) / 2` for real numbers `x`. @@ -331,7 +343,7 @@ Return `x` if `x` is negative, else zero. In other words, return end """ - stolarsky_mean(x, y, gamma) + Trixi.stolarsky_mean(x::Real, y:Real, gamma::Real) Compute an instance of a weighted Stolarsky mean of the form @@ -382,15 +394,18 @@ Given ε = 1.0e-4, we use the following algorithm. for Intel, AMD, and VIA CPUs. [https://www.agner.org/optimize/instruction_tables.pdf](https://www.agner.org/optimize/instruction_tables.pdf) """ -@inline function stolarsky_mean(x, y, gamma) - epsilon_f2 = 1.0e-4 +@inline stolarsky_mean(x::Real, y::Real, gamma::Real) = stolarsky_mean(promote(x, y, + gamma)...) + +@inline function stolarsky_mean(x::RealT, y::RealT, gamma::RealT) where {RealT <: Real} + epsilon_f2 = convert(RealT, 1.0e-4) f2 = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y) # f2 = f^2 if f2 < epsilon_f2 # convenience coefficients - c1 = (1 / 3) * (gamma - 2) - c2 = -(1 / 15) * (gamma + 1) * (gamma - 3) * c1 - c3 = -(1 / 21) * (2 * gamma * (gamma - 2) - 9) * c2 - return 0.5 * (x + y) * @evalpoly(f2, 1, c1, c2, c3) + c1 = convert(RealT, 1 / 3) * (gamma - 2) + c2 = convert(RealT, -1 / 15) * (gamma + 1) * (gamma - 3) * c1 + c3 = convert(RealT, -1 / 21) * (2 * gamma * (gamma - 2) - 9) * c2 + return 0.5f0 * (x + y) * @evalpoly(f2, 1, c1, c2, c3) else return (gamma - 1) / gamma * (y^gamma - x^gamma) / (y^(gamma - 1) - x^(gamma - 1)) diff --git a/src/solvers/dgmulti/shock_capturing.jl b/src/solvers/dgmulti/shock_capturing.jl index d224e5ed03..99261b82ed 100644 --- a/src/solvers/dgmulti/shock_capturing.jl +++ b/src/solvers/dgmulti/shock_capturing.jl @@ -160,10 +160,14 @@ function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, mesh::DGMultiMesh, dg::DGMulti) empty!(element_ids_dg) empty!(element_ids_dgfv) + # For `Float64`, this gives 1.8189894035458565e-12 + # For `Float32`, this gives 1.1920929f-5 + RealT = eltype(alpha) + atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0)) for element in eachelement(mesh, dg) # Clip blending factor for values close to zero (-> pure DG) - dg_only = isapprox(alpha[element], 0, atol = 1e-12) + dg_only = isapprox(alpha[element], 0, atol = atol) if dg_only push!(element_ids_dg, element) else diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 29e94f4370..17b9af0446 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -379,7 +379,7 @@ end # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, node_index, direction_index, element_index] = flux_[v] + - 0.5 * + 0.5f0 * noncons_[v] end end diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 0993b3c9b8..febfb0b121 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -24,10 +24,14 @@ function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, cache) empty!(element_ids_dg) empty!(element_ids_dgfv) + # For `Float64`, this gives 1.8189894035458565e-12 + # For `Float32`, this gives 1.1920929f-5 + RealT = eltype(alpha) + atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0)) for element in eachelement(dg, cache) # Clip blending factor for values close to zero (-> pure DG) - dg_only = isapprox(alpha[element], 0, atol = 1e-12) + dg_only = isapprox(alpha[element], 0, atol = atol) if dg_only push!(element_ids_dg, element) else diff --git a/src/solvers/dgsem_unstructured/dg_2d.jl b/src/solvers/dgsem_unstructured/dg_2d.jl index 645bf05f07..a531ded91c 100644 --- a/src/solvers/dgsem_unstructured/dg_2d.jl +++ b/src/solvers/dgsem_unstructured/dg_2d.jl @@ -20,6 +20,8 @@ function create_cache(mesh::UnstructuredMesh2D, equations, # perform a check on the sufficient metric identities condition for free-stream preservation # and halt computation if it fails + # For `Float64`, this gives 1.8189894035458565e-12 + # For `Float32`, this gives 1.1920929f-5 atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0)) if !isapprox(max_discrete_metric_identities(dg, cache), 0, atol = atol) error("metric terms fail free-stream preservation check with maximum error $(max_discrete_metric_identities(dg, cache))") diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index f56dc9cd76..95c69b8695 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -221,6 +221,34 @@ end end end +@trixi_testset "elixir_euler_shockcapturing_ec_float32.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_shockcapturing_ec_float32.jl"), + l2=[ + 0.09539953f0, + 0.10563527f0, + 0.105637245f0, + 0.3507514f0, + ], + linf=[ + 0.2942562f0, + 0.4079147f0, + 0.3972956f0, + 1.0810697f0, + ], + tspan=(0.0f0, 1.0f0), + rtol=10 * sqrt(eps(Float32)), # to make CI pass + RealT=Float32) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + @trixi_testset "elixir_euler_sedov.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"), l2=[ diff --git a/test/test_type.jl b/test/test_type.jl index d13b626b06..bd553ee89c 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -11,6 +11,20 @@ isdir(outdir) && rm(outdir, recursive = true) # Run unit tests for various equations @testset "Test Type Stability" begin + @timed_testset "mean values" begin + for RealT1 in (Float32, Float64), RealT2 in (Float32, Float64) + RealT = promote_type(RealT1, RealT2) + @test typeof(@inferred Trixi.ln_mean(RealT1(1), RealT2(2))) == RealT + @test typeof(@inferred Trixi.inv_ln_mean(RealT1(1), RealT2(2))) == RealT + for RealT3 in (Float32, Float64) + RealT = promote_type(RealT1, RealT2, RealT3) + @test typeof(@inferred Trixi.stolarsky_mean(RealT1(1), RealT2(2), + RealT3(3))) == + RealT + end + end + end + @timed_testset "Acoustic Perturbation 2D" begin for RealT in (Float32, Float64) v_mean_global = (zero(RealT), zero(RealT)) @@ -121,25 +135,11 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux_kennedy_gruber(u_ll, u_rr, 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_chandrashekar(u_ll, u_rr, orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == - RealT - end - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == - RealT - end + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, equations))) == @@ -228,26 +228,11 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux_hllc(u_ll, u_rr, normal_direction, equations)) == RealT - if RealT == Float32 - # check `ln_mean` (test broken) - @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, - normal_direction, - equations)) == - RealT - else - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, - equations)) == - RealT - end - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, - equations)) == - RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, - equations)) == RealT - end + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, + equations)) == + RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT @test eltype(eltype(@inferred splitting_drikakis_tsangaris(u, normal_direction, equations))) == RealT @@ -282,26 +267,11 @@ isdir(outdir) && rm(outdir, recursive = true) 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_chandrashekar(u_ll, u_rr, - orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == - RealT - end - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == - RealT - end + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, equations))) == @@ -351,7 +321,7 @@ isdir(outdir) && rm(outdir, recursive = true) @timed_testset "Compressible Euler 3D" begin for RealT in (Float32, Float64) - # set gamma = 2 for the coupling convergence test + # set gamma = 2 for the coupling convergence test equations = @inferred CompressibleEulerEquations3D(RealT(2)) x = SVector(zero(RealT), zero(RealT), zero(RealT)) @@ -400,24 +370,10 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux_hllc(u_ll, u_rr, normal_direction, equations)) == RealT - if RealT == Float32 - # check `ln_mean` (test broken) - @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, - normal_direction, - equations)) == RealT - else - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, - equations)) == RealT - end - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, - equations)) == - RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, - equations)) == RealT - end + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, normal_direction, + equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == @@ -440,24 +396,10 @@ isdir(outdir) && rm(outdir, recursive = true) 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_chandrashekar(u_ll, u_rr, - orientation, - equations)) == RealT - else - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == RealT - end - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == - RealT - end + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT @test eltype(eltype(@inferred splitting_steger_warming(u, orientation, equations))) == @@ -509,23 +451,10 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux(u, orientation, equations)) == RealT - if RealT == Float32 - # check `ln_mean` (test broken) - @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == RealT - else - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == RealT - end - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == - RealT - end + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -565,34 +494,15 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, - equations)) == RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, - equations)) == RealT - end + @test eltype(@inferred flux_ranocha(u_ll, u_rr, normal_direction, + equations)) == RealT for orientation in orientations @test eltype(@inferred flux(u, orientation, equations)) == RealT - if RealT == Float32 - # check `ln_mean` (test broken) - @test_broken eltype(@inferred flux_chandrashekar(u_ll, u_rr, - orientation, - equations)) == RealT - else - @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, - equations)) == RealT - end - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, - equations)) == RealT - else - @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == - RealT - end + @test eltype(@inferred flux_chandrashekar(u_ll, u_rr, orientation, + equations)) == RealT + @test eltype(@inferred flux_ranocha(u_ll, u_rr, orientation, equations)) == + RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == @@ -637,15 +547,8 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, normal_rr, equations)) == RealT - if RealT == Float32 - # check `ln_mean` and `inv_ln_mean` (test broken) - @test_broken eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, equations)) == - RealT - end + @test eltype(@inferred flux_chan_etal(u_ll, u_rr, orientation, equations)) == + RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -989,21 +892,11 @@ isdir(outdir) && rm(outdir, recursive = true) @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 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 @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -1073,15 +966,8 @@ isdir(outdir) && rm(outdir, recursive = true) 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 eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, + equations)) == RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @@ -1101,22 +987,10 @@ isdir(outdir) && rm(outdir, recursive = true) 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 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 for nonconservative_term in nonconservative_terms @test eltype(@inferred flux_nonconservative_powell_local_symmetric(u_ll, orientation, @@ -1210,15 +1084,8 @@ isdir(outdir) && rm(outdir, recursive = true) 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 eltype(@inferred flux_hindenlang_gassner(u_ll, u_rr, normal_direction, + equations)) == RealT @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @@ -1233,22 +1100,10 @@ isdir(outdir) && rm(outdir, recursive = true) @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 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 @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -1318,21 +1173,10 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux(u, 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 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 @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == RealT @@ -1381,22 +1225,11 @@ isdir(outdir) && rm(outdir, recursive = true) @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 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 @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, equations)) == @@ -2011,17 +1844,9 @@ isdir(outdir) && rm(outdir, recursive = true) RealT @test eltype(@inferred flux(u, normal_direction, equations)) == RealT - if RealT == Float32 - # check `ln_mean` and `stolarsky_mean` (test broken) - @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr, - normal_direction, - equations)) == - RealT - else - @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction, - equations)) == - RealT - end + @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction, + equations)) == + RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction, equations)) == RealT @@ -2034,17 +1859,9 @@ isdir(outdir) && rm(outdir, recursive = true) for orientation in orientations @test eltype(@inferred flux(u, orientation, equations)) == RealT - if RealT == Float32 - # check `ln_mean` and `stolarsky_mean` (test broken) - @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr, - orientation, - equations)) == - RealT - else - @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation, - equations)) == - RealT - end + @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation, + equations)) == + RealT @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation, equations)) == RealT @@ -2124,7 +1941,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred dissipation(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred numflux(u_ll, u_rr, orientation, equations)) == RealT - # no matching method + # no matching method # @test eltype(@inferred numflux(u_ll, u_rr, normal_direction, equations)) == RealT @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, orientation, equations)) == RealT From e6ba83857eb4447aaccc6574b7a3675bbd21f9ae Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Thu, 22 Aug 2024 18:39:06 -1000 Subject: [PATCH 6/6] Add BC tests for type stability (#2045) * start * BC tests for Navier Stokes 1D * BC tests for Navier Stokes 2D * BC tests for Navier Stokes 3D * BC tests for Laplace Diffusion 1D --- test/test_type.jl | 279 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 261 insertions(+), 18 deletions(-) diff --git a/test/test_type.jl b/test/test_type.jl index bd553ee89c..d507cded95 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -315,7 +315,13 @@ isdir(outdir) && rm(outdir, recursive = true) @test typeof(@inferred Trixi.entropy_thermodynamic(cons, equations)) == RealT @test typeof(@inferred energy_internal(cons, equations)) == RealT - # TODO: test `gradient_conservative`, not necessary but good to have + @test eltype(@inferred Trixi.gradient_conservative(pressure, u, equations)) == + RealT + @test eltype(@inferred Trixi.gradient_conservative(Trixi.entropy_math, u, + equations)) == RealT + @test eltype(@inferred Trixi.gradient_conservative(Trixi.entropy_guermond_etal, + u, + equations)) == RealT end end @@ -577,11 +583,33 @@ isdir(outdir) && rm(outdir, recursive = true) Prandtl = prandtl_number, gradient_variables = GradientVariablesEntropy()) - u = u_transformed = SVector(one(RealT), zero(RealT), - zero(RealT)) + x = SVector(zero(RealT)) + t = zero(RealT) + u = u_inner = u_transformed = flux_inner = SVector(one(RealT), zero(RealT), + zero(RealT)) orientation = 1 + directions = [1, 2] gradients = SVector(RealT(0.1), RealT(0.1), RealT(0.1)) + operator_gradient = Trixi.Gradient() + operator_divergence = Trixi.Divergence() + + # For BC tests + function initial_condition_navier_stokes_convergence_test(x, t, equations) + RealT_local = eltype(x) + A = 0.5f0 + c = 2 + + pi_x = convert(RealT_local, pi) * x[1] + pi_t = convert(RealT_local, pi) * t + + rho = c + A * cos(pi_x) * cos(pi_t) + v1 = log(x[1] + 2) * (1 - exp(-A * (x[1] - 1))) * cos(pi_t) + p = rho^2 + + return prim2cons(SVector(rho, v1, p), equations) + end + for equations_parabolic in (equations_parabolic_primitive, equations_parabolic_entropy) @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == @@ -599,10 +627,49 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred Trixi.convert_derivative_to_primitive(u, gradients, equations_parabolic)) == RealT - end - # TODO: BC tests for GradientVariablesPrimitive - # TODO: BC tests for GradientVariablesEntropy + # For BC tests + velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, + t, + equations)[2]) + heat_bc_left = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, + t, + equations), + equations_parabolic)) + heat_bc_right = Adiabatic((x, t, equations) -> oftype(t, 0)) + + boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right, + heat_bc_left) + boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, + heat_bc_right) + # BC tests + for direction in directions + @test eltype(@inferred boundary_condition_right(flux_inner, u_inner, + orientation, direction, + x, + t, operator_gradient, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_right(flux_inner, u_inner, + orientation, direction, + x, + t, operator_divergence, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_left(flux_inner, u_inner, + orientation, direction, + x, + t, operator_gradient, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_left(flux_inner, u_inner, + orientation, direction, + x, + t, operator_divergence, + equations_parabolic)) == + RealT + end + end end end @@ -620,11 +687,37 @@ isdir(outdir) && rm(outdir, recursive = true) Prandtl = prandtl_number, gradient_variables = GradientVariablesEntropy()) - u = u_transformed = SVector(one(RealT), zero(RealT), zero(RealT), zero(RealT)) + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = w_inner = u_transformed = flux_inner = normal = 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) + operator_gradient = Trixi.Gradient() + operator_divergence = Trixi.Divergence() + + # For BC tests + function initial_condition_navier_stokes_convergence_test(x, t, equations) + RealT_local = eltype(x) + A = 0.5f0 + c = 2 + + pi_x = convert(RealT_local, pi) * x[1] + pi_y = convert(RealT_local, pi) * x[2] + pi_t = convert(RealT_local, pi) * t + + rho = c + A * sin(pi_x) * cos(pi_y) * cos(pi_t) + v1 = sin(pi_x) * log(x[2] + 2) * (1 - exp(-A * (x[2] - 1))) * cos(pi_t) + v2 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, p), equations) + end + for equations_parabolic in (equations_parabolic_primitive, equations_parabolic_entropy) for orientation in orientations @@ -648,10 +741,50 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred Trixi.convert_derivative_to_primitive(u, gradient, equations_parabolic)) == RealT - end - # TODO: BC tests for GradientVariablesPrimitive - # TODO: BC tests for GradientVariablesEntropy + # For BC tests + velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, + t, + equations)[2:3]) + heat_bc_left = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, + t, + equations), + equations_parabolic)) + heat_bc_right = Adiabatic((x, t, equations) -> oftype(t, 0)) + + boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right, + heat_bc_left) + boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, + heat_bc_right) + + # BC tests + @test eltype(@inferred boundary_condition_right(flux_inner, w_inner, + normal, + x, + t, + operator_gradient, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_right(flux_inner, w_inner, + normal, + x, + t, + operator_divergence, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_left(flux_inner, w_inner, + normal, + x, + t, operator_gradient, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_left(flux_inner, w_inner, + normal, + x, + t, operator_divergence, + equations_parabolic)) == + RealT + end end end @@ -669,12 +802,43 @@ isdir(outdir) && rm(outdir, recursive = true) Prandtl = prandtl_number, gradient_variables = GradientVariablesEntropy()) - u = u_transformed = SVector(one(RealT), zero(RealT), zero(RealT), zero(RealT), - zero(RealT)) + x = SVector(zero(RealT), zero(RealT), zero(RealT)) + t = zero(RealT) + u = w_inner = u_transformed = flux_inner = normal = 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) + operator_gradient = Trixi.Gradient() + operator_divergence = Trixi.Divergence() + + # For BC tests + function initial_condition_navier_stokes_convergence_test(x, t, equations) + RealT_local = eltype(x) + c = 2 + A1 = 0.5f0 + A2 = 1 + A3 = 0.5f0 + + pi_x = convert(RealT_local, pi) * x[1] + pi_y = convert(RealT_local, pi) * x[2] + pi_z = convert(RealT_local, pi) * x[3] + pi_t = convert(RealT_local, pi) * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2) * (1 - exp(-A3 * (x[2] - 1))) * + sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) + end + for equations_parabolic in (equations_parabolic_primitive, equations_parabolic_entropy) for orientation in orientations @@ -698,10 +862,50 @@ isdir(outdir) && rm(outdir, recursive = true) @test eltype(@inferred Trixi.convert_derivative_to_primitive(u, gradient, equations_parabolic)) == RealT - end - # TODO: BC tests for GradientVariablesPrimitive - # TODO: BC tests for GradientVariablesEntropy + # For BC tests + velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, + t, + equations)[2:4]) + heat_bc_left = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x, + t, + equations), + equations_parabolic)) + heat_bc_right = Adiabatic((x, t, equations) -> oftype(t, 0)) + + boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right, + heat_bc_left) + boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right, + heat_bc_right) + + # BC tests + @test eltype(@inferred boundary_condition_right(flux_inner, w_inner, + normal, + x, + t, + operator_gradient, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_right(flux_inner, w_inner, + normal, + x, + t, + operator_divergence, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_left(flux_inner, w_inner, + normal, + x, + t, operator_gradient, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_left(flux_inner, w_inner, + normal, + x, + t, operator_divergence, + equations_parabolic)) == + RealT + end end end @@ -1299,14 +1503,53 @@ isdir(outdir) && rm(outdir, recursive = true) x = SVector(zero(RealT)) t = zero(RealT) - u = gradients = SVector(one(RealT)) + u = u_inner = flux_inner = normal = gradients = SVector(one(RealT)) orientation = 1 + operator_gradient = Trixi.Gradient() + operator_divergence = Trixi.Divergence() + @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) == RealT - # TODO: BC tests for BoundaryConditionDirichlet - # TODO: BC tests for BoundaryConditionNeumann + # For BC tests + function initial_condition_convergence_test(x, t, + equation::LaplaceDiffusion1D) + RealT_local = eltype(x) + x_trans = x[1] - equation.diffusivity * t + + c = 1 + A = 0.5f0 + L = 2 + f = 1.0f0 / L + omega = 2 * convert(RealT_local, pi) * f + scalar = c + A * sin(omega * sum(x_trans)) + return SVector(scalar) + end + + boundary_condition_dirichlet = BoundaryConditionDirichlet(initial_condition_convergence_test) + boundary_condition_neumann = BoundaryConditionNeumann((x, t, equations) -> oftype(t, + 0)) + + # BC tests + @test eltype(@inferred boundary_condition_dirichlet(flux_inner, u_inner, normal, + x, t, + operator_gradient, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_dirichlet(flux_inner, u_inner, normal, + x, t, + operator_divergence, + equations_parabolic)) == + RealT + @test eltype(@inferred boundary_condition_neumann(flux_inner, u_inner, normal, + x, t, + operator_gradient, + equations_parabolic)) == RealT + @test eltype(@inferred boundary_condition_neumann(flux_inner, u_inner, normal, + x, t, + operator_divergence, + equations_parabolic)) == RealT end end