Add 3d support for FV using T8codeMesh (#130)

* Add first version of 3d support

* Complete 3d support

* Add 3d elixirs and tests

* Adapt name in elixir

* Reduce initial refinement level of one test

* Simplify NDIMS=3

* Add 3d mpi tests

* Remove initial_condition from rhs; fmt

* Fix parameter in `x_trans_periodic_3d`

* Remove comment

* Some last changes

* fmt

examples/t8code_2d_fv/elixir_advection_basic.jl

@@ -16,6 +16,8 @@ solver = FV(order = 2, extended_reconstruction_stencil = false,
 # 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.

examples/t8code_2d_fv/elixir_advection_gauss.jl

@@ -11,9 +11,9 @@ initial_condition = initial_condition_gauss
 solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
 
 initial_refinement_level = 4
+# Note: The initial_condition is set up for a domain [-5,5]^2.
+# This becomes important when using a non-periodic mesh.
 cmesh = Trixi.cmesh_new_periodic_hybrid()
-# Note: A non-periodic run with the tri mesh is unstable. Same as in `elixir_advection_nonperiodic.jl`
-# cmesh = Trixi.cmesh_new_tri(periodicity = (false, false))
 
 mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
 

examples/t8code_2d_fv/elixir_advection_nonperiodic.jl

@@ -23,7 +23,7 @@ solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
 initial_refinement_level = 5
 cmesh = Trixi.cmesh_new_quad(periodicity = (false, false))
 
-# **Note**: A non-periodic run with the tri mesh is unstable. (Same happens for a non-periodic run with `elixir_advection_gauss.jl`)
+# **Note**: A non-periodic run with the tri mesh is unstable.
 # - This only happens with **2nd order**
 # - When increasing refinement_level to 6 and lower CFL to 0.4, it seems like the simulation is stable again (except of some smaller noises at the corners)
 # - With a lower resolution (5) I cannot get the simulation stable. Even with a cfl of 0.01 etc.

examples/t8code_3d_fv/elixir_advection_basic.jl

@@ -0,0 +1,96 @@
+using OrdinaryDiffEq
+using Trixi
+using T8code
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(order = 2, extended_reconstruction_stencil = false,
+            surface_flux = flux_lax_friedrichs)
+
+# 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 trixi_t8_mapping(cmesh, gtreeid, ref_coords, num_coords, out_coords,
+                          tree_data, user_data)
+    ltreeid = t8_cmesh_get_local_id(cmesh, gtreeid)
+    eclass = t8_cmesh_get_tree_class(cmesh, ltreeid)
+    T8code.t8_geom_compute_linear_geometry(eclass, tree_data,
+                                           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
+
+function trixi_t8_mapping_c()
+    @cfunction($trixi_t8_mapping, Cvoid,
+               (t8_cmesh_t, t8_gloidx_t, Ptr{Cdouble}, Csize_t,
+                Ptr{Cdouble}, Ptr{Cvoid}, Ptr{Cvoid}))
+end
+
+trees_per_dimension = (2, 2, 2)
+
+eclass = T8_ECLASS_HEX
+mesh = T8codeMesh(trees_per_dimension, eclass;
+                  mapping = trixi_t8_mapping_c(),
+                  # Plan is to use either
+                  # coordinates_max = coordinates_max, coordinates_min = coordinates_min,
+                  # or at least
+                  # mapping = mapping,
+                  initial_refinement_level = 5)
+
+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()

examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl

@@ -0,0 +1,42 @@
+using OrdinaryDiffEq
+using Trixi
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(order = 2, extended_reconstruction_stencil = false,
+            surface_flux = flux_lax_friedrichs)
+
+# TODO: There are no other cmesh functions implemented yet in 3d.
+cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true))
+mesh = T8codeMesh(cmesh, solver,
+                  initial_refinement_level = 4)
+
+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 = 1,
+                                     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()

examples/t8code_3d_fv/elixir_advection_gauss.jl

@@ -0,0 +1,47 @@
+using OrdinaryDiffEq
+using Trixi
+
+####################################################
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_gauss
+
+solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
+
+initial_refinement_level = 4
+
+# TODO: There are no other cmesh functions implemented yet in 3d.
+cmesh = Trixi.cmesh_new_quad_3d(coordinates_min = (-5.0, -5.0, -5.0),
+                                coordinates_max = (5.0, 5.0, 5.0),
+                                periodicity = (true, true, true))
+mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.7)
+
+callbacks = CallbackSet(summary_callback, save_solution, analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),#Euler(),
+            dt = 5.0e-2, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback()

examples/t8code_3d_fv/elixir_advection_nonperiodic.jl

@@ -0,0 +1,53 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the linear advection equation.
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:all => boundary_condition)
+# Problem: T8codeMesh interface with parameter cmesh cannot distinguish between boundaries
+# boundary_conditions = Dict(:x_neg => boundary_condition,
+#                            :x_pos => boundary_condition,
+#                            :y_neg => boundary_condition_periodic,
+#                            :y_pos => boundary_condition_periodic)
+
+solver = FV(order = 2, surface_flux = flux_lax_friedrichs)
+
+# TODO: When using mesh construction as in elixir_advection_basic.jl boundary Symbol :all is not defined
+initial_refinement_level = 5
+cmesh = Trixi.cmesh_new_quad_3d(periodicity = (false, false, false))
+mesh = T8codeMesh(cmesh, solver, initial_refinement_level = initial_refinement_level)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# Run the simulation.
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # Solve needs some value here but it will be overwritten by the stepsize_callback.
+            save_everystep = false, callback = callbacks);
+summary_callback()

examples/t8code_3d_fv/elixir_euler_source_terms.jl

@@ -0,0 +1,51 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+equations = CompressibleEulerEquations3D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+# boundary_condition = BoundaryConditionDirichlet(initial_condition)
+# boundary_conditions = Dict(:all => boundary_condition)
+
+source_terms = source_terms_convergence_test
+
+solver = FV(order = 2, extended_reconstruction_stencil = false,
+            surface_flux = flux_lax_friedrichs)
+
+# TODO: There are no other cmesh functions implemented yet in 3d.
+cmesh = Trixi.cmesh_new_quad_3d(periodicity = (true, true, true))
+mesh = T8codeMesh(cmesh, solver,
+                  initial_refinement_level = 3)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms)
+# boundary_conditions = boundary_conditions)
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

src/equations/linear_scalar_advection_3d.jl

@@ -64,14 +64,24 @@ function initial_condition_convergence_test(x, t,
     return SVector(scalar)
 end
 
+# Calculates translated coordinates `x` for a periodic domain
+function x_trans_periodic_3d(x, domain_length = SVector(10, 10, 10),
+                             center = SVector(0, 0, 0))
+    x_normalized = x .- center
+    x_shifted = x_normalized .% domain_length
+    x_offset = ((x_shifted .< -0.5f0 * domain_length) -
+                (x_shifted .> 0.5f0 * domain_length)) .* domain_length
+    return center + x_shifted + x_offset
+end
+
 """
-    initial_condition_gauss(x, t, equations::LinearScalarAdvectionEquation1D)
+    initial_condition_gauss(x, t, equations::LinearScalarAdvectionEquation3D)
 
 A Gaussian pulse.
 """
 function initial_condition_gauss(x, t, equation::LinearScalarAdvectionEquation3D)
     # Store translated coordinate for easy use of exact solution
-    x_trans = x - equation.advection_velocity * t
+    x_trans = x_trans_periodic_3d(x - equation.advection_velocity * t)
 
     scalar = exp(-(x_trans[1]^2 + x_trans[2]^2 + x_trans[3]^2))
     return SVector(scalar)