diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl index 72dbe2c4256..7dfe4430244 100644 --- a/src/solvers/dgmulti/dg_parabolic.jl +++ b/src/solvers/dgmulti/dg_parabolic.jl @@ -62,9 +62,10 @@ end function transform_variables!(u_transformed, u, mesh, equations_parabolic::AbstractEquationsParabolic, dg::DGMulti, parabolic_scheme, cache, cache_parabolic) + transformation = gradient_variable_transformation(equations_parabolic) + @threaded for i in eachindex(u) - u_transformed[i] = gradient_variable_transformation(equations_parabolic)(u[i], - equations_parabolic) + u_transformed[i] = transformation(u[i], equations_parabolic) end end diff --git a/src/solvers/dgsem_tree/dg_1d_parabolic.jl b/src/solvers/dgsem_tree/dg_1d_parabolic.jl index c2aa75388c8..7602331d7c8 100644 --- a/src/solvers/dgsem_tree/dg_1d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_1d_parabolic.jl @@ -105,12 +105,13 @@ end function transform_variables!(u_transformed, u, mesh::TreeMesh{1}, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) + transformation = gradient_variable_transformation(equations_parabolic) + @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for i in eachnode(dg) u_node = get_node_vars(u, equations_parabolic, dg, i, element) - u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node, - equations_parabolic) + u_transformed_node = transformation(u_node, equations_parabolic) set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg, i, element) end @@ -147,16 +148,18 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) @unpack interfaces = cache_parabolic + @unpack neighbor_ids = interfaces + interfaces_u = interfaces.u @threaded for interface in eachinterface(dg, cache) - left_element = interfaces.neighbor_ids[1, interface] - right_element = interfaces.neighbor_ids[2, interface] + left_element = neighbor_ids[1, interface] + right_element = neighbor_ids[2, interface] # interface in x-direction for v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - interfaces.u[1, v, interface] = flux_viscous[v, nnodes(dg), left_element] - interfaces.u[2, v, interface] = flux_viscous[v, 1, right_element] + # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! + interfaces_u[1, v, interface] = flux_viscous[v, nnodes(dg), left_element] + interfaces_u[2, v, interface] = flux_viscous[v, 1, right_element] end end @@ -204,21 +207,22 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) @unpack boundaries = cache_parabolic - @unpack neighbor_sides = boundaries + @unpack neighbor_sides, neighbor_ids = boundaries + boundaries_u = boundaries.u @threaded for boundary in eachboundary(dg, cache_parabolic) - element = boundaries.neighbor_ids[boundary] + element = neighbor_ids[boundary] if neighbor_sides[boundary] == 1 # element in -x direction of boundary for v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[1, v, boundary] = flux_viscous[v, nnodes(dg), element] + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[1, v, boundary] = flux_viscous[v, nnodes(dg), element] end else # Element in +x direction of boundary for v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[2, v, boundary] = flux_viscous[v, 1, element] + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[2, v, boundary] = flux_viscous[v, 1, element] end end end @@ -552,8 +556,10 @@ end # where f(u) is the inviscid flux and g(u) is the viscous flux. function apply_jacobian_parabolic!(du, mesh::TreeMesh{1}, equations::AbstractEquationsParabolic, dg::DG, cache) + @unpack inverse_jacobian = cache.elements + @threaded for element in eachelement(dg, cache) - factor = cache.elements.inverse_jacobian[element] + factor = inverse_jacobian[element] for i in eachnode(dg) for v in eachvariable(equations) diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 0da25230380..3dbc55412ad 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -118,12 +118,13 @@ end function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) + transformation = gradient_variable_transformation(equations_parabolic) + @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations_parabolic, dg, i, j, element) - u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node, - equations_parabolic) + u_transformed_node = transformation(u_node, equations_parabolic) set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg, i, j, element) end @@ -168,30 +169,31 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) @unpack interfaces = cache_parabolic - @unpack orientations = interfaces + @unpack orientations, neighbor_ids = interfaces + interfaces_u = interfaces.u flux_viscous_x, flux_viscous_y = flux_viscous @threaded for interface in eachinterface(dg, cache) - left_element = interfaces.neighbor_ids[1, interface] - right_element = interfaces.neighbor_ids[2, interface] + left_element = neighbor_ids[1, interface] + right_element = neighbor_ids[2, interface] if orientations[interface] == 1 # interface in x-direction for j in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - interfaces.u[1, v, j, interface] = flux_viscous_x[v, nnodes(dg), j, + # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! + interfaces_u[1, v, j, interface] = flux_viscous_x[v, nnodes(dg), j, left_element] - interfaces.u[2, v, j, interface] = flux_viscous_x[v, 1, j, + interfaces_u[2, v, j, interface] = flux_viscous_x[v, 1, j, right_element] end else # if orientations[interface] == 2 # interface in y-direction for i in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - interfaces.u[1, v, i, interface] = flux_viscous_y[v, i, nnodes(dg), + # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! + interfaces_u[1, v, i, interface] = flux_viscous_y[v, i, nnodes(dg), left_element] - interfaces.u[2, v, i, interface] = flux_viscous_y[v, i, 1, + interfaces_u[2, v, i, interface] = flux_viscous_y[v, i, 1, right_element] end end @@ -244,25 +246,26 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) @unpack boundaries = cache_parabolic - @unpack orientations, neighbor_sides = boundaries + @unpack orientations, neighbor_sides, neighbor_ids = boundaries + boundaries_u = boundaries.u flux_viscous_x, flux_viscous_y = flux_viscous @threaded for boundary in eachboundary(dg, cache_parabolic) - element = boundaries.neighbor_ids[boundary] + element = neighbor_ids[boundary] if orientations[boundary] == 1 # boundary in x-direction if neighbor_sides[boundary] == 1 # element in -x direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[1, v, l, boundary] = flux_viscous_x[v, nnodes(dg), l, + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[1, v, l, boundary] = flux_viscous_x[v, nnodes(dg), l, element] end else # Element in +x direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[2, v, l, boundary] = flux_viscous_x[v, 1, l, element] + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[2, v, l, boundary] = flux_viscous_x[v, 1, l, element] end end else # if orientations[boundary] == 2 @@ -270,15 +273,15 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, if neighbor_sides[boundary] == 1 # element in -y direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[1, v, l, boundary] = flux_viscous_y[v, l, nnodes(dg), + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[1, v, l, boundary] = flux_viscous_y[v, l, nnodes(dg), element] end else # element in +y direction of boundary for l in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[2, v, l, boundary] = flux_viscous_y[v, l, 1, element] + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[2, v, l, boundary] = flux_viscous_y[v, l, 1, element] end end end @@ -608,7 +611,7 @@ function prolong2mortars!(cache, flux_viscous::Tuple{AbstractArray, AbstractArra end # NOTE: Use analogy to "calc_mortar_flux!" for hyperbolic eqs with no nonconservative terms. -# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for +# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for # hyperbolic terms with conserved terms only, i.e., no nonconservative terms. function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{2}, @@ -934,8 +937,10 @@ end # where f(u) is the inviscid flux and g(u) is the viscous flux. function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}}, equations::AbstractEquationsParabolic, dg::DG, cache) + @unpack inverse_jacobian = cache.elements + @threaded for element in eachelement(dg, cache) - factor = cache.elements.inverse_jacobian[element] + factor = inverse_jacobian[element] for j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations) diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 2745d312b37..9817e0e5f0e 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -118,12 +118,13 @@ end function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) + transformation = gradient_variable_transformation(equations_parabolic) + @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) u_node = get_node_vars(u, equations_parabolic, dg, i, j, k, element) - u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node, - equations_parabolic) + u_transformed_node = transformation(u_node, equations_parabolic) set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg, i, j, k, element) end @@ -175,43 +176,44 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) @unpack interfaces = cache_parabolic - @unpack orientations = interfaces + @unpack orientations, neighbor_ids = interfaces + interfaces_u = interfaces.u flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @threaded for interface in eachinterface(dg, cache) - left_element = interfaces.neighbor_ids[1, interface] - right_element = interfaces.neighbor_ids[2, interface] + left_element = neighbor_ids[1, interface] + right_element = neighbor_ids[2, interface] if orientations[interface] == 1 # interface in x-direction for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - interfaces.u[1, v, j, k, interface] = flux_viscous_x[v, nnodes(dg), j, + # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! + interfaces_u[1, v, j, k, interface] = flux_viscous_x[v, nnodes(dg), j, k, left_element] - interfaces.u[2, v, j, k, interface] = flux_viscous_x[v, 1, j, k, + interfaces_u[2, v, j, k, interface] = flux_viscous_x[v, 1, j, k, right_element] end elseif orientations[interface] == 2 # interface in y-direction for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - interfaces.u[1, v, i, k, interface] = flux_viscous_y[v, i, nnodes(dg), + # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! + interfaces_u[1, v, i, k, interface] = flux_viscous_y[v, i, nnodes(dg), k, left_element] - interfaces.u[2, v, i, k, interface] = flux_viscous_y[v, i, 1, k, + interfaces_u[2, v, i, k, interface] = flux_viscous_y[v, i, 1, k, right_element] end else # if orientations[interface] == 3 # interface in z-direction for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - interfaces.u[1, v, i, j, interface] = flux_viscous_z[v, i, j, + # OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*! + interfaces_u[1, v, i, j, interface] = flux_viscous_z[v, i, j, nnodes(dg), left_element] - interfaces.u[2, v, i, j, interface] = flux_viscous_z[v, i, j, 1, + interfaces_u[2, v, i, j, interface] = flux_viscous_z[v, i, j, 1, right_element] end end @@ -265,11 +267,12 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) @unpack boundaries = cache_parabolic - @unpack orientations, neighbor_sides = boundaries + @unpack orientations, neighbor_sides, neighbor_ids = boundaries + boundaries_u = boundaries.u flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @threaded for boundary in eachboundary(dg, cache_parabolic) - element = boundaries.neighbor_ids[boundary] + element = neighbor_ids[boundary] if orientations[boundary] == 1 # boundary in x-direction @@ -277,15 +280,15 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, # element in -x direction of boundary for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[1, v, j, k, boundary] = flux_viscous_x[v, nnodes(dg), + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[1, v, j, k, boundary] = flux_viscous_x[v, nnodes(dg), j, k, element] end else # Element in +x direction of boundary for k in eachnode(dg), j in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[2, v, j, k, boundary] = flux_viscous_x[v, 1, j, k, + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[2, v, j, k, boundary] = flux_viscous_x[v, 1, j, k, element] end end @@ -295,8 +298,8 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, # element in -y direction of boundary for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[1, v, i, k, boundary] = flux_viscous_y[v, i, + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[1, v, i, k, boundary] = flux_viscous_y[v, i, nnodes(dg), k, element] end @@ -304,8 +307,8 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, # element in +y direction of boundary for k in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[2, v, i, k, boundary] = flux_viscous_y[v, i, 1, k, + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[2, v, i, k, boundary] = flux_viscous_y[v, i, 1, k, element] end end @@ -315,8 +318,8 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, # element in -z direction of boundary for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[1, v, i, j, boundary] = flux_viscous_z[v, i, j, + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[1, v, i, j, boundary] = flux_viscous_z[v, i, j, nnodes(dg), element] end @@ -324,8 +327,8 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, # element in +z direction of boundary for j in eachnode(dg), i in eachnode(dg), v in eachvariable(equations_parabolic) - # OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*! - boundaries.u[2, v, i, j, boundary] = flux_viscous_z[v, i, j, 1, + # OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*! + boundaries_u[2, v, i, j, boundary] = flux_viscous_z[v, i, j, 1, element] end end @@ -820,7 +823,7 @@ function prolong2mortars!(cache, end # NOTE: Use analogy to "calc_mortar_flux!" for hyperbolic eqs with no nonconservative terms. -# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for +# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for # hyperbolic terms with conserved terms only, i.e., no nonconservative terms. function calc_mortar_flux!(surface_flux_values, mesh::TreeMesh{3}, @@ -1124,8 +1127,10 @@ end # where f(u) is the inviscid flux and g(u) is the viscous flux. function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations::AbstractEquationsParabolic, dg::DG, cache) + @unpack inverse_jacobian = cache.elements + @threaded for element in eachelement(dg, cache) - factor = cache.elements.inverse_jacobian[element] + factor = inverse_jacobian[element] for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) for v in eachvariable(equations)