Modularize calc_boundary_flux! for overset grid method (#71)

Author: JoshuaLampert

src/semidiscretization/semidiscretization_overset_grid.jl

@@ -138,49 +138,73 @@ function rhs!(du, u, t, mesh::OversetGridMesh, equations, initial_condition,
     end
 end
 
-function calc_boundary_flux!(surface_flux_values, u, t, boundary_conditions,
-                             mesh::OversetGridMesh, equations,
-                             integral::Tuple, solver, cache)
-    surface_flux_values_left, surface_flux_values_right = surface_flux_values
+function calc_boundary_flux_left!(surface_flux_values_left, u, t, x_neg,
+                                  equations, mesh, integral_left, solver,
+                                  cache)
     u_left, u_right = u
-    (; x_neg, x_pos) = boundary_conditions
-    mesh_left, mesh_right = mesh.mesh_left, mesh.mesh_right
-    integral_left, integral_right = integral
-    solver_left, solver_right = solver
-    cache_left, cache_right = cache
-    l_left, l_right = cache.l_left, cache.l_right
-    e_M_left, e_M_right = cache.e_M_left, cache.e_M_right
+    solver_left, _ = solver
+    mesh_left = mesh.mesh_left
+    cache_left, _ = cache
+    l_right = cache.l_right
+    e_M_right = cache.e_M_right
 
-    # Left boundary condition of left mesh
+    # Left boundary condition of left mesh (at a)
     e_left_L = get_projection_operator(cache_left.e_left, solver_left, 1)
     u_ll = x_neg(u, xmin(mesh), t, mesh, equations, solver, true, cache)
     u_rr = get_multiplied_node_vars(u_left, equations, e_left_L', :, 1)
     f = integral_left.surface_flux_boundary(u_ll, u_rr, equations)
     set_node_vars!(surface_flux_values_left, f, equations, 1, 1)
 
-    # Right boundary condition of left mesh
+    # Right boundary condition of left mesh (at c)
     e_right_L = get_projection_operator(cache_left.e_right, solver_left,
                                         nelements(mesh_left))
     u_ll = get_multiplied_node_vars(u_left, equations, e_right_L', :, nelements(mesh_left))
     u_rr = get_multiplied_node_vars(u_right, equations, e_M_right', :, l_right)
     f = integral_left.surface_flux_boundary(u_ll, u_rr, equations)
     set_node_vars!(surface_flux_values_left, f, equations, 2, nelements(mesh_left))
+    return nothing
+end
 
-    # Left boundary condition of right mesh
+function calc_boundary_flux_right!(surface_flux_values_right, u, t, x_pos,
+                                   equations, mesh, integral_right, solver,
+                                   cache)
+    u_left, u_right = u
+    _, solver_right = solver
+    mesh_right = mesh.mesh_right
+    _, cache_right = cache
+    l_left = cache.l_left
+    e_M_left = cache.e_M_left
+
+    # Left boundary condition of right mesh (at b)
     e_left_R = get_projection_operator(cache_right.e_left, solver_right, 1)
     u_ll = get_multiplied_node_vars(u_left, equations, e_M_left', :, l_left)
     u_rr = get_multiplied_node_vars(u_right, equations, e_left_R', :, 1)
     f = integral_right.surface_flux_boundary(u_ll, u_rr, equations)
     set_node_vars!(surface_flux_values_right, f, equations, 1, 1)
 
-    # Right boundary condition of right mesh
+    # Right boundary condition of right mesh (at d)
     e_right_R = get_projection_operator(cache_right.e_right, solver_right,
                                         nelements(mesh_right))
     u_ll = get_multiplied_node_vars(u_right, equations, e_right_R', :,
                                     nelements(mesh_right))
     u_rr = x_pos(u, xmax(mesh), t, mesh, equations, solver, false, cache)
     f = integral_right.surface_flux_boundary(u_ll, u_rr, equations)
-    return set_node_vars!(surface_flux_values_right, f, equations, 2, nelements(mesh_right))
+    set_node_vars!(surface_flux_values_right, f, equations, 2, nelements(mesh_right))
+    return nothing
+end
+
+function calc_boundary_flux!(surface_flux_values, u, t, boundary_conditions,
+                             mesh::OversetGridMesh, equations,
+                             integral::Tuple, solver, cache)
+    surface_flux_values_left, surface_flux_values_right = surface_flux_values
+    (; x_neg, x_pos) = boundary_conditions
+    integral_left, integral_right = integral
+
+    calc_boundary_flux_left!(surface_flux_values_left, u, t, x_neg, equations,
+                             mesh, integral_left, solver, cache)
+    calc_boundary_flux_right!(surface_flux_values_right, u, t, x_pos, equations,
+                              mesh, integral_right, solver, cache)
+    return nothing
 end
 
 # This method is for integrating a vector quantity for all variables over the entire domain,