diff --git a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl index 74a6ba74fb2..0221817e1cf 100644 --- a/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl +++ b/examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl @@ -56,12 +56,12 @@ save_solution = SaveSolutionCallback(dt = 0.1, # Construct embedded order paired explicit Runge-Kutta method with 10 stages and 7 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) +ode_algorithm = Trixi.EmbeddedPairedRK3(16, 5, 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 = 0.75) # Warisa: This number is quite small in contrast the other one from optimizing A +stepsize_callback = StepsizeCallback() # 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, @@ -83,9 +83,9 @@ summary_callback() function construct_b_vector(b_unknown, num_stages_embedded, num_stage_evals_embedded) # Construct the b vector b = [ - 1 - sum(b_unknown), + b_unknown[1], zeros(Float64, num_stages_embedded - num_stage_evals_embedded)..., - b_unknown..., + b_unknown[2:end]..., 0 ] return b diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl index 0b00b1bf364..d859014d712 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_embedded_PERK3.jl @@ -7,6 +7,45 @@ using DelimitedFiles: readdlm @muladd begin #! format: noindent +function compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c) + # Different scheme for c? + + A = zeros(num_stage_evals - 1, num_stage_evals - 1) + b_embedded = zeros(num_stage_evals - 1) + rhs = [1, 1/2, embedded_monomial_coeffs...] + + # sum(b) = 1 + A[1, :] .= 1 + + # The second order constraint: dot(b,c) = 1/2 + for i in 2:num_stage_evals - 1 + A[2, i] = c[num_stages - num_stage_evals + i] + end + + # Fill the A matrix + for i in 3:(num_stage_evals - 1) + # z^i + for j in i: (num_stage_evals - 1) + println("i = ", i, ", j = ", j) + println("[num_stages - num_stage_evals + j - 1] = ", num_stages - num_stage_evals + j - 1) + A[i,j] = c[num_stages - num_stage_evals + j - 1] + # number of times a_unknown should be multiplied in each power of z + for k in 1: i-2 + # so we want to get from a[k] * ... i-2 times (1 time is already accounted for by c-coeff) + # j-1 - k + 1 = j - k + println("a_unknown at index: ", j - k) + A[i, j] *= a_unknown[j - k] # a[k] is in the same stage as b[k-1] -> since we also store b_1 + end + end + #rhs[i] /= factorial(i) + end + + display(A) + + b_embedded = A \ rhs + return b_embedded +end + # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 3 # using a list of eigenvalues function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, tspan, @@ -42,7 +81,7 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, dtmax, dteps, eig_vals; verbose) - embedded_monomial_coeffs = undo_normalization!(embedded_monomial_coeffs, consistency_order, num_stage_evals-1) + embedded_monomial_coeffs = undo_normalization!(embedded_monomial_coeffs, consistency_order-1, num_stage_evals-1) # Solve the nonlinear system of equations from monomial coefficient and # Butcher array abscissae c to find Butcher matrix A @@ -54,13 +93,14 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals, println("dt_opt_b = ", dt_opt_b) - b = Trixi.solve_b_butcher_coeffs_unknown!(b, num_stages, num_stage_evals, embedded_monomial_coeffs, c, a_unknown; verbose) + b = compute_b_embedded_coeffs(num_stage_evals, num_stages, embedded_monomial_coeffs, a_unknown, c) + println("b: ", b) println("Embedded monomial after normalization: ", embedded_monomial_coeffs) println("dt_opt_a = ", dt_opt_a) - #error("Stability polynomial optimization process is complete.") + #error("b found.") end # Fill A-matrix in P-ERK style a_matrix = zeros(num_stage_evals - 2, 2) @@ -309,7 +349,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator) @threaded for i in eachindex(integrator.u) # add the contribution of the first stage - integrator.u[i] += (1 - sum(alg.b)) * integrator.k1[i] + integrator.u[i] += alg.b[1] * integrator.k1[i] end # Higher stages where the weight of b in the butcher tableau is zero @@ -345,7 +385,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator) @threaded for i in eachindex(integrator.u) integrator.u[i] += integrator.k_higher[i] * - alg.b[1] + alg.b[2] end # Higher stages after num_stage_evals where b is non-zero @@ -370,10 +410,11 @@ function step!(integrator::EmbeddedPairedRK3Integrator) @threaded for i in eachindex(integrator.u) integrator.u[i] += integrator.k_higher[i] * - alg.b[stage - alg.num_stages + alg.num_stage_evals - 1] + alg.b[stage - alg.num_stages + alg.num_stage_evals] end end + # Last stage @threaded for i in eachindex(integrator.du) integrator.u_tmp[i] = integrator.u[i] + @@ -382,6 +423,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator) alg.a_matrix[alg.num_stage_evals - 2, 2] * integrator.k_higher[i] end + integrator.f(integrator.du, integrator.u_tmp, prob.p, integrator.t + alg.c[alg.num_stages] * integrator.dt)