text:
using Oceananigans
using Statistics
using KernelAbstractions: @kernel, @index, Event
using CUDA
using Test
using Printf
using Test
using Oceananigans.TimeSteppers: QuasiAdamsBashforth2TimeStepper, RungeKutta3TimeStepper, update_state!

import Oceananigans.Fields: interior

test_architectures() = CUDA.has_cuda() ? tuple(GPU()) : tuple(CPU())

function summarize_regression_test(fields, correct_fields)
    for (field_name, φ, φ_c) in zip(keys(fields), fields, correct_fields)
        Δ = φ .- φ_c

        Δ_min      = minimum(Δ)
        Δ_max      = maximum(Δ)
        Δ_mean     = mean(Δ)
        Δ_abs_mean = mean(abs, Δ)
        Δ_std      = std(Δ)

        matching    = sum(φ .≈ φ_c)
        grid_points = length(φ_c)

        @info @sprintf("Δ%s: min=%+.6e, max=%+.6e, mean=%+.6e, absmean=%+.6e, std=%+.6e (%d/%d matching grid points)",
                       field_name, Δ_min, Δ_max, Δ_mean, Δ_abs_mean, Δ_std, matching, grid_points)
    end
end

#####
##### Useful kernels
#####

@kernel function ∇²!(∇²f, grid, f)
    i, j, k = @index(Global, NTuple)
    @inbounds ∇²f[i, j, k] = ∇²ᶜᶜᶜ(i, j, k, grid, f)
end

@kernel function divergence!(grid, u, v, w, div)
    i, j, k = @index(Global, NTuple)
    @inbounds div[i, j, k] = divᶜᶜᶜ(i, j, k, grid, u, v, w)
end

function compute_∇²!(∇²ϕ, ϕ, arch, grid)
    fill_halo_regions!(ϕ, arch)
    child_arch = child_architecture(arch)
    event = launch!(child_arch, grid, :xyz, ∇²!, ∇²ϕ, grid, ϕ, dependencies=Event(device(child_arch)))
    wait(device(child_arch), event)
    fill_halo_regions!(∇²ϕ, arch)
    return nothing
end

#####
##### Useful utilities
#####

interior(a, grid) = view(a, grid.Hx+1:grid.Nx+grid.Hx,
                            grid.Hy+1:grid.Ny+grid.Hy,
                            grid.Hz+1:grid.Nz+grid.Hz)

datatuple(A) = NamedTuple{propertynames(A)}(Array(data(a)) for a in A)
datatuple(args, names) = NamedTuple{names}(a.data for a in args)

function get_model_field(field_name, model)
    if field_name ∈ (:u, :v, :w)
        return getfield(model.velocities, field_name)
    else
        return getfield(model.tracers, field_name)
    end
end

function get_output_tuple(output, iter, tuplename)
    file = jldopen(output.filepath, "r")
    output_tuple = file["timeseries/$tuplename/$iter"]
    close(file)
    return output_tuple
end

function run_script(replace_strings, script_name, script_filepath, module_suffix="")
    file_content = read(script_filepath, String)
    test_script_filepath = script_name * "_test.jl"

    for strs in replace_strings
        new_file_content = replace(file_content, strs[1] => strs[2])
        if new_file_content == file_content
            error("$(strs[1]) => $(strs[2]) replacement not found in $script_filepath. " *
                  "Make sure the script has not changed, otherwise the test might take a long time.")
            return false
        else
            file_content = new_file_content
        end
    end

    open(test_script_filepath, "w") do f
        # Wrap test script inside module to avoid polluting namespaces
        write(f, "module _Test_$script_name" * "_$module_suffix\n")
        write(f, file_content)
        write(f, "\nend # module")
    end

    try
        include(test_script_filepath)
    catch err
        @error sprint(showerror, err)

        # Print the content of the file to the test log, with line numbers, for debugging
        test_file_content = read(test_script_filepath, String)
        delineated_file_content = split(test_file_content, '\n')
        for (number, line) in enumerate(delineated_file_content)
            @printf("% 3d %s\n", number, line)
        end

        rm(test_script_filepath)
        return false
    end

    # Delete the test script (if it hasn't been deleted already)
    rm(test_script_filepath)

    return true
end

#####
##### Boundary condition utils
#####

discrete_func(i, j, grid, clock, model_fields) = - model_fields.u[i, j, grid.Nz]
parameterized_discrete_func(i, j, grid, clock, model_fields, p) = - p.μ * model_fields.u[i, j, grid.Nz]

parameterized_fun(ξ, η, t, p) = p.μ * cos(p.ω * t)
field_dependent_fun(ξ, η, t, u, v, w) = - w * sqrt(u^2 + v^2)
exploding_fun(ξ, η, t, T, S, p) = - p.μ * cosh(S - p.S0) * exp((T - p.T0) / p.λ)

# Many bc. Very many
                 integer_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, 1)
                   float_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, FT(π))
              irrational_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, π)
                   array_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, ArrayType(rand(FT, 1, 1)))
         simple_function_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, (ξ, η, t) -> exp(ξ) * cos(η) * sin(t))
  parameterized_function_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, parameterized_fun, parameters=(μ=0.1, ω=2π))
field_dependent_function_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, field_dependent_fun, field_dependencies=(:u, :v, :w))
       discrete_function_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, discrete_func, discrete_form=true)

       parameterized_discrete_function_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, parameterized_discrete_func, discrete_form=true, parameters=(μ=0.1,))
parameterized_field_dependent_function_bc(C, FT=Float64, ArrayType=Array) = BoundaryCondition(C, exploding_fun, field_dependencies=(:T, :S), parameters=(S0=35, T0=100, μ=2π, λ=FT(2)))


text:
using WaterLily
using GLMakie
using FileIO
using MeshIO

# For some reason MeshIO tries to load this binary stl as STL_ASCII DataFormat
path = abspath(joinpath(@__DIR__, "assets", "Sleigh.stl"))
sleigh = load(path)

ext: jl


text:
module Atmo

export atmo

# To think about:
# - How to handle parameters for different case? Dictionaries?

# FIXME: Be consistent with tuple assignments (either with or without parens)
#
# FIXME: Add logging

# FIXME: Add link to https://github.com/paranumal/libparanumal here and in
# advection (also update the license)

include(joinpath(@__DIR__,"vtk.jl"))
using MPI
using Canary
using Printf: @sprintf
using Utils: @hascuda

# {{{ reshape for CuArray
@hascuda function Base.reshape(A::CuArray, dims::NTuple{N, Int}) where {N}
  @assert prod(dims) == prod(size(A))
  CuArray{eltype(A), length(dims)}(dims, A.buf)
end
# }}}

# {{{ constants
# note the order of the fields below is also assumed in the code.
const _nstate = 5
const _U, _V, _W, _ρ, _E = 1:_nstate
const stateid = (U = _U, V = _V, W = _W, ρ = _ρ, E = _E)

const _nvgeo = 14
const _ξx, _ηx, _ζx, _ξy, _ηy, _ζy, _ξz, _ηz, _ζz, _MJ, _MJI,
       _x, _y, _z = 1:_nvgeo
const vgeoid = (ξx = _ξx, ηx = _ηx, ζx = _ζx,
                ξy = _ξy, ηy = _ηy, ζy = _ζy,
                ξz = _ξz, ηz = _ηz, ζz = _ζz,
                MJ = _MJ, MJI = _MJI,
                 x = _x,   y = _y,   z = _z)

const _nsgeo = 5
const _nx, _ny, _nz, _sMJ, _vMJI = 1:_nsgeo
const sgeoid = (nx = _nx, ny = _ny, nz = _nz, sMJ = _sMJ, vMJI = _vMJI)

const _γ = 14  // 10
const _p0 = 100000
const _R_gas = 28717 // 100
const _c_p = 100467 // 100
const _c_v = 7175 // 10
const _gravity = 10

# }}}

# {{{ cfl
function cfl(::Val{dim}, ::Val{N}, vgeo, Q, mpicomm) where {dim, N}
  DFloat = eltype(Q)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  Np = (N+1)^dim
  (~, ~, nelem) = size(Q)

  dt = [floatmax(DFloat)]

  if dim == 1
    @inbounds for e = 1:nelem, n = 1:Np
      ρ, U = Q[n, _ρ, e], Q[n, _U, e]
      E = Q[n, _E, e]
      P = p0 * (R_gas * E / p0)^(c_p / c_v)
      ξx = vgeo[n, _ξx, e]

      loc_dt = 2ρ / (abs(U * ξx) + ρ * sqrt(γ * P / ρ))
      dt[1] = min(dt[1], loc_dt)
    end
  end

  if dim == 2
    @inbounds for e = 1:nelem, n = 1:Np
      ρ, U, V = Q[n, _ρ, e], Q[n, _U, e], Q[n, _V, e]
      E = Q[n, _E, e]
      P = p0 * (R_gas * E / p0)^(c_p / c_v)
      ξx, ξy, ηx, ηy = vgeo[n, _ξx, e], vgeo[n, _ξy, e],
                       vgeo[n, _ηx, e], vgeo[n, _ηy, e]

      loc_dt = 2ρ / max(abs(U * ξx + V * ξy) + ρ * sqrt(γ * P / ρ),
                        abs(U * ηx + V * ηy) + ρ * sqrt(γ * P / ρ))
      dt[1] = min(dt[1], loc_dt)
    end
  end

  if dim == 3
    @inbounds for e = 1:nelem, n = 1:Np
      ρ, U, V, W = Q[n, _ρ, e], Q[n, _U, e], Q[n, _V, e], Q[n, _W, e]
      E = Q[n, _E, e]
      P = p0 * (R_gas * E / p0)^(c_p / c_v)

      ξx, ξy, ξz = vgeo[n, _ξx, e], vgeo[n, _ξy, e], vgeo[n, _ξz, e]
      ηx, ηy, ηz = vgeo[n, _ηx, e], vgeo[n, _ηy, e], vgeo[n, _ηz, e]
      ζx, ζy, ζz = vgeo[n, _ζx, e], vgeo[n, _ζy, e], vgeo[n, _ζz, e]

      loc_dt = 2ρ / max(abs(U * ξx + V * ξy + W * ξz) + ρ * sqrt(γ * P / ρ),
                        abs(U * ηx + V * ηy + W * ηz) + ρ * sqrt(γ * P / ρ),
                        abs(U * ζx + V * ζy + W * ζz) + ρ * sqrt(γ * P / ρ))
      dt[1] = min(dt[1], loc_dt)
    end
  end

  MPI.Allreduce(dt[1], MPI.MIN, mpicomm)
end
# }}}

# {{{ compute geometry
function computegeometry(::Val{dim}, mesh, D, ξ, ω, meshwarp, vmapM) where dim
  # Compute metric terms
  Nq = size(D, 1)
  DFloat = eltype(D)

  (nface, nelem) = size(mesh.elemtoelem)

  crd = creategrid(Val(dim), mesh.elemtocoord, ξ)

  vgeo = zeros(DFloat, Nq^dim, _nvgeo, nelem)
  sgeo = zeros(DFloat, _nsgeo, Nq^(dim-1), nface, nelem)

  (ξx, ηx, ζx, ξy, ηy, ζy, ξz, ηz, ζz, MJ, MJI, x, y, z) =
      ntuple(j->(@view vgeo[:, j, :]), _nvgeo)
  J = similar(x)
  (nx, ny, nz, sMJ, vMJI) = ntuple(j->(@view sgeo[ j, :, :, :]), _nsgeo)
  sJ = similar(sMJ)

  X = ntuple(j->(@view vgeo[:, _x+j-1, :]), dim)
  creategrid!(X..., mesh.elemtocoord, ξ)

  @inbounds for j = 1:length(x)
    (x[j], y[j], z[j]) = meshwarp(x[j], y[j], z[j])
  end

  # Compute the metric terms
  if dim == 1
    computemetric!(x, J, ξx, sJ, nx, D)
  elseif dim == 2
    computemetric!(x, y, J, ξx, ηx, ξy, ηy, sJ, nx, ny, D)
  elseif dim == 3
    computemetric!(x, y, z, J, ξx, ηx, ζx, ξy, ηy, ζy, ξz, ηz, ζz, sJ,
                   nx, ny, nz, D)
  end

  M = kron(1, ntuple(j->ω, dim)...)
  MJ .= M .* J
  MJI .= 1 ./ MJ
  vMJI .= MJI[vmapM]

  sM = dim > 1 ? kron(1, ntuple(j->ω, dim-1)...) : one(DFloat)
  sMJ .= sM .* sJ

  (vgeo, sgeo)
end
# }}}

# {{{ CPU Kernels
# {{{ 2-D
# Volume RHS for 2-D
function volumerhs!(::Val{2}, ::Val{N}, rhs::Array, Q, vgeo, D, elems) where N
  DFloat = eltype(Q)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  Nq = N + 1

  nelem = size(Q)[end]

  Q = reshape(Q, Nq, Nq, _nstate, nelem)
  rhs = reshape(rhs, Nq, Nq, _nstate, nelem)
  vgeo = reshape(vgeo, Nq, Nq, _nvgeo, nelem)

  s_F = Array{DFloat}(undef, Nq, Nq, _nstate)
  s_G = Array{DFloat}(undef, Nq, Nq, _nstate)

  @inbounds for e in elems
    for j = 1:Nq, i = 1:Nq
      MJ = vgeo[i, j, _MJ, e]
      ξx, ξy = vgeo[i,j,_ξx,e], vgeo[i,j,_ξy,e]
      ηx, ηy = vgeo[i,j,_ηx,e], vgeo[i,j,_ηy,e]

      U, V = Q[i, j, _U, e], Q[i, j, _V, e]
      ρ, E = Q[i, j, _ρ, e], Q[i, j, _E, e]

      P = p0 * (R_gas * E / p0)^(c_p / c_v)

      ρinv = 1 / ρ
      fluxρ_x = U
      fluxU_x = ρinv * U * U + P
      fluxV_x = ρinv * V * U
      fluxE_x = ρinv * U * E

      fluxρ_y = V
      fluxU_y = ρinv * U * V
      fluxV_y = ρinv * V * V + P
      fluxE_y = ρinv * V * E

      s_F[i, j, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y)
      s_F[i, j, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y)
      s_F[i, j, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y)
      s_F[i, j, _W] = 0
      s_F[i, j, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y)

      s_G[i, j, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y)
      s_G[i, j, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y)
      s_G[i, j, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y)
      s_G[i, j, _W] = 0
      s_G[i, j, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y)

      # buoyancy term
      rhs[i, j, _V, e] -= MJ * ρ * gravity
    end

    # loop of ξ-grid lines
    for s = 1:_nstate, j = 1:Nq, i = 1:Nq, n = 1:Nq
      rhs[i, j, s, e] += D[n, i] * s_F[n, j, s]
    end
    # loop of η-grid lines
    for s = 1:_nstate, j = 1:Nq, i = 1:Nq, n = 1:Nq
      rhs[i, j, s, e] += D[n, j] * s_G[i, n, s]
    end
  end
end

# Face RHS for 2-D
function facerhs!(::Val{2}, ::Val{N}, rhs::Array, Q, sgeo, elems, vmapM,
                  vmapP, elemtobndy) where N
  DFloat = eltype(Q)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  Np = (N+1)^2
  Nfp = N+1
  nface = 4

  @inbounds for e in elems
    for f = 1:nface
      for n = 1:Nfp
        nxM, nyM, sMJ = sgeo[_nx, n, f, e], sgeo[_ny, n, f, e], sgeo[_sMJ, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        ρM = Q[vidM, _ρ, eM]
        UM = Q[vidM, _U, eM]
        VM = Q[vidM, _V, eM]
        EM = Q[vidM, _E, eM]

        bc = elemtobndy[f, e]
        PM = p0 * (R_gas * EM / p0)^(c_p / c_v)
        if bc == 0
          ρP = Q[vidP, _ρ, eP]
          UP = Q[vidP, _U, eP]
          VP = Q[vidP, _V, eP]
          EP = Q[vidP, _E, eP]
          PP = p0 * (R_gas * EP / p0)^(c_p / c_v)
        elseif bc == 1
          UnM = nxM * UM + nyM * VM
          UP = UM - 2 * UnM * nxM
          VP = VM - 2 * UnM * nyM
          ρP = ρM
          EP = EM
          PP = PM
        else
          error("Invalid boundary conditions $bc on face $f of element $e")
        end

        ρMinv = 1 / ρM
        fluxρM_x = UM
        fluxUM_x = ρMinv * UM * UM + PM
        fluxVM_x = ρMinv * VM * UM
        fluxEM_x = ρMinv * UM * EM

        fluxρM_y = VM
        fluxUM_y = ρMinv * UM * VM
        fluxVM_y = ρMinv * VM * VM + PM
        fluxEM_y = ρMinv * VM * EM

        ρPinv = 1 / ρP
        fluxρP_x = UP
        fluxUP_x = ρPinv * UP * UP + PP
        fluxVP_x = ρPinv * VP * UP
        fluxEP_x = ρPinv * UP * EP

        fluxρP_y = VP
        fluxUP_y = ρPinv * UP * VP
        fluxVP_y = ρPinv * VP * VP + PP
        fluxEP_y = ρPinv * VP * EP

        λM = ρMinv * abs(nxM * UM + nyM * VM) + sqrt(ρMinv * γ * PM)
        λP = ρPinv * abs(nxM * UP + nyM * VP) + sqrt(ρPinv * γ * PP)
        λ  =  max(λM, λP)

        #Compute Numerical Flux and Update
        fluxρS = (nxM * (fluxρM_x + fluxρP_x) + nyM * (fluxρM_y + fluxρP_y) +
                  - λ * (ρP - ρM)) / 2
        fluxUS = (nxM * (fluxUM_x + fluxUP_x) + nyM * (fluxUM_y + fluxUP_y) +
                  - λ * (UP - UM)) / 2
        fluxVS = (nxM * (fluxVM_x + fluxVP_x) + nyM * (fluxVM_y + fluxVP_y) +
                  - λ * (VP - VM)) / 2
        fluxES = (nxM * (fluxEM_x + fluxEP_x) + nyM * (fluxEM_y + fluxEP_y) +
                  - λ * (EP - EM)) / 2

        #Update RHS
        rhs[vidM, _ρ, eM] -= sMJ * fluxρS
        rhs[vidM, _U, eM] -= sMJ * fluxUS
        rhs[vidM, _V, eM] -= sMJ * fluxVS
        rhs[vidM, _E, eM] -= sMJ * fluxES
      end
    end
  end
end
# }}}

# {{{ 3-D
# Volume RHS for 3-D
function volumerhs!(::Val{3}, ::Val{N}, rhs::Array, Q, vgeo, D, elems) where N
  DFloat = eltype(Q)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  Nq = N + 1

  nelem = size(Q)[end]

  Q = reshape(Q, Nq, Nq, Nq, _nstate, nelem)
  rhs = reshape(rhs, Nq, Nq, Nq, _nstate, nelem)
  vgeo = reshape(vgeo, Nq, Nq, Nq, _nvgeo, nelem)

  s_F = Array{DFloat}(undef, Nq, Nq, Nq, _nstate)
  s_G = Array{DFloat}(undef, Nq, Nq, Nq, _nstate)
  s_H = Array{DFloat}(undef, Nq, Nq, Nq, _nstate)

  @inbounds for e in elems
    for k = 1:Nq, j = 1:Nq, i = 1:Nq
      MJ = vgeo[i, j, k, _MJ, e]
      ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
      ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
      ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]

      U, V, W = Q[i, j, k, _U, e], Q[i, j, k, _V, e], Q[i, j, k, _W, e]
      ρ, E = Q[i, j, k, _ρ, e], Q[i, j, k, _E, e]

      P = p0 * (R_gas * E / p0)^(c_p / c_v)

      ρinv = 1 / ρ
      fluxρ_x = U
      fluxU_x = ρinv * U * U + P
      fluxV_x = ρinv * V * U
      fluxW_x = ρinv * W * U
      fluxE_x = E * ρinv * U

      fluxρ_y = V
      fluxU_y = ρinv * U * V
      fluxV_y = ρinv * V * V + P
      fluxW_y = ρinv * W * V
      fluxE_y = E * ρinv * V

      fluxρ_z = W
      fluxU_z = ρinv * U * W
      fluxV_z = ρinv * V * W
      fluxW_z = ρinv * W * W + P
      fluxE_z = E * ρinv * W

      s_F[i, j, k, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y + ξz * fluxρ_z)
      s_F[i, j, k, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y + ξz * fluxU_z)
      s_F[i, j, k, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y + ξz * fluxV_z)
      s_F[i, j, k, _W] = MJ * (ξx * fluxW_x + ξy * fluxW_y + ξz * fluxW_z)
      s_F[i, j, k, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y + ξz * fluxE_z)

      s_G[i, j, k, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y + ηz * fluxρ_z)
      s_G[i, j, k, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y + ηz * fluxU_z)
      s_G[i, j, k, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y + ηz * fluxV_z)
      s_G[i, j, k, _W] = MJ * (ηx * fluxW_x + ηy * fluxW_y + ηz * fluxW_z)
      s_G[i, j, k, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y + ηz * fluxE_z)

      s_H[i, j, k, _ρ] = MJ * (ζx * fluxρ_x + ζy * fluxρ_y + ζz * fluxρ_z)
      s_H[i, j, k, _U] = MJ * (ζx * fluxU_x + ζy * fluxU_y + ζz * fluxU_z)
      s_H[i, j, k, _V] = MJ * (ζx * fluxV_x + ζy * fluxV_y + ζz * fluxV_z)
      s_H[i, j, k, _W] = MJ * (ζx * fluxW_x + ζy * fluxW_y + ζz * fluxW_z)
      s_H[i, j, k, _E] = MJ * (ζx * fluxE_x + ζy * fluxE_y + ζz * fluxE_z)

      # buoyancy term
      rhs[i, j, k, _W, e] -= MJ * ρ * gravity
    end

    # loop of ξ-grid lines
    for s = 1:_nstate, k = 1:Nq, j = 1:Nq, i = 1:Nq, n = 1:Nq
      rhs[i, j, k, s, e] += D[n, i] * s_F[n, j, k, s]
    end
    # loop of η-grid lines
    for s = 1:_nstate, k = 1:Nq, j = 1:Nq, i = 1:Nq, n = 1:Nq
      rhs[i, j, k, s, e] += D[n, j] * s_G[i, n, k, s]
    end
    # loop of ζ-grid lines
    for s = 1:_nstate, k = 1:Nq, j = 1:Nq, i = 1:Nq, n = 1:Nq
      rhs[i, j, k, s, e] += D[n, k] * s_H[i, j, n, s]
    end
  end
end

# Face RHS for 3-D
function facerhs!(::Val{3}, ::Val{N}, rhs::Array, Q, sgeo, elems, vmapM,
                  vmapP, elemtobndy) where N
  DFloat = eltype(Q)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  Np = (N+1)^3
  Nfp = (N+1)^2
  nface = 6

  @inbounds for e in elems
    for f = 1:nface
      for n = 1:Nfp
        (nxM, nyM, nzM, sMJ, ~) = sgeo[:, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        ρM = Q[vidM, _ρ, eM]
        UM = Q[vidM, _U, eM]
        VM = Q[vidM, _V, eM]
        WM = Q[vidM, _W, eM]
        EM = Q[vidM, _E, eM]

        bc = elemtobndy[f, e]
        PM = p0 * (R_gas * EM / p0)^(c_p / c_v)
        if bc == 0
          ρP = Q[vidP, _ρ, eP]
          UP = Q[vidP, _U, eP]
          VP = Q[vidP, _V, eP]
          WP = Q[vidP, _W, eP]
          EP = Q[vidP, _E, eP]
          PP = p0 * (R_gas * EP / p0)^(c_p / c_v)
        elseif bc == 1
          UnM = nxM * UM + nyM * VM + nzM * WM
          UP = UM - 2 * UnM * nxM
          VP = VM - 2 * UnM * nyM
          WP = WM - 2 * UnM * nzM
          ρP = ρM
          EP = EM
          PP = PM
        else
          error("Invalid boundary conditions $bc on face $f of element $e")
        end

        ρMinv = 1 / ρM
        fluxρM_x = UM
        fluxUM_x = ρMinv * UM * UM + PM
        fluxVM_x = ρMinv * VM * UM
        fluxWM_x = ρMinv * WM * UM
        fluxEM_x = ρMinv * UM * EM

        fluxρM_y = VM
        fluxUM_y = ρMinv * UM * VM
        fluxVM_y = ρMinv * VM * VM + PM
        fluxWM_y = ρMinv * WM * VM
        fluxEM_y = ρMinv * VM * EM

        fluxρM_z = WM
        fluxUM_z = ρMinv * UM * WM
        fluxVM_z = ρMinv * VM * WM
        fluxWM_z = ρMinv * WM * WM + PM
        fluxEM_z = ρMinv * WM * EM

        ρPinv = 1 / ρP
        fluxρP_x = UP
        fluxUP_x = ρPinv * UP * UP + PP
        fluxVP_x = ρPinv * VP * UP
        fluxWP_x = ρPinv * WP * UP
        fluxEP_x = ρPinv * UP * EP

        fluxρP_y = VP
        fluxUP_y = ρPinv * UP * VP
        fluxVP_y = ρPinv * VP * VP + PP
        fluxWP_y = ρPinv * WP * VP
        fluxEP_y = ρPinv * VP * EP

        fluxρP_z = WP
        fluxUP_z = ρPinv * UP * WP
        fluxVP_z = ρPinv * VP * WP
        fluxWP_z = ρPinv * WP * WP + PP
        fluxEP_z = ρPinv * WP * EP

        λM = ρMinv * abs(nxM * UM + nyM * VM + nzM * WM) + sqrt(ρMinv * γ * PM)
        λP = ρPinv * abs(nxM * UP + nyM * VP + nzM * WP) + sqrt(ρPinv * γ * PP)
        λ  =  max(λM, λP)

        #Compute Numerical Flux and Update
        fluxρS = (nxM * (fluxρM_x + fluxρP_x) + nyM * (fluxρM_y + fluxρP_y) +
                  nzM * (fluxρM_z + fluxρP_z) - λ * (ρP - ρM)) / 2
        fluxUS = (nxM * (fluxUM_x + fluxUP_x) + nyM * (fluxUM_y + fluxUP_y) +
                  nzM * (fluxUM_z + fluxUP_z) - λ * (UP - UM)) / 2
        fluxVS = (nxM * (fluxVM_x + fluxVP_x) + nyM * (fluxVM_y + fluxVP_y) +
                  nzM * (fluxVM_z + fluxVP_z) - λ * (VP - VM)) / 2
        fluxWS = (nxM * (fluxWM_x + fluxWP_x) + nyM * (fluxWM_y + fluxWP_y) +
                  nzM * (fluxWM_z + fluxWP_z) - λ * (WP - WM)) / 2
        fluxES = (nxM * (fluxEM_x + fluxEP_x) + nyM * (fluxEM_y + fluxEP_y) +
                  nzM * (fluxEM_z + fluxEP_z) - λ * (EP - EM)) / 2

        #Update RHS
        rhs[vidM, _ρ, eM] -= sMJ * fluxρS
        rhs[vidM, _U, eM] -= sMJ * fluxUS
        rhs[vidM, _V, eM] -= sMJ * fluxVS
        rhs[vidM, _W, eM] -= sMJ * fluxWS
        rhs[vidM, _E, eM] -= sMJ * fluxES
      end
    end
  end
end
# }}}

# {{{ Update solution (for all dimensions)
function updatesolution!(::Val{dim}, ::Val{N}, rhs::Array, Q, vgeo, elems, rka,
                         rkb, dt) where {dim, N}
  @inbounds for e = elems, s = 1:_nstate, i = 1:(N+1)^dim
    Q[i, s, e] += rkb * dt * rhs[i, s, e] * vgeo[i, _MJI, e]
    rhs[i, s, e] *= rka
  end
end

# }}}
# }}}

# {{{ improved GPU kernles

# {{{ Volume RHS for 2-D
@hascuda function knl_volumerhs!(::Val{2}, ::Val{N}, rhs, Q, vgeo, D, nelem) where N
  DFloat = eltype(D)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  Nq = N + 1

  (i, j, k) = threadIdx()
  e = blockIdx().x

  s_D = @cuStaticSharedMem(eltype(D), (Nq, Nq))
  s_F = @cuStaticSharedMem(eltype(Q), (Nq, Nq, _nstate))
  s_G = @cuStaticSharedMem(eltype(Q), (Nq, Nq, _nstate))

  rhsU = rhsV = rhsρ = rhsE = zero(eltype(rhs))
  if i <= Nq && j <= Nq && k == 1 && e <= nelem
    # Load derivative into shared memory
    if k == 1
      s_D[i, j] = D[i, j]
    end

    # Load values will need into registers
    MJ = vgeo[i, j, _MJ, e]
    ξx, ξy = vgeo[i, j, _ξx, e], vgeo[i, j, _ξy, e]
    ηx, ηy = vgeo[i, j, _ηx, e], vgeo[i, j, _ηy, e]
    U, V = Q[i, j, _U, e], Q[i, j, _V, e]
    ρ, E = Q[i, j, _ρ, e], Q[i, j, _E, e]
    rhsU, rhsV = rhs[i, j, _U, e], rhs[i, j, _V, e]
    rhsρ, rhsE = rhs[i, j, _ρ, e], rhs[i, j, _E, e]

    P = p0 * CUDAnative.pow(R_gas * E / p0, c_p / c_v)

    ρinv = 1 / ρ
    fluxρ_x = U
    fluxU_x = ρinv * U * U + P
    fluxV_x = ρinv * V * U
    fluxE_x = ρinv * U * E

    fluxρ_y = V
    fluxU_y = ρinv * U * V
    fluxV_y = ρinv * V * V + P
    fluxE_y = ρinv * V * E

    s_F[i, j, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y)
    s_F[i, j, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y)
    s_F[i, j, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y)
    s_F[i, j, _W] = 0
    s_F[i, j, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y)

    s_G[i, j, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y)
    s_G[i, j, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y)
    s_G[i, j, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y)
    s_G[i, j, _W] = 0
    s_G[i, j, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y)

    # buoyancy term
    rhsV -= MJ * ρ * gravity
  end

  sync_threads()

  @inbounds if i <= Nq && j <= Nq && k == 1 && e <= nelem
    # loop of ξ-grid lines
    for n = 1:Nq
      Dni = s_D[n, i]
      Dnj = s_D[n, j]

      rhsρ += Dni * s_F[n, j, _ρ]
      rhsρ += Dnj * s_G[i, n, _ρ]

      rhsU += Dni * s_F[n, j, _U]
      rhsU += Dnj * s_G[i, n, _U]

      rhsV += Dni * s_F[n, j, _V]
      rhsV += Dnj * s_G[i, n, _V]

      rhsE += Dni * s_F[n, j, _E]
      rhsE += Dnj * s_G[i, n, _E]
    end

    rhs[i, j, _U, e] = rhsU
    rhs[i, j, _V, e] = rhsV
    rhs[i, j, _ρ, e] = rhsρ
    rhs[i, j, _E, e] = rhsE
  end
  nothing
end
# }}}

# {{{ Volume RHS for 3-D
@hascuda function knl_volumerhs!(::Val{3}, ::Val{N}, rhs, Q, vgeo, D, nelem) where N
  DFloat = eltype(D)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  Nq = N + 1

  (i, j, k) = threadIdx()
  e = blockIdx().x

  s_D = @cuStaticSharedMem(eltype(D), (Nq, Nq))
  s_F = @cuStaticSharedMem(eltype(Q), (Nq, Nq, Nq, _nstate))
  s_G = @cuStaticSharedMem(eltype(Q), (Nq, Nq, Nq, _nstate))
  s_H = @cuStaticSharedMem(eltype(Q), (Nq, Nq, Nq, _nstate))

  rhsU = rhsV = rhsW = rhsρ = rhsE = zero(eltype(rhs))
  @inbounds if i <= Nq && j <= Nq && k <= Nq && e <= nelem
    # Load derivative into shared memory
    if k == 1
      s_D[i, j] = D[i, j]
    end

    # Load values will need into registers
    MJ = vgeo[i, j, k, _MJ, e]
    ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
    ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
    ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]

    U, V, W = Q[i, j, k, _U, e], Q[i, j, k, _V, e], Q[i, j, k, _W, e]
    ρ, E = Q[i, j, k, _ρ, e], Q[i, j, k, _E, e]

    P = p0 * CUDAnative.pow(R_gas * E / p0, c_p / c_v)

    ρinv = 1 / ρ
    fluxρ_x = U
    fluxU_x = ρinv * U * U + P
    fluxV_x = ρinv * V * U
    fluxW_x = ρinv * W * U
    fluxE_x = E * ρinv * U

    fluxρ_y = V
    fluxU_y = ρinv * U * V
    fluxV_y = ρinv * V * V + P
    fluxW_y = ρinv * W * V
    fluxE_y = E * ρinv * V

    fluxρ_z = W
    fluxU_z = ρinv * U * W
    fluxV_z = ρinv * V * W
    fluxW_z = ρinv * W * W + P
    fluxE_z = E * ρinv * W

    s_F[i, j, k, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y + ξz * fluxρ_z)
    s_F[i, j, k, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y + ξz * fluxU_z)
    s_F[i, j, k, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y + ξz * fluxV_z)
    s_F[i, j, k, _W] = MJ * (ξx * fluxW_x + ξy * fluxW_y + ξz * fluxW_z)
    s_F[i, j, k, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y + ξz * fluxE_z)

    s_G[i, j, k, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y + ηz * fluxρ_z)
    s_G[i, j, k, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y + ηz * fluxU_z)
    s_G[i, j, k, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y + ηz * fluxV_z)
    s_G[i, j, k, _W] = MJ * (ηx * fluxW_x + ηy * fluxW_y + ηz * fluxW_z)
    s_G[i, j, k, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y + ηz * fluxE_z)

    s_H[i, j, k, _ρ] = MJ * (ζx * fluxρ_x + ζy * fluxρ_y + ζz * fluxρ_z)
    s_H[i, j, k, _U] = MJ * (ζx * fluxU_x + ζy * fluxU_y + ζz * fluxU_z)
    s_H[i, j, k, _V] = MJ * (ζx * fluxV_x + ζy * fluxV_y + ζz * fluxV_z)
    s_H[i, j, k, _W] = MJ * (ζx * fluxW_x + ζy * fluxW_y + ζz * fluxW_z)
    s_H[i, j, k, _E] = MJ * (ζx * fluxE_x + ζy * fluxE_y + ζz * fluxE_z)

    rhsU, rhsV, rhsW = (rhs[i, j, k, _U, e],
                        rhs[i, j, k, _V, e],
                        rhs[i, j, k, _W, e])
    rhsρ, rhsE = rhs[i, j, k, _ρ, e], rhs[i, j, k, _E, e]

    # buoyancy term
    rhsW -= MJ * ρ * gravity
  end

  sync_threads()

  @inbounds if i <= Nq && j <= Nq && k <= Nq && e <= nelem
    # loop of ξ-grid lines
    for n = 1:Nq
      Dni = s_D[n, i]
      Dnj = s_D[n, j]
      Dnk = s_D[n, k]

      rhsρ += Dni * s_F[n, j, k, _ρ]
      rhsρ += Dnj * s_G[i, n, k, _ρ]
      rhsρ += Dnk * s_H[i, j, n, _ρ]

      rhsU += Dni * s_F[n, j, k, _U]
      rhsU += Dnj * s_G[i, n, k, _U]
      rhsU += Dnk * s_H[i, j, n, _U]

      rhsV += Dni * s_F[n, j, k, _V]
      rhsV += Dnj * s_G[i, n, k, _V]
      rhsV += Dnk * s_H[i, j, n, _V]

      rhsW += Dni * s_F[n, j, k, _W]
      rhsW += Dnj * s_G[i, n, k, _W]
      rhsW += Dnk * s_H[i, j, n, _W]

      rhsE += Dni * s_F[n, j, k, _E]
      rhsE += Dnj * s_G[i, n, k, _E]
      rhsE += Dnk * s_H[i, j, n, _E]
    end

    rhs[i, j, k, _U, e] = rhsU
    rhs[i, j, k, _V, e] = rhsV
    rhs[i, j, k, _W, e] = rhsW
    rhs[i, j, k, _ρ, e] = rhsρ
    rhs[i, j, k, _E, e] = rhsE
  end
  nothing
end
# }}}

# {{{ Face RHS (all dimensions)
@hascuda function knl_facerhs!(::Val{dim}, ::Val{N}, rhs, Q, sgeo, nelem, vmapM,
                               vmapP, elemtobndy) where {dim, N}
  DFloat = eltype(Q)
  γ::DFloat       = _γ
  p0::DFloat      = _p0
  R_gas::DFloat   = _R_gas
  c_p::DFloat     = _c_p
  c_v::DFloat     = _c_v
  gravity::DFloat = _gravity

  if dim == 1
    Np = (N+1)
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    nface = 6
  end

  (i, j, k) = threadIdx()
  e = blockIdx().x

  Nq = N+1
  half = convert(eltype(Q), 0.5)

  @inbounds if i <= Nq && j <= Nq && k == 1 && e <= nelem
    n = i + (j-1) * Nq
    for lf = 1:2:nface
      for f = lf:lf+1
        (nxM, nyM) = (sgeo[_nx, n, f, e], sgeo[_ny, n, f, e])
        (nzM, sMJ) = (sgeo[_nz, n, f, e], sgeo[_sMJ, n, f, e])

        (idM, idP) = (vmapM[n, f, e], vmapP[n, f, e])

        (eM, eP) = (e, ((idP - 1) ÷ Np) + 1)
        (vidM, vidP) = (((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1)

        ρM = Q[vidM, _ρ, eM]
        UM = Q[vidM, _U, eM]
        VM = Q[vidM, _V, eM]
        WM = Q[vidM, _W, eM]
        EM = Q[vidM, _E, eM]

        bc = elemtobndy[f, e]
        PM = p0 * CUDAnative.pow(R_gas * EM / p0, c_p / c_v)
        ρP = UP = VP = WP = EP = PP = zero(eltype(Q))
        if bc == 0
          ρP = Q[vidP, _ρ, eP]
          UP = Q[vidP, _U, eP]
          VP = Q[vidP, _V, eP]
          WP = Q[vidP, _W, eP]
          EP = Q[vidP, _E, eP]
          PP = p0 * CUDAnative.pow(R_gas * EP / p0, c_p / c_v)
        elseif bc == 1
          UnM = nxM * UM + nyM * VM + nzM * WM
          UP = UM - 2 * UnM * nxM
          VP = VM - 2 * UnM * nyM
          WP = WM - 2 * UnM * nzM
          ρP = ρM
          EP = EM
          PP = PM
        end

        ρMinv = 1 / ρM
        fluxρM_x = UM
        fluxUM_x = ρMinv * UM * UM + PM
        fluxVM_x = ρMinv * VM * UM
        fluxWM_x = ρMinv * WM * UM
        fluxEM_x = ρMinv * UM * EM

        fluxρM_y = VM
        fluxUM_y = ρMinv * UM * VM
        fluxVM_y = ρMinv * VM * VM + PM
        fluxWM_y = ρMinv * WM * VM
        fluxEM_y = ρMinv * VM * EM

        fluxρM_z = WM
        fluxUM_z = ρMinv * UM * WM
        fluxVM_z = ρMinv * VM * WM
        fluxWM_z = ρMinv * WM * WM + PM
        fluxEM_z = ρMinv * WM * EM

        ρPinv = 1 / ρP
        fluxρP_x = UP
        fluxUP_x = ρPinv * UP * UP + PP
        fluxVP_x = ρPinv * VP * UP
        fluxWP_x = ρPinv * WP * UP
        fluxEP_x = ρPinv * UP * EP

        fluxρP_y = VP
        fluxUP_y = ρPinv * UP * VP
        fluxVP_y = ρPinv * VP * VP + PP
        fluxWP_y = ρPinv * WP * VP
        fluxEP_y = ρPinv * VP * EP

        fluxρP_z = WP
        fluxUP_z = ρPinv * UP * WP
        fluxVP_z = ρPinv * VP * WP
        fluxWP_z = ρPinv * WP * WP + PP
        fluxEP_z = ρPinv * WP * EP

        λM = ρMinv * abs(nxM * UM + nyM * VM + nzM * WM) + CUDAnative.sqrt(ρMinv * γ * PM)
        λP = ρPinv * abs(nxM * UP + nyM * VP + nzM * WP) + CUDAnative.sqrt(ρPinv * γ * PP)
        λ  =  max(λM, λP)

        #Compute Numerical Flux and Update
        fluxρS = (nxM * (fluxρM_x + fluxρP_x) + nyM * (fluxρM_y + fluxρP_y) +
                  nzM * (fluxρM_z + fluxρP_z) - λ * (ρP - ρM)) / 2
        fluxUS = (nxM * (fluxUM_x + fluxUP_x) + nyM * (fluxUM_y + fluxUP_y) +
                  nzM * (fluxUM_z + fluxUP_z) - λ * (UP - UM)) / 2
        fluxVS = (nxM * (fluxVM_x + fluxVP_x) + nyM * (fluxVM_y + fluxVP_y) +
                  nzM * (fluxVM_z + fluxVP_z) - λ * (VP - VM)) / 2
        fluxWS = (nxM * (fluxWM_x + fluxWP_x) + nyM * (fluxWM_y + fluxWP_y) +
                  nzM * (fluxWM_z + fluxWP_z) - λ * (WP - WM)) / 2
        fluxES = (nxM * (fluxEM_x + fluxEP_x) + nyM * (fluxEM_y + fluxEP_y) +
                  nzM * (fluxEM_z + fluxEP_z) - λ * (EP - EM)) / 2

        #Update RHS
        rhs[vidM, _ρ, eM] -= sMJ * fluxρS
        rhs[vidM, _U, eM] -= sMJ * fluxUS
        rhs[vidM, _V, eM] -= sMJ * fluxVS
        rhs[vidM, _W, eM] -= sMJ * fluxWS
        rhs[vidM, _E, eM] -= sMJ * fluxES
      end
      sync_threads()
    end
  end
  nothing
end
# }}}

# {{{ Update solution (for all dimensions)
@hascuda function knl_updatesolution!(::Val{dim}, ::Val{N}, rhs, Q, vgeo, nelem, rka,
                             rkb, dt) where {dim, N}
  (i, j, k) = threadIdx()
  e = blockIdx().x

  Nq = N+1
  @inbounds if i <= Nq && j <= Nq && k <= Nq && e <= nelem
    n = i + (j-1) * Nq + (k-1) * Nq * Nq
    MJI = vgeo[n, _MJI, e]
    for s = 1:_nstate
      Q[n, s, e] += rkb * dt * rhs[n, s, e] * MJI
      rhs[n, s, e] *= rka
    end
  end
  nothing
end
# }}}

# }}}

# {{{ Fill sendQ on device with Q (for all dimensions)
@hascuda function knl_fillsendQ!(::Val{dim}, ::Val{N}, sendQ, Q,
                        sendelems) where {N, dim}
  Nq = N + 1
  (i, j, k) = threadIdx()
  e = blockIdx().x

  @inbounds if i <= Nq && j <= Nq && k <= Nq && e <= length(sendelems)
    n = i + (j-1) * Nq + (k-1) * Nq * Nq
    re = sendelems[e]
    for s = 1:_nstate
      sendQ[n, s, e] = Q[n, s, re]
    end
  end
  nothing
end
# }}}

# {{{ Fill Q on device with recvQ (for all dimensions)
@hascuda function knl_transferrecvQ!(::Val{dim}, ::Val{N}, Q, recvQ, nelem,
                            nrealelem) where {N, dim}
  Nq = N + 1
  (i, j, k) = threadIdx()
  e = blockIdx().x

  @inbounds if i <= Nq && j <= Nq && k <= Nq && e <= nelem
    n = i + (j-1) * Nq + (k-1) * Nq * Nq
    for s = 1:_nstate
      Q[n, s, nrealelem + e] = recvQ[n, s, e]
    end
  end
  nothing
end
# }}}

# {{{ MPI Buffer handling
function fillsendQ!(::Val{dim}, ::Val{N}, sendQ, d_sendQ::Array, Q,
                    sendelems) where {dim, N}
  sendQ[:, :, :] .= Q[:, :, sendelems]
end

@hascuda function fillsendQ!(::Val{dim}, ::Val{N}, sendQ, d_sendQ::CuArray,
                             d_QL, d_sendelems) where {dim, N}
  nsendelem = length(d_sendelems)
  if nsendelem > 0
    @cuda(threads=ntuple(j->N+1, dim), blocks=nsendelem,
          knl_fillsendQ!(Val(dim), Val(N), d_sendQ, d_QL, d_sendelems))
    sendQ .= d_sendQ
  end
end

@hascuda function transferrecvQ!(::Val{dim}, ::Val{N}, d_recvQ::CuArray, recvQ,
                                 d_QL, nrealelem) where {dim, N}
  nrecvelem = size(recvQ)[end]
  if nrecvelem > 0
    d_recvQ .= recvQ
    @cuda(threads=ntuple(j->N+1, dim), blocks=nrecvelem,
          knl_transferrecvQ!(Val(dim), Val(N), d_QL, d_recvQ, nrecvelem,
                             nrealelem))
  end
end

function transferrecvQ!(::Val{dim}, ::Val{N}, d_recvQ::Array, recvQ, Q,
                        nrealelem) where {dim, N}
  Q[:, :, nrealelem+1:end] .= recvQ[:, :, :]
end
# }}}

# {{{ GPU kernel wrappers
@hascuda function volumerhs!(::Val{dim}, ::Val{N}, d_rhsC::CuArray, d_QC,
                             d_vgeoC, d_D, elems) where {dim, N}
  nelem = length(elems)
  @cuda(threads=ntuple(j->N+1, dim), blocks=nelem,
        knl_volumerhs!(Val(dim), Val(N), d_rhsC, d_QC, d_vgeoC, d_D, nelem))
end

@hascuda function facerhs!(::Val{dim}, ::Val{N}, d_rhsL::CuArray, d_QL, d_sgeo,
                           elems, d_vmapM, d_vmapP, d_elemtobndy) where {dim, N}
  nelem = length(elems)
  @cuda(threads=(ntuple(j->N+1, dim-1)..., 1), blocks=nelem,
        knl_facerhs!(Val(dim), Val(N), d_rhsL, d_QL, d_sgeo, nelem, d_vmapM,
                     d_vmapP, d_elemtobndy))
end

@hascuda function updatesolution!(::Val{dim}, ::Val{N}, d_rhsL::CuArray, d_QL,
                                  d_vgeoL, elems, rka, rkb, dt) where {dim, N}
  nelem = length(elems)
  @cuda(threads=ntuple(j->N+1, dim), blocks=nelem,
        knl_updatesolution!(Val(dim), Val(N), d_rhsL, d_QL, d_vgeoL, nelem, rka,
                            rkb, dt))
end
# }}}

# {{{ L2 Energy (for all dimensions)
function L2energysquared(::Val{dim}, ::Val{N}, Q, vgeo, elems) where {dim, N}
  DFloat = eltype(Q)
  Np = (N+1)^dim
  (~, nstate, nelem) = size(Q)

  energy = zero(DFloat)

  @inbounds for e = elems, q = 1:nstate, i = 1:Np
    energy += vgeo[i, _MJ, e] * Q[i, q, e]^2
  end

  energy
end
# }}}

# {{{ RK loop
function lowstorageRK(::Val{dim}, ::Val{N}, mesh, vgeo, sgeo, Q, rhs, D,
                      dt, nsteps, tout, vmapM, vmapP, mpicomm;
                      ArrType=ArrType, plotstep=0) where {dim, N}
  DFloat = eltype(Q)
  mpirank = MPI.Comm_rank(mpicomm)

  # Fourth-order, low-storage, Runge–Kutta scheme of Carpenter and Kennedy
  # (1994) ((5,4) 2N-Storage RK scheme.
  #
  # Ref:
  # @TECHREPORT{CarpenterKennedy1994,
  #   author = {M.~H. Carpenter and C.~A. Kennedy},
  #   title = {Fourth-order {2N-storage} {Runge-Kutta} schemes},
  #   institution = {National Aeronautics and Space Administration},
  #   year = {1994},
  #   number = {NASA TM-109112},
  #   address = {Langley Research Center, Hampton, VA},
  # }
  RKA = (DFloat(0),
         DFloat(-567301805773)  / DFloat(1357537059087),
         DFloat(-2404267990393) / DFloat(2016746695238),
         DFloat(-3550918686646) / DFloat(2091501179385),
         DFloat(-1275806237668) / DFloat(842570457699 ))

  RKB = (DFloat(1432997174477) / DFloat(9575080441755 ),
         DFloat(5161836677717) / DFloat(13612068292357),
         DFloat(1720146321549) / DFloat(2090206949498 ),
         DFloat(3134564353537) / DFloat(4481467310338 ),
         DFloat(2277821191437) / DFloat(14882151754819))

  RKC = (DFloat(0),
         DFloat(1432997174477) / DFloat(9575080441755),
         DFloat(2526269341429) / DFloat(6820363962896),
         DFloat(2006345519317) / DFloat(3224310063776),
         DFloat(2802321613138) / DFloat(2924317926251))

  # Create send and recv request array
  nnabr = length(mesh.nabrtorank)
  sendreq = fill(MPI.REQUEST_NULL, nnabr)
  recvreq = fill(MPI.REQUEST_NULL, nnabr)

  # Create send and recv buffer
  sendQ = zeros(DFloat, (N+1)^dim, size(Q,2), length(mesh.sendelems))
  recvQ = zeros(DFloat, (N+1)^dim, size(Q,2), length(mesh.ghostelems))

  nrealelem = length(mesh.realelems)
  nsendelem = length(mesh.sendelems)
  nrecvelem = length(mesh.ghostelems)
  nelem = length(mesh.elems)

  d_QL, d_rhsL = ArrType(Q), ArrType(rhs)
  d_vgeoL, d_sgeo = ArrType(vgeo), ArrType(sgeo)
  d_vmapM, d_vmapP = ArrType(vmapM), ArrType(vmapP)
  d_sendelems, d_elemtobndy = ArrType(mesh.sendelems), ArrType(mesh.elemtobndy)
  d_sendQ, d_recvQ = ArrType(sendQ), ArrType(recvQ)
  d_D = ArrType(D)

  Qshape    = (fill(N+1, dim)..., size(Q, 2), size(Q, 3))
  vgeoshape = (fill(N+1, dim)..., _nvgeo, size(Q, 3))

  d_QC = reshape(d_QL, Qshape)
  d_rhsC = reshape(d_rhsL, Qshape...)
  d_vgeoC = reshape(d_vgeoL, vgeoshape)

  start_time = t1 = time_ns()
  for step = 1:nsteps
    for s = 1:length(RKA)
      # post MPI receives
      for n = 1:nnabr
        recvreq[n] = MPI.Irecv!((@view recvQ[:, :, mesh.nabrtorecv[n]]),
                                mesh.nabrtorank[n], 777, mpicomm)
      end

      # wait on (prior) MPI sends
      MPI.Waitall!(sendreq)

      # pack data in send buffer
      fillsendQ!(Val(dim), Val(N), sendQ, d_sendQ, d_QL, d_sendelems)

      # post MPI sends
      for n = 1:nnabr
        sendreq[n] = MPI.Isend((@view sendQ[:, :, mesh.nabrtosend[n]]),
                               mesh.nabrtorank[n], 777, mpicomm)
      end

      # volume RHS computation
      volumerhs!(Val(dim), Val(N), d_rhsC, d_QC, d_vgeoC, d_D, mesh.realelems)

      # wait on MPI receives
      MPI.Waitall!(recvreq)

      # copy data to state vectors
      transferrecvQ!(Val(dim), Val(N), d_recvQ, recvQ, d_QL, nrealelem)

      # face RHS computation
      facerhs!(Val(dim), Val(N), d_rhsL, d_QL, d_sgeo, mesh.realelems, d_vmapM,
               d_vmapP, d_elemtobndy)

      # update solution and scale RHS
      updatesolution!(Val(dim), Val(N), d_rhsL, d_QL, d_vgeoL, mesh.realelems,
                      RKA[s%length(RKA)+1], RKB[s], dt)
    end
    if step == 1
      @hascuda synchronize()
      start_time = time_ns()
    end
    if mpirank == 0 && (time_ns() - t1)*1e-9 > tout
      @hascuda synchronize()
      t1 = time_ns()
      avg_stage_time = (time_ns() - start_time) * 1e-9 / ((step-1) * length(RKA))
      @show (step, nsteps, avg_stage_time)
    end
    # TODO: Fix VTK for 1-D
    if dim > 1 && plotstep > 0 && step % plotstep == 0
      Q .= d_QL
      X = ntuple(j->reshape((@view vgeo[:, _x+j-1, :]), ntuple(j->N+1,dim)...,
                            nelem), dim)
      ρ = reshape((@view Q[:, _ρ, :]), ntuple(j->(N+1),dim)..., nelem)
      U = reshape((@view Q[:, _U, :]), ntuple(j->(N+1),dim)..., nelem)
      V = reshape((@view Q[:, _V, :]), ntuple(j->(N+1),dim)..., nelem)
      W = reshape((@view Q[:, _W, :]), ntuple(j->(N+1),dim)..., nelem)
      E = reshape((@view Q[:, _E, :]), ntuple(j->(N+1),dim)..., nelem)
      writemesh(@sprintf("viz/euler%dD_%s_rank_%04d_step_%05d",
                         dim, ArrType, mpirank, step), X...;
                fields=(("ρ", ρ), ("U", U), ("V", V), ("W", W), ("E", E)),
                realelems=mesh.realelems)
    end
  end
  if mpirank == 0
    avg_stage_time = (time_ns() - start_time) * 1e-9 / ((nsteps-1) * length(RKA))
    @show (nsteps, avg_stage_time)
  end
  Q .= d_QL
  rhs .= d_rhsL
end
# }}}

# {{{ atmo driver
function atmo(::Val{dim}, ::Val{N}, mpicomm, ic, mesh, tend;
               meshwarp=(x...)->identity(x),
               tout = 1, ArrType=Array, plotstep=0) where {dim, N}
  DFloat = typeof(tend)

  mpirank = MPI.Comm_rank(mpicomm)
  mpisize = MPI.Comm_size(mpicomm)

  # Partion the mesh using a Hilbert curve based partitioning
  mpirank == 0 && println("[CPU] partiting mesh...")
  mesh = partition(mpicomm, mesh...)

  # Connect the mesh in parallel
  mpirank == 0 && println("[CPU] connecting mesh...")
  mesh = connectmesh(mpicomm, mesh...)

  # Get the vmaps
  mpirank == 0 && println("[CPU] computing mappings...")
  (vmapM, vmapP) = mappings(N, mesh.elemtoelem, mesh.elemtoface,
                            mesh.elemtoordr)

  # Create 1-D operators
  (ξ, ω) = lglpoints(DFloat, N)
  D = spectralderivative(ξ)

  # Compute the geometry
  mpirank == 0 && println("[CPU] computing metrics...")
  (vgeo, sgeo) = computegeometry(Val(dim), mesh, D, ξ, ω, meshwarp, vmapM)
  (nface, nelem) = size(mesh.elemtoelem)

  # Storage for the solution, rhs, and error
  mpirank == 0 && println("[CPU] creating fields (CPU)...")
  Q = zeros(DFloat, (N+1)^dim, _nstate, nelem)
  rhs = zeros(DFloat, (N+1)^dim, _nstate, nelem)

  # setup the initial condition
  mpirank == 0 && println("[CPU] computing initial conditions (CPU)...")
  @inbounds for e = 1:nelem, i = 1:(N+1)^dim
    x, y, z = vgeo[i, _x, e], vgeo[i, _y, e], vgeo[i, _z, e]
    ρ, U, V, W, E = ic(x, y, z)
    Q[i, _ρ, e] = ρ
    Q[i, _U, e] = U
    Q[i, _V, e] = V
    Q[i, _W, e] = W
    Q[i, _E, e] = E
  end

  # Compute time step
  mpirank == 0 && println("[CPU] computing dt (CPU)...")
  base_dt = cfl(Val(dim), Val(N), vgeo, Q, mpicomm) / N^√2
  mpirank == 0 && @show base_dt

  nsteps = ceil(Int64, tend / base_dt)
  dt = tend / nsteps
  mpirank == 0 && @show (dt, nsteps, dt * nsteps, tend)

  # Do time stepping
  stats = zeros(DFloat, 2)
  mpirank == 0 && println("[CPU] computing initial energy...")
  stats[1] = L2energysquared(Val(dim), Val(N), Q, vgeo, mesh.realelems)

  # plot the initial condition
  mkpath("viz")
  # TODO: Fix VTK for 1-D
  if dim > 1
    X = ntuple(j->reshape((@view vgeo[:, _x+j-1, :]), ntuple(j->N+1,dim)...,
                          nelem), dim)
    ρ = reshape((@view Q[:, _ρ, :]), ntuple(j->(N+1),dim)..., nelem)
    U = reshape((@view Q[:, _U, :]), ntuple(j->(N+1),dim)..., nelem)
    V = reshape((@view Q[:, _V, :]), ntuple(j->(N+1),dim)..., nelem)
    W = reshape((@view Q[:, _W, :]), ntuple(j->(N+1),dim)..., nelem)
    E = reshape((@view Q[:, _E, :]), ntuple(j->(N+1),dim)..., nelem)
    writemesh(@sprintf("viz/euler%dD_%s_rank_%04d_step_%05d",
                       dim, ArrType, mpirank, 0), X...;
              fields=(("ρ", ρ), ("U", U), ("V", V), ("W", W), ("E", E)),
              realelems=mesh.realelems)
  end

  mpirank == 0 && println("[DEV] starting time stepper...")
  lowstorageRK(Val(dim), Val(N), mesh, vgeo, sgeo, Q, rhs, D, dt, nsteps, tout,
               vmapM, vmapP, mpicomm; ArrType=ArrType, plotstep=plotstep)

  # TODO: Fix VTK for 1-D
  if dim > 1
    X = ntuple(j->reshape((@view vgeo[:, _x+j-1, :]), ntuple(j->N+1,dim)...,
                          nelem), dim)
    ρ = reshape((@view Q[:, _ρ, :]), ntuple(j->(N+1),dim)..., nelem)
    U = reshape((@view Q[:, _U, :]), ntuple(j->(N+1),dim)..., nelem)
    V = reshape((@view Q[:, _V, :]), ntuple(j->(N+1),dim)..., nelem)
    W = reshape((@view Q[:, _W, :]), ntuple(j->(N+1),dim)..., nelem)
    E = reshape((@view Q[:, _E, :]), ntuple(j->(N+1),dim)..., nelem)
    writemesh(@sprintf("viz/euler%dD_%s_rank_%04d_step_%05d",
                       dim, ArrType, mpirank, nsteps), X...;
              fields=(("ρ", ρ), ("U", U), ("V", V), ("W", W), ("E", E)),
              realelems=mesh.realelems)
  end

  mpirank == 0 && println("[CPU] computing final energy...")
  stats[2] = L2energysquared(Val(dim), Val(N), Q, vgeo, mesh.realelems)

  stats = sqrt.(MPI.allreduce(stats, MPI.SUM, mpicomm))

  if  mpirank == 0
    @show eng0 = stats[1]
    @show engf = stats[2]
    @show Δeng = engf - eng0
  end
end
# }}}

end # module


text:
module TestDefaultlogger

using Historic.Internal: DefaultLogger
using Logging
using Test

function with_defaultlogger(body)
    buffer = (main = IOBuffer(), fallback = IOBuffer())
    context = (:displaysize => (5, 20),)
    io = (
        main = IOContext(buffer.main, context...),
        fallback = IOContext(buffer.fallback, context...),
    )
    logger = DefaultLogger(io.main, ConsoleLogger(io.fallback))
    ans = with_logger(body, logger)
    output = (main = String(take!(buffer.main)), fallback = String(take!(buffer.fallback)))
    return (; output, logger, io, buffer, ans)
end

function test()
    (; output) = with_defaultlogger() do
        @info(:msg1, a = 1, b = 2, c = 3)
        @info(:msg2, a = 1, very_very_very_very_very_very_long_key = 2)
        try
            error("an error to be caught")
        catch err
            @error(:msg3, a = 1, exception = (err, catch_backtrace()))
        end
    end
    @test output.main == """
        msg1 a=1 b=2 c=3
        msg2 a=1...
        Error: msg3 a=1
        """
    @test occursin("an error to be caught", output.fallback)
end

end  # module


text:
include("acs.jl")
include("web.jl")


text:
signal_spec =
[
    (:name, AbstractString, nothing),
    (:init, Any, nothing),
    (:verbose, Bool, false),
    (:expr, AbstractString, nothing),
    (:scale, Any, nothing),
    (:streams, Any, nothing),

]

primitivefactory(:VegaSignal, signal_spec)

max_stars_repo_name: jpthor/vega


text:
######
###### `inv`
######

function frule((_, Δx), ::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    return Ω, -Ω * Δx * Ω
end

function rrule(::typeof(inv), x::AbstractArray)
    Ω = inv(x)
    function inv_pullback(ΔΩ)
        return NoTangent(), -Ω' * ΔΩ * Ω'
    end
    return Ω, inv_pullback
end

#####
##### `*`
#####

frule((_, ΔA, ΔB), ::typeof(*), A, B) = A * B, muladd(ΔA, B, A * ΔB)

frule((_, ΔA, ΔB, ΔC), ::typeof(*), A, B, C) = A*B*C, ΔA*B*C + A*ΔB*C + A*B*ΔC

function rrule(
    ::typeof(*),
    A::AbstractVecOrMat{<:CommutativeMulNumber},
    B::AbstractVecOrMat{<:CommutativeMulNumber},
)
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        dA = @thunk(project_A(Ȳ * B'))
        dB = @thunk(project_B(A' * Ȳ))
        return NoTangent(), dA, dB
    end
    return A * B, times_pullback
end

# Optimized case for StridedMatrixes
# no need to project as already dense, and we are allowed to use InplaceableThunk because
# we know the destination will also be dense. TODO workout how to apply this generally:
# https://github.com/JuliaDiff/ChainRulesCore.jl/issues/411
function rrule(
    ::typeof(*),
    A::StridedMatrix{<:CommutativeMulNumber},
    B::StridedVecOrMat{<:CommutativeMulNumber},
)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        dA = InplaceableThunk(
            X̄ -> mul!(X̄, Ȳ, B', true, true),
            @thunk(Ȳ * B'),
        )
        dB = InplaceableThunk(
            X̄ -> mul!(X̄, A', Ȳ, true, true),
            @thunk(A' * Ȳ),
        )
        return NoTangent(), dA, dB
    end
    return A * B, times_pullback
end

function rrule(
    ::typeof(*),
    A::AbstractVector{<:CommutativeMulNumber},
    B::AbstractMatrix{<:CommutativeMulNumber},
)
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        @assert size(B, 1) === 1   # otherwise primal would have failed.
        return (
            NoTangent(),
            InplaceableThunk(
                X̄ -> mul!(X̄, Ȳ, vec(B'), true, true),
                @thunk(project_A(Ȳ * vec(B'))),
            ),
            InplaceableThunk(
                X̄ -> mul!(X̄, A', Ȳ, true, true),
                @thunk(project_B(A' * Ȳ)),
            )
        )
    end
    return A * B, times_pullback
end

#####
##### `*` matrix-scalar_rule
#####

function rrule(
   ::typeof(*), A::CommutativeMulNumber, B::AbstractArray{<:CommutativeMulNumber}
)
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        return (
            NoTangent(),
            @thunk(project_A(dot(Ȳ, B)')),
            InplaceableThunk(
                X̄ -> mul!(X̄, conj(A), Ȳ, true, true),
                @thunk(project_B(A' * Ȳ)),
            )
        )
    end
    return A * B, times_pullback
end

function rrule(
    ::typeof(*), B::AbstractArray{<:CommutativeMulNumber}, A::CommutativeMulNumber
)
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)
    function times_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        return (
            NoTangent(),
            InplaceableThunk(
                X̄ -> mul!(X̄, conj(A), Ȳ, true, true),
                @thunk(project_B(A' * Ȳ)),
            ),
            @thunk(project_A(dot(Ȳ, B)')),
        )
    end
    return A * B, times_pullback
end

#####
##### fused 3-argument *
#####

if VERSION > v"1.7.0-"

    @eval using LinearAlgebra: mat_mat_scalar, mat_vec_scalar, StridedMaybeAdjOrTransMat

    function rrule(
        ::typeof(mat_mat_scalar),
        A::StridedMaybeAdjOrTransMat{<:CommutativeMulNumber},
        B::StridedMaybeAdjOrTransMat{<:CommutativeMulNumber},
        γ::CommutativeMulNumber
    )
        project_A = ProjectTo(A)
        project_B = ProjectTo(B)
        project_γ = ProjectTo(γ)
        C = mat_mat_scalar(A, B, γ)
        function mat_mat_scalar_back(Ȳraw)
            Ȳ = unthunk(Ȳraw)
            Athunk = InplaceableThunk(
                dA -> mul!(dA, Ȳ, B', conj(γ), true),
                @thunk(project_A(mat_mat_scalar(Ȳ, B', conj(γ)))),
            )
            Bthunk = InplaceableThunk(
                dB -> mul!(dB, A', Ȳ, conj(γ), true),
                @thunk(project_B(mat_mat_scalar(A', Ȳ, conj(γ)))),
            )
            γthunk = @thunk if iszero(γ)
                # Could save A*B on the forward pass, but it's messy.
                # This ought to be rare, should guarantee the same type:
                project_γ(dot(mat_mat_scalar(A, B, oneunit(γ)), Ȳ) / one(γ))
            else
                project_γ(dot(C, Ȳ) / conj(γ))
            end
            return (NoTangent(), Athunk, Bthunk, γthunk)
        end
        return C, mat_mat_scalar_back
    end

    function rrule(
        ::typeof(mat_vec_scalar),
        A::StridedMaybeAdjOrTransMat{<:CommutativeMulNumber},
        b::StridedVector{<:CommutativeMulNumber},
        γ::CommutativeMulNumber
    )
        project_A = ProjectTo(A)
        project_b = ProjectTo(b)
        project_γ = ProjectTo(γ)
        y = mat_vec_scalar(A, b, γ)
        function mat_vec_scalar_back(dy_raw)
            dy = unthunk(dy_raw)
            Athunk = InplaceableThunk(
                dA -> mul!(dA, dy, b', conj(γ), true),
                @thunk(project_A(*(dy, b', conj(γ)))),
            )
            Bthunk = InplaceableThunk(
                db -> mul!(db, A', dy, conj(γ), true),
                @thunk(project_b(*(A', dy, conj(γ)))),
            )
            γthunk = @thunk if iszero(γ)
                project_γ(dot(mat_vec_scalar(A, b, oneunit(γ)), dy))
            else
                project_γ(dot(y, dy) / conj(γ))
            end
            return (NoTangent(), Athunk, Bthunk, γthunk)
        end
        return y, mat_vec_scalar_back
    end

end # VERSION

#####
##### `muladd`
#####

function frule((_, ΔA, ΔB, Δz), ::typeof(muladd), A, B, z)
    Ω = muladd(A, B, z)
    return Ω, ΔA * B .+ A * ΔB .+ Δz
end

function rrule(
        ::typeof(muladd),
        A::AbstractMatrix{<:CommutativeMulNumber},
        B::AbstractVecOrMat{<:CommutativeMulNumber},
        z::Union{CommutativeMulNumber, AbstractVecOrMat{<:CommutativeMulNumber}},
    )
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)
    project_z = ProjectTo(z)

    # The useful case, mul! fused with +
    function muladd_pullback_1(ȳ)
        Ȳ = unthunk(ȳ)
        matmul = (
            InplaceableThunk(
                dA -> mul!(dA, Ȳ, B', true, true),
                @thunk(project_A(Ȳ * B')),
            ),
            InplaceableThunk(
                dB -> mul!(dB, A', Ȳ, true, true),
                @thunk(project_B(A' * Ȳ)),
            )
        )
        addon = if z isa Bool
            NoTangent()
        elseif z isa Number
            @thunk(project_z(sum(Ȳ)))
        else
            InplaceableThunk(
                dz -> sum!(dz, Ȳ; init=false),
                @thunk(project_z(sum!(similar(z, eltype(Ȳ)), Ȳ))),
            )
        end
        (NoTangent(), matmul..., addon)
    end
    return muladd(A, B, z), muladd_pullback_1
end

function rrule(
        ::typeof(muladd),
        ut::LinearAlgebra.AdjOrTransAbsVec{<:CommutativeMulNumber},
        v::AbstractVector{<:CommutativeMulNumber},
        z::CommutativeMulNumber,
    )
    project_ut = ProjectTo(ut)
    project_v = ProjectTo(v)
    project_z = ProjectTo(z)

    # This case is dot(u,v)+z, but would also match signature above.
    function muladd_pullback_2(ȳ)
        dy = unthunk(ȳ)
        ut_thunk = InplaceableThunk(
            dut -> dut .+= v' .* dy,
            @thunk(project_ut(v' .* dy)),
        )
        v_thunk = InplaceableThunk(
            dv -> dv .+= ut' .* dy,
            @thunk(project_v(ut' .* dy)),
        )
        (NoTangent(), ut_thunk, v_thunk, z isa Bool ? NoTangent() : project_z(dy))
    end
    return muladd(ut, v, z), muladd_pullback_2
end

function rrule(
        ::typeof(muladd),
        u::AbstractVector{<:CommutativeMulNumber},
        vt::LinearAlgebra.AdjOrTransAbsVec{<:CommutativeMulNumber},
        z::Union{CommutativeMulNumber, AbstractVecOrMat{<:CommutativeMulNumber}},
    )
    project_u = ProjectTo(u)
    project_vt = ProjectTo(vt)
    project_z = ProjectTo(z)

    # Outer product, just broadcasting
    function muladd_pullback_3(ȳ)
        Ȳ = unthunk(ȳ)
        proj = (
            @thunk(project_u(vec(sum(Ȳ .* conj.(vt), dims=2)))),
            @thunk(project_vt(vec(sum(u .* conj.(Ȳ), dims=1))')),
        )
        addon = if z isa Bool
            NoTangent()
        elseif z isa Number
            @thunk(project_z(sum(Ȳ)))
        else
            InplaceableThunk(
                dz -> sum!(dz, Ȳ; init=false),
                @thunk(project_z(sum!(similar(z, eltype(Ȳ)), Ȳ))),
            )
        end
        (NoTangent(), proj..., addon)
    end
    return muladd(u, vt, z), muladd_pullback_3
end

#####
##### `/`
#####

function rrule(::typeof(/), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    Aᵀ, dA_pb = rrule(adjoint, A)
    Bᵀ, dB_pb = rrule(adjoint, B)
    Cᵀ, dS_pb = rrule(\, Bᵀ, Aᵀ)
    C, dC_pb = rrule(adjoint, Cᵀ)
    function slash_pullback(Ȳ)
        # Optimization note: dAᵀ, dBᵀ, dC are calculated no matter which partial you want
        _, dC = dC_pb(Ȳ)
        _, dBᵀ, dAᵀ = dS_pb(unthunk(dC))

        ∂A = last(dA_pb(unthunk(dAᵀ)))
        ∂B = last(dB_pb(unthunk(dBᵀ)))

        (NoTangent(), ∂A, ∂B)
    end
    return C, slash_pullback
end

#####
##### `\`
#####

function rrule(::typeof(\), A::AbstractVecOrMat{<:Real}, B::AbstractVecOrMat{<:Real})
    project_A = ProjectTo(A)
    project_B = ProjectTo(B)

    Y = A \ B
    function backslash_pullback(ȳ)
        Ȳ = unthunk(ȳ)
        ∂A = @thunk begin
            B̄ = A' \ Ȳ
            Ā = -B̄ * Y'
            Ā = add!!(Ā, (B - A * Y) * B̄' / A')
            Ā = add!!(Ā, A' \ Y * (Ȳ' - B̄'A))
            project_A(Ā)
        end
        ∂B = @thunk project_B(A' \ Ȳ)
        return NoTangent(), ∂A, ∂B
    end
    return Y, backslash_pullback

end

#####
##### `\`, `/` matrix-scalar_rule
#####

function frule((_, ΔA, Δb), ::typeof(/), A::AbstractArray{<:CommutativeMulNumber}, b::CommutativeMulNumber)
    return A/b, ΔA/b - A*(Δb/b^2)
end
function frule((_, Δa, ΔB), ::typeof(\), a::CommutativeMulNumber, B::AbstractArray{<:CommutativeMulNumber})
    return B/a, ΔB/a - B*(Δa/a^2)
end

function rrule(::typeof(/), A::AbstractArray{<:CommutativeMulNumber}, b::CommutativeMulNumber)
    Y = A/b
    function slash_pullback_scalar(ȳ)
        Ȳ = unthunk(ȳ)
        Athunk = InplaceableThunk(
            dA -> dA .+= Ȳ ./ conj(b),
            @thunk(Ȳ / conj(b)),
        )
        bthunk = @thunk(-dot(A,Ȳ) / conj(b^2))
        return (NoTangent(), Athunk, bthunk)
    end
    return Y, slash_pullback_scalar
end

function rrule(::typeof(\), b::CommutativeMulNumber, A::AbstractArray{<:CommutativeMulNumber})
    Y, back = rrule(/, A, b)
    function backslash_pullback(dY)  # just reverses the arguments!
        d0, dA, db = back(dY)
        return (d0, db, dA)
    end
    return Y, backslash_pullback
end

#####
##### Negation (Unary -)
#####

frule((_, ΔA), ::typeof(-), A::AbstractArray) = -A, -ΔA

function rrule(::typeof(-), x::AbstractArray)
    function negation_pullback(ȳ)
        return NoTangent(), InplaceableThunk(ā -> ā .-= ȳ, @thunk(-ȳ))
    end
    return -x, negation_pullback
end

#####
##### Addition (Multiarg `+`)
#####

frule((_, ΔAs...), ::typeof(+), As::AbstractArray...) = +(As...), +(ΔAs...)

function rrule(::typeof(+), arrs::AbstractArray...)
    y = +(arrs...)
    arr_axs = map(axes, arrs)
    function add_pullback(dy)
        return (NoTangent(), map(ax -> reshape(dy, ax), arr_axs)...)
    end
    return y, add_pullback
end


text:
"""
    timezone_names() -> Array{AbstractString}

Returns a sorted list of all of the valid names for constructing a `TimeZone`.
"""
function timezone_names()
    # Note: IANA time zone names are typically encoded only in ASCII.
    names = AbstractString[]
    check = Tuple{AbstractString,AbstractString}[(COMPILED_DIR, "")]

    for (dir, partial) in check
        for filename in readdir(dir)
            startswith(filename, ".") && continue

            path = joinpath(dir, filename)
            name = partial == "" ? filename : join([partial, filename], "/")

            if isdir(path)
                push!(check, (path, name))
            else
                push!(names, name)
            end
        end
    end

    return sort!(names)
end

"""
    all_timezones() -> Array{TimeZone}

Returns all pre-computed `TimeZone`s.
"""
function all_timezones()
    results = TimeZone[]
    for name in timezone_names()
        push!(results, TimeZone(name))
    end
    return results
end

"""
    all_timezones(criteria::Function) -> Array{TimeZone}

Returns `TimeZone`s that match the given `criteria` function. The `criteria` function takes
two parameters: UTC transition (`DateTime`) and transition zone (`FixedTimeZone`).

## Examples

Find all time zones which contain an absolute UTC offset greater than 15 hours:

```julia
all_timezones() do dt, zone
    abs(zone.offset.std) > Dates.Second(Dates.Hour(15))
end
```

Determine all time zones which have a non-hourly daylight saving time offset:

```julia
all_timezones() do dt, zone
    zone.offset.dst % Dates.Second(Dates.Hour(1)) != 0
end
```
"""
function all_timezones(criteria::Function)
    results = TimeZone[]
    for tz in all_timezones()
        if isa(tz, FixedTimeZone)
            criteria(typemin(DateTime), tz) && push!(results, tz)
        else
            for t in tz.transitions
                if criteria(t.utc_datetime, t.zone)
                    push!(results, tz)
                    break
                end
            end
        end
    end
    return results
end

"""
    timezones_from_abbr(abbr) -> Array{TimeZone}

Returns all `TimeZone`s that have the specified abbrevation
"""
function timezones_from_abbr end

function timezones_from_abbr(abbr::Symbol)
    results = TimeZone[]
    for tz in all_timezones()
        if isa(tz, FixedTimeZone)
            tz.name == abbr && push!(results, tz)
        else
            for t in tz.transitions
                if t.zone.name == abbr
                    push!(results, tz)
                    break
                end
            end
        end
    end
    return results
end
timezones_from_abbr(abbr::AbstractString) = timezones_from_abbr(Symbol(abbr))

"""
    timezone_abbrs -> Array{AbstractString}

Returns a sorted list of all pre-computed time zone abbrevations.
"""
function timezone_abbrs()
    abbrs = Set{AbstractString}()
    for tz in all_timezones()
        if isa(tz, FixedTimeZone)
            push!(abbrs, string(tz.name))
        else
            for t in tz.transitions
                push!(abbrs, string(t.zone.name))
            end
        end
    end
    return sort!(collect(abbrs))
end


text:
"""Simplicial Complex Type

Create a simplicial complex with simplices with a value type `P`.
"""
mutable struct SimplicialComplex <: AbstractComplex
    cells::Dict{Int,Vector{AbstractSimplex}}       # cells per dimension
    order::Vector{Dict{UInt64,Int}}  # total order of cells per dimension
end
show(io::IO, cplx::SimplicialComplex) = print(io, "SimplicialComplex($(size(cplx)))")

#
# Constructors
#
SimplicialComplex() = SimplicialComplex(Dict{Int,Vector{AbstractSimplex}}(), Dict{UInt64,Int}[])
# (::Type{SimplicialComplex{S}})() where {S<:AbstractSimplex} = SimplicialComplex(S)

function SimplicialComplex(splxs::AbstractSimplex...)
    cplx = SimplicialComplex()
    added = Set{AbstractSimplex}()
    toprocess = Vector{AbstractSimplex}()
    for splx in splxs
        push!(toprocess, splx) # add simplex to process its faces
        while(length(toprocess) > 0)
            tmp = pop!(toprocess) # processing simplex faces
            tmp in added && continue # skip if already processed
            addsimplex!(cplx, tmp) # add the simplex to the complex
            push!(added, tmp) # mark as processed
            for f in faces(tmp) # add simplex faces for processing
                push!(toprocess, f)
            end
        end
    end

    return cplx
end

SimplicialComplex(splxs::AbstractVector{AbstractSimplex}) = SimplicialComplex(splxs...)

# ---------------
# Private Methods
# ---------------

"""Add a simplex to the complex (as is) and return a complex size, a dimension of simplex and its index in it

This function **doesn't** add missing faces of the simplex to the complex. Use `addsimplices!` function for instead.
"""
function addsimplex!(cplx::SimplicialComplex, splx::AbstractSimplex)
    d = dim(splx)
    d < 0 && return nothing
    if !haskey(cplx.cells, d)
        cplx.cells[d] = AbstractSimplex[]
    end
    if length(cplx.order) <= d
        # println("Missing dimensions: $(length(cplx.order)) need to be $(d+1)")
        for i in length(cplx.order):d
            push!(cplx.order, Dict{UInt64,Int}())
        end
    end
    push!(cplx.cells[d], splx)
    cplx.order[d+1][hash(splx)] = length(cplx.cells[d]) # length(cplx.order)+1
    return splx
end

"""Add a simplex to the complex and all of its faces recursivly and return a complex size, a dimension of simplex and its index in it"""
function addsimplices!(cplx::SimplicialComplex, splx::AbstractSimplex)
    added = Set{AbstractSimplex}() # already cached
    toprocess = AbstractSimplex[]
    ret = AbstractSimplex[]
    addedidxs = NTuple{3,Integer}[]

    # add simplex to complex
    s = addsimplex!(cplx, splx)
    push!(addedidxs, (sum(size(cplx)), dim(s), hash(s)))
    push!(ret, s)
    for f in faces(splx) # add simplex faces for processing
        push!(toprocess, f)
    end

    # process simplex faces and add missing simplexes
    while(length(toprocess) > 0)
        tmp = pop!(toprocess) # processing simplex faces
        d = dim(tmp)

        tmp ∈ cplx  && continue # skip if already in complex
        tmp ∈ added && continue # skip if already processed

        s = addsimplex!(cplx, tmp) # add simplex to the complex
        push!(addedidxs, (sum(size(cplx)), dim(s), hash(s)))
        push!(added, tmp) # mark as processed
        push!(ret, s)

        dim(tmp) == 0 && continue # do not create faces for 0-simplexes
        for f in faces(tmp) # add simplex faces for processing
            push!(toprocess, f)
        end
    end

    return ret
end

#
# AbstractComplex Interface
#
eltype(cplx::SimplicialComplex) = Simplex

function cells(cplx::SimplicialComplex)
    length(cplx.cells) == 0 && return Vector{AbstractSimplex}[] # no cell in complex
    dims = maximum(keys(cplx.cells))
    return [haskey(cplx.cells, d) ? cplx.cells[d] : AbstractSimplex[] for d in 0:dims]
end

cells(cplx::SimplicialComplex, d::Int) = get(cplx.cells, d,  AbstractSimplex[])

function faces(cplx::SimplicialComplex, ci::IX) where {IX<:Integer}
    idxs = IX[]
    d = findfirst(d->ci ∈ keys(d), cplx.order)
    d === nothing && return idxs
    splx = cplx[ci, d-1]
    for f in faces(splx)
        push!(idxs, hash(f))
    end
    return idxs
end

function cofaces(cplx::SimplicialComplex, ci::IX) where {IX<:Integer}
    idxs = IX[]
    d = findfirst(d->ci ∈ keys(d), cplx.order)
    d === nothing || d > dim(cplx) && return idxs
    for cid in keys(cplx.order[d+1])
        ci ∈ faces(cplx, cid) && push!(idxs, cid)
    end
    return idxs
end

function boundary(cplx::SimplicialComplex, idx::IX, d::Int, ::Type{R}) where {R, IX<:Integer}
    ch = Chain(d-1, IX, R)
    d == 0 && return ch

    splx = cplx[idx, d]
    splx === nothing && return ch

    pos = one(R)
    neg = -one(R)
    sgn = true
    for face in faces(splx)
        push!(ch, hash(face)=>(sgn ? pos : neg))
        sgn = !sgn
    end

    return ch
end

function coboundary(cplx::SimplicialComplex, idx::IX, d::Int, ::Type{R}) where {R, IX<:Integer}
    ch = Chain(d+1, IX, R)
    d == 0 && return ch

    # (δϕ)(c) = ϕ(∂c)
    for c in cells(cplx, d+1)
        i = hash(c)
        for (elem, coef) in boundary(cplx, i, d+1, R)
            elem == idx && push!(ch, i=>coef)
        end
    end

    return ch
end

push!(cplx::SimplicialComplex, splx::AbstractSimplex; recursive=false) =
    recursive ? addsimplices!(cplx, splx) : [addsimplex!(cplx, splx)]

function position(cplx::SimplicialComplex, idx::Integer, d::Int)
    length(cplx.order) <= d && return 0
    return get(cplx.order[d+1], idx, 0)
end

#
# Miscellaneous
#

# similar(cplx::SimplicialComplex) = SimplicialComplex()

function read(io::IO, s::S, ::Type{T}) where {S<:AbstractComplex, T<:AbstractSimplex}
    while !eof(io)
        l = chomp(readline(io))
        push!(s, parse(T, l), recursive=false)
    end
    return s
end

function write(io::IO, cplx::SimplicialComplex)
    zsplxs = Dict{Int,UInt}()
    for s in cells(cplx, 0)
        idx = position(cplx, s)
        zsplxs[first(indices(s))] = idx
        write(io, "$idx")
        write(io, 0x0A)
    end
    for d in 1:dim(cplx)
        for s in cplx.cells[d]
            for (j,i) in enumerate(indices(s))
                write(io, "$(zsplxs[i])")
                if j>d
                    write(io, 0x0A)
                else
                    write(io, ' ')
                end
            end
        end
    end
    return
end


text:
"""
    MatAcoustFluidModule

Module for acoustic-fluid  material.
"""
module MatAcoustFluidModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.MatModule: AbstractMat

"""
    MatAcoustFluid <: AbstractMat

Type for acoustic fluid material.
"""
struct MatAcoustFluid <: AbstractMat
	bulk_modulus::FFlt;# Bulk modulus
	mass_density::FFlt;# Mass density
end

"""
    bulkmodulus(self::MatAcoustFluid)

Return the bulk modulus.
"""
function bulkmodulus(self::MatAcoustFluid)
	return self.bulk_modulus
end

end


text:
"""
    direction1Dp(sym)

Return true if `sym` is 1D direction or false otherwise.

The only possible 1D direction is `:x`.
"""
function direction1Dp end

"""
    direction2Dp(sym)

Return true if `sym` is 2D direction or false otherwise.

 Known directions are:

* `:x` and `:y`. Correlation functions are computed along unit vectors
   `(1, 0)` and `(0, 1)` respectively.
* `:xy` and `:yx`. Correlation functions are computed in diagonal
   directions `(1, 1)` and `(-1, 1)` respectively.
"""
function direction2Dp end

"""
    direction3Dp(sym)

Return true if `sym` is 3D direction or false otherwise.

 Known directions are:

* `:x`, `:y` and `:z`. Correlation functions are computed along unit
   vectors `(1, 0, 0)`, `(0, 1, 0)` and `(0, 0, 1)` respectively.
* `:xy`, `:xz` and `:yz`. Correlation functions are computed in
   diagonal directions `(1, 1, 0)`, `(1, 0, 1)` and `(0, 1, 1)`
   respectively.
* `:yx`, `:zx` and `:zy`. Correlation functions are computed in
   diagonal directions `(-1, 1, 0)`, `(-1, 0, 1)` and `(0, -1, 1)`
   respectively.
* `xyz`, `yxz`, `xzy` and `zyx`. Corresponding directions are
   `(1, 1, 1)`, `(-1, 1, 1)`, `(1, -1, 1)` and `(1, 1, -1)`.
"""
function direction3Dp end

macro def_direction_predicate(name)
    exp =
        quote
            $name(x :: Symbol) = x |> Val |> $name
            $name(  :: Any) = false
        end
    return esc(exp)
end

macro def_direction(predicate, sym)
    return :($predicate(:: Val{$sym}) = true)
end

@def_direction_predicate direction1Dp
@def_direction_predicate direction2Dp
@def_direction_predicate direction3Dp

@def_direction(direction1Dp, :x)

@def_direction(direction2Dp, :x)
@def_direction(direction2Dp, :y)
@def_direction(direction2Dp, :xy)
@def_direction(direction2Dp, :yx)

@def_direction(direction3Dp, :x)
@def_direction(direction3Dp, :y)
@def_direction(direction3Dp, :z)
@def_direction(direction3Dp, :xy)
@def_direction(direction3Dp, :yz)
@def_direction(direction3Dp, :xz)
@def_direction(direction3Dp, :yx)
@def_direction(direction3Dp, :zy)
@def_direction(direction3Dp, :zx)

@def_direction(direction3Dp, :xyz)
@def_direction(direction3Dp, :yxz)
@def_direction(direction3Dp, :xzy)
@def_direction(direction3Dp, :zyx)

"""
    default_directions(array)

Get default direction in which correlation functions are calculated
for the given array.
"""
function default_directions end

default_directions(x :: AbstractArray) = x |> ndims |> Val |> default_directions
default_directions(:: Val{3}) = [:x, :y, :z]
default_directions(:: Val{2}) = [:x, :y]
default_directions(:: Val{1}) = [:x]

"""
    direction_predicate(ndims)

Get direction predicate for specified number of dimensions `ndims`.
"""
function direction_predicate end

direction_predicate(x :: Integer) = x |> Val |> direction_predicate
direction_predicate(:: Val{3}) = direction3Dp
direction_predicate(:: Val{2}) = direction2Dp
direction_predicate(:: Val{1}) = direction1Dp
direction_predicate(:: Any) = error("Wrong number of dimensions")

function check_directions(directions :: Vector{Symbol},
                          shape      :: Tuple,
                          periodic   :: Bool)
    ndims = length(shape)
    predicate = direction_predicate(ndims)
    directions = unique(directions)

    if !all(predicate, directions)
        error("Unknown directions found.")
    end

    cubic = all(x -> x == shape[1], shape)
    axial = all(x -> x ∈ [:x, :y, :z], directions)
    if periodic && !axial && !cubic
        error("Periodic diagonals for non-cubic arrays are not supported")
    end

    return directions
end

"""
    unit_length(direction)

Return length of one unit (pixel or voxel) in chosen direction.

# Examples
```jldoctest
julia> unit_length(:x)
1.0
```

```jldoctest
julia> unit_length(:xy)
1.4142135623730951
```

```jldoctest
julia> unit_length(:xyz)
1.7320508075688772
```
"""
function unit_length end

macro def_unit_length(sym, expr)
    return :(unit_length(:: Val{$sym}) = $expr)
end

unit_length(direction :: Symbol) = direction |> Val |> unit_length
@def_unit_length :x sqrt(1.0)
@def_unit_length :y sqrt(1.0)
@def_unit_length :z sqrt(1.0)

@def_unit_length :xy sqrt(2.0)
@def_unit_length :yx sqrt(2.0)
@def_unit_length :yz sqrt(2.0)
@def_unit_length :zy sqrt(2.0)
@def_unit_length :xz sqrt(2.0)
@def_unit_length :zx sqrt(2.0)

@def_unit_length :xyz sqrt(3.0)
@def_unit_length :xzy sqrt(3.0)
@def_unit_length :zyx sqrt(3.0)
@def_unit_length :yxz sqrt(3.0)


text:
module StylegroundWindow

using Gtk, Gtk.ShortNames, Gtk.GConstants
using ..Ahorn, Maple

stylegroundWindow = nothing
effectSectionGrid = nothing
effectOptions = nothing

effectFieldOrder = String[
    "name", "only", "exclude", "tag",
    "flag", "notflag"
]

parallaxFieldOrder = String[
    "texture", "only", "exclude", "tag",
    "flag", "notflag", "blendmode", "color",
    "x", "y", "scrollx", "scrolly",
    "speedx", "speedy", "alpha"
]

function drawPreview(canvas::Gtk.GtkCanvas, textureOption::Ahorn.Form.Option, colorOption::Ahorn.Form.Option)
    ctx = Gtk.getgc(canvas)
    texture = Ahorn.Form.getValue(textureOption)
    sprite = Ahorn.getSprite(texture, "Gameplay")
    rawColor = Ahorn.Form.getValue(colorOption)

    width, height = sprite.realWidth, sprite.realWidth

    if width > 0 && height > 0
        try
            @assert length(rawColor) == 6
            color = Ahorn.argb32ToRGBATuple(parse(Int, "0x" * rawColor) + 255 << 24) ./ 255.0

            offsetX, offsetY = sprite.offsetX, sprite.offsetY

            set_gtk_property!(canvas, :width_request, width)
            set_gtk_property!(canvas, :height_request, height)

            Ahorn.clearSurface(ctx)
            Ahorn.drawImage(ctx, sprite, -offsetX, -offsetY, tint=color)

        catch e
            Ahorn.clearSurface(ctx)
            Ahorn.drawCenteredText(ctx, "Color is invalid.\nExpected hex triplet without #.", 0, 0, 320, 180)
        end

    else
        Ahorn.clearSurface(ctx)
        Ahorn.drawCenteredText(ctx, "Unable to preview backdrop.", 0, 0, 320, 180)
    end
end

function spritesToBackgroundTextures(sprites::Dict{String, Ahorn.SpriteHolder})
    Ahorn.loadExternalSprites!()

    res = String[]

    for (path, spriteHolder) in sprites
        if startswith(path, "bgs/")
            if spriteHolder.sprite.width > 0 && spriteHolder.sprite.height > 0
                push!(res, path)
            end
        end
    end

    push!(res, "mist")
    push!(res, "purplesunset")
    push!(res, "northernlights")
    push!(res, "darkswamp")

    sort!(res)

    return res
end

parallaxFields = Dict{String, Any}(
    "x" => 0,
    "y" => 0,

    "scrollx" => 1,
    "scrolly" => 1,
    "speedx" => 0,
    "speedy" => 0,

    "alpha" => 1,
    "color" => "FFFFFF",

    "only" => "*",
    "exclude" => "",

    "texture" => "",

    "flipx" => false,
    "flipy" => false,
    "loopx" => true,
    "loopy" => true,

    "flag" => "",
    "notflag" => "",

    "blendmode" => "additive",
    "instantIn" => false,
    "instantOut" => false,
    "fadeIn" => false,

    "tag" => ""
)

# Additional fields to add when creating an Effect
effectFields = Dict{String, Any}(
    "flag" => "",
    "notflag" => "",
    "tag" => "",

    "fg" => true
)

# Fake fields that are not set in the data, but makes sense to have options for
parallaxFakeFields = Dict{String, Any}(
    "fg" => true
)

effectTemplates = Dict{String, Dict{String, Any}}(

)

function updateEffectTemplates()
    for effect in Ahorn.effectPlacements
        e = effect()
        effectTemplates[e.name] = merge(effectFields, e.data)
    end
end

function getParallaxListRow(backdrop::Maple.Parallax, fg::Bool=true)
    return (
        get(backdrop, "texture", ""),
        fg,
        get(backdrop, "x", 0),
        get(backdrop, "y", 0),
        get(backdrop, "only", "*"),
        get(backdrop, "exclude", ""),
    )
end

function getEffectListRow(backdrop::Maple.Effect, fg::Bool=true)
    return (
        backdrop.name,
        fg,
        get(backdrop, "only", "*"),
        get(backdrop, "exclude", "")
    )
end

function getParallaxListRows(style::Maple.Style)
    Maple.expandStylegroundApplies!(style)
    rows = Tuple{String, Bool, Int, Int, String, String, Int}[]

    for (fg, backdrops) in ((true, style.foregrounds), (false, style.backgrounds))
        for (i, backdrop) in enumerate(backdrops)
            if isa(backdrop, Maple.Parallax)
                push!(rows, (getParallaxListRow(backdrop, fg)..., i))
            end
        end
    end

    return rows
end

function getEffectListRows(style::Maple.Style)
    Maple.expandStylegroundApplies!(style)
    rows = Tuple{String, Bool, String, String, Int}[]

    for (fg, backdrops) in ((true, style.foregrounds), (false, style.backgrounds))
        for (i, backdrop) in enumerate(backdrops)
            if isa(backdrop, Maple.Effect)
                push!(rows, (getEffectListRow(backdrop, fg)..., i))
            end
        end
    end

    return rows
end

function getBackdrop(map::Maple.Map, fg::Bool, index::Int)
    backdrops = fg ? map.style.foregrounds : map.style.backgrounds

    return backdrops[index]
end

function getListIndices(list::Ahorn.ListContainer, fgCol::Int, indexCol::Int)
    fgIndices = Int[]
    bgIndices = Int[]

    for row in list.data
        fg = row[fgCol]
        index = row[indexCol]

        push!(fg ? fgIndices : bgIndices, index)
    end

    return fgIndices, bgIndices
end

function moveBackdrop!(target::Array{Maple.Backdrop, 1}, indices::Array{Int, 1}, index::Int, offset::Int)
    indicesIndex = findfirst(isequal(index), indices)
    newIndex = indices[clamp(indicesIndex + offset, 1, length(indices))]

    if index != newIndex && 1 <= index <= length(target)
        value = target[index]
        deleteat!(target, index)
        insert!(target, newIndex, value)

        return true
    end

    return false
end

function getParallaxOptions(fields::Dict{String, Any}, langdata::Ahorn.LangData)
    options = Ahorn.Form.Option[]

    # Make sure we get a new list of valid textures
    dropdownOptions = Dict{String, Any}(
        "texture" => sort(spritesToBackgroundTextures(Ahorn.getAtlas("Gameplay"))),
        "blendmode" => String[
            "additive", "alphablend"
        ]
    )

    names = get(langdata, :names)
    tooltips = get(langdata, :tooltips)

    for (dataName, value) in fields
        symbolDataName = Symbol(dataName)
        keyOptions = get(dropdownOptions, dataName, nothing)
        displayName = haskey(names, symbolDataName) ? names[symbolDataName] : Ahorn.humanizeVariableName(dataName)
        tooltip = Ahorn.expandTooltipText(get(tooltips, symbolDataName, ""))

        push!(options, Ahorn.Form.suggestOption(displayName, value, tooltip=tooltip, dataName=dataName, choices=keyOptions, editable=true))
    end

    return options
end

function getEffectOptions(effect::Maple.Effect, langdata::Ahorn.LangData, fg::Bool=true)
    updateEffectTemplates()

    options = Ahorn.Form.Option[]

    # Make sure we get a new list of valid names
    nameOptions = Dict{String, Any}(
        "name" => sort(collect(keys(effectTemplates))),
    )

    names = get(langdata, :names)
    tooltips = get(langdata, :tooltips)
    dropdownOptions = merge(
        Ahorn.editingOptions(effect),
        nameOptions
    )

    data = effect.data
    data["name"] = effect.name

    canFg, canBg = Ahorn.canFgBg(effect)

    if canFg && canBg
        data["fg"] = fg
    end

    for (dataName, value) in merge(effectFields, data)
        symbolDataName = Symbol(dataName)
        keyOptions = get(dropdownOptions, dataName, nothing)
        displayName = haskey(names, symbolDataName) ? names[symbolDataName] : Ahorn.humanizeVariableName(dataName)
        tooltip = Ahorn.expandTooltipText(get(tooltips, symbolDataName, ""))

        push!(options, Ahorn.Form.suggestOption(displayName, value, tooltip=tooltip, dataName=dataName, choices=keyOptions, editable=true))
    end

    return options
end

function updateEffect(effect::Maple.Effect, data::Dict{String, Any}, fg::Bool=true)
    updateEffectTemplates()

    newName = data["name"]
    res = Effect(newName, Dict{String, Any}())

    if effect.name != newName
        for (attr, value) in effectTemplates[newName]
            res.data[attr] = get(effect.data, attr, value)
        end

        res.data = merge(effectTemplates[newName], effect.data)

    else
        res.data = data
    end

    delete!(res.data, "fg")

    return res
end

function moveIllegalFgBg(map::Maple.Map)
    moved = false
    style = map.style

    for (fg, backdrops) in ((true, style.foregrounds), (false, style.backgrounds))
        for i in length(backdrops):-1:1
            backdrop = backdrops[i]

            if isa(backdrop, Maple.Effect)
                canFg, canBg = Ahorn.canFgBg(backdrop)

                if fg && !canFg
                    moved = true
                    deleteat!(style.foregrounds, i)
                    push!(style.backgrounds, backdrop)

                elseif !fg && !canBg
                    moved = true
                    deleteat!(style.backgrounds, i)
                    push!(style.foregrounds, backdrop)
                end
            end
        end
    end

    return moved
end

function filterFakeFields!(data::Dict{String, Any}, fake::Dict{String, Any})
    return filter(pair -> !(haskey(fake, pair[1])), data)
end

function updateEffectSectionGrid(grid::Gtk.GtkGrid, backdrop::Maple.Effect, fg::Bool=true)
    if effectSectionGrid !== nothing
        Gtk.destroy(effectSectionGrid)
    end

    langdataEffects = get(Ahorn.langdata, [:styleground_window, :effects])
    langdataEffect = get(Ahorn.langdata, [:placements, :effects, Symbol(backdrop.name)])

    langdataCombined = Ahorn.LangData(Dict{Symbol, Any}(
        :names => merge(get(langdataEffect, :names), get(langdataEffects, :names)),
        :tooltips => merge(get(langdataEffect, :tooltips), get(langdataEffects, :tooltips))
    ))

    global effectOptions = getEffectOptions(backdrop, langdataCombined, fg)
    section = Ahorn.Form.Section("Effect", effectOptions, fieldOrder=effectFieldOrder)
    global effectSectionGrid = Ahorn.Form.generateSectionGrid(section, columns=8)

    set_gtk_property!(effectSectionGrid, :hexpand, true)
    Ahorn.Form.setGtkProperty!.(effectOptions, :hexpand, true)

    grid[1:8, 2] = effectSectionGrid

    nameOption = Ahorn.Form.getOptionByDataName(effectOptions, "name")
    fgOption = Ahorn.Form.getOptionByDataName(effectOptions, "fg")

    signal_connect(nameOption.combobox, "changed") do widget
        newName = Ahorn.Form.getValue(nameOption)
        backdrop = Effect(newName, effectTemplates[newName])

        fg = fgOption !== nothing ? Ahorn.Form.getValue(fgOption) : fg

        updateEffectSectionGrid(grid, backdrop, fg)
    end

    showall(grid)
end

function getParallaxGrid(map::Maple.Map)
    grid = Grid()

    parallaxList = Ahorn.generateTreeView(
        ("Backdrop", "Foreground", "X", "Y", "Rooms", "Exclude", "Index"),
        getParallaxListRows(map.style),
        sortable=false,
        visible=[
            true, true, true, true,
            true, true, false
        ]
    )

    scrollableList = ScrolledWindow(vexpand=true, hexpand=true, hscrollbar_policy=Gtk.GtkPolicyType.NEVER)
    push!(scrollableList, parallaxList.tree)

    fields = merge(parallaxFields, parallaxFakeFields)
    langdataParallax = get(Ahorn.langdata, [:styleground_window, :parallax])
    options = getParallaxOptions(fields, langdataParallax)
    section = Ahorn.Form.Section("Parallax", options, fieldOrder=parallaxFieldOrder)
    sectionGrid = Ahorn.Form.generateSectionGrid(section, columns=8)

    set_gtk_property!(sectionGrid, :hexpand, true)
    Ahorn.Form.setGtkProperty!.(options, :hexpand, true)

    textureOption = Ahorn.Form.getOptionByDataName(options, "texture")
    colorOption = Ahorn.Form.getOptionByDataName(options, "color")
    fgOption = Ahorn.Form.getOptionByDataName(options, "fg")

    preview = Canvas(320, 180)

    addButton = Button("Add")
    removeButton = Button("Remove")
    updateButton = Button("Update")
    upButton = Button("↑")
    downButton = Button("↓")

    function parallaxRowHandler(list, row)
        index = row[7]
        fg = row[2]

        backdrop = getBackdrop(map, fg, index)
        Ahorn.Form.setOptionsData!(options, merge(parallaxFields, backdrop.data))
        Ahorn.Form.setValue!(fgOption, fg)
    end

    @guarded signal_connect(addButton, "clicked") do widget
        fg = Ahorn.Form.getValue(fgOption)

        backdrops = fg ? map.style.foregrounds : map.style.backgrounds
        backdrop = Parallax(filterFakeFields!(Ahorn.Form.getOptionsData(options), parallaxFakeFields))

        if hasselection(parallaxList.selection)
            row = parallaxList.data[Ahorn.getSelected(parallaxList, 1)]
            index = row[7]
            targetFg = row[2]

            if targetFg == fg
                insert!(backdrops, index + 1, backdrop)

            else
                push!(backdrops, backdrop)
            end

        else
            push!(backdrops, backdrop)
        end

        Ahorn.updateTreeView!(parallaxList, getParallaxListRows(map.style), Ahorn.currentRow(parallaxList) + 1, updateByReplacement=true)
    end

    @guarded signal_connect(removeButton, "clicked") do widget
        if hasselection(parallaxList.selection)
            row = parallaxList.data[Ahorn.getSelected(parallaxList)]
            index = row[7]
            fg = row[2]

            backdrops = fg ? map.style.foregrounds : map.style.backgrounds
            deleteat!(backdrops, index)

            Ahorn.updateTreeView!(parallaxList, getParallaxListRows(map.style), max(Ahorn.currentRow(parallaxList) - 1, 1), updateByReplacement=true)
        end
    end

    @guarded signal_connect(updateButton, "clicked") do widget
        if hasselection(parallaxList.selection)
            row = parallaxList.data[Ahorn.getSelected(parallaxList)]
            index = row[7]
            fg = row[2]
            fgNew = Ahorn.Form.getValue(fgOption)

            if fg == fgNew
                backdrops = fg ? map.style.foregrounds : map.style.backgrounds
                backdrop = backdrops[index]
                backdrop.data = Ahorn.Form.getOptionsData(options)

            else
                backdrops = fg ? map.style.foregrounds : map.style.backgrounds
                backdropsNew = fgNew ? map.style.foregrounds : map.style.backgrounds
                backdrop = backdrops[index]
                backdrop.data = Ahorn.Form.getOptionsData(options)

                deleteat!(backdrops, index)
                push!(backdropsNew, backdrop)
            end

            Ahorn.updateTreeView!(parallaxList, getParallaxListRows(map.style), Ahorn.currentRow(parallaxList), updateByReplacement=true)
        end
    end

    @guarded signal_connect(upButton, "clicked") do widget
        if hasselection(parallaxList.selection)
            row = parallaxList.data[Ahorn.getSelected(parallaxList)]
            index = row[7]
            fg = row[2]

            fgIndices, bgIndices = getListIndices(parallaxList, 2, 7)
            indices = fg ? fgIndices : bgIndices
            backdrops = fg ? map.style.foregrounds : map.style.backgrounds

            moved = moveBackdrop!(backdrops, indices, index, -1)
            select = moved ? Ahorn.currentRow(parallaxList) - 1 : nothing

            Ahorn.updateTreeView!(parallaxList, getParallaxListRows(map.style), select, updateByReplacement=true)
        end
    end

    @guarded signal_connect(downButton, "clicked") do widget
        if hasselection(parallaxList.selection)
            row = parallaxList.data[Ahorn.getSelected(parallaxList)]
            index = row[7]
            fg = row[2]

            fgIndices, bgIndices = getListIndices(parallaxList, 2, 7)
            indices = fg ? fgIndices : bgIndices
            backdrops = fg ? map.style.foregrounds : map.style.backgrounds

            moved = moveBackdrop!(backdrops, indices, index, 1)
            select = moved ? Ahorn.currentRow(parallaxList) + 1 : nothing

            Ahorn.updateTreeView!(parallaxList, getParallaxListRows(map.style), select, updateByReplacement=true)
        end
    end

    Ahorn.connectChanged(parallaxRowHandler, parallaxList)

    signal_connect(widget -> draw(preview), textureOption.combobox, "changed")
    signal_connect(widget -> draw(preview), colorOption.entry, "changed")

    @guarded draw(preview) do widget
        drawPreview(preview, textureOption, colorOption)
    end

    grid[1:6, 1] = scrollableList
    grid[7:8, 1] = preview
    grid[1:8, 2] = sectionGrid
    grid[1:2, 3] = addButton
    grid[3:4, 3] = removeButton
    grid[5:6, 3] = updateButton
    grid[7, 3] = upButton
    grid[8, 3] = downButton

    return grid
end

function getEffectGrid(map::Maple.Map)
    grid = Grid()

    effectList = Ahorn.generateTreeView(
        ("Backdrop", "Foreground", "Rooms", "Exclude", "Index"),
        getEffectListRows(map.style),
        sortable=false,
        visible=[
            true, true, true, true, false
        ]
    )

    scrollableList = ScrolledWindow(vexpand=true, hexpand=true, hscrollbar_policy=Gtk.GtkPolicyType.NEVER)
    push!(scrollableList, effectList.tree)

    global effectSectionGrid = nothing
    options = []

    addButton = Button("Add")
    removeButton = Button("Remove")
    updateButton = Button("Update")
    upButton = Button("↑")
    downButton = Button("↓")

    # Placeholder effect
    updateEffectTemplates()
    templateEffectName = Ahorn.effectPlacements[1]().name
    templateEffectData = effectTemplates[templateEffectName]
    templateEffect = Effect(templateEffectName, templateEffectData)
    updateEffectSectionGrid(grid, templateEffect, true)

    function effectRowHandler(list, row)
        index = row[5]
        fg = row[2]

        backdrop = getBackdrop(map, fg, index)
        updateEffectSectionGrid(grid, backdrop, fg)
    end

    @guarded signal_connect(addButton, "clicked") do widget
        if hasselection(effectList.selection)
            fgOption = Ahorn.Form.getOptionByDataName(effectOptions, "fg")
            row = effectList.data[Ahorn.getSelected(effectList, 1)]
            index = row[5]
            fg = row[2]
            fgNew = fgOption !== nothing ? Ahorn.Form.getValue(fgOption) : fg

            backdrops = fgNew ? map.style.foregrounds : map.style.backgrounds
            backdrop = Effect(row[1], Dict{String, Any}())
            backdrop = updateEffect(backdrop, Ahorn.Form.getOptionsData(effectOptions), fgNew)

            updateEffectSectionGrid(grid, backdrop, fgNew)

            if fg == fgNew
                insert!(backdrops, index + 1, backdrop)

            else
                push!(backdrops, backdrop)
            end

        else
            nameOption = Ahorn.Form.getOptionByDataName(effectOptions, "name")
            name = nameOption !== nothing ? Ahorn.Form.getValue(nameOption) : Ahorn.effectPlacements[1]().name
            backdrop = Effect(name, effectTemplates[name])

            push!(map.style.foregrounds, backdrop)
        end

        select = moveIllegalFgBg(map) ? 0 : something(Ahorn.getSelected(effectList), 0) + 1
        Ahorn.updateTreeView!(effectList, getEffectListRows(map.style), select, updateByReplacement=true)
    end

    @guarded signal_connect(removeButton, "clicked") do widget
        if hasselection(effectList.selection)
            row = effectList.data[Ahorn.getSelected(effectList)]
            index = row[5]
            fg = row[2]

            backdrops = fg ? map.style.foregrounds : map.style.backgrounds
            deleteat!(backdrops, index)

            Ahorn.updateTreeView!(effectList, getEffectListRows(map.style), max(Ahorn.currentRow(effectList) - 1, 1), updateByReplacement=true)
        end
    end

    @guarded signal_connect(updateButton, "clicked") do widget
        if hasselection(effectList.selection)
            fgOption = Ahorn.Form.getOptionByDataName(effectOptions, "fg")

            row = effectList.data[Ahorn.getSelected(effectList)]
            index = row[5]
            fg = row[2]
            fgNew = fgOption !== nothing ? Ahorn.Form.getValue(fgOption) : fg

            if fg == fgNew
                backdrops = fg ? map.style.foregrounds : map.style.backgrounds
                backdrop = backdrops[index]
                backdrops[index] = updateEffect(backdrop, Ahorn.Form.getOptionsData(effectOptions), fgNew)

            else
                backdrops = fg ? map.style.foregrounds : map.style.backgrounds
                backdrop = backdrops[index]
                backdrop = updateEffect(backdrop, Ahorn.Form.getOptionsData(effectOptions), fgNew)

                backdropsNew = fgNew ? map.style.foregrounds : map.style.backgrounds

                deleteat!(backdrops, index)
                push!(backdropsNew, backdrop)
            end

            select = moveIllegalFgBg(map) ? 0 : something(Ahorn.getSelected(effectList), 0)
            Ahorn.updateTreeView!(effectList, getEffectListRows(map.style), select, updateByReplacement=true)
        end
    end

    @guarded signal_connect(upButton, "clicked") do widget
        if hasselection(effectList.selection)
            row = effectList.data[Ahorn.getSelected(effectList)]
            index = row[5]
            fg = row[2]

            fgIndices, bgIndices = getListIndices(effectList, 2, 5)
            indices = fg ? fgIndices : bgIndices
            backdrops = fg ? map.style.foregrounds : map.style.backgrounds

            moved = moveBackdrop!(backdrops, indices, index, -1)
            select = moved ? Ahorn.currentRow(effectList) - 1 : nothing

            Ahorn.updateTreeView!(effectList, getEffectListRows(map.style), select, updateByReplacement=true)
        end
    end

    @guarded signal_connect(downButton, "clicked") do widget
        if hasselection(effectList.selection)
            row = effectList.data[Ahorn.getSelected(effectList)]
            index = row[5]
            fg = row[2]

            fgIndices, bgIndices = getListIndices(effectList, 2, 5)
            indices = fg ? fgIndices : bgIndices
            backdrops = fg ? map.style.foregrounds : map.style.backgrounds

            moved = moveBackdrop!(backdrops, indices, index, 1)
            select = moved ? Ahorn.currentRow(effectList) : nothing

            Ahorn.updateTreeView!(effectList, getEffectListRows(map.style), select, updateByReplacement=true)
        end
    end

    Ahorn.connectChanged(effectRowHandler, effectList)

    grid[1:8, 1] = scrollableList
    grid[1:2, 3] = addButton
    grid[3:4, 3] = removeButton
    grid[5:6, 3] = updateButton
    grid[7, 3] = upButton
    grid[8, 3] = downButton

    return grid
end

function createStylegroundWindow(widget::Ahorn.MenuItemsTypes)
    if Ahorn.loadedState.map === nothing
        info_dialog("No map is currently loaded.", Ahorn.window)

    else
        @Ahorn.catchall begin
            map = Ahorn.loadedState.map

            Maple.expandStylegroundApplies!(map.style)

            box = Box(:v)

            notebook = Notebook()

            parallaxGrid = getParallaxGrid(map)
            parallaxBox = Box(:v)

            effectGrid = getEffectGrid(map)
            effectBox = Box(:v)

            push!(parallaxBox, parallaxGrid)
            push!(notebook, parallaxBox, "Parallax")

            push!(effectBox, effectGrid)
            push!(notebook, effectBox, "Effect")

            push!(box, notebook)

            if stylegroundWindow !== nothing
                Gtk.destroy(stylegroundWindow)
            end

            window = Window("$(Ahorn.baseTitle) - Stylegrounds", -1, -1, true, icon=Ahorn.windowIcon)
            push!(window, box)

            showall(window)

            global stylegroundWindow = window
        end
    end
end

end

max_stars_repo_name: zyowo/Ahorn

max_issues_repo_name: zyowo/Ahorn

max_forks_repo_name: zyowo/Ahorn


text:
"""
Split the work into `Task`s that are spawned on the local machine.

!!! warning

    Currently (Julia version 1.1), all tasks are scheduled onto the same thread.
    Check out [the docs](https://github.com/JuliaLang/julia/blame/80516ca20297a67b996caa08c38786332379b6a5/doc/src/manual/parallel-computing.md#L213-L217)
    or the help of `asyncmap` for more information.

    However, I personally hope/expect that in a future version of Julia this
    issue will be resolved. Check out [this PR on GitHub](https://github.com/JuliaLang/julia/pull/22631)
    and the [references therein](https://github.com/JuliaLang/julia/pull/22631#issuecomment-476472243).

"""
struct LocalTasks <: DistributionStrategy end

"""
There are several constructs in Julia to use coroutines. The example example
here uses a data-pipeline-ish style to demonstrate the use of `Channel`s.
For a more implicit/pragmatic implementation, take a look at `step!` using
`DArray`s where each `@async` creates and schedules a new `Task`.
"""
function step!(new, old, ::LocalTasks)
    m, n = size(old)

    Indices = UnitRange{Int}
    Tile = Tuple{Indices,Indices}
    tiles = Channel{Tile}(Sys.CPU_THREADS)

    # Define the work to be done asynchronously:
    @async begin
        for c in 2:TILESIZE:n-1
            for r in 2:TILESIZE:m-1
                rows = r:min(r+TILESIZE-1, m-1)
                cols = c:min(c+TILESIZE-1, n-1)
                I = (rows,cols)
                put!(tiles, I)
            end
        end
        close(tiles)
    end

    prepare!(old)

    # Start working on it right away, but wait until all are finished.
    # This will be ensured by the channel being closed.
    @sync for _ in 1:Sys.CPU_THREADS
        @async for I in tiles
            life_rule!(new, old, I, I)
        end
    end
end


text:
"""
    struct Sema <: Any
Hold a pointer to a `clang::Sema` object.
"""
struct Sema
    ptr::CXSema
end

Base.unsafe_convert(::Type{CXSema}, x::Sema) = x.ptr
Base.cconvert(::Type{CXSema}, x::Sema) = x