diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index a7f3345168f..299f2f6140a 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -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) diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index 0bb97c7af02..83d663809a7 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -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, diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 3083ae30680..b1c27343999 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -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) diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 9ad28c6aa8e..ee0e7c6b069 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -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?