diff --git a/examples/elixir_euler_spherical_advection_cartesian.jl b/examples/elixir_euler_spherical_advection_cartesian.jl deleted file mode 100644 index 9a0e7b5..0000000 --- a/examples/elixir_euler_spherical_advection_cartesian.jl +++ /dev/null @@ -1,153 +0,0 @@ - -using OrdinaryDiffEq -using Trixi -using TrixiAtmo - -############################################################################### -# semidiscretization of the linear advection equation - -# We use the Euler equations structure but modify the rhs! function to convert it to a -# variable-coefficient advection equation -equations = CompressibleEulerEquations3D(1.4) - -# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux -polydeg = 3 -solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs) - -# Initial condition for a Gaussian density profile with constant pressure -# and the velocity of a rotating solid body -function initial_condition_advection_sphere(x, t, equations::CompressibleEulerEquations3D) - # Gaussian density - rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) - # Constant pressure - p = 1.0 - - # Spherical coordinates for the point x - if sign(x[2]) == 0.0 - signy = 1.0 - else - signy = sign(x[2]) - end - # Co-latitude - colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) - # Latitude (auxiliary variable) - lat = -colat + 0.5 * pi - # Longitude - r_xy = sqrt(x[1]^2 + x[2]^2) - if r_xy == 0.0 - phi = pi / 2 - else - phi = signy * acos(x[1] / r_xy) - end - - # Compute the velocity of a rotating solid body - # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system) - v0 = 1.0 # Velocity at the "equator" - alpha = 0.0 #pi / 4 - v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) - v_lat = -v0 * sin(phi) * sin(alpha) - - # Transform to Cartesian coordinate system - v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long - v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long - v3 = sin(colat) * v_lat - - return prim2cons(SVector(rho, v1, v2, v3, p), equations) -end - -# Source term function to transform the Euler equations into the linear advection equations with variable advection velocity -function source_terms_convert_to_linear_advection(u, du, x, t, - equations::CompressibleEulerEquations3D) - v1 = u[2] / u[1] - v2 = u[3] / u[1] - v3 = u[4] / u[1] - - s2 = du[1] * v1 - du[2] - s3 = du[1] * v2 - du[3] - s4 = du[1] * v3 - du[4] - s5 = 0.5 * (s2 * v1 + s3 * v2 + s4 * v3) - du[5] - - return SVector(0.0, s2, s3, s4, s5) -end - -# Custom RHS that applies a source term that depends on du (used to convert the 3D Euler equations into the 3D linear advection -# equations with position-dependent advection speed) -function rhs_semi_custom!(du_ode, u_ode, semi, t) - # Compute standard Trixi RHS - Trixi.rhs!(du_ode, u_ode, semi, t) - - # Now apply the custom source term - Trixi.@trixi_timeit Trixi.timer() "custom source term" begin - @unpack solver, equations, cache = semi - @unpack node_coordinates = cache.elements - - # Wrap the solution and RHS - u = Trixi.wrap_array(u_ode, semi) - du = Trixi.wrap_array(du_ode, semi) - - Trixi.@threaded for element in eachelement(solver, cache) - for j in eachnode(solver), i in eachnode(solver) - u_local = Trixi.get_node_vars(u, equations, solver, i, j, element) - du_local = Trixi.get_node_vars(du, equations, solver, i, j, element) - x_local = Trixi.get_node_coords(node_coordinates, equations, solver, - i, j, element) - source = source_terms_convert_to_linear_advection(u_local, du_local, - x_local, t, equations) - Trixi.add_to_node_vars!(du, source, equations, solver, i, j, element) - end - end - end -end - -initial_condition = initial_condition_advection_sphere - -mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, - initial_refinement_level = 0) - -# A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span from 0.0 to π -tspan = (0.0, pi) -ode = semidiscretize(semi, tspan) - -# Create custom discretization that runs with the custom RHS -ode_semi_custom = ODEProblem(rhs_semi_custom!, - ode.u0, - ode.tspan, - semi) - -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers -summary_callback = SummaryCallback() - -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback = AnalysisCallback(semi, interval = 10, - save_analysis = true, - extra_analysis_errors = (:conservation_error,), - extra_analysis_integrals = (Trixi.density,)) - -# The SaveSolutionCallback allows to save the solution to a file in regular intervals -save_solution = SaveSolutionCallback(interval = 10, - solution_variables = cons2prim) - -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl = 0.7) - -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, - stepsize_callback) - -############################################################################### -# run the simulation - -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode_semi_custom, 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); - -# Print the timer summary -summary_callback() diff --git a/examples/elixir_spherical_advection_cartesian.jl b/examples/elixir_spherical_advection_cartesian.jl index 3d7755e..9a0e7b5 100644 --- a/examples/elixir_spherical_advection_cartesian.jl +++ b/examples/elixir_spherical_advection_cartesian.jl @@ -1,13 +1,13 @@ -############################################################################### -# DGSEM for the linear advection equation on the cubed using Lagrange multiplier approach -############################################################################### using OrdinaryDiffEq using Trixi +using TrixiAtmo ############################################################################### # semidiscretization of the linear advection equation +# We use the Euler equations structure but modify the rhs! function to convert it to a +# variable-coefficient advection equation equations = CompressibleEulerEquations3D(1.4) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux @@ -101,7 +101,8 @@ end initial_condition = initial_condition_advection_sphere -mesh = P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, initial_refinement_level = 0) +mesh = TrixiAtmo.P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, + initial_refinement_level = 0) # A semidiscretization collects data structures and functions for the spatial discretization semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) @@ -126,6 +127,7 @@ summary_callback = SummaryCallback() # The AnalysisCallback allows to analyse the solution in regular intervals and prints the results analysis_callback = AnalysisCallback(semi, interval = 10, save_analysis = true, + extra_analysis_errors = (:conservation_error,), extra_analysis_integrals = (Trixi.density,)) # The SaveSolutionCallback allows to save the solution to a file in regular intervals @@ -133,17 +135,19 @@ save_solution = SaveSolutionCallback(interval = 10, solution_variables = cons2prim) # The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -# stepsize_callback = StepsizeCallback(cfl = 1.4) +stepsize_callback = StepsizeCallback(cfl = 0.7) # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) +callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, + stepsize_callback) ############################################################################### # run the simulation # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks sol = solve(ode_semi_custom, CarpenterKennedy2N54(williamson_condition = false), - dt = pi * 1e-3, adaptive = false, save_everystep = false, callback = callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); # Print the timer summary summary_callback() diff --git a/examples/elixir_spherical_advection_covariant.jl b/examples/elixir_spherical_advection_covariant.jl index 63ba43a..b1b6f02 100644 --- a/examples/elixir_spherical_advection_covariant.jl +++ b/examples/elixir_spherical_advection_covariant.jl @@ -8,7 +8,7 @@ using OrdinaryDiffEq, Trixi, TrixiAtmo # Problem definition function initial_condition_advection_sphere(x, t, ::CovariantLinearAdvectionEquation2D) - + # Gaussian density rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) diff --git a/examples/elixir_spherical_advection_covariant_flux_differencing.jl b/examples/elixir_spherical_advection_covariant_flux_differencing.jl deleted file mode 100644 index d9e853b..0000000 --- a/examples/elixir_spherical_advection_covariant_flux_differencing.jl +++ /dev/null @@ -1,100 +0,0 @@ -############################################################################### -# DGSEM for the linear advection equation on the cubed sphere -############################################################################### - -using OrdinaryDiffEq, Trixi, TrixiAtmo - -############################################################################### -# Problem definition - -function initial_condition_advection_sphere(x, t, ::CovariantLinearAdvectionEquation2D) - # Gaussian density - - rho = 1.0 + exp(-20 * (x[1]^2 + x[3]^2)) - - # Spherical coordinates for the point x - if sign(x[2]) == 0.0 - signy = 1.0 - else - signy = sign(x[2]) - end - # Co-latitude - colat = acos(x[3] / sqrt(x[1]^2 + x[2]^2 + x[3]^2)) - # Latitude (auxiliary variable) - lat = -colat + 0.5 * pi - # Longitude - r_xy = sqrt(x[1]^2 + x[2]^2) - if r_xy == 0.0 - phi = pi / 2 - else - phi = signy * acos(x[1] / r_xy) - end - - # Compute the velocity of a rotating solid body - # (alpha is the angle between the rotation axis and the polar axis of the spherical coordinate system) - v0 = 1.0 # Velocity at the "equator" - alpha = 0.0 #pi / 4 - v_long = v0 * (cos(lat) * cos(alpha) + sin(lat) * cos(phi) * sin(alpha)) - v_lat = -v0 * sin(phi) * sin(alpha) - - # Transform to Cartesian coordinate system - v1 = -cos(colat) * cos(phi) * v_lat - sin(phi) * v_long - v2 = -cos(colat) * sin(phi) * v_lat + cos(phi) * v_long - v3 = sin(colat) * v_lat - - return SVector(rho, v1, v2, v3) -end - -############################################################################### -# Spatial discretization - -polydeg = 3 -cells_per_dimension = 2 - -mesh = P4estMeshCubedSphere2D(5, 1.0, polydeg = polydeg, initial_refinement_level = 0) - -equations = CovariantLinearAdvectionEquation2D() - -initial_condition = initial_condition_advection_sphere - -# Create DG solver with polynomial degree = p and a local Lax-Friedrichs flux -solver = DGSEM(polydeg = polydeg, surface_flux = flux_lax_friedrichs, - volume_integral = VolumeIntegralFluxDifferencing(flux_central)) - -# A semidiscretization collects data structures and functions for the spatial discretization -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solvers, callbacks etc. - -# Create ODE problem with time span from 0 to T -ode = semidiscretize(semi, (0.0, pi)) - -# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup -# and resets the timers -summary_callback = SummaryCallback() - -# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results -analysis_callback = AnalysisCallback(semi, interval = 10, - save_analysis = true, - extra_analysis_errors = (:conservation_error,)) - -# The SaveSolutionCallback allows to save the solution to a file in regular intervals -save_solution = SaveSolutionCallback(interval = 10, - solution_variables = cons2cons) - -# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step -stepsize_callback = StepsizeCallback(cfl = 0.4) - -# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver -callbacks = CallbackSet(summary_callback, analysis_callback, save_solution) - -############################################################################### -# run the simulation - -# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks -sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), - dt = pi * 1e-3, adaptive = false, save_everystep = false, callback = callbacks); - -# Print the timer summary -summary_callback() diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index 6d5e47a..9c2212e 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -17,6 +17,8 @@ using StrideArrays: StrideArray, StaticInt, PtrArray using LinearAlgebra: norm, dot using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect +@reexport using StaticArrays: SVector + include("auxiliary/auxiliary.jl") include("equations/equations.jl") include("meshes/meshes.jl") diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index ecefa87..36421ec 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -30,7 +30,8 @@ end """ Compute a given contravariant flux component """ -@inline function Trixi.flux(u, orientation::Integer, ::CovariantLinearAdvectionEquation2D) +@inline function Trixi.flux(u, orientation::Integer, + ::CovariantLinearAdvectionEquation2D) z = zero(eltype(u)) return SVector(u[orientation + 1] * u[1], z, z) end @@ -62,13 +63,12 @@ end return abs(u[2]), abs(u[3]) end - """ Convert from 4-vector of scalar plus 3 Cartesian components to 3-vector of scalar plus 2 contravariant velocity/momentum components """ -@inline function cartesian2contravariant(u_cartesian, ::CovariantLinearAdvectionEquation2D, +@inline function cartesian2contravariant(u_cartesian, + ::CovariantLinearAdvectionEquation2D, i, j, element, cache) - (; contravariant_vectors, inverse_jacobian) = cache.elements Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, j, @@ -89,7 +89,6 @@ end """ @inline function contravariant2cartesian(u_node, ::CovariantLinearAdvectionEquation2D, i, j, element, cache) - A11, A21, A31, A12, A22, A32 = view(cache.elements.jacobian_matrix, :, :, i, j, element) return SVector(u_node[1], @@ -97,5 +96,4 @@ end A21 * u_node[2] + A22 * u_node[3], A31 * u_node[2] + A32 * u_node[3]) end - end # @muladd diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index 687ae35..1ba9142 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -5,7 +5,8 @@ @muladd begin #! format: noindent -struct CovariantShallowWaterEquations2D{RealT <: Real} <: AbstractCovariantEquations2D{6} +struct CovariantShallowWaterEquations2D{RealT <: Real} <: + AbstractCovariantEquations2D{6} gravitational_acceleration::RealT rotation_rate::RealT end @@ -14,8 +15,7 @@ function Trixi.varnames(::typeof(cons2cons), ::CovariantShallowWaterEquations2D) return ("h", "M_con_1", "M_con_2", "G_cov_11", "G_cov_12", "G_cov_22") end -function Trixi.cons2entropy(u, equations::CovariantShallowWaterEquations2D) - +function Trixi.cons2entropy(u, equations::CovariantShallowWaterEquations2D) z = zero(eltype(u)) h, M_con_1, M_con_2, G_cov_11, G_cov_12, G_cov_22 = u (; gravitational_acceleration) = equations @@ -26,33 +26,32 @@ function Trixi.cons2entropy(u, equations::CovariantShallowWaterEquations2D) v_cov_1 = G_cov_11 * v_con_1 + G_cov_12 * v_con_2 v_cov_2 = G_cov_12 * v_con_1 + G_cov_22 * v_con_2 - w_1 = gravitational_acceleration * h - - 0.5 * (v_cov_1 * v_con_1 + v_cov_2 * v_con_2) + w_1 = gravitational_acceleration * h - + 0.5 * (v_cov_1 * v_con_1 + v_cov_2 * v_con_2) return SVector(w_1, v_cov_1, v_cov_2, z, z, z) - end """ Custom dissipation to ensure no flux is applied to metric coefficients """ @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, - normal_direction::AbstractVector, - equations::CovariantShallowWaterEquations2D) + normal_direction::AbstractVector, + equations::CovariantShallowWaterEquations2D) z = zero(eltype(u_ll)) λ = dissipation.max_abs_speed(u_ll, u_rr, normal_direction, equations) - return -0.5f0 * λ * SVector(u_rr[1] - u_ll[1], - u_rr[2] - u_ll[2], - u_rr[3] - u_ll[3], - z, z, z) + return -0.5f0 * λ * + SVector(u_rr[1] - u_ll[1], + u_rr[2] - u_ll[2], + u_rr[3] - u_ll[3], + z, z, z) end """ Compute a given contravariant flux component """ -@inline function Trixi.flux(u, orientation::Integer, +@inline function Trixi.flux(u, orientation::Integer, equations::CovariantShallowWaterEquations2D) - z = zero(eltype(u)) h, M_con_1, M_con_2, G_cov_11, G_cov_12, G_cov_22 = u @@ -64,14 +63,14 @@ end G_con_12 = -G_cov_12 / G G_con_22 = G_cov_11 / G - if orientation == 1 + if orientation == 1 M_con_j = M_con_1 - T_con_1j = (M_con_1 * M_con_1 / h) + G_con_11 * half_g_h_squared - T_con_2j = (M_con_2 * M_con_1 / h) + G_con_12 * half_g_h_squared + T_con_1j = (M_con_1 * M_con_1 / h) + G_con_11 * half_g_h_squared + T_con_2j = (M_con_2 * M_con_1 / h) + G_con_12 * half_g_h_squared else # orientation == 2 M_con_j = M_con_2 - T_con_1j = (M_con_2 * M_con_1 / h) + G_con_12 * half_g_h_squared - T_con_2j = (M_con_2 * M_con_2 / h) + G_con_22 * half_g_h_squared + T_con_1j = (M_con_2 * M_con_1 / h) + G_con_12 * half_g_h_squared + T_con_2j = (M_con_2 * M_con_2 / h) + G_con_22 * half_g_h_squared end return SVector(M_con_j, T_con_1j, T_con_2j, z, z, z) @@ -82,14 +81,12 @@ end """ @inline function Trixi.flux(u, normal_direction::AbstractVector, equations::CovariantShallowWaterEquations2D) - - return flux(u, 1, equations) * normal_direction[1] + + return flux(u, 1, equations) * normal_direction[1] + flux(u, 2, equations) * normal_direction[2] end @inline function source_terms_coriolis_sphere(u, x, t, equations::CovariantShallowWaterEquations2D) - _, M_con_1, M_con_2, G_cov_11, G_cov_12, G_cov_22 = u z = zero(eltype(u)) @@ -100,27 +97,24 @@ end G_con_22 = G_cov_11 / G fG = 2 * equations.rotation_rate * x[3] / norm(x) * G - + return SVector(z, fG * (G_con_12 * M_con_1 - G_con_11 * M_con_2), - fG * (G_con_22 * M_con_1 - G_con_12 * M_con_2), - z, z, z) + fG * (G_con_22 * M_con_1 - G_con_12 * M_con_2), + z, z, z) end - @inline function flux_nonconservative_naive(u_ll, u_rr, orientation::Integer, equations::CovariantShallowWaterEquations2D) - _, _, _, G_cov_11_rr, G_cov_12_rr, G_cov_22_rr = u_rr z = zero(eltype(u_ll)) - + _, T_con_11_ll, _ = flux(u_ll, orientation, equations) _, T_con_12_ll, T_con_22_ll = flux(u_ll, orientation, equations) G_ll = G_cov_11_ll * G_cov_22_ll - G_cov_12_ll^2 if orientation == 1 - a11_ll = T_con_11_ll * G_cov_22_ll / G_ll a12_ll = 2 * T_con_11_ll * G_cov_12_ll / G_ll a13_ll = (T_con_11_ll * G_cov_11 - 2 * T_con_12_ll * G_cov_12_ll - @@ -131,8 +125,8 @@ end a23_ll = (T_con_22_ll * G_cov_12 + 3 * T_con_12 * G_cov_11) / (2 * G_ll) else # orientation == 2 - - a11_ll = (T_con_11_ll * G_cov_12_ll + 3 * T_con_12_ll * G_cov_22_ll) / (2 * G_ll) + a11_ll = (T_con_11_ll * G_cov_12_ll + 3 * T_con_12_ll * G_cov_22_ll) / + (2 * G_ll) a12_ll = (T_con_22_ll * G_cov_22_ll - T_con_12_ll * G_cov_12_ll) / G_ll a13_ll = (T_con_12_ll * G_cov_11_ll - T_con_22_ll * G_cov_12_ll) / (2 * G_ll) @@ -141,18 +135,18 @@ end a22_ll = -2 * T_con_22_ll * G_cov_12_ll / G_ll a23_ll = T_con_22_ll * G_cov_11_ll / G_ll end - - return SVector(z, a11_ll * G_cov_11_rr + a12_ll * G_cov_12_rr + a13_ll * G_cov_22_rr, - a21_ll * G_cov_21_rr + a22_ll * G_cov_12_rr + a23_ll * G_cov_22_rr, - z, z, z) + + return SVector(z, + a11_ll * G_cov_11_rr + a12_ll * G_cov_12_rr + a13_ll * G_cov_22_rr, + a21_ll * G_cov_21_rr + a22_ll * G_cov_12_rr + a23_ll * G_cov_22_rr, + z, z, z) end """ Maximum directional wave speed for Lax-Friedrichs dissipation """ @inline function Trixi.max_abs_speed_naive(u_ll, u_rr, normal_direction, - ::CovariantShallowWaterEquations2D) - + ::CovariantShallowWaterEquations2D) h_ll, M_con_1_ll, M_con_2_ll, G_cov_11, G_cov_12, G_cov_22 = u_ll h_rr, M_con_1_rr, M_con_2_rr = u_rr @@ -162,16 +156,16 @@ end G_con_11 = G_cov_22 / G G_con_12 = -G_cov_12 / G G_con_22 = G_cov_11 / G - - G_nn = n_con_1 * G_con_11 * n_con_1 + n_con_1 * G_con_12 * n_con_2 + - n_con_2 * G_con_12 * n_con_1 + n_con_2 * G_con_22 * n_con_2 + + G_nn = n_con_1 * G_con_11 * n_con_1 + n_con_1 * G_con_12 * n_con_2 + + n_con_2 * G_con_12 * n_con_1 + n_con_2 * G_con_22 * n_con_2 v_n_ll = (M_con_1_ll * n_con_1 + M_con_2_ll * n_con_2) / h_ll v_n_rr = (M_con_1_rr * n_con_1 + M_con_2_rr * n_con_2) / h_rr gh_ll = max(h_ll * equations.gravitational_acceleration, 0) gh_rr = max(h_rr * equations.gravitational_acceleration, 0) - + return max(abs(v_n_ll) + sqrt(G_nn * gh_ll), abs(v_n_rr) + sqrt(G_nn * gh_rr)) end @@ -179,7 +173,6 @@ end Maximum wave speeds with respect to the contravariant basis """ @inline function Trixi.max_abs_speeds(u, ::CovariantShallowWaterEquations2D) - h, M_con_1, M_con_2, G_cov_11, G_cov_12, G_cov_22 = u G = G_cov_11 * G_cov_22 - G_cov_12^2 @@ -188,15 +181,15 @@ end gh = max(h * equations.gravitational_acceleration, 0) - v_con_1 = M_con_1 / h + v_con_1 = M_con_1 / h v_con_2 = M_con_2 / h return abs(v_con_1) + sqrt(G_con_11 * gh), abs(v_con_2) + sqrt(G_con_22 * gh) end -@inline function cartesian2contravariant(u_cartesian, ::CovariantShallowWaterEquations2D, - i, j, element, cache) - +@inline function cartesian2contravariant(u_cartesian, + ::CovariantShallowWaterEquations2D, + i, j, element, cache) (; contravariant_vectors, inverse_jacobian) = cache.elements A11, A21, A31, A12, A22, A32 = view(cache.elements.jacobian_matrix, :, :, i, j, @@ -207,22 +200,21 @@ end G_cov_22 = A12 * A12 + A22 * A22 + A32 * A32 Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, j, - element) + element) Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, i, j, - element) + element) return SVector(u_cartesian[1], - inverse_jacobian[i, j, element] * (Ja11 * u_cartesian[2] + + inverse_jacobian[i, j, element] * (Ja11 * u_cartesian[2] + Ja12 * u_cartesian[3] + Ja13 * u_cartesian[4]), - inverse_jacobian[i, j, element] * (Ja21 * u_cartesian[2] + + inverse_jacobian[i, j, element] * (Ja21 * u_cartesian[2] + Ja22 * u_cartesian[3] + Ja23 * u_cartesian[4]), - G_cov_11, G_cov_12, G_cov_22) + G_cov_11, G_cov_12, G_cov_22) end @inline function contravariant2cartesian(u_node, ::CovariantShallowWaterEquations2D, - i, j, element, cache) - + i, j, element, cache) A11, A21, A31, A12, A22, A32 = view(cache.elements.jacobian_matrix, :, :, i, j, element) return SVector(u_node[1], @@ -230,5 +222,4 @@ end A21 * u_node[2] + A22 * u_node[3], A31 * u_node[2] + A32 * u_node[3]) end - end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 052ffe7..8e54cd8 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -10,7 +10,7 @@ const SECONDS_PER_DAY = 8.64f4 abstract type AbstractCovariantEquations2D{NVARS} <: AbstractEquations{2, NVARS} end abstract type AbstractCompressibleMoistEulerEquations{NDIMS, NVARS} <: - AbstractEquations{NDIMS, NVARS} end + AbstractEquations{NDIMS, NVARS} end export CovariantLinearAdvectionEquation2D include("covariant_advection.jl") diff --git a/src/meshes/p4est_cubed_sphere.jl b/src/meshes/p4est_cubed_sphere.jl index 1de97a8..07f0d9b 100644 --- a/src/meshes/p4est_cubed_sphere.jl +++ b/src/meshes/p4est_cubed_sphere.jl @@ -309,4 +309,4 @@ function calc_tree_node_coordinates_cubed_sphere_2D!(node_coordinates::AbstractA end end end -end # @muladd \ No newline at end of file +end # @muladd diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl index 68a0b72..e27af63 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl @@ -44,7 +44,8 @@ end # This function dispatches on the dimensions of the mesh and the equation (AbstractEquations{3}) function Trixi.init_elements(mesh::Union{P4estMesh{2, RealT}, T8codeMesh{2, RealT}}, - equations::Union{AbstractEquations{3}, AbstractCovariantEquations2D}, + equations::Union{AbstractEquations{3}, + AbstractCovariantEquations2D}, basis, ::Type{uEltype}) where {RealT <: Real, uEltype <: Real} diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl index f2ce452..0d3e8e7 100644 --- a/src/solvers/dgsem_p4est/dg.jl +++ b/src/solvers/dgsem_p4est/dg.jl @@ -5,4 +5,4 @@ static2val(::Trixi.StaticInt{N}) where {N} = Val{N}() indices...) SVector(ntuple(@inline(dim->contravariant_vectors[dim, index, indices...]), static2val(static_size(contravariant_vectors, Trixi.StaticInt(1))))) -end \ No newline at end of file +end diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl index dfb6b9d..54c36c4 100644 --- a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl @@ -1,7 +1,11 @@ +@muladd begin +#! format: noindent + # Weak-form kernel for 3D equations solved in 2D manifolds @inline function Trixi.weak_form_kernel!(du, u, element, - mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, + mesh::Union{StructuredMesh{2}, + UnstructuredMesh2D, P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations::AbstractEquations{3}, @@ -20,7 +24,8 @@ # Compute the contravariant flux by taking the scalar product of the # first contravariant vector Ja^1 and the flux vector - Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, j, + Ja11, Ja12, Ja13 = Trixi.get_contravariant_vector(1, contravariant_vectors, i, + j, element) contravariant_flux1 = Ja11 * flux1 + Ja12 * flux2 + Ja13 * flux3 for ii in eachnode(dg) @@ -31,7 +36,8 @@ # Compute the contravariant flux by taking the scalar product of the # second contravariant vector Ja^2 and the flux vector - Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, i, j, + Ja21, Ja22, Ja23 = Trixi.get_contravariant_vector(2, contravariant_vectors, i, + j, element) contravariant_flux2 = Ja21 * flux1 + Ja22 * flux2 + Ja23 * flux3 for jj in eachnode(dg) @@ -43,3 +49,4 @@ return nothing end +end # muladd diff --git a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl index d0cb092..ad8cbb8 100644 --- a/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl +++ b/src/solvers/dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl @@ -75,52 +75,6 @@ end return nothing end -""" - Weak form kernel for fluxes already in contravariant components -""" -@inline function Trixi.flux_differencing_kernel!(du, u, - element, - mesh::Union{StructuredMesh{2}, - StructuredMeshView{2}, - UnstructuredMesh2D, P4estMesh{2}, - T8codeMesh{2}}, - nonconservative_terms::False, - equations::AbstractCovariantEquations2D, - volume_flux, dg::DGSEM, cache, - alpha = true) - (; derivative_split) = dg.basis - (; inverse_jacobian) = cache.elements - - for j in eachnode(dg), i in eachnode(dg) - u_node = Trixi.get_node_vars(u, equations, dg, i, j, element) - J_node = 1 / inverse_jacobian[i, j, element] - - # x direction - for ii in (i + 1):nnodes(dg) - u_node_ii = Trixi.get_node_vars(u, equations, dg, ii, j, element) - J_avg = 0.5 * (J_node + 1 / inverse_jacobian[ii, j, element]) - flux1 = J_avg * volume_flux(u_node, u_node_ii, 1, equations) - Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[i, ii], flux1, - equations, dg, i, j, element) - Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[ii, i], flux1, - equations, dg, ii, j, element) - end - - # y direction - for jj in (j + 1):nnodes(dg) - u_node_jj = Trixi.get_node_vars(u, equations, dg, i, jj, element) - J_avg = 0.5 * (J_node + 1 / inverse_jacobian[i, jj, element]) - flux2 = J_avg * volume_flux(u_node, u_node_jj, 2, equations) - Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[j, jj], flux2, - equations, dg, i, j, element) - Trixi.multiply_add_to_node_vars!(du, alpha * derivative_split[jj, j], flux2, - equations, dg, i, jj, element) - end - end - - return nothing -end - """ Interface flux for the covariant form, where the numerical flux is computed in the direction of the 2D reference normal """ @@ -175,7 +129,8 @@ function Trixi.calc_interface_flux!(surface_flux_values, for node in eachnode(dg) Trixi.calc_interface_flux!(surface_flux_values, mesh, nonconservative_terms, - equations, surface_integral, dg, cache, interface, + equations, surface_integral, dg, cache, + interface, i_primary, j_primary, i_secondary, j_secondary, node, primary_direction, primary_element, node_secondary, secondary_direction, diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index 5977385..b32b2ee 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -3,4 +3,4 @@ include("dgsem_p4est/dg_2d_manifold_in_3d_cartesian.jl") include("dgsem_p4est/containers_2d_manifold_in_3d.jl") export cartesian2contravariant, covariant2cartesian -include("dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl") \ No newline at end of file +include("dgsem_p4est/dg_2d_manifold_in_3d_covariant.jl") diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 31f7a29..9d6b862 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -7,9 +7,9 @@ include("test_trixiatmo.jl") EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") -@trixiatmo_testset "elixir_euler_spherical_advection_cartesian" begin +@trixiatmo_testset "elixir_spherical_advection_cartesian" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_euler_spherical_advection_cartesian.jl"), + "elixir_spherical_advection_cartesian.jl"), l2=[ 8.44505073e-03, 8.23414117e-03, @@ -34,4 +34,27 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") end end +@trixiatmo_testset "elixir_spherical_advection_covariant" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_spherical_advection_covariant.jl"), + l2=[ + 9.71730194e-03, + 0.00000000e+00, + 0.00000000e+00, + ], + linf=[ + 6.98678150e-02, + 0.00000000e+00, + 0.00000000e+00, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + end # module