From 3256ecc87a13ea9c2713dfae049db024001ce680 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 8 Oct 2024 11:55:24 +0200 Subject: [PATCH 1/3] Thread-parallelize src term computation Euler Gravity --- .../semidiscretization_euler_gravity.jl | 41 ++++++++++++------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index 3ca429f63b6..f7a48b66e88 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -232,22 +232,29 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t) # compute gravitational potential and forces @trixi_timeit timer() "gravity solver" update_gravity!(semi, u_ode) + n_elements = size(u_euler)[end] # add gravitational source source_terms to the Euler part if ndims(semi_euler) == 1 - @views @. du_euler[2, .., :] -= u_euler[1, .., :] * u_gravity[2, .., :] - @views @. du_euler[3, .., :] -= u_euler[2, .., :] * u_gravity[2, .., :] + @threaded for i in 1:n_elements + @views @. du_euler[2, .., i] -= u_euler[1, .., i] * u_gravity[2, .., i] + @views @. du_euler[3, .., i] -= u_euler[2, .., i] * u_gravity[2, .., i] + end elseif ndims(semi_euler) == 2 - @views @. du_euler[2, .., :] -= u_euler[1, .., :] * u_gravity[2, .., :] - @views @. du_euler[3, .., :] -= u_euler[1, .., :] * u_gravity[3, .., :] - @views @. du_euler[4, .., :] -= (u_euler[2, .., :] * u_gravity[2, .., :] + - u_euler[3, .., :] * u_gravity[3, .., :]) + @threaded for i in 1:n_elements + @views @. du_euler[2, .., i] -= u_euler[1, .., i] * u_gravity[2, .., i] + @views @. du_euler[3, .., i] -= u_euler[1, .., i] * u_gravity[3, .., i] + @views @. du_euler[4, .., i] -= (u_euler[2, .., i] * u_gravity[2, .., i] + + u_euler[3, .., i] * u_gravity[3, .., i]) + end elseif ndims(semi_euler) == 3 - @views @. du_euler[2, .., :] -= u_euler[1, .., :] * u_gravity[2, .., :] - @views @. du_euler[3, .., :] -= u_euler[1, .., :] * u_gravity[3, .., :] - @views @. du_euler[4, .., :] -= u_euler[1, .., :] * u_gravity[4, .., :] - @views @. du_euler[5, .., :] -= (u_euler[2, .., :] * u_gravity[2, .., :] + - u_euler[3, .., :] * u_gravity[3, .., :] + - u_euler[4, .., :] * u_gravity[4, .., :]) + @threaded for i in 1:n_elements + @views @. du_euler[2, .., i] -= u_euler[1, .., i] * u_gravity[2, .., i] + @views @. du_euler[3, .., i] -= u_euler[1, .., i] * u_gravity[3, .., i] + @views @. du_euler[4, .., i] -= u_euler[1, .., i] * u_gravity[4, .., i] + @views @. du_euler[5, .., i] -= (u_euler[2, .., i] * u_gravity[2, .., i] + + u_euler[3, .., i] * u_gravity[3, .., i] + + u_euler[4, .., i] * u_gravity[4, .., i]) + end else error("Number of dimensions $(ndims(semi_euler)) not supported.") end @@ -335,7 +342,10 @@ function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gr # Source term: Jeans instability OR coupling convergence test OR blast wave # put in gravity source term proportional to Euler density # OBS! subtract off the background density ρ_0 (spatial mean value) - @views @. du_gravity[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0) + n_elements = size(u_euler)[end] + @threaded for i in 1:n_elements + @views @. du_gravity[1, .., i] += grav_scale * (u_euler[1, .., i] - rho0) + end a_stage = a[stage] b_stage_dt = b[stage] * dt @@ -390,7 +400,10 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, # Source term: Jeans instability OR coupling convergence test OR blast wave # put in gravity source term proportional to Euler density # OBS! subtract off the background density ρ_0 around which the Jeans instability is perturbed - @views @. du_gravity[1, .., :] += grav_scale * (u_euler[1, .., :] - rho0) + n_elements = size(u_euler)[end] + @threaded for i in 1:n_elements + @views @. du_gravity[1, .., i] += grav_scale * (u_euler[1, .., i] - rho0) + end delta_stage = delta[stage] gamma1_stage = gamma1[stage] From 9b22f46308f00b9c83ea4d4ae4444cf20af0b96e Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Wed, 9 Oct 2024 10:42:19 +0200 Subject: [PATCH 2/3] element --- .../semidiscretization_euler_gravity.jl | 52 ++++++++++++------- 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index f7a48b66e88..71e3e47fd2a 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -235,25 +235,37 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t) n_elements = size(u_euler)[end] # add gravitational source source_terms to the Euler part if ndims(semi_euler) == 1 - @threaded for i in 1:n_elements - @views @. du_euler[2, .., i] -= u_euler[1, .., i] * u_gravity[2, .., i] - @views @. du_euler[3, .., i] -= u_euler[2, .., i] * u_gravity[2, .., i] + @threaded for element in 1:n_elements + @views @. du_euler[2, .., element] -= u_euler[1, .., element] * + u_gravity[2, .., element] + @views @. du_euler[3, .., element] -= u_euler[2, .., element] * + u_gravity[2, .., element] end elseif ndims(semi_euler) == 2 - @threaded for i in 1:n_elements - @views @. du_euler[2, .., i] -= u_euler[1, .., i] * u_gravity[2, .., i] - @views @. du_euler[3, .., i] -= u_euler[1, .., i] * u_gravity[3, .., i] - @views @. du_euler[4, .., i] -= (u_euler[2, .., i] * u_gravity[2, .., i] + - u_euler[3, .., i] * u_gravity[3, .., i]) + @threaded for element in 1:n_elements + @views @. du_euler[2, .., element] -= u_euler[1, .., element] * + u_gravity[2, .., element] + @views @. du_euler[3, .., element] -= u_euler[1, .., element] * + u_gravity[3, .., element] + @views @. du_euler[4, .., element] -= (u_euler[2, .., element] * + u_gravity[2, .., element] + + u_euler[3, .., element] * + u_gravity[3, .., element]) end elseif ndims(semi_euler) == 3 - @threaded for i in 1:n_elements - @views @. du_euler[2, .., i] -= u_euler[1, .., i] * u_gravity[2, .., i] - @views @. du_euler[3, .., i] -= u_euler[1, .., i] * u_gravity[3, .., i] - @views @. du_euler[4, .., i] -= u_euler[1, .., i] * u_gravity[4, .., i] - @views @. du_euler[5, .., i] -= (u_euler[2, .., i] * u_gravity[2, .., i] + - u_euler[3, .., i] * u_gravity[3, .., i] + - u_euler[4, .., i] * u_gravity[4, .., i]) + @threaded for element in 1:n_elements + @views @. du_euler[2, .., element] -= u_euler[1, .., element] * + u_gravity[2, .., element] + @views @. du_euler[3, .., element] -= u_euler[1, .., element] * + u_gravity[3, .., element] + @views @. du_euler[4, .., element] -= u_euler[1, .., element] * + u_gravity[4, .., element] + @views @. du_euler[5, .., element] -= (u_euler[2, .., element] * + u_gravity[2, .., element] + + u_euler[3, .., element] * + u_gravity[3, .., element] + + u_euler[4, .., element] * + u_gravity[4, .., element]) end else error("Number of dimensions $(ndims(semi_euler)) not supported.") @@ -343,8 +355,9 @@ function timestep_gravity_2N!(cache, u_euler, t, dt, gravity_parameters, semi_gr # put in gravity source term proportional to Euler density # OBS! subtract off the background density ρ_0 (spatial mean value) n_elements = size(u_euler)[end] - @threaded for i in 1:n_elements - @views @. du_gravity[1, .., i] += grav_scale * (u_euler[1, .., i] - rho0) + @threaded for element in 1:n_elements + @views @. du_gravity[1, .., element] += grav_scale * + (u_euler[1, .., element] - rho0) end a_stage = a[stage] @@ -401,8 +414,9 @@ function timestep_gravity_3Sstar!(cache, u_euler, t, dt, gravity_parameters, # put in gravity source term proportional to Euler density # OBS! subtract off the background density ρ_0 around which the Jeans instability is perturbed n_elements = size(u_euler)[end] - @threaded for i in 1:n_elements - @views @. du_gravity[1, .., i] += grav_scale * (u_euler[1, .., i] - rho0) + @threaded for element in 1:n_elements + @views @. du_gravity[1, .., element] += grav_scale * + (u_euler[1, .., element] - rho0) end delta_stage = delta[stage] From d1d307af55f83e11a803e5b6b04022af177a7722 Mon Sep 17 00:00:00 2001 From: Daniel_Doehring Date: Tue, 7 Jan 2025 11:47:51 +0100 Subject: [PATCH 3/3] revisit --- .../semidiscretization_euler_gravity.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl index 260fefccc20..3f8093822a5 100644 --- a/src/semidiscretization/semidiscretization_euler_gravity.jl +++ b/src/semidiscretization/semidiscretization_euler_gravity.jl @@ -227,6 +227,7 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t) u_euler = wrap_array(u_ode, semi_euler) du_euler = wrap_array(du_ode, semi_euler) u_gravity = wrap_array(cache.u_ode, semi_gravity) + n_elements = size(u_euler)[end] time_start = time_ns() @@ -236,7 +237,6 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerGravity, t) # compute gravitational potential and forces @trixi_timeit timer() "gravity solver" update_gravity!(semi, u_ode) - n_elements = size(u_euler)[end] # add gravitational source source_terms to the Euler part if ndims(semi_euler) == 1 @threaded for element in 1:n_elements @@ -347,8 +347,11 @@ function timestep_gravity_2N!(cache, u_euler, tau, dtau, gravity_parameters, # Note that `u_ode` is `u_gravity` in `rhs!` above @unpack u_ode, du_ode, u_tmp1_ode = cache + n_elements = size(u_euler)[end] + u_tmp1_ode .= zero(eltype(u_tmp1_ode)) du_gravity = wrap_array(du_ode, semi_gravity) + for stage in eachindex(c) tau_stage = tau + dtau * c[stage] @@ -360,7 +363,6 @@ function timestep_gravity_2N!(cache, u_euler, tau, dtau, gravity_parameters, # Source term: Jeans instability OR coupling convergence test OR blast wave # put in gravity source term proportional to Euler density # OBS! subtract off the background density ρ_0 (spatial mean value) - n_elements = size(u_euler)[end] # Note: Adding to `du_gravity` is essentially adding to `du_ode`! @threaded for element in 1:n_elements @views @. du_gravity[1, .., element] += grav_scale * @@ -407,9 +409,12 @@ function timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters, # Note that `u_ode` is `u_gravity` in `rhs!` above @unpack u_ode, du_ode, u_tmp1_ode, u_tmp2_ode = cache + n_elements = size(u_euler)[end] + u_tmp1_ode .= zero(eltype(u_tmp1_ode)) u_tmp2_ode .= u_ode du_gravity = wrap_array(du_ode, semi_gravity) + for stage in eachindex(c) tau_stage = tau + dtau * c[stage] @@ -421,7 +426,6 @@ function timestep_gravity_3Sstar!(cache, u_euler, tau, dtau, gravity_parameters, # Source term: Jeans instability OR coupling convergence test OR blast wave # put in gravity source term proportional to Euler density # OBS! subtract off the background density ρ_0 around which the Jeans instability is perturbed - n_elements = size(u_euler)[end] # Note: Adding to `du_gravity` is essentially adding to `du_ode`! @threaded for element in 1:n_elements @views @. du_gravity[1, .., element] += grav_scale *