diff --git a/Project.toml b/Project.toml index 612b65ec6f0..012b64bcf35 100644 --- a/Project.toml +++ b/Project.toml @@ -103,7 +103,7 @@ StaticArrays = "1.5" StrideArrays = "0.1.26" StructArrays = "0.6.11" SummationByPartsOperators = "0.5.41" -T8code = "0.5" +T8code = "0.7" TimerOutputs = "0.5.7" Triangulate = "2.2" TriplotBase = "0.1" diff --git a/examples/t8code_2d_fv/elixir_advection_basic.jl b/examples/t8code_2d_fv/elixir_advection_basic.jl index fde4e3d16f8..c26c6dc85db 100644 --- a/examples/t8code_2d_fv/elixir_advection_basic.jl +++ b/examples/t8code_2d_fv/elixir_advection_basic.jl @@ -1,5 +1,6 @@ using OrdinaryDiffEq using Trixi +using T8code advection_velocity = (0.2, -0.7) equations = LinearScalarAdvectionEquation2D(advection_velocity) @@ -9,14 +10,73 @@ 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) # minimum coordinates (min(x), min(y)) -coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y)) - -trees_per_dimension = (8, 8) - -mesh = T8codeMesh(trees_per_dimension, - coordinates_min = coordinates_min, coordinates_max = coordinates_max, - initial_refinement_level = 4) +# Note: +# For now, it is completely irrelevant that coordinates_max/min are. +# The used t8code routine creates the mesh on [0, nx] x [0, ny], where (nx, ny) = trees_per_dimension. +# Afterwards and only inside Trixi, `tree_node_coordinates` are mapped back to [-1, 1]^2. +# But, this variable is not used for the FV method. +# That's why I use the cmesh interface in all other elixirs. + +# Option 1: coordinates +coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y)) +coordinates_max = (8.0, 8.0) # maximum coordinates (max(x), max(y)) +# Note and TODO: The plan is to move the auxiliary routine f and the macro to a different place. +# Then, somehow, I get SegFaults when using this `mapping_coordinates` or (equally) when +# using `coordinates_min/max` and then use the `coordinates2mapping` within the constructor. +# With both other mappings I don't get that. +mapping_coordinates = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +# Option 2: faces +f1(s) = SVector(-1.0, s - 1.0) +f2(s) = SVector(1.0, s + 1.0) +f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s)) +f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s)) +faces = (f1, f2, f3, f4) +mapping_faces = Trixi.transfinite_mapping(faces) + +# Option 3: classic mapping +function mapping(xi, eta) + x = 4 * (xi + 1) + y = 4 * (eta + 1) + + return SVector(x, y) +end + +# Note and TODO: +# Normally, this should be put somewhere else. For now, that doesn't properly. +# See note in `src/auxiliary/t8code.jl` +function f(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, + ref_coords, num_coords, out_coords) + + for i in 1:num_coords + offset_3d = 3 * (i - 1) + 1 + + xi = unsafe_load(out_coords, offset_3d) + eta = unsafe_load(out_coords, offset_3d + 1) + # xy = mapping_coordinates(xi, eta) + # xy = mapping_faces(xi, eta) + xy = mapping(xi, eta) + + unsafe_store!(out_coords, xy[1], offset_3d) + unsafe_store!(out_coords, xy[2], offset_3d + 1) + end + + return nothing +end + +trees_per_dimension = (2, 2) + +eclass = T8_ECLASS_QUAD +mesh = T8codeMesh(trees_per_dimension, eclass, + mapping = @t8_analytic_callback(f), + # Plan is to use either + # coordinates_max = coordinates_max, coordinates_min = coordinates_min, + # or at least + # mapping = mapping, + initial_refinement_level = 6) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_2d_fv/elixir_advection_basic_hybrid.jl index e51d4f86623..41470fefdc3 100644 --- a/examples/t8code_2d_fv/elixir_advection_basic_hybrid.jl +++ b/examples/t8code_2d_fv/elixir_advection_basic_hybrid.jl @@ -10,8 +10,8 @@ solver = FV(order = 2, extended_reconstruction_stencil = false, surface_flux = flux_lax_friedrichs) cmesh = Trixi.cmesh_new_periodic_hybrid() -# cmesh = Trixi.cmesh_quad(periodicity = (true, true)) -# cmesh = Trixi.cmesh_new_periodic_tri() +# cmesh = Trixi.cmesh_new_quad(periodicity = (true, true)) +# cmesh = Trixi.cmesh_new_tri(periodicity = (true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = 3) diff --git a/examples/t8code_2d_fv/elixir_advection_gauss.jl b/examples/t8code_2d_fv/elixir_advection_gauss.jl index 14c75953f6c..177fce928c6 100644 --- a/examples/t8code_2d_fv/elixir_advection_gauss.jl +++ b/examples/t8code_2d_fv/elixir_advection_gauss.jl @@ -12,6 +12,9 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs) initial_refinement_level = 4 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) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl b/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl index 23748694991..7e956fe1754 100644 --- a/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl +++ b/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl @@ -21,7 +21,16 @@ 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(periodicity = (false, false)) +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`) +# - 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. +# -> That can't be expected. +# -> Maybe, the reconstruction is just not fitted for this example/mesh/resolution?! +# cmesh = Trixi.cmesh_new_tri(periodicity = (false, false)) + mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, @@ -39,10 +48,13 @@ 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.8) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - stepsize_callback) + save_solution, stepsize_callback) ############################################################################### # Run the simulation. diff --git a/examples/t8code_2d_fv/elixir_euler_blast_wave.jl b/examples/t8code_2d_fv/elixir_euler_blast_wave.jl index 7779448420f..93277d97838 100644 --- a/examples/t8code_2d_fv/elixir_euler_blast_wave.jl +++ b/examples/t8code_2d_fv/elixir_euler_blast_wave.jl @@ -30,7 +30,7 @@ solver = FV(order = 2, extended_reconstruction_stencil = true, initial_refinement_level = 4 # cmesh = Trixi.cmesh_new_periodic_hybrid2() -cmesh = Trixi.cmesh_quad(periodicity = (true, true)) +cmesh = Trixi.cmesh_new_quad(periodicity = (true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_2d_fv/elixir_euler_free_stream.jl b/examples/t8code_2d_fv/elixir_euler_free_stream.jl index 3c0ac37ef68..b6abb8aa5f8 100644 --- a/examples/t8code_2d_fv/elixir_euler_free_stream.jl +++ b/examples/t8code_2d_fv/elixir_euler_free_stream.jl @@ -32,7 +32,8 @@ end mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", joinpath(@__DIR__, "square_unstructured_1.inp")) -# TODO: Including mapping does not work yet for FV. +# Note: +# Including mapping does not work yet for FV. mesh = T8codeMesh(mesh_file, 2; polydeg = 0, mapping = mapping, initial_refinement_level = 2) diff --git a/examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl b/examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl index 659e55701b6..7adf83fc47f 100644 --- a/examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl +++ b/examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl @@ -36,8 +36,8 @@ solver = FV(order = 2, extended_reconstruction_stencil = false, initial_refinement_level = 4 # cmesh = Trixi.cmesh_new_periodic_hybrid() -# cmesh = Trixi.cmesh_quad(periodicity = (true, true)) -# cmesh = Trixi.cmesh_new_periodic_tri() +# cmesh = Trixi.cmesh_new_quad(periodicity = (true, true)) +# cmesh = Trixi.cmesh_new_tri(periodicity = (true, true)) cmesh = Trixi.cmesh_new_periodic_tri2() mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) diff --git a/examples/t8code_2d_fv/elixir_euler_source_terms.jl b/examples/t8code_2d_fv/elixir_euler_source_terms.jl index 63f22b05a36..74cc1aae7ee 100644 --- a/examples/t8code_2d_fv/elixir_euler_source_terms.jl +++ b/examples/t8code_2d_fv/elixir_euler_source_terms.jl @@ -15,8 +15,8 @@ solver = FV(order = 2, extended_reconstruction_stencil = false, surface_flux = flux_lax_friedrichs) cmesh = Trixi.cmesh_new_periodic_hybrid() -# cmesh = Trixi.cmesh_quad(periodicity = (true, true)) -# cmesh = Trixi.cmesh_new_periodic_tri() +# cmesh = Trixi.cmesh_new_quad(periodicity = (true, true)) +# cmesh = Trixi.cmesh_new_tri(periodicity = (true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = 3) diff --git a/examples/t8code_3d_fv/elixir_advection_basic.jl b/examples/t8code_3d_fv/elixir_advection_basic.jl index 071ef23dc8d..d5937ae7cbd 100644 --- a/examples/t8code_3d_fv/elixir_advection_basic.jl +++ b/examples/t8code_3d_fv/elixir_advection_basic.jl @@ -1,5 +1,6 @@ using OrdinaryDiffEq using Trixi +using T8code advection_velocity = (0.2, -0.7, 0.5) equations = LinearScalarAdvectionEquation3D(advection_velocity) @@ -9,14 +10,55 @@ 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) +# Option 1: coordinates +# For all problems see the 2D file... +coordinates_min = (0.0, 0.0, 0.0) # minimum coordinates (min(x), min(y), min(z)) +coordinates_max = (8.0, 8.0, 8.0) # maximum coordinates (max(x), max(y), max(z)) + +mapping_coordinates = Trixi.coordinates2mapping(coordinates_min, coordinates_max) + +# Option 3: classic mapping +function mapping(xi, eta, zeta) + x = 4 * (xi + 1) + y = 4 * (eta + 1) + z = 4 * (zeta + 1) + + return SVector(x, y, z) +end + +function f(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, + ref_coords, num_coords, out_coords) + + for i in 1:num_coords + offset_3d = 3 * (i - 1) + 1 + + xi = unsafe_load(out_coords, offset_3d) + eta = unsafe_load(out_coords, offset_3d + 1) + zeta = unsafe_load(out_coords, offset_3d + 2) + # xyz = mapping_coordinates(xi, eta, zeta) + xyz = mapping(xi, eta, zeta) + + unsafe_store!(out_coords, xyz[1], offset_3d) + unsafe_store!(out_coords, xyz[2], offset_3d + 1) + unsafe_store!(out_coords, xyz[3], offset_3d + 2) + end + + return nothing +end + +trees_per_dimension = (1, 1, 1) + +eclass = T8_ECLASS_HEX +mesh = T8codeMesh(trees_per_dimension, eclass; + mapping = @t8_analytic_callback(f), + # Plan is to use either + # coordinates_max = coordinates_max, coordinates_min = coordinates_min, + # or at least + # mapping = mapping, + initial_refinement_level = 6) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl index f07ae5bed3d..8d6dddca1e0 100644 --- a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl +++ b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl @@ -10,7 +10,7 @@ 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)) +cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = 4) diff --git a/examples/t8code_3d_fv/elixir_advection_gauss.jl b/examples/t8code_3d_fv/elixir_advection_gauss.jl index ef223dd6aeb..1f9c051f63d 100644 --- a/examples/t8code_3d_fv/elixir_advection_gauss.jl +++ b/examples/t8code_3d_fv/elixir_advection_gauss.jl @@ -13,7 +13,7 @@ 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_quad_3d(periodicity = (true, true, true)) +cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) diff --git a/examples/t8code_3d_fv/elixir_advection_nonperiodic.jl b/examples/t8code_3d_fv/elixir_advection_nonperiodic.jl index 22ba4037898..414fcbe4a8f 100644 --- a/examples/t8code_3d_fv/elixir_advection_nonperiodic.jl +++ b/examples/t8code_3d_fv/elixir_advection_nonperiodic.jl @@ -21,7 +21,7 @@ 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)) +cmesh = Trixi.cmesh_new_quad_3d(periodicity = (false, false, false)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level) semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, diff --git a/examples/t8code_3d_fv/elixir_euler_source_terms.jl b/examples/t8code_3d_fv/elixir_euler_source_terms.jl index 8568a1dac96..e588e24674f 100644 --- a/examples/t8code_3d_fv/elixir_euler_source_terms.jl +++ b/examples/t8code_3d_fv/elixir_euler_source_terms.jl @@ -15,7 +15,7 @@ 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)) +cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true)) mesh = T8codeMesh(cmesh, solver, initial_refinement_level = 3) diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl index 7c1399fc803..8327a618734 100644 --- a/src/auxiliary/t8code.jl +++ b/src/auxiliary/t8code.jl @@ -190,3 +190,34 @@ function trixi_t8_adapt!(mesh, indicators) return differences end + +# Note and TODO: +# This routine seems to work for most of the mappings. +# Somehow, when using `coordinates_min/max` and then `coordinates2mapping`, +# I get SegFault errors. +# This happens when running this in a new julia session or after some runs when I ran simulations +# with other mappings before. +# Cannot figure out why. For now, I will leave this auxiliary mapping within the elixir +# and comment this one out. +# function trixi_t8_mapping_c(mapping) +# function f(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, ref_coords, num_coords, out_coords) + +# for i in 1:num_coords +# offset_3d = 3 * (i - 1) + 1 + +# xi = unsafe_load(out_coords, offset_3d) +# eta = unsafe_load(out_coords, offset_3d + 1) +# xy = mapping(xi, eta) + +# unsafe_store!(out_coords, xy[1], offset_3d) +# unsafe_store!(out_coords, xy[2], offset_3d + 1) +# end + +# return nothing +# end + +# return @cfunction($f, Cvoid, (t8_cmesh_t, t8_gloidx_t, Ptr{Cdouble}, Csize_t, Ptr{Cdouble}, Ptr{Cvoid}, Ptr{Cvoid})) +# end diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index aa83fceaafb..7cb2f7e86b4 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -348,6 +348,264 @@ function T8codeMesh(trees_per_dimension; polydeg = 1, mapping = mapping_) end +# New interface for FV solver +# Parameter eclass decides which class of element are used. +function T8codeMesh(trees_per_dimension, eclass; + mapping = nothing, faces = nothing, coordinates_min = nothing, + coordinates_max = nothing, + RealT = Float64, initial_refinement_level = 0, + periodicity = true) + @assert ((coordinates_min === nothing)===(coordinates_max === nothing)) "Either both or none of coordinates_min and coordinates_max must be specified" + + @assert count(i -> i !== nothing, + (mapping, faces, coordinates_min))==1 "Exactly one of mapping, faces and coordinates_min/max must be specified" + + # Extract mapping + if faces !== nothing + validate_faces(faces) + mapping = transfinite_mapping(faces) + elseif coordinates_min !== nothing + mapping = coordinates2mapping(coordinates_min, coordinates_max) + end + + NDIMS = length(trees_per_dimension) + @assert (NDIMS == 2||NDIMS == 3) "NDIMS should be 2 or 3." + + # Convert periodicity to a Tuple of a Bool for every dimension + if all(periodicity) + # Also catches case where periodicity = true + periodicity = ntuple(_ -> true, NDIMS) + elseif !any(periodicity) + # Also catches case where periodicity = false + periodicity = ntuple(_ -> false, NDIMS) + else + # Default case if periodicity is an iterable + periodicity = Tuple(periodicity) + end + + do_partition = 0 + if NDIMS == 2 + 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) + 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 + vertices[index] = coordinates_tree_x[itree_x] + vertices[index + 1] = coordinates_tree_y[itree_y] + vertices[index + 2] = 0.0 + + index += 3 + vertices[index] = coordinates_tree_x[itree_x + 1] + vertices[index + 1] = coordinates_tree_y[itree_y] + vertices[index + 2] = 0.0 + + index += 3 + vertices[index] = coordinates_tree_x[itree_x] + vertices[index + 1] = coordinates_tree_y[itree_y + 1] + vertices[index + 2] = 0.0 + + index += 3 + vertices[index] = coordinates_tree_x[itree_x + 1] + 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[2] * 3 * 8 * (itree_z - 1) + + trees_per_dimension[1] * 3 * 8 * (itree_y - 1) + + 3 * 8 * (itree_x - 1) + 1 + + vertices[index] = coordinates_tree_x[itree_x] + vertices[index + 1] = coordinates_tree_y[itree_y] + vertices[index + 2] = coordinates_tree_z[itree_z] + + index += 3 + vertices[index] = coordinates_tree_x[itree_x + 1] + vertices[index + 1] = coordinates_tree_y[itree_y] + vertices[index + 2] = coordinates_tree_z[itree_z] + + index += 3 + vertices[index] = coordinates_tree_x[itree_x] + vertices[index + 1] = coordinates_tree_y[itree_y + 1] + vertices[index + 2] = coordinates_tree_z[itree_z] + + index += 3 + vertices[index] = coordinates_tree_x[itree_x + 1] + vertices[index + 1] = coordinates_tree_y[itree_y + 1] + vertices[index + 2] = coordinates_tree_z[itree_z] + + # itree_z + 1 + index += 3 + vertices[index] = coordinates_tree_x[itree_x] + vertices[index + 1] = coordinates_tree_y[itree_y] + vertices[index + 2] = coordinates_tree_z[itree_z + 1] + + index += 3 + vertices[index] = coordinates_tree_x[itree_x + 1] + vertices[index + 1] = coordinates_tree_y[itree_y] + vertices[index + 2] = coordinates_tree_z[itree_z + 1] + + index += 3 + vertices[index] = coordinates_tree_x[itree_x] + vertices[index + 1] = coordinates_tree_y[itree_y + 1] + vertices[index + 2] = coordinates_tree_z[itree_z + 1] + + index += 3 + vertices[index] = coordinates_tree_x[itree_x + 1] + vertices[index + 1] = coordinates_tree_y[itree_y + 1] + vertices[index + 2] = coordinates_tree_z[itree_z + 1] + 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 + + # 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 + if periodicity[3] + t8_cmesh_set_join(cmesh, 0, 0, 4, 5, 0) + end + elseif num_trees == 8 + error("Not supported") + t8_cmesh_set_join_by_stash(cmesh, C_NULL, 0) + else + error("Not supported trees_per_dimension") + end + + t8_cmesh_commit(cmesh, mpi_comm()) + end + + do_face_ghost = mpi_isparallel() + scheme = t8_scheme_new_default_cxx() + forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost, + mpi_comm()) + + # Non-periodic boundaries. + boundary_names = fill(Symbol("---"), 2 * NDIMS, prod(trees_per_dimension)) + + for itree in 1:t8_forest_get_num_global_trees(forest) + if !periodicity[1] + boundary_names[1, itree] = :x_neg + boundary_names[2, itree] = :x_pos + end + + if !periodicity[2] + boundary_names[3, itree] = :y_neg + boundary_names[4, itree] = :y_pos + end + + if NDIMS > 2 + if !periodicity[3] + boundary_names[5, itree] = :z_neg + boundary_names[6, itree] = :z_pos + end + end + end + + return T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = 1) +end + """ T8codeMesh(cmesh::Ptr{t8_cmesh}, mapping=nothing, polydeg=1, RealT=Float64, @@ -1473,7 +1731,7 @@ function cmesh_new_periodic_hybrid(; comm = mpi_comm())::t8_cmesh_t 1.0, 1.0, 0, 0, 0, 0, # tree 5, triangle 1.0, 1.0, 0, - 0, 1.0, 0, + 0, 1.0, 0 ] # Generally, one can define other geometries. But besides linear the other @@ -1534,25 +1792,9 @@ function cmesh_new_periodic_hybrid(; comm = mpi_comm())::t8_cmesh_t return cmesh end -function cmesh_quad(; comm = mpi_comm(), periodicity = (true, true))::t8_cmesh_t - n_dims = 2 - vertices = [ # Just all vertices of all trees. partly duplicated - -1.0, -1.0, 0, # tree 0, quad - 1.0, -1.0, 0, - -1.0, 1.0, 0, - 1.0, 1.0, 0, - - # rotated: - # -1.0, 0.0, 0, # tree 0, quad - # 0.0, -1.0, 0, - # 0.0, 1.0, 0, - # 1.0, 0.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) - +function cmesh_new_quad(; trees_per_dimension = (1, 1), + coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0), + periodicity = (true, true))::t8_cmesh_t # This is how the cmesh looks like. The numbers are the tree numbers: # # +---+ @@ -1562,45 +1804,29 @@ function cmesh_quad(; comm = mpi_comm(), periodicity = (true, true))::t8_cmesh_t # +---+ # - 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_QUAD) - - t8_cmesh_set_tree_vertices(cmesh, 0, @views(vertices[(1 + 0):end]), 4) - - 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 - - t8_cmesh_commit(cmesh, comm) + polygons_x, polygons_y = trees_per_dimension + periodic_x, periodic_y = periodicity + boundary = [coordinates_min[1], coordinates_min[2], 0.0, + coordinates_max[1], coordinates_min[2], 0.0, + coordinates_min[1], coordinates_max[2], 0.0, + coordinates_max[1], coordinates_max[2], 0.0] + + comm = mpi_comm() + set_partition = 0 + offset = 0.0 + use_axis_aligned = 0 + eclass = T8_ECLASS_QUAD + cmesh = t8_cmesh_new_hypercube_pad_ext(eclass, comm, boundary, + polygons_x, polygons_y, 0, + periodic_x, periodic_y, 0, + use_axis_aligned, set_partition, offset) 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) - +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), + periodicity = (true, true, true))::t8_cmesh_t # This is how the cmesh looks like. The numbers are the tree numbers: # +---+ # / /| @@ -1611,54 +1837,35 @@ function cmesh_quad_3d(; comm = mpi_comm(), periodicity = (true, true, true))::t # +---+ # - 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) + polygons_x, polygons_y, polygons_z = trees_per_dimension + periodic_x, periodic_y, periodic_z = periodicity + boundary = [coordinates_min[1], coordinates_min[2], coordinates_min[3], + coordinates_max[1], coordinates_min[2], coordinates_min[3], + coordinates_min[1], coordinates_max[2], coordinates_min[3], + coordinates_max[1], coordinates_max[2], coordinates_min[3], + coordinates_min[1], coordinates_min[2], coordinates_max[3], + coordinates_max[1], coordinates_min[2], coordinates_max[3], + coordinates_min[1], coordinates_max[2], coordinates_max[3], + coordinates_max[1], coordinates_max[2], coordinates_max[3]] + + comm = mpi_comm() + set_partition = 0 + offset = 0.0 + use_axis_aligned = 0 + eclass = T8_ECLASS_HEX + cmesh = t8_cmesh_new_hypercube_pad_ext(eclass, comm, boundary, + polygons_x, polygons_y, polygons_z, + periodic_x, periodic_y, periodic_z, + use_axis_aligned, set_partition, offset) 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 - -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, - - # rotated: - # -1.0, 0, 0, # tree 0, triangle - # 0, -1.0, 0, - # 1.0, 0, 0, - # -1.0, 0, 0, # tree 1, triangle - # 1.0, 0, 0, - # 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) - +# TODO: The structure of `cmesh_new_quad` and `cmesh_new_tri` is equal and only differs for `eclass`. +# Use only one routine! +function cmesh_new_tri(; trees_per_dimension = (1, 1), + coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0), + periodicity = (true, true))::t8_cmesh_t # This is how the cmesh looks like. The numbers are the tree numbers: # # +---+ @@ -1668,27 +1875,79 @@ function cmesh_new_periodic_tri(; comm = mpi_comm())::t8_cmesh_t # +---+ # - 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) - t8_cmesh_set_join(cmesh, 0, 1, 0, 1, 0) - t8_cmesh_set_join(cmesh, 0, 1, 2, 0, 1) - - t8_cmesh_commit(cmesh, comm) + # Note: + # Number of trees = prod(trees_per_dimension) * 2 + # since each quad is divided into 2 triangles + polygons_x, polygons_y = trees_per_dimension + periodic_x, periodic_y = periodicity + boundary = [coordinates_min[1], coordinates_min[2], 0.0, + coordinates_max[1], coordinates_min[2], 0.0, + coordinates_min[1], coordinates_max[2], 0.0, + coordinates_max[1], coordinates_max[2], 0.0] + + comm = mpi_comm() + set_partition = 0 + offset = 0.0 + use_axis_aligned = 0 + eclass = T8_ECLASS_TRIANGLE + cmesh = t8_cmesh_new_hypercube_pad_ext(eclass, comm, boundary, + polygons_x, polygons_y, 0, + periodic_x, periodic_y, 0, + use_axis_aligned, set_partition, offset) 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 @@ -1715,7 +1974,7 @@ function cmesh_new_periodic_tri2(; comm = mpi_comm())::t8_cmesh_t 0, 1.0, 0, 0, 1.0, 0, # tree 7, triangle 1.0, 0, 0, - 1.0, 1.0, 0, + 1.0, 1.0, 0 ] # Generally, one can define other geometries. But besides linear the other @@ -1802,7 +2061,7 @@ function cmesh_new_periodic_hybrid2(; comm = mpi_comm())::t8_cmesh_t 0, -2.0, 0, # tree 4, quad 2.0, 0, 0, -2.0, 0, 0, - 0, 2.0, 0, + 0, 2.0, 0 ] # This is how the cmesh looks like. The numbers are the tree numbers: diff --git a/test/test_mpi_t8code_fv.jl b/test/test_mpi_t8code_fv.jl index 80e1aa60984..5a2d9525558 100644 --- a/test/test_mpi_t8code_fv.jl +++ b/test/test_mpi_t8code_fv.jl @@ -30,7 +30,6 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") end @trixi_testset "second-order FV" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), - order=2, l2=[0.008142380494734171], linf=[0.018687916234976898]) # Ensure that we do not have excessive memory allocations @@ -42,8 +41,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end - # The extended reconstruction stencil currently is not mpi parallel. - # The current version runs through but an error occurs on some rank (somehow not printed to the terminal). + # 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, @@ -93,6 +92,110 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") 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, + initial_refinement_level=2, + cmesh=Trixi.cmesh_new_periodic_hybrid(), + l2=[0.2253867410593706], + linf=[0.34092690256865166]) + # 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 "first-order FV - triangles" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_basic_hybrid.jl"), + order=1, + initial_refinement_level=1, + cmesh=Trixi.cmesh_new_tri(trees_per_dimension = (2, 2)), + l2=[0.29924666807083133], + linf=[0.4581996753014146]) + # 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 "first-order FV - hybrid2" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_basic_hybrid.jl"), + order=1, + initial_refinement_level=2, + cmesh=Trixi.cmesh_new_periodic_hybrid2(), + l2=[0.20740154468889108], + linf=[0.4659917007721659]) + # 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"), + initial_refinement_level=2, + cmesh=Trixi.cmesh_new_periodic_hybrid(), + l2=[0.1296561675517274], + linf=[0.25952934874433753]) + # 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 + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_nonperiodic.jl"), + order=1, + l2=[0.07215018673798403], + linf=[0.12087525707243896]) + # 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_nonperiodic.jl"), + l2=[0.017076631443535124], + linf=[0.05613089948002803]) + # 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 + # TODO: Somehow, this non-periodic run with a triangular mesh is unstable. + # When fixed, also add here. + 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"), @@ -101,13 +204,13 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") 0.059376731961878, 0.019737707470047838, 0.019737707470047747, - 0.09982550390697936, + 0.09982550390697936 ], linf=[ 0.08501451493301548, 0.029105783468157398, 0.029105783468157842, - 0.1451756151490775, + 0.1451756151490775 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -125,13 +228,13 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") 0.031971635416962525, 0.016630983283552267, 0.016630873316960327, - 0.04813244762272231, + 0.04813244762272231 ], linf=[ 0.055105205670854085, 0.03647221946045942, 0.036470504033139894, - 0.0811201478913759, + 0.0811201478913759 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -152,13 +255,13 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") 0.5733341919395405, 0.11399976571202448, 0.1139997657120245, - 1.3548613737038315, + 1.3548613737038315 ], linf=[ 1.7328363346781415, 0.27645456051355827, 0.27645456051355827, - 2.6624886901791407, + 2.6624886901791407 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -179,13 +282,13 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") 0.25420413805862135, 0.22153054262689362, 0.11870842058617848, - 0.03626117501911353, + 0.03626117501911353 ], linf=[ 0.5467894727797227, 0.4156157752497065, 0.26176691230685767, - 0.0920609123083227, + 0.0920609123083227 ], tspan=(0.0, 1.0)) # Ensure that we do not have excessive memory allocations @@ -201,18 +304,18 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"), order=2, - cmesh=Trixi.cmesh_quad(periodicity = (true, true)), + cmesh=Trixi.cmesh_new_quad(periodicity = (true, true)), l2=[ 0.2307479238046326, 0.19300139957275295, 0.1176326315506721, - 0.020439850138732837, + 0.020439850138732837 ], linf=[ 0.5069212604421109, 0.365579474379667, 0.24226411409222004, - 0.049201093470609525, + 0.049201093470609525 ], tspan=(0.0, 1.0)) # Ensure that we do not have excessive memory allocations @@ -232,13 +335,13 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") 0.16181566308374545, 0.10090964843765918, 0.15229553744179888, - 0.0395037796064376, + 0.0395037796064376 ], linf=[ 0.6484515918189779, 0.3067327488921227, 0.34771083375083534, - 0.10713502930441887, + 0.10713502930441887 ], tspan=(0.0, 1.0)) # Ensure that we do not have excessive memory allocations diff --git a/test/test_t8code_fv_2d.jl b/test/test_t8code_fv_2d.jl index c4b1a0639ee..9ddef51a9d1 100644 --- a/test/test_t8code_fv_2d.jl +++ b/test/test_t8code_fv_2d.jl @@ -125,6 +125,7 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), order=1, initial_refinement_level=2, + cmesh=Trixi.cmesh_new_periodic_hybrid(), l2=[0.2253867410593706], linf=[0.34092690256865166]) # Ensure that we do not have excessive memory allocations @@ -136,9 +137,42 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + @trixi_testset "first-order FV - triangles" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + order=1, + initial_refinement_level=1, + cmesh=Trixi.cmesh_new_tri(trees_per_dimension = (2, 2)), + l2=[0.29924666807083133], + linf=[0.4581996753014146]) + # 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 "first-order FV - hybrid2" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_hybrid.jl"), + order=1, + initial_refinement_level=2, + cmesh=Trixi.cmesh_new_periodic_hybrid2(), + l2=[0.20740154468889108], + linf=[0.4659917007721659]) + # 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"), initial_refinement_level=2, + cmesh=Trixi.cmesh_new_periodic_hybrid(), l2=[0.1296561675517274], linf=[0.25952934874433753]) # Ensure that we do not have excessive memory allocations @@ -153,18 +187,49 @@ end end @trixi_testset "elixir_advection_nonperiodic.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"), - order=1, - l2=[0.07215018673798403], - linf=[0.12087525707243896]) - # 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 + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"), + order=1, + l2=[0.07215018673798403], + linf=[0.12087525707243896]) + # 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_nonperiodic.jl"), + l2=[0.017076631443535124], + linf=[0.05613089948002803]) + # 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 + # TODO: Somehow, this non-periodic run with a triangular mesh is unstable. + # When fixed, also add to mpi test file. + # @trixi_testset "second-order FV - triangles" begin + # @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"), + # cmesh=Trixi.cmesh_new_tri(periodicity = (false, false)), + # l2=[0.0], + # linf=[0.0]) + # # 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_euler_source_terms.jl" begin @@ -175,13 +240,13 @@ end 0.059376731961878, 0.019737707470047838, 0.019737707470047747, - 0.09982550390697936, + 0.09982550390697936 ], linf=[ 0.08501451493301548, 0.029105783468157398, 0.029105783468157842, - 0.1451756151490775, + 0.1451756151490775 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -199,13 +264,13 @@ end 0.031971635993647315, 0.016631028330554957, 0.016630833188111448, - 0.04813246238398825, + 0.04813246238398825 ], linf=[ 0.055105654108624336, 0.03647317645079773, 0.03647020577993976, - 0.08112180586875883, + 0.08112180586875883 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -224,13 +289,13 @@ end 0.058781550112786296, 0.02676026477370241, 0.02673090935979779, - 0.08033279155463603, + 0.08033279155463603 ], linf=[ 0.09591263836601072, 0.05351985245787505, 0.05264935415308125, - 0.14318629241962988, + 0.14318629241962988 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -252,13 +317,13 @@ end 0.5733341919395403, 0.11399976571202451, 0.11399976571202453, - 1.3548613737038324, + 1.3548613737038324 ], linf=[ 1.732836334678142, 0.27645456051355827, 0.27645456051355827, - 2.6624886901791407, + 2.6624886901791407 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -276,13 +341,13 @@ end 0.7331398527938754, 0.15194349346989244, 0.1519434934698924, - 1.299914264830515, + 1.299914264830515 ], linf=[ 2.2864127304524726, 0.3051023693829176, 0.30510236938291757, - 2.6171402581107936, + 2.6171402581107936 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -304,13 +369,13 @@ end 0.25420413805862135, 0.22153054262689362, 0.11870842058617848, - 0.03626117501911353, + 0.03626117501911353 ], linf=[ 0.5467894727797227, 0.4156157752497065, 0.26176691230685767, - 0.0920609123083227, + 0.0920609123083227 ], tspan=(0.0, 1.0)) # Ensure that we do not have excessive memory allocations @@ -326,18 +391,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"), order=2, - cmesh=Trixi.cmesh_quad(periodicity = (true, true)), + cmesh=Trixi.cmesh_new_quad(periodicity = (true, true)), l2=[ 0.2307479238046326, 0.19300139957275295, 0.1176326315506721, - 0.020439850138732837, + 0.020439850138732837 ], linf=[ 0.5069212604421109, 0.365579474379667, 0.24226411409222004, - 0.049201093470609525, + 0.049201093470609525 ], tspan=(0.0, 1.0)) # Ensure that we do not have excessive memory allocations @@ -357,13 +422,13 @@ end 0.16181566308374545, 0.10090964843765918, 0.15229553744179888, - 0.0395037796064376, + 0.0395037796064376 ], linf=[ 0.6484515918189779, 0.3067327488921227, 0.34771083375083534, - 0.10713502930441887, + 0.10713502930441887 ], tspan=(0.0, 1.0)) # 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 ec980288af1..320aefd4607 100644 --- a/test/test_t8code_fv_3d.jl +++ b/test/test_t8code_fv_3d.jl @@ -44,11 +44,12 @@ mkdir(outdir) # end # end +# NOTE: Since I use 1x1x1 tree instead of 8x8x8, I need to increase the resolution three times by the factor of 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, + initial_refinement_level=2+3, l2=[0.2848617953369851], linf=[0.3721898718954475]) # Ensure that we do not have excessive memory allocations @@ -62,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, + initial_refinement_level=2+3, l2=[0.10381089565603231], linf=[0.13787405651527007]) # Ensure that we do not have excessive memory allocations @@ -76,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=2, + initial_refinement_level=2+3, extended_reconstruction_stencil=true, l2=[0.24826100771542928], linf=[0.3799973662927152])