diff --git a/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl index 35272fec96..f9bf18e977 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_blast_wave_amr.jl @@ -1,4 +1,5 @@ +using Downloads: download using OrdinaryDiffEq using Trixi using Plots @@ -56,7 +57,7 @@ isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/a07 mesh = P4estMesh{2}(mesh_file, polydeg=3, initial_refinement_level=1) -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=Dict( :all => BoundaryConditionDirichlet(initial_condition) )) diff --git a/examples/p4est_2d_dgsem/elixir_euler_sedov.jl b/examples/p4est_2d_dgsem/elixir_euler_sedov.jl index e2459dc6d6..9f5247e8c4 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_sedov.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi diff --git a/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec.jl b/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec.jl index 49f1197aaa..6ef551c486 100644 --- a/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec.jl +++ b/examples/p4est_2d_dgsem/elixir_euler_shockcapturing_ec.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -60,7 +59,7 @@ save_solution = SaveSolutionCallback(interval=100, stepsize_callback = StepsizeCallback(cfl=1.0) callbacks = CallbackSet(summary_callback, - analysis_callback, + analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl b/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl index f932dd21c6..3ce4f8da13 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_baroclinic_instability.jl @@ -16,49 +16,103 @@ using Trixi using LinearAlgebra ############################################################################### -# semidiscretization of the compressible Euler equations +# Setup for the baroclinic instability test gamma = 1.4 equations = CompressibleEulerEquations3D(gamma) # Initial condition for an idealized baroclinic instability test # https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A function initial_condition_baroclinic_instability(x, t, equations::CompressibleEulerEquations3D) - # Parameters from Table 1 in the paper - # Corresponding names in the paper are commented - radius_earth = 6.371229e6 # a - half_width_parameter = 2 # b - gravity_constant = 9.80616 # g - k = 3 # k - p0 = 1e5 # p₀ - gas_constant = 287 # R - temperature0e = 310.0 # T₀ᴱ - temperature0p = 240.0 # T₀ᴾ - lapse = 0.005 # Γ - omega_earth = 7.29212e-5 # Ω + lon, lat, r = cartesian_to_sphere(x) + radius_earth = 6.371229e6 + # Make sure that the r is not smaller than radius_earth + z = max(r - radius_earth, 0.0) + + # Unperturbated basic state + rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + + # Stream function type perturbation + u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z) + + u += u_perturbation + v = v_perturbation + # Convert spherical velocity to Cartesian + v1 = -sin(lon) * u - sin(lat) * cos(lon) * v + v2 = cos(lon) * u - sin(lat) * sin(lon) * v + v3 = cos(lat) * v + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +# Steady state for RHS correction below +function steady_state_baroclinic_instability(x, t, equations::CompressibleEulerEquations3D) lon, lat, r = cartesian_to_sphere(x) + radius_earth = 6.371229e6 # Make sure that the r is not smaller than radius_earth z = max(r - radius_earth, 0.0) + + # Unperturbated basic state + rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + + # Convert spherical velocity to Cartesian + v1 = -sin(lon) * u + v2 = cos(lon) * u + v3 = 0.0 + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end + +function cartesian_to_sphere(x) + r = norm(x) + lambda = atan(x[2], x[1]) + if lambda < 0 + lambda += 2 * pi + end + phi = asin(x[3] / r) + + return lambda, phi, r +end + +# Unperturbated balanced steady-state. +# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p). +# The other velocity components are zero. +function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + radius_earth = 6.371229e6 # a + half_width_parameter = 2 # b + gravitational_acceleration = 9.80616 # g + k = 3 # k + surface_pressure = 1e5 # p₀ + gas_constant = 287 # R + surface_equatorial_temperature = 310.0 # T₀ᴱ + surface_polar_temperature = 240.0 # T₀ᴾ + lapse_rate = 0.005 # Γ + angular_velocity = 7.29212e-5 # Ω + + # Distance to the center of the Earth r = z + radius_earth # In the paper: T₀ - temperature0 = 0.5 * (temperature0e + temperature0p) + temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature) # In the paper: A, B, C, H - const_a = 1 / lapse - const_b = (temperature0 - temperature0p) / (temperature0 * temperature0p) - const_c = 0.5 * (k + 2) * (temperature0e - temperature0p) / (temperature0e * temperature0p) - const_h = gas_constant * temperature0 / gravity_constant + const_a = 1 / lapse_rate + const_b = (temperature0 - surface_polar_temperature) / + (temperature0 * surface_polar_temperature) + const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) / + (surface_equatorial_temperature * surface_polar_temperature) + const_h = gas_constant * temperature0 / gravitational_acceleration # In the paper: (r - a) / bH scaled_z = z / (half_width_parameter * const_h) # Temporary variables - temp1 = exp(lapse/temperature0 * z) + temp1 = exp(lapse_rate/temperature0 * z) temp2 = exp(-scaled_z^2) # In the paper: ̃τ₁, ̃τ₂ - tau1 = const_a * lapse / temperature0 * temp1 + const_b * (1 - 2 * scaled_z^2) * temp2 + tau1 = const_a * lapse_rate / temperature0 * temp1 + const_b * (1 - 2 * scaled_z^2) * temp2 tau2 = const_c * (1 - 2 * scaled_z^2) * temp2 # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr' @@ -73,95 +127,89 @@ function initial_condition_baroclinic_instability(x, t, equations::CompressibleE temperature = 1 / ((r/radius_earth)^2 * (tau1 - tau2 * temp4)) # In the paper: U, u (zonal wind, first component of spherical velocity) - big_u = gravity_constant/radius_earth * k * temperature * inttau2 * (temp3^(k-1) - temp3^(k+1)) + big_u = gravitational_acceleration/radius_earth * k * temperature * inttau2 * (temp3^(k-1) - temp3^(k+1)) temp5 = radius_earth * cos(lat) - u = -omega_earth * temp5 + sqrt(omega_earth^2 * temp5^2 + temp5 * big_u) + u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u) # Hydrostatic pressure - p = p0 * exp(-gravity_constant/gas_constant * (inttau1 - inttau2 * temp4)) - - # Perturbation - u += evaluate_exponential(lon,lat,z) - - # Convert spherical velocity to Cartesian - v1 = -sin(lon) * u - v2 = cos(lon) * u - v3 = 0.0 + p = surface_pressure * exp(-gravitational_acceleration/gas_constant * (inttau1 - inttau2 * temp4)) # Density (via ideal gas law) rho = p / (gas_constant * temperature) - return prim2cons(SVector(rho, v1, v2, v3, p), equations) + return rho, u, p end -function cartesian_to_sphere(x) - r = norm(x) - lambda = atan(x[2], x[1]) - if lambda < 0 - lambda += 2 * pi +# Perturbation as in Equations 25 and 26 of the paper (analytical derivative) +function perturbation_stream_function(lon, lat, z) + # Parameters from Table 1 in the paper + # Corresponding names in the paper are commented + perturbation_radius = 1/6 # d₀ / a + perturbed_wind_amplitude = 1.0 # Vₚ + perturbation_lon = pi/9 # Longitude of perturbation location + perturbation_lat = 2 * pi/9 # Latitude of perturbation location + pertz = 15000 # Perturbation height cap + + # Great circle distance (d in the paper) divided by a (radius of the Earth) + # because we never actually need d without dividing by a + great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) + + cos(perturbation_lat) * cos(lat) * + cos(lon - perturbation_lon)) + + # In the first case, the vertical taper function is per definition zero. + # In the second case, the stream function is per definition zero. + if z > pertz || great_circle_distance_by_a > perturbation_radius + return 0.0, 0.0 end - phi = asin(x[3] / r) - return lambda, phi, r -end + # Vertical tapering of stream function + perttaper = 1.0 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 -# Exponential perturbation function, taken from Fortran initialisation routine provided with -# https://doi.org/10.1002/qj.2241 -function evaluate_exponential(lon, lat, z) - pertup = 1.0 - pertexpr = 0.1 - pertlon = pi/9.0 - pertlat = 2.0*pi/9.0 - pertz = 15000.0 + # sin/cos(pi * d / (2 * d_0)) in the paper + sin_, cos_ = sincos(0.5 * pi * great_circle_distance_by_a / perturbation_radius) - # Great circle distance - greatcircler = 1/pertexpr * acos( - sin(pertlat) * sin(lat) + cos(pertlat) * cos(lat) * cos(lon - pertlon)) + # Common factor for both u and v + factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_ - # Vertical tapering of stream function - if z < pertz - perttaper = 1.0 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3 - else - perttaper = 0.0 - end + u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) + + cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon) + ) / sin(great_circle_distance_by_a) - # Zonal velocity perturbation - if greatcircler < 1 - result = pertup * perttaper * exp(- greatcircler^2) - else - result = 0.0 - end + v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) / + sin(great_circle_distance_by_a) - return result + return u_perturbation, v_perturbation end @inline function source_terms_baroclinic_instability(u, x, t, equations::CompressibleEulerEquations3D) - radius_earth = 6.371229e6 # a - gravity_constant = 9.80616 # g - omega_earth = 7.29212e-5 # Ω + radius_earth = 6.371229e6 # a + gravitational_acceleration = 9.80616 # g + angular_velocity = 7.29212e-5 # Ω r = norm(x) - # Make sure that the r is not smaller than radius_earth + # Make sure that r is not smaller than radius_earth z = max(r - radius_earth, 0.0) r = z + radius_earth du1 = zero(eltype(u)) # Gravity term - temp = -gravity_constant * radius_earth^2 / r^3 + temp = -gravitational_acceleration * radius_earth^2 / r^3 du2 = temp * u[1] * x[1] du3 = temp * u[1] * x[2] du4 = temp * u[1] * x[3] du5 = temp * (u[2] * x[1] + u[3] * x[2] + u[4] * x[3]) - # Coriolis term, -2Ω × ρv = -2 * omega_earth * (0, 0, 1) × u[2:4] - du2 -= -2 * omega_earth * u[3] - du3 -= 2 * omega_earth * u[2] + # Coriolis term, -2Ω × ρv = -2 * angular_velocity * (0, 0, 1) × u[2:4] + du2 -= -2 * angular_velocity * u[3] + du3 -= 2 * angular_velocity * u[2] return SVector(du1, du2, du3, du4, du5) end +############################################################################### +# Start of the actual elixir, semidiscretization of the problem initial_condition = initial_condition_baroclinic_instability @@ -190,7 +238,35 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, # ODE solvers, callbacks etc. tspan = (0.0, 10 * 24 * 60 * 60.0) # time in seconds for 10 days -ode = semidiscretize(semi, tspan) + +# Save RHS of the steady state and subtract it in every RHS evaluation. +# This trick preserves the steady state exactly (to machine rounding errors, of course). +# Otherwise, this elixir produces entirely unusable results for a resolution of 8x8x4 cells +# per cube face with a polydeg of 3. +# With this trick, even the polydeg 3 simulation produces usable (although badly resolved) results, +# and most of the grid imprinting in higher polydeg simulation is eliminated. +# +# See https://github.com/trixi-framework/Trixi.jl/issues/980 for more information. +u_steady_state = compute_coefficients(steady_state_baroclinic_instability, tspan[1], semi) +# Use a `let` block for performance (otherwise du_steady_state will be a global variable) +let du_steady_state = similar(u_steady_state) + # Save RHS of the steady state + Trixi.rhs!(du_steady_state, u_steady_state, semi, tspan[1]) + + global function corrected_rhs!(du, u, semi, t) + # Normal RHS evaluation + Trixi.rhs!(du, u, semi, t) + # Correct by subtracting the steady-state RHS + Trixi.@trixi_timeit Trixi.timer() "rhs correction" begin + # Use Trixi.@threaded for threaded performance + Trixi.@threaded for i in eachindex(du) + du[i] -= du_steady_state[i] + end + end + end +end +u0 = compute_coefficients(tspan[1], semi) +ode = ODEProblem(corrected_rhs!, u0, tspan, semi) summary_callback = SummaryCallback() diff --git a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl index 5f5d895c08..00da413285 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_sedov.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_sedov.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -27,7 +26,7 @@ function initial_condition_medium_sedov_blast_wave(x, t, equations::Compressible r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) E = 1.0 p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) - p0_outer = 1.0e-3 + p0_outer = 1.0e-3 # Calculate primitive variables rho = 1.0 @@ -42,7 +41,7 @@ end initial_condition = initial_condition_medium_sedov_blast_wave surface_flux = flux_lax_friedrichs -volume_flux = flux_ranocha +volume_flux = flux_ranocha polydeg = 5 basis = LobattoLegendreBasis(polydeg) indicator_sc = IndicatorHennemannGassner(equations, basis, @@ -53,8 +52,8 @@ indicator_sc = IndicatorHennemannGassner(equations, basis, volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg=volume_flux, volume_flux_fv=surface_flux) - -solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral) + +solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volume_integral) coordinates_min = (-1.0, -1.0, -1.0) coordinates_max = ( 1.0, 1.0, 1.0) diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl index 67194a521a..dcc2bb1c7c 100644 --- a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonconforming_unstructured_curved.jl @@ -1,4 +1,5 @@ +using Downloads: download using OrdinaryDiffEq using Trixi diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl index 555c907098..69f4e71077 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl index 8e6d9f331a..f2841514ab 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl index 6d7309ad0b..f2433e1555 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi diff --git a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl index b072178627..3d5a391bd9 100644 --- a/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl +++ b/examples/unstructured_2d_dgsem/elixir_euler_sedov.jl @@ -1,5 +1,4 @@ -using Downloads: download using OrdinaryDiffEq using Trixi @@ -56,7 +55,8 @@ solver = DGSEM(polydeg=polydeg, surface_flux=surface_flux, volume_integral=volum # Get the curved quad mesh from a file default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh") - +isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh", + default_mesh_file) mesh_file = default_mesh_file mesh = UnstructuredMesh2D(mesh_file, periodicity=true) diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 22dfba764b..84660bedcc 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -122,8 +122,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_euler_baroclinic_instability.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_baroclinic_instability.jl"), - l2 = [2.09897807904154e-6, 0.0008762268315931632, 0.002292983258753421, 0.00014890463331164162, 0.44731850924464284], - linf = [0.00018844037951093462, 0.14135647414705824, 0.41606363387428025, 0.011495022607007977, 54.019175442983396], + l2 = [6.725065410642336e-7, 0.00021710117340245454, 0.000438679759422352, 0.00020836356588024185, 0.07602006689579247], + linf = [1.9101671995258585e-5, 0.029803626911022396, 0.04847630924006063, 0.022001371349740104, 4.847761006938526], tspan = (0.0, 1e2), # Decrease tolerance of adaptive time stepping to get similar results across different systems abstol=1.0e-9, reltol=1.0e-9,