From 1ee929709dbca08d3dbd458aca7f01c353751493 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 6 Aug 2024 16:53:48 +0200 Subject: [PATCH 01/12] Add first version of 3d support --- .../elixir_advection_basic_hybrid.jl | 43 +++++ src/meshes/t8code_mesh.jl | 51 ++++++ src/solvers/fv_t8code/containers.jl | 4 +- src/solvers/fv_t8code/fv.jl | 159 +++++++++++++++++- 4 files changed, 252 insertions(+), 5 deletions(-) create mode 100644 examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl diff --git a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl new file mode 100644 index 00000000000..67c664471d5 --- /dev/null +++ b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl @@ -0,0 +1,43 @@ +using OrdinaryDiffEq +using Trixi + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +initial_condition = initial_condition_convergence_test + +solver = FV(order = 2, extended_reconstruction_stencil = false, + surface_flux = flux_lax_friedrichs) + +# cmesh = Trixi.cmesh_new_periodic_hybrid() +cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true)) +# cmesh = Trixi.cmesh_new_periodic_tri() +mesh = T8codeMesh(cmesh, solver, + initial_refinement_level = 5) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +ode = semidiscretize(semi, (0.0, 10.0)); + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.3) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback, save_solution) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 8f6ec2a93da..0dafdf0af9e 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -1579,6 +1579,57 @@ function cmesh_quad(; comm = mpi_comm(), periodicity = (true, true))::t8_cmesh_t return cmesh end +function cmesh_quad_3d(; comm = mpi_comm(), periodicity = (true, true, true))::t8_cmesh_t + n_dims = 3 + vertices = [ # Just all vertices of all trees. partly duplicated + -1.0, -1.0, -1.0, # tree 0, quad + 1.0, -1.0, -1.0, + -1.0, 1.0, -1.0, + 1.0, 1.0, -1.0, + -1.0, -1.0, 1.0, + 1.0, -1.0, 1.0, + -1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, + ] + + # Generally, one can define other geometries. But besides linear the other + # geometries in t8code do not have C interface yet. + linear_geom = t8_geometry_linear_new(n_dims) + + # This is how the cmesh looks like. The numbers are the tree numbers: + # + # +---+ + # | | + # | 0 | + # | | + # +---+ + # + + cmesh_ref = Ref(t8_cmesh_t()) + t8_cmesh_init(cmesh_ref) + cmesh = cmesh_ref[] + + # Use linear geometry + t8_cmesh_register_geometry(cmesh, linear_geom) + t8_cmesh_set_tree_class(cmesh, 0, T8_ECLASS_HEX) + + t8_cmesh_set_tree_vertices(cmesh, 0, @views(vertices[(1 + 0):end]), 8) + + if periodicity[1] + t8_cmesh_set_join(cmesh, 0, 0, 0, 1, 0) + end + if periodicity[2] + t8_cmesh_set_join(cmesh, 0, 0, 2, 3, 0) + end + if periodicity[3] + t8_cmesh_set_join(cmesh, 0, 0, 4, 5, 0) + end + + t8_cmesh_commit(cmesh, comm) + + return cmesh +end + function cmesh_new_periodic_tri(; comm = mpi_comm())::t8_cmesh_t n_dims = 2 vertices = [ # Just all vertices of all trees. partly duplicated diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl index 5694e6e5dd2..44dcd19a5e5 100644 --- a/src/solvers/fv_t8code/containers.jl +++ b/src/solvers/fv_t8code/containers.jl @@ -248,7 +248,7 @@ function init_elements!(elements, mesh::T8codeMesh, solver::FV) return nothing end -@inline function init_extended_reconstruction_stencil!(corners, elements, solver) +@inline function init_extended_reconstruction_stencil!(corners, elements, solver::FV) if solver.order != 2 return nothing end @@ -327,7 +327,7 @@ end function init_reconstruction_stencil!(elements, interfaces, boundaries, communication_data, - mesh, equations, solver) + mesh, equations, solver::FV) if solver.order != 2 return nothing end diff --git a/src/solvers/fv_t8code/fv.jl b/src/solvers/fv_t8code/fv.jl index 08c82df056b..6445a2d7b23 100644 --- a/src/solvers/fv_t8code/fv.jl +++ b/src/solvers/fv_t8code/fv.jl @@ -216,7 +216,7 @@ function rhs!(du, u, t, mesh::T8codeMesh, equations, return nothing end -function calc_gradient_reconstruction!(u, mesh, equations, solver, cache) +function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver, cache) if solver.order == 1 return nothing elseif solver.order > 2 @@ -237,7 +237,7 @@ function calc_gradient_reconstruction!(u, mesh, equations, solver, cache) # (A^T A)^-1 = determinant_factor * [a3 -a2; -a2, a1] # A^T b = [d1; d2] - # Note: A^T b depends on the variable v. Using structure [d/e, v] + # Note: A^T b depends on the variable v. Using structure [1/2, v] a = zeros(eltype(u), 3) # [a1, a2, a3] d = zeros(eltype(u), size(u, 1), 2) # [d1(v), d2(v)] @@ -322,7 +322,7 @@ end @inline function calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d, reconstruction_stencil, reconstruction_distance, - solution_data, equations) + solution_data, equations::AbstractEquations{2}) for i in 1:n_stencil_neighbors # stencil=1 contains information from all neighbors # stencil=2,...,n_stencil_neighbors+1 is computed without (stencil+1)-th neighbor's information @@ -397,6 +397,159 @@ end return nothing end +function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver, cache) + if solver.order == 1 + return nothing + elseif solver.order > 2 + error("Order $(solver.order) is not supported yet!") + end + + (; elements, communication_data) = cache + (; reconstruction_stencil, reconstruction_distance, reconstruction_corner_elements, reconstruction_gradient, reconstruction_gradient_limited) = elements + (; solution_data) = communication_data + + # A N x 3 Matrix, where N is the number of stencil neighbors + # A^T A 3 x 3 Matrix + # b N Vector + # A^T b 3 Vector + + # Matrix/vector notation + # A^T A = [a11 a12 a13; a12 a22 a23; a13 a23 a33] + # det(A^T A) = a11 * (a22*a33 - a23^2) - a12 * (a12*a33 - a13*a23) + a13 * (a12*a23 - a22*a13) + # (A^T A)^-1 = 1/det(A^T A) * [(a33*a22-a23^2) -(a33*a12-a23*a13) (a23*a12-a22*a13) ; + # *** (a11*a33-a13^2) -(a23*a11-a12*a13); + # *** *** (a11*a22-a12^2) ] + + # A^T b = [d1; d2; d3] + # Note: A^T b depends on the variable v. Using structure [1/2/3, v] + a = zeros(eltype(u), 6) # [a11, a12, a13, a22, a23, a33] + d = zeros(eltype(u), size(u, 1), 3) # [d1(v), d2(v), d3(v)] + + # Parameter for limiting weights + lambda = [0.0, 1.0] + # lambda = [1.0, 0.0] # No limiting + r = 4 + epsilon = 1.0e-13 + + for element in eachelement(mesh, solver, cache) + if solver.extended_reconstruction_stencil + n_stencil_neighbors = elements.num_faces[element] + error("not yet supported!") + else + n_stencil_neighbors = length(reconstruction_stencil[element]) + end + + for stencil in 1:(n_stencil_neighbors + 1) + # Reset variables + a .= zero(eltype(u)) + # A^T b = [d1(v), d2(v), d3(v)] + # Note: A^T b depends on the variable v. Using vectors with d[v, 1/2/3] + d .= zero(eltype(u)) + + if solver.extended_reconstruction_stencil + calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d, + reconstruction_stencil, reconstruction_distance, + reconstruction_corner_elements, solution_data, + equations, solver, cache) + error("Not supported yet") + else + calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d, + reconstruction_stencil, reconstruction_distance, + solution_data, equations) + end + + # Divide by determinant + # Determinant = a11 * (a22 * a33 - a23^2 ) - a12 * (a12 * a33 - a13 * a23 ) + a13 * (a12 * a23 - a22 * a13 ) + AT_A_determinant = a[1] * (a[4] * a[6] - a[5]^2) - a[2] * (a[2] * a[6] - a[3] * a[5]) + a[3] * (a[2] * a[5] - a[4] * a[3]) + if isapprox(AT_A_determinant, 0.0) + for v in eachvariable(equations) + reconstruction_gradient[:, v, stencil, element] .= 0.0 + end + continue + end + + # det(A^T A) = a11 * (a22*a33 - a23^2) - a12 * (a12*a33 - a13*a23) + a13 * (a12*a23 - a22*a13) + # (A^T A)^-1 = 1/det(A^T A) * [(a33*a22-a23^2) -(a33*a12-a23*a13) (a23*a12-a22*a13) ; + # *** (a11*a33-a13^2) -(a23*a11-a12*a13); + # *** *** (a11*a22-a12^2) ] + # = [m11 m12 m13; m12 m22 m23; m13 m23 m33] + m11 = a[6] * a[4] - a[5]^2 + m12 = -(a[6] * a[2] - a[5] * a[3]) + m13 = (a[5] * a[2] - a[4] * a[3]) + m22 = (a[1] * a[6] - a[3]^2) + m23 = -(a[5] * a[1] - a[2] * a[3]) + m33 = (a[1] * a[4] - a[2]^2) + + # Solving least square problem + for v in eachvariable(equations) + reconstruction_gradient[1, v, stencil, element] = 1 / AT_A_determinant * (m11 * d[v, 1] + m12 * d[v, 2] + m13 * d[v, 3]) + reconstruction_gradient[2, v, stencil, element] = 1 / AT_A_determinant * (m12 * d[v, 1] + m22 * d[v, 2] + m23 * d[v, 3]) + reconstruction_gradient[3, v, stencil, element] = 1 / AT_A_determinant * (m13 * d[v, 1] + m23 * d[v, 2] + m33 * d[v, 3]) + end + end + + # Get limited reconstruction gradient by weighting the just computed + weight_sum = zero(eltype(reconstruction_gradient)) + for v in eachvariable(equations) + reconstruction_gradient_limited[:, v, element] .= zero(eltype(reconstruction_gradient_limited)) + weight_sum = zero(eltype(reconstruction_gradient)) + for j in 1:(n_stencil_neighbors + 1) + gradient = get_node_coords(reconstruction_gradient, equations, solver, + v, j, element) + indicator = sum(gradient .^ 2) + lambda_j = (j == 1) ? lambda[1] : lambda[2] + weight = (lambda_j / (indicator + epsilon))^r + for x in axes(reconstruction_gradient_limited, 1) + reconstruction_gradient_limited[x, v, element] += weight * + gradient[x] + end + weight_sum += weight + end + for x in axes(reconstruction_gradient_limited, 1) + reconstruction_gradient_limited[x, v, element] /= weight_sum + end + end + end + + exchange_gradient_data!(reconstruction_gradient_limited, mesh, equations, solver, + cache) + + return nothing +end + +@inline function calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d, + reconstruction_stencil, reconstruction_distance, + solution_data, equations::AbstractEquations{3}) + for i in 1:n_stencil_neighbors + # stencil=1 contains information from all neighbors + # stencil=2,...,n_stencil_neighbors+1 is computed without (stencil+1)-th neighbor's information + if i + 1 != stencil + neighbor = reconstruction_stencil[element][i] + distance = reconstruction_distance[element][i] + a[1] += distance[1]^2 # = a11 + a[2] += distance[1] * distance[2] # = a12 + a[3] += distance[1] * distance[3] + a[4] += distance[2]^2 + a[5] += distance[2] * distance[3] + a[6] += distance[3]^2 + + for v in eachvariable(equations) + d[v, 1] += distance[1] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + d[v, 2] += distance[2] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + d[v, 3] += distance[3] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + end + end + end + + return nothing +end + function prolong2interfaces!(cache, mesh::T8codeMesh, equations, solver::FV) (; interfaces, communication_data) = cache (; solution_data, domain_data, gradient_data) = communication_data From df555d212b10d7ae7c741b52efc531b203162c4f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 19 Aug 2024 15:54:38 +0200 Subject: [PATCH 02/12] Complete 3d support --- .../t8code_3d_fv/elixir_advection_basic.jl | 47 +++++++ .../elixir_advection_basic_hybrid.jl | 6 +- .../t8code_3d_fv/elixir_advection_gauss.jl | 44 +++++++ src/equations/linear_scalar_advection_3d.jl | 14 ++- src/meshes/t8code_mesh.jl | 11 +- src/solvers/fv_t8code/containers.jl | 1 - src/solvers/fv_t8code/fv.jl | 117 +++++++++++++++--- 7 files changed, 209 insertions(+), 31 deletions(-) create mode 100644 examples/t8code_3d_fv/elixir_advection_basic.jl create mode 100644 examples/t8code_3d_fv/elixir_advection_gauss.jl diff --git a/examples/t8code_3d_fv/elixir_advection_basic.jl b/examples/t8code_3d_fv/elixir_advection_basic.jl new file mode 100644 index 00000000000..071ef23dc8d --- /dev/null +++ b/examples/t8code_3d_fv/elixir_advection_basic.jl @@ -0,0 +1,47 @@ +using OrdinaryDiffEq +using Trixi + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +initial_condition = initial_condition_convergence_test + +solver = FV(order = 2, extended_reconstruction_stencil = false, + surface_flux = flux_lax_friedrichs) + +coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z)) + +trees_per_dimension = (8, 8, 8) + +mesh = T8codeMesh(trees_per_dimension, + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = 3) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback, save_solution) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() diff --git a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl index 67c664471d5..a9528182320 100644 --- a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl +++ b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl @@ -13,11 +13,11 @@ solver = FV(order = 2, extended_reconstruction_stencil = false, cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true)) # cmesh = Trixi.cmesh_new_periodic_tri() mesh = T8codeMesh(cmesh, solver, - initial_refinement_level = 5) + initial_refinement_level = 3) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) -ode = semidiscretize(semi, (0.0, 10.0)); +ode = semidiscretize(semi, (0.0, 1.0)); summary_callback = SummaryCallback() @@ -29,7 +29,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval) save_solution = SaveSolutionCallback(interval = 1, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.3) +stepsize_callback = StepsizeCallback(cfl = 0.9) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback, save_solution) diff --git a/examples/t8code_3d_fv/elixir_advection_gauss.jl b/examples/t8code_3d_fv/elixir_advection_gauss.jl new file mode 100644 index 00000000000..5c22d4fd6f9 --- /dev/null +++ b/examples/t8code_3d_fv/elixir_advection_gauss.jl @@ -0,0 +1,44 @@ +using OrdinaryDiffEq +using Trixi + +#################################################### + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +initial_condition = initial_condition_gauss + +solver = FV(order = 2, surface_flux = flux_lax_friedrichs) + +initial_refinement_level = 4 +# cmesh = Trixi.cmesh_new_periodic_hybrid() +cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true)) +mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +tspan = (0.0, 5.0) +ode = semidiscretize(semi, tspan); + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.7) + +callbacks = CallbackSet(summary_callback, save_solution, analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),#Euler(), + dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index 088f934cc3e..006179c1fc9 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/src/equations/linear_scalar_advection_3d.jl @@ -64,14 +64,24 @@ function initial_condition_convergence_test(x, t, return SVector(scalar) end +# Calculates translated coordinates `x` for a periodic domain +function x_trans_periodic_3d(x, domain_length = SVector(2, 2, 2), + center = SVector(0, 0, 0)) + x_normalized = x .- center + x_shifted = x_normalized .% domain_length + x_offset = ((x_shifted .< -0.5f0 * domain_length) - + (x_shifted .> 0.5f0 * domain_length)) .* domain_length + return center + x_shifted + x_offset +end + """ - initial_condition_gauss(x, t, equations::LinearScalarAdvectionEquation1D) + initial_condition_gauss(x, t, equations::LinearScalarAdvectionEquation3D) A Gaussian pulse. """ function initial_condition_gauss(x, t, equation::LinearScalarAdvectionEquation3D) # Store translated coordinate for easy use of exact solution - x_trans = x - equation.advection_velocity * t + x_trans = x_trans_periodic_3d(x - equation.advection_velocity * t) scalar = exp(-(x_trans[1]^2 + x_trans[2]^2 + x_trans[3]^2)) return SVector(scalar) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 0dafdf0af9e..65ad8377fe3 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -1597,11 +1597,12 @@ function cmesh_quad_3d(; comm = mpi_comm(), periodicity = (true, true, true))::t linear_geom = t8_geometry_linear_new(n_dims) # This is how the cmesh looks like. The numbers are the tree numbers: - # - # +---+ - # | | - # | 0 | - # | | + # +---+ + # / /| + # +---+ | + # | | | + # | 0 | + + # | |/ # +---+ # diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl index 44dcd19a5e5..635a391b4cb 100644 --- a/src/solvers/fv_t8code/containers.jl +++ b/src/solvers/fv_t8code/containers.jl @@ -371,7 +371,6 @@ function init_reconstruction_stencil!(elements, interfaces, boundaries, return nothing end - # Container data structure (structure-of-arrays style) for FV interfaces mutable struct T8codeFVInterfaceContainer{uEltype <: Real} <: AbstractContainer u::Array{uEltype, 3} # [primary/secondary, variable, interface] diff --git a/src/solvers/fv_t8code/fv.jl b/src/solvers/fv_t8code/fv.jl index 6445a2d7b23..1e4646c89a3 100644 --- a/src/solvers/fv_t8code/fv.jl +++ b/src/solvers/fv_t8code/fv.jl @@ -248,9 +248,12 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver epsilon = 1.0e-13 for element in eachelement(mesh, solver, cache) + # The actual number of used stencils is `n_stencil_neighbors + 1`, since the full stencil is additionally used once. if solver.extended_reconstruction_stencil + # Number of faces = Number of corners n_stencil_neighbors = elements.num_faces[element] else + # Number of direct (edge) neighbors n_stencil_neighbors = length(reconstruction_stencil[element]) end @@ -262,7 +265,7 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver d .= zero(eltype(u)) if solver.extended_reconstruction_stencil - calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d, + calc_gradient_extended!(stencil, element, a, d, reconstruction_stencil, reconstruction_distance, reconstruction_corner_elements, solution_data, equations, solver, cache) @@ -302,14 +305,14 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver indicator = sum(gradient .^ 2) lambda_j = (j == 1) ? lambda[1] : lambda[2] weight = (lambda_j / (indicator + epsilon))^r - for x in axes(reconstruction_gradient_limited, 1) - reconstruction_gradient_limited[x, v, element] += weight * - gradient[x] + for dim in axes(reconstruction_gradient_limited, 1) + reconstruction_gradient_limited[dim, v, element] += weight * + gradient[dim] end weight_sum += weight end - for x in axes(reconstruction_gradient_limited, 1) - reconstruction_gradient_limited[x, v, element] /= weight_sum + for dim in axes(reconstruction_gradient_limited, 1) + reconstruction_gradient_limited[dim, v, element] /= weight_sum end end end @@ -347,11 +350,11 @@ end return nothing end -@inline function calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d, +@inline function calc_gradient_extended!(stencil, element, a, d, reconstruction_stencil, reconstruction_distance, reconstruction_corner_elements, solution_data, - equations, solver, cache) + equations::AbstractEquations{2}, solver, cache) if stencil == 1 # Full stencil for i in eachindex(reconstruction_stencil[element]) @@ -432,10 +435,12 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver epsilon = 1.0e-13 for element in eachelement(mesh, solver, cache) + # The actual number of used stencils is `n_stencil_neighbors + 1`, since the full stencil is additionally used once. if solver.extended_reconstruction_stencil + # Number of faces = Number of corners n_stencil_neighbors = elements.num_faces[element] - error("not yet supported!") else + # Number of direct (surface) neighbors n_stencil_neighbors = length(reconstruction_stencil[element]) end @@ -447,11 +452,10 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver d .= zero(eltype(u)) if solver.extended_reconstruction_stencil - calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d, + calc_gradient_extended!(stencil, element, a, d, reconstruction_stencil, reconstruction_distance, reconstruction_corner_elements, solution_data, equations, solver, cache) - error("Not supported yet") else calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d, reconstruction_stencil, reconstruction_distance, @@ -460,7 +464,9 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver # Divide by determinant # Determinant = a11 * (a22 * a33 - a23^2 ) - a12 * (a12 * a33 - a13 * a23 ) + a13 * (a12 * a23 - a22 * a13 ) - AT_A_determinant = a[1] * (a[4] * a[6] - a[5]^2) - a[2] * (a[2] * a[6] - a[3] * a[5]) + a[3] * (a[2] * a[5] - a[4] * a[3]) + AT_A_determinant = a[1] * (a[4] * a[6] - a[5]^2) - + a[2] * (a[2] * a[6] - a[3] * a[5]) + + a[3] * (a[2] * a[5] - a[4] * a[3]) if isapprox(AT_A_determinant, 0.0) for v in eachvariable(equations) reconstruction_gradient[:, v, stencil, element] .= 0.0 @@ -482,9 +488,18 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver # Solving least square problem for v in eachvariable(equations) - reconstruction_gradient[1, v, stencil, element] = 1 / AT_A_determinant * (m11 * d[v, 1] + m12 * d[v, 2] + m13 * d[v, 3]) - reconstruction_gradient[2, v, stencil, element] = 1 / AT_A_determinant * (m12 * d[v, 1] + m22 * d[v, 2] + m23 * d[v, 3]) - reconstruction_gradient[3, v, stencil, element] = 1 / AT_A_determinant * (m13 * d[v, 1] + m23 * d[v, 2] + m33 * d[v, 3]) + reconstruction_gradient[1, v, stencil, element] = 1 / AT_A_determinant * + (m11 * d[v, 1] + + m12 * d[v, 2] + + m13 * d[v, 3]) + reconstruction_gradient[2, v, stencil, element] = 1 / AT_A_determinant * + (m12 * d[v, 1] + + m22 * d[v, 2] + + m23 * d[v, 3]) + reconstruction_gradient[3, v, stencil, element] = 1 / AT_A_determinant * + (m13 * d[v, 1] + + m23 * d[v, 2] + + m33 * d[v, 3]) end end @@ -499,14 +514,14 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver indicator = sum(gradient .^ 2) lambda_j = (j == 1) ? lambda[1] : lambda[2] weight = (lambda_j / (indicator + epsilon))^r - for x in axes(reconstruction_gradient_limited, 1) - reconstruction_gradient_limited[x, v, element] += weight * - gradient[x] + for dim in axes(reconstruction_gradient_limited, 1) + reconstruction_gradient_limited[dim, v, element] += weight * + gradient[dim] end weight_sum += weight end - for x in axes(reconstruction_gradient_limited, 1) - reconstruction_gradient_limited[x, v, element] /= weight_sum + for dim in axes(reconstruction_gradient_limited, 1) + reconstruction_gradient_limited[dim, v, element] /= weight_sum end end end @@ -550,6 +565,68 @@ end return nothing end +@inline function calc_gradient_extended!(stencil, element, a, d, + reconstruction_stencil, + reconstruction_distance, + reconstruction_corner_elements, solution_data, + equations::AbstractEquations{3}, solver, cache) + if stencil == 1 + # Full stencil + for i in eachindex(reconstruction_stencil[element]) + neighbor = reconstruction_stencil[element][i] + distance = reconstruction_distance[element][i] + a[1] += distance[1]^2 # = a11 + a[2] += distance[1] * distance[2] # = a12 + a[3] += distance[1] * distance[3] + a[4] += distance[2]^2 + a[5] += distance[2] * distance[3] + a[6] += distance[3]^2 + + for v in eachvariable(equations) + d[v, 1] += distance[1] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + d[v, 2] += distance[2] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + d[v, 3] += distance[3] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + end + end + else + # Partial stencils + midpoint_element = get_node_coords(cache.elements.midpoint, equations, solver, + element) + for neighbor in reconstruction_corner_elements[stencil - 1, element] + midpoint_neighbor = get_node_coords(cache.elements.midpoint, equations, + solver, neighbor) + distance = midpoint_neighbor .- midpoint_element + + a[1] += distance[1]^2 # = a11 + a[2] += distance[1] * distance[2] # = a12 + a[3] += distance[1] * distance[3] + a[4] += distance[2]^2 + a[5] += distance[2] * distance[3] + a[6] += distance[3]^2 + + for v in eachvariable(equations) + d[v, 1] += distance[1] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + d[v, 2] += distance[2] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + d[v, 3] += distance[3] * + (solution_data[neighbor].u[v] - + solution_data[element].u[v]) + end + end + end + + return nothing +end + function prolong2interfaces!(cache, mesh::T8codeMesh, equations, solver::FV) (; interfaces, communication_data) = cache (; solution_data, domain_data, gradient_data) = communication_data From 55a30ca92d8fd289732054925357fe4eaca43df4 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 20 Aug 2024 11:31:17 +0200 Subject: [PATCH 03/12] Add 3d elixirs and tests --- .../elixir_advection_basic_hybrid.jl | 5 +- .../t8code_3d_fv/elixir_advection_gauss.jl | 3 +- .../elixir_advection_nonperiodic.jl | 53 ++++ .../t8code_3d_fv/elixir_euler_source_terms.jl | 51 ++++ src/solvers/fv_t8code/containers.jl | 2 +- test/runtests.jl | 3 +- ...test_t8code_fv.jl => test_t8code_fv_2d.jl} | 0 test/test_t8code_fv_3d.jl | 256 ++++++++++++++++++ 8 files changed, 367 insertions(+), 6 deletions(-) create mode 100644 examples/t8code_3d_fv/elixir_advection_nonperiodic.jl create mode 100644 examples/t8code_3d_fv/elixir_euler_source_terms.jl rename test/{test_t8code_fv.jl => test_t8code_fv_2d.jl} (100%) create mode 100644 test/test_t8code_fv_3d.jl diff --git a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl index a9528182320..f07ae5bed3d 100644 --- a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl +++ b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl @@ -9,11 +9,10 @@ initial_condition = initial_condition_convergence_test solver = FV(order = 2, extended_reconstruction_stencil = false, surface_flux = flux_lax_friedrichs) -# cmesh = Trixi.cmesh_new_periodic_hybrid() +# TODO: There are no other cmesh functions implemented yet in 3d. cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true)) -# cmesh = Trixi.cmesh_new_periodic_tri() mesh = T8codeMesh(cmesh, solver, - initial_refinement_level = 3) + initial_refinement_level = 4) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_3d_fv/elixir_advection_gauss.jl b/examples/t8code_3d_fv/elixir_advection_gauss.jl index 5c22d4fd6f9..ef223dd6aeb 100644 --- a/examples/t8code_3d_fv/elixir_advection_gauss.jl +++ b/examples/t8code_3d_fv/elixir_advection_gauss.jl @@ -11,7 +11,8 @@ initial_condition = initial_condition_gauss solver = FV(order = 2, surface_flux = flux_lax_friedrichs) initial_refinement_level = 4 -# cmesh = Trixi.cmesh_new_periodic_hybrid() + +# TODO: There are no other cmesh functions implemented yet in 3d. cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) diff --git a/examples/t8code_3d_fv/elixir_advection_nonperiodic.jl b/examples/t8code_3d_fv/elixir_advection_nonperiodic.jl new file mode 100644 index 00000000000..22ba4037898 --- /dev/null +++ b/examples/t8code_3d_fv/elixir_advection_nonperiodic.jl @@ -0,0 +1,53 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the linear advection equation. + +advection_velocity = (0.2, -0.7, 0.5) +equations = LinearScalarAdvectionEquation3D(advection_velocity) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = Dict(:all => boundary_condition) +# Problem: T8codeMesh interface with parameter cmesh cannot distinguish between boundaries +# boundary_conditions = Dict(:x_neg => boundary_condition, +# :x_pos => boundary_condition, +# :y_neg => boundary_condition_periodic, +# :y_pos => boundary_condition_periodic) + +solver = FV(order = 2, surface_flux = flux_lax_friedrichs) + +# TODO: When using mesh construction as in elixir_advection_basic.jl boundary Symbol :all is not defined +initial_refinement_level = 5 +cmesh = Trixi.cmesh_quad_3d(periodicity = (false, false, false)) +mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# Run the simulation. + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # Solve needs some value here but it will be overwritten by the stepsize_callback. + save_everystep = false, callback = callbacks); +summary_callback() diff --git a/examples/t8code_3d_fv/elixir_euler_source_terms.jl b/examples/t8code_3d_fv/elixir_euler_source_terms.jl new file mode 100644 index 00000000000..8568a1dac96 --- /dev/null +++ b/examples/t8code_3d_fv/elixir_euler_source_terms.jl @@ -0,0 +1,51 @@ + +using OrdinaryDiffEq +using Trixi + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +# boundary_condition = BoundaryConditionDirichlet(initial_condition) +# boundary_conditions = Dict(:all => boundary_condition) + +source_terms = source_terms_convergence_test + +solver = FV(order = 2, extended_reconstruction_stencil = false, + surface_flux = flux_lax_friedrichs) + +# TODO: There are no other cmesh functions implemented yet in 3d. +cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true)) +mesh = T8codeMesh(cmesh, solver, + initial_refinement_level = 3) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms) +# boundary_conditions = boundary_conditions) + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl index 635a391b4cb..52278e0b69a 100644 --- a/src/solvers/fv_t8code/containers.jl +++ b/src/solvers/fv_t8code/containers.jl @@ -352,7 +352,7 @@ function init_reconstruction_stencil!(elements, interfaces, boundaries, face_midpoint_element1 = domain_data[element1].face_midpoints[face_element1] face_midpoint_element2 = domain_data[element2].face_midpoints[face_element2] - # How to handle periodic boundaries? + # TODO: How to handle periodic boundaries? if isapprox(face_midpoint_element1, face_midpoint_element2) distance = midpoint_element2 .- midpoint_element1 else diff --git a/test/runtests.jl b/test/runtests.jl index 1ecb13c7fab..75bd210d0eb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -94,7 +94,8 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) end @time if TRIXI_TEST == "all" || TRIXI_TEST == "t8code_fv" - include("test_t8code_fv.jl") + include("test_t8code_fv_2d.jl") + include("test_t8code_fv_3d.jl") end @time if TRIXI_TEST == "all" || TRIXI_TEST == "unstructured_dgmulti" diff --git a/test/test_t8code_fv.jl b/test/test_t8code_fv_2d.jl similarity index 100% rename from test/test_t8code_fv.jl rename to test/test_t8code_fv_2d.jl diff --git a/test/test_t8code_fv_3d.jl b/test/test_t8code_fv_3d.jl new file mode 100644 index 00000000000..ec980288af1 --- /dev/null +++ b/test/test_t8code_fv_3d.jl @@ -0,0 +1,256 @@ +module TestExamplesT8codeMesh2D + +using Test +using Trixi + +include("test_trixi.jl") + +# I added this temporary test file for constantly testing while developing. +# The tests have to be adapted at the end. +EXAMPLES_DIR = joinpath(examples_dir(), "t8code_3d_fv") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) +mkdir(outdir) + +@testset "T8codeMesh3D" begin +#! format: noindent + +# @trixi_testset "test save_mesh_file" begin +# @test_throws Exception begin +# # Save mesh file support will be added in the future. The following +# # lines of code are here for satisfying code coverage. + +# # Create dummy mesh. +# mesh = T8codeMesh((1, 1), polydeg = 1, +# mapping = Trixi.coordinates2mapping((-1.0, -1.0), (1.0, 1.0)), +# initial_refinement_level = 1) + +# # This call throws an error. +# Trixi.save_mesh_file(mesh, "dummy") +# end +# end + +# @trixi_testset "test check_for_negative_volumes" begin +# @test_warn "Discovered negative volumes" begin +# # Unstructured mesh with six cells which have left-handed node ordering. +# mesh_file = Trixi.download("https://gist.githubusercontent.com/jmark/bfe0d45f8e369298d6cc637733819013/raw/cecf86edecc736e8b3e06e354c494b2052d41f7a/rectangle_with_negative_volumes.msh", +# joinpath(EXAMPLES_DIR, +# "rectangle_with_negative_volumes.msh")) + +# # This call should throw a warning about negative volumes detected. +# mesh = T8codeMesh(mesh_file, 2) +# end +# end + +@trixi_testset "elixir_advection_basic.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + order=1, + initial_refinement_level=2, + l2=[0.2848617953369851], + linf=[0.3721898718954475]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + initial_refinement_level=2, + l2=[0.10381089565603231], + linf=[0.13787405651527007]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV, extended reconstruction stencil" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + initial_refinement_level=2, + extended_reconstruction_stencil=true, + l2=[0.24826100771542928], + linf=[0.3799973662927152]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + +@trixi_testset "elixir_advection_gauss.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_gauss.jl"), + order=1, + l2=[0.1515258539168874], + linf=[0.43164936150417055]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_gauss.jl"), + l2=[0.04076672839289378], + linf=[0.122537463101035582]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + +@trixi_testset "elixir_advection_basic_hybrid.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + order=1, + l2=[0.20282363730327146], + linf=[0.28132446651281295]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + l2=[0.02153993127089835], + linf=[0.039109618097251886]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end + +@trixi_testset "elixir_advection_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"), + l2=[0.022202106950138526], + linf=[0.0796166790338586]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_euler_source_terms.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + order=1, + l2=[ + 0.050763790354290725, + 0.0351299673616484, + 0.0351299673616484, + 0.03512996736164839, + 0.1601847269543808, + ], + linf=[ + 0.07175521415072939, + 0.04648499338897771, + 0.04648499338897816, + 0.04648499338897816, + 0.2235470564880404, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + order=2, + l2=[ + 0.012308219704695382, + 0.010791416898840429, + 0.010791416898840464, + 0.010791416898840377, + 0.036995680347196136, + ], + linf=[ + 0.01982294164697862, + 0.01840725612418126, + 0.01840725612418148, + 0.01840725612418148, + 0.05736595182767079, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV, extended reconstruction stencil" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + order=2, + extended_reconstruction_stencil=true, + l2=[ + 0.05057867333486591, + 0.03596196296013507, + 0.03616867188152877, + 0.03616867188152873, + 0.14939041550302212, + ], + linf=[ + 0.07943789383956079, + 0.06389365911606859, + 0.06469291944863809, + 0.0646929194486372, + 0.23507781748792533, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end +end + +# Clean up afterwards: delete Trixi.jl output directory +@test_nowarn rm(outdir, recursive = true) + +end # module From 5d7ea73e980f744b7e42b9d5a96471ba53d8f5b0 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 14 Oct 2024 17:41:35 +0200 Subject: [PATCH 04/12] Adapt name in elixir --- examples/t8code_3d_fv/elixir_advection_basic.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/t8code_3d_fv/elixir_advection_basic.jl b/examples/t8code_3d_fv/elixir_advection_basic.jl index 1b98ad11524..a9e08bf897e 100644 --- a/examples/t8code_3d_fv/elixir_advection_basic.jl +++ b/examples/t8code_3d_fv/elixir_advection_basic.jl @@ -26,7 +26,7 @@ function mapping(xi, eta, zeta) return SVector(x, y, z) end -function f(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data) +function trixi_t8_mapping(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data) ltreeid = t8_cmesh_get_local_id(cmesh, gtreeid) eclass = t8_cmesh_get_tree_class(cmesh, ltreeid) T8code.t8_geom_compute_linear_geometry(eclass, tree_data, @@ -49,8 +49,8 @@ function f(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_d return nothing end -function f_c() - @cfunction($f, Cvoid, +function trixi_t8_mapping_c() + @cfunction($trixi_t8_mapping, Cvoid, (t8_cmesh_t, t8_gloidx_t, Ptr{Cdouble}, Csize_t, Ptr{Cdouble}, Ptr{Cvoid}, Ptr{Cvoid})) end @@ -59,7 +59,7 @@ trees_per_dimension = (2, 2, 2) eclass = T8_ECLASS_HEX mesh = T8codeMesh(trees_per_dimension, eclass; - mapping = f_c(), + mapping = trixi_t8_mapping_c(), # Plan is to use either # coordinates_max = coordinates_max, coordinates_min = coordinates_min, # or at least From d6a5a6d3d4e5dd602425de8b11090ec0a8a01a32 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 14 Oct 2024 17:48:04 +0200 Subject: [PATCH 05/12] Reduce initial refinement level of one test --- test/test_t8code_fv_3d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_t8code_fv_3d.jl b/test/test_t8code_fv_3d.jl index ef14531d91f..5f09902198d 100644 --- a/test/test_t8code_fv_3d.jl +++ b/test/test_t8code_fv_3d.jl @@ -77,10 +77,10 @@ mkdir(outdir) end @trixi_testset "second-order FV, extended reconstruction stencil" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - initial_refinement_level=2+2, + initial_refinement_level=1+2, extended_reconstruction_stencil=true, - l2=[0.24826100771542928], - linf=[0.3799973662927152]) + l2=[0.3282177575292713], + linf=[0.39002345444858333]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 9b09a4bdaa9fc17706b872f1849c95f0618ffe6c Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 14 Oct 2024 18:18:16 +0200 Subject: [PATCH 06/12] Simplify NDIMS=3 --- src/meshes/t8code_mesh.jl | 160 +++++++++++++++----------------------- 1 file changed, 63 insertions(+), 97 deletions(-) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 469dc62841b..ca1bf70396a 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -384,22 +384,23 @@ function T8codeMesh(trees_per_dimension, eclass; end do_partition = 0 - if NDIMS == 2 - cmesh_ref = Ref(t8_cmesh_t()) - t8_cmesh_init(cmesh_ref) - cmesh = cmesh_ref[] + cmesh_ref = Ref(t8_cmesh_t()) + t8_cmesh_init(cmesh_ref) + cmesh = cmesh_ref[] - @assert(allequal(trees_per_dimension), + @assert(allequal(trees_per_dimension), "Different trees per dimensions are not supported for quad mesh. `trees_per_dimension`: $trees_per_dimension") - num_trees = prod(trees_per_dimension) * 1 # 1 = number of trees for single hypercube with quads + num_trees = prod(trees_per_dimension) * 1 # 1 = number of trees for single hypercube with quads + vertices_per_tree = 2^NDIMS # number of vertices (=corners) in one tree + vertices = Vector{Cdouble}(undef, 3 * vertices_per_tree * num_trees) # 3 (dimensions) * vertices_per_tree) * num_trees + if NDIMS == 2 coordinates_tree_x = range(-1.0, 1.0, length = trees_per_dimension[1] + 1) coordinates_tree_y = range(-1.0, 1.0, length = trees_per_dimension[2] + 1) - vertices = Vector{Cdouble}(undef, 3 * 4 * num_trees) # 3 (dimensions) * 4 (vertices per tree) * num_trees for itree_y in 1:trees_per_dimension[2], itree_x in 1:trees_per_dimension[1] - index = trees_per_dimension[1] * 3 * 4 * (itree_y - 1) + - 3 * 4 * (itree_x - 1) + 1 + index = trees_per_dimension[1] * 3 * vertices_per_tree * (itree_y - 1) + + 3 * vertices_per_tree * (itree_x - 1) + 1 vertices[index] = coordinates_tree_x[itree_x] vertices[index + 1] = coordinates_tree_y[itree_y] vertices[index + 2] = 0.0 @@ -419,74 +420,15 @@ function T8codeMesh(trees_per_dimension, eclass; vertices[index + 1] = coordinates_tree_y[itree_y + 1] vertices[index + 2] = 0.0 end - - # analytical = trixi_t8_mapping_c(mapping) - analytical = mapping - name = "" - user_data = vertices - # user_data = C_NULL - - jacobian = C_NULL # type: t8_geom_analytic_jacobian_fn - # TODO: Since @t8_load_tree_data is not yet included into T8code.jl, precompiling Trixi fails because of this line. - # load_tree_data = @t8_load_tree_data(t8_geom_load_tree_data_vertices) # type: t8_geom_load_tree_data_fn - # For now, I remove it and pass C_NULL. Then, `elixir_advection_basic.jl` will fail with a SegFault. - load_tree_data = C_NULL - tree_negative_volume = C_NULL # type: t8_geom_tree_negative_volume_fn - - geometry = t8_geometry_analytic_new(NDIMS, name, analytical, jacobian, - load_tree_data, tree_negative_volume, user_data) - - t8_cmesh_register_geometry(cmesh, geometry) - - for itree in 1:num_trees - offset_vertices = 3 * 4 * (itree - 1) - t8_cmesh_set_tree_class(cmesh, itree - 1, eclass) - t8_cmesh_set_tree_vertices(cmesh, itree - 1, - @views(vertices[(1 + offset_vertices):end]), 4) - end - - # Note and TODO: - # Only hardcoded to `trees_per_dimension = (1, 1) or (2, 2)` - if num_trees == 1 - if periodicity[1] - t8_cmesh_set_join(cmesh, 0, 0, 0, 1, 0) - end - if periodicity[2] - t8_cmesh_set_join(cmesh, 0, 0, 2, 3, 0) - end - elseif num_trees == 4 - if periodicity[1] - t8_cmesh_set_join(cmesh, 0, 1, 0, 1, 0) - t8_cmesh_set_join(cmesh, 2, 3, 0, 1, 0) - end - if periodicity[2] - t8_cmesh_set_join(cmesh, 0, 2, 2, 3, 0) - t8_cmesh_set_join(cmesh, 1, 3, 2, 3, 0) - end - t8_cmesh_set_join_by_stash(cmesh, C_NULL, 0) - else - error("Not supported trees_per_dimension") - end - - t8_cmesh_commit(cmesh, mpi_comm()) elseif NDIMS == 3 - cmesh_ref = Ref(t8_cmesh_t()) - t8_cmesh_init(cmesh_ref) - cmesh = cmesh_ref[] - - @assert(allequal(trees_per_dimension), - "Different trees per dimensions are not supported for quad mesh. `trees_per_dimension`: $trees_per_dimension") - num_trees = prod(trees_per_dimension) * 1 # 1 = number of trees for single hypercube with quads - coordinates_tree_x = range(-1.0, 1.0, length = trees_per_dimension[1] + 1) coordinates_tree_y = range(-1.0, 1.0, length = trees_per_dimension[2] + 1) coordinates_tree_z = range(-1.0, 1.0, length = trees_per_dimension[3] + 1) - vertices = Vector{Cdouble}(undef, 3 * 8 * num_trees) # 3 (dimensions) * 8 (vertices per tree) * num_trees for itree_z in 1:trees_per_dimension[3], itree_y in 1:trees_per_dimension[2], itree_x in 1:trees_per_dimension[1] - index = trees_per_dimension[1] * trees_per_dimension[2] * 3 * 8 * (itree_z - 1) + - trees_per_dimension[1] * 3 * 8 * (itree_y - 1) + - 3 * 8 * (itree_x - 1) + 1 + index = trees_per_dimension[1] * trees_per_dimension[2] * 3 * vertices_per_tree * (itree_z - 1) + + trees_per_dimension[1] * 3 * vertices_per_tree * (itree_y - 1) + + 3 * vertices_per_tree * (itree_x - 1) + 1 vertices[index] = coordinates_tree_x[itree_x] vertices[index + 1] = coordinates_tree_y[itree_y] @@ -528,34 +470,59 @@ function T8codeMesh(trees_per_dimension, eclass; vertices[index + 1] = coordinates_tree_y[itree_y + 1] vertices[index + 2] = coordinates_tree_z[itree_z + 1] end + end - # analytical = trixi_t8_mapping_c(mapping) - analytical = mapping - name = "" - user_data = vertices - # user_data = C_NULL - - jacobian = C_NULL # type: t8_geom_analytic_jacobian_fn - # TODO: Since @t8_load_tree_data is not yet included into T8code.jl, precompiling Trixi fails because of this line. - load_tree_data = @t8_load_tree_data(t8_geom_load_tree_data_vertices) # type: t8_geom_load_tree_data_fn - # For now, I remove it and pass C_NULL. Then, `elixir_advection_basic.jl` will fail with a SegFault. - # load_tree_data = C_NULL - tree_negative_volume = C_NULL # type: t8_geom_tree_negative_volume_fn - - geometry = t8_geometry_analytic_new(NDIMS, name, analytical, jacobian, - load_tree_data, tree_negative_volume, user_data) - - t8_cmesh_register_geometry(cmesh, geometry) - - for itree in 1:num_trees - offset_vertices = 3 * 8 * (itree - 1) - t8_cmesh_set_tree_class(cmesh, itree - 1, eclass) - t8_cmesh_set_tree_vertices(cmesh, itree - 1, - @views(vertices[(1 + offset_vertices):end]), 8) - end + # analytical = trixi_t8_mapping_c(mapping) + analytical = mapping + name = "" + user_data = vertices + # user_data = C_NULL + + jacobian = C_NULL # type: t8_geom_analytic_jacobian_fn + # TODO: Since @t8_load_tree_data is not yet included into T8code.jl, precompiling Trixi fails because of this line. + # load_tree_data = @t8_load_tree_data(t8_geom_load_tree_data_vertices) # type: t8_geom_load_tree_data_fn + # For now, I remove it and pass C_NULL. Then, `elixir_advection_basic.jl` will fail with a SegFault. + load_tree_data = C_NULL + tree_negative_volume = C_NULL # type: t8_geom_tree_negative_volume_fn + + geometry = t8_geometry_analytic_new(NDIMS, name, analytical, jacobian, + load_tree_data, tree_negative_volume, user_data) + + t8_cmesh_register_geometry(cmesh, geometry) + + for itree in 1:num_trees + offset_vertices = 3 * vertices_per_tree * (itree - 1) + t8_cmesh_set_tree_class(cmesh, itree - 1, eclass) + t8_cmesh_set_tree_vertices(cmesh, itree - 1, + @views(vertices[(1 + offset_vertices):end]), vertices_per_tree) + end + if NDIMS == 2 # Note and TODO: # Only hardcoded to `trees_per_dimension = (1, 1) or (2, 2)` + if num_trees == 1 + if periodicity[1] + t8_cmesh_set_join(cmesh, 0, 0, 0, 1, 0) + end + if periodicity[2] + t8_cmesh_set_join(cmesh, 0, 0, 2, 3, 0) + end + elseif num_trees == 4 + if periodicity[1] + t8_cmesh_set_join(cmesh, 0, 1, 0, 1, 0) + t8_cmesh_set_join(cmesh, 2, 3, 0, 1, 0) + end + if periodicity[2] + t8_cmesh_set_join(cmesh, 0, 2, 2, 3, 0) + t8_cmesh_set_join(cmesh, 1, 3, 2, 3, 0) + end + t8_cmesh_set_join_by_stash(cmesh, C_NULL, 0) + else + error("Not supported trees_per_dimension") + end + elseif NDIMS == 3 + # Note and TODO: + # Only hardcoded to `trees_per_dimension = (1, 1, 1) or (2, 2, 2)` if num_trees == 1 if periodicity[1] t8_cmesh_set_join(cmesh, 0, 0, 0, 1, 0) @@ -589,9 +556,8 @@ function T8codeMesh(trees_per_dimension, eclass; else error("Not supported trees_per_dimension") end - - t8_cmesh_commit(cmesh, mpi_comm()) end + t8_cmesh_commit(cmesh, mpi_comm()) do_face_ghost = mpi_isparallel() scheme = t8_scheme_new_default_cxx() From 32dab35d07939db1242e4e397c2ed27a36d8faaf Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 15 Oct 2024 10:36:34 +0200 Subject: [PATCH 07/12] Add 3d mpi tests --- test/test_mpi.jl | 3 +- ..._t8code_fv.jl => test_mpi_t8code_fv_2d.jl} | 15 --- test/test_mpi_t8code_fv_3d.jl | 95 +++++++++++++++++++ 3 files changed, 97 insertions(+), 16 deletions(-) rename test/{test_mpi_t8code_fv.jl => test_mpi_t8code_fv_2d.jl} (95%) create mode 100644 test/test_mpi_t8code_fv_3d.jl diff --git a/test/test_mpi.jl b/test/test_mpi.jl index f36b669d3c3..4f9eb65ecd8 100644 --- a/test/test_mpi.jl +++ b/test/test_mpi.jl @@ -23,10 +23,11 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows() # P4estMesh and T8codeMesh tests include("test_mpi_p4est_2d.jl") include("test_mpi_t8code_2d.jl") - include("test_mpi_t8code_fv.jl") + include("test_mpi_t8code_fv_2d.jl") if !CI_ON_WINDOWS # see comment on `CI_ON_WINDOWS` above include("test_mpi_p4est_3d.jl") include("test_mpi_t8code_3d.jl") + include("test_mpi_t8code_fv_3d.jl") end end # MPI diff --git a/test/test_mpi_t8code_fv.jl b/test/test_mpi_t8code_fv_2d.jl similarity index 95% rename from test/test_mpi_t8code_fv.jl rename to test/test_mpi_t8code_fv_2d.jl index 5a2d9525558..54b4cd83333 100644 --- a/test/test_mpi_t8code_fv.jl +++ b/test/test_mpi_t8code_fv_2d.jl @@ -43,21 +43,6 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") end # The extended reconstruction stencil is currently not mpi parallel. # The current version runs through but an error occurs on some rank. - # @trixi_testset "second-order FV, extended reconstruction stencil" begin - # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - # order=2, - # extended_reconstruction_stencil=true, - # l2=[0.020331012873518642], - # linf=[0.05571209803860677]) - # # Ensure that we do not have excessive memory allocations - # # (e.g., from type instabilities) - # let - # t = sol.t[end] - # u_ode = sol.u[end] - # du_ode = similar(u_ode) - # @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - # end - # end end @trixi_testset "elixir_advection_gauss.jl" begin diff --git a/test/test_mpi_t8code_fv_3d.jl b/test/test_mpi_t8code_fv_3d.jl new file mode 100644 index 00000000000..5818b777b5a --- /dev/null +++ b/test/test_mpi_t8code_fv_3d.jl @@ -0,0 +1,95 @@ +module TestExamplesMPIT8codeMesh2D + +using Test +using Trixi + +include("test_trixi.jl") + +const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") + +@testset "T8codeMesh MPI FV" begin +#! format: noindent + +# Run basic tests +@testset "Examples 3D" begin + @trixi_testset "elixir_advection_basic.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + order=1, + initial_refinement_level=2+2, + l2=[0.2848617953369851], + linf=[0.3721898718954475]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + initial_refinement_level=2+2, + l2=[0.10381089565603231], + linf=[0.13787405651527007]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + # The extended reconstruction stencil is currently not mpi parallel. + # The current version runs through but an error occurs on some rank. + end + + @trixi_testset "elixir_advection_basic_hybrid.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + order=1, + l2=[0.20282363730327146], + linf=[0.28132446651281295]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + l2=[0.02153993127089835], + linf=[0.039109618097251886]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + end + + @trixi_testset "elixir_advection_nonperiodic.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"), + l2=[0.022202106950138526], + linf=[0.0796166790338586]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end +end # T8codeMesh MPI + +end # module From 8d02b1ec83e39303f727365f1354ebf579ac2b1c Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 15 Oct 2024 11:19:51 +0200 Subject: [PATCH 08/12] Remove initial_condition from rhs; fmt --- .../t8code_2d_fv/elixir_advection_basic.jl | 11 ++++++----- .../t8code_3d_fv/elixir_advection_basic.jl | 3 ++- src/meshes/t8code_mesh.jl | 15 ++++++++++----- src/solvers/fv_t8code/fv.jl | 2 +- test/test_mpi_t8code_fv_3d.jl | 11 +++++++---- test/test_t8code_fv_3d.jl | 18 +++++++++--------- 6 files changed, 35 insertions(+), 25 deletions(-) diff --git a/examples/t8code_2d_fv/elixir_advection_basic.jl b/examples/t8code_2d_fv/elixir_advection_basic.jl index 4e6c13f1d08..e16fc5a7ff1 100644 --- a/examples/t8code_2d_fv/elixir_advection_basic.jl +++ b/examples/t8code_2d_fv/elixir_advection_basic.jl @@ -34,10 +34,10 @@ mapping_coordinates = Trixi.coordinates2mapping(coordinates_min, coordinates_max # f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) # [0,8]^2 -f1(s) = SVector(0.0, 4*(s+1)) -f2(s) = SVector(8.0, 4*(s+1)) -f3(s) = SVector(4*(s+1), 0.0) -f4(s) = SVector(4*(s+1), 8.0) +f1(s) = SVector(0.0, 4 * (s + 1)) +f2(s) = SVector(8.0, 4 * (s + 1)) +f3(s) = SVector(4 * (s + 1), 0.0) +f4(s) = SVector(4 * (s + 1), 8.0) faces = (f1, f2, f3, f4) Trixi.validate_faces(faces) mapping_faces = Trixi.transfinite_mapping(faces) @@ -57,7 +57,8 @@ end # directly within this elixir (e.g. mapping = trixi_t8_mapping_c(mapping)), we get the SegFault error. # - Is the function called with the correct parameters? Memory location correct? It seems so, yes. # - Life time issue for the GC tracked Julia object used in C? **Yes, see gc deactivation in elixir.** -function trixi_t8_mapping(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data) +function trixi_t8_mapping(cmesh, gtreeid, ref_coords, num_coords, out_coords, + tree_data, user_data) ltreeid = t8_cmesh_get_local_id(cmesh, gtreeid) eclass = t8_cmesh_get_tree_class(cmesh, ltreeid) T8code.t8_geom_compute_linear_geometry(eclass, tree_data, diff --git a/examples/t8code_3d_fv/elixir_advection_basic.jl b/examples/t8code_3d_fv/elixir_advection_basic.jl index a9e08bf897e..4fe2f147fd4 100644 --- a/examples/t8code_3d_fv/elixir_advection_basic.jl +++ b/examples/t8code_3d_fv/elixir_advection_basic.jl @@ -26,7 +26,8 @@ function mapping(xi, eta, zeta) return SVector(x, y, z) end -function trixi_t8_mapping(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data) +function trixi_t8_mapping(cmesh, gtreeid, ref_coords, num_coords, out_coords, + tree_data, user_data) ltreeid = t8_cmesh_get_local_id(cmesh, gtreeid) eclass = t8_cmesh_get_tree_class(cmesh, ltreeid) T8code.t8_geom_compute_linear_geometry(eclass, tree_data, diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index ca1bf70396a..164a816aa15 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -389,7 +389,7 @@ function T8codeMesh(trees_per_dimension, eclass; cmesh = cmesh_ref[] @assert(allequal(trees_per_dimension), - "Different trees per dimensions are not supported for quad mesh. `trees_per_dimension`: $trees_per_dimension") + "Different trees per dimensions are not supported for quad mesh. `trees_per_dimension`: $trees_per_dimension") num_trees = prod(trees_per_dimension) * 1 # 1 = number of trees for single hypercube with quads vertices_per_tree = 2^NDIMS # number of vertices (=corners) in one tree @@ -425,8 +425,11 @@ function T8codeMesh(trees_per_dimension, eclass; coordinates_tree_y = range(-1.0, 1.0, length = trees_per_dimension[2] + 1) coordinates_tree_z = range(-1.0, 1.0, length = trees_per_dimension[3] + 1) - for itree_z in 1:trees_per_dimension[3], itree_y in 1:trees_per_dimension[2], itree_x in 1:trees_per_dimension[1] - index = trees_per_dimension[1] * trees_per_dimension[2] * 3 * vertices_per_tree * (itree_z - 1) + + for itree_z in 1:trees_per_dimension[3], itree_y in 1:trees_per_dimension[2], + itree_x in 1:trees_per_dimension[1] + + index = trees_per_dimension[1] * trees_per_dimension[2] * 3 * + vertices_per_tree * (itree_z - 1) + trees_per_dimension[1] * 3 * vertices_per_tree * (itree_y - 1) + 3 * vertices_per_tree * (itree_x - 1) + 1 @@ -494,7 +497,8 @@ function T8codeMesh(trees_per_dimension, eclass; offset_vertices = 3 * vertices_per_tree * (itree - 1) t8_cmesh_set_tree_class(cmesh, itree - 1, eclass) t8_cmesh_set_tree_vertices(cmesh, itree - 1, - @views(vertices[(1 + offset_vertices):end]), vertices_per_tree) + @views(vertices[(1 + offset_vertices):end]), + vertices_per_tree) end if NDIMS == 2 @@ -1808,7 +1812,8 @@ function cmesh_new_quad(; trees_per_dimension = (1, 1), end function cmesh_new_quad_3d(; trees_per_dimension = (1, 1, 1), - coordinates_min = (-1.0, -1.0, -1.0), coordinates_max = (1.0, 1.0, 1.0), + coordinates_min = (-1.0, -1.0, -1.0), + coordinates_max = (1.0, 1.0, 1.0), periodicity = (true, true, true))::t8_cmesh_t # This is how the cmesh looks like. The numbers are the tree numbers: # +---+ diff --git a/src/solvers/fv_t8code/fv.jl b/src/solvers/fv_t8code/fv.jl index 1e4646c89a3..e531c6ee8a5 100644 --- a/src/solvers/fv_t8code/fv.jl +++ b/src/solvers/fv_t8code/fv.jl @@ -160,7 +160,7 @@ function create_cache(mesh::T8codeMesh, equations::AbstractEquations, solver::FV end function rhs!(du, u, t, mesh::T8codeMesh, equations, - initial_condition, boundary_conditions, source_terms::Source, + boundary_conditions, source_terms::Source, solver::FV, cache) where {Source} # Reset du @trixi_timeit timer() "reset ∂u/∂t" du.=zero(eltype(du)) diff --git a/test/test_mpi_t8code_fv_3d.jl b/test/test_mpi_t8code_fv_3d.jl index 5818b777b5a..ec3316c754b 100644 --- a/test/test_mpi_t8code_fv_3d.jl +++ b/test/test_mpi_t8code_fv_3d.jl @@ -12,11 +12,12 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") # Run basic tests @testset "Examples 3D" begin + # NOTE: Since I use 2x2x2 tree instead of 8x8x8, I need to increase the resolution 2 times by the factor of 2 -> +2 @trixi_testset "elixir_advection_basic.jl" begin @trixi_testset "first-order FV" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), order=1, - initial_refinement_level=2+2, + initial_refinement_level=2 + 2, l2=[0.2848617953369851], linf=[0.3721898718954475]) # Ensure that we do not have excessive memory allocations @@ -30,7 +31,7 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") end @trixi_testset "second-order FV" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - initial_refinement_level=2+2, + initial_refinement_level=2 + 2, l2=[0.10381089565603231], linf=[0.13787405651527007]) # Ensure that we do not have excessive memory allocations @@ -48,7 +49,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") @trixi_testset "elixir_advection_basic_hybrid.jl" begin @trixi_testset "first-order FV" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_basic_hybrid.jl"), order=1, l2=[0.20282363730327146], linf=[0.28132446651281295]) @@ -62,7 +64,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") end end @trixi_testset "second-order FV" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_basic_hybrid.jl"), l2=[0.02153993127089835], linf=[0.039109618097251886]) # Ensure that we do not have excessive memory allocations diff --git a/test/test_t8code_fv_3d.jl b/test/test_t8code_fv_3d.jl index 5f09902198d..e2e1225abde 100644 --- a/test/test_t8code_fv_3d.jl +++ b/test/test_t8code_fv_3d.jl @@ -49,7 +49,7 @@ mkdir(outdir) @trixi_testset "first-order FV" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), order=1, - initial_refinement_level=2+2, + initial_refinement_level=2 + 2, l2=[0.2848617953369851], linf=[0.3721898718954475]) # Ensure that we do not have excessive memory allocations @@ -63,7 +63,7 @@ mkdir(outdir) end @trixi_testset "second-order FV" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - initial_refinement_level=2+2, + initial_refinement_level=2 + 2, l2=[0.10381089565603231], linf=[0.13787405651527007]) # Ensure that we do not have excessive memory allocations @@ -77,7 +77,7 @@ mkdir(outdir) end @trixi_testset "second-order FV, extended reconstruction stencil" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - initial_refinement_level=1+2, + initial_refinement_level=1 + 2, extended_reconstruction_stencil=true, l2=[0.3282177575292713], linf=[0.39002345444858333]) @@ -177,14 +177,14 @@ end 0.0351299673616484, 0.0351299673616484, 0.03512996736164839, - 0.1601847269543808, + 0.1601847269543808 ], linf=[ 0.07175521415072939, 0.04648499338897771, 0.04648499338897816, 0.04648499338897816, - 0.2235470564880404, + 0.2235470564880404 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -203,14 +203,14 @@ end 0.010791416898840429, 0.010791416898840464, 0.010791416898840377, - 0.036995680347196136, + 0.036995680347196136 ], linf=[ 0.01982294164697862, 0.01840725612418126, 0.01840725612418148, 0.01840725612418148, - 0.05736595182767079, + 0.05736595182767079 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -230,14 +230,14 @@ end 0.03596196296013507, 0.03616867188152877, 0.03616867188152873, - 0.14939041550302212, + 0.14939041550302212 ], linf=[ 0.07943789383956079, 0.06389365911606859, 0.06469291944863809, 0.0646929194486372, - 0.23507781748792533, + 0.23507781748792533 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From bb5bceb2cbc8f90b61e5c0653b0d4ea1995e0415 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 15 Oct 2024 14:00:44 +0200 Subject: [PATCH 09/12] Fix parameter in `x_trans_periodic_3d` --- examples/t8code_2d_fv/elixir_advection_gauss.jl | 4 ++-- examples/t8code_2d_fv/elixir_advection_nonperiodic.jl | 2 +- examples/t8code_3d_fv/elixir_advection_gauss.jl | 4 +++- src/equations/linear_scalar_advection_3d.jl | 2 +- test/test_t8code_fv_3d.jl | 8 ++++---- 5 files changed, 11 insertions(+), 9 deletions(-) diff --git a/examples/t8code_2d_fv/elixir_advection_gauss.jl b/examples/t8code_2d_fv/elixir_advection_gauss.jl index 177fce928c6..d4394451df6 100644 --- a/examples/t8code_2d_fv/elixir_advection_gauss.jl +++ b/examples/t8code_2d_fv/elixir_advection_gauss.jl @@ -11,9 +11,9 @@ initial_condition = initial_condition_gauss solver = FV(order = 2, surface_flux = flux_lax_friedrichs) initial_refinement_level = 4 +# Note: The initial_condition is set up for a domain [-5,5]^2. +# This becomes important when using a non-periodic mesh. cmesh = Trixi.cmesh_new_periodic_hybrid() -# Note: A non-periodic run with the tri mesh is unstable. Same as in `elixir_advection_nonperiodic.jl` -# cmesh = Trixi.cmesh_new_tri(periodicity = (false, false)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) diff --git a/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl b/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl index 7e956fe1754..fbe7940617d 100644 --- a/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl +++ b/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl @@ -23,7 +23,7 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs) initial_refinement_level = 5 cmesh = Trixi.cmesh_new_quad(periodicity = (false, false)) -# **Note**: A non-periodic run with the tri mesh is unstable. (Same happens for a non-periodic run with `elixir_advection_gauss.jl`) +# **Note**: A non-periodic run with the tri mesh is unstable. # - This only happens with **2nd order** # - When increasing refinement_level to 6 and lower CFL to 0.4, it seems like the simulation is stable again (except of some smaller noises at the corners) # - With a lower resolution (5) I cannot get the simulation stable. Even with a cfl of 0.01 etc. diff --git a/examples/t8code_3d_fv/elixir_advection_gauss.jl b/examples/t8code_3d_fv/elixir_advection_gauss.jl index 1f9c051f63d..71e6d82eca6 100644 --- a/examples/t8code_3d_fv/elixir_advection_gauss.jl +++ b/examples/t8code_3d_fv/elixir_advection_gauss.jl @@ -13,7 +13,9 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs) initial_refinement_level = 4 # TODO: There are no other cmesh functions implemented yet in 3d. -cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true)) +cmesh = Trixi.cmesh_new_quad_3d(coordinates_min = (-5.0, -5.0, -5.0), + coordinates_max = (5.0, 5.0, 5.0), + periodicity = (true, true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index 006179c1fc9..9c68ce6cd47 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/src/equations/linear_scalar_advection_3d.jl @@ -65,7 +65,7 @@ function initial_condition_convergence_test(x, t, end # Calculates translated coordinates `x` for a periodic domain -function x_trans_periodic_3d(x, domain_length = SVector(2, 2, 2), +function x_trans_periodic_3d(x, domain_length = SVector(10, 10, 10), center = SVector(0, 0, 0)) x_normalized = x .- center x_shifted = x_normalized .% domain_length diff --git a/test/test_t8code_fv_3d.jl b/test/test_t8code_fv_3d.jl index e2e1225abde..b29c09c99b7 100644 --- a/test/test_t8code_fv_3d.jl +++ b/test/test_t8code_fv_3d.jl @@ -97,8 +97,8 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_gauss.jl"), order=1, - l2=[0.1515258539168874], - linf=[0.43164936150417055]) + l2=[0.03422041780251932], + linf=[0.7655163504661142]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -111,8 +111,8 @@ end @trixi_testset "second-order FV" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_gauss.jl"), - l2=[0.04076672839289378], - linf=[0.122537463101035582]) + l2=[0.02129708543319383], + linf=[0.5679262915623222]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 27cb9f603bac0de6e047cdefb2ef9b6c08b9f407 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 15 Oct 2024 15:01:36 +0200 Subject: [PATCH 10/12] Remove comment --- test/test_mpi_t8code_fv_3d.jl | 5 ++--- test/test_t8code_fv_3d.jl | 7 +++---- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/test/test_mpi_t8code_fv_3d.jl b/test/test_mpi_t8code_fv_3d.jl index 5dbd83353b5..575baba084d 100644 --- a/test/test_mpi_t8code_fv_3d.jl +++ b/test/test_mpi_t8code_fv_3d.jl @@ -12,12 +12,11 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") # Run basic tests @testset "Examples 3D" begin - # NOTE: Since I use 2x2x2 tree instead of 8x8x8, I need to increase the resolution 2 times by the factor of 2 -> +2 # @trixi_testset "elixir_advection_basic.jl" begin # @trixi_testset "first-order FV" begin # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), # order=1, - # initial_refinement_level=2 + 2, + # initial_refinement_level=4, # l2=[0.2848617953369851], # linf=[0.3721898718954475]) # # Ensure that we do not have excessive memory allocations @@ -31,7 +30,7 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") # end # @trixi_testset "second-order FV" begin # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - # initial_refinement_level=2 + 2, + # initial_refinement_level=4, # l2=[0.10381089565603231], # linf=[0.13787405651527007]) # # Ensure that we do not have excessive memory allocations diff --git a/test/test_t8code_fv_3d.jl b/test/test_t8code_fv_3d.jl index 1a8c5da9e63..18a5dd2cb05 100644 --- a/test/test_t8code_fv_3d.jl +++ b/test/test_t8code_fv_3d.jl @@ -44,12 +44,11 @@ mkdir(outdir) # end # end -# NOTE: Since I use 2x2x2 tree instead of 8x8x8, I need to increase the resolution 2 times by the factor of 2 -> +2 # @trixi_testset "elixir_advection_basic.jl" begin # @trixi_testset "first-order FV" begin # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), # order=1, -# initial_refinement_level=2 + 2, +# initial_refinement_level=4, # l2=[0.2848617953369851], # linf=[0.3721898718954475]) # # Ensure that we do not have excessive memory allocations @@ -63,7 +62,7 @@ mkdir(outdir) # end # @trixi_testset "second-order FV" begin # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# initial_refinement_level=2 + 2, +# initial_refinement_level=4, # l2=[0.10381089565603231], # linf=[0.13787405651527007]) # # Ensure that we do not have excessive memory allocations @@ -77,7 +76,7 @@ mkdir(outdir) # end # @trixi_testset "second-order FV, extended reconstruction stencil" begin # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# initial_refinement_level=1 + 2, +# initial_refinement_level=3, # extended_reconstruction_stencil=true, # l2=[0.3282177575292713], # linf=[0.39002345444858333]) From ffe857d2bd00fb84d052f609fd711a5ad3d332b5 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 15 Oct 2024 15:42:51 +0200 Subject: [PATCH 11/12] Some last changes --- src/meshes/t8code_mesh.jl | 50 ------------------- test/test_mpi_t8code_fv_2d.jl | 62 +++++++++++------------ test/test_mpi_t8code_fv_3d.jl | 66 ++++++++++++------------- test/test_t8code_fv_2d.jl | 93 +++++++++++++++++++---------------- test/test_t8code_fv_3d.jl | 92 +++++++++++++++++----------------- 5 files changed, 160 insertions(+), 203 deletions(-) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 164a816aa15..47d0eb0304f 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -1886,56 +1886,6 @@ function cmesh_new_tri(; trees_per_dimension = (1, 1), return cmesh end -# # Old version if triangular mesh -# function cmesh_new_periodic_tri(; periodicity = (true, true), comm = mpi_comm())::t8_cmesh_t -# n_dims = 2 -# vertices = [ # Just all vertices of all trees. partly duplicated -# -1.0, -1.0, 0, # tree 0, triangle -# 1.0, -1.0, 0, -# 1.0, 1.0, 0, -# -1.0, -1.0, 0, # tree 1, triangle -# 1.0, 1.0, 0, -# -1.0, 1.0, 0, -# ] - -# # Generally, one can define other geometries. But besides linear the other -# # geometries in t8code do not have C interface yet. -# linear_geom = t8_geometry_linear_new(n_dims) - -# # This is how the cmesh looks like. The numbers are the tree numbers: -# # -# # +---+ -# # |1 /| -# # | / | -# # |/0 | -# # +---+ -# # - -# cmesh_ref = Ref(t8_cmesh_t()) -# t8_cmesh_init(cmesh_ref) -# cmesh = cmesh_ref[] - -# # Use linear geometry -# t8_cmesh_register_geometry(cmesh, linear_geom) -# t8_cmesh_set_tree_class(cmesh, 0, T8_ECLASS_TRIANGLE) -# t8_cmesh_set_tree_class(cmesh, 1, T8_ECLASS_TRIANGLE) - -# t8_cmesh_set_tree_vertices(cmesh, 0, @views(vertices[(1 + 0):end]), 3) -# t8_cmesh_set_tree_vertices(cmesh, 1, @views(vertices[(1 + 9):end]), 3) - -# t8_cmesh_set_join(cmesh, 0, 1, 1, 2, 0) -# if periodicity[1] -# t8_cmesh_set_join(cmesh, 0, 1, 0, 1, 0) -# end -# if periodicity[2] -# t8_cmesh_set_join(cmesh, 0, 1, 2, 0, 1) -# end - -# t8_cmesh_commit(cmesh, comm) - -# return cmesh -# end - function cmesh_new_periodic_tri2(; comm = mpi_comm())::t8_cmesh_t n_dims = 2 vertices = [ # Just all vertices of all trees. partly duplicated diff --git a/test/test_mpi_t8code_fv_2d.jl b/test/test_mpi_t8code_fv_2d.jl index 14d918794d2..54b4cd83333 100644 --- a/test/test_mpi_t8code_fv_2d.jl +++ b/test/test_mpi_t8code_fv_2d.jl @@ -13,37 +13,37 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") # Run basic tests @testset "Examples 2D" begin # Linear scalar advection - # @trixi_testset "elixir_advection_basic.jl" begin - # @trixi_testset "first-order FV" begin - # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - # order=1, - # l2=[0.08551397247817498], - # linf=[0.12087467695430498]) - # # Ensure that we do not have excessive memory allocations - # # (e.g., from type instabilities) - # let - # t = sol.t[end] - # u_ode = sol.u[end] - # du_ode = similar(u_ode) - # @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - # end - # end - # @trixi_testset "second-order FV" begin - # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - # l2=[0.008142380494734171], - # linf=[0.018687916234976898]) - # # Ensure that we do not have excessive memory allocations - # # (e.g., from type instabilities) - # let - # t = sol.t[end] - # u_ode = sol.u[end] - # du_ode = similar(u_ode) - # @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - # end - # end - # # The extended reconstruction stencil is currently not mpi parallel. - # # The current version runs through but an error occurs on some rank. - # end + @trixi_testset "elixir_advection_basic.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + order=1, + l2=[0.08551397247817498], + linf=[0.12087467695430498]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + l2=[0.008142380494734171], + linf=[0.018687916234976898]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + # The extended reconstruction stencil is currently not mpi parallel. + # The current version runs through but an error occurs on some rank. + end @trixi_testset "elixir_advection_gauss.jl" begin @trixi_testset "first-order FV" begin diff --git a/test/test_mpi_t8code_fv_3d.jl b/test/test_mpi_t8code_fv_3d.jl index 575baba084d..89faeb6ac57 100644 --- a/test/test_mpi_t8code_fv_3d.jl +++ b/test/test_mpi_t8code_fv_3d.jl @@ -12,39 +12,39 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_fv") # Run basic tests @testset "Examples 3D" begin - # @trixi_testset "elixir_advection_basic.jl" begin - # @trixi_testset "first-order FV" begin - # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - # order=1, - # initial_refinement_level=4, - # l2=[0.2848617953369851], - # linf=[0.3721898718954475]) - # # Ensure that we do not have excessive memory allocations - # # (e.g., from type instabilities) - # let - # t = sol.t[end] - # u_ode = sol.u[end] - # du_ode = similar(u_ode) - # @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - # end - # end - # @trixi_testset "second-order FV" begin - # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - # initial_refinement_level=4, - # l2=[0.10381089565603231], - # linf=[0.13787405651527007]) - # # Ensure that we do not have excessive memory allocations - # # (e.g., from type instabilities) - # let - # t = sol.t[end] - # u_ode = sol.u[end] - # du_ode = similar(u_ode) - # @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - # end - # end - # # The extended reconstruction stencil is currently not mpi parallel. - # # The current version runs through but an error occurs on some rank. - # end + @trixi_testset "elixir_advection_basic.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + order=1, + initial_refinement_level=4, + l2=[0.2848617953369851], + linf=[0.3721898718954475]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + initial_refinement_level=4, + l2=[0.10381089565603231], + linf=[0.13787405651527007]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + # The extended reconstruction stencil is currently not mpi parallel. + # The current version runs through but an error occurs on some rank. + end @trixi_testset "elixir_advection_basic_hybrid.jl" begin @trixi_testset "first-order FV" begin diff --git a/test/test_t8code_fv_2d.jl b/test/test_t8code_fv_2d.jl index 46330d2d98e..22eb5fd49dd 100644 --- a/test/test_t8code_fv_2d.jl +++ b/test/test_t8code_fv_2d.jl @@ -14,6 +14,13 @@ outdir = "out" isdir(outdir) && rm(outdir, recursive = true) mkdir(outdir) +# The are some reasons why the FV tests fail on the CI runs: +# - T8code.jl is no dependency of the Trixi test environment. So `using T8code` within the `advection_basic` elixir throws an error. +# (But: When removing this line, e.g. code like `trixi_t8_mapping_c` (with unknown types) or unknown eclasses failes) +# - `load_tree_data = @t8_load_tree_data(t8_geom_load_tree_data_vertices)` is required. +# Sadly, the macro is not yet in the latest T8code.jl release (but for instance in `jmark/bumb-t8code-3.0.0` or `bennibolm/t8-trixi-fv-scheme`). + + @testset "T8codeMesh2D" begin #! format: noindent @@ -44,49 +51,49 @@ mkdir(outdir) # end # end -# @trixi_testset "elixir_advection_basic.jl" begin -# @trixi_testset "first-order FV" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# order=1, -# l2=[0.08551397247817498], -# linf=[0.12087467695430498]) -# # Ensure that we do not have excessive memory allocations -# # (e.g., from type instabilities) -# let -# t = sol.t[end] -# u_ode = sol.u[end] -# du_ode = similar(u_ode) -# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 -# end -# end -# @trixi_testset "second-order FV" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# l2=[0.008142380494734171], -# linf=[0.018687916234976898]) -# # Ensure that we do not have excessive memory allocations -# # (e.g., from type instabilities) -# let -# t = sol.t[end] -# u_ode = sol.u[end] -# du_ode = similar(u_ode) -# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 -# end -# end -# @trixi_testset "second-order FV, extended reconstruction stencil" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# extended_reconstruction_stencil=true, -# l2=[0.028436326251639936], -# linf=[0.08696815845435057]) -# # Ensure that we do not have excessive memory allocations -# # (e.g., from type instabilities) -# let -# t = sol.t[end] -# u_ode = sol.u[end] -# du_ode = similar(u_ode) -# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 -# end -# end -# end +@trixi_testset "elixir_advection_basic.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + order=1, + l2=[0.08551397247817498], + linf=[0.12087467695430498]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + l2=[0.008142380494734171], + linf=[0.018687916234976898]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV, extended reconstruction stencil" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + extended_reconstruction_stencil=true, + l2=[0.028436326251639936], + linf=[0.08696815845435057]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end @trixi_testset "elixir_advection_gauss.jl" begin @trixi_testset "first-order FV" begin diff --git a/test/test_t8code_fv_3d.jl b/test/test_t8code_fv_3d.jl index 18a5dd2cb05..a3150a3ab29 100644 --- a/test/test_t8code_fv_3d.jl +++ b/test/test_t8code_fv_3d.jl @@ -44,52 +44,52 @@ mkdir(outdir) # end # end -# @trixi_testset "elixir_advection_basic.jl" begin -# @trixi_testset "first-order FV" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# order=1, -# initial_refinement_level=4, -# l2=[0.2848617953369851], -# linf=[0.3721898718954475]) -# # Ensure that we do not have excessive memory allocations -# # (e.g., from type instabilities) -# let -# t = sol.t[end] -# u_ode = sol.u[end] -# du_ode = similar(u_ode) -# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 -# end -# end -# @trixi_testset "second-order FV" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# initial_refinement_level=4, -# l2=[0.10381089565603231], -# linf=[0.13787405651527007]) -# # Ensure that we do not have excessive memory allocations -# # (e.g., from type instabilities) -# let -# t = sol.t[end] -# u_ode = sol.u[end] -# du_ode = similar(u_ode) -# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 -# end -# end -# @trixi_testset "second-order FV, extended reconstruction stencil" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), -# initial_refinement_level=3, -# extended_reconstruction_stencil=true, -# l2=[0.3282177575292713], -# linf=[0.39002345444858333]) -# # Ensure that we do not have excessive memory allocations -# # (e.g., from type instabilities) -# let -# t = sol.t[end] -# u_ode = sol.u[end] -# du_ode = similar(u_ode) -# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 -# end -# end -# end +@trixi_testset "elixir_advection_basic.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + order=1, + initial_refinement_level=4, + l2=[0.2848617953369851], + linf=[0.3721898718954475]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + initial_refinement_level=4, + l2=[0.10381089565603231], + linf=[0.13787405651527007]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "second-order FV, extended reconstruction stencil" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + initial_refinement_level=3, + extended_reconstruction_stencil=true, + l2=[0.3282177575292713], + linf=[0.39002345444858333]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end @trixi_testset "elixir_advection_gauss.jl" begin @trixi_testset "first-order FV" begin From 59401502f99b13726ed7cc4b55eebd72546fbb34 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 15 Oct 2024 15:44:30 +0200 Subject: [PATCH 12/12] fmt --- test/test_t8code_fv_2d.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_t8code_fv_2d.jl b/test/test_t8code_fv_2d.jl index 22eb5fd49dd..47f607252ad 100644 --- a/test/test_t8code_fv_2d.jl +++ b/test/test_t8code_fv_2d.jl @@ -16,11 +16,10 @@ mkdir(outdir) # The are some reasons why the FV tests fail on the CI runs: # - T8code.jl is no dependency of the Trixi test environment. So `using T8code` within the `advection_basic` elixir throws an error. -# (But: When removing this line, e.g. code like `trixi_t8_mapping_c` (with unknown types) or unknown eclasses failes) +# (But: When removing this line, e.g. code like `trixi_t8_mapping_c` (with unknown types) or unknown eclasses fails) # - `load_tree_data = @t8_load_tree_data(t8_geom_load_tree_data_vertices)` is required. # Sadly, the macro is not yet in the latest T8code.jl release (but for instance in `jmark/bumb-t8code-3.0.0` or `bennibolm/t8-trixi-fv-scheme`). - @testset "T8codeMesh2D" begin #! format: noindent