Merge branch 't8codemesh-fv' into t8codemesh-fv-3d

Project.toml

@@ -103,7 +103,7 @@ StaticArrays = "1.5"
 StrideArrays = "0.1.26"
 StructArrays = "0.6.11"
 SummationByPartsOperators = "0.5.41"
-T8code = "0.5"
+T8code = "0.7"
 TimerOutputs = "0.5.7"
 Triangulate = "2.2"
 TriplotBase = "0.1"

examples/t8code_2d_fv/elixir_advection_basic.jl

@@ -1,5 +1,6 @@
 using OrdinaryDiffEq
 using Trixi
+using T8code
 
 advection_velocity = (0.2, -0.7)
 equations = LinearScalarAdvectionEquation2D(advection_velocity)
@@ -9,14 +10,73 @@ initial_condition = initial_condition_convergence_test
 solver = FV(order = 2, extended_reconstruction_stencil = false,
             surface_flux = flux_lax_friedrichs)
 
-coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
-coordinates_max = (1.0, 1.0) # maximum coordinates (max(x), max(y))
-
-trees_per_dimension = (8, 8)
-
-mesh = T8codeMesh(trees_per_dimension,
-                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
-                  initial_refinement_level = 4)
+# Note:
+# For now, it is completely irrelevant that coordinates_max/min are.
+# The used t8code routine creates the mesh on [0, nx] x [0, ny], where (nx, ny) = trees_per_dimension.
+# Afterwards and only inside Trixi, `tree_node_coordinates` are mapped back to [-1, 1]^2.
+# But, this variable is not used for the FV method.
+# That's why I use the cmesh interface in all other elixirs.
+
+# Option 1: coordinates
+coordinates_min = (0.0, 0.0) # minimum coordinates (min(x), min(y))
+coordinates_max = (8.0, 8.0) # maximum coordinates (max(x), max(y))
+# Note and TODO: The plan is to move the auxiliary routine f and the macro to a different place.
+# Then, somehow, I get SegFaults when using this `mapping_coordinates` or (equally) when
+# using `coordinates_min/max` and then use the `coordinates2mapping` within the constructor.
+# With both other mappings I don't get that.
+mapping_coordinates = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
+
+# Option 2: faces
+f1(s) = SVector(-1.0, s - 1.0)
+f2(s) = SVector(1.0, s + 1.0)
+f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
+f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s))
+faces = (f1, f2, f3, f4)
+mapping_faces = Trixi.transfinite_mapping(faces)
+
+# Option 3: classic mapping
+function mapping(xi, eta)
+    x = 4 * (xi + 1)
+    y = 4 * (eta + 1)
+
+    return SVector(x, y)
+end
+
+# Note and TODO:
+# Normally, this should be put somewhere else. For now, that doesn't properly.
+# See note in `src/auxiliary/t8code.jl`
+function f(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data)
+    ltreeid = t8_cmesh_get_local_id(cmesh, gtreeid)
+    eclass = t8_cmesh_get_tree_class(cmesh, ltreeid)
+    T8code.t8_geom_compute_linear_geometry(eclass, tree_data,
+                                           ref_coords, num_coords, out_coords)
+
+    for i in 1:num_coords
+        offset_3d = 3 * (i - 1) + 1
+
+        xi = unsafe_load(out_coords, offset_3d)
+        eta = unsafe_load(out_coords, offset_3d + 1)
+        # xy = mapping_coordinates(xi, eta)
+        # xy = mapping_faces(xi, eta)
+        xy = mapping(xi, eta)
+
+        unsafe_store!(out_coords, xy[1], offset_3d)
+        unsafe_store!(out_coords, xy[2], offset_3d + 1)
+    end
+
+    return nothing
+end
+
+trees_per_dimension = (2, 2)
+
+eclass = T8_ECLASS_QUAD
+mesh = T8codeMesh(trees_per_dimension, eclass,
+                  mapping = @t8_analytic_callback(f),
+                  # Plan is to use either
+                  # coordinates_max = coordinates_max, coordinates_min = coordinates_min,
+                  # or at least
+                  # mapping = mapping,
+                  initial_refinement_level = 6)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 

examples/t8code_2d_fv/elixir_advection_basic_hybrid.jl

@@ -10,8 +10,8 @@ solver = FV(order = 2, extended_reconstruction_stencil = false,
             surface_flux = flux_lax_friedrichs)
 
 cmesh = Trixi.cmesh_new_periodic_hybrid()
-# cmesh = Trixi.cmesh_quad(periodicity = (true, true))
-# cmesh = Trixi.cmesh_new_periodic_tri()
+# cmesh = Trixi.cmesh_new_quad(periodicity = (true, true))
+# cmesh = Trixi.cmesh_new_tri(periodicity = (true, true))
 mesh = T8codeMesh(cmesh, solver,
                   initial_refinement_level = 3)
 

examples/t8code_2d_fv/elixir_advection_gauss.jl

@@ -12,6 +12,9 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
 
 initial_refinement_level = 4
 cmesh = Trixi.cmesh_new_periodic_hybrid()
+# Note: A non-periodic run with the tri mesh is unstable. Same as in `elixir_advection_nonperiodic.jl`
+# cmesh = Trixi.cmesh_new_tri(periodicity = (false, false))
+
 mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/t8code_2d_fv/elixir_advection_nonperiodic.jl

@@ -21,7 +21,16 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
 
 # TODO: When using mesh construction as in elixir_advection_basic.jl boundary Symbol :all is not defined
 initial_refinement_level = 5
-cmesh = Trixi.cmesh_quad(periodicity = (false, false))
+cmesh = Trixi.cmesh_new_quad(periodicity = (false, false))
+
+# **Note**: A non-periodic run with the tri mesh is unstable. (Same happens for a non-periodic run with `elixir_advection_gauss.jl`)
+# - This only happens with **2nd order**
+# - When increasing refinement_level to 6 and lower CFL to 0.4, it seems like the simulation is stable again (except of some smaller noises at the corners)
+# - With a lower resolution (5) I cannot get the simulation stable. Even with a cfl of 0.01 etc.
+# -> That can't be expected.
+# -> Maybe, the reconstruction is just not fitted for this example/mesh/resolution?!
+# cmesh = Trixi.cmesh_new_tri(periodicity = (false, false))
+
 mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
@@ -39,10 +48,13 @@ analysis_interval = 100
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = cons2prim)
+
 stepsize_callback = StepsizeCallback(cfl = 0.8)
 
 callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
-                        stepsize_callback)
+                        save_solution, stepsize_callback)
 
 ###############################################################################
 # Run the simulation.

examples/t8code_2d_fv/elixir_euler_blast_wave.jl

@@ -30,7 +30,7 @@ solver = FV(order = 2, extended_reconstruction_stencil = true,
 
 initial_refinement_level = 4
 # cmesh = Trixi.cmesh_new_periodic_hybrid2()
-cmesh = Trixi.cmesh_quad(periodicity = (true, true))
+cmesh = Trixi.cmesh_new_quad(periodicity = (true, true))
 mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/t8code_2d_fv/elixir_euler_free_stream.jl

@@ -32,7 +32,8 @@ end
 mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/a075f8ec39a67fa9fad8f6f84342cbca/raw/a7206a02ed3a5d3cadacd8d9694ac154f9151db7/square_unstructured_1.inp",
                            joinpath(@__DIR__, "square_unstructured_1.inp"))
 
-# TODO: Including mapping does not work yet for FV.
+# Note:
+# Including mapping does not work yet for FV.
 mesh = T8codeMesh(mesh_file, 2; polydeg = 0,
                   mapping = mapping,
                   initial_refinement_level = 2)

examples/t8code_2d_fv/elixir_euler_kelvin_helmholtz_instability.jl

@@ -36,8 +36,8 @@ solver = FV(order = 2, extended_reconstruction_stencil = false,
 
 initial_refinement_level = 4
 # cmesh = Trixi.cmesh_new_periodic_hybrid()
-# cmesh = Trixi.cmesh_quad(periodicity = (true, true))
-# cmesh = Trixi.cmesh_new_periodic_tri()
+# cmesh = Trixi.cmesh_new_quad(periodicity = (true, true))
+# cmesh = Trixi.cmesh_new_tri(periodicity = (true, true))
 cmesh = Trixi.cmesh_new_periodic_tri2()
 mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
 

examples/t8code_2d_fv/elixir_euler_source_terms.jl

@@ -15,8 +15,8 @@ solver = FV(order = 2, extended_reconstruction_stencil = false,
             surface_flux = flux_lax_friedrichs)
 
 cmesh = Trixi.cmesh_new_periodic_hybrid()
-# cmesh = Trixi.cmesh_quad(periodicity = (true, true))
-# cmesh = Trixi.cmesh_new_periodic_tri()
+# cmesh = Trixi.cmesh_new_quad(periodicity = (true, true))
+# cmesh = Trixi.cmesh_new_tri(periodicity = (true, true))
 mesh = T8codeMesh(cmesh, solver,
                   initial_refinement_level = 3)
 

examples/t8code_3d_fv/elixir_advection_basic.jl

@@ -1,5 +1,6 @@
 using OrdinaryDiffEq
 using Trixi
+using T8code
 
 advection_velocity = (0.2, -0.7, 0.5)
 equations = LinearScalarAdvectionEquation3D(advection_velocity)
@@ -9,14 +10,55 @@ initial_condition = initial_condition_convergence_test
 solver = FV(order = 2, extended_reconstruction_stencil = false,
             surface_flux = flux_lax_friedrichs)
 
-coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z))
-coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z))
-
-trees_per_dimension = (8, 8, 8)
-
-mesh = T8codeMesh(trees_per_dimension,
-                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
-                  initial_refinement_level = 3)
+# Option 1: coordinates
+# For all problems see the 2D file...
+coordinates_min = (0.0, 0.0, 0.0) # minimum coordinates (min(x), min(y), min(z))
+coordinates_max = (8.0, 8.0, 8.0) # maximum coordinates (max(x), max(y), max(z))
+
+mapping_coordinates = Trixi.coordinates2mapping(coordinates_min, coordinates_max)
+
+# Option 3: classic mapping
+function mapping(xi, eta, zeta)
+    x = 4 * (xi + 1)
+    y = 4 * (eta + 1)
+    z = 4 * (zeta + 1)
+
+    return SVector(x, y, z)
+end
+
+function f(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data)
+    ltreeid = t8_cmesh_get_local_id(cmesh, gtreeid)
+    eclass = t8_cmesh_get_tree_class(cmesh, ltreeid)
+    T8code.t8_geom_compute_linear_geometry(eclass, tree_data,
+                                           ref_coords, num_coords, out_coords)
+
+    for i in 1:num_coords
+        offset_3d = 3 * (i - 1) + 1
+
+        xi = unsafe_load(out_coords, offset_3d)
+        eta = unsafe_load(out_coords, offset_3d + 1)
+        zeta = unsafe_load(out_coords, offset_3d + 2)
+        # xyz = mapping_coordinates(xi, eta, zeta)
+        xyz = mapping(xi, eta, zeta)
+
+        unsafe_store!(out_coords, xyz[1], offset_3d)
+        unsafe_store!(out_coords, xyz[2], offset_3d + 1)
+        unsafe_store!(out_coords, xyz[3], offset_3d + 2)
+    end
+
+    return nothing
+end
+
+trees_per_dimension = (1, 1, 1)
+
+eclass = T8_ECLASS_HEX
+mesh = T8codeMesh(trees_per_dimension, eclass;
+                  mapping = @t8_analytic_callback(f),
+                  # Plan is to use either
+                  # coordinates_max = coordinates_max, coordinates_min = coordinates_min,
+                  # or at least
+                  # mapping = mapping,
+                  initial_refinement_level = 6)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 

examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl

@@ -10,7 +10,7 @@ solver = FV(order = 2, extended_reconstruction_stencil = false,
             surface_flux = flux_lax_friedrichs)
 
 # TODO: There are no other cmesh functions implemented yet in 3d.
-cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true))
+cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true))
 mesh = T8codeMesh(cmesh, solver,
                   initial_refinement_level = 4)
 

examples/t8code_3d_fv/elixir_advection_gauss.jl

@@ -13,7 +13,7 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
 initial_refinement_level = 4
 
 # TODO: There are no other cmesh functions implemented yet in 3d.
-cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true))
+cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true))
 mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

examples/t8code_3d_fv/elixir_advection_nonperiodic.jl

@@ -21,7 +21,7 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
 
 # TODO: When using mesh construction as in elixir_advection_basic.jl boundary Symbol :all is not defined
 initial_refinement_level = 5
-cmesh = Trixi.cmesh_quad_3d(periodicity = (false, false, false))
+cmesh = Trixi.cmesh_new_quad_3d(periodicity = (false, false, false))
 mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,

examples/t8code_3d_fv/elixir_euler_source_terms.jl

@@ -15,7 +15,7 @@ solver = FV(order = 2, extended_reconstruction_stencil = false,
             surface_flux = flux_lax_friedrichs)
 
 # TODO: There are no other cmesh functions implemented yet in 3d.
-cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true))
+cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true))
 mesh = T8codeMesh(cmesh, solver,
                   initial_refinement_level = 3)
 

src/auxiliary/t8code.jl

@@ -190,3 +190,34 @@ function trixi_t8_adapt!(mesh, indicators)
 
     return differences
 end
+
+# Note and TODO:
+# This routine seems to work for most of the mappings.
+# Somehow, when using `coordinates_min/max` and then `coordinates2mapping`,
+# I get SegFault errors.
+# This happens when running this in a new julia session or after some runs when I ran simulations
+# with other mappings before.
+# Cannot figure out why. For now, I will leave this auxiliary mapping within the elixir
+# and comment this one out.
+# function trixi_t8_mapping_c(mapping)
+#     function f(cmesh, gtreeid, ref_coords, num_coords, out_coords, tree_data, user_data)
+#         ltreeid = t8_cmesh_get_local_id(cmesh, gtreeid)
+#         eclass = t8_cmesh_get_tree_class(cmesh, ltreeid)
+#         T8code.t8_geom_compute_linear_geometry(eclass, tree_data, ref_coords, num_coords, out_coords)
+
+#         for i in 1:num_coords
+#             offset_3d = 3 * (i - 1) + 1
+
+#             xi = unsafe_load(out_coords, offset_3d)
+#             eta = unsafe_load(out_coords, offset_3d + 1)
+#             xy = mapping(xi, eta)
+
+#             unsafe_store!(out_coords, xy[1], offset_3d)
+#             unsafe_store!(out_coords, xy[2], offset_3d + 1)
+#         end
+
+#         return nothing
+#     end
+
+#     return @cfunction($f, Cvoid, (t8_cmesh_t, t8_gloidx_t, Ptr{Cdouble}, Csize_t, Ptr{Cdouble}, Ptr{Cvoid}, Ptr{Cvoid}))
+# end