Tendencies.jl

module Tendencies



abstract type Term end
# rows: equations of each variable
abstract type PrognosticQuantity <: Term end
struct Mass <: PrognosticQuantity end
struct Momentum <: PrognosticQuantity end

struct Divergence <: Term
    operand
end
struct Gradient <: Term
    operand
end

# define operators
struct Grad end
const ∇ = Grad()
(::Grad)(operand) = Gradient(operand)
(⋅)(::Grad, operand) = Divergence(operand)

struct TermSum <: Term
    operands
end
Base.(:+)(t::Term...) = TermSum(t)


islinear(::PrognosticQuantity) = true
islinear(d::Divergence) = islinear(d.operand)
islinear(d::Gradient) = islinear(d.operand)
islinear(d::TermSum) = all(islinear, d.operands)

isvertical(::Momentum) = false
isvertical(::VericalProjection) = true

struct Pressure <: DiagnosticQuantity
end

islinear(::Pressure) = false

const ρ = Mass()
const ρu = Momentum()

u = ρu / ρ
p = Pressure()

∂t(ρ) ~ ∇ ⋅ ρu + s(ρ)
S ~ (∇(u) + ∇(u)')/2
τ = -2*ν .* S

ρu_euler = ∇⋅(u ⊗ ρu + p * I)
ρu_diffusive = ∇⋅(ρ * τ)

mooddeell ==  ∇⋅(u ⊗ ρu + p * I)∇⋅(u ⊗ ρu + p * I)


ffoo==∂t(ρu) ~  ∇⋅(u ⊗ ρu + p * I)) + ∇⋅(ρ * τ; numflux = Central())
ffoo['eeuullerr']

# challenges
# - how to "name" subexpressions
#   - numerical fluxes
#   - boundary conditions
#   - time rates


@something begin
    G = g(Qstate, Qaux)
    ∇G = ∇DG(G; numflux=Central())
    H = h(∇G)
end

# fuse these
@something begin
    F1 = f1(Qstate, Qaux)
    filter(F1, ExponentialFilter())
    D1 = ∇DG⋅(F1; numflux=Rusanov())
    F2 = f2(Qstate, H, Qaux)
    D2 = ∇DG⋅(F2; numflux=Central())
    S = s(Qstate, Qaux)
    D1 + D2 + S
end


[[ ∇⋅(u ⊗ ρu + p * I)

]]




# linear
∂t(Mass()) ~ vertical(∇ ⋅ (Momentum()))
∂t(Momentum()) ~ 0



# remainder
∂t(Mass()) ~
∂t(Momentum()) ~ 0



# column: tendencies
abstract type Tendency
end

struct Advection <: Tendency
end
struct PressureGradient <: Tendency
end

# third dimension: splitting
struct IMEXAcoustic{Tendency,Rate} <: Tendency
end

# Advection
## Mass
flux_firstorder(::Mass, ::Advection, model, state) = momentum(model, state)
flux_firstorder(::Mass, ::IMEXAcoustic{Advection,VerticalImplicit}, model, state) = project_vertical(momentum(model, state))
flux_firstorder(::Mass, ::IMEXAcoustic{Advection,Explicit}, model, state) = 0


flux(::Mass, ::FirstOrder, ::Acoustic, ::VerticallyImplicit{rate}, ::CentralNumericalFlux, ::Antialiasing, model, state) = func()

∇⋅(ρu)



## Momentum
flux_firstorder(::Momentum, ::Advection, ...) =
    momentum(state) .* velocity(state)'
flux_firstorder(::Momentum, ::IMEXAcoustic{Advection,Implicit}, ...) =
    0
flux_firstorder(::Momentum, ::IMEXAcoustic{Advection,Explicit}, ...) =
    momentum(state) .* velocity(state)'

## Energy
flux_firstorder(::Energy, ::Advection, ...) =
    momentum(state) * energy_per_mass(state)
flux_firstorder(::Energy, ::IMEXAcoustic{Advection,Implicit}, ...) =
    momentum(state) * ref_energy_per_mass(state)
flux_firstorder(::Energy, ::IMEXAcoustic{Advection,Explicit}, ...) =
    momentum(state) * (energy_per_mass(state) - ref_energy_per_mass(state))

flux_firstorder(::AbstractTracer, ::Advection, ...) =
    momentum(state) * quantity_per_mass(state)'



# PressureGradient
## Momentum
flux_firstorder(::Momentum, ::PressureGradientPlusGravity, ...) =
    (pressure(...) - reference_pressure(...))*I
source(::Momentum, ::PressureGradientPlusGravity, ...) =
    -(state.ρ - reference_mass(...)) * aux.orientation.∇Φ

flux_firstorder(::Momentum, ::IMEXAcoustic{PressureGradient,Implicit},...) =
    linear_pressure(....) * I
flux_firstorder(::Momentum, ::IMEXAcoustic{PressureGradient,Explicit},...) =
    (pressure(...) - linear_pressure(...))*I

flux_firstorder(::Momentum, ::PressureGradient, ::Linear,...) = linear_pressure(....) * I

## Energy
flux_firstorder(::Energy, ::PressureGradient, ...) = u * pressure(...)

flux_firstorder(::Energy, ::Advection, ::Linear, ...) = (ref.ρe + ref.p) / ref.ρ * state.ρu
flux_firstorder(::Energy, ::PressureGradient, ::Linear, ...) = e_pot * state.ρu

abstract type Term
end




end # module