diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0c636ee8b0..39d9f532a4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -71,6 +71,7 @@ jobs: - p4est_part2 - t8code_part1 - t8code_part2 + - t8code_fv - unstructured_dgmulti - parabolic - paper_self_gravitating_gas_dynamics diff --git a/Project.toml b/Project.toml index 612b65ec6f..f7e203158b 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.6, 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 new file mode 100644 index 0000000000..49d32d9706 --- /dev/null +++ b/examples/t8code_2d_fv/elixir_advection_basic.jl @@ -0,0 +1,112 @@ +using OrdinaryDiffEq +using Trixi +using T8code + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +initial_condition = initial_condition_convergence_test + +solver = FV(order = 2, extended_reconstruction_stencil = false, + surface_flux = flux_lax_friedrichs) + +# 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. +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 + +function f_c() + @cfunction($f, Cvoid, + (t8_cmesh_t, t8_gloidx_t, Ptr{Cdouble}, Csize_t, + Ptr{Cdouble}, Ptr{Cvoid}, Ptr{Cvoid})) +end + +trees_per_dimension = (2, 2) + +eclass = T8_ECLASS_QUAD +mesh = T8codeMesh(trees_per_dimension, eclass, + # mapping = Trixi.trixi_t8_mapping_c(mapping), + mapping = f_c(), + # 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) + +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_2d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_2d_fv/elixir_advection_basic_hybrid.jl new file mode 100644 index 0000000000..41470fefdc --- /dev/null +++ b/examples/t8code_2d_fv/elixir_advection_basic_hybrid.jl @@ -0,0 +1,43 @@ +using OrdinaryDiffEq +using Trixi + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(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_new_quad(periodicity = (true, true)) +# cmesh = Trixi.cmesh_new_tri(periodicity = (true, true)) +mesh = T8codeMesh(cmesh, solver, + initial_refinement_level = 3) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +ode = semidiscretize(semi, (0.0, 1.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 = 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_2d_fv/elixir_advection_gauss.jl b/examples/t8code_2d_fv/elixir_advection_gauss.jl new file mode 100644 index 0000000000..177fce928c --- /dev/null +++ b/examples/t8code_2d_fv/elixir_advection_gauss.jl @@ -0,0 +1,45 @@ +using OrdinaryDiffEq +using Trixi + +#################################################### + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(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() +# 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) + +ode = semidiscretize(semi, (0.0, 20.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 = 10, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +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/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl b/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl new file mode 100644 index 0000000000..7e956fe175 --- /dev/null +++ b/examples/t8code_2d_fv/elixir_advection_nonperiodic.jl @@ -0,0 +1,65 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the linear advection equation. + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(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_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, + 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) + +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + save_solution, 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_2d_fv/elixir_euler_blast_wave.jl b/examples/t8code_2d_fv/elixir_euler_blast_wave.jl new file mode 100644 index 0000000000..93277d9783 --- /dev/null +++ b/examples/t8code_2d_fv/elixir_euler_blast_wave.jl @@ -0,0 +1,61 @@ +using OrdinaryDiffEq +using Trixi + +#################################################### + +equations = CompressibleEulerEquations2D(1.4) + +function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_blast_wave + +solver = FV(order = 2, extended_reconstruction_stencil = true, + surface_flux = flux_lax_friedrichs) + +initial_refinement_level = 4 +# cmesh = Trixi.cmesh_new_periodic_hybrid2() +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) + +ode = semidiscretize(semi, (0.0, 2.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 = 10, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.9) + +callbacks = CallbackSet(summary_callback, save_solution, analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + 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/examples/t8code_2d_fv/elixir_euler_free_stream.jl b/examples/t8code_2d_fv/elixir_euler_free_stream.jl new file mode 100644 index 0000000000..b6abb8aa5f --- /dev/null +++ b/examples/t8code_2d_fv/elixir_euler_free_stream.jl @@ -0,0 +1,73 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the compressible Euler equations. + +equations = CompressibleEulerEquations2D(1.4) + +initial_condition = initial_condition_constant + +solver = FV(order = 2, surface_flux = flux_lax_friedrichs) + +# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D +function mapping(xi_, eta_) + # Transform input variables between -1 and 1 onto [0,3] + xi = 1.5 * xi_ + 1.5 + eta = 1.5 * eta_ + 1.5 + + y = eta + 3 / 8 * (cos(1.5 * pi * (2 * xi - 3) / 3) * + cos(0.5 * pi * (2 * eta - 3) / 3)) + + x = xi + 3 / 8 * (cos(0.5 * pi * (2 * xi - 3) / 3) * + cos(2 * pi * (2 * y - 3) / 3)) + + return SVector(x, y) +end + +############################################################################### +# Get the uncurved mesh from a file (downloads the file if not available locally) + +# Unstructured mesh with 48 cells of the square domain [-1, 1]^n +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp", + joinpath(@__DIR__, "square_unstructured_1.inp")) + +# Note: +# Including mapping does not work yet for FV. +mesh = T8codeMesh(mesh_file, 2; polydeg = 0, + mapping = mapping, + initial_refinement_level = 2) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = Dict(:all => BoundaryConditionDirichlet(initial_condition))) + +############################################################################### +# 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.9) + +save_solution = SaveSolutionCallback(interval = 10, + solution_variables = cons2prim) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + 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() # print the timer summary diff --git a/examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl b/examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl new file mode 100644 index 0000000000..7adf83fc47 --- /dev/null +++ b/examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl @@ -0,0 +1,74 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +gamma = 1.4 +equations = CompressibleEulerEquations2D(gamma) + +""" + initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D) + +A version of the classical Kelvin-Helmholtz instability based on +- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021) + A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations + of the Euler Equations + [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017) +""" +function initial_condition_kelvin_helmholtz_instability(x, t, + equations::CompressibleEulerEquations2D) + # change discontinuity to tanh + # typical resolution 128^2, 256^2 + # domain size is [-1,+1]^2 + slope = 15 + amplitude = 0.02 + B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5) + rho = 0.5 + 0.75 * B + v1 = 0.5 * (B - 1) + v2 = 0.1 * sin(2 * pi * x[1]) + p = 1.0 + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_kelvin_helmholtz_instability + +solver = FV(order = 2, extended_reconstruction_stencil = false, + surface_flux = flux_lax_friedrichs) + +initial_refinement_level = 4 +# cmesh = Trixi.cmesh_new_periodic_hybrid() +# 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) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +tspan = (0.0, 3.7) +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 = 40, + 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) + +############################################################################### +# 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() # print the timer summary diff --git a/examples/t8code_2d_fv/elixir_euler_source_terms.jl b/examples/t8code_2d_fv/elixir_euler_source_terms.jl new file mode 100644 index 0000000000..74cc1aae7e --- /dev/null +++ b/examples/t8code_2d_fv/elixir_euler_source_terms.jl @@ -0,0 +1,52 @@ + +using OrdinaryDiffEq +using Trixi + +equations = CompressibleEulerEquations2D(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) + +cmesh = Trixi.cmesh_new_periodic_hybrid() +# cmesh = Trixi.cmesh_new_quad(periodicity = (true, true)) +# cmesh = Trixi.cmesh_new_tri(periodicity = (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 = 1.0) + +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/Trixi.jl b/src/Trixi.jl index d9b9245918..a1cf7722b7 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -245,6 +245,8 @@ export DG, export VolumeIntegralSubcellLimiting, BoundsCheckCallback, SubcellLimiterIDP, SubcellLimiterIDPCorrection +export FV + export nelements, nnodes, nvariables, eachelement, eachnode, eachvariable diff --git a/src/auxiliary/t8code.jl b/src/auxiliary/t8code.jl index 7c1399fc80..5200d483f7 100644 --- a/src/auxiliary/t8code.jl +++ b/src/auxiliary/t8code.jl @@ -190,3 +190,42 @@ function trixi_t8_adapt!(mesh, indicators) return differences end + +# Note and TODO: +# This routine seems to work with a "standard" mapping, but not with `coordinates2mapping`. +# Now, even when called directly within the elixir. +# - Is the function called with the correct parameters? Memory location correct? +# - Life time issue for the GC tracked Julia object used in C? +# function trixi_t8_mapping_c(mapping) +# function f(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data) +# @T8_ASSERT(cmesh isa Ptr{t8_cmesh}) +# @T8_ASSERT(gtreeid isa t8_gloidx_t) +# @T8_ASSERT(ref_coords isa Ptr{Cdouble}) +# @T8_ASSERT(num_coords isa Csize_t) +# @T8_ASSERT(out_coords isa Ptr{Cdouble}) +# @T8_ASSERT(tree_data isa Ptr{Cvoid}) +# @T8_ASSERT(user_data isa Ptr{Cvoid}) + +# 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/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 4cd92ce7c5..7251661d1c 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -180,6 +180,51 @@ function calc_error_norms(func, u, t, analyzer, return l2_error, linf_error end +function calc_error_norms(func, u, t, analyzer, + mesh::T8codeMesh, equations, + initial_condition, solver::FV, cache, cache_analysis) + # Set up data structures + l2_error = zero(func(get_node_vars(u, equations, solver, 1), equations)) + linf_error = copy(l2_error) + total_volume = zero(real(mesh)) + + # Iterate over all elements for error calculations + for element in eachelement(mesh, solver, cache) + midpoint = get_node_coords(cache.elements.midpoint, equations, solver, element) + volume = cache.elements.volume[element] + + u_exact = initial_condition(midpoint, t, equations) + diff = func(u_exact, equations) - + func(get_node_vars(u, equations, solver, element), equations) + l2_error += diff .^ 2 * volume + linf_error = @. max(linf_error, abs(diff)) + total_volume += volume + end + + # Accumulate local results on root process + if mpi_isparallel() + global_l2_error = Vector(l2_error) + global_linf_error = Vector(linf_error) + MPI.Reduce!(global_l2_error, +, mpi_root(), mpi_comm()) + MPI.Reduce!(global_linf_error, max, mpi_root(), mpi_comm()) + total_volume_ = MPI.Reduce(total_volume, +, mpi_root(), mpi_comm()) + if mpi_isroot() + l2_error = convert(typeof(l2_error), global_l2_error) + linf_error = convert(typeof(linf_error), global_linf_error) + # For L2 error, divide by total volume + l2_error = @. sqrt(l2_error / total_volume_) + else + l2_error = convert(typeof(l2_error), NaN * global_l2_error) + linf_error = convert(typeof(linf_error), NaN * global_linf_error) + end + else + # For L2 error, divide by total volume + l2_error = @. sqrt(l2_error / total_volume) + end + + return l2_error, linf_error +end + function integrate_via_indices(func::Func, u, mesh::TreeMesh{2}, equations, dg::DGSEM, cache, args...; normalize = true) where {Func} @@ -264,6 +309,54 @@ function integrate(func::Func, u, end end +function integrate_via_indices(func::Func, u, + mesh::T8codeMesh, equations, + solver::FV, cache, args...; + normalize = true) where {Func} + # Initialize integral with zeros of the right shape + integral = zero(func(u, 1, equations, solver, args...)) + total_volume = zero(real(mesh)) + + # Use quadrature to numerically integrate over entire domain + for element in eachelement(mesh, solver, cache) + volume = cache.elements.volume[element] + integral += volume * func(u, element, equations, solver, args...) + total_volume += volume + end + + if mpi_isparallel() + global_integral = MPI.Reduce!(Ref(integral), +, mpi_root(), mpi_comm()) + total_volume_ = MPI.Reduce(total_volume, +, mpi_root(), mpi_comm()) + if mpi_isroot() + integral = convert(typeof(integral), global_integral[]) + # Normalize with total volume + if normalize + integral = integral / total_volume_ + end + else + integral = convert(typeof(integral), NaN * integral) + total_volume_ = total_volume # non-root processes receive nothing from reduce -> overwrite + end + else + # Normalize with total volume + if normalize + integral = integral / total_volume + end + end + + return integral +end + +function integrate(func::Func, u, + mesh, + equations, solver::FV, cache; normalize = true) where {Func} + integrate_via_indices(u, mesh, equations, solver, cache; + normalize = normalize) do u, element, equations, solver + u_local = get_node_vars(u, equations, solver, element) + return func(u_local, equations) + end +end + function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, StructuredMeshView{2}, UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, @@ -277,6 +370,18 @@ function analyze(::typeof(entropy_timederivative), du, u, t, end end +function analyze(::typeof(entropy_timederivative), du, u, t, + mesh::T8codeMesh, + equations, solver::FV, cache) + # Calculate ∫(∂S/∂u ⋅ ∂u/∂t)dΩ + integrate_via_indices(u, mesh, equations, solver, cache, + du) do u, element, equations, solver, du + u_node = get_node_vars(u, equations, solver, element) + du_node = get_node_vars(du, equations, solver, element) + dot(cons2entropy(u_node, equations), du_node) + end +end + function analyze(::Val{:l2_divb}, du, u, t, mesh::TreeMesh{2}, equations::IdealGlmMhdEquations2D, dg::DGSEM, cache) diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl index 3f9abef7f3..58fedbb2d1 100644 --- a/src/callbacks_step/save_solution.jl +++ b/src/callbacks_step/save_solution.jl @@ -246,5 +246,32 @@ end # TODO: Taal refactor, move save_mesh_file? # function save_mesh_file(mesh::TreeMesh, output_directory, timestep=-1) in io/io.jl +# Temporary and experimental routine for SaveSolutionCallback. +# output_data_to_vtu directly creates .vtu files. +function save_solution_file(u, time, dt, timestep, + mesh::T8codeMesh, + equations, solver::FV, cache, + solution_callback, + element_variables = Dict{Symbol, Any}(), + node_variables = Dict{Symbol, Any}(); + system = "") + @unpack output_directory, solution_variables = solution_callback + + # Filename based on current time step + if isempty(system) + filename = joinpath(output_directory, @sprintf("solution_%09d.h5", timestep)) + else + filename = joinpath(output_directory, + @sprintf("solution_%s_%09d.h5", system, timestep)) + end + MPI.Barrier(MPI.COMM_WORLD) + Trixi.exchange_solution_data!(u, mesh, equations, solver, cache) + Trixi.output_data_to_vtu(mesh, equations, solver, + cache.communication_data.solution_data, filename, + solution_variables) + + return filename +end + include("save_solution_dg.jl") end # @muladd diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index c7922cecc6..7dc0a735fd 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -212,4 +212,40 @@ function max_dt(u, t, mesh::ParallelT8codeMesh{2}, return dt end + +function max_dt(u, t, mesh, + constant_speed::False, equations, solver::FV, cache) + dt = typemax(eltype(u)) + + for element in eachelement(mesh, solver, cache) + u_node = get_node_vars(u, equations, solver, element) + lambda1, lambda2 = max_abs_speeds(u_node, equations) + dx = cache.elements.dx[element] + dt = min(dt, dx / (lambda1 + lambda2)) + end + + if mpi_isparallel() + dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[] + end + + return dt +end + +function max_dt(u, t, mesh, + constant_speed::True, equations, solver::FV, cache) + dx_min = typemax(eltype(u)) + + max_lambda1, max_lambda2 = max_abs_speeds(equations) + for element in eachelement(mesh, solver, cache) + dx = cache.elements.dx[element] + dx_min = min(dx_min, dx) + end + dt = dx_min / (max_lambda1 + max_lambda2) + + if mpi_isparallel() + dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[] + end + + return dt +end end # @muladd diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index 127bc420f6..e060159dbc 100644 --- a/src/meshes/mesh_io.jl +++ b/src/meshes/mesh_io.jl @@ -227,7 +227,8 @@ end # TODO: Implement this function as soon as there is support for this in `t8code`. function save_mesh_file(mesh::T8codeMesh, output_directory, timestep, mpi_parallel) - error("Mesh file output not supported yet for `T8codeMesh`.") + # TODO: commented following line. SaveSolutionCallback doesn't work for FV with this error. + # error("Mesh file output not supported yet for `T8codeMesh`.") return joinpath(output_directory, "dummy_mesh.h5") end diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 231e296566..51078f8eec 100644 --- a/src/meshes/t8code_mesh.jl +++ b/src/meshes/t8code_mesh.jl @@ -18,10 +18,12 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: boundary_names :: Array{Symbol, 2} # [face direction, tree] current_filename :: String + unsaved_changes :: Bool - ninterfaces :: Int - nmortars :: Int - nboundaries :: Int + max_number_faces :: Int + ninterfaces :: Int + nmortars :: Int + nboundaries :: Int nmpiinterfaces :: Int nmpimortars :: Int @@ -31,13 +33,17 @@ mutable struct T8codeMesh{NDIMS, RealT <: Real, IsParallel, NDIMSP2, NNODES} <: current_filename) where {NDIMS} is_parallel = mpi_isparallel() ? True() : False() + @assert t8_forest_get_local_num_elements(forest)>0 "Too many ranks to properly partition the mesh!" + mesh = new{NDIMS, Float64, typeof(is_parallel), NDIMS + 2, length(nodes)}(forest, is_parallel) mesh.nodes = nodes mesh.boundary_names = boundary_names mesh.current_filename = current_filename + mesh.unsaved_changes = false # TODO: if set to `true`, the mesh will be saved to a mesh file. mesh.tree_node_coordinates = tree_node_coordinates + mesh.max_number_faces = 2 * NDIMS # TODO: How to automatically adapt for other element types, e.g. triangles? finalizer(mesh) do mesh # When finalizing `mesh.forest`, `mesh.scheme` and `mesh.cmesh` are @@ -342,6 +348,161 @@ 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(trees_per_dimension[1]==trees_per_dimension[2], + "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 + error("3D not supported") + 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, @@ -1268,9 +1429,553 @@ function fill_mesh_info!(mesh::T8codeMesh, interfaces, mortars, boundaries, return nothing end +function fill_mesh_info_fv!(mesh::T8codeMesh, interfaces, boundaries, + boundary_names; mpi_mesh_info = nothing) + @assert t8_forest_is_committed(mesh.forest) != 0 + + num_local_elements = t8_forest_get_local_num_elements(mesh.forest) + num_local_trees = t8_forest_get_num_local_trees(mesh.forest) + + # Process-local index of the current element in the space-filling curve. + current_index = t8_locidx_t(0) + + # Increment counters for the different interface/mortar/boundary types. + local_num_conform = 0 + local_num_mortars = 0 + local_num_boundary = 0 + + local_num_mpi_conform = 0 + local_num_mpi_mortars = 0 + + # Helper variables to compute unique global MPI interface/mortar ids. + max_level = t8_forest_get_maxlevel(mesh.forest) #UInt64 + max_tree_num_elements = UInt64(2^ndims(mesh))^max_level + + # Loop over all local trees. + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(mesh.forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(mesh.forest, tree_class) + + num_elements_in_tree = t8_forest_get_tree_num_elements(mesh.forest, itree) + + global_itree = t8_forest_global_tree_id(mesh.forest, itree) + + # Loop over all local elements of the current local tree. + for ielement in 0:(num_elements_in_tree - 1) + element = t8_forest_get_element_in_tree(mesh.forest, itree, ielement) + + level = t8_element_level(eclass_scheme, element) + + num_faces = t8_element_num_faces(eclass_scheme, element) + + # Loop over all faces of the current local element. + for iface in 0:(num_faces - 1) + pelement_indices_ref = Ref{Ptr{t8_locidx_t}}() + pneighbor_leaves_ref = Ref{Ptr{Ptr{t8_element}}}() + pneigh_scheme_ref = Ref{Ptr{t8_eclass_scheme}}() + + dual_faces_ref = Ref{Ptr{Cint}}() + num_neighbors_ref = Ref{Cint}() + + forest_is_balanced = Cint(1) + + # Query neighbor information from t8code. + t8_forest_leaf_face_neighbors(mesh.forest, itree, element, + pneighbor_leaves_ref, iface, dual_faces_ref, + num_neighbors_ref, + pelement_indices_ref, pneigh_scheme_ref, + forest_is_balanced) + + num_neighbors = num_neighbors_ref[] + dual_faces = unsafe_wrap(Array, dual_faces_ref[], num_neighbors) + neighbor_ielements = unsafe_wrap(Array, pelement_indices_ref[], + num_neighbors) + neighbor_leaves = unsafe_wrap(Array, pneighbor_leaves_ref[], num_neighbors) + neighbor_scheme = pneigh_scheme_ref[] + + # Now we check for the different cases. The nested if-structure is as follows: + # + # if `boundary`: + # + # + # else: // It must be an interface or mortar. + # + # if `local interface`: + # + # elseif `local mortar from larger element point of view`: + # Not supported yet + # else: // `local mortar from smaller elements point of view` + # // We only count local mortars once. + # + # // end + + # Domain boundary. + if num_neighbors == 0 + local_num_boundary += 1 + boundary_id = local_num_boundary + + # One-based indexing. + boundaries.neighbor_ids[boundary_id] = current_index + 1 + boundaries.faces[boundary_id] = iface + 1 + boundaries.name[boundary_id] = boundary_names[iface + 1, itree + 1] + else # Interface or mortar. + neighbor_level = t8_element_level(neighbor_scheme, neighbor_leaves[1]) + + # Local interface: The second condition ensures we only visit the interface once. + if level == neighbor_level && current_index <= neighbor_ielements[1] + local_num_conform += 1 + interfaces.neighbor_ids[1, local_num_conform] = current_index + + 1 + interfaces.neighbor_ids[2, local_num_conform] = neighbor_ielements[1] + + 1 + + interfaces.faces[1, local_num_conform] = iface + 1 + interfaces.faces[2, local_num_conform] = dual_faces[1] + 1 + + # Local mortar. + elseif level < neighbor_level + error("Mortars are not supported yet!") + end + end + + t8_free(dual_faces_ref[]) + t8_free(pneighbor_leaves_ref[]) + t8_free(pelement_indices_ref[]) + end # for iface + + current_index += 1 + end # for ielement + end # for itree + + return nothing +end + #! format: off @deprecate T8codeMesh{2}(conn::Ptr{p4est_connectivity}; kwargs...) T8codeMesh(conn::Ptr{p4est_connectivity}; kwargs...) @deprecate T8codeMesh{3}(conn::Ptr{p8est_connectivity}; kwargs...) T8codeMesh(conn::Ptr{p8est_connectivity}; kwargs...) @deprecate T8codeMesh{2}(meshfile::String; kwargs...) T8codeMesh(meshfile::String, 2; kwargs...) @deprecate T8codeMesh{3}(meshfile::String; kwargs...) T8codeMesh(meshfile::String, 3; kwargs...) #! format: on + +# Write the forest as vtu and also write the element's volumes in the file. +# +# t8code supports writing element based data to vtu as long as its stored +# as doubles. Each of the data fields to write has to be provided in its own +# array of length num_local_elements. +# t8code supports two types: T8_VTK_SCALAR - One double per element. +# and T8_VTK_VECTOR - Three doubles per element. +function output_data_to_vtu(mesh::T8codeMesh, equations, solver, + u_tmp, out, solution_variables) + vars = varnames(solution_variables, equations) + + vtk_data = Vector{t8_vtk_data_field_t}(undef, nvariables(equations)) + + data = Array{Float64}(undef, ncells(mesh), nvariables(equations)) + for element in 1:ncells(mesh) + variables = solution_variables(u_tmp[element].u, equations) + for v in eachvariable(equations) + data[element, v] = variables[v] + end + end + + GC.@preserve data begin + for v in eachvariable(equations) + data_ptr = pointer(@views(data[:, v])) + vtk_data[v] = t8_vtk_data_field_t(T8_VTK_SCALAR, + NTuple{8192, Cchar}(rpad("$(vars[v])\0", + 8192, ' ')), + data_ptr) + end + + # The number of user defined data fields to write. + num_data = length(vtk_data) + + # Write user defined data to vtu file. + write_treeid = 1 + write_mpirank = 1 + write_level = 1 + write_element_id = 1 + write_ghosts = 0 + t8_forest_write_vtk_ext(mesh.forest, out, write_treeid, write_mpirank, + write_level, write_element_id, write_ghosts, + 0, 0, num_data, pointer(vtk_data)) + end +end + +# Simple meshes +# Temporary routines to create simple `cmesh`s by hand + +# Directly ported from: `src/t8_cmesh/t8_cmesh_examples.c: t8_cmesh_new_periodic_hybrid`. +function cmesh_new_periodic_hybrid(; 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 + 0, -1.0, 0, + 0, 0, 0, + -1.0, -1.0, 0, # tree 1, triangle + 0, 0, 0, + -1.0, 0, 0, + 0, -1.0, 0, # tree 2, quad + 1.0, -1.0, 0, + 0, 0, 0, + 1.0, 0, 0, + -1.0, 0, 0, # tree 3, quad + 0, 0, 0, + -1.0, 1.0, 0, + 0, 1.0, 0, + 0, 0, 0, # tree 4, triangle + 1.0, 0, 0, + 1.0, 1.0, 0, + 0, 0, 0, # tree 5, triangle + 1.0, 1.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) + + # This is how the cmesh looks like. The numbers are the tree numbers: + # Domain size [-1,1]^2 + # + # +---+---+ + # | |5 /| + # | 3 | / | + # | |/ 4| + # +---+---+ + # |1 /| | + # | / | 2 | + # |/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_class(cmesh, 2, T8_ECLASS_QUAD) + t8_cmesh_set_tree_class(cmesh, 3, T8_ECLASS_QUAD) + t8_cmesh_set_tree_class(cmesh, 4, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 5, 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_tree_vertices(cmesh, 2, @views(vertices[(1 + 18):end]), 4) + t8_cmesh_set_tree_vertices(cmesh, 3, @views(vertices[(1 + 30):end]), 4) + t8_cmesh_set_tree_vertices(cmesh, 4, @views(vertices[(1 + 42):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 5, @views(vertices[(1 + 51):end]), 3) + + t8_cmesh_set_join(cmesh, 0, 1, 1, 2, 0) + t8_cmesh_set_join(cmesh, 0, 2, 0, 0, 0) + t8_cmesh_set_join(cmesh, 0, 3, 2, 3, 0) + + t8_cmesh_set_join(cmesh, 1, 3, 0, 2, 1) + t8_cmesh_set_join(cmesh, 1, 2, 1, 1, 0) + + t8_cmesh_set_join(cmesh, 2, 4, 3, 2, 0) + t8_cmesh_set_join(cmesh, 2, 5, 2, 0, 1) + + t8_cmesh_set_join(cmesh, 3, 5, 1, 1, 0) + t8_cmesh_set_join(cmesh, 3, 4, 0, 0, 0) + + t8_cmesh_set_join(cmesh, 4, 5, 1, 2, 0) + + t8_cmesh_commit(cmesh, comm) + + return cmesh +end + +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: + # + # +---+ + # | | + # | 0 | + # | | + # +---+ + # + + 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 + +# 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: + # + # +---+ + # |1 /| + # | / | + # |/0 | + # +---+ + # + + # 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 + -1.0, -1.0, 0, # tree 0, triangle + 0, -1.0, 0, + 0, 0, 0, + -1.0, -1.0, 0, # tree 1, triangle + 0, 0, 0, + -1.0, 0, 0, + 0, -1.0, 0, # tree 2, triangle + 1.0, -1.0, 0, + 1.0, 0, 0, + 0, -1.0, 0, # tree 3, triangle + 1.0, 0, 0, + 0, 0, 0, + -1.0, 0, 0, # tree 4, triangle + 0, 0, 0, + -1.0, 1.0, 0, + -1.0, 1.0, 0, # tree 5, triangle + 0, 0, 0, + 0, 1.0, 0, + 0, 0, 0, # tree 6, triangle + 1.0, 0, 0, + 0, 1.0, 0, + 0, 1.0, 0, # tree 7, triangle + 1.0, 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: + # + # +---+---+ + # |\ 5|\ 7| + # | \ | \ | + # |4 \| 6\| + # +---+---+ + # |1 /|3 /| + # | / | / | + # |/0 |/ 2| + # +---+---+ + # + + 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_class(cmesh, 2, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 3, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 4, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 5, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 6, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 7, 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_tree_vertices(cmesh, 2, @views(vertices[(1 + 18):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 3, @views(vertices[(1 + 27):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 4, @views(vertices[(1 + 36):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 5, @views(vertices[(1 + 45):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 6, @views(vertices[(1 + 54):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 7, @views(vertices[(1 + 63):end]), 3) + + t8_cmesh_set_join(cmesh, 0, 1, 1, 2, 0) + t8_cmesh_set_join(cmesh, 0, 3, 0, 1, 0) + t8_cmesh_set_join(cmesh, 0, 5, 2, 1, 0) + + t8_cmesh_set_join(cmesh, 1, 4, 0, 2, 1) + t8_cmesh_set_join(cmesh, 1, 2, 1, 0, 0) + + t8_cmesh_set_join(cmesh, 2, 3, 1, 2, 0) + t8_cmesh_set_join(cmesh, 2, 7, 2, 1, 0) + + t8_cmesh_set_join(cmesh, 3, 6, 0, 2, 1) + + t8_cmesh_set_join(cmesh, 4, 5, 0, 2, 1) + t8_cmesh_set_join(cmesh, 4, 7, 1, 0, 0) + + t8_cmesh_set_join(cmesh, 5, 6, 0, 1, 0) + + t8_cmesh_set_join(cmesh, 6, 7, 0, 2, 1) + + t8_cmesh_commit(cmesh, comm) + + return cmesh +end + +function cmesh_new_periodic_hybrid2(; comm = mpi_comm())::t8_cmesh_t + n_dims = 2 + vertices = [ # Just all vertices of all trees. partly duplicated + -2.0, -2.0, 0, # tree 0, triangle + 0, -2.0, 0, + -2.0, 0, 0, + -2.0, 2.0, 0, # tree 1, triangle + -2.0, 0, 0, + 0, 2.0, 0, + 2.0, -2.0, 0, # tree 2, triangle + 2.0, 0, 0, + 0, -2.0, 0, + 2.0, 2.0, 0, # tree 3, triangle + 0, 2.0, 0, + 2.0, 0, 0, + 0, -2.0, 0, # tree 4, quad + 2.0, 0, 0, + -2.0, 0, 0, + 0, 2.0, 0 + ] + + # This is how the cmesh looks like. The numbers are the tree numbers: + # Domain size [-2,2]^2 + # + # +----------+ + # | 1 /\ 3 | + # | / \ | + # | / \ | + # | / \ | + # |/ 4 \| + # |\ /| + # | \ / | + # | \ / | + # | 0 \ / 2 | + # | \/ | + # +----------+ + # + + # 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) + + 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_class(cmesh, 2, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 3, T8_ECLASS_TRIANGLE) + t8_cmesh_set_tree_class(cmesh, 4, T8_ECLASS_QUAD) + + 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_tree_vertices(cmesh, 2, @views(vertices[(1 + 18):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 3, @views(vertices[(1 + 27):end]), 3) + t8_cmesh_set_tree_vertices(cmesh, 4, @views(vertices[(1 + 36):end]), 4) + + t8_cmesh_set_join(cmesh, 0, 4, 0, 0, 0) + t8_cmesh_set_join(cmesh, 0, 2, 1, 2, 0) + t8_cmesh_set_join(cmesh, 0, 1, 2, 1, 0) + + t8_cmesh_set_join(cmesh, 1, 4, 0, 3, 0) + t8_cmesh_set_join(cmesh, 1, 3, 2, 1, 0) + + t8_cmesh_set_join(cmesh, 2, 4, 0, 2, 1) + t8_cmesh_set_join(cmesh, 2, 3, 1, 2, 0) + + t8_cmesh_set_join(cmesh, 3, 4, 0, 1, 1) + + t8_cmesh_commit(cmesh, comm) + + return cmesh +end diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl new file mode 100644 index 0000000000..18007d493d --- /dev/null +++ b/src/solvers/fv_t8code/containers.jl @@ -0,0 +1,614 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# Container data structure (structure-of-arrays style) for DG elements +mutable struct T8codeFVElementContainer{NDIMS, RealT <: Real, uEltype <: Real} + level::Vector{Cint} # [element] + volume::Vector{RealT} # [element] + midpoint::Matrix{RealT} # [dimension, element] + dx::Vector{RealT} # [element] - characteristic length (for CFL condition). + + num_faces::Vector{Cint} # [element] + face_midpoints::Array{RealT, 3} # [dimension, face, element] + face_areas::Matrix{RealT} # [face, element] + face_normals::Array{RealT, 3} # [dimension, face, element] + + reconstruction_stencil::Vector{Vector{Int}} # Reconstruction stencil vector with neighbors per [element] + reconstruction_distance::Vector{Vector{Vector{RealT}}} # Reconstruction stencil vector with distances per [element] + reconstruction_corner_elements::Array{Vector{Int}, 2} # Reconstruction stencil array with neighbor elements at corner per [corner, element] + reconstruction_gradient::Array{RealT, 4} # [dimension, variable, slope_stencils, element], slope_stencil: i = use all neighbors exclusive i, n_neighbors = use all neighbors + reconstruction_gradient_limited::Array{RealT, 3} # [dimension, variable, element] + + # internal `resize!`able storage + _midpoint::Vector{RealT} + _face_midpoints::Vector{RealT} + _face_areas::Vector{RealT} + _face_normals::Vector{RealT} + _reconstruction_corner_elements::Vector{Vector{Int}} + _reconstruction_gradient::Vector{RealT} + _reconstruction_gradient_limited::Vector{RealT} +end + +@inline Base.ndims(::T8codeFVElementContainer{NDIMS}) where {NDIMS} = NDIMS +@inline function Base.eltype(::T8codeFVElementContainer{NDIMS, RealT, uEltype}) where { + NDIMS, + RealT, + uEltype + } + uEltype +end + +@inline is_ghost_cell(element, mesh) = element > ncells(mesh) + +# See explanation of Base.resize! for the element container +function Base.resize!(elements::T8codeFVElementContainer, capacity) + (; _midpoint, _face_midpoints, _face_areas, _face_normals, _reconstruction_gradient, _reconstruction_gradient_limited, _reconstruction_corner_elements) = elements + + n_dims = ndims(elements) + n_variables = size(elements.reconstruction_gradient, 2) + max_number_faces = size(elements.face_midpoints, 2) + + resize!(elements.level, capacity) + resize!(elements.volume, capacity) + resize!(elements.dx, capacity) + resize!(elements.num_faces, capacity) + + resize!(_midpoint, n_dims * capacity) + elements.midpoint = unsafe_wrap(Array, pointer(_midpoint), (ndims, capacity)) + + resize!(_face_midpoints, n_dims * max_number_faces * capacity) + elements.face_midpoints = unsafe_wrap(Array, pointer(_face_midpoints), + (ndims, max_number_faces, capacity)) + + resize!(_face_areas, max_number_faces * capacity) + elements.face_areas = unsafe_wrap(Array, pointer(_face_areas), + (max_number_faces, capacity)) + + resize!(_face_normals, n_dims * max_number_faces * capacity) + elements.face_normals = unsafe_wrap(Array, pointer(_face_normals), + (ndims, max_number_faces, capacity)) + + resize!(elements.reconstruction_stencil, capacity) + resize!(elements.reconstruction_distance, capacity) + + resize!(_reconstruction_corner_elements, max_number_faces * capacity) + elements.reconstruction_corner_elements = unsafe_wrap(Array, + pointer(_reconstruction_corner_elements), + (max_number_faces, capacity)) + + resize!(_reconstruction_gradient, + n_dims * n_variables * (max_number_faces + 1) * capacity) + elements.reconstruction_gradient = unsafe_wrap(Array, + pointer(_reconstruction_gradient), + (ndims, n_variables, + (max_number_faces + 1), capacity)) + + resize!(_reconstruction_gradient_limited, n_dims * n_variables * capacity) + elements.reconstruction_gradient_limited = unsafe_wrap(Array, + pointer(_reconstruction_gradient_limited), + (ndims, n_variables, + capacity)) + + return nothing +end + +# Create element container and initialize element data +function init_elements(mesh::T8codeMesh{NDIMS, RealT}, + equations, + solver::FV, + ::Type{uEltype}) where {NDIMS, RealT <: Real, uEltype <: Real} + (; forest) = mesh + # Check that the forest is a committed. + @assert(t8_forest_is_committed(forest)==1) + + nelements = ncells(mesh) + n_variables = nvariables(equations) + (; max_number_faces) = mesh + + level = Vector{Cint}(undef, nelements) + volume = Vector{RealT}(undef, nelements) + dx = Vector{RealT}(undef, nelements) + num_faces = Vector{Cint}(undef, nelements) + + _midpoint = Vector{RealT}(undef, NDIMS * nelements) + midpoint = unsafe_wrap(Array, pointer(_midpoint), (NDIMS, nelements)) + + _face_midpoints = Vector{RealT}(undef, NDIMS * max_number_faces * nelements) + face_midpoints = unsafe_wrap(Array, pointer(_face_midpoints), + (NDIMS, max_number_faces, nelements)) + + _face_areas = Vector{RealT}(undef, max_number_faces * nelements) + face_areas = unsafe_wrap(Array, pointer(_face_areas), (max_number_faces, nelements)) + + _face_normals = Vector{RealT}(undef, NDIMS * max_number_faces * nelements) + face_normals = unsafe_wrap(Array, pointer(_face_normals), + (NDIMS, max_number_faces, nelements)) + + reconstruction_stencil = Vector{Vector{Int}}(undef, nelements) + reconstruction_distance = Vector{Vector{Array{RealT, NDIMS}}}(undef, nelements) + + _reconstruction_corner_elements = Vector{Vector{Int}}(undef, + max_number_faces * nelements) + reconstruction_corner_elements = unsafe_wrap(Array, + pointer(_reconstruction_corner_elements), + (max_number_faces, nelements)) + + _reconstruction_gradient = Vector{RealT}(undef, + NDIMS * n_variables * + (max_number_faces + 1) * nelements) + reconstruction_gradient = unsafe_wrap(Array, pointer(_reconstruction_gradient), + (NDIMS, n_variables, max_number_faces + 1, + nelements)) + + _reconstruction_gradient_limited = Vector{RealT}(undef, + NDIMS * n_variables * nelements) + reconstruction_gradient_limited = unsafe_wrap(Array, + pointer(_reconstruction_gradient_limited), + (NDIMS, n_variables, nelements)) + + elements = T8codeFVElementContainer{NDIMS, RealT, uEltype}(level, volume, midpoint, + dx, num_faces, + face_midpoints, + face_areas, face_normals, + reconstruction_stencil, + reconstruction_distance, + reconstruction_corner_elements, + reconstruction_gradient, + reconstruction_gradient_limited, + _midpoint, + _face_midpoints, + _face_areas, + _face_normals, + _reconstruction_corner_elements, + _reconstruction_gradient, + _reconstruction_gradient_limited) + + init_elements!(elements, mesh, solver) + + return elements +end + +function init_elements!(elements, mesh::T8codeMesh, solver::FV) + (; forest) = mesh + + n_dims = ndims(mesh) + (; level, volume, dx, num_faces, face_areas, face_midpoints, face_normals) = elements + + midpoint = Vector{Cdouble}(undef, n_dims) + + face_midpoint = Vector{Cdouble}(undef, 3) # Need NDIMS=3 for t8code API + face_normal = Vector{Cdouble}(undef, 3) # Need NDIMS=3 for t8code API + corners = Array{Cdouble, 3}(undef, n_dims, mesh.max_number_faces, length(volume)) + corner_node = Vector{Cdouble}(undef, n_dims) + + num_local_trees = t8_forest_get_num_local_trees(forest) + + # Loop over all local trees in the forest. + current_index = 0 + for itree in 0:(num_local_trees - 1) + tree_class = t8_forest_get_tree_class(forest, itree) + eclass_scheme = t8_forest_get_eclass_scheme(forest, tree_class) + + # Get the number of elements of this tree. + num_elements_in_tree = t8_forest_get_tree_num_elements(forest, itree) + + # Loop over all local elements in the tree. + for ielement in 0:(num_elements_in_tree - 1) + current_index += 1 # Note: Julia has 1-based indexing, while C/C++ starts with 0. + + element = t8_forest_get_element_in_tree(forest, itree, ielement) + + level[current_index] = t8_element_level(eclass_scheme, element) + volume[current_index] = t8_forest_element_volume(forest, itree, element) + + t8_forest_element_centroid(forest, itree, element, pointer(midpoint)) + for dim in 1:n_dims + elements.midpoint[dim, current_index] = midpoint[dim] + end + + # Characteristic length of the element. It is an approximation since only + # for the element type `lines` this can be exact. + dx[current_index] = t8_forest_element_diam(forest, itree, element) + + # Loop over all faces of an element. + num_faces[current_index] = t8_element_num_faces(eclass_scheme, element) + + for iface in 1:num_faces[current_index] + # C++ is zero-indexed + t8_forest_element_coordinate(forest, itree, element, iface - 1, + @views(corner_node)) + + t8_forest_element_face_centroid(forest, itree, element, iface - 1, + face_midpoint) + + face_area = t8_forest_element_face_area(forest, itree, element, + iface - 1) + face_areas[iface, current_index] = face_area + + t8_forest_element_face_normal(forest, itree, element, iface - 1, + face_normal) + for dim in 1:n_dims + corners[dim, iface, current_index] = corner_node[dim] + face_midpoints[dim, iface, current_index] = face_midpoint[dim] + face_normals[dim, iface, current_index] = face_normal[dim] + end + end + end + end + + # Init stencil for reconstruction + if solver.extended_reconstruction_stencil + init_extended_reconstruction_stencil!(corners, elements, solver) + end + + return nothing +end + +@inline function init_extended_reconstruction_stencil!(corners, elements, solver) + if solver.order != 2 + return nothing + end + (; reconstruction_stencil, reconstruction_distance, reconstruction_corner_elements) = elements + (; volume, num_faces) = elements + + # Create empty vectors for every element + for element in eachindex(volume) + reconstruction_stencil[element] = [] + reconstruction_distance[element] = [] + for corner in 1:num_faces[element] + reconstruction_corner_elements[corner, element] = [] + end + end + + # Add all stencil neighbors to list; including doubled elements + for element in eachindex(volume) + # loop over all elements with higher index + for possible_stencil_neighbor in (element + 1):length(volume) + # loop over all corners of `element` + for corner in 1:num_faces[element] + corner_coords = view(corners, :, corner, element) + # loop over all corners of `possible_stencil_neighbor` + for possible_corner in 1:num_faces[possible_stencil_neighbor] + possible_corner_coords = view(corners, :, possible_corner, + possible_stencil_neighbor) + if corner_coords == possible_corner_coords + neighbor = possible_stencil_neighbor + + midpoint_element = view(elements.midpoint, :, element) + midpoint_neighbor = view(elements.midpoint, :, neighbor) + + distance = midpoint_neighbor .- midpoint_element + append!(reconstruction_stencil[element], neighbor) + push!(reconstruction_distance[element], distance) + append!(reconstruction_stencil[neighbor], element) + push!(reconstruction_distance[neighbor], -distance) + append!(reconstruction_corner_elements[corner, element], + neighbor) + append!(reconstruction_corner_elements[possible_corner, + neighbor], element) + + # elseif # TODO: Handle periodic boundaries; Something like: + # distance = (face_midpoint_element .- midpoint_element) .+ + # (midpoint_neighbor .- face_midpoint_neighbor) + end + end + end + end + end + + # Remove all doubled elements from vectors + for element in eachindex(volume) + for i in length(reconstruction_stencil[element]):-1:1 + neighbor = reconstruction_stencil[element][i] + if neighbor in reconstruction_stencil[element][1:(i - 1)] + popat!(reconstruction_stencil[element], i) + popat!(reconstruction_distance[element], i) + end + end + end + + for element in eachindex(volume) + for corner in 1:num_faces[element] + for i in length(reconstruction_corner_elements[corner, element]):-1:1 + neighbor = reconstruction_corner_elements[corner, element][i] + if neighbor in reconstruction_corner_elements[corner, element][1:(i - 1)] + popat!(reconstruction_corner_elements[corner, element], i) + end + end + end + end + + return nothing +end + +function init_reconstruction_stencil!(elements, interfaces, boundaries, + communication_data, + mesh, equations, solver) + if solver.order != 2 + return nothing + end + (; reconstruction_stencil, reconstruction_distance) = elements + (; neighbor_ids, faces) = interfaces + (; domain_data) = communication_data + + # Create empty vectors for every element + for element in eachindex(reconstruction_stencil) + reconstruction_stencil[element] = [] + reconstruction_distance[element] = [] + end + + for interface in axes(neighbor_ids, 2) + element1 = neighbor_ids[1, interface] + element2 = neighbor_ids[2, interface] + face_element1 = faces[1, interface] + face_element2 = faces[2, interface] + + midpoint_element1 = domain_data[element1].midpoint + midpoint_element2 = domain_data[element2].midpoint + 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? + if isapprox(face_midpoint_element1, face_midpoint_element2) + distance = midpoint_element2 .- midpoint_element1 + else + distance = (face_midpoint_element1 .- midpoint_element1) .+ + (midpoint_element2 .- face_midpoint_element2) + end + append!(reconstruction_stencil[element1], element2) + push!(reconstruction_distance[element1], distance) + # only if element2 is local element + if element2 <= ncells(mesh) + append!(reconstruction_stencil[element2], element1) + push!(reconstruction_distance[element2], -distance) + end + end + + 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] + neighbor_ids::Matrix{Int} # [primary/secondary, interface] + faces::Matrix{Int} # [primary/secondary, interface] + + # internal `resize!`able storage + _u::Vector{uEltype} + _neighbor_ids::Vector{Int} + _faces::Vector{Int} +end + +@inline function ninterfaces(interfaces::T8codeFVInterfaceContainer) + size(interfaces.neighbor_ids, 2) +end + +# See explanation of Base.resize! for the element container +function Base.resize!(interfaces::T8codeFVInterfaceContainer, capacity) + (; _u, _neighbor_ids, _faces) = interfaces + + n_variables = size(interfaces.u, 2) + + resize!(_u, 2 * n_variables * capacity) + interfaces.u = unsafe_wrap(Array, pointer(_u), + (2, n_variables, capacity)) + + resize!(_neighbor_ids, 2 * capacity) + interfaces.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, capacity)) + + resize!(_faces, 2 * capacity) + interfaces.faces = unsafe_wrap(Array, pointer(_faces), (2, capacity)) + + return nothing +end + +# Create interface container and initialize interface data. +function init_interfaces(mesh::T8codeMesh, equations, solver::FV, uEltype) + # Initialize container + n_interfaces = count_required_surfaces(mesh).interfaces + if mpi_parallel(mesh) == true + n_interfaces += count_required_surfaces(mesh).mpi_interfaces + end + + _u = Vector{uEltype}(undef, + 2 * nvariables(equations) * n_interfaces) + u = unsafe_wrap(Array, pointer(_u), + (2, nvariables(equations), n_interfaces)) + + _neighbor_ids = Vector{Int}(undef, 2 * n_interfaces) + neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, n_interfaces)) + + _faces = Vector{Int}(undef, 2 * n_interfaces) + faces = unsafe_wrap(Array, pointer(_faces), (2, n_interfaces)) + + interfaces = T8codeFVInterfaceContainer{uEltype}(u, neighbor_ids, faces, + _u, _neighbor_ids, _faces) + + return interfaces +end + +mutable struct T8codeFVBoundaryContainer{uEltype <: Real} <: AbstractContainer + u::Array{uEltype, 2} # [variable, boundary] + neighbor_ids::Vector{Int} # [boundary] + faces::Vector{Int} # [boundary] + name::Vector{Symbol} # [boundary] + + # internal `resize!`able storage + _u::Vector{uEltype} +end + +@inline function nboundaries(boundaries::T8codeFVBoundaryContainer) + length(boundaries.neighbor_ids) +end + +# See explanation of Base.resize! for the element container +function Base.resize!(boundaries::T8codeFVBoundaryContainer, capacity) + (; _u, neighbor_ids, faces, name) = boundaries + + n_variables = size(boundaries.u, 1) + + resize!(_u, n_variables * capacity) + boundaries.u = unsafe_wrap(Array, pointer(_u), + (n_variables, capacity)) + + resize!(neighbor_ids, capacity) + + resize!(faces, capacity) + + resize!(name, capacity) + + return nothing +end + +# Create interface container and initialize interface data in `elements`. +function init_boundaries(mesh::T8codeMesh, equations, solver::FV, uEltype) + # Initialize container + n_boundaries = count_required_surfaces(mesh).boundaries + + _u = Vector{uEltype}(undef, + nvariables(equations) * n_boundaries) + u = unsafe_wrap(Array, pointer(_u), + (nvariables(equations), n_boundaries)) + + neighbor_ids = Vector{Int}(undef, n_boundaries) + faces = Vector{Int}(undef, n_boundaries) + names = Vector{Symbol}(undef, n_boundaries) + + boundaries = T8codeFVBoundaryContainer{uEltype}(u, neighbor_ids, faces, names, _u) + + return boundaries +end + +function init_communication_data!(mesh::T8codeMesh, equations) + (; forest) = mesh + # Check that the forest is a committed. + @assert(t8_forest_is_committed(forest)==1) + + # Get the number of local elements of forest. + num_local_elements = ncells(mesh) + # Get the number of ghost elements of forest. + num_ghost_elements = t8_forest_get_num_ghosts(forest) + + # Build an array of our data that is as long as the number of elements plus + # the number of ghosts. + solution_data = Vector{T8codeSolutionContainer{nvariables(equations)}}(undef, + num_local_elements + + num_ghost_elements) + + domain_data = Vector{T8codeReconstructionContainer{ndims(equations), + mesh.max_number_faces}}(undef, + num_local_elements + + num_ghost_elements) + + gradient_data = Vector{T8codeGradientContainer{ndims(equations), + nvariables(equations)}}(undef, + num_local_elements + + num_ghost_elements) + + return (; solution_data, domain_data, gradient_data) +end + +# Each process has computed the data entries for its local elements. In order +# to get the values for the ghost elements, we use +# t8_forest_ghost_exchange_data. Calling this function will fill all the ghost +# entries of our element data array with the value on the process that owns the +# corresponding element. +function exchange_ghost_data(mesh, container) + # t8_forest_ghost_exchange_data expects an sc_array (of length num_local_elements + num_ghosts). + # We wrap our data array to an sc_array. + sc_array_wrapper = T8code.Libt8.sc_array_new_data(pointer(container), + sizeof(typeof(container[1])), + length(container)) + + # Carry out the data exchange. The entries with indices > num_local_elements will get overwritten. + t8_forest_ghost_exchange_data(mesh.forest, sc_array_wrapper) + + # Destroy the wrapper array. This will not free the data memory since we used sc_array_new_data. + T8code.Libt8.sc_array_destroy(sc_array_wrapper) +end + +struct T8codeSolutionContainer{NVARS} + u::NTuple{NVARS, Cdouble} + + function T8codeSolutionContainer(u) + new{length(u)}(u) + end +end + +function exchange_solution_data!(u, mesh, equations, solver, cache) + (; solution_data) = cache.communication_data + for element in eachelement(mesh, solver, cache) + solution_data[element] = T8codeSolutionContainer(Tuple(get_node_vars(u, + equations, + solver, + element))) + end + exchange_ghost_data(mesh, solution_data) + + return nothing +end + +struct T8codeGradientContainer{NDIMS, NVARS} + reconstruction_gradient_limited::NTuple{NVARS, SVector{NDIMS, Cdouble}} + + function T8codeGradientContainer(reconstruction_gradient_limited) + new{length(reconstruction_gradient_limited[1]), + length(reconstruction_gradient_limited)}(reconstruction_gradient_limited) + end +end + +function exchange_gradient_data!(reconstruction_gradient_limited, + mesh, equations, solver, cache) + (; gradient_data) = cache.communication_data + for element in eachelement(mesh, solver, cache) + gradient_data[element] = T8codeGradientContainer(ntuple(v -> get_node_coords(reconstruction_gradient_limited, + equations, + solver, + v, + element), + Val(nvariables(equations)))) + end + exchange_ghost_data(mesh, gradient_data) + + return nothing +end + +struct T8codeReconstructionContainer{NDIMS, NFACES} + midpoint::SVector{NDIMS, Cdouble} + face_midpoints::NTuple{NFACES, SVector{NDIMS, Cdouble}} + + function T8codeReconstructionContainer(midpoint, face_midpoints) + new{length(midpoint), length(face_midpoints)}(midpoint, face_midpoints) + end +end + +function exchange_domain_data!(communication_data, elements, mesh, equations, solver) + (; domain_data) = communication_data + (; midpoint, face_midpoints, num_faces) = elements + for element in 1:ncells(mesh) + face_midpoints_tuple = [] + for face in 1:(num_faces[element]) + push!(face_midpoints_tuple, + get_node_coords(face_midpoints, equations, solver, face, element)) + end + for i in (num_faces[element] + 1):(mesh.max_number_faces) + push!(face_midpoints_tuple, + SVector{ndims(equations)}(zeros(eltype(midpoint), ndims(equations)))) + end + face_midpoints_tuple_ = NTuple{mesh.max_number_faces, + eltype(face_midpoints_tuple)}(face_midpoints_tuple) + domain_data[element] = T8codeReconstructionContainer(get_node_coords(midpoint, + equations, + solver, + element), + face_midpoints_tuple_) + end + exchange_ghost_data(mesh, domain_data) + + return nothing +end +end # @muladd diff --git a/src/solvers/fv_t8code/fv.jl b/src/solvers/fv_t8code/fv.jl new file mode 100644 index 0000000000..08c82df056 --- /dev/null +++ b/src/solvers/fv_t8code/fv.jl @@ -0,0 +1,639 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +struct FV{RealT <: Real, SurfaceFlux} + order::Integer + extended_reconstruction_stencil::Bool + surface_flux::SurfaceFlux + + function FV(; order = 1, extended_reconstruction_stencil = false, RealT = Float64, + surface_flux = flux_central) + new{RealT, typeof(surface_flux)}(order, extended_reconstruction_stencil, + surface_flux) + end +end + +function Base.show(io::IO, solver::FV) + @nospecialize solver # reduce precompilation time + + print(io, "FV(") + print(io, "order $(solver.order)") + print(io, ", ", solver.surface_flux) + print(io, ")") +end + +function Base.show(io::IO, mime::MIME"text/plain", solver::FV) + @nospecialize solver # reduce precompilation time + + if get(io, :compact, false) + show(io, solver) + else + summary_header(io, "FV{" * string(real(solver)) * "}") + summary_line(io, "order", solver.order) + summary_line(io, "surface flux", solver.surface_flux) + summary_footer(io) + end +end + +Base.summary(io::IO, solver::FV) = print(io, "FV(order=$(solver.order))") + +@inline Base.real(solver::FV{RealT}) where {RealT} = RealT + +@inline ndofs(mesh, solver::FV, cache) = ncells(mesh) + +@inline nelements(mesh::T8codeMesh, solver::FV, cache) = ncells(mesh) +@inline function ndofsglobal(mesh, solver::FV, cache) + nelementsglobal(mesh, solver, cache) +end + +@inline function eachelement(mesh, solver::FV, cache) + Base.OneTo(nelements(mesh, solver, cache)) +end + +@inline eachinterface(solver::FV, cache) = Base.OneTo(ninterfaces(solver, cache)) +@inline eachboundary(solver::FV, cache) = Base.OneTo(nboundaries(solver, cache)) + +@inline function nelementsglobal(mesh, solver::FV, cache) + if mpi_isparallel() + Int(t8_forest_get_global_num_elements(mesh.forest)) + else + nelements(mesh, solver, cache) + end +end + +@inline ninterfaces(solver::FV, cache) = ninterfaces(cache.interfaces) +@inline nboundaries(solver::FV, cache) = nboundaries(cache.boundaries) + +@inline function get_node_coords(x, equations, solver::FV, indices...) + SVector(ntuple(@inline(idx->x[idx, indices...]), Val(ndims(equations)))) +end + +@inline function get_node_vars(u, equations, solver::FV, element) + SVector(ntuple(@inline(v->u[v, element]), Val(nvariables(equations)))) +end + +@inline function set_node_vars!(u, u_node, equations, solver::FV, element) + for v in eachvariable(equations) + u[v, element] = u_node[v] + end + return nothing +end + +@inline function get_surface_node_vars(u, equations, solver::FV, indices...) + # There is a cut-off at `n == 10` inside of the method + # `ntuple(f::F, n::Integer) where F` in Base at ntuple.jl:17 + # in Julia `v1.5`, leading to type instabilities if + # more than ten variables are used. That's why we use + # `Val(...)` below. + u_ll = SVector(ntuple(@inline(v->u[1, v, indices...]), Val(nvariables(equations)))) + u_rr = SVector(ntuple(@inline(v->u[2, v, indices...]), Val(nvariables(equations)))) + return u_ll, u_rr +end + +function allocate_coefficients(mesh::T8codeMesh, equations, solver::FV, cache) + # We must allocate a `Vector` in order to be able to `resize!` it (AMR). + # cf. wrap_array + zeros(eltype(cache.elements), + nvariables(equations) * nelements(mesh, solver, cache)) +end + +# General fallback +@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, + solver::FV, cache) + wrap_array_native(u_ode, mesh, equations, solver, cache) +end + +# Like `wrap_array`, but guarantees to return a plain `Array`, which can be better +# for interfacing with external C libraries (MPI, HDF5, visualization), +# writing solution files etc. +@inline function wrap_array_native(u_ode::AbstractVector, mesh::AbstractMesh, equations, + solver::FV, cache) + @boundscheck begin + @assert length(u_ode) == + nvariables(equations) * nelements(mesh, solver, cache) + end + unsafe_wrap(Array{eltype(u_ode), 2}, pointer(u_ode), + (nvariables(equations), nelements(mesh, solver, cache))) +end + +function compute_coefficients!(u, func, t, mesh::T8codeMesh, + equations, solver::FV, cache) + for element in eachelement(mesh, solver, cache) + x_node = get_node_coords(cache.elements.midpoint, equations, solver, element) + u_node = func(x_node, t, equations) + set_node_vars!(u, u_node, equations, solver, element) + end +end + +function create_cache(mesh::T8codeMesh, equations::AbstractEquations, solver::FV, ::Any, + ::Type{uEltype}) where {uEltype <: Real} + count_required_surfaces!(mesh) + + # After I saved some data (e.g. normal) in the interfaces and boundaries, + # the element data structure is not used anymore after this `create_cache` routine. + # Possible to remove it and directly save the data in interface, boundars (and mortar) data structure? + elements = init_elements(mesh, equations, solver, uEltype) + interfaces = init_interfaces(mesh, equations, solver, uEltype) + boundaries = init_boundaries(mesh, equations, solver, uEltype) + # mortars = init_mortars(mesh, equations, solver, uEltype) + + fill_mesh_info_fv!(mesh, interfaces, boundaries, + mesh.boundary_names) + + # Data structure for exchange between MPI ranks. + communication_data = init_communication_data!(mesh, equations) + exchange_domain_data!(communication_data, elements, mesh, equations, solver) + + # Initialize reconstruction stencil + if !solver.extended_reconstruction_stencil + init_reconstruction_stencil!(elements, interfaces, boundaries, + communication_data, mesh, equations, solver) + end + + cache = (; elements, interfaces, boundaries, communication_data) + + return cache +end + +function rhs!(du, u, t, mesh::T8codeMesh, equations, + initial_condition, boundary_conditions, source_terms::Source, + solver::FV, cache) where {Source} + # Reset du + @trixi_timeit timer() "reset ∂u/∂t" du.=zero(eltype(du)) + + # Exchange solution between MPI ranks + @trixi_timeit timer() "exchange_solution_data!" exchange_solution_data!(u, mesh, + equations, + solver, + cache) + + @trixi_timeit timer() "gradient reconstruction" calc_gradient_reconstruction!(u, + mesh, + equations, + solver, + cache) + + # Prolong solution to interfaces + @trixi_timeit timer() "prolong2interfaces" begin + prolong2interfaces!(cache, mesh, equations, solver) + end + + # Calculate interface fluxes + @trixi_timeit timer() "interface flux" begin + calc_interface_flux!(du, mesh, have_nonconservative_terms(equations), equations, + solver, cache) + end + + # Prolong solution to boundaries + @trixi_timeit timer() "prolong2boundaries!" begin + prolong2boundaries!(cache, mesh, equations, solver) + end + + # Calculate boundary fluxes + @trixi_timeit timer() "calc boundary flux" begin + calc_boundary_flux!(du, cache, t, boundary_conditions, mesh, + equations, solver) + end + + @trixi_timeit timer() "volume" begin + for element in eachelement(mesh, solver, cache) + volume = cache.elements.volume[element] + for v in eachvariable(equations) + du[v, element] = (1 / volume) * du[v, element] + end + end + end + + # Calculate source terms + @trixi_timeit timer() "source terms" begin + calc_sources!(du, u, t, source_terms, mesh, equations, solver, cache) + end + + return nothing +end + +function calc_gradient_reconstruction!(u, mesh, 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 2 Matrix, where N is the number of stencil neighbors + # A^T A 2 x 2 Matrix + # b N Vector + # A^T b 2 Vector + + # Matrix/vector notation + # A^T A = [a1 a2; a2 a3] + # (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] + a = zeros(eltype(u), 3) # [a1, a2, a3] + d = zeros(eltype(u), size(u, 1), 2) # [d1(v), d2(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] + 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)] + # Note: A^T b depends on the variable v. Using vectors with d[v, 1/2] + 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) + else + calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d, + reconstruction_stencil, reconstruction_distance, + solution_data, equations) + end + + # Divide by determinant + AT_A_determinant = a[1] * a[3] - a[2]^2 + if isapprox(AT_A_determinant, 0.0) + for v in eachvariable(equations) + reconstruction_gradient[:, v, stencil, element] .= 0.0 + end + continue + end + a .*= 1 / AT_A_determinant + + # Solving least square problem + for v in eachvariable(equations) + reconstruction_gradient[1, v, stencil, element] = a[3] * d[v, 1] - + a[2] * d[v, 2] + reconstruction_gradient[2, v, stencil, element] = -a[2] * d[v, 1] + + a[1] * d[v, 2] + 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) + 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 + a[2] += distance[1] * distance[2] + a[3] += distance[2]^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]) + end + end + end + + return nothing +end + +@inline function calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d, + reconstruction_stencil, + reconstruction_distance, + reconstruction_corner_elements, solution_data, + equations, 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 + a[2] += distance[1] * distance[2] + a[3] += distance[2]^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]) + 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 + a[2] += distance[1] * distance[2] + a[3] += distance[2]^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]) + 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 + + for interface in eachinterface(solver, cache) + element = interfaces.neighbor_ids[1, interface] + neighbor = interfaces.neighbor_ids[2, interface] + if solver.order == 1 + for v in eachvariable(equations) + interfaces.u[1, v, interface] = solution_data[element].u[v] + interfaces.u[2, v, interface] = solution_data[neighbor].u[v] + end + elseif solver.order == 2 + face_element = interfaces.faces[1, interface] + face_neighbor = interfaces.faces[2, interface] + + face_midpoint_element = domain_data[element].face_midpoints[face_element] + face_midpoint_neighbor = domain_data[neighbor].face_midpoints[face_neighbor] + + midpoint_element = domain_data[element].midpoint + midpoint_neighbor = domain_data[neighbor].midpoint + + vector_element = face_midpoint_element .- midpoint_element + vector_neighbor = face_midpoint_neighbor .- midpoint_neighbor + for v in eachvariable(equations) + gradient_v_element = gradient_data[element].reconstruction_gradient_limited[v] + gradient_v_neighbor = gradient_data[neighbor].reconstruction_gradient_limited[v] + interfaces.u[1, v, interface] = solution_data[element].u[v] + + dot(gradient_v_element, vector_element) + interfaces.u[2, v, interface] = solution_data[neighbor].u[v] + + dot(gradient_v_neighbor, + vector_neighbor) + end + else + error("Order $(solver.order) is not supported.") + end + end + + return nothing +end + +function calc_interface_flux!(du, mesh::T8codeMesh, + nonconservative_terms::False, equations, + solver::FV, cache) + (; surface_flux) = solver + (; elements, interfaces) = cache + + for interface in eachinterface(solver, cache) + element = interfaces.neighbor_ids[1, interface] + neighbor = interfaces.neighbor_ids[2, interface] + face = interfaces.faces[1, interface] + + # TODO: Save normal and face_areas in interface? + normal = get_node_coords(elements.face_normals, equations, solver, + face, element) + u_ll, u_rr = get_surface_node_vars(interfaces.u, equations, solver, + interface) + flux = surface_flux(u_ll, u_rr, normal, equations) + + for v in eachvariable(equations) + flux_ = elements.face_areas[face, element] * flux[v] + du[v, element] -= flux_ + if !is_ghost_cell(neighbor, mesh) + du[v, neighbor] += flux_ + end + end + end + + return nothing +end + +function prolong2boundaries!(cache, mesh::T8codeMesh, equations, solver::FV) + (; elements, boundaries, communication_data) = cache + (; solution_data) = communication_data + (; midpoint, face_midpoints, reconstruction_gradient_limited) = elements + + for boundary in eachboundary(solver, cache) + element = boundaries.neighbor_ids[boundary] + if solver.order == 1 + for v in eachvariable(equations) + boundaries.u[v, boundary] = solution_data[element].u[v] + end + elseif solver.order == 2 + face_element = boundaries.faces[boundary] + + face_midpoint_element = get_node_coords(face_midpoints, equations, solver, + face_element, element) + + midpoint_element = get_node_coords(midpoint, equations, solver, element) + + vector_element = face_midpoint_element .- midpoint_element + for v in eachvariable(equations) + gradient_v_element = get_node_coords(reconstruction_gradient_limited, + equations, solver, v, element) + boundaries.u[v, boundary] = solution_data[element].u[v] + + dot(gradient_v_element, vector_element) + end + else + error("Order $(solver.order) is not supported.") + end + end + + return nothing +end + +function calc_boundary_flux!(du, cache, t, + boundary_condition::BoundaryConditionPeriodic, + mesh::T8codeMesh, + equations, solver::FV) + @assert isempty(eachboundary(solver, cache)) +end + +# Function barrier for type stability +function calc_boundary_flux!(du, cache, t, boundary_conditions, + mesh::T8codeMesh, + equations, solver::FV) + @unpack boundary_condition_types, boundary_indices = boundary_conditions + + calc_boundary_flux_by_type!(du, cache, t, boundary_condition_types, + boundary_indices, mesh, equations, solver) + return nothing +end + +# Iterate over tuples of boundary condition types and associated indices +# in a type-stable way using "lispy tuple programming". +function calc_boundary_flux_by_type!(du, cache, t, BCs::NTuple{N, Any}, + BC_indices::NTuple{N, Vector{Int}}, + mesh::T8codeMesh, + equations, solver::FV) where {N} + # Extract the boundary condition type and index vector + boundary_condition = first(BCs) + boundary_condition_indices = first(BC_indices) + # Extract the remaining types and indices to be processed later + remaining_boundary_conditions = Base.tail(BCs) + remaining_boundary_condition_indices = Base.tail(BC_indices) + + # process the first boundary condition type + calc_boundary_flux!(du, cache, t, boundary_condition, boundary_condition_indices, + mesh, equations, solver) + + # recursively call this method with the unprocessed boundary types + calc_boundary_flux_by_type!(du, cache, t, remaining_boundary_conditions, + remaining_boundary_condition_indices, + mesh, equations, solver) + + return nothing +end + +# terminate the type-stable iteration over tuples +function calc_boundary_flux_by_type!(du, cache, t, BCs::Tuple{}, BC_indices::Tuple{}, + mesh::T8codeMesh, + equations, solver::FV) + nothing +end + +function calc_boundary_flux!(du, cache, t, boundary_condition::BC, boundary_indexing, + mesh::T8codeMesh, + equations, solver::FV) where {BC} + (; elements, boundaries) = cache + (; surface_flux) = solver + + for local_index in eachindex(boundary_indexing) + # Use the local index to get the global boundary index from the pre-sorted list + boundary = boundary_indexing[local_index] + + # Get information on the adjacent element, compute the surface fluxes, + # and store them + element = boundaries.neighbor_ids[boundary] + face = boundaries.faces[boundary] + + # TODO: Save normal and face_areas in interface? + normal = get_node_coords(cache.elements.face_normals, equations, solver, + face, element) + + u_inner = get_node_vars(boundaries.u, equations, solver, boundary) + + # Coordinates at boundary node + face_midpoint = get_node_coords(cache.elements.face_midpoints, equations, + solver, face, element) + + flux = boundary_condition(u_inner, normal, face_midpoint, t, surface_flux, + equations) + for v in eachvariable(equations) + flux_ = elements.face_areas[face, element] * flux[v] + du[v, element] -= flux_ + end + end + + return nothing +end + +function calc_sources!(du, u, t, source_terms::Nothing, mesh::T8codeMesh, + equations::AbstractEquations, solver::FV, cache) + return nothing +end + +function calc_sources!(du, u, t, source_terms, mesh::T8codeMesh, + equations::AbstractEquations, solver::FV, cache) + @threaded for element in eachelement(mesh, solver, cache) + u_local = get_node_vars(u, equations, solver, element) + x_local = get_node_coords(cache.elements.midpoint, equations, solver, element) + du_local = source_terms(u_local, x_local, t, equations) + for v in eachvariable(equations) + du[v, element] += du_local[v] + end + end + + return nothing +end + +function get_element_variables!(element_variables, u, + mesh, equations, + solver::FV, cache) + return nothing +end + +function get_node_variables!(node_variables, mesh, + equations, solver::FV, cache) + return nothing +end + +function SolutionAnalyzer(solver::FV; kwargs...) +end + +function create_cache_analysis(analyzer, mesh, + equations, solver::FV, cache, + RealT, uEltype) +end + +function T8codeMesh(cmesh::Ptr{t8_cmesh}, solver::DG; kwargs...) + T8codeMesh(cmesh; kwargs...) +end + +function T8codeMesh(cmesh::Ptr{t8_cmesh}, solver::FV; kwargs...) + T8codeMesh(cmesh; polydeg = 0, kwargs...) +end + +# Container data structures +include("containers.jl") +end # @muladd diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index a39f7cb175..e2e47467b4 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -10,4 +10,5 @@ include("solvers_parabolic.jl") include("dg.jl") include("dgmulti.jl") +include("fv_t8code/fv.jl") end # @muladd diff --git a/test/runtests.jl b/test/runtests.jl index 836488d0d8..1ecb13c7fa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -93,6 +93,10 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) include("test_t8code_3d.jl") end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "t8code_fv" + include("test_t8code_fv.jl") + end + @time if TRIXI_TEST == "all" || TRIXI_TEST == "unstructured_dgmulti" include("test_unstructured_2d.jl") include("test_dgmulti_1d.jl") diff --git a/test/test_mpi.jl b/test/test_mpi.jl index 001d9bff86..f36b669d3c 100644 --- a/test/test_mpi.jl +++ b/test/test_mpi.jl @@ -23,6 +23,7 @@ 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") if !CI_ON_WINDOWS # see comment on `CI_ON_WINDOWS` above include("test_mpi_p4est_3d.jl") include("test_mpi_t8code_3d.jl") diff --git a/test/test_mpi_t8code_fv.jl b/test/test_mpi_t8code_fv.jl new file mode 100644 index 0000000000..5a2d952555 --- /dev/null +++ b/test/test_mpi_t8code_fv.jl @@ -0,0 +1,360 @@ +module TestExamplesMPIT8codeMesh2D + +using Test +using Trixi + +include("test_trixi.jl") + +const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_fv") + +@testset "T8codeMesh MPI FV" begin +#! format: noindent + +# 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. + # @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 + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_gauss.jl"), + order=1, + l2=[0.5598148317954682], + linf=[0.6301130236005371]) + # 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.5899077806567905], + linf=[0.8972489222157533]) + # 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, + 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"), + order=1, + l2=[ + 0.059376731961878, + 0.019737707470047838, + 0.019737707470047747, + 0.09982550390697936 + ], + linf=[ + 0.08501451493301548, + 0.029105783468157398, + 0.029105783468157842, + 0.1451756151490775 + ]) + # 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.031971635416962525, + 0.016630983283552267, + 0.016630873316960327, + 0.04813244762272231 + ], + linf=[ + 0.055105205670854085, + 0.03647221946045942, + 0.036470504033139894, + 0.0811201478913759 + ]) + # 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_blast_wave.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave.jl"), + order=1, + l2=[ + 0.5733341919395405, + 0.11399976571202448, + 0.1139997657120245, + 1.3548613737038315 + ], + linf=[ + 1.7328363346781415, + 0.27645456051355827, + 0.27645456051355827, + 2.6624886901791407 + ]) + # 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_kelvin_helmholtz_instability.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability.jl"), + order=1, + l2=[ + 0.25420413805862135, + 0.22153054262689362, + 0.11870842058617848, + 0.03626117501911353 + ], + linf=[ + 0.5467894727797227, + 0.4156157752497065, + 0.26176691230685767, + 0.0920609123083227 + ], + tspan=(0.0, 1.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 + @trixi_testset "second-order FV - quads" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability.jl"), + order=2, + cmesh=Trixi.cmesh_new_quad(periodicity = (true, true)), + l2=[ + 0.2307479238046326, + 0.19300139957275295, + 0.1176326315506721, + 0.020439850138732837 + ], + linf=[ + 0.5069212604421109, + 0.365579474379667, + 0.24226411409222004, + 0.049201093470609525 + ], + tspan=(0.0, 1.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 + @trixi_testset "second-order FV - triangles" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability.jl"), + order=2, + l2=[ + 0.16181566308374545, + 0.10090964843765918, + 0.15229553744179888, + 0.0395037796064376 + ], + linf=[ + 0.6484515918189779, + 0.3067327488921227, + 0.34771083375083534, + 0.10713502930441887 + ], + tspan=(0.0, 1.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 +end +end # T8codeMesh MPI + +end # module diff --git a/test/test_t8code_fv.jl b/test/test_t8code_fv.jl new file mode 100644 index 0000000000..9ddef51a9d --- /dev/null +++ b/test/test_t8code_fv.jl @@ -0,0 +1,449 @@ +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_2d_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 "T8codeMesh2D" 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, + 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 + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_gauss.jl"), + order=1, + l2=[0.5598148317954682], + linf=[0.6301130236005371]) + # 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.5899012906928289], + linf=[0.8970164922705812]) + # 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, + 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 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 + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"), + order=1, + l2=[ + 0.059376731961878, + 0.019737707470047838, + 0.019737707470047747, + 0.09982550390697936 + ], + linf=[ + 0.08501451493301548, + 0.029105783468157398, + 0.029105783468157842, + 0.1451756151490775 + ]) + # 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.031971635993647315, + 0.016631028330554957, + 0.016630833188111448, + 0.04813246238398825 + ], + linf=[ + 0.055105654108624336, + 0.03647317645079773, + 0.03647020577993976, + 0.08112180586875883 + ]) + # 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.058781550112786296, + 0.02676026477370241, + 0.02673090935979779, + 0.08033279155463603 + ], + linf=[ + 0.09591263836601072, + 0.05351985245787505, + 0.05264935415308125, + 0.14318629241962988 + ]) + # 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_blast_wave.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_blast_wave.jl"), + order=1, + l2=[ + 0.5733341919395403, + 0.11399976571202451, + 0.11399976571202453, + 1.3548613737038324 + ], + linf=[ + 1.732836334678142, + 0.27645456051355827, + 0.27645456051355827, + 2.6624886901791407 + ]) + # 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_blast_wave.jl"), + l2=[ + 0.7331398527938754, + 0.15194349346989244, + 0.1519434934698924, + 1.299914264830515 + ], + linf=[ + 2.2864127304524726, + 0.3051023693829176, + 0.30510236938291757, + 2.6171402581107936 + ]) + # 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_kelvin_helmholtz_instability.jl" begin + @trixi_testset "first-order FV" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability.jl"), + order=1, + l2=[ + 0.25420413805862135, + 0.22153054262689362, + 0.11870842058617848, + 0.03626117501911353 + ], + linf=[ + 0.5467894727797227, + 0.4156157752497065, + 0.26176691230685767, + 0.0920609123083227 + ], + tspan=(0.0, 1.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 + @trixi_testset "second-order FV - quads" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability.jl"), + order=2, + cmesh=Trixi.cmesh_new_quad(periodicity = (true, true)), + l2=[ + 0.2307479238046326, + 0.19300139957275295, + 0.1176326315506721, + 0.020439850138732837 + ], + linf=[ + 0.5069212604421109, + 0.365579474379667, + 0.24226411409222004, + 0.049201093470609525 + ], + tspan=(0.0, 1.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 + @trixi_testset "second-order FV - triangles" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_kelvin_helmholtz_instability.jl"), + order=2, + l2=[ + 0.16181566308374545, + 0.10090964843765918, + 0.15229553744179888, + 0.0395037796064376 + ], + linf=[ + 0.6484515918189779, + 0.3067327488921227, + 0.34771083375083534, + 0.10713502930441887 + ], + tspan=(0.0, 1.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 +end + +# Clean up afterwards: delete Trixi.jl output directory +@test_nowarn rm(outdir, recursive = true) + +end # module