Skip to content

Commit

Permalink
Remove duplicate code (rhs_parabolic!) for 2D, 3D (#1810)
Browse files Browse the repository at this point in the history 
* Remove duplicate code

* comment
  • Loading branch information
@DanielDoehring
DanielDoehring authored Jan 22, 2024
1 parent 1946f9d commit 43dac2c
Show file tree
Hide file tree
Showing 4 changed files with 24 additions and 224 deletions.
26 changes: 23 additions & 3 deletions src/solvers/dgsem_p4est/dg_2d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,29 @@ function create_cache_parabolic(mesh::P4estMesh{2}, equations_hyperbolic::Abstra
return cache
end

# TODO: Remove in favor of the implementation for the TreeMesh
# once the P4estMesh can handle mortars as well
function rhs_parabolic!(du, u, t, mesh::P4estMesh{2},
#=
Reusing `rhs_parabolic!` for `TreeMesh`es is not easily possible as
for `P4estMesh`es we call
```
prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic,
dg.mortar, dg.surface_integral, dg)
calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values,
mesh, equations_parabolic, dg.mortar,
dg.surface_integral, dg, cache)
```
instead of
```
prolong2mortars!(cache, flux_viscous, mesh, equations_parabolic,
dg.mortar, dg.surface_integral, dg)
calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, mesh,
equations_parabolic,
dg.mortar, dg.surface_integral, dg, cache)
```
=#
function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{2}, P4estMesh{3}},
equations_parabolic::AbstractEquationsParabolic,
initial_condition, boundary_conditions_parabolic, source_terms,
dg::DG, parabolic_scheme, cache, cache_parabolic)
Expand Down
112 changes: 0 additions & 112 deletions src/solvers/dgsem_p4est/dg_3d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,118 +19,6 @@ function create_cache_parabolic(mesh::P4estMesh{3}, equations_hyperbolic::Abstra
return cache
end

# This file collects all methods that have been updated to work with parabolic systems of equations
#
# assumptions: parabolic terms are of the form div(f(u, grad(u))) and
# will be discretized first order form as follows:
# 1. compute grad(u)
# 2. compute f(u, grad(u))
# 3. compute div(f(u, grad(u))) (i.e., the "regular" rhs! call)
# boundary conditions will be applied to both grad(u) and div(f(u, grad(u))).
# TODO: Remove in favor of the implementation for the TreeMesh
# once the P4estMesh can handle mortars as well
function rhs_parabolic!(du, u, t, mesh::P4estMesh{3},
equations_parabolic::AbstractEquationsParabolic,
initial_condition, boundary_conditions_parabolic, source_terms,
dg::DG, parabolic_scheme, cache, cache_parabolic)
@unpack viscous_container = cache_parabolic
@unpack u_transformed, gradients, flux_viscous = viscous_container

# Convert conservative variables to a form more suitable for viscous flux calculations
@trixi_timeit timer() "transform variables" begin
transform_variables!(u_transformed, u, mesh, equations_parabolic,
dg, parabolic_scheme, cache, cache_parabolic)
end

# Compute the gradients of the transformed variables
@trixi_timeit timer() "calculate gradient" begin
calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic,
boundary_conditions_parabolic, dg, cache, cache_parabolic)
end

# Compute and store the viscous fluxes
@trixi_timeit timer() "calculate viscous fluxes" begin
calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh,
equations_parabolic, dg, cache, cache_parabolic)
end

# The remainder of this function is essentially a regular rhs! for parabolic
# equations (i.e., it computes the divergence of the viscous fluxes)
#
# OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have
# been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the
# `interfaces` and `boundaries` containers in `cache_parabolic` to interpolate and store the
# *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it
# is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values*
# and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we
# do not need to recreate the existing data structure only with a different name, and c) we do not
# need to interpolate solutions *and* gradients to the surfaces.

# TODO: parabolic; reconsider current data structure reuse strategy

# Reset du
@trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache)

# Calculate volume integral
@trixi_timeit timer() "volume integral" begin
calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache)
end

# Prolong solution to interfaces
@trixi_timeit timer() "prolong2interfaces" begin
prolong2interfaces!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
dg.surface_integral, dg, cache)
end

# Calculate interface fluxes
@trixi_timeit timer() "interface flux" begin
calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh,
equations_parabolic, dg, cache_parabolic)
end

# Prolong solution to boundaries
@trixi_timeit timer() "prolong2boundaries" begin
prolong2boundaries!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
dg.surface_integral, dg, cache)
end

# Calculate boundary fluxes
@trixi_timeit timer() "boundary flux" begin
calc_boundary_flux_divergence!(cache_parabolic, t,
boundary_conditions_parabolic,
mesh, equations_parabolic,
dg.surface_integral, dg)
end

# Prolong solution to mortars (specialized for AbstractEquationsParabolic)
# !!! NOTE: we reuse the hyperbolic cache here since it contains "mortars" and "u_threaded"
# !!! Is this OK?
@trixi_timeit timer() "prolong2mortars" begin
prolong2mortars_divergence!(cache, flux_viscous, mesh, equations_parabolic,
dg.mortar, dg.surface_integral, dg)
end

# Calculate mortar fluxes (specialized for AbstractEquationsParabolic)
@trixi_timeit timer() "mortar flux" begin
calc_mortar_flux_divergence!(cache_parabolic.elements.surface_flux_values,
mesh, equations_parabolic, dg.mortar,
dg.surface_integral, dg, cache)
end

# Calculate surface integrals
@trixi_timeit timer() "surface integral" begin
calc_surface_integral!(du, u, mesh, equations_parabolic,
dg.surface_integral, dg, cache_parabolic)
end

# Apply Jacobian from mapping to reference element
@trixi_timeit timer() "Jacobian" begin
apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic)
end

return nothing
end

function calc_gradient!(gradients, u_transformed, t,
mesh::P4estMesh{3}, equations_parabolic,
boundary_conditions_parabolic, dg::DG,
Expand Down
2 changes: 1 addition & 1 deletion src/solvers/dgsem_tree/dg_2d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
# 2. compute f(u, grad(u))
# 3. compute div(f(u, grad(u))) (i.e., the "regular" rhs! call)
# boundary conditions will be applied to both grad(u) and div(f(u, grad(u))).
function rhs_parabolic!(du, u, t, mesh::TreeMesh{2},
function rhs_parabolic!(du, u, t, mesh::Union{TreeMesh{2}, TreeMesh{3}},
equations_parabolic::AbstractEquationsParabolic,
initial_condition, boundary_conditions_parabolic, source_terms,
dg::DG, parabolic_scheme, cache, cache_parabolic)
Expand Down
108 changes: 0 additions & 108 deletions src/solvers/dgsem_tree/dg_3d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,114 +5,6 @@
@muladd begin
#! format: noindent

# This file collects all methods that have been updated to work with parabolic systems of equations
#
# assumptions: parabolic terms are of the form div(f(u, grad(u))) and
# will be discretized first order form as follows:
# 1. compute grad(u)
# 2. compute f(u, grad(u))
# 3. compute div(f(u, grad(u))) (i.e., the "regular" rhs! call)
# boundary conditions will be applied to both grad(u) and div(f(u, grad(u))).
function rhs_parabolic!(du, u, t, mesh::TreeMesh{3},
equations_parabolic::AbstractEquationsParabolic,
initial_condition, boundary_conditions_parabolic, source_terms,
dg::DG, parabolic_scheme, cache, cache_parabolic)
@unpack viscous_container = cache_parabolic
@unpack u_transformed, gradients, flux_viscous = viscous_container

# Convert conservative variables to a form more suitable for viscous flux calculations
@trixi_timeit timer() "transform variables" begin
transform_variables!(u_transformed, u, mesh, equations_parabolic,
dg, parabolic_scheme, cache, cache_parabolic)
end

# Compute the gradients of the transformed variables
@trixi_timeit timer() "calculate gradient" begin
calc_gradient!(gradients, u_transformed, t, mesh, equations_parabolic,
boundary_conditions_parabolic, dg, cache, cache_parabolic)
end

# Compute and store the viscous fluxes
@trixi_timeit timer() "calculate viscous fluxes" begin
calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh,
equations_parabolic, dg, cache, cache_parabolic)
end

# The remainder of this function is essentially a regular rhs! for parabolic
# equations (i.e., it computes the divergence of the viscous fluxes)
#
# OBS! In `calc_viscous_fluxes!`, the viscous flux values at the volume nodes of each element have
# been computed and stored in `fluxes_viscous`. In the following, we *reuse* (abuse) the
# `interfaces` and `boundaries` containers in `cache_parabolic` to interpolate and store the
# *fluxes* at the element surfaces, as opposed to interpolating and storing the *solution* (as it
# is done in the hyperbolic operator). That is, `interfaces.u`/`boundaries.u` store *viscous flux values*
# and *not the solution*. The advantage is that a) we do not need to allocate more storage, b) we
# do not need to recreate the existing data structure only with a different name, and c) we do not
# need to interpolate solutions *and* gradients to the surfaces.

# TODO: parabolic; reconsider current data structure reuse strategy

# Reset du
@trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache)

# Calculate volume integral
@trixi_timeit timer() "volume integral" begin
calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache)
end

# Prolong solution to interfaces
@trixi_timeit timer() "prolong2interfaces" begin
prolong2interfaces!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
dg.surface_integral, dg, cache)
end

# Calculate interface fluxes
@trixi_timeit timer() "interface flux" begin
calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh,
equations_parabolic, dg, cache_parabolic)
end

# Prolong solution to boundaries
@trixi_timeit timer() "prolong2boundaries" begin
prolong2boundaries!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
dg.surface_integral, dg, cache)
end

# Calculate boundary fluxes
@trixi_timeit timer() "boundary flux" begin
calc_boundary_flux_divergence!(cache_parabolic, t,
boundary_conditions_parabolic,
mesh, equations_parabolic,
dg.surface_integral, dg)
end

# Prolong solution to mortars
@trixi_timeit timer() "prolong2mortars" begin
prolong2mortars!(cache, flux_viscous, mesh, equations_parabolic,
dg.mortar, dg.surface_integral, dg)
end

# Calculate mortar fluxes
@trixi_timeit timer() "mortar flux" begin
calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, mesh,
equations_parabolic,
dg.mortar, dg.surface_integral, dg, cache)
end

# Calculate surface integrals
@trixi_timeit timer() "surface integral" begin
calc_surface_integral!(du, u, mesh, equations_parabolic,
dg.surface_integral, dg, cache_parabolic)
end

# Apply Jacobian from mapping to reference element
@trixi_timeit timer() "Jacobian" begin
apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic)
end

return nothing
end

# Transform solution variables prior to taking the gradient
# (e.g., conservative to primitive variables). Defaults to doing nothing.
# TODO: can we avoid copying data?
Expand Down

0 comments on commit 43dac2c

Please sign in to comment.