diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl index 98b581ae45a..f9453d0c2c6 100644 --- a/src/solvers/fv_t8code/containers.jl +++ b/src/solvers/fv_t8code/containers.jl @@ -70,7 +70,8 @@ function Base.resize!(elements::T8codeFVElementContainer, capacity) resize!(elements.reconstruction_stencil, capacity) resize!(_reconstruction_gradient, n_dims * n_variables * capacity) - elements.reconstruction_gradient = unsafe_wrap(Array, pointer(_face_normals), (ndims, n_variables, capacity)) + elements.reconstruction_gradient = unsafe_wrap(Array, pointer(_face_normals), + (ndims, n_variables, capacity)) return nothing end @@ -109,15 +110,19 @@ function init_elements(mesh::T8codeMesh{NDIMS, RealT}, reconstruction_stencil = Vector{Vector{Int}}(undef, nelements) _reconstruction_gradient = Vector{RealT}(undef, NDIMS * n_variables * nelements) - reconstruction_gradient = unsafe_wrap(Array, pointer(_reconstruction_gradient), (NDIMS, n_variables, nelements)) + reconstruction_gradient = unsafe_wrap(Array, pointer(_reconstruction_gradient), + (NDIMS, n_variables, nelements)) elements = T8codeFVElementContainer{NDIMS, RealT, uEltype}(level, volume, midpoint, dx, num_faces, face_midpoints, face_areas, face_normals, - reconstruction_stencil, reconstruction_gradient, - _midpoint, _face_midpoints, - _face_areas, _face_normals, + reconstruction_stencil, + reconstruction_gradient, + _midpoint, + _face_midpoints, + _face_areas, + _face_normals, _reconstruction_gradient) init_elements!(elements, mesh, solver) @@ -219,9 +224,12 @@ end corner_coords = view(corners, :, corner, element) # loop over all corners of `possible_stencil_neighbor` for possible_corner in 1:num_faces[possible_stencil_neighbor] - if corner_coords == view(corners, :, possible_corner, possible_stencil_neighbor) - append!(reconstruction_stencil[element], possible_stencil_neighbor) - append!(reconstruction_stencil[possible_stencil_neighbor], element) + if corner_coords == + view(corners, :, possible_corner, possible_stencil_neighbor) + append!(reconstruction_stencil[element], + possible_stencil_neighbor) + append!(reconstruction_stencil[possible_stencil_neighbor], + element) end end end diff --git a/src/solvers/fv_t8code/fv.jl b/src/solvers/fv_t8code/fv.jl index 2795fd5b329..3e97a713caa 100644 --- a/src/solvers/fv_t8code/fv.jl +++ b/src/solvers/fv_t8code/fv.jl @@ -159,7 +159,11 @@ function rhs!(du, u, t, mesh::T8codeMesh, equations, @trixi_timeit timer() "exchange_solution!" exchange_solution!(u, mesh, equations, solver, cache) - @trixi_timeit timer() "calc gradient reconstruction" calc_gradient_reconstruction!(u, mesh, equations, solver, cache) + @trixi_timeit timer() "gradient reconstruction" calc_gradient_reconstruction!(u, + mesh, + equations, + solver, + cache) # Prolong solution to interfaces @trixi_timeit timer() "prolong2interfaces" begin @@ -221,7 +225,8 @@ function calc_gradient_reconstruction!(u, mesh, equations, solver, cache) for element in eachelement(mesh, solver, cache) n_stencil_neighbors = length(reconstruction_stencil[element]) - coordinates_element = get_node_coords(elements.midpoint, equations, solver, element) + coordinates_element = get_node_coords(elements.midpoint, equations, solver, + element) # Reset variables a = zero(eltype(u)) @@ -233,15 +238,19 @@ function calc_gradient_reconstruction!(u, mesh, equations, solver, cache) e .= zero(eltype(u)) for i in 1:n_stencil_neighbors neighbor = reconstruction_stencil[element][i] - coordinates_neighbor = get_node_coords(elements.midpoint, equations, solver, neighbor) + coordinates_neighbor = get_node_coords(elements.midpoint, equations, solver, + neighbor) a += (coordinates_neighbor[1] - coordinates_element[1])^2 - b += (coordinates_neighbor[1] - coordinates_element[1]) * (coordinates_neighbor[2] - coordinates_element[2]) + b += (coordinates_neighbor[1] - coordinates_element[1]) * + (coordinates_neighbor[2] - coordinates_element[2]) c += (coordinates_neighbor[2] - coordinates_element[2])^2 for v in eachvariable(equations) - d[v] += (coordinates_neighbor[1] - coordinates_element[1]) * (u[v, neighbor] - u[v, element]) - e[v] += (coordinates_neighbor[2] - coordinates_element[2]) * (u[v, neighbor] - u[v, element]) + d[v] += (coordinates_neighbor[1] - coordinates_element[1]) * + (u[v, neighbor] - u[v, element]) + e[v] += (coordinates_neighbor[2] - coordinates_element[2]) * + (u[v, neighbor] - u[v, element]) end end @@ -263,7 +272,7 @@ end function prolong2interfaces!(cache, mesh::T8codeMesh, equations, solver::FV) (; elements, interfaces, u_tmp) = cache - (; reconstruction_gradient) = elements + (; midpoint, face_midpoints, reconstruction_gradient) = elements for interface in eachinterface(solver, cache) element = interfaces.neighbor_ids[1, interface] @@ -277,21 +286,26 @@ function prolong2interfaces!(cache, mesh::T8codeMesh, equations, solver::FV) face_element = interfaces.faces[1, interface] face_neighbor = interfaces.faces[2, interface] - face_midpoint_element = get_node_coords(elements.face_midpoints, equations, solver, face_element, element) - face_midpoint_neighbor = get_node_coords(elements.face_midpoints, equations, solver, face_neighbor, neighbor) + face_midpoint_element = get_node_coords(face_midpoints, equations, solver, + face_element, element) + face_midpoint_neighbor = get_node_coords(face_midpoints, equations, solver, + face_neighbor, neighbor) - midpoint_element = get_node_coords(elements.midpoint, equations, solver, element) - midpoint_neighbor = get_node_coords(elements.midpoint, equations, solver, neighbor) + midpoint_element = get_node_coords(midpoint, equations, solver, element) + midpoint_neighbor = get_node_coords(midpoint, equations, solver, neighbor) vector_element = face_midpoint_element .- midpoint_element vector_neighbor = face_midpoint_neighbor .- midpoint_neighbor for v in eachvariable(equations) - gradient_v_element = get_node_coords(reconstruction_gradient, equations, solver, v, element) - gradient_v_neighbor = get_node_coords(reconstruction_gradient, equations, solver, v, neighbor) + gradient_v_element = get_node_coords(reconstruction_gradient, equations, + solver, v, element) + gradient_v_neighbor = get_node_coords(reconstruction_gradient, + equations, solver, v, neighbor) interfaces.u[1, v, interface] = u_tmp[element].u[v] + dot(gradient_v_element, vector_element) interfaces.u[2, v, interface] = u_tmp[neighbor].u[v] + - dot(gradient_v_neighbor, vector_neighbor) + dot(gradient_v_neighbor, + vector_neighbor) end else error("Order $(solver.order) is not supported.")