Complete 3d support

examples/t8code_3d_fv/elixir_advection_basic.jl

@@ -0,0 +1,47 @@
+using OrdinaryDiffEq
+using Trixi
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(order = 2, extended_reconstruction_stencil = false,
+            surface_flux = flux_lax_friedrichs)
+
+coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z))
+coordinates_max = (1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z))
+
+trees_per_dimension = (8, 8, 8)
+
+mesh = T8codeMesh(trees_per_dimension,
+                  coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                  initial_refinement_level = 3)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.9)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback()

examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl

@@ -13,11 +13,11 @@ solver = FV(order = 2, extended_reconstruction_stencil = false,
 cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true))
 # cmesh = Trixi.cmesh_new_periodic_tri()
 mesh = T8codeMesh(cmesh, solver,
-                  initial_refinement_level = 5)
+                  initial_refinement_level = 3)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
-ode = semidiscretize(semi, (0.0, 10.0));
+ode = semidiscretize(semi, (0.0, 1.0));
 
 summary_callback = SummaryCallback()
 
@@ -29,7 +29,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = 1,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.3)
+stepsize_callback = StepsizeCallback(cfl = 0.9)
 
 callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
                         stepsize_callback, save_solution)

examples/t8code_3d_fv/elixir_advection_gauss.jl

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

src/equations/linear_scalar_advection_3d.jl

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

src/meshes/t8code_mesh.jl

@@ -1597,11 +1597,12 @@ function cmesh_quad_3d(; comm = mpi_comm(), periodicity = (true, true, true))::t
     linear_geom = t8_geometry_linear_new(n_dims)
 
     # This is how the cmesh looks like. The numbers are the tree numbers:
-    #
-    #   +---+
-    #   |   |
-    #   | 0 |
-    #   |   |
+    #     +---+
+    #    /   /|
+    #   +---+ |
+    #   |   | |
+    #   | 0 | +
+    #   |   |/
     #   +---+
     #
 

src/solvers/fv_t8code/containers.jl

@@ -371,7 +371,6 @@ function init_reconstruction_stencil!(elements, interfaces, boundaries,
     return nothing
 end
 
-
 # Container data structure (structure-of-arrays style) for FV interfaces
 mutable struct T8codeFVInterfaceContainer{uEltype <: Real} <: AbstractContainer
     u::Array{uEltype, 3}                # [primary/secondary, variable, interface]

src/solvers/fv_t8code/fv.jl

@@ -248,9 +248,12 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver
     epsilon = 1.0e-13
 
     for element in eachelement(mesh, solver, cache)
+        # The actual number of used stencils is `n_stencil_neighbors + 1`, since the full stencil is additionally used once.
         if solver.extended_reconstruction_stencil
+            # Number of faces = Number of corners
             n_stencil_neighbors = elements.num_faces[element]
         else
+            # Number of direct (edge) neighbors
             n_stencil_neighbors = length(reconstruction_stencil[element])
         end
 
@@ -262,7 +265,7 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver
             d .= zero(eltype(u))
 
             if solver.extended_reconstruction_stencil
-                calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d,
+                calc_gradient_extended!(stencil, element, a, d,
                                         reconstruction_stencil, reconstruction_distance,
                                         reconstruction_corner_elements, solution_data,
                                         equations, solver, cache)
@@ -302,14 +305,14 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver
                 indicator = sum(gradient .^ 2)
                 lambda_j = (j == 1) ? lambda[1] : lambda[2]
                 weight = (lambda_j / (indicator + epsilon))^r
-                for x in axes(reconstruction_gradient_limited, 1)
-                    reconstruction_gradient_limited[x, v, element] += weight *
-                                                                      gradient[x]
+                for dim in axes(reconstruction_gradient_limited, 1)
+                    reconstruction_gradient_limited[dim, v, element] += weight *
+                                                                        gradient[dim]
                 end
                 weight_sum += weight
             end
-            for x in axes(reconstruction_gradient_limited, 1)
-                reconstruction_gradient_limited[x, v, element] /= weight_sum
+            for dim in axes(reconstruction_gradient_limited, 1)
+                reconstruction_gradient_limited[dim, v, element] /= weight_sum
             end
         end
     end
@@ -347,11 +350,11 @@ end
     return nothing
 end
 
-@inline function calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d,
+@inline function calc_gradient_extended!(stencil, element, a, d,
                                          reconstruction_stencil,
                                          reconstruction_distance,
                                          reconstruction_corner_elements, solution_data,
-                                         equations, solver, cache)
+                                         equations::AbstractEquations{2}, solver, cache)
     if stencil == 1
         # Full stencil
         for i in eachindex(reconstruction_stencil[element])
@@ -432,10 +435,12 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver
     epsilon = 1.0e-13
 
     for element in eachelement(mesh, solver, cache)
+        # The actual number of used stencils is `n_stencil_neighbors + 1`, since the full stencil is additionally used once.
         if solver.extended_reconstruction_stencil
+            # Number of faces = Number of corners
             n_stencil_neighbors = elements.num_faces[element]
-            error("not yet supported!")
         else
+            # Number of direct (surface) neighbors
             n_stencil_neighbors = length(reconstruction_stencil[element])
         end
 
@@ -447,11 +452,10 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver
             d .= zero(eltype(u))
 
             if solver.extended_reconstruction_stencil
-                calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d,
+                calc_gradient_extended!(stencil, element, a, d,
                                         reconstruction_stencil, reconstruction_distance,
                                         reconstruction_corner_elements, solution_data,
                                         equations, solver, cache)
-                error("Not supported yet")
             else
                 calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d,
                                       reconstruction_stencil, reconstruction_distance,
@@ -460,7 +464,9 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver
 
             # Divide by determinant
             # Determinant    = a11  * (a22  * a33  -  a23^2 ) - a12  * (a12  * a33  - a13  * a23 ) + a13  * (a12  * a23  - a22  * a13 )
-            AT_A_determinant = a[1] * (a[4] * a[6] -  a[5]^2) - a[2] * (a[2] * a[6] - a[3] * a[5]) + a[3] * (a[2] * a[5] - a[4] * a[3])
+            AT_A_determinant = a[1] * (a[4] * a[6] - a[5]^2) -
+                               a[2] * (a[2] * a[6] - a[3] * a[5]) +
+                               a[3] * (a[2] * a[5] - a[4] * a[3])
             if isapprox(AT_A_determinant, 0.0)
                 for v in eachvariable(equations)
                     reconstruction_gradient[:, v, stencil, element] .= 0.0
@@ -482,9 +488,18 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver
 
             # Solving least square problem
             for v in eachvariable(equations)
-                reconstruction_gradient[1, v, stencil, element] = 1 / AT_A_determinant * (m11 * d[v, 1] + m12 * d[v, 2] + m13 * d[v, 3])
-                reconstruction_gradient[2, v, stencil, element] = 1 / AT_A_determinant * (m12 * d[v, 1] + m22 * d[v, 2] + m23 * d[v, 3])
-                reconstruction_gradient[3, v, stencil, element] = 1 / AT_A_determinant * (m13 * d[v, 1] + m23 * d[v, 2] + m33 * d[v, 3])
+                reconstruction_gradient[1, v, stencil, element] = 1 / AT_A_determinant *
+                                                                  (m11 * d[v, 1] +
+                                                                   m12 * d[v, 2] +
+                                                                   m13 * d[v, 3])
+                reconstruction_gradient[2, v, stencil, element] = 1 / AT_A_determinant *
+                                                                  (m12 * d[v, 1] +
+                                                                   m22 * d[v, 2] +
+                                                                   m23 * d[v, 3])
+                reconstruction_gradient[3, v, stencil, element] = 1 / AT_A_determinant *
+                                                                  (m13 * d[v, 1] +
+                                                                   m23 * d[v, 2] +
+                                                                   m33 * d[v, 3])
             end
         end
 
@@ -499,14 +514,14 @@ function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver
                 indicator = sum(gradient .^ 2)
                 lambda_j = (j == 1) ? lambda[1] : lambda[2]
                 weight = (lambda_j / (indicator + epsilon))^r
-                for x in axes(reconstruction_gradient_limited, 1)
-                    reconstruction_gradient_limited[x, v, element] += weight *
-                                                                      gradient[x]
+                for dim in axes(reconstruction_gradient_limited, 1)
+                    reconstruction_gradient_limited[dim, v, element] += weight *
+                                                                        gradient[dim]
                 end
                 weight_sum += weight
             end
-            for x in axes(reconstruction_gradient_limited, 1)
-                reconstruction_gradient_limited[x, v, element] /= weight_sum
+            for dim in axes(reconstruction_gradient_limited, 1)
+                reconstruction_gradient_limited[dim, v, element] /= weight_sum
             end
         end
     end
@@ -550,6 +565,68 @@ end
     return nothing
 end
 
+@inline function calc_gradient_extended!(stencil, element, a, d,
+                                         reconstruction_stencil,
+                                         reconstruction_distance,
+                                         reconstruction_corner_elements, solution_data,
+                                         equations::AbstractEquations{3}, solver, cache)
+    if stencil == 1
+        # Full stencil
+        for i in eachindex(reconstruction_stencil[element])
+            neighbor = reconstruction_stencil[element][i]
+            distance = reconstruction_distance[element][i]
+            a[1] += distance[1]^2 # = a11
+            a[2] += distance[1] * distance[2] # = a12
+            a[3] += distance[1] * distance[3]
+            a[4] += distance[2]^2
+            a[5] += distance[2] * distance[3]
+            a[6] += distance[3]^2
+
+            for v in eachvariable(equations)
+                d[v, 1] += distance[1] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+                d[v, 2] += distance[2] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+                d[v, 3] += distance[3] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+            end
+        end
+    else
+        # Partial stencils
+        midpoint_element = get_node_coords(cache.elements.midpoint, equations, solver,
+                                           element)
+        for neighbor in reconstruction_corner_elements[stencil - 1, element]
+            midpoint_neighbor = get_node_coords(cache.elements.midpoint, equations,
+                                                solver, neighbor)
+            distance = midpoint_neighbor .- midpoint_element
+
+            a[1] += distance[1]^2 # = a11
+            a[2] += distance[1] * distance[2] # = a12
+            a[3] += distance[1] * distance[3]
+            a[4] += distance[2]^2
+            a[5] += distance[2] * distance[3]
+            a[6] += distance[3]^2
+
+            for v in eachvariable(equations)
+                d[v, 1] += distance[1] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+                d[v, 2] += distance[2] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+                d[v, 3] += distance[3] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+            end
+        end
+    end
+
+    return nothing
+end
+
 function prolong2interfaces!(cache, mesh::T8codeMesh, equations, solver::FV)
     (; interfaces, communication_data) = cache
     (; solution_data, domain_data, gradient_data) = communication_data