diff --git a/Project.toml b/Project.toml index a0e461f9..63bd61ab 100644 --- a/Project.toml +++ b/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" diff --git a/examples/advection/covariant/elixir_tri_icosahedron.jl b/examples/advection/covariant/elixir_tri_icosahedron.jl new file mode 100644 index 00000000..3a362aa2 --- /dev/null +++ b/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) diff --git a/examples/shallow_water/covariant/elixir_tri_barotropic_instability.jl b/examples/shallow_water/covariant/elixir_tri_barotropic_instability.jl new file mode 100644 index 00000000..c52bd096 --- /dev/null +++ b/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) diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index f8ea19e4..909804eb 100644 --- a/src/TrixiAtmo.jl +++ b/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 diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index 57f95386..06904218 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/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 diff --git a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl index c55514e9..baabb9d7 100644 --- a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl +++ b/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) diff --git a/src/callbacks_step/save_solution_covariant.jl b/src/callbacks_step/save_solution_covariant.jl index db2e9c1b..677ff8a4 100644 --- a/src/callbacks_step/save_solution_covariant.jl +++ b/src/callbacks_step/save_solution_covariant.jl @@ -33,6 +33,30 @@ 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 = 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) + if applicable(solution_variables, u, mesh, equations, dg, cache, i) + data[i] = solution_variables(u, mesh, equations, dg, cache, i) + else + u_node = u[:, i] + aux_node = aux_values[i] + data[i] = solution_variables(u_node, aux_node, equations) + end + 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 @@ -49,6 +73,22 @@ end return SVector(primitive_global..., h_s, relative_vorticity) end +@inline function cons2prim_and_vorticity(u, mesh::DGMultiMesh, + equations::AbstractCovariantEquations{2}, + dg::DGMulti, cache, i) + (; aux_values) = cache + u_node = u[:, i] + aux_node = aux_values[i] + element = (i - 1) ÷ nnodes(dg) + 1 + local_node_index = (i - 1) % nnodes(dg) + 1 + relative_vorticity = calc_vorticity_node(u, equations, dg, cache, local_node_index, + element) + h_s = bottom_topography(aux_node, equations) + primitive_global = contravariant2global(cons2prim(u_node, aux_node, equations), + aux_node, equations) + return SVector(primitive_global..., h_s, relative_vorticity) +end + # Calculate relative vorticity ζ = (∂₁v₂ - ∂₂v₁)/J for equations in covariant form @inline function calc_vorticity_node(u, equations::AbstractCovariantEquations{2}, dg::DGSEM, cache, i, j, element) diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index 0d9c888d..b5e5069c 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/src/callbacks_step/stepsize_dg2d.jl @@ -69,4 +69,30 @@ 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} + (; 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 = Trixi.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 diff --git a/src/equations/covariant_advection.jl b/src/equations/covariant_advection.jl index 84e998b5..32fac313 100644 --- a/src/equations/covariant_advection.jl +++ b/src/equations/covariant_advection.jl @@ -96,6 +96,15 @@ end return SVector(J * u[1] * vcon[orientation], z, z) end +@inline function 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 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) diff --git a/src/equations/covariant_shallow_water.jl b/src/equations/covariant_shallow_water.jl index 616cc749..0be3c270 100644 --- a/src/equations/covariant_shallow_water.jl +++ b/src/equations/covariant_shallow_water.jl @@ -186,6 +186,32 @@ end return SVector(J * h_vcon[orientation], J * momentum_flux_1, J * momentum_flux_2) end +# Flux as a function of the state vector u, as well as the auxiliary variables aux_vars, +# which contain the geometric information required for the covariant form +@inline function flux(u, aux_vars, normal_direction::AbstractVector, + equations::CovariantShallowWaterEquations2D) + # Geometric variables + Gcon = metric_contravariant(aux_vars, equations) + J = area_element(aux_vars, equations) + + # Physical variables + h = waterheight(u, equations) + h_vcon = momentum_contravariant(u, equations) + + # Compute and store the velocity and gravitational terms + vcon = dot(h_vcon, normal_direction) / h + gravitational_term = 0.5f0 * equations.gravity * h^2 + + # Compute the momentum flux components in the desired orientation + momentum_flux_1 = h_vcon[1] * vcon + + dot(Gcon[1, :], normal_direction) * gravitational_term + momentum_flux_2 = h_vcon[2] * vcon + + dot(Gcon[2, :], normal_direction) * gravitational_term + + return SVector(J * dot(h_vcon, normal_direction), J * momentum_flux_1, + J * momentum_flux_2) +end + # Standard geometric and Coriolis source terms for a rotating sphere @inline function source_terms_geometric_coriolis(u, x, t, aux_vars, equations::CovariantShallowWaterEquations2D) @@ -218,7 +244,7 @@ end # Maximum wave speed along the normal direction in reference space @inline function max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr, - orientation, + orientation::Integer, equations::AbstractCovariantShallowWaterEquations2D) # Geometric variables Gcon_ll = metric_contravariant(aux_vars_ll, equations) @@ -236,6 +262,28 @@ end sqrt(Gcon_rr[orientation, orientation] * h_rr * equations.gravity)) end +# Maximum wave speed along the normal direction in reference space +@inline function max_abs_speed(u_ll, u_rr, aux_vars_ll, aux_vars_rr, + normal_direction::AbstractVector, + equations::AbstractCovariantShallowWaterEquations2D) + # Geometric variables + Gcon_ll = metric_contravariant(aux_vars_ll, equations) + Gcon_rr = metric_contravariant(aux_vars_rr, equations) + + # Physical variables + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) + h_vcon_ll = momentum_contravariant(u_ll, equations) + h_vcon_rr = momentum_contravariant(u_rr, equations) + + return max(abs(dot(normal_direction, h_vcon_ll) / h_ll) + + sqrt(dot(normal_direction, Gcon_ll * normal_direction) * h_ll * + equations.gravity), + abs(dot(normal_direction, h_vcon_rr) / h_rr) + + sqrt(dot(normal_direction, Gcon_rr * normal_direction) * h_rr * + equations.gravity)) +end + # Maximum wave speeds with respect to the covariant basis @inline function max_abs_speeds(u, aux_vars, equations::AbstractCovariantShallowWaterEquations2D) diff --git a/src/meshes/dgmulti_icosahedron_tri.jl b/src/meshes/dgmulti_icosahedron_tri.jl new file mode 100644 index 00000000..f8895bb7 --- /dev/null +++ b/src/meshes/dgmulti_icosahedron_tri.jl @@ -0,0 +1,244 @@ +function DGMultiMeshTriIcosahedron2D(dg::DGMulti{NDIMS, <:Tri}, radius; + initial_refinement = 3, + is_on_boundary = nothing) where {NDIMS} + NDIMS_AMBIENT = 3 + + vertex_coordinates = calc_node_coordinates_icosahedron_vertices(radius) + + Vxyz = ntuple(n -> vertex_coordinates[n, :], NDIMS_AMBIENT) + + EToV = zeros(Int, 20, 3) + + for i in 1:size(EToV, 1) + EToV[i, :] = icosahedron_triangle_vertices_idx_map[i] + end + + for j in 1:initial_refinement + EToV_old = EToV + Vxyz_old = ntuple(n -> copy(Vxyz[n]), NDIMS_AMBIENT) + old_to_new = Dict{Int, Int}() + edge_to_new = Dict{Tuple{Int, Int}, Int}() + EToV = zeros(Int, (size(EToV_old, 1) * 4, 3)) + + Vxyz = ntuple(n -> Vector{eltype(Vxyz[1])}(), NDIMS_AMBIENT) + + for i in 1:size(EToV_old, 1) + idx_old = EToV_old[i, :] + + for k in idx_old + if !haskey(old_to_new, k) + old_to_new[k] = length(Vxyz[1]) + 1 + for n in 1:NDIMS_AMBIENT + push!(Vxyz[n], Vxyz_old[n][k]) + end + end + end + + for k in idx_old, l in idx_old + if k < l + edge = (k, l) + if !haskey(edge_to_new, edge) + edge_to_new[edge] = length(Vxyz[1]) + 1 + vk = ntuple(n -> Vxyz_old[n][k], NDIMS_AMBIENT) + vl = ntuple(n -> Vxyz_old[n][l], NDIMS_AMBIENT) + midpoint = 0.5 .* (vk .+ vl) + midpoint = midpoint .* (EARTH_RADIUS / norm(midpoint)) # Normalize to outer radius + for n in 1:NDIMS_AMBIENT + push!(Vxyz[n], midpoint[n]) + end + end + end + end + id1 = old_to_new[idx_old[1]] + id2 = edge_to_new[Tuple(sort([idx_old[1], idx_old[2]]))] + id3 = old_to_new[idx_old[2]] + id4 = edge_to_new[Tuple(sort([idx_old[2], idx_old[3]]))] + id5 = old_to_new[idx_old[3]] + id6 = edge_to_new[Tuple(sort([idx_old[3], idx_old[1]]))] + + ids = [id1, id2, id3, id4, id5, id6] + + # Fill EToV for the 4 new triangles + for (sub_i, vertex_map) in enumerate(icosahedron_tri_vertices_idx_map) + EToV[(i - 1) * 4 + sub_i, :] = getindex.(Ref(ids), vertex_map) + end + end + end + + md = StartUpDG.MeshData(Vxyz, EToV, dg.basis) + md = spherify_meshdata!(md, dg, EToV, Vxyz, radius) + boundary_faces = StartUpDG.tag_boundary_faces(md, is_on_boundary) + return DGMultiMesh(dg, Trixi.GeometricTermsType(Trixi.Curved(), dg), md, boundary_faces) +end + +function spherify_meshdata!(md::MeshData, dg::DGMulti{NDIMS}, EToV, Vxyz, + radius) where {NDIMS} + rd = dg.basis + (; xyz, xyzq, xyzf) = md + for e in 1:size(EToV, 1) + v1, v2, v3 = ntuple(n -> SVector(ntuple(d -> Vxyz[d][EToV[e, n]], 3)), 3) + for j in 1:size(rd.rst[1], 1) + r, s = rd.rst[1][j], rd.rst[2][j] + # Bilinear mapping from reference square to physical space + x_node = local_mapping(r, s, v1, v2, v3, radius) + + for n in 1:3 + xyz[n][j, e] = x_node[n] + end + end + + for j in 1:size(rd.rstq[1], 1) + r, s = rd.rstq[1][j], rd.rstq[2][j] + # Bilinear mapping from reference square to physical space + x_node = local_mapping(r, s, v1, v2, v3, radius) + + for n in 1:3 + xyzq[n][j, e] = x_node[n] + end + end + + for j in 1:size(rd.rstf[1], 1) + r, s = rd.rstf[1][j], rd.rstf[2][j] + # Bilinear mapping from reference square to physical space + x_node = local_mapping(r, s, v1, v2, v3, radius) + + for n in 1:3 + xyzf[n][j, e] = x_node[n] + end + end + end + md = @set md.xyz = xyz + md = @set md.xyzq = xyzq + md = @set md.xyzf = xyzf + return md +end + +# We use a local numbering to obtain the triangle vertices of each triangular face +# +# Fig: Local quad vertex numbering for a triangular face (corner vertices of the triangular face in parenthesis) +# 5 (3) +# /\ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# / \ +# 6/ 4\ +# /⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺⎺\ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# / \ / \ +# 1/_________________\_/_________________3\ +# (1) 2 (2) + +# Index map for the vertices of each triangle on the triangular faces of the icosahedron (see Fig 4) +const icosahedron_tri_vertices_idx_map = ([1, 2, 6], # Tri 1 + [2, 3, 4], # Tri 2 + [4, 5, 6], # Tri 3 + [2, 4, 6]) # Tri 4 + +# Local mapping from the reference triangle to a spherical patch based on three vertex +# positions on the sphere, provided in Cartesian coordinates +@inline function local_mapping(r, s, v1, v2, v3, radius) + + # Construct a bilinear mapping based on the four corner vertices + xe = 0.5f0 * (-(r + s) * v1 + + (1 + r) * v2 + + (1 + s) * v3) + + # Project the mapped local coordinates onto the sphere + return radius * xe / norm(xe) +end + +function StartUpDG.geometric_factors(x, y, z, Dr, Ds; Filters = (I, I, I)) + xr, xs = Dr * x, Ds * x + yr, ys = Dr * y, Ds * y + zr, zs = Dr * z, Ds * z + + a1 = Dr * x, Dr * y, Dr * z + a2 = Ds * x, Dr * y, Ds * z + + a3 = @. a1[2] * a2[3] - a1[3] * a2[2], + a1[3] * a2[1] - a1[1] * a2[3], + a1[1] * a2[2] - a1[2] * a2[1] + + norm_a3 = sqrt.(a3[1] .^ 2 + a3[2] .^ 2 + a3[3] .^ 2) + + a3 = @. a3[1] / norm_a3, a3[2] / norm_a3, a3[3] / norm_a3 + + rxJ, ryJ, rzJ = @. a2[2] * a3[3] - a2[3] * a3[2], + a2[3] * a3[1] - a2[1] * a3[3], + a2[1] * a3[2] - a2[2] * a3[1] + sxJ, syJ, szJ = @. a3[2] * a1[3] - a3[3] * a1[2], + a3[3] * a1[1] - a3[1] * a1[3], + a3[1] * a1[2] - a3[2] * a1[1] + + J = @. sqrt((xr^2 + yr^2 + zr^2) * (xs^2 + ys^2 + zs^2) - + (xr * xs + yr * ys + zr * zs)^2) + return rxJ, sxJ, ryJ, syJ, rzJ, szJ, J +end + +# physical outward normals are computed via Nanson's formula: G * nhatJ, where +# G = matrix of J-scaled geometric terms. Here, Vf is a face interpolation matrix +# which maps interpolation nodes to face nodes. +function StartUpDG.compute_normals(geo::SMatrix{Dim, DimAmbient}, Vf, + nrstJ...) where {Dim, DimAmbient} + nxyzJ = ntuple(x -> zeros(size(Vf, 1), size(first(geo), 2)), DimAmbient) + for i in 1:Dim, j in 1:DimAmbient + nxyzJ[i] .+= (Vf * geo[i, j]) .* nrstJ[j] + end + Jf = sqrt.(sum(map(x -> x .^ 2, nxyzJ))) + return nxyzJ..., Jf +end + +function StartUpDG.MeshData(VX, VY, VZ, EToV, rd::RefElemData{2}; + is_periodic = (false, false, false)) + (; fv) = rd + FToF = StartUpDG.connect_mesh(EToV, fv) + Nfaces, K = size(FToF) + + #Construct global coordinates + (; V1) = rd + x, y, z = (x -> V1 * x[transpose(EToV)]).((VX, VY, VZ)) + + #Compute connectivity maps: uP = exterior value used in DG numerical fluxes + (; r, s, Vf) = rd + xf, yf, zf = (x -> Vf * x).((x, y, z)) + mapM, mapP, mapB = StartUpDG.build_node_maps(FToF, (xf, yf, zf)) + mapM = reshape(mapM, :, K) + mapP = reshape(mapP, :, K) + + #Compute geometric factors and surface normals + (; Dr, Ds) = rd + rxJ, sxJ, ryJ, syJ, rzJ, szJ, J = StartUpDG.geometric_factors(x, y, z, Dr, Ds) + rstxyzJ = SMatrix{2, 3}(rxJ, ryJ, rzJ, sxJ, syJ, szJ) + nrstJ = (rd.nrstJ..., zeros(size(rd.nrstJ[1]))) + (; Vq, wq) = rd + xq, yq, zq = (x -> Vq * x).((x, y, z)) + wJq = diagm(wq) * (Vq * J) + + nxJ, nyJ, nzJ, Jf = StartUpDG.compute_normals(rstxyzJ, rd.Vf, nrstJ...) + + periodicity = (false, false, false) + md = MeshData(StartUpDG.VertexMappedMesh(rd.element_type, tuple(VX, VY, VZ), EToV), + FToF, tuple(x, y, z), tuple(xf, yf, zf), tuple(xq, yq, zq), wJq, + mapM, mapP, mapB, + rstxyzJ, J, tuple(nxJ, nyJ, nzJ), Jf, + periodicity) + + if any(is_periodic) + # loosen the tolerance if N >> 1 + tol = length(rd.r) * 100 * eps() + md = make_periodic(md, is_periodic; tol) + end + + return md +end diff --git a/src/meshes/meshes.jl b/src/meshes/meshes.jl index 8bd9c4b0..17ca7777 100644 --- a/src/meshes/meshes.jl +++ b/src/meshes/meshes.jl @@ -1,2 +1,3 @@ include("p4est_cubed_sphere.jl") include("p4est_icosahedron_quads.jl") +include("dgmulti_icosahedron_tri.jl") diff --git a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl index 7c33a001..64d9f0c9 100644 --- a/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl +++ b/src/semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl @@ -40,7 +40,39 @@ function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{2}, performance_counter) end -# Constructor for SemidiscretizationHyperbolic for the covariant form. Requires +function Trixi.SemidiscretizationHyperbolic(mesh::DGMultiMesh, + equations::AbstractCovariantEquations, + initial_condition, + solver; + source_terms = nothing, + boundary_conditions = boundary_condition_periodic, + # `RealT` is used as real type for node locations etc. + # while `uEltype` is used as element type of solutions etc. + RealT = real(solver), uEltype = RealT, + initial_cache = NamedTuple(), + metric_terms = MetricTermsCrossProduct(), + auxiliary_field = nothing) + cache = (; + Trixi.create_cache(mesh, equations, solver, RealT, metric_terms, + auxiliary_field, uEltype)..., initial_cache...) + _boundary_conditions = Trixi.digest_boundary_conditions(boundary_conditions, mesh, + solver, + cache) + + performance_counter = Trixi.PerformanceCounter() + + SemidiscretizationHyperbolic{typeof(mesh), typeof(equations), + typeof(initial_condition), + typeof(_boundary_conditions), typeof(source_terms), + typeof(solver), typeof(cache)}(mesh, equations, + initial_condition, + _boundary_conditions, + source_terms, solver, + cache, + performance_counter) +end + +# Constructor for SemidiscretizationHyperbolic for the covariant form. Requires # compatibility between the mesh and equations (i.e. the same `NDIMS` and `NDIMS_AMBIENT`) # and sets the default metric terms to MetricTermsCovariantSphere. function Trixi.SemidiscretizationHyperbolic(mesh::P4estMesh{NDIMS, NDIMS_AMBIENT}, diff --git a/src/solvers/dgmulti/containers_manifold_covariant.jl b/src/solvers/dgmulti/containers_manifold_covariant.jl new file mode 100644 index 00000000..98f07f67 --- /dev/null +++ b/src/solvers/dgmulti/containers_manifold_covariant.jl @@ -0,0 +1,175 @@ +@muladd begin + #! format: noindnent + + function init_auxiliary_node_variables!(aux_values, mesh::DGMultiMesh, + equations::AbstractCovariantEquations{2, 3}, + dg::DGMulti{<:Any, <:Tri}, + bottom_topography) + rd = dg.basis + (; xyz) = mesh.md + md = mesh.md + n_aux = n_aux_node_vars(equations) + + # Identify the vertices corresponding to the corners of the reference element. rd.V1 is not useful, + # since the physical nodes are projected onto the sphere and thus the corner nodes would be slightly off. + # Instead, we identify the corner nodes by their reference coordinates and build a mask to access them directly + VMask = [] + for corner in [(-1.0, -1.0), (1.0, -1.0), (-1.0, 1.0)] + for j in 1:size(rd.rst[1], 1) + r, s = rd.rst[1][j], rd.rst[2][j] + if all(isapprox.((r, s), corner)) + push!(VMask, j) + end + end + end + + # Compute the radius of the sphere from the first element's fourth vertex, such that we can use it + # throughout the computation. We assume that each Wedge element's last three corner vertices lie + # on the simulated sphere. + VX, VY, VZ = map(coords -> coords[VMask, 1], xyz) + v_outer = getindex.([VX, VY, VZ], 1) + radius = norm(v_outer) + + for element in eachelement(mesh, dg) + # Compute corner vertices of the element + VX, VY, VZ = map(coords -> coords[VMask, element], xyz) + v1, v2, v3 = map(i -> getindex.([VX, VY, VZ], i), 1:3) + + aux_node = Vector{eltype(aux_values[1, 1])}(undef, n_aux) + + # Compute the auxiliary metric information at each node + for i in 1:Trixi.nnodes(dg) + # Covariant basis in the desired global coordinate system as columns of a matrix + basis_covariant = calc_basis_covariant(v1, v2, v3, + rd.rst[1][i], rd.rst[2][i], + radius, + equations.global_coordinate_system) + + aux_node[1:6] = SVector(basis_covariant) + + # Covariant metric tensor G := basis_covariant' * basis_covariant + metric_covariant = basis_covariant' * basis_covariant + + # Contravariant metric tensor inv(G) + metric_contravariant = inv(metric_covariant) + + # Contravariant basis vectors as rows of a matrix + basis_contravariant = metric_contravariant * basis_covariant' + + aux_node[7:12] = SVector(basis_contravariant) + # Area element + aux_node[13] = sqrt(det(metric_covariant)) + + # Covariant metric tensor components + aux_node[14:16] = SVector(metric_covariant[1, 1], + metric_covariant[1, 2], + metric_covariant[2, 2]) + + # Contravariant metric tensor components + aux_node[17:19] = SVector(metric_contravariant[1, 1], + metric_contravariant[1, 2], + metric_contravariant[2, 2]) + # Bottom topography + if !isnothing(bottom_topography) + x_node = map(coords -> coords[i, element], xyz) + aux_node[20] = bottom_topography(x_node) + else + aux_node[20] = zero(eltype(aux_node)) + end + aux_values[i, element] = SVector{n_aux}(aux_node) + end + # Christoffel symbols of the second kind (aux_values[21:26, :, :, element]) + calc_christoffel_symbols!(aux_values, mesh, equations, dg, element) + end + + return nothing + end + + # Analytically compute the transformation matrix A, such that G = AᵀA is the + # covariant metric tensor and a_i = A[1,i] * e_x + A[2,i] * e_y + A[3,i] * e_z denotes + # the covariant tangent basis, where e_x, e_y, and e_z are the Cartesian unit basis vectors. + @inline function calc_basis_covariant(v1, v2, v3, xi1, xi2, radius, + ::GlobalCartesianCoordinates) + + # Construct a bilinear mapping based on the three corner vertices + xe = 0.5f0 * (-(xi1 + xi2) * v1 + (1 + xi1) * v2 + + (1 + xi2) * v3) + # Derivatives of bilinear map with respect to reference coordinates xi1, xi2 + dxedxi1 = 0.5f0 * + (-v1 + v2) + dxedxi2 = 0.5f0 * + (-v1 + v3) + + # Use product/quotient rule on the projection + norm_xe = norm(xe) + dxdxi1 = radius / norm_xe * (dxedxi1 - dot(xe, dxedxi1) / norm_xe^2 * xe) + dxdxi2 = radius / norm_xe * (dxedxi2 - dot(xe, dxedxi2) / norm_xe^2 * xe) + + return SMatrix{3, 2}(dxdxi1[1], dxdxi1[2], dxdxi1[3], + dxdxi2[1], dxdxi2[2], dxdxi2[3]) + end + + function calc_christoffel_symbols!(aux_values, mesh::DGMultiMesh, + equations::AbstractCovariantEquations{2, 3}, dg, + element) + rd = dg.basis + (; Vq, Drst) = rd + Dr, Ds = Drst + + for i in 1:Trixi.nnodes(dg) + + # Numerically differentiate covariant metric components with respect to ξ¹ + dG11dxi1 = zero(eltype(aux_values[i, element])) + dG12dxi1 = zero(eltype(aux_values[i, element])) + dG22dxi1 = zero(eltype(aux_values[i, element])) + for jj in 1:nnodes(dg) + aux_node_jj = aux_values[jj, element] + Gcov_jj = metric_covariant(aux_node_jj, equations) + dG11dxi1 += Dr[i, jj] * Gcov_jj[1, 1] + dG12dxi1 += Dr[i, jj] * Gcov_jj[1, 2] + dG22dxi1 += Dr[i, jj] * Gcov_jj[2, 2] + end + + # Numerically differentiate covariant metric components with respect to ξ² + dG11dxi2 = zero(eltype(aux_values[i, element])) + dG12dxi2 = zero(eltype(aux_values[i, element])) + dG22dxi2 = zero(eltype(aux_values[i, element])) + for jj in 1:nnodes(dg) + aux_node_jj = aux_values[jj, element] + Gcov_jj = metric_covariant(aux_node_jj, equations) + dG11dxi2 += Ds[i, jj] * Gcov_jj[1, 1] + dG12dxi2 += Ds[i, jj] * Gcov_jj[1, 2] + dG22dxi2 += Ds[i, jj] * Gcov_jj[2, 2] + end + + # Compute Christoffel symbols of the first kind + christoffel_firstkind_1 = SMatrix{2, 2}(0.5f0 * dG11dxi1, + 0.5f0 * dG11dxi2, + 0.5f0 * dG11dxi2, + dG12dxi2 - 0.5f0 * dG22dxi1) + christoffel_firstkind_2 = SMatrix{2, 2}(dG12dxi1 - 0.5f0 * dG11dxi2, + 0.5f0 * dG22dxi1, + 0.5f0 * dG22dxi1, + 0.5f0 * dG22dxi2) + + # Raise indices to get Christoffel symbols of the second kind + aux_node = Vector(aux_values[i, element]) + Gcon = metric_contravariant(aux_node, equations) + aux_node[21] = Gcon[1, 1] * christoffel_firstkind_1[1, 1] + + Gcon[1, 2] * christoffel_firstkind_2[1, 1] + aux_node[22] = Gcon[1, 1] * christoffel_firstkind_1[1, 2] + + Gcon[1, 2] * christoffel_firstkind_2[1, 2] + aux_node[23] = Gcon[1, 1] * christoffel_firstkind_1[2, 2] + + Gcon[1, 2] * christoffel_firstkind_2[2, 2] + + aux_node[24] = Gcon[2, 1] * christoffel_firstkind_1[1, 1] + + Gcon[2, 2] * christoffel_firstkind_2[1, 1] + aux_node[25] = Gcon[2, 1] * christoffel_firstkind_1[1, 2] + + Gcon[2, 2] * christoffel_firstkind_2[1, 2] + aux_node[26] = Gcon[2, 1] * christoffel_firstkind_1[2, 2] + + Gcon[2, 2] * christoffel_firstkind_2[2, 2] + + aux_values[i, element] = SVector{n_aux_node_vars(equations)}(aux_node) + end + end +end # @muladd diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl new file mode 100644 index 00000000..65fa3e4f --- /dev/null +++ b/src/solvers/dgmulti/dg.jl @@ -0,0 +1,64 @@ +# Constructs cache variables including auxiliary variables for covariant equations and DGMultiMeshes +function Trixi.create_cache(mesh::DGMultiMesh{NDIMS}, equations::AbstractCovariantEquations, + dg::Trixi.DGMultiWeakForm, + RealT, metric_terms, auxiliary_field, + uEltype) where {NDIMS} + rd = dg.basis + md = mesh.md + + # volume quadrature weights, volume interpolation matrix, mass matrix, differentiation matrices + @unpack wq, Vq, M, Drst = rd + + # ∫f(u) * dv/dx_i = ∑_j (Vq*Drst[i])'*diagm(wq)*(rstxyzJ[i,j].*f(Vq*u)) + weak_differentiation_matrices = map(D -> -M \ ((Vq * D)' * Diagonal(wq)), Drst) + + nvars = nvariables(equations) + naux = n_aux_node_vars(equations) + + # storage for volume quadrature values, face quadrature values, flux values + u_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xq), dg) + u_face_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xf), dg) + flux_face_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xf), dg) + aux_values = Trixi.allocate_nested_array(uEltype, naux, size(md.x), dg) + aux_quad_values = Trixi.allocate_nested_array(uEltype, naux, size(md.xq), dg) + aux_face_values = Trixi.allocate_nested_array(uEltype, naux, size(md.xf), dg) + + if typeof(rd.approximation_type) <: + Union{SBP, Trixi.AbstractNonperiodicDerivativeOperator} + lift_scalings = rd.wf ./ rd.wq[rd.Fmask] # lift scalings for diag-norm SBP operators + else + lift_scalings = nothing + end + + # local storage for volume integral and source computations + local_values_threaded = [Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:Threads.nthreads()] + + # For curved meshes, we interpolate geometric terms from nodal points to quadrature points. + # For affine meshes, we just access one element of this interpolated data. + dxidxhatj = map(x -> rd.Vq * x, md.rstxyzJ) + + init_auxiliary_node_variables!(aux_values, mesh, equations, dg, auxiliary_field) + + # Interpolate auxiliary variables to quadrature and face points + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), aux_quad_values, aux_values) + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vf), aux_face_values, aux_values) + + # interpolate J to quadrature points for weight-adjusted DG (WADG) + invJ = inv.(area_element.(aux_quad_values, equations)) + + # for scaling by curved geometric terms (not used by affine DGMultiMesh) + flux_threaded = [[Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:NDIMS] for _ in 1:Threads.nthreads()] + rotated_flux_threaded = [Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg) + for _ in 1:Threads.nthreads()] + + cache = (; md, weak_differentiation_matrices, lift_scalings, invJ, dxidxhatj, + u_values, u_face_values, flux_face_values, + aux_values, aux_quad_values, aux_face_values, + local_values_threaded, flux_threaded, rotated_flux_threaded) + return cache +end + +include("containers_manifold_covariant.jl") +include("dg_manifold_covariant.jl") diff --git a/src/solvers/dgmulti/dg_manifold_covariant.jl b/src/solvers/dgmulti/dg_manifold_covariant.jl new file mode 100644 index 00000000..cf197ae4 --- /dev/null +++ b/src/solvers/dgmulti/dg_manifold_covariant.jl @@ -0,0 +1,112 @@ +# Compute coefficients for an initial condition that uses auxiliary variables +function Trixi.compute_coefficients!(::Nothing, u, initial_condition, t, + mesh::DGMultiMesh, + equations::AbstractCovariantEquations, + dg::DGMulti, cache) + md = mesh.md + rd = dg.basis + @unpack u_values, aux_quad_values = cache + # evaluate the initial condition at quadrature points + Trixi.@threaded for i in Trixi.each_quad_node_global(mesh, dg, cache) + x_node = SVector(getindex.(md.xyzq, i)) + aux_node = aux_quad_values[i] + u_values[i] = initial_condition(x_node, t, aux_node, equations) + end + + # multiplying by Pq computes the L2 projection + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Pq), u, u_values) +end + +# uses quadrature + projection to compute source terms. +function Trixi.calc_sources!(du, u, t, source_terms, + mesh, equations::AbstractCovariantEquations, dg::DGMulti, + cache) + rd = dg.basis + md = mesh.md + @unpack Pq = rd + @unpack u_values, aux_quad_values, local_values_threaded = cache + Trixi.@threaded for e in Trixi.eachelement(mesh, dg, cache) + source_values = local_values_threaded[Threads.threadid()] + + u_e = view(u_values, :, e) # u_values should already be computed from volume integral + aux_e = view(aux_quad_values, :, e) + + for i in Trixi.each_quad_node(mesh, dg, cache) + source_values[i] = source_terms(u_e[i], SVector(getindex.(md.xyzq, i, e)), + t, aux_e[i], equations) + end + Trixi.apply_to_each_field(Trixi.mul_by_accum!(Pq), view(du, :, e), source_values) + end + + return nothing +end + +function Trixi.calc_sources!(du, u, t, source_term::Nothing, + mesh, equations::AbstractCovariantEquations, dg::DGMulti, + cache) + nothing +end + +# version for covariant equations on DGMultiMeshes +function Trixi.calc_volume_integral!(du, u, + mesh::DGMultiMesh{NDIMS_AMBIENT, <:Trixi.NonAffine}, + have_nonconservative_terms::False, + equations::AbstractCovariantEquations{NDIMS}, + volume_integral::VolumeIntegralWeakForm, dg::DGMulti, + cache) where {NDIMS_AMBIENT, NDIMS} + rd = dg.basis + md = mesh.md + (; weak_differentiation_matrices, u_values, aux_quad_values, local_values_threaded) = cache + + # interpolate to quadrature points + Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u) + + Trixi.@threaded for e in Trixi.eachelement(mesh, dg, cache) + flux_values = local_values_threaded[Threads.threadid()] + for i in 1:NDIMS + for j in Trixi.eachindex(flux_values) + u_node = u_values[j, e] + aux_node = aux_quad_values[j, e] + area_elem = area_element(aux_node, equations) + flux_values[j] = flux(u_node, aux_node, i, equations) + end + + Trixi.apply_to_each_field(Trixi.mul_by_accum!(weak_differentiation_matrices[i]), + view(du, :, e), flux_values) + end + end +end + +function Trixi.calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm, + mesh::DGMultiMesh, + have_nonconservative_terms::False, + equations::AbstractCovariantEquations{NDIMS}, + dg::DGMulti{NDIMS_AMBIENT, <:Tri}) where {NDIMS_AMBIENT, + NDIMS} + @unpack surface_flux = surface_integral + md = mesh.md + rd = dg.basis + @unpack mapM, mapP, Jf = md + @unpack nrstJ, Nfq = rd + @unpack u_face_values, flux_face_values, aux_face_values = cache + Trixi.@threaded for face_node_index in Trixi.each_face_node_global(mesh, dg, cache) + + # inner (idM -> minus) and outer (idP -> plus) indices + idM, idP = mapM[face_node_index], mapP[face_node_index] + uM = u_face_values[idM] + uP = u_face_values[idP] + auxM = aux_face_values[idM] + auxP = aux_face_values[idP] + # Transform uP to the same coordinate system as uM + uP_global = contravariant2global(uP, auxP, equations) + uP_transformed_to_M = global2contravariant(uP_global, auxM, equations) + + # Compute ref_index from face_node_index in order to access covariant normal vector + # in reference element coordinates. + ref_index = mod(face_node_index - 1, Nfq) + 1 + normal = SVector(ntuple(k -> nrstJ[k][ref_index], NDIMS)) + + flux_face_values[idM] = surface_flux(uM, uP_transformed_to_M, auxM, auxM, normal, + equations) + end +end diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index cfcf2cf9..9a852cd7 100644 --- a/src/solvers/solvers.jl +++ b/src/solvers/solvers.jl @@ -1 +1,2 @@ include("dgsem_p4est/dg.jl") +include("dgmulti/dg.jl") diff --git a/test/test_2d_shallow_water_covariant.jl b/test/test_2d_shallow_water_covariant.jl index f5c670c8..a175a2ac 100644 --- a/test/test_2d_shallow_water_covariant.jl +++ b/test/test_2d_shallow_water_covariant.jl @@ -100,4 +100,17 @@ end @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100) end +@trixi_testset "elixir_tri_barotropic_instability" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_barotropic_instability.jl"), + l2=[41.05018196765347, 0.04598801953369521, 0.03324228006147076], + linf=[202.17195189961058, 0.2046503536574818, 0.14813768215260187]), + polydeg = 3, + initial_refinement_level = 1, + tspan = (0.0, 1.0 * SECONDS_PER_DAY) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100) +end + end # module diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 6d1c418e..d763afdf 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -128,4 +128,22 @@ end @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100) end +@trixi_testset "Spherical advection on icosahedral grid, covariant weak form, LLF surface flux" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "advection/covariant", + "elixir_tri_icosahedron.jl"), + l2=[ + 0.00037441291907931983, + 1.0851024650858329e-13, + 1.0760196451424737e-13 + ], + linf=[ + 0.010071157495758598, + 5.18004530328936e-13, + 5.133704062982442e-13 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(TrixiAtmo.Trixi.rhs!, semi, sol, 100) +end + end # module