Skip to content

Commit

Permalink
change with the integrator
Browse files Browse the repository at this point in the history
  • Loading branch information
warisa-r committed Oct 25, 2024
1 parent b2358e5 commit 2af55dd
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 20 deletions.
4 changes: 2 additions & 2 deletions examples/tree_1d_dgsem/elixir_advection_embedded_perk3.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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, 7, tspan, semi)
ode_algorithm = Trixi.EmbeddedPairedRK3(10, 9, 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() # Warisa: This number is quite small in contrast the other one from optimizing A
stepsize_callback = StepsizeCallback()
# I've tried using cfl of 1.5 and the error is very similar.

callbacks = CallbackSet(summary_callback,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ function compute_EmbeddedPairedRK3_butcher_tableau(num_stages, num_stage_evals,
println("dt_opt_a = ", dt_opt_a)


error("b found.")
#error("b found.")
end
# Fill A-matrix in P-ERK style
a_matrix = zeros(num_stage_evals - 2, 2)
Expand Down Expand Up @@ -370,7 +370,7 @@ function step!(integrator::EmbeddedPairedRK3Integrator)
end

# Higher stages where the weight of b in the butcher tableau is zero
for stage in 2:(alg.num_stages - alg.num_stage_evals + 1)
for stage in 2:(alg.num_stages - alg.num_stage_evals)
# Construct current state
@threaded for i in eachindex(integrator.du)
integrator.u_tmp[i] = integrator.u[i] + alg.c[stage] * integrator.k1[i]
Expand All @@ -384,25 +384,28 @@ function step!(integrator::EmbeddedPairedRK3Integrator)
end
end

# k_(s-e+2)
# #k_(s-e+1) and k_(s-e+2)
# Construct current state
@threaded for i in eachindex(integrator.du)
integrator.u_tmp[i] = integrator.u[i] +
alg.c[alg.num_stages - alg.num_stage_evals + 2] *
integrator.k1[i]
end

integrator.f(integrator.du, integrator.u_tmp, prob.p,
integrator.t +
alg.c[alg.num_stages - alg.num_stage_evals + 2] * integrator.dt)
for j in 1:2
@threaded for i in eachindex(integrator.du)
integrator.u_tmp[i] = integrator.u[i] +
alg.c[alg.num_stages - alg.num_stage_evals + 2] *
integrator.k1[i]
end

@threaded for i in eachindex(integrator.du)
integrator.k_higher[i] = integrator.du[i] * integrator.dt
end
integrator.f(integrator.du, integrator.u_tmp, prob.p,
integrator.t +
alg.c[alg.num_stages - alg.num_stage_evals + j] * integrator.dt)

@threaded for i in eachindex(integrator.u)
integrator.u[i] += integrator.k_higher[i] *
alg.b[2]
@threaded for i in eachindex(integrator.du)
integrator.k_higher[i] = integrator.du[i] * integrator.dt
end

@threaded for i in eachindex(integrator.u)
integrator.u[i] += integrator.k_higher[i] *
alg.b[j+1]
end
end

# Higher stages after num_stage_evals where b is non-zero
Expand All @@ -427,7 +430,7 @@ 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]
alg.b[stage - alg.num_stages + alg.num_stage_evals+1]
end
end

Expand Down

0 comments on commit 2af55dd

Please sign in to comment.