diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl new file mode 100644 index 0000000000..8b06e06974 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl @@ -0,0 +1,277 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal compressible Navier-Stokes equations + +prandtl_number() = 0.72 +mu() = 0.01 + +equations = CompressibleEulerEquations3D(1.4) +equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu=mu(), Prandtl=prandtl_number(), + gradient_variables=GradientVariablesPrimitive()) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs, + volume_integral=VolumeIntegralWeakForm()) + +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 = (2, 2, 2) + +mesh = P4estMesh(trees_per_dimension, polydeg=3, + coordinates_min=coordinates_min, coordinates_max=coordinates_max, + periodicity=(true, false, true), initial_refinement_level=2) + +# Note: the initial condition cannot be specialized to `CompressibleNavierStokesDiffusion3D` +# since it is called by both the parabolic solver (which passes in `CompressibleNavierStokesDiffusion3D`) +# and by the initial condition (which passes in `CompressibleEulerEquations3D`). +# This convergence test setup was originally derived by Andrew Winters (@andrewwinters5000) +function initial_condition_navier_stokes_convergence_test(x, t, equations) + # Constants. OBS! Must match those in `source_terms_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + v1 = A2 * sin(pi_x) * log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) * sin(pi_z) * cos(pi_t) + v2 = v1 + v3 = v1 + p = rho^2 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +@inline function source_terms_navier_stokes_convergence_test(u, x, t, equations) + # TODO: parabolic + # we currently need to hardcode these parameters until we fix the "combined equation" issue + # see also https://github.com/trixi-framework/Trixi.jl/pull/1160 + inv_gamma_minus_one = inv(equations.gamma - 1) + Pr = prandtl_number() + mu_ = mu() + + # Constants. OBS! Must match those in `initial_condition_navier_stokes_convergence_test` + c = 2.0 + A1 = 0.5 + A2 = 1.0 + A3 = 0.5 + + # Convenience values for trig. functions + pi_x = pi * x[1] + pi_y = pi * x[2] + pi_z = pi * x[3] + pi_t = pi * t + + # Define auxiliary functions for the strange function of the y variable + # to make expressions easier to read + g = log(x[2] + 2.0) * (1.0 - exp(-A3 * (x[2] - 1.0))) + g_y = ( A3 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) + + (1.0 - exp(-A3 * (x[2] - 1.0))) / (x[2] + 2.0) ) + g_yy = ( 2.0 * A3 * exp(-A3 * (x[2] - 1.0)) / (x[2] + 2.0) + - (1.0 - exp(-A3 * (x[2] - 1.0))) / ((x[2] + 2.0)^2) + - A3^2 * log(x[2] + 2.0) * exp(-A3 * (x[2] - 1.0)) ) + + # Density and its derivatives + rho = c + A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_t = -pi * A1 * sin(pi_x) * cos(pi_y) * sin(pi_z) * sin(pi_t) + rho_x = pi * A1 * cos(pi_x) * cos(pi_y) * sin(pi_z) * cos(pi_t) + rho_y = -pi * A1 * sin(pi_x) * sin(pi_y) * sin(pi_z) * cos(pi_t) + rho_z = pi * A1 * sin(pi_x) * cos(pi_y) * cos(pi_z) * cos(pi_t) + rho_xx = -pi^2 * (rho - c) + rho_yy = -pi^2 * (rho - c) + rho_zz = -pi^2 * (rho - c) + + # Velocities and their derivatives + # v1 terms + v1 = A2 * sin(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_t = -pi * A2 * sin(pi_x) * g * sin(pi_z) * sin(pi_t) + v1_x = pi * A2 * cos(pi_x) * g * sin(pi_z) * cos(pi_t) + v1_y = A2 * sin(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_z = pi * A2 * sin(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_xx = -pi^2 * v1 + v1_yy = A2 * sin(pi_x) * g_yy * sin(pi_z) * cos(pi_t) + v1_zz = -pi^2 * v1 + v1_xy = pi * A2 * cos(pi_x) * g_y * sin(pi_z) * cos(pi_t) + v1_xz = pi^2 * A2 * cos(pi_x) * g * cos(pi_z) * cos(pi_t) + v1_yz = pi * A2 * sin(pi_x) * g_y * cos(pi_z) * cos(pi_t) + # v2 terms (simplifies from ansatz) + v2 = v1 + v2_t = v1_t + v2_x = v1_x + v2_y = v1_y + v2_z = v1_z + v2_xx = v1_xx + v2_yy = v1_yy + v2_zz = v1_zz + v2_xy = v1_xy + v2_yz = v1_yz + # v3 terms (simplifies from ansatz) + v3 = v1 + v3_t = v1_t + v3_x = v1_x + v3_y = v1_y + v3_z = v1_z + v3_xx = v1_xx + v3_yy = v1_yy + v3_zz = v1_zz + v3_xz = v1_xz + v3_yz = v1_yz + + # Pressure and its derivatives + p = rho^2 + p_t = 2.0 * rho * rho_t + p_x = 2.0 * rho * rho_x + p_y = 2.0 * rho * rho_y + p_z = 2.0 * rho * rho_z + + # Total energy and its derivatives; simiplifies from ansatz that v2 = v1 and v3 = v1 + E = p * inv_gamma_minus_one + 1.5 * rho * v1^2 + E_t = p_t * inv_gamma_minus_one + 1.5 * rho_t * v1^2 + 3.0 * rho * v1 * v1_t + E_x = p_x * inv_gamma_minus_one + 1.5 * rho_x * v1^2 + 3.0 * rho * v1 * v1_x + E_y = p_y * inv_gamma_minus_one + 1.5 * rho_y * v1^2 + 3.0 * rho * v1 * v1_y + E_z = p_z * inv_gamma_minus_one + 1.5 * rho_z * v1^2 + 3.0 * rho * v1 * v1_z + + # Divergence of Fick's law ∇⋅∇q = kappa ∇⋅∇T; simplifies because p = rho², so T = p/rho = rho + kappa = equations.gamma * inv_gamma_minus_one / Pr + q_xx = kappa * rho_xx # kappa T_xx + q_yy = kappa * rho_yy # kappa T_yy + q_zz = kappa * rho_zz # kappa T_zz + + # Stress tensor and its derivatives (exploit symmetry) + tau11 = 4.0 / 3.0 * v1_x - 2.0 / 3.0 * (v2_y + v3_z) + tau12 = v1_y + v2_x + tau13 = v1_z + v3_x + tau22 = 4.0 / 3.0 * v2_y - 2.0 / 3.0 * (v1_x + v3_z) + tau23 = v2_z + v3_y + tau33 = 4.0 / 3.0 * v3_z - 2.0 / 3.0 * (v1_x + v2_y) + + tau11_x = 4.0 / 3.0 * v1_xx - 2.0 / 3.0 * (v2_xy + v3_xz) + tau12_x = v1_xy + v2_xx + tau13_x = v1_xz + v3_xx + + tau12_y = v1_yy + v2_xy + tau22_y = 4.0 / 3.0 * v2_yy - 2.0 / 3.0 * (v1_xy + v3_yz) + tau23_y = v2_yz + v3_yy + + tau13_z = v1_zz + v3_xz + tau23_z = v2_zz + v3_yz + tau33_z = 4.0 / 3.0 * v3_zz - 2.0 / 3.0 * (v1_xz + v2_yz) + + # Compute the source terms + # Density equation + du1 = ( rho_t + rho_x * v1 + rho * v1_x + + rho_y * v2 + rho * v2_y + + rho_z * v3 + rho * v3_z ) + # x-momentum equation + du2 = ( rho_t * v1 + rho * v1_t + p_x + rho_x * v1^2 + + 2.0 * rho * v1 * v1_x + + rho_y * v1 * v2 + + rho * v1_y * v2 + + rho * v1 * v2_y + + rho_z * v1 * v3 + + rho * v1_z * v3 + + rho * v1 * v3_z + - mu_ * (tau11_x + tau12_y + tau13_z) ) + # y-momentum equation + du3 = ( rho_t * v2 + rho * v2_t + p_y + rho_x * v1 * v2 + + rho * v1_x * v2 + + rho * v1 * v2_x + + rho_y * v2^2 + + 2.0 * rho * v2 * v2_y + + rho_z * v2 * v3 + + rho * v2_z * v3 + + rho * v2 * v3_z + - mu_ * (tau12_x + tau22_y + tau23_z) ) + # z-momentum equation + du4 = ( rho_t * v3 + rho * v3_t + p_z + rho_x * v1 * v3 + + rho * v1_x * v3 + + rho * v1 * v3_x + + rho_y * v2 * v3 + + rho * v2_y * v3 + + rho * v2 * v3_y + + rho_z * v3^2 + + 2.0 * rho * v3 * v3_z + - mu_ * (tau13_x + tau23_y + tau33_z) ) + # Total energy equation + du5 = ( E_t + v1_x * (E + p) + v1 * (E_x + p_x) + + v2_y * (E + p) + v2 * (E_y + p_y) + + v3_z * (E + p) + v3 * (E_z + p_z) + # stress tensor and temperature gradient from x-direction + - mu_ * ( q_xx + v1_x * tau11 + v2_x * tau12 + v3_x * tau13 + + v1 * tau11_x + v2 * tau12_x + v3 * tau13_x) + # stress tensor and temperature gradient terms from y-direction + - mu_ * ( q_yy + v1_y * tau12 + v2_y * tau22 + v3_y * tau23 + + v1 * tau12_y + v2 * tau22_y + v3 * tau23_y) + # stress tensor and temperature gradient terms from z-direction + - mu_ * ( q_zz + v1_z * tau13 + v2_z * tau23 + v3_z * tau33 + + v1 * tau13_z + v2 * tau23_z + v3 * tau33_z) ) + + return SVector(du1, du2, du3, du4, du5) +end + +initial_condition = initial_condition_navier_stokes_convergence_test + +# BC types +velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x, t, equations)[2:4]) +heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0) +boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom, heat_bc_top_bottom) + +# define inviscid boundary conditions +# boundary_conditions = Dict( +# :y_neg => boundary_condition_slip_wall, +# :y_pos => boundary_condition_slip_wall +# ) + +boundary_conditions = Dict( :x_neg => boundary_condition_periodic, + :x_pos => boundary_condition_periodic, + :y_neg => boundary_condition_slip_wall, + :y_pos => boundary_condition_slip_wall, + :z_neg => boundary_condition_periodic, + :z_pos => boundary_condition_periodic) + +# define viscous boundary conditions +# boundary_conditions_parabolic = Dict( +# :y_neg => boundary_condition_top_bottom, +# :y_pos => boundary_condition_top_bottom +# ) + +boundary_conditions_parabolic = Dict( :x_neg => boundary_condition_periodic, + :x_pos => boundary_condition_periodic, + :y_neg => boundary_condition_top_bottom, + :y_pos => boundary_condition_top_bottom, + :z_neg => boundary_condition_periodic, + :z_pos => boundary_condition_periodic) + +semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, solver; + boundary_conditions=(boundary_conditions, boundary_conditions_parabolic), + source_terms=source_terms_navier_stokes_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span `tspan` +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() +alive_callback = AliveCallback(alive_interval=10) +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval) +callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback) + +############################################################################### +# run the simulation + +time_int_tol = 1e-8 +sol = solve(ode, RDPK3SpFSAL49(); abstol=time_int_tol, reltol=time_int_tol, dt = 1e-5, + ode_default_options()..., callback=callbacks) +summary_callback() # print the timer summary + diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index 73ac47ed1e..7e90a83a9c 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -1,7 +1,7 @@ # This method is called when a SemidiscretizationHyperbolicParabolic is constructed. # It constructs the basic `cache` used throughout the simulation to compute # the RHS etc. -function create_cache_parabolic(mesh::P4estMesh, equations_hyperbolic::AbstractEquations, +function create_cache_parabolic(mesh::P4estMesh{2}, equations_hyperbolic::AbstractEquations, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, RealT, uEltype) balance!(mesh) diff --git a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl index b0a79d69de..3040c330f5 100644 --- a/src/solvers/dgsem_p4est/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_3d_parabolic.jl @@ -1,26 +1,26 @@ # This method is called when a SemidiscretizationHyperbolicParabolic is constructed. # It constructs the basic `cache` used throughout the simulation to compute # the RHS etc. -# function create_cache_parabolic(mesh::P4estMesh, equations_hyperbolic::AbstractEquations, -# equations_parabolic::AbstractEquationsParabolic, -# dg::DG, parabolic_scheme, RealT, uEltype) -# balance!(mesh) +function create_cache_parabolic(mesh::P4estMesh{3}, equations_hyperbolic::AbstractEquations, + equations_parabolic::AbstractEquationsParabolic, + dg::DG, parabolic_scheme, RealT, uEltype) + balance!(mesh) -# elements = init_elements(mesh, equations_hyperbolic, dg.basis, uEltype) -# interfaces = init_interfaces(mesh, equations_hyperbolic, dg.basis, elements) -# boundaries = init_boundaries(mesh, equations_hyperbolic, dg.basis, elements) + elements = init_elements(mesh, equations_hyperbolic, dg.basis, uEltype) + interfaces = init_interfaces(mesh, equations_hyperbolic, dg.basis, elements) + boundaries = init_boundaries(mesh, equations_hyperbolic, dg.basis, elements) -# n_vars = nvariables(equations_hyperbolic) -# n_elements = nelements(elements) -# n_nodes = nnodes(dg.basis) # nodes in one direction -# u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_elements) -# gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) -# flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) + n_vars = nvariables(equations_hyperbolic) + n_elements = nelements(elements) + n_nodes = nnodes(dg.basis) # nodes in one direction + u_transformed = Array{uEltype}(undef, n_vars, n_nodes, n_nodes, n_nodes, n_elements) + gradients = ntuple(_ -> similar(u_transformed), ndims(mesh)) + flux_viscous = ntuple(_ -> similar(u_transformed), ndims(mesh)) -# cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) + cache = (; elements, interfaces, boundaries, gradients, flux_viscous, u_transformed) -# return cache -# end + return cache +end function calc_gradient!(gradients, u_transformed, t, mesh::P4estMesh{3}, equations_parabolic, @@ -64,34 +64,45 @@ function calc_gradient!(gradients, u_transformed, t, u_node, equations_parabolic, dg, i, j, kk, element) end - + end # now that the reference coordinate gradients are computed, transform them node-by-node to physical gradients # using the contravariant vectors - for j in eachnode(dg), i in eachnode(dg) - Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j,k, element) - Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, + Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j,k, + element) + Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j,k, element) gradients_reference_1 = get_node_vars(gradients_x, equations_parabolic, dg, - i, j, element) + i, j, k, element) gradients_reference_2 = get_node_vars(gradients_y, equations_parabolic, dg, - i, j, element) + i, j, k, element) + gradients_reference_3 = get_node_vars(gradients_z, equations_parabolic, dg, + i, j, k, element) # note that the contravariant vectors are transposed compared with computations of flux # divergences in `calc_volume_integral!`. See # https://github.com/trixi-framework/Trixi.jl/pull/1490#discussion_r1213345190 # for a more detailed discussion. gradient_x_node = Ja11 * gradients_reference_1 + - Ja21 * gradients_reference_2 + Ja21 * gradients_reference_2 + + Ja31 * gradients_reference_3 gradient_y_node = Ja12 * gradients_reference_1 + - Ja22 * gradients_reference_2 + Ja22 * gradients_reference_2 + + Ja32 * gradients_reference_3 + gradient_z_node = Ja13 * gradients_reference_1 + + Ja23 * gradients_reference_2 + + Ja33 * gradients_reference_3 - set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, i, j, + set_node_vars!(gradients_x, gradient_x_node, equations_parabolic, dg, i, j, k, + element) + set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, i, j, k, element) - set_node_vars!(gradients_y, gradient_y_node, equations_parabolic, dg, i, j, + set_node_vars!(gradients_z, gradient_z_node, equations_parabolic, dg, i, j, k, element) end end @@ -138,76 +149,57 @@ function calc_gradient!(gradients, u_transformed, t, factor_1 = boundary_interpolation[1, 1] factor_2 = boundary_interpolation[nnodes(dg), 2] @threaded for element in eachelement(dg, cache) - for l in eachnode(dg) + for m in eachnode(dg), l in eachnode(dg) for v in eachvariable(equations_parabolic) - - # Compute x-component of gradients - - # surface at -x - normal_direction_x, _ = get_normal_direction(1, contravariant_vectors, - 1, l, element) - gradients_x[v, 1, l, element] = (gradients_x[v, 1, l, element] + - surface_flux_values[v, l, 1, element] * - factor_1 * normal_direction_x) - - # surface at +x - normal_direction_x, _ = get_normal_direction(2, contravariant_vectors, - nnodes(dg), l, element) - gradients_x[v, nnodes(dg), l, element] = (gradients_x[v, nnodes(dg), l, - element] + - surface_flux_values[v, l, 2, - element] * - factor_2 * normal_direction_x) - - # surface at -y - normal_direction_x, _ = get_normal_direction(3, contravariant_vectors, - l, 1, element) - gradients_x[v, l, 1, element] = (gradients_x[v, l, 1, element] + - surface_flux_values[v, l, 3, element] * - factor_1 * normal_direction_x) - - # surface at +y - normal_direction_x, _ = get_normal_direction(4, contravariant_vectors, - l, nnodes(dg), element) - gradients_x[v, l, nnodes(dg), element] = (gradients_x[v, l, nnodes(dg), - element] + - surface_flux_values[v, l, 4, - element] * - factor_2 * normal_direction_x) - - # Compute y-component of gradients - - # surface at -x - _, normal_direction_y = get_normal_direction(1, contravariant_vectors, - 1, l, element) - gradients_y[v, 1, l, element] = (gradients_y[v, 1, l, element] + - surface_flux_values[v, l, 1, element] * - factor_1 * normal_direction_y) - - # surface at +x - _, normal_direction_y = get_normal_direction(2, contravariant_vectors, - nnodes(dg), l, element) - gradients_y[v, nnodes(dg), l, element] = (gradients_y[v, nnodes(dg), l, - element] + - surface_flux_values[v, l, 2, - element] * - factor_2 * normal_direction_y) - - # surface at -y - _, normal_direction_y = get_normal_direction(3, contravariant_vectors, - l, 1, element) - gradients_y[v, l, 1, element] = (gradients_y[v, l, 1, element] + - surface_flux_values[v, l, 3, element] * - factor_1 * normal_direction_y) - - # surface at +y - _, normal_direction_y = get_normal_direction(4, contravariant_vectors, - l, nnodes(dg), element) - gradients_y[v, l, nnodes(dg), element] = (gradients_y[v, l, nnodes(dg), - element] + - surface_flux_values[v, l, 4, - element] * - factor_2 * normal_direction_y) + for dim in 1:3 + # surface at -x + normal_direction = get_normal_direction(1, contravariant_vectors, + 1, l, m, element) + gradients[dim][v, 1, l, m, element] = (gradients[dim][v, 1, l, m, element] + + surface_flux_values[v, l, m, 1, element] * + factor_1 * normal_direction[dim]) + + # surface at +x + normal_direction = get_normal_direction(2, contravariant_vectors, + nnodes(dg), l, m, element) + gradients[dim][v, nnodes(dg), l, m, element] = (gradients[dim][v, nnodes(dg), l, m, + element] + + surface_flux_values[v, l, m, 2, + element] * + factor_2 * normal_direction[dim]) + + # surface at -y + normal_direction = get_normal_direction(3, contravariant_vectors, + l, m, 1, element) + gradients[dim][v, l, m, 1, element] = (gradients[dim][v, l, m, 1, element] + + surface_flux_values[v, l, m, 3, element] * + factor_1 * normal_direction[dim]) + + # surface at +y + normal_direction = get_normal_direction(4, contravariant_vectors, + l, nnodes(dg), m, element) + gradients[dim][v, l, nnodes(dg), m, element] = (gradients[dim][v, l, nnodes(dg), m, + element] + + surface_flux_values[v, l, m, 4, + element] * + factor_2 * normal_direction[dim]) + + # surface at -z + normal_direction = get_normal_direction(5, contravariant_vectors, + l, m, 1, element) + gradients[dim][v, l, m, 1, element] = (gradients[dim][v, l, m, 1, element] + + surface_flux_values[v, l, m, 5, element] * + factor_1 * normal_direction[dim]) + + # surface at +z + normal_direction = get_normal_direction(6, contravariant_vectors, + l, m, nnodes(dg), element) + gradients[dim][v, l, m, nnodes(dg), element] = (gradients[dim][v, l, m, nnodes(dg), + element] + + surface_flux_values[v, l, m, 6, + element] * + factor_2 * normal_direction[dim]) + end end end end @@ -219,25 +211,28 @@ function calc_gradient!(gradients, u_transformed, t, cache_parabolic) apply_jacobian_parabolic!(gradients_y, mesh, equations_parabolic, dg, cache_parabolic) + apply_jacobian_parabolic!(gradients_z, mesh, equations_parabolic, dg, + cache_parabolic) end return nothing end # This version is used for parabolic gradient computations -@inline function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{2}, +@inline function calc_interface_flux!(surface_flux_values, mesh::P4estMesh{3}, nonconservative_terms::False, equations::AbstractEquationsParabolic, surface_integral, dg::DG, cache, interface_index, normal_direction, - primary_node_index, primary_direction_index, - primary_element_index, - secondary_node_index, secondary_direction_index, + primary_i_node_index, primary_j_node_index, + primary_direction_index, primary_element_index, + secondary_i_node_index, secondary_j_node_index, + secondary_direction_index, secondary_element_index) @unpack u = cache.interfaces @unpack surface_flux = surface_integral - u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_node_index, + u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_i_node_index, primary_j_node_index, interface_index) flux_ = 0.5 * (u_ll + u_rr) # we assume that the gradient computations utilize a central flux @@ -245,42 +240,52 @@ end # Note that we don't flip the sign on the secondondary flux. This is because for parabolic terms, # the normals are not embedded in `flux_` for the parabolic gradient computations. for v in eachvariable(equations) - surface_flux_values[v, primary_node_index, primary_direction_index, primary_element_index] = flux_[v] - surface_flux_values[v, secondary_node_index, secondary_direction_index, secondary_element_index] = flux_[v] + surface_flux_values[v, primary_i_node_index, primary_j_node_index, primary_direction_index, primary_element_index] = flux_[v] + surface_flux_values[v, secondary_i_node_index, secondary_j_node_index, secondary_direction_index, secondary_element_index] = flux_[v] end end # This is the version used when calculating the divergence of the viscous fluxes function calc_volume_integral!(du, flux_viscous, - mesh::P4estMesh{2}, + mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, dg::DGSEM, cache) (; derivative_dhat) = dg.basis (; contravariant_vectors) = cache.elements - flux_viscous_x, flux_viscous_y = flux_viscous + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @threaded for element in eachelement(dg, cache) # Calculate volume terms in one element - for j in eachnode(dg), i in eachnode(dg) - flux1 = get_node_vars(flux_viscous_x, equations_parabolic, dg, i, j, element) - flux2 = get_node_vars(flux_viscous_y, equations_parabolic, dg, i, j, element) + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + flux1 = get_node_vars(flux_viscous_x, equations_parabolic, dg, i, j, k, element) + flux2 = get_node_vars(flux_viscous_y, equations_parabolic, dg, i, j, k, element) + flux3 = get_node_vars(flux_viscous_z, equations_parabolic, dg, i, j, k, element) # Compute the contravariant flux by taking the scalar product of the # first contravariant vector Ja^1 and the flux vector - Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j, element) - contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j, k, element) + contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 for ii in eachnode(dg) multiply_add_to_node_vars!(du, derivative_dhat[ii, i], contravariant_flux1, - equations_parabolic, dg, ii, j, element) + equations_parabolic, dg, ii, j, k, element) end # Compute the contravariant flux by taking the scalar product of the # second contravariant vector Ja^2 and the flux vector - Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element) - contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j, k, element) + contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 for jj in eachnode(dg) multiply_add_to_node_vars!(du, derivative_dhat[jj, j], contravariant_flux2, - equations_parabolic, dg, i, jj, element) + equations_parabolic, dg, i, jj, k, element) + end + + # Compute the contravariant flux by taking the scalar product of the + # second contravariant vector Ja^2 and the flux vector + Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j, k, element) + contravariant_flux3 = Ja31 * flux1 + Ja32 * flux2 + Ja33 * flux3 + for kk in eachnode(dg) + multiply_add_to_node_vars!(du, derivative_dhat[kk, k], contravariant_flux3, + equations_parabolic, dg, i, j, kk, element) end end end @@ -291,13 +296,13 @@ end # This is the version used when calculating the divergence of the viscous fluxes # We pass the `surface_integral` argument solely for dispatch function prolong2interfaces!(cache_parabolic, flux_viscous, - mesh::P4estMesh{2}, + mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) (; interfaces) = cache_parabolic (; contravariant_vectors) = cache_parabolic.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y = flux_viscous + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @threaded for interface in eachinterface(dg, cache) # Copy solution data from the primary element using "delayed indexing" with @@ -309,31 +314,42 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, primary_indices = interfaces.node_indices[1, interface] primary_direction = indices2direction(primary_indices) - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) + i_primary_start, i_primary_step_i, i_primary_step_j = index_to_start_step_3d(primary_indices[1], + index_range) + j_primary_start, j_primary_step_i, j_primary_step_j = index_to_start_step_3d(primary_indices[2], + index_range) + k_primary_start, k_primary_step_i, k_primary_step_j = index_to_start_step_3d(primary_indices[3], + index_range) i_primary = i_primary_start j_primary = j_primary_start - for i in eachnode(dg) + k_primary = k_primary_start - # this is the outward normal direction on the primary element - normal_direction = get_normal_direction(primary_direction, + for j in eachnode(dg) + for i in eachnode(dg) + # this is the outward normal direction on the primary element + normal_direction = get_normal_direction(primary_direction, contravariant_vectors, - i_primary, j_primary, primary_element) - - for v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, - primary_element], - flux_viscous_y[v, i_primary, j_primary, - primary_element]) + i_primary, j_primary, k_primary, primary_element) - interfaces.u[1, v, i, interface] = dot(flux_viscous, normal_direction) + for v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + flux_viscous = SVector(flux_viscous_x[v, i_primary, j_primary, k_primary, + primary_element], + flux_viscous_y[v, i_primary, j_primary, k_primary, + primary_element], + flux_viscous_z[v, i_primary, j_primary, k_primary, + primary_element]) + + interfaces.u[1, v, i, j, interface] = dot(flux_viscous, normal_direction) + end + i_primary += i_primary_step_i + j_primary += j_primary_step_i + k_primary += k_primary_step_i end - i_primary += i_primary_step - j_primary += j_primary_step + i_primary += i_primary_step_j + j_primary += j_primary_step_j + k_primary += k_primary_step_j end # Copy solution data from the secondary element using "delayed indexing" with @@ -342,33 +358,44 @@ function prolong2interfaces!(cache_parabolic, flux_viscous, secondary_indices = interfaces.node_indices[2, interface] secondary_direction = indices2direction(secondary_indices) - i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], - index_range) - j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], - index_range) + i_secondary_start, i_secondary_step_i, i_secondary_step_j = index_to_start_step_3d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step_i, j_secondary_step_j = index_to_start_step_3d(secondary_indices[2], + index_range) + k_secondary_start, k_secondary_step_i, k_secondary_step_j = index_to_start_step_3d(secondary_indices[3], + index_range) i_secondary = i_secondary_start j_secondary = j_secondary_start - for i in eachnode(dg) - # This is the outward normal direction on the secondary element. - # Here, we assume that normal_direction on the secondary element is - # the negative of normal_direction on the primary element. - normal_direction = get_normal_direction(secondary_direction, - contravariant_vectors, - i_secondary, j_secondary, - secondary_element) - - for v in eachvariable(equations_parabolic) - # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! - flux_viscous = SVector(flux_viscous_x[v, i_secondary, j_secondary, - secondary_element], - flux_viscous_y[v, i_secondary, j_secondary, - secondary_element]) - # store the normal flux with respect to the primary normal direction - interfaces.u[2, v, i, interface] = -dot(flux_viscous, normal_direction) + k_secondary = k_secondary_start + for j in eachnode(dg) + for i in eachnode(dg) + # This is the outward normal direction on the secondary element. + # Here, we assume that normal_direction on the secondary element is + # the negative of normal_direction on the primary element. + normal_direction = get_normal_direction(secondary_direction, + contravariant_vectors, + i_secondary, j_secondary, k_secondary, + secondary_element) + + for v in eachvariable(equations_parabolic) + # OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*! + flux_viscous = SVector(flux_viscous_x[v, i_secondary, j_secondary, k_secondary, + secondary_element], + flux_viscous_y[v, i_secondary, j_secondary, k_secondary, + secondary_element], + flux_viscous_z[v, i_secondary, j_secondary, k_secondary, + secondary_element]) + # store the normal flux with respect to the primary normal direction + interfaces.u[2, v, i, j, interface] = -dot(flux_viscous, normal_direction) + end + i_secondary += i_secondary_step_i + j_secondary += j_secondary_step_i + k_secondary += k_secondary_step_i end - i_secondary += i_secondary_step - j_secondary += j_secondary_step + i_secondary += i_secondary_step_j + j_secondary += j_secondary_step_j + k_secondary += k_secondary_step_j end end @@ -377,12 +404,10 @@ end # This version is used for divergence flux computations function calc_interface_flux!(surface_flux_values, - mesh::P4estMesh{2}, equations_parabolic, + mesh::P4estMesh{3}, equations_parabolic, dg::DG, cache_parabolic) (; neighbor_ids, node_indices) = cache_parabolic.interfaces - (; contravariant_vectors) = cache_parabolic.elements index_range = eachnode(dg) - index_end = last(index_range) @threaded for interface in eachinterface(dg, cache_parabolic) # Get element and side index information on the primary element @@ -390,51 +415,67 @@ function calc_interface_flux!(surface_flux_values, primary_indices = node_indices[1, interface] primary_direction_index = indices2direction(primary_indices) - # Create the local i,j indexing on the primary element used to pull normal direction information - i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], - index_range) - j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], - index_range) + i_primary_start, i_primary_step_i, i_primary_step_j = index_to_start_step_3d(primary_indices[1], + index_range) + j_primary_start, j_primary_step_i, j_primary_step_j = index_to_start_step_3d(primary_indices[2], + index_range) + k_primary_start, k_primary_step_i, k_primary_step_j = index_to_start_step_3d(primary_indices[3], + index_range) i_primary = i_primary_start j_primary = j_primary_start + k_primary = k_primary_start # Get element and side index information on the secondary element secondary_element = neighbor_ids[2, interface] secondary_indices = node_indices[2, interface] secondary_direction_index = indices2direction(secondary_indices) + secondary_surface_indices = surface_indices(secondary_indices) # Initiate the secondary index to be used in the surface for loop. # This index on the primary side will always run forward but # the secondary index might need to run backwards for flipped sides. - if :i_backward in secondary_indices - node_secondary = index_end - node_secondary_step = -1 - else - node_secondary = 1 - node_secondary_step = 1 - end + # Get the surface indexing on the secondary element. + # Note that the indices of the primary side will always run forward but + # the secondary indices might need to run backwards for flipped sides. + i_secondary_start, i_secondary_step_i, i_secondary_step_j = index_to_start_step_3d(secondary_surface_indices[1], + index_range) + j_secondary_start, j_secondary_step_i, j_secondary_step_j = index_to_start_step_3d(secondary_surface_indices[2], + index_range) + i_secondary = i_secondary_start + j_secondary = j_secondary_start - for node in eachnode(dg) - # We prolong the viscous flux dotted with respect the outward normal on the - # primary element. We assume a BR-1 type of flux. - viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache_parabolic.interfaces.u, - equations_parabolic, - dg, node, - interface) + for j in eachnode(dg) + for i in eachnode(dg) + # We prolong the viscous flux dotted with respect the outward normal on the + # primary element. We assume a BR-1 type of flux. + viscous_flux_normal_ll, viscous_flux_normal_rr = get_surface_node_vars(cache_parabolic.interfaces.u, + equations_parabolic, + dg, i, j, + interface) - flux = 0.5 * (viscous_flux_normal_ll + viscous_flux_normal_rr) + flux = 0.5 * (viscous_flux_normal_ll + viscous_flux_normal_rr) - for v in eachvariable(equations_parabolic) - surface_flux_values[v, node, primary_direction_index, primary_element] = flux[v] - surface_flux_values[v, node_secondary, secondary_direction_index, secondary_element] = -flux[v] - end + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, primary_direction_index, primary_element] = flux[v] + surface_flux_values[v, i_secondary, j_secondary, secondary_direction_index, secondary_element] = -flux[v] + end - # Increment primary element indices to pull the normal direction - i_primary += i_primary_step - j_primary += j_primary_step - # Increment the surface node index along the secondary element - node_secondary += node_secondary_step + # Increment the primary element indices + i_primary += i_primary_step_i + j_primary += j_primary_step_i + k_primary += k_primary_step_i + # Increment the secondary element surface indices + i_secondary += i_secondary_step_i + j_secondary += j_secondary_step_i + end + # Increment the primary element indices + i_primary += i_primary_step_j + j_primary += j_primary_step_j + k_primary += k_primary_step_j + # Increment the secondary element surface indices + i_secondary += i_secondary_step_j + j_secondary += j_secondary_step_j end end @@ -443,14 +484,14 @@ end # TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` function prolong2boundaries!(cache_parabolic, flux_viscous, - mesh::P4estMesh{2}, + mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG, cache) (; boundaries) = cache_parabolic (; contravariant_vectors) = cache_parabolic.elements index_range = eachnode(dg) - flux_viscous_x, flux_viscous_y = flux_viscous + flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous @threaded for boundary in eachboundary(dg, cache_parabolic) # Copy solution data from the element using "delayed indexing" with @@ -459,24 +500,37 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, node_indices = boundaries.node_indices[boundary] direction = indices2direction(node_indices) - i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) - j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + i_node_start, i_node_step_i, i_node_step_j = index_to_start_step_3d(node_indices[1], + index_range) + j_node_start, j_node_step_i, j_node_step_j = index_to_start_step_3d(node_indices[2], + index_range) + k_node_start, k_node_step_i, k_node_step_j = index_to_start_step_3d(node_indices[3], + index_range) i_node = i_node_start j_node = j_node_start - for i in eachnode(dg) - # this is the outward normal direction on the primary element - normal_direction = get_normal_direction(direction, contravariant_vectors, - i_node, j_node, element) + k_node = k_node_start + + for j in eachnode(dg) + for i in eachnode(dg) + # this is the outward normal direction on the primary element + normal_direction = get_normal_direction(direction, contravariant_vectors, + i_node, j_node, k_node, element) - for v in eachvariable(equations_parabolic) - flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, element], - flux_viscous_y[v, i_node, j_node, element]) + for v in eachvariable(equations_parabolic) + flux_viscous = SVector(flux_viscous_x[v, i_node, j_node, k_node, element], + flux_viscous_y[v, i_node, j_node, k_node, element], + flux_viscous_z[v, i_node, j_node, k_node, element]) - boundaries.u[v, i, boundary] = dot(flux_viscous, normal_direction) + boundaries.u[v, i, j, boundary] = dot(flux_viscous, normal_direction) + end + i_node += i_node_step_i + j_node += j_node_step_i + k_node += k_node_step_i end - i_node += i_node_step - j_node += j_node_step + i_node += i_node_step_j + j_node += j_node_step_j + k_node += k_node_step_j end end return nothing @@ -546,7 +600,7 @@ end function calc_boundary_flux!(cache, t, boundary_condition_parabolic, # works with Dict types boundary_condition_indices, - operator_type, mesh::P4estMesh{2}, + operator_type, mesh::P4estMesh{3}, equations_parabolic::AbstractEquationsParabolic, surface_integral, dg::DG) (; boundaries) = cache @@ -564,39 +618,51 @@ function calc_boundary_flux!(cache, t, node_indices = boundaries.node_indices[boundary_index] direction_index = indices2direction(node_indices) - i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range) - j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range) + i_node_start, i_node_step_i, i_node_step_j = index_to_start_step_3d(node_indices[1], + index_range) + j_node_start, j_node_step_i, j_node_step_j = index_to_start_step_3d(node_indices[2], + index_range) + k_node_start, k_node_step_i, k_node_step_j = index_to_start_step_3d(node_indices[3], + index_range) i_node = i_node_start j_node = j_node_start - for node_index in eachnode(dg) - # Extract solution data from boundary container - u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, node_index, - boundary_index) - - # Outward-pointing normal direction (not normalized) - normal_direction = get_normal_direction(direction_index, contravariant_vectors, - i_node, j_node, element) - - # TODO: revisit if we want more general boundary treatments. - # This assumes the gradient numerical flux at the boundary is the gradient variable, - # which is consistent with BR1, LDG. - flux_inner = u_inner - - # Coordinates at boundary node - x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, - element) - - flux_ = boundary_condition_parabolic(flux_inner, u_inner, normal_direction, - x, t, operator_type, equations_parabolic) - - # Copy flux to element storage in the correct orientation - for v in eachvariable(equations_parabolic) - surface_flux_values[v, node_index, direction_index, element] = flux_[v] - end + k_node = k_node_start - i_node += i_node_step - j_node += j_node_step + for j in eachnode(dg) + for i in eachnode(dg) + # Extract solution data from boundary container + u_inner = get_node_vars(boundaries.u, equations_parabolic, dg, i, j, + boundary_index) + + # Outward-pointing normal direction (not normalized) + normal_direction = get_normal_direction(direction_index, contravariant_vectors, + i_node, j_node, k_node, element) + + # TODO: revisit if we want more general boundary treatments. + # This assumes the gradient numerical flux at the boundary is the gradient variable, + # which is consistent with BR1, LDG. + flux_inner = u_inner + + # Coordinates at boundary node + x = get_node_coords(node_coordinates, equations_parabolic, dg, i_node, j_node, k_node, + element) + + flux_ = boundary_condition_parabolic(flux_inner, u_inner, normal_direction, + x, t, operator_type, equations_parabolic) + + # Copy flux to element storage in the correct orientation + for v in eachvariable(equations_parabolic) + surface_flux_values[v, i, j, direction_index, element] = flux_[v] + end + + i_node += i_node_step_i + j_node += j_node_step_i + k_node += k_node_step_i + end + i_node += i_node_step_j + j_node += j_node_step_j + k_node += k_node_step_j end end end diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index d6d7463702..e089fd866a 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -13,7 +13,7 @@ # 2. compute f(u, grad(u)) # 3. compute div(f(u, grad(u))) (i.e., the "regular" rhs! call) # boundary conditions will be applied to both grad(u) and div(f(u, grad(u))). -function rhs_parabolic!(du, u, t, mesh::TreeMesh{3}, +function rhs_parabolic!(du, u, t, mesh::Union{P4estMesh{3},TreeMesh{3}}, equations_parabolic::AbstractEquationsParabolic, initial_condition, boundary_conditions_parabolic, source_terms, dg::DG, parabolic_scheme, cache, cache_parabolic) @@ -105,7 +105,7 @@ end # Transform solution variables prior to taking the gradient # (e.g., conservative to primitive variables). Defaults to doing nothing. # TODO: can we avoid copying data? -function transform_variables!(u_transformed, u, mesh::TreeMesh{3}, +function transform_variables!(u_transformed, u, mesh::Union{P4estMesh{3},TreeMesh{3}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, parabolic_scheme, cache, cache_parabolic) @threaded for element in eachelement(dg, cache) @@ -325,7 +325,7 @@ function prolong2boundaries!(cache_parabolic, flux_viscous, return nothing end -function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::TreeMesh{3}, +function calc_viscous_fluxes!(flux_viscous, gradients, u_transformed, mesh::Union{P4estMesh{3},TreeMesh{3}}, equations_parabolic::AbstractEquationsParabolic, dg::DG, cache, cache_parabolic) gradients_x, gradients_y, gradients_z = gradients @@ -806,7 +806,7 @@ end # This is because the parabolic fluxes are assumed to be of the form # `du/dt + df/dx = dg/dx + source(x,t)`, # where f(u) is the inviscid flux and g(u) is the viscous flux. -function apply_jacobian_parabolic!(du, mesh::TreeMesh{3}, +function apply_jacobian_parabolic!(du, mesh::Union{P4estMesh{3},TreeMesh{3}}, equations::AbstractEquationsParabolic, dg::DG, cache) @threaded for element in eachelement(dg, cache) factor = cache.elements.inverse_jacobian[element] diff --git a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl index cad5542aae..fc79abf3d1 100644 --- a/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl +++ b/src/solvers/dgsem_unstructured/sort_boundary_conditions.jl @@ -67,6 +67,7 @@ function initialize!(boundary_types_container::UnstructuredSortedBoundaryTypes{N else for key in keys(boundary_dictionary) if !(key in unique_names) + @show unique_names error("Key $(repr(key)) is not a valid boundary name") end end