Arbitrary Precision LGL Basis (#2128)

* test higher precision * try to get higher order convergence * type stability * fmt * conservation * correct * working version * changes * rm ODE * rm unused * Apply suggestions from code review Co-authored-by: Joshua Lampert <[email protected]> * RealT for TreeMesh * undo tree * Revert "undo tree" This reverts commit 8103d0a. * ones, zeros * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> Co-authored-by: Joshua Lampert <[email protected]> * reprod test vals f32 * LV suggestion * typo * shorten * data safety * typesafety * fix * docstring * comments * Type general projection * error * syntax * bring back conv * tolerance always dynamic * convert pi * tighter tolerances * warning only * update some float test vals * remove unused/NON_TESTED * restructure arguments * docstrings * remove random 1000 * old behaviour * compliance with deprecated dgsem * Avoid breaking change * avoid ambiguity * pass datatype through * do not cast to float * remvoe comment * more type generality * do not mix up things * move example * tests * News etc * fmt * update test vals to different basis * try to disable warnings * fmt * CI hazzle * fun with Ci * more warnings * Compat entries * Tolerances for quad precision types * Comment * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * comment * shorter code for analysis * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> --------- Co-authored-by: Joshua Lampert <[email protected]> Co-authored-by: Michael Schlottke-Lakemper <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]>
trixi-framework · Dec 2, 2024 · 5db379b · 5db379b
1 parent 50d30e9
commit 5db379b
Show file tree

Hide file tree

Showing 22 changed files with 330 additions and 139 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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