Covariant solver in DGMulti

Project.toml

@@ -13,6 +13,8 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+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 +31,7 @@ NLsolve = "4.5.1"
 Printf = "1"
 QuadGK = "2"
 Reexport = "1.0"
+Setfield = "1.1.2"
 Static = "0.8, 1"
 StaticArrayInterface = "1.5.1"
 StaticArrays = "1"

examples/advection/covariant/elixir_tri_icosahedron.jl

@@ -0,0 +1,75 @@
+###############################################################################
+# 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 = 3
+
+dg = DGMulti(element_type = Tri(),
+             approximation_type = Polynomial(),
+             surface_flux = flux_lax_friedrichs,
+             polydeg = tensor_polydeg)
+
+###############################################################################
+# Build mesh.
+
+initial_refinement_level = 4
+
+mesh = DGMultiMeshTriIcosahedron2D(dg, 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 = 100,
+                                     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 = 100,
+                                     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 = false, callback = callbacks)

examples/shallow_water/covariant/elixir_tri_barotropic_instability.jl

@@ -0,0 +1,77 @@
+###############################################################################
+# DGSEM for the linear advection equation on a prismed icosahedral grid
+###############################################################################
+
+using OrdinaryDiffEq, Trixi, TrixiAtmo
+
+###############################################################################
+# Spatial discretization
+
+initial_condition = initial_condition_barotropic_instability
+
+equations = CovariantShallowWaterEquations2D(EARTH_GRAVITATIONAL_ACCELERATION,
+                                             EARTH_ROTATION_RATE,
+                                             global_coordinate_system = GlobalCartesianCoordinates())
+
+###############################################################################
+# Build DG solver.
+
+polydeg = 4
+
+dg = DGMulti(element_type = Tri(),
+             approximation_type = Polynomial(),
+             surface_flux = flux_lax_friedrichs,
+             polydeg = polydeg)
+
+###############################################################################
+# Build mesh.
+
+initial_refinement_level = 4
+
+mesh = DGMultiMeshTriIcosahedron2D(dg, 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,
+                                    source_terms = source_terms_geometric_coriolis)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0 to T
+tspan = (0.0, 12.0 * SECONDS_PER_DAY)
+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 = 100,
+                                     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 = 300,
+                                     solution_variables = cons2prim_and_vorticity)
+
+# 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 = false, callback = callbacks)

src/TrixiAtmo.jl

@@ -15,12 +15,13 @@ using Printf: @sprintf
 using Static: True, False
 using StrideArrays: PtrArray
 using StaticArrayInterface: static_size
-using LinearAlgebra: cross, norm, dot, det
+using LinearAlgebra: Diagonal, I, cross, norm, dot, det, diagm
 using Reexport: @reexport
 using LoopVectorization: @turbo
 using QuadGK: quadgk
 using ForwardDiff: derivative
 using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
+using Setfield
 
 @reexport using StaticArrays: SVector, SMatrix
 @reexport import Trixi: waterheight, varnames, cons2cons, cons2prim,
@@ -35,6 +36,9 @@ using HDF5: HDF5, h5open, attributes, create_dataset, datatype, dataspace
 
 using Trixi: ln_mean, stolarsky_mean, inv_ln_mean
 
+# DGMulti Solvers
+using StartUpDG: MeshData, RefElemData
+
 include("auxiliary/auxiliary.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")
@@ -74,9 +78,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, DGMultiMeshTriIcosahedron2D
+MetricTermsCrossProduct, MetricTermsInvariantCurl, MetricTermsCovariantSphere,
+ChristoffelSymbolsAutodiff, ChristoffelSymbolsCollocationDerivative
 
 export EARTH_RADIUS, EARTH_GRAVITATIONAL_ACCELERATION,
        EARTH_ROTATION_RATE, SECONDS_PER_DAY

src/callbacks_step/analysis_covariant.jl

@@ -21,6 +21,30 @@ function Trixi.integrate(func::Func, u,
     end
 end
 
+function Trixi.integrate(func::Func, u,
+                         mesh::DGMultiMesh,
+                         equations::AbstractCovariantEquations,
+                         dg::DGMulti, cache; normalize = true) where {Func}
+    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)
+
+    integral = zero(func(u_values[1], equations))
+    total_volume = zero(sum(rd.wq))
+    for element in Trixi.eachelement(dg, cache)
+        weights = area_element.(aux_quad_values[:, element], equations) .* rd.wq
+        integral += sum(weights .* func.(u_values[:, element], equations))
+        total_volume += sum(weights)
+    end
+    if normalize == true
+        integral /= total_volume
+    end
+    return integral
+end
+
 # For the covariant form, we want to integrate using the exact area element 
 # J = √G = (det(AᵀA))^(1/2), which is stored in cache.auxiliary_variables, not the approximate 
 # area element used in the Cartesian formulation, which stored in cache.elements
@@ -78,6 +102,33 @@ 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)
+    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))
+        ref_index = mod(i - 1, rd.Nq) + 1
+        node_weight = rd.wq[ref_index] * area_element(aux_quad_values[i], equations)
+        dS_dt += dot(cons2entropy(u_values[i], aux_quad_values[i], equations),
+                     du_values[i]) * node_weight
+    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},
@@ -125,4 +176,34 @@ 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_AMBIENT},
+                                equations::AbstractCovariantEquations,
+                                initial_condition,
+                                dg::DGMulti{NDIMS}, cache,
+                                cache_analysis) where {NDIMS, NDIMS_AMBIENT}
+    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))
+    total_volume = zero(eltype(u_values[1]))
+    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)
+        ref_index = mod(i - 1, rd.Nq) + 1
+        node_weight = rd.wq[ref_index] * area_element(aux_quad_values[i], equations)
+        component_l2_errors += node_weight * error_at_node .^ 2
+        component_linf_errors = max.(component_linf_errors, abs.(error_at_node))
+        total_volume += node_weight
+    end
+    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 = true)
+end
+
 # Calculate the primitive variables and the relative vorticity at a given node
 @inline function cons2prim_and_vorticity(u, normal_vector,
                                          mesh::P4estMesh{2},
@@ -177,6 +211,28 @@ end
     return SVector(dv3dy - dv2dz, dv1dz - dv3dx, dv2dx - dv1dy)
 end
 
+@inline function calc_vorticity_node(u, equations::AbstractEquations,
+                                     dg::DGMulti, cache, i, element)
+    rd = dg.basis
+    (; aux_values) = cache
+    Dr, Ds = rd.Drst
+
+    dv2dxi1 = dv1dxi2 = zero(eltype(u))
+    for ii in eachnode(dg)
+        index = (element - 1) * nnodes(dg) + ii
+        u_node_ii = u[:, index]
+        aux_node_ii = aux_values[ii, element]
+        vcov = metric_covariant(aux_node_ii, equations) *
+               velocity_contravariant(u_node_ii, equations)
+        dv2dxi1 = dv2dxi1 + Dr[i, ii] * vcov[2]
+        dv1dxi2 = dv1dxi2 + Ds[i, ii] * vcov[1]
+    end
+
+    # compute the relative vorticity
+    aux_node = aux_values[i, element]
+    return (dv2dxi1 - dv1dxi2) / area_element(aux_node, equations)
+end
+
 # Variable names for cons2prim_and_vorticity
 function Trixi.varnames(::typeof(cons2prim_and_vorticity),
                         equations::ShallowWaterEquations3D)