From 1a9d315945178ac109b32a37b8325ba62889eda8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 21 Nov 2023 15:29:34 +0100 Subject: [PATCH] First working version of AMR on a 2D cubed sphere --- .../elixir_euler_cubed_sphere_shell_amr.jl | 127 ++++++++++++++++++ src/solvers/dgsem_p4est/containers.jl | 14 +- src/solvers/dgsem_tree/indicators_2d.jl | 13 ++ 3 files changed, 148 insertions(+), 6 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell_amr.jl diff --git a/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell_amr.jl new file mode 100644 index 0000000000..f2818bd816 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_euler_cubed_sphere_shell_amr.jl @@ -0,0 +1,127 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear 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 + +initial_condition = initial_condition_advection_sphere + +mesh = Trixi.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, + source_terms = source_terms_convert_to_linear_advection) + +############################################################################### +# ODE solvers, callbacks etc. + +# Create ODE problem with time span from 0.0 to π +tspan = (0.0, pi) +ode = semidiscretize(semi, tspan) + +# 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_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) + +# AMR callback +amr_indicator = IndicatorLöhner(semi, variable = Trixi.density) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.02, + max_level = 2, max_threshold = 0.03) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 1, + adapt_initial_condition = true, + adapt_initial_condition_only_refine = true) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, + stepsize_callback, amr_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 = 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/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index ffb6eea5e9..2371805d9a 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -49,22 +49,24 @@ function Base.resize!(elements::P4estElementContainer, capacity) _inverse_jacobian, _surface_flux_values = elements n_dims = ndims(elements) + ndims_spa = size(elements.node_coordinates, 1) n_nodes = size(elements.node_coordinates, 2) n_variables = size(elements.surface_flux_values, 1) - resize!(_node_coordinates, n_dims * n_nodes^n_dims * capacity) + resize!(_node_coordinates, ndims_spa * n_nodes^n_dims * capacity) elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates), - (n_dims, ntuple(_ -> n_nodes, n_dims)..., + (ndims_spa, ntuple(_ -> n_nodes, n_dims)..., capacity)) - resize!(_jacobian_matrix, n_dims^2 * n_nodes^n_dims * capacity) + resize!(_jacobian_matrix, ndims_spa * n_dims * n_nodes^n_dims * capacity) elements.jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix), - (n_dims, n_dims, + (ndims_spa, n_dims, ntuple(_ -> n_nodes, n_dims)..., capacity)) - resize!(_contravariant_vectors, length(_jacobian_matrix)) + resize!(_contravariant_vectors, ndims_spa^2 * n_nodes^n_dims * capacity) elements.contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors), - size(elements.jacobian_matrix)) + (ndims_spa, ndims_spa, + ntuple(_ -> n_nodes, n_dims)..., capacity)) resize!(_inverse_jacobian, n_nodes^n_dims * capacity) elements.inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian), diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index 8333bb515d..86652d9526 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -260,6 +260,19 @@ function create_cache(typ::Type{IndicatorLöhner}, mesh, equations::AbstractEqua create_cache(typ, equations, dg.basis) end +# this method is used when the indicator is constructed as for AMR +function create_cache(typ::Type{IndicatorLöhner}, mesh::AbstractMesh{2}, equations::AbstractEquations{3}, + dg::DGSEM, cache) + @unpack basis = dg + alpha = Vector{real(basis)}() + + A = Array{real(basis), ndims(mesh)} + indicator_threaded = [A(undef, nnodes(basis), nnodes(basis)) + for _ in 1:Threads.nthreads()] + + return (; alpha, indicator_threaded) +end + function (löhner::IndicatorLöhner)(u::AbstractArray{<:Any, 4}, mesh, equations, dg::DGSEM, cache; kwargs...)