Covariant Advection in DGMulti

Project.toml

@@ -13,6 +13,7 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+StartUpDG = "472ebc20-7c99-4d4b-9470-8fde4e9faa0f"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -29,6 +30,7 @@ NLsolve = "4.5.1"
 Printf = "1"
 QuadGK = "2"
 Reexport = "1.0"
+StartUpDG = "1.2.0"
 Static = "0.8, 1"
 StaticArrayInterface = "1.5.1"
 StaticArrays = "1"

examples/elixir_spherical_advection_covariant_prismed_icosahedron.jl

@@ -0,0 +1,77 @@
+###############################################################################
+# DGSEM for the linear advection equation on a prismed icosahedral grid
+###############################################################################
+# To run a convergence test, use
+# convergence_test("../examples/elixir_spherical_advection_covariant_prismed_icosahedron.jl", 4, initial_refinement_level = 1)
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Spatial discretization
+
+initial_condition = initial_condition_gaussian
+
+equations = CovariantLinearAdvectionEquation2D(global_coordinate_system = GlobalCartesianCoordinates())
+
+###############################################################################
+# Build DG solver.
+
+tensor_polydeg = (2, 1)
+
+dg = DGMulti(element_type = Wedge(),
+             approximation_type = Polynomial(),
+             surface_flux = flux_central,
+             polydeg = tensor_polydeg)
+
+###############################################################################
+# Build mesh.
+
+initial_refinement_level = 3
+
+mesh = DGMultiMeshPrismIcosahedron(dg;
+    inner_radius = 0.999 * EARTH_RADIUS,
+    outer_radius = EARTH_RADIUS,
+    initial_refinement = initial_refinement_level)
+
+# Transform the initial condition to the proper set of conservative variables
+initial_condition_transformed = transform_initial_condition(initial_condition, equations)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_transformed, dg)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0 to T
+ode = semidiscretize(semi, (0.0, 12 * SECONDS_PER_DAY))
+
+# 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,),
+                                     uEltype = real(dg))
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 10,
+                                     solution_variables = contravariant2global)
+
+# 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, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, save_everystep = true, callback = callbacks)

src/TrixiAtmo.jl

@@ -15,7 +15,7 @@ using Printf: @sprintf
 using Static: True, False
 using StrideArrays: PtrArray
 using StaticArrayInterface: static_size
-using LinearAlgebra: cross, norm, dot, det
+using LinearAlgebra: Diagonal, cross, norm, dot, det
 using Reexport: @reexport
 using LoopVectorization: @turbo
 using QuadGK: quadgk
@@ -35,6 +35,9 @@ using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
 
 using Trixi: ln_mean, stolarsky_mean, inv_ln_mean
 
+# DGMulti solvers
+using StartUpDG: RefElemData, MeshData, AbstractElemShape
+
 include("auxiliary/auxiliary.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")
@@ -69,9 +72,9 @@ export source_terms_lagrange_multiplier, clean_solution_lagrange_multiplier!
 
 export cons2prim_and_vorticity, contravariant2global
 
-export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, MetricTermsCrossProduct,
-       MetricTermsInvariantCurl, MetricTermsCovariantSphere, ChristoffelSymbolsAutodiff,
-       ChristoffelSymbolsCollocationDerivative
+export P4estMeshCubedSphere2D, P4estMeshQuadIcosahedron2D, DGMultiMeshPrismIcosahedron,
+       MetricTermsCrossProduct, MetricTermsInvariantCurl, MetricTermsCovariantSphere,
+       ChristoffelSymbolsAutodiff, ChristoffelSymbolsCollocationDerivative
 
 export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
        EARTH_ROTATION_RATE, SECONDS_PER_DAY

src/callbacks_step/analysis_covariant.jl

@@ -77,6 +77,29 @@ function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t,
     end
 end
 
+# Entropy time derivative for cons2entropy function which depends on auxiliary variables
+function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t,
+                 mesh::DGMultiMesh, equations::AbstractCovariantEquations, dg::DGMulti, cache)
+    rd = dg.basis
+    md = mesh.md
+    @unpack u_values, aux_quad_values = cache
+
+    # interpolate u, du to quadrature points
+    du_values = similar(u_values) # TODO: DGMulti. Can we move this to the analysis cache somehow?
+    Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), du_values, du)
+    Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u)
+
+    # compute ∫v(u) * du/dt = ∫dS/dt. We can directly compute v(u) instead of computing the entropy
+    # projection here, since the RHS will be projected to polynomials of degree N and testing with
+    # the L2 projection of v(u) would be equivalent to testing with v(u) due to the moment-preserving
+    # property of the L2 projection.
+    dS_dt = zero(eltype(first(du)))
+    for i in Base.OneTo(length(md.wJq))
+        dS_dt += dot(cons2entropy(u_values[i], aux_quad_values[i], equations), du_values[i]) * md.wJq[i]
+    end
+    return dS_dt
+end
+
 # L2 and Linf error calculation for the covariant form
 function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
                                 equations::AbstractCovariantEquations{2},
@@ -124,4 +147,27 @@ function Trixi.calc_error_norms(func, u, t, analyzer, mesh::P4estMesh{2},
 
     return l2_error, linf_error
 end
+
+# L2 and Linf error calculation for the covariant form
+function Trixi.calc_error_norms(func, u, t, analyzer,
+                                mesh::DGMultiMesh{NDIMS}, equations::AbstractCovariantEquations, initial_condition,
+                                dg::DGMulti{NDIMS}, cache, cache_analysis) where {NDIMS}
+    rd = dg.basis
+    md = mesh.md
+    @unpack u_values, aux_quad_values = cache
+
+    # interpolate u to quadrature points
+    Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u)
+
+    component_l2_errors = zero(eltype(u_values))
+    component_linf_errors = zero(eltype(u_values))
+    for i in Trixi.each_quad_node_global(mesh, dg, cache)
+        u_exact = initial_condition(SVector(getindex.(md.xyzq, i)), t, aux_quad_values[i], equations)
+        error_at_node = func(u_values[i], equations) - func(u_exact, equations)
+        component_l2_errors += md.wJq[i] * error_at_node .^ 2
+        component_linf_errors = max.(component_linf_errors, abs.(error_at_node))
+    end
+    total_volume = sum(md.wJq)
+    return sqrt.(component_l2_errors ./ total_volume), component_linf_errors
+end
 end # @muladd

src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl

@@ -118,6 +118,40 @@ function Trixi.save_solution_file(u, time, dt, timestep,
     return filename
 end
 
+function Trixi.save_solution_file(u, time, dt, timestep,
+                                  mesh::DGMultiMesh,
+                                  equations::AbstractCovariantEquations,
+                                  dg::DG, cache,
+                                  solution_callback,
+                                  element_variables = Dict{Symbol, Any}(),
+                                  node_variables = Dict{Symbol, Any}();
+                                  system = "")
+    @unpack output_directory, solution_variables = solution_callback
+
+    # Filename based on current time step.
+    if isempty(system)
+        filename = joinpath(output_directory, @sprintf("solution_%09d.vtu", timestep))
+    else
+        filename = joinpath(output_directory,
+                            @sprintf("solution_%s_%09d.vtu", system, timestep))
+    end
+
+    # Convert to different set of variables if requested.
+    if solution_variables === cons2cons
+        data = u
+        n_vars = nvariables(equations)
+    else
+        data = convert_variables(u, solution_variables, mesh, equations, dg, cache)
+        # Find out variable count by looking at output from `solution_variables` function.
+        n_vars = length(data[1])
+    end
+    # Create a dictionary with variable names and data
+    var_names = Trixi.varnames(solution_variables, equations)
+    var_dict = Dict(var_names[i] => getindex.(data, i) for i in 1:n_vars)
+    StartUpDG.export_to_vtk(dg.basis, mesh.md, var_dict, filename;
+                            equi_dist_nodes = false)
+end
+
 # Calculate the primitive variables and the relative vorticity at a given node
 @inline function cons2prim_and_vorticity(u, normal_vector,
                                          mesh::P4estMesh{2},

src/callbacks_step/save_solution_covariant.jl

@@ -33,6 +33,27 @@ function convert_variables(u, solution_variables, mesh::P4estMesh{2},
     return data
 end
 
+# Convert to another set of variables using a solution_variables function
+function convert_variables(u, solution_variables, mesh::DGMultiMesh,
+                           equations::AbstractCovariantEquations, dg, cache)
+    (; aux_values) = cache
+    # Extract the number of solution variables to be output 
+    # (may be different than the number of conservative variables) 
+    n_vars = length(Trixi.varnames(solution_variables, equations))
+
+    # Allocate storage for output, an Np x n_elements array of nvars arrays
+    # data = Array{eltype(u)}(undef, n_vars, (size(aux_node_vars)[2:end])...)
+    data = map(_ -> SVector{n_vars, eltype(u)}(zeros(n_vars)), u)
+
+    # Loop over all nodes and convert variables, passing in auxiliary variables
+    for i in Trixi.each_dof_global(mesh, dg)
+        u_node = u[:, i]
+        aux_node = aux_values[i]
+        data[i] = solution_variables(u_node, aux_node, equations)
+    end
+    return data
+end
+
 # Calculate the primitive variables and relative vorticity at a given node. The velocity
 # components in the global coordinate system and the bottom topography are returned, such 
 # that the outputs for CovariantShallowWaterEquations2D and ShallowWaterEquations3D are the 

src/callbacks_step/stepsize_dg2d.jl

@@ -69,4 +69,31 @@ function Trixi.max_dt(u, t, mesh::P4estMesh{2}, constant_speed::False,
     end
     return 2 / (nnodes(dg) * max_scaled_speed)
 end
+
+function Trixi.max_dt(u, t, mesh::DGMultiMesh,
+                constant_speed::False, equations::AbstractCovariantEquations{NDIMS},
+                dg::DGMulti, cache) where {NDIMS}
+    @unpack md = mesh
+    rd = dg.basis
+    (; aux_values) = cache
+
+    dt_min = Inf
+    for e in eachelement(mesh, dg, cache)
+        max_speeds = ntuple(_ -> nextfloat(zero(t)), NDIMS)
+        for i in 1:nnodes(dg)
+            u_node = u[i, e]
+            aux_node = aux_values[i, e]
+            detg = area_element(aux_node, equations)
+            lambda_i = max_abs_speeds(u_node, aux_node, equations)
+            max_speeds = max.(max_speeds, lambda_i)
+        end
+        dt_min = min(dt_min, 1 / sum(max_speeds))
+    end
+
+    # This mimics `max_dt` for `TreeMesh`, except that `nnodes(dg)` is replaced by
+    # `polydeg+1`. This is because `nnodes(dg)` returns the total number of
+    # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns
+    # the number of 1D nodes for `DGSEM` solvers.
+    return 2 * dt_min * Trixi.dt_polydeg_scaling(dg)
+end
 end # @muladd

src/equations/covariant_advection.jl

@@ -96,6 +96,15 @@ end
     return SVector(J * u[1] * vcon[orientation], z, z)
 end
 
+@inline function Trixi.flux(u, aux_vars, normal_direction::AbstractVector,
+                      equations::CovariantLinearAdvectionEquation2D)
+    z = zero(eltype(u))
+    J = area_element(aux_vars, equations)
+    vcon = velocity_contravariant(u, equations)
+    a = dot(vcon, normal_direction) # velocity in normal direction
+    return SVector(J * u[1] * a, z, z)
+end
+
 # Local Lax-Friedrichs dissipation which is not applied to the contravariant velocity 
 # components, as they should remain unchanged in time
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, aux_vars_ll,
@@ -118,6 +127,17 @@ end
     return max(abs(vcon_ll[orientation]), abs(vcon_rr[orientation]))
 end
 
+@inline function Trixi.max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr,
+                                     normal_direction::AbstractVector,
+                                     equations::CovariantLinearAdvectionEquation2D)
+    vcon_ll = velocity_contravariant(u_ll, equations)  # Contravariant components on left side
+    vcon_rr = velocity_contravariant(u_rr, equations)  # Contravariant components on right side
+    # Calculate the velocity in the normal direction
+    a_ll = abs(dot(vcon_ll, normal_direction))
+    a_rr = abs(dot(vcon_rr, normal_direction))
+    return max(a_ll, a_rr)
+end
+
 # Maximum wave speeds in each direction for CFL calculation
 @inline function max_abs_speeds(u, aux_vars,
                                 equations::CovariantLinearAdvectionEquation2D)