diff --git a/NEWS.md b/NEWS.md index 2f71c45e57e..c187f229519 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,8 @@ for human readability. - New time integrator `PairedExplicitRK3`, implementing the third-order paired explicit Runge-Kutta method with [Convex.jl](https://github.com/jump-dev/Convex.jl), [ECOS.jl](https://github.com/jump-dev/ECOS.jl), and [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) ([#2008]) +- `LobattoLegendreBasis` and related datastructures made fully floating-type general, + enabling calculations with higher than double (`Float64`) precision ([#2128]) ## Changes when updating to v0.9 from v0.8.x diff --git a/README.md b/README.md index 70c28799f31..5b58612bb09 100644 --- a/README.md +++ b/README.md @@ -33,6 +33,7 @@ installation and postprocessing procedures. Its features include: * Hierarchical quadtree/octree grid with adaptive mesh refinement * Forests of quadtrees/octrees with [p4est](https://github.com/cburstedde/p4est) via [P4est.jl](https://github.com/trixi-framework/P4est.jl) * High-order accuracy in space and time + * Arbitrary floating-point precision * Discontinuous Galerkin methods * Kinetic energy-preserving and entropy-stable methods based on flux differencing * Entropy-stable shock capturing diff --git a/docs/src/index.md b/docs/src/index.md index 88752cc7f11..40191e97b9a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -26,6 +26,7 @@ installation and postprocessing procedures. Its features include: * Hierarchical quadtree/octree grid with adaptive mesh refinement * Forests of quadtrees/octrees with [p4est](https://github.com/cburstedde/p4est) via [P4est.jl](https://github.com/trixi-framework/P4est.jl) * High-order accuracy in space and time + * Arbitrary floating-point precision * Discontinuous Galerkin methods * Kinetic energy-preserving and entropy-stable methods based on flux differencing * Entropy-stable shock capturing diff --git a/examples/structured_1d_dgsem/elixir_advection_float128.jl b/examples/structured_1d_dgsem/elixir_advection_float128.jl new file mode 100644 index 00000000000..7f72eba58e3 --- /dev/null +++ b/examples/structured_1d_dgsem/elixir_advection_float128.jl @@ -0,0 +1,60 @@ + +using OrdinaryDiffEq +using Trixi + +using Quadmath + +############################################################################### +# semidiscretization of the linear advection equation + +# See https://github.com/JuliaMath/Quadmath.jl +RealT = Float128 + +advection_velocity = 4 / 3 # Does not need to be in higher precision +equations = LinearScalarAdvectionEquation1D(advection_velocity) + +solver = DGSEM(RealT = RealT, polydeg = 13, surface_flux = flux_lax_friedrichs) + +# CARE: Important to use higher precision datatype for coordinates +# as these are used for type promotion of the mesh (points etc.) +coordinates_min = (-one(RealT),) +coordinates_max = (one(RealT),) +cells_per_dimension = (1,) + +# `StructuredMesh` infers datatype from coordinates +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# CARE: Important to use higher precision datatype in specification of final time +tspan = (zero(RealT), one(RealT)) + +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() + +analysis_callback = AnalysisCallback(semi, interval = 100, + extra_analysis_errors = (:conservation_error,)) + +# cfl does not need to be in higher precision +stepsize_callback = StepsizeCallback(cfl = 0.25) + +callbacks = CallbackSet(summary_callback, + stepsize_callback, + analysis_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, Feagin14(), + # Turn off adaptivity to avoid setting very small tolerances + adaptive = false, + dt = 42, # `dt` does not need to be in higher precision + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/examples/tree_1d_dgsem/elixir_advection_doublefloat.jl b/examples/tree_1d_dgsem/elixir_advection_doublefloat.jl new file mode 100644 index 00000000000..8405d0125d0 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_advection_doublefloat.jl @@ -0,0 +1,63 @@ + +using OrdinaryDiffEq +using Trixi + +using DoubleFloats + +############################################################################### +# semidiscretization of the linear advection equation + +# See https://github.com/JuliaMath/DoubleFloats.jl +RealT = Double64 + +advection_velocity = 4 / 3 # Does not need to be in higher precision +equations = LinearScalarAdvectionEquation1D(advection_velocity) + +solver = DGSEM(RealT = RealT, polydeg = 7, surface_flux = flux_lax_friedrichs) + +# CARE: Important to use higher precision datatype for coordinates +# as these are used for type promotion of the mesh (points etc.) +coordinates_min = -one(RealT) # minimum coordinate +coordinates_max = one(RealT) # maximum coordinate + +# For `TreeMesh` the datatype has to be specified explicitly, +# i.e., is not inferred from the coordinates. +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 30_000, + RealT = RealT) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver) + +############################################################################### +# ODE solvers, callbacks etc. + +# CARE: Important to use higher precision datatype in specification of final time +tspan = (zero(RealT), one(RealT)) + +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() + +analysis_callback = AnalysisCallback(semi, interval = 100, + extra_analysis_errors = (:conservation_error,)) + +# cfl does not need to be in higher precision +stepsize_callback = StepsizeCallback(cfl = 1.4) + +callbacks = CallbackSet(summary_callback, + stepsize_callback, + analysis_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, DP8(), + # Turn off adaptivity to avoid setting very small tolerances + adaptive = false, + dt = 42, # `dt` does not need to be in higher precision + save_everystep = false, callback = callbacks); + +# Print the timer summary +summary_callback() diff --git a/src/auxiliary/precompile.jl b/src/auxiliary/precompile.jl index dfda8ece687..5a1de036808 100644 --- a/src/auxiliary/precompile.jl +++ b/src/auxiliary/precompile.jl @@ -329,8 +329,10 @@ function _precompile_manual_() @assert Base.precompile(Tuple{typeof(Trixi.gauss_nodes_weights), Int}) @assert Base.precompile(Tuple{typeof(Trixi.calc_forward_upper), Int}) @assert Base.precompile(Tuple{typeof(Trixi.calc_forward_lower), Int}) - @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_upper), Int, Val{:gauss}}) - @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_lower), Int, Val{:gauss}}) + @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_upper), Int, + Val{:gauss}}) + @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_lower), Int, + Val{:gauss}}) @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_upper), Int, Val{:gauss_lobatto}}) @assert Base.precompile(Tuple{typeof(Trixi.calc_reverse_lower), Int, diff --git a/src/auxiliary/special_elixirs.jl b/src/auxiliary/special_elixirs.jl index d71a27aa96a..3e960945e7a 100644 --- a/src/auxiliary/special_elixirs.jl +++ b/src/auxiliary/special_elixirs.jl @@ -6,10 +6,11 @@ #! format: noindent """ - convergence_test([mod::Module=Main,] elixir::AbstractString, iterations; kwargs...) + convergence_test([mod::Module=Main,] elixir::AbstractString, iterations, RealT = Float64; kwargs...) Run `iterations` Trixi.jl simulations using the setup given in `elixir` and compute the experimental order of convergence (EOC) in the ``L^2`` and ``L^\\infty`` norm. +Use `RealT` as the data type to represent the errors. In each iteration, the resolution of the respective mesh will be doubled. Additional keyword arguments `kwargs...` and the optional module `mod` are passed directly to [`trixi_include`](@ref). @@ -18,12 +19,14 @@ This function assumes that the spatial resolution is set via the keywords `initial_refinement_level` (an integer) or `cells_per_dimension` (a tuple of integers, one per spatial dimension). """ -function convergence_test(mod::Module, elixir::AbstractString, iterations; kwargs...) +function convergence_test(mod::Module, elixir::AbstractString, iterations, + RealT = Float64; + kwargs...) @assert(iterations>1, "Number of iterations must be bigger than 1 for a convergence analysis") # Types of errors to be calculated - errors = Dict(:l2 => Float64[], :linf => Float64[]) + errors = Dict(:l2 => RealT[], :linf => RealT[]) initial_resolution = extract_initial_resolution(elixir, kwargs) @@ -105,7 +108,7 @@ function analyze_convergence(errors, iterations, println("") # Print mean EOCs - mean_values = zeros(nvariables) + mean_values = zeros(eltype(errors[:l2]), nvariables) for v in 1:nvariables mean_values[v] = sum(eocs[kind][:, v]) ./ length(eocs[kind][:, v]) @printf("%-10s", "mean") @@ -119,8 +122,9 @@ function analyze_convergence(errors, iterations, return eoc_mean_values end -function convergence_test(elixir::AbstractString, iterations; kwargs...) - convergence_test(Main, elixir::AbstractString, iterations; kwargs...) +function convergence_test(elixir::AbstractString, iterations, RealT = Float64; + kwargs...) + convergence_test(Main, elixir::AbstractString, iterations, RealT; kwargs...) end # Helper methods used in the functions defined above diff --git a/src/callbacks_step/analysis_dg1d.jl b/src/callbacks_step/analysis_dg1d.jl index d2613c325be..b4d49a39d5a 100644 --- a/src/callbacks_step/analysis_dg1d.jl +++ b/src/callbacks_step/analysis_dg1d.jl @@ -61,16 +61,17 @@ function calc_error_norms(func, u, t, analyzer, inv.(view(inverse_jacobian, :, element))) # Calculate errors at each analysis node - @. jacobian_local = abs(jacobian_local) - for i in eachnode(analyzer) u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t, equations) diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i), equations) - l2_error += diff .^ 2 * (weights[i] * jacobian_local[i]) + # We take absolute value as we need the Jacobian here for the volume calculation + abs_jacobian_local_i = abs(jacobian_local[i]) + + l2_error += diff .^ 2 * (weights[i] * abs_jacobian_local_i) linf_error = @. max(linf_error, abs(diff)) - total_volume += weights[i] * jacobian_local[i] + total_volume += weights[i] * abs_jacobian_local_i end end diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index e089133fa17..e33533738ec 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -161,16 +161,17 @@ function calc_error_norms(func, u, t, analyzer, jacobian_tmp1) # Calculate errors at each analysis node - @. jacobian_local = abs(jacobian_local) - for j in eachnode(analyzer), i in eachnode(analyzer) u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j), t, equations) diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j), equations) - l2_error += diff .^ 2 * (weights[i] * weights[j] * jacobian_local[i, j]) + # We take absolute value as we need the Jacobian here for the volume calculation + abs_jacobian_local_ij = abs(jacobian_local[i, j]) + + l2_error += diff .^ 2 * (weights[i] * weights[j] * abs_jacobian_local_ij) linf_error = @. max(linf_error, abs(diff)) - total_volume += weights[i] * weights[j] * jacobian_local[i, j] + total_volume += weights[i] * weights[j] * abs_jacobian_local_ij end end diff --git a/src/callbacks_step/analysis_dg2d_parallel.jl b/src/callbacks_step/analysis_dg2d_parallel.jl index 5b3ae858ab7..10b6e30f95c 100644 --- a/src/callbacks_step/analysis_dg2d_parallel.jl +++ b/src/callbacks_step/analysis_dg2d_parallel.jl @@ -114,16 +114,17 @@ function calc_error_norms(func, u, t, analyzer, jacobian_tmp1) # Calculate errors at each analysis node - @. jacobian_local = abs(jacobian_local) - for j in eachnode(analyzer), i in eachnode(analyzer) u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j), t, equations) diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j), equations) - l2_error += diff .^ 2 * (weights[i] * weights[j] * jacobian_local[i, j]) + # We take absolute value as we need the Jacobian here for the volume calculation + abs_jacobian_local_ij = abs(jacobian_local[i, j]) + + l2_error += diff .^ 2 * (weights[i] * weights[j] * abs_jacobian_local_ij) linf_error = @. max(linf_error, abs(diff)) - volume += weights[i] * weights[j] * jacobian_local[i, j] + volume += weights[i] * weights[j] * abs_jacobian_local_ij end end diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 3c9c3a69f5e..efc1dc3bb3c 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -187,18 +187,19 @@ function calc_error_norms(func, u, t, analyzer, jacobian_tmp1, jacobian_tmp2) # Calculate errors at each analysis node - @. jacobian_local = abs(jacobian_local) - for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer) u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j, k), t, equations) diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j, k), equations) + # We take absolute value as we need the Jacobian here for the volume calculation + abs_jacobian_local_ijk = abs(jacobian_local[i, j, k]) + l2_error += diff .^ 2 * - (weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k]) + (weights[i] * weights[j] * weights[k] * abs_jacobian_local_ijk) linf_error = @. max(linf_error, abs(diff)) - total_volume += weights[i] * weights[j] * weights[k] * - jacobian_local[i, j, k] + total_volume += (weights[i] * weights[j] * weights[k] * + abs_jacobian_local_ijk) end end diff --git a/src/callbacks_step/analysis_dg3d_parallel.jl b/src/callbacks_step/analysis_dg3d_parallel.jl index 285f0f62de6..ab9ba6a0055 100644 --- a/src/callbacks_step/analysis_dg3d_parallel.jl +++ b/src/callbacks_step/analysis_dg3d_parallel.jl @@ -31,17 +31,18 @@ function calc_error_norms(func, u, t, analyzer, jacobian_tmp1, jacobian_tmp2) # Calculate errors at each analysis node - @. jacobian_local = abs(jacobian_local) - for k in eachnode(analyzer), j in eachnode(analyzer), i in eachnode(analyzer) u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j, k), t, equations) diff = func(u_exact, equations) - func(get_node_vars(u_local, equations, dg, i, j, k), equations) + # We take absolute value as we need the Jacobian here for the volume calculation + abs_jacobian_local_ijk = abs(jacobian_local[i, j, k]) + l2_error += diff .^ 2 * - (weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k]) + (weights[i] * weights[j] * weights[k] * abs_jacobian_local_ijk) linf_error = @. max(linf_error, abs(diff)) - volume += weights[i] * weights[j] * weights[k] * jacobian_local[i, j, k] + volume += weights[i] * weights[j] * weights[k] * abs_jacobian_local_ijk end end diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl index cac1dba9c74..3ea668c0264 100644 --- a/src/solvers/dgsem/basis_lobatto_legendre.jl +++ b/src/solvers/dgsem/basis_lobatto_legendre.jl @@ -6,7 +6,7 @@ #! format: noindent """ - LobattoLegendreBasis([RealT=Float64,] polydeg::Integer) + LobattoLegendreBasis([RealT = Float64,] polydeg::Integer) Create a nodal Lobatto-Legendre basis for polynomials of degree `polydeg`. @@ -37,15 +37,14 @@ end function LobattoLegendreBasis(RealT, polydeg::Integer) nnodes_ = polydeg + 1 - # compute everything using `Float64` by default - nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_) + nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT) inverse_weights_ = inv.(weights_) - _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_) + _, inverse_vandermonde_legendre_ = vandermonde_legendre(nodes_, RealT) - boundary_interpolation_ = zeros(nnodes_, 2) - boundary_interpolation_[:, 1] = calc_lhat(-1.0, nodes_, weights_) - boundary_interpolation_[:, 2] = calc_lhat(1.0, nodes_, weights_) + boundary_interpolation_ = zeros(RealT, nnodes_, 2) + boundary_interpolation_[:, 1] = calc_lhat(-one(RealT), nodes_, weights_) + boundary_interpolation_[:, 2] = calc_lhat(one(RealT), nodes_, weights_) derivative_matrix_ = polynomial_derivative_matrix(nodes_) derivative_split_ = calc_dsplit(nodes_, weights_) @@ -79,7 +78,6 @@ function LobattoLegendreBasis(RealT, polydeg::Integer) derivative_split_transpose, derivative_dhat) end - LobattoLegendreBasis(polydeg::Integer) = LobattoLegendreBasis(Float64, polydeg) function Base.show(io::IO, basis::LobattoLegendreBasis) @@ -165,11 +163,10 @@ function MortarL2(basis::LobattoLegendreBasis) RealT = real(basis) nnodes_ = nnodes(basis) - # compute everything using `Float64` by default - forward_upper_ = calc_forward_upper(nnodes_) - forward_lower_ = calc_forward_lower(nnodes_) - reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss)) - reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss)) + forward_upper_ = calc_forward_upper(nnodes_, RealT) + forward_lower_ = calc_forward_lower(nnodes_, RealT) + reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss), RealT) + reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss), RealT) # type conversions to get the requested real type and enable possible # optimizations of runtime performance and latency @@ -228,10 +225,10 @@ end # end # function MortarEC(basis::LobattoLegendreBasis{RealT}, surface_flux) -# forward_upper = calc_forward_upper(n_nodes) -# forward_lower = calc_forward_lower(n_nodes) -# l2reverse_upper = calc_reverse_upper(n_nodes, Val(:gauss_lobatto)) -# l2reverse_lower = calc_reverse_lower(n_nodes, Val(:gauss_lobatto)) +# forward_upper = calc_forward_upper(n_nodes, RealT) +# forward_lower = calc_forward_lower(n_nodes, RealT) +# l2reverse_upper = calc_reverse_upper(n_nodes, Val(:gauss_lobatto), RealT) +# l2reverse_lower = calc_reverse_lower(n_nodes, Val(:gauss_lobatto), RealT) # # type conversions to make use of StaticArrays etc. # nnodes_ = nnodes(basis) @@ -263,8 +260,7 @@ function SolutionAnalyzer(basis::LobattoLegendreBasis; RealT = real(basis) nnodes_ = analysis_polydeg + 1 - # compute everything using `Float64` by default - nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_) + nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT) vandermonde_ = polynomial_interpolation_matrix(get_nodes(basis), nodes_) # type conversions to get the requested real type and enable possible @@ -322,11 +318,10 @@ end function AdaptorL2(basis::LobattoLegendreBasis{RealT}) where {RealT} nnodes_ = nnodes(basis) - # compute everything using `Float64` by default - forward_upper_ = calc_forward_upper(nnodes_) - forward_lower_ = calc_forward_lower(nnodes_) - reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss)) - reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss)) + forward_upper_ = calc_forward_upper(nnodes_, RealT) + forward_lower_ = calc_forward_lower(nnodes_, RealT) + reverse_upper_ = calc_reverse_upper(nnodes_, Val(:gauss), RealT) + reverse_lower_ = calc_reverse_lower(nnodes_, Val(:gauss), RealT) # type conversions to get the requested real type and enable possible # optimizations of runtime performance and latency @@ -408,7 +403,7 @@ end # This implements algorithm 37 "PolynomialDerivativeMatrix" from Kopriva's book. function polynomial_derivative_matrix(nodes) n_nodes = length(nodes) - d = zeros(n_nodes, n_nodes) + d = zeros(eltype(nodes), n_nodes, n_nodes) wbary = barycentric_weights(nodes) for i in 1:n_nodes, j in 1:n_nodes @@ -481,7 +476,7 @@ For details, see (especially Section 3) """ function barycentric_weights(nodes) n_nodes = length(nodes) - weights = ones(n_nodes) + weights = ones(eltype(nodes), n_nodes) for j in 2:n_nodes, k in 1:(j - 1) weights[k] *= nodes[k] - nodes[j] @@ -528,7 +523,7 @@ For details, see e.g. Section 2 of """ function lagrange_interpolating_polynomials(x, nodes, wbary) n_nodes = length(nodes) - polynomials = zeros(n_nodes) + polynomials = zeros(eltype(nodes), n_nodes) for i in 1:n_nodes # Avoid division by zero when `x` is close to node by using @@ -553,7 +548,7 @@ function lagrange_interpolating_polynomials(x, nodes, wbary) end """ - gauss_lobatto_nodes_weights(n_nodes::Integer) + gauss_lobatto_nodes_weights(n_nodes::Integer, RealT = Float64) Computes nodes ``x_j`` and weights ``w_j`` for the (Legendre-)Gauss-Lobatto quadrature. This implements algorithm 25 "GaussLobattoNodesAndWeights" from the book @@ -563,15 +558,14 @@ This implements algorithm 25 "GaussLobattoNodesAndWeights" from the book Algorithms for scientists and engineers. [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5) """ -# From FLUXO (but really from blue book by Kopriva) -function gauss_lobatto_nodes_weights(n_nodes::Integer) - # From Kopriva's book - n_iterations = 10 - tolerance = 1e-15 +function gauss_lobatto_nodes_weights(n_nodes::Integer, RealT = Float64) + # From FLUXO (but really from blue book by Kopriva) + n_iterations = 20 + tolerance = 2 * eps(RealT) # Relative tolerance for Newton iteration # Initialize output - nodes = zeros(n_nodes) - weights = zeros(n_nodes) + nodes = zeros(RealT, n_nodes) + weights = zeros(RealT, n_nodes) # Special case for polynomial degree zero (first order finite volume) if n_nodes == 1 @@ -584,15 +578,15 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer) N = n_nodes - 1 # Calculate values at boundary - nodes[1] = -1.0 - nodes[end] = 1.0 - weights[1] = 2 / (N * (N + 1)) + nodes[1] = -1 + nodes[end] = 1 + weights[1] = RealT(2) / (N * (N + 1)) weights[end] = weights[1] # Calculate interior values if N > 1 - cont1 = pi / N - cont2 = 3 / (8 * N * pi) + cont1 = convert(RealT, pi) / N + cont2 = 3 / (8 * N * convert(RealT, pi)) # Use symmetry -> only left side is computed for i in 1:(div(N + 1, 2) - 1) @@ -608,6 +602,10 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer) if abs(dx) < tolerance * abs(nodes[i + 1]) break end + + if k == n_iterations + @warn "`gauss_lobatto_nodes_weights` Newton iteration did not converge" + end end # Calculate weight @@ -622,8 +620,8 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer) # If odd number of nodes, set center node to origin (= 0.0) and calculate weight if n_nodes % 2 == 1 - _, _, L = calc_q_and_l(N, 0) - nodes[div(N, 2) + 1] = 0.0 + _, _, L = calc_q_and_l(N, zero(RealT)) + nodes[div(N, 2) + 1] = 0 weights[div(N, 2) + 1] = weights[1] / L^2 end @@ -631,11 +629,13 @@ function gauss_lobatto_nodes_weights(n_nodes::Integer) end # From FLUXO (but really from blue book by Kopriva, algorithm 24) -function calc_q_and_l(N::Integer, x::Float64) - L_Nm2 = 1.0 +function calc_q_and_l(N::Integer, x::Real) + RealT = typeof(x) + + L_Nm2 = one(RealT) L_Nm1 = x - Lder_Nm2 = 0.0 - Lder_Nm1 = 1.0 + Lder_Nm2 = zero(RealT) + Lder_Nm1 = one(RealT) local L for i in 2:N @@ -652,10 +652,9 @@ function calc_q_and_l(N::Integer, x::Float64) return q, qder, L end -calc_q_and_l(N::Integer, x::Real) = calc_q_and_l(N, convert(Float64, x)) """ - gauss_nodes_weights(n_nodes::Integer) + gauss_nodes_weights(n_nodes::Integer, RealT = Float64) Computes nodes ``x_j`` and weights ``w_j`` for the Gauss-Legendre quadrature. This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book @@ -665,31 +664,30 @@ This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book Algorithms for scientists and engineers. [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5) """ -function gauss_nodes_weights(n_nodes::Integer) - # From Kopriva's book - n_iterations = 10 - tolerance = 1e-15 +function gauss_nodes_weights(n_nodes::Integer, RealT = Float64) + n_iterations = 20 + tolerance = 2 * eps(RealT) # Relative tolerance for Newton iteration # Initialize output - nodes = ones(n_nodes) * 1000 - weights = zeros(n_nodes) + nodes = ones(RealT, n_nodes) + weights = zeros(RealT, n_nodes) # Get polynomial degree for convenience N = n_nodes - 1 if N == 0 - nodes .= 0.0 - weights .= 2.0 + nodes .= 0 + weights .= 2 return nodes, weights elseif N == 1 - nodes[1] = -sqrt(1 / 3) + nodes[1] = -sqrt(one(RealT) / 3) nodes[end] = -nodes[1] - weights .= 1.0 + weights .= 1 return nodes, weights else # N > 1 # Use symmetry property of the roots of the Legendre polynomials for i in 0:(div(N + 1, 2) - 1) # Starting guess for Newton method - nodes[i + 1] = -cos(pi / (2 * N + 2) * (2 * i + 1)) + nodes[i + 1] = -cos(convert(RealT, pi) / (2 * N + 2) * (2 * i + 1)) # Newton iteration to find root of Legendre polynomial (= integration node) for k in 0:n_iterations @@ -699,6 +697,10 @@ function gauss_nodes_weights(n_nodes::Integer) if abs(dx) < tolerance * abs(nodes[i + 1]) break end + + if k == n_iterations + @warn "`gauss_nodes_weights` Newton iteration did not converge" + end end # Calculate weight @@ -712,8 +714,8 @@ function gauss_nodes_weights(n_nodes::Integer) # If odd number of nodes, set center node to origin (= 0.0) and calculate weight if n_nodes % 2 == 1 - poly, deriv = legendre_polynomial_and_derivative(N + 1, 0.0) - nodes[div(N, 2) + 1] = 0.0 + poly, deriv = legendre_polynomial_and_derivative(N + 1, zero(RealT)) + nodes[div(N, 2) + 1] = 0 weights[div(N, 2) + 1] = (2 * N + 3) / deriv^2 end @@ -733,20 +735,21 @@ This implements algorithm 22 "LegendrePolynomialAndDerivative" from the book [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5) """ function legendre_polynomial_and_derivative(N::Int, x::Real) + RealT = typeof(x) if N == 0 - poly = 1.0 - deriv = 0.0 + poly = one(RealT) + deriv = zero(RealT) elseif N == 1 - poly = convert(Float64, x) - deriv = 1.0 + poly = x + deriv = one(RealT) else - poly_Nm2 = 1.0 - poly_Nm1 = convert(Float64, x) - deriv_Nm2 = 0.0 - deriv_Nm1 = 1.0 + poly_Nm2 = one(RealT) + poly_Nm1 = copy(x) + deriv_Nm2 = zero(RealT) + deriv_Nm1 = one(RealT) - poly = 0.0 - deriv = 0.0 + poly = zero(RealT) + deriv = zero(RealT) for i in 2:N poly = ((2 * i - 1) * x * poly_Nm1 - (i - 1) * poly_Nm2) / i deriv = deriv_Nm2 + (2 * i - 1) * poly_Nm1 @@ -765,10 +768,10 @@ function legendre_polynomial_and_derivative(N::Int, x::Real) end # Calculate Legendre vandermonde matrix and its inverse -function vandermonde_legendre(nodes, N) +function vandermonde_legendre(nodes, N::Integer, RealT = Float64) n_nodes = length(nodes) n_modes = N + 1 - vandermonde = zeros(n_nodes, n_modes) + vandermonde = zeros(RealT, n_nodes, n_modes) for i in 1:n_nodes for m in 1:n_modes @@ -779,5 +782,7 @@ function vandermonde_legendre(nodes, N) inverse_vandermonde = inv(vandermonde) return vandermonde, inverse_vandermonde end -vandermonde_legendre(nodes) = vandermonde_legendre(nodes, length(nodes) - 1) +vandermonde_legendre(nodes, RealT = Float64) = vandermonde_legendre(nodes, + length(nodes) - 1, + RealT) end # @muladd diff --git a/src/solvers/dgsem/l2projection.jl b/src/solvers/dgsem/l2projection.jl index 0bb46f5ca15..436c1b92032 100644 --- a/src/solvers/dgsem/l2projection.jl +++ b/src/solvers/dgsem/l2projection.jl @@ -23,9 +23,9 @@ # Calculate forward projection matrix for discrete L2 projection from large to upper # # Note: This is actually an interpolation. -function calc_forward_upper(n_nodes) +function calc_forward_upper(n_nodes, RealT = Float64) # Calculate nodes, weights, and barycentric weights - nodes, weights = gauss_lobatto_nodes_weights(n_nodes) + nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT) wbary = barycentric_weights(nodes) # Calculate projection matrix (actually: interpolation) @@ -43,9 +43,9 @@ end # Calculate forward projection matrix for discrete L2 projection from large to lower # # Note: This is actually an interpolation. -function calc_forward_lower(n_nodes) +function calc_forward_lower(n_nodes, RealT = Float64) # Calculate nodes, weights, and barycentric weights - nodes, weights = gauss_lobatto_nodes_weights(n_nodes) + nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT) wbary = barycentric_weights(nodes) # Calculate projection matrix (actually: interpolation) @@ -64,9 +64,9 @@ end # # Note: To not make the L2 projection exact, first convert to Gauss nodes, # perform projection, and convert back to Gauss-Lobatto. -function calc_reverse_upper(n_nodes, ::Val{:gauss}) +function calc_reverse_upper(n_nodes, ::Val{:gauss}, RealT = Float64) # Calculate nodes, weights, and barycentric weights for Legendre-Gauss - gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes) + gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes, RealT) gauss_wbary = barycentric_weights(gauss_nodes) # Calculate projection matrix (actually: discrete L2 projection with errors) @@ -80,7 +80,7 @@ function calc_reverse_upper(n_nodes, ::Val{:gauss}) end # Calculate Vandermondes - lobatto_nodes, lobatto_weights = gauss_lobatto_nodes_weights(n_nodes) + lobatto_nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT) gauss2lobatto = polynomial_interpolation_matrix(gauss_nodes, lobatto_nodes) lobatto2gauss = polynomial_interpolation_matrix(lobatto_nodes, gauss_nodes) @@ -91,9 +91,9 @@ end # # Note: To not make the L2 projection exact, first convert to Gauss nodes, # perform projection, and convert back to Gauss-Lobatto. -function calc_reverse_lower(n_nodes, ::Val{:gauss}) +function calc_reverse_lower(n_nodes, ::Val{:gauss}, RealT = Float64) # Calculate nodes, weights, and barycentric weights for Legendre-Gauss - gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes) + gauss_nodes, gauss_weights = gauss_nodes_weights(n_nodes, RealT) gauss_wbary = barycentric_weights(gauss_nodes) # Calculate projection matrix (actually: discrete L2 projection with errors) @@ -107,7 +107,7 @@ function calc_reverse_lower(n_nodes, ::Val{:gauss}) end # Calculate Vandermondes - lobatto_nodes, lobatto_weights = gauss_lobatto_nodes_weights(n_nodes) + lobatto_nodes, _ = gauss_lobatto_nodes_weights(n_nodes, RealT) gauss2lobatto = polynomial_interpolation_matrix(gauss_nodes, lobatto_nodes) lobatto2gauss = polynomial_interpolation_matrix(lobatto_nodes, gauss_nodes) @@ -116,9 +116,9 @@ end # Calculate reverse projection matrix for discrete L2 projection from upper to large (Gauss-Lobatto # version) -function calc_reverse_upper(n_nodes, ::Val{:gauss_lobatto}) +function calc_reverse_upper(n_nodes, ::Val{:gauss_lobatto}, RealT = Float64) # Calculate nodes, weights, and barycentric weights - nodes, weights = gauss_lobatto_nodes_weights(n_nodes) + nodes, weights = gauss_lobatto_nodes_weights(n_nodes, RealT) wbary = barycentric_weights(nodes) # Calculate projection matrix (actually: discrete L2 projection with errors) @@ -135,9 +135,9 @@ end # Calculate reverse projection matrix for discrete L2 projection from lower to large (Gauss-Lobatto # version) -function calc_reverse_lower(n_nodes, ::Val{:gauss_lobatto}) +function calc_reverse_lower(n_nodes, ::Val{:gauss_lobatto}, RealT = Float64) # Calculate nodes, weights, and barycentric weights - nodes, weights = gauss_lobatto_nodes_weights(n_nodes) + nodes, weights = gauss_lobatto_nodes_weights(n_nodes, RealT) wbary = barycentric_weights(nodes) # Calculate projection matrix (actually: discrete L2 projection with errors) diff --git a/test/Project.toml b/test/Project.toml index 64c20ea4c37..f9a6241421b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Convex = "f65535da-76fb-5f13-bab9-19810c17039a" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" @@ -14,6 +15,7 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -23,6 +25,7 @@ Aqua = "0.8" CairoMakie = "0.10" Convex = "0.16" DelimitedFiles = "1" +DoubleFloats = "1.4.0" Downloads = "1" ECOS = "1.1.2" ExplicitImports = "1.0.1" @@ -34,6 +37,7 @@ NLsolve = "4.5.1" OrdinaryDiffEq = "6.49.1" Plots = "1.19" Printf = "1" +Quadmath = "0.5.10" Random = "1" StableRNGs = "1.0.2" Test = "1" diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 7a1a3f3ac09..0e70282e77c 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -231,10 +231,10 @@ end 0.3507514f0 ], linf=[ - 0.2942562f0, + 0.2930063f0, 0.4079147f0, 0.3972956f0, - 1.0810697f0 + 1.0764117f0 ], tspan=(0.0f0, 1.0f0), rtol=10 * sqrt(eps(Float32)), # to make CI pass diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl index ccb0c8cacac..a5b561ed8fe 100644 --- a/test/test_structured_1d.jl +++ b/test/test_structured_1d.jl @@ -58,6 +58,20 @@ end end end +@trixi_testset "elixir_advection_float128.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_float128.jl"), + l2=Float128[6.49879312655540217059228636803492411e-09], + linf=Float128[5.35548407857266390181158920649552284e-08]) + # 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 + # Testing the third-order paired explicit Runge-Kutta (PERK) method with its optimal CFL number @trixi_testset "elixir_burgers_perk3.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_perk3.jl"), diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl index f061e2e1c30..0665a1be105 100644 --- a/test/test_tree_1d_advection.jl +++ b/test/test_tree_1d_advection.jl @@ -142,6 +142,20 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000 end end + +@trixi_testset "elixir_advection_doublefloat.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_doublefloat.jl"), + l2=Double64[6.80895929885700039832943251427357703e-11], + linf=Double64[5.82834770064525291688100323411704252e-10]) + # 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 end # module diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index b7152cb8b75..5d52f4a5d42 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -150,24 +150,24 @@ end @trixi_testset "elixir_mhd_ec_float32.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_float32.jl"), - l2=Float32[0.036355644, + l2=Float32[0.03635566, + 0.042947732, 0.042947736, - 0.04294775, - 0.025748005, - 0.16211236, - 0.017452478, - 0.017452475, - 0.026877576, - 2.0951437f-7], - linf=Float32[0.22100884, - 0.28798944, - 0.2879896, - 0.20858234, - 0.8812654, - 0.09208202, - 0.09208143, - 0.14795381, - 2.0912607f-6]) + 0.025748001, + 0.16211228, + 0.01745248, + 0.017452491, + 0.026877586, + 2.417227f-7], + linf=Float32[0.2210092, + 0.28798974, + 0.28799006, + 0.20858109, + 0.8812673, + 0.09208107, + 0.09208131, + 0.14795369, + 2.2078211f-6]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_trixi.jl b/test/test_trixi.jl index be3521ed16a..024f3654ea0 100644 --- a/test/test_trixi.jl +++ b/test/test_trixi.jl @@ -15,6 +15,7 @@ are compared approximately against these reference values, using `atol, rtol` as absolute/relative tolerance. """ macro test_trixi_include(elixir, args...) + # Note: The variables below are just Symbols, not actual errors/types local l2 = get_kwarg(args, :l2, nothing) local linf = get_kwarg(args, :linf, nothing) local RealT = get_kwarg(args, :RealT, :Float64) @@ -24,6 +25,12 @@ macro test_trixi_include(elixir, args...) elseif RealT === :Float32 atol_default = 500 * eps(Float32) rtol_default = sqrt(eps(Float32)) + elseif RealT === :Float128 + atol_default = 500 * eps(Float128) + rtol_default = sqrt(eps(Float128)) + elseif RealT === :Double64 + atol_default = 500 * eps(Double64) + rtol_default = sqrt(eps(Double64)) end local atol = get_kwarg(args, :atol, atol_default) local rtol = get_kwarg(args, :rtol, rtol_default) @@ -167,7 +174,10 @@ macro test_nowarn_mod(expr, additional_ignore_content = String[]) r"WARNING: Method definition .* in module .* at .* overwritten .*.\n", # Warnings from third party packages r"┌ Warning: Problem status ALMOST_INFEASIBLE; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n", - r"┌ Warning: Problem status ALMOST_OPTIMAL; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n"] + r"┌ Warning: Problem status ALMOST_OPTIMAL; solution may be inaccurate.\n└ @ Convex ~/.julia/packages/Convex/.*\n", + # Warnings for higher-precision floating data types + r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/interpolation.jl:118 =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/tIJUA/src/condense_loopset.jl:1166\n", + r"┌ Warning: #= /home/runner/work/Trixi.jl/Trixi.jl/src/solvers/dgsem/interpolation.jl:136 =#:\n│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\n│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.\n└ @ Trixi ~/.julia/packages/LoopVectorization/tIJUA/src/condense_loopset.jl:1166\n"] append!(ignore_content, $additional_ignore_content) for pattern in ignore_content stderr_content = replace(stderr_content, pattern => "") diff --git a/test/test_unit.jl b/test/test_unit.jl index cd219050555..97e32726ca7 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -268,6 +268,12 @@ end @timed_testset "interpolation" begin @testset "nodes and weights" begin @test Trixi.gauss_nodes_weights(1) == ([0.0], [2.0]) + + @test Trixi.gauss_nodes_weights(2)[1] ≈ [-1 / sqrt(3), 1 / sqrt(3)] + @test Trixi.gauss_nodes_weights(2)[2] == [1.0, 1.0] + + @test Trixi.gauss_nodes_weights(3)[1] ≈ [-sqrt(3 / 5), 0.0, sqrt(3 / 5)] + @test Trixi.gauss_nodes_weights(3)[2] ≈ [5 / 9, 8 / 9, 5 / 9] end @testset "multiply_dimensionwise" begin diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index c433338a238..fed71e49f0b 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -335,8 +335,8 @@ end ], linf=[ Float32(2.776782342826098), - 3.2162943f0, # this needs to be adapted - 3.6683278f0, # this needed to be adapted + 3.2116454f0, # this needs to be adapted + 3.6616623f0, # this needed to be adapted Float32(2.052861364219655) ], tspan=(0.0f0, 0.25f0),