Skip to content

Instantly share code, notes, and snippets.

@simonbyrne
Created June 10, 2020 21:48
Show Gist options
  • Save simonbyrne/2c0a121d9b2fd9744b7eae2ff8a9d9f5 to your computer and use it in GitHub Desktop.
Save simonbyrne/2c0a121d9b2fd9744b7eae2ff8a9d9f5 to your computer and use it in GitHub Desktop.
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment