From df555d212b10d7ae7c741b52efc531b203162c4f Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 19 Aug 2024 15:54:38 +0200 Subject: [PATCH] Complete 3d support --- .../t8code_3d_fv/elixir_advection_basic.jl | 47 +++++++ .../elixir_advection_basic_hybrid.jl | 6 +- .../t8code_3d_fv/elixir_advection_gauss.jl | 44 +++++++ src/equations/linear_scalar_advection_3d.jl | 14 ++- src/meshes/t8code_mesh.jl | 11 +- src/solvers/fv_t8code/containers.jl | 1 - src/solvers/fv_t8code/fv.jl | 117 +++++++++++++++--- 7 files changed, 209 insertions(+), 31 deletions(-) create mode 100644 examples/t8code_3d_fv/elixir_advection_basic.jl create mode 100644 examples/t8code_3d_fv/elixir_advection_gauss.jl diff --git a/examples/t8code_3d_fv/elixir_advection_basic.jl b/examples/t8code_3d_fv/elixir_advection_basic.jl new file mode 100644 index 00000000000..071ef23dc8d --- /dev/null +++ b/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() diff --git a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl index 67c664471d5..a9528182320 100644 --- a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl +++ b/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) diff --git a/examples/t8code_3d_fv/elixir_advection_gauss.jl b/examples/t8code_3d_fv/elixir_advection_gauss.jl new file mode 100644 index 00000000000..5c22d4fd6f9 --- /dev/null +++ b/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() diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl index 088f934cc3e..006179c1fc9 100644 --- a/src/equations/linear_scalar_advection_3d.jl +++ b/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) diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl index 0dafdf0af9e..65ad8377fe3 100644 --- a/src/meshes/t8code_mesh.jl +++ b/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 | + + # | |/ # +---+ # diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl index 44dcd19a5e5..635a391b4cb 100644 --- a/src/solvers/fv_t8code/containers.jl +++ b/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] diff --git a/src/solvers/fv_t8code/fv.jl b/src/solvers/fv_t8code/fv.jl index 6445a2d7b23..1e4646c89a3 100644 --- a/src/solvers/fv_t8code/fv.jl +++ b/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