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 <[email protected]> * remove redundant using Downloads in examples * more comments * format --------- Co-authored-by: Andrew Winters <[email protected]>
trixi-framework · Aug 21, 2024 · 20ab97b · 20ab97b
1 parent afd8060
commit 20ab97b
Show file tree

Hide file tree

Showing 49 changed files with 737 additions and 287 deletions.
diff --git 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
@@ -1,7 +1,6 @@
 
 using Trixi
 using OrdinaryDiffEq
-using Downloads: download
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
@@ -62,24 +61,24 @@ 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)
 
 # Mesh generated from the following gmsh geometry input file:
 # 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]
 
 mesh = P4estMesh{2}(mesh_file, polydeg = polydeg, boundary_symbols = boundary_symbols)
 
 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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