Arbitrary Precision LGL Basis

examples/structured_1d_dgsem/elixir_advection_double_float.jl

@@ -0,0 +1,58 @@
+
+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 = (-RealT(1.0),) # minimum coordinate
+coordinates_max = (RealT(1.0),) # maximum coordinate
+cells_per_dimension = (256,)
+
+# NOTE: StructuredMesh supports higher precision 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), RealT(1) / 1000)
+
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 1000,
+                                     extra_analysis_errors = (:conservation_error,))
+
+# cfl does not need to be in higher precision
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        stepsize_callback,
+                        analysis_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(),
+            dt = 42.0, # `dt` does not need to be in higher precision
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

examples/structured_1d_dgsem/elixir_advection_float128.jl

@@ -0,0 +1,58 @@
+
+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 = 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 = (-RealT(1.0),) # minimum coordinate
+coordinates_max = (RealT(1.0),) # maximum coordinate
+cells_per_dimension = (256,)
+
+# NOTE: StructuredMesh supports higher precision 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), RealT(1) / 1000)
+
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+analysis_callback = AnalysisCallback(semi, interval = 1000,
+                                     extra_analysis_errors = (:conservation_error,))
+
+# cfl does not need to be in higher precision
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        stepsize_callback,
+                        analysis_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(),
+            dt = 42.0, # `dt` does not need to be in higher precision
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

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

src/callbacks_step/analysis_dg1d.jl

@@ -61,7 +61,10 @@ function calc_error_norms(func, u, t, analyzer,
                                        inv.(view(inverse_jacobian, :, element)))
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
+        #@. jacobian_local = abs(jacobian_local) # Does not work with LoopVectorization and higher precision datatypes
+        for (index, value) in enumerate(jacobian_local)
+            jacobian_local[index] = abs(value)
+        end
 
         for i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t,

src/callbacks_step/analysis_dg2d.jl

@@ -161,7 +161,10 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
+        #@. jacobian_local = abs(jacobian_local) # Does not work with LoopVectorization and higher precision datatypes
+        for (index, value) in enumerate(jacobian_local)
+            jacobian_local[index] = abs(value)
+        end
 
         for j in eachnode(analyzer), i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j),

src/callbacks_step/analysis_dg2d_parallel.jl

@@ -114,7 +114,10 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
+        #@. jacobian_local = abs(jacobian_local) # Does not work with LoopVectorization and higher precision datatypes
+        for (index, value) in enumerate(jacobian_local)
+            jacobian_local[index] = abs(value)
+        end
 
         for j in eachnode(analyzer), i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i, j),

src/callbacks_step/analysis_dg3d.jl

@@ -187,7 +187,10 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1, jacobian_tmp2)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
+        #@. jacobian_local = abs(jacobian_local) # Does not work with LoopVectorization and higher precision datatypes
+        for (index, value) in enumerate(jacobian_local)
+            jacobian_local[index] = abs(value)
+        end
 
         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,

src/callbacks_step/analysis_dg3d_parallel.jl

@@ -31,7 +31,10 @@ function calc_error_norms(func, u, t, analyzer,
                                        jacobian_tmp1, jacobian_tmp2)
 
         # Calculate errors at each analysis node
-        @. jacobian_local = abs(jacobian_local)
+        #@. jacobian_local = abs(jacobian_local) # Does not work with LoopVectorization and higher precision datatypes
+        for (index, value) in enumerate(jacobian_local)
+            jacobian_local[index] = abs(value)
+        end
 
         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,