Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

computations using only Float32 and no Float64 #2042

Merged
merged 21 commits into from
Aug 21, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using Downloads: download

using OrdinaryDiffEq
using Trixi

Expand Down
9 changes: 4 additions & 5 deletions examples/p4est_2d_dgsem/elixir_euler_NACA6412airfoil_mach2.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@

using Trixi
using OrdinaryDiffEq
using Downloads: download

###############################################################################
# semidiscretization of the compressible Euler equations
Expand Down Expand Up @@ -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,
Expand Down
2 changes: 1 addition & 1 deletion examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using Downloads: download

using OrdinaryDiffEq
using Trixi

Expand Down
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using Downloads: download

using OrdinaryDiffEq
using Trixi

Expand All @@ -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)

Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@

using Downloads: download
using OrdinaryDiffEq
using Trixi

Expand All @@ -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]
Expand Down
Original file line number Diff line number Diff line change
@@ -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
61 changes: 61 additions & 0 deletions examples/structured_2d_dgsem/elixir_advection_float32.jl
Original file line number Diff line number Diff line change
@@ -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()
124 changes: 124 additions & 0 deletions examples/unstructured_2d_dgsem/elixir_shallowwater_ec_float32.jl
Original file line number Diff line number Diff line change
@@ -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
Loading
Loading