partitioned rhs structured & tree

DanielDoehring · Mar 21, 2024 · 20acb37 · 20acb37
1 parent a424eb1
commit 20acb37
Show file tree

Hide file tree

Showing 3 changed files with 119 additions and 6 deletions.
diff --git a/src/solvers/dgsem_structured/dg_1d.jl b/src/solvers/dgsem_structured/dg_1d.jl
@@ -9,7 +9,6 @@ function rhs!(du, u, t,
               mesh::StructuredMesh{1}, equations,
               initial_condition, boundary_conditions, source_terms::Source,
               dg::DG, cache) where {Source}
-
     element_range = eachelement(dg, cache)
 
     # Reset du
@@ -52,6 +51,55 @@ function rhs!(du, u, t,
     return nothing
 end
 
+# RHS for PERK integrator
+function rhs!(du, u, t,
+              mesh::StructuredMesh{1}, equations,
+              initial_condition, boundary_conditions, source_terms::Source,
+              dg::DG, cache,
+              element_range,
+              # Interfaces, boundaries, boundary_orientations, mortars not present for structured mesh
+              _, _, _, _) where {Source}
+
+    # Reset du
+    @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache, element_range)
+
+    # Calculate volume integral
+    @trixi_timeit timer() "volume integral" begin
+        calc_volume_integral!(du, u, mesh,
+                              have_nonconservative_terms(equations), equations,
+                              dg.volume_integral, dg, cache, element_range)
+    end
+
+    # Calculate interface and boundary fluxes
+    @trixi_timeit timer() "interface flux" begin
+        calc_interface_flux!(cache, u, mesh, equations, dg.surface_integral, dg,
+                             element_range)
+    end
+
+    # Calculate boundary fluxes
+    @trixi_timeit timer() "boundary flux" begin
+        calc_boundary_flux!(cache, u, t, boundary_conditions, mesh, equations,
+                            dg.surface_integral, dg)
+    end
+
+    # Calculate surface integrals
+    @trixi_timeit timer() "surface integral" begin
+        calc_surface_integral!(du, u, mesh, equations,
+                               dg.surface_integral, dg, cache, element_range)
+    end
+
+    # Apply Jacobian from mapping to reference element
+    @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache,
+                                                     element_range)
+
+    # Calculate source terms
+    @trixi_timeit timer() "source terms" begin
+        calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range)
+    end
+
+    return nothing
+end
+
 function calc_interface_flux!(cache, u, mesh::StructuredMesh{1},
                               equations, surface_integral, dg::DG,
                               element_range)

diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
@@ -78,7 +78,6 @@ function rhs!(du, u, t,
               mesh::TreeMesh{1}, equations,
               initial_condition, boundary_conditions, source_terms::Source,
               dg::DG, cache) where {Source}
-
     element_range = eachelement(dg, cache)
     interface_range = eachinterface(dg, cache)
     boundary_range = eachboundary(dg, cache)
@@ -136,6 +135,69 @@ function rhs!(du, u, t,
     return nothing
 end
 
+# RHS for PERK integrator
+function rhs!(du, u, t,
+              mesh::TreeMesh{1}, equations,
+              initial_condition, boundary_conditions, source_terms::Source,
+              dg::DG, cache,
+              element_range, interface_range, boundary_range,
+              boundary_orientation_range,
+              # Mortars not appearant in 1D
+              _) where {Source}
+
+    # Reset du
+    @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, element_range)
+
+    # Calculate volume integral
+    @trixi_timeit timer() "volume integral" begin
+        calc_volume_integral!(du, u, mesh,
+                              have_nonconservative_terms(equations), equations,
+                              dg.volume_integral, dg, cache, element_range)
+    end
+
+    # Prolong solution to interfaces, i.e., reconstruct interface/trace values
+    @trixi_timeit timer() "prolong2interfaces" begin
+        prolong2interfaces!(cache, u, mesh, equations,
+                            dg.surface_integral, dg, interface_range)
+    end
+
+    # Calculate interface fluxes
+    @trixi_timeit timer() "interface flux" begin
+        calc_interface_flux!(cache.elements.surface_flux_values, mesh,
+                             have_nonconservative_terms(equations), equations,
+                             dg.surface_integral, dg, cache, interface_range)
+    end
+
+    # Prolong solution to boundaries
+    @trixi_timeit timer() "prolong2boundaries" begin
+        prolong2boundaries!(cache, u, mesh, equations,
+                            dg.surface_integral, dg, boundary_range)
+    end
+
+    # Calculate boundary fluxes
+    @trixi_timeit timer() "boundary flux" begin
+        calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations,
+                            dg.surface_integral, dg, boundary_orientation_range)
+    end
+
+    # Calculate surface integrals
+    @trixi_timeit timer() "surface integral" begin
+        calc_surface_integral!(du, u, mesh, equations,
+                               dg.surface_integral, dg, cache, element_range)
+    end
+
+    # Apply Jacobian from mapping to reference element
+    @trixi_timeit timer() "Jacobian" apply_jacobian!(du, mesh, equations, dg, cache,
+                                                     element_range)
+
+    # Calculate source terms
+    @trixi_timeit timer() "source terms" begin
+        calc_sources!(du, u, t, source_terms, equations, dg, cache, element_range)
+    end
+
+    return nothing
+end
+
 function calc_volume_integral!(du, u,
                                mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                                nonconservative_terms, equations,
@@ -613,7 +675,8 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
 end
 
 function calc_surface_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}},
-                                equations, surface_integral, dg::DGSEM, cache, element_range)
+                                equations, surface_integral, dg::DGSEM, cache,
+                                element_range)
     @unpack boundary_interpolation = dg.basis
     @unpack surface_flux_values = cache.elements
 

diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl
@@ -64,7 +64,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1},
 
     # Calculate volume integral
     @trixi_timeit timer() "volume integral" begin
-        calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache, element_range)
+        calc_volume_integral!(du, flux_viscous, mesh, equations_parabolic, dg, cache,
+                              element_range)
     end
 
     # Prolong solution to interfaces
@@ -101,7 +102,8 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1},
 
     # Apply Jacobian from mapping to reference element
     @trixi_timeit timer() "Jacobian" begin
-        apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic, element_range)
+        apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic,
+                                  element_range)
     end
 
     return nothing
@@ -112,7 +114,7 @@ end
 # TODO: can we avoid copying data?
 function transform_variables!(u_transformed, u, mesh::TreeMesh{1},
                               equations_parabolic::AbstractEquationsParabolic,
-                              dg::DG, parabolic_scheme, cache, cache_parabolic, 
+                              dg::DG, parabolic_scheme, cache, cache_parabolic,
                               element_range)
     transformation = gradient_variable_transformation(equations_parabolic)