add some comments

examples/structured_1d_dgsem/elixir_euler_source_terms_embedded_perk3.jl

@@ -44,7 +44,9 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-
+# Construct embedded order paired explicit Runge-Kutta method with 10 stages and 6 evaluation stages for given simulation setup.
+# Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used 
+# in calculating the polynomial coefficients in the ODE algorithm.
 ode_algorithm = Trixi.EmbeddedPairedRK3(10, 6, tspan, semi)
 cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
 
@@ -64,4 +66,16 @@ sol = Trixi.solve(ode, ode_algorithm,
             dt = 1.0, # 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
+
+# Some function defined so that I can check if the second order condition is met. This will be removed later.
+function construct_b_vector(b_unknown, num_stages_embedded, num_stage_evals_embedded)
+    # Construct the b vector
+    b = [1 - sum(b_unknown), zeros(Float64, num_stages_embedded - num_stage_evals_embedded)..., b_unknown..., 0]
+    return b
+end
+
+b = construct_b_vector(ode_algorithm.b, ode_algorithm.num_stages - 1, ode_algorithm.num_stage_evals - 1)
+println("dot(b, c) = ", dot(b, ode_algorithm.c))
+println("sum(b) = ", sum(b))
+
 println("cfl_number = ", cfl_number)

examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl

@@ -62,7 +62,7 @@ ode_algorithm = Trixi.EmbeddedPairedRK3(10, 6, tspan, semi)
 # Calculate the CFL number for the given ODE algorithm and ODE problem (cfl_number calculate from dt_opt of the optimization of
 # b values in the Butcher tableau of the ODE algorithm).
 cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
-stepsize_callback = StepsizeCallback(cfl = cfl_number) # Warisa: This number is very small in contrast the other one from optimizing A
+stepsize_callback = StepsizeCallback(cfl = cfl_number) # Warisa: This number is quite small in contrast the other one from optimizing A
                                                        # I've tried using cfl of 1.5 and the error is very similar.
 
 callbacks = CallbackSet(summary_callback,
@@ -80,6 +80,7 @@ sol = Trixi.solve(ode, ode_algorithm,
 # Print the timer summary
 summary_callback()
 
+# Some function defined so that I can check if the second order condition is met. This will be removed later.
 function construct_b_vector(b_unknown, num_stages_embedded, num_stage_evals_embedded)
     # Construct the b vector
     b = [1 - sum(b_unknown), zeros(Float64, num_stages_embedded - num_stage_evals_embedded)..., b_unknown..., 0]

src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl

@@ -51,52 +51,59 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
     a_matrix[:, 2] = a_unknown
 
 
-    #TODO: something wrong when s = e
+    # Find the optimal b coeffficient from the Butcher tableau and its time step
     b, dt_opt_b = solve_b_butcher_coeffs_unknown(num_eig_vals, eig_vals,
                                            num_stages, num_stage_evals,
                                            num_stages - 1, # num_stages_embedded = num_stages - 1	
                                            num_stage_evals - 1, # num_stage_evals_embedded = num_stage_evals - 1
                                            a_unknown, c, dtmax, dteps)
 
 
-    return a_matrix, b, c, dt_opt_b
+    return a_matrix, b, c, dt_opt_b # Return the optimal time step from the b coefficients for testing purposes
 end
 
 # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3
 # using provided values of coefficients a in A-matrix of Butcher tableau
-function compute_EmbeddedPairedRK3_butcher_tableau(num_stages,
-                                                   base_path_a_coeffs::AbstractString;
+function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
+                                                   base_path_coeffs::AbstractString;
                                                    cS2)
 
     # Initialize array of c
     c = compute_c_coeffs(num_stages, cS2)
 
     # - 2 Since First entry of A is always zero (explicit method) and second is given by c_2 (consistency)
-    a_coeffs_max = num_stages - 2
+    coeffs_max = num_stage_evals - 2
 
-    #TODO: Read the values from the file
-    b = zeros(num_stages)
-
-    a_matrix = zeros(a_coeffs_max, 2)
+    a_matrix = zeros(coeffs_max, 2)
     a_matrix[:, 1] = c[3:end]
 
-    path_a_coeffs = joinpath(base_path_a_coeffs,
-                             "a_" * string(num_stages) * ".txt")
+    b = zeros(coeffs_max)
+
+    path_a_coeffs = joinpath(base_path_coeffs,
+                             "a_" * string(num_stages) * "_" * string(num_stage_evals) * ".txt")
+
+    path_b_coeffs = joinpath(base_path_coeffs,
+                             "b_" * string(num_stages) * "_" * string(num_stage_evals) * ".txt")
 
     @assert isfile(path_a_coeffs) "Couldn't find file $path_a_coeffs"
     a_coeffs = readdlm(path_a_coeffs, Float64)
     num_a_coeffs = size(a_coeffs, 1)
 
-    @assert num_a_coeffs == a_coeffs_max
+    @assert num_a_coeffs == coeffs_max
     # Fill A-matrix in P-ERK style
     a_matrix[:, 1] -= a_coeffs
     a_matrix[:, 2] = a_coeffs
 
+    @assert isfile(path_b_coeffs) "Couldn't find file $path_b_coeffs"
+    b = readdlm(path_b_coeffs, Float64)
+    num_b_coeffs = size(b, 1)
+    @assert num_b_coeffs == coeffs_max
+
     return a_matrix, b, c
 end
 
 @doc raw"""
-    EmbeddedPairedRK3(num_stages, base_path_a_coeffs::AbstractString, dt_opt;
+    EmbeddedPairedRK3(num_stages, base_path_coeffs::AbstractString, dt_opt;
                       cS2 = 1.0f0)
     EmbeddedPairedRK3(num_stages, tspan, semi::AbstractSemidiscretization;
                       verbose = false, cS2 = 1.0f0)
@@ -105,8 +112,8 @@ end
 
     Parameters:
     - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method.
-    - `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in 
-      the Butcher tableau of the Runge Kutta method.
+    - `base_path_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix and a file constaining 
+      in some coefficients in the b vector of the Butcher tableau of the Runge Kutta method.
       The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks.
     - `dt_opt` (`Float64`): Optimal time step size for the simulation setup.
     - `tspan`: Time span of the simulation.
@@ -141,11 +148,10 @@ end # struct EmbeddedPairedRK3
 
 # Constructor for previously computed A Coeffs
 function EmbeddedPairedRK3(num_stages, num_stage_evals,
-                           base_path_a_coeffs::AbstractString, dt_opt;
+                           base_path_coeffs::AbstractString, dt_opt;
                            cS2 = 1.0f0)
-    #TODO: Update this constructor to have num_stage_evals as well and also return correct things
-    a_matrix, b, c = compute_EmbeddedPairedRK3_butcher_tableau(num_stages,
-                                                               base_path_a_coeffs;
+    a_matrix, b, c = compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
+                                                               base_path_coeffs;
                                                                cS2)
 
     return EmbeddedPairedRK3(num_stages, num_stage_evals, a_matrix, b, c, dt_opt)