-
Notifications
You must be signed in to change notification settings - Fork 0
/
compressible_operators.jl
222 lines (177 loc) · 15.3 KB
/
compressible_operators.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
#####
##### Legend:
##### ρ -> density fields
##### ρũ -> momentum fields
##### ρc̃ -> conservative tracer fields
##### K̃ -> diffusivity fields
##### ũ -> velocity fields
##### c̃ -> tracer fields
#####
#####
##### Convert conservative variables to primitive variables
#####
@inline U_over_ρ(i, j, k, grid, ρ, ρu) = @inbounds ρu[i, j, k] / ℑxᶠᵃᵃ(i, j, k, grid, ρ)
@inline V_over_ρ(i, j, k, grid, ρ, ρv) = @inbounds ρv[i, j, k] / ℑyᵃᶠᵃ(i, j, k, grid, ρ)
@inline W_over_ρ(i, j, k, grid, ρ, ρw) = @inbounds ρw[i, j, k] / ℑzᵃᵃᶠ(i, j, k, grid, ρ)
@inline C_over_ρ(i, j, k, grid, ρ, ρc) = @inbounds ρc[i, j, k] / ρ[i, j, k]
#####
##### Kinetic energy
#####
@inline kinetic_energy(i, j, k, grid::AbstractGrid{FT}, ρ, ρũ) where FT =
@inbounds FT(0.5) * ( (ℑxᶜᵃᵃ(i, j, k, grid, ρũ.ρu) / ρ[i, j, k])^2
+ (ℑyᵃᶜᵃ(i, j, k, grid, ρũ.ρv) / ρ[i, j, k])^2
+ (ℑzᵃᵃᶜ(i, j, k, grid, ρũ.ρw) / ρ[i, j, k])^2)
#####
##### Pressure work
#####
@inline p∂u∂x(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
( ℑxᶠᵃᵃ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
* Axᵃᵃᶠ(i, j, k, grid) * U_over_ρ(i, j, k, grid, ρ, ρũ.ρu))
@inline p∂v∂y(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
( ℑyᵃᶠᵃ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
* Ayᵃᵃᶠ(i, j, k, grid) * V_over_ρ(i, j, k, grid, ρ, ρũ.ρv))
@inline p∂w∂z(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
( ℑzᵃᵃᶠ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
* Azᵃᵃᵃ(i, j, k, grid) * W_over_ρ(i, j, k, grid, ρ, ρũ.ρw))
@inline ∂ⱼpuⱼ(i, j, k, grid, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, p∂u∂x, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
+ δyᵃᶜᵃ(i, j, k, grid, p∂v∂y, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)
+ δzᵃᵃᶜ(i, j, k, grid, p∂w∂z, diagnose_p, tvar, gases, gravity, ρ, ρũ, ρc̃)))
#####
##### Tracer advection
#####
@inline advective_tracer_flux_x(i, j, k, grid, ρ, ρu, ρc) = Ax_ψᵃᵃᶠ(i, j, k, grid, ρu) * ℑxᶠᵃᵃ(i, j, k, grid, ρc) / ℑxᶠᵃᵃ(i, j, k, grid, ρ)
@inline advective_tracer_flux_y(i, j, k, grid, ρ, ρv, ρc) = Ay_ψᵃᵃᶠ(i, j, k, grid, ρv) * ℑyᵃᶠᵃ(i, j, k, grid, ρc) / ℑyᵃᶠᵃ(i, j, k, grid, ρ)
@inline advective_tracer_flux_z(i, j, k, grid, ρ, ρw, ρc) = Az_ψᵃᵃᵃ(i, j, k, grid, ρw) * ℑzᵃᵃᶠ(i, j, k, grid, ρc) / ℑzᵃᵃᶠ(i, j, k, grid, ρ)
@inline div_uc(i, j, k, grid, ρ, ρũ, ρc) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, advective_tracer_flux_x, ρ, ρũ.ρu, ρc)
+ δyᵃᶜᵃ(i, j, k, grid, advective_tracer_flux_y, ρ, ρũ.ρv, ρc)
+ δzᵃᵃᶜ(i, j, k, grid, advective_tracer_flux_z, ρ, ρũ.ρw, ρc)))
#####
##### Diffusion
##### FIXME: need a better way to enforce boundary conditions on diffusive fluxes at rigid boundaries
#####
@inline diffusive_flux_x(i, j, k, grid, κᶠᶜᶜ, ρ, ρc) =
ℑxᶠᵃᵃ(i, j, k, grid, ρ) * κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, C_over_ρ, ρ, ρc)
@inline diffusive_flux_y(i, j, k, grid, κᶜᶠᶜ, ρ, ρc) =
ℑyᵃᶠᵃ(i, j, k, grid, ρ) * κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, C_over_ρ, ρ, ρc)
@inline diffusive_flux_z(i, j, k, grid, κᶜᶜᶠ, ρ, ρc) =
ℑzᵃᵃᶠ(i, j, k, grid, ρ) * κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, C_over_ρ, ρ, ρc)
@inline diffusive_tvar_flux_x(i, j, k, grid, κᶠᶜᶜ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc) =
(ℑxᶠᵃᵃ(i, j, k, grid, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃)
* κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid) * ∂xᶠᵃᵃ(i, j, k, grid, C_over_ρ, ρ, ρc))
@inline diffusive_tvar_flux_y(i, j, k, grid, κᶜᶠᶜ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc) =
(ℑyᵃᶠᵃ(i, j, k, grid, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃)
* κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid) * ∂yᵃᶠᵃ(i, j, k, grid, C_over_ρ, ρ, ρc))
@inline diffusive_tvar_flux_z(i, j, k, grid, κᶜᶜᶠ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc) =
(ℑzᵃᵃᶠ(i, j, k, grid, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃)
* κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * ∂zᵃᵃᶠ(i, j, k, grid, C_over_ρ, ρ, ρc))
@inline diffusive_pressure_flux_x(i, j, k, grid, κᶠᶜᶜ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃) =
(ℑxᶠᵃᵃ(i, j, k, grid, ρ) * κᶠᶜᶜ * Axᵃᵃᶠ(i, j, k, grid)
* ∂xᶠᵃᵃ(i, j, k, grid, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃))
@inline diffusive_pressure_flux_y(i, j, k, grid, κᶜᶠᶜ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃) =
(ℑyᵃᶠᵃ(i, j, k, grid, ρ) * κᶜᶠᶜ * Ayᵃᵃᶠ(i, j, k, grid)
* ∂yᵃᶠᵃ(i, j, k, grid, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃))
@inline function diffusive_pressure_flux_z(i, j, k, grid::AbstractGrid{FT}, κᶜᶜᶠ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃) where FT
(k <= 1 || k > grid.Nz) && return zero(FT)
return (ℑzᵃᵃᶠ(i, j, k, grid, ρ) * κᶜᶜᶠ * Azᵃᵃᵃ(i, j, k, grid)
* ∂zᵃᵃᶠ(i, j, k, grid, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃))
end
@inline function ∂ⱼDᶜⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, tracer_index, ρ, ρc, args...)
@inbounds κ = closure.κ[tracer_index]
return (1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, diffusive_flux_x, κ, ρ, ρc)
+ δyᵃᶜᵃ(i, j, k, grid, diffusive_flux_y, κ, ρ, ρc)
+ δzᵃᵃᶜ(i, j, k, grid, diffusive_flux_z, κ, ρ, ρc)))
end
@inline function ∂ⱼtᶜDᶜⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, tvar_diag,
tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc, args...)
@inbounds κ = closure.κ[tracer_index]
return (1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, diffusive_tvar_flux_x, κ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)
+ δyᵃᶜᵃ(i, j, k, grid, diffusive_tvar_flux_y, κ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)
+ δzᵃᵃᶜ(i, j, k, grid, diffusive_tvar_flux_z, κ, tvar_diag, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)))
end
@inline function ∂ⱼDᵖⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, tracer_index,
p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃, args ...)
@inbounds κ = closure.κ[tracer_index]
return (1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, diffusive_pressure_flux_x, κ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃)
+ δyᵃᶜᵃ(i, j, k, grid, diffusive_pressure_flux_y, κ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃)
+ δzᵃᵃᶜ(i, j, k, grid, diffusive_pressure_flux_z, κ, p_over_ρ_diag, tvar, gases, gravity, ρ, ρũ, ρc̃)))
end
#####
##### Momentum advection
#####
@inline momentum_flux_ρuu(i, j, k, grid, ρ, ρu) = @inbounds ℑxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, ρu) * ℑxᶜᵃᵃ(i, j, k, grid, ρu) / ρ[i, j, k]
@inline momentum_flux_ρvv(i, j, k, grid, ρ, ρv) = @inbounds ℑyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, ρv) * ℑyᵃᶜᵃ(i, j, k, grid, ρv) / ρ[i, j, k]
@inline momentum_flux_ρww(i, j, k, grid, ρ, ρw) = @inbounds ℑzᵃᵃᶜ(i, j, k, grid, Az_ψᵃᵃᵃ, ρw) * ℑzᵃᵃᶜ(i, j, k, grid, ρw) / ρ[i, j, k]
@inline momentum_flux_ρuv(i, j, k, grid, ρ, ρu, ρv) = ℑxᶠᵃᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, ρv) * ℑyᵃᶠᵃ(i, j, k, grid, ρu) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
@inline momentum_flux_ρuw(i, j, k, grid, ρ, ρu, ρw) = ℑxᶠᵃᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, ρw) * ℑzᵃᵃᶠ(i, j, k, grid, ρu) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
@inline momentum_flux_ρvu(i, j, k, grid, ρ, ρu, ρv) = ℑyᵃᶠᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, ρu) * ℑxᶠᵃᵃ(i, j, k, grid, ρv) / ℑxyᶠᶠᵃ(i, j, k, grid, ρ)
@inline momentum_flux_ρvw(i, j, k, grid, ρ, ρv, ρw) = ℑyᵃᶠᵃ(i, j, k, grid, Az_ψᵃᵃᵃ, ρw) * ℑzᵃᵃᶠ(i, j, k, grid, ρv) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)
@inline momentum_flux_ρwu(i, j, k, grid, ρ, ρu, ρw) = ℑzᵃᵃᶠ(i, j, k, grid, Ax_ψᵃᵃᶠ, ρu) * ℑxᶠᵃᵃ(i, j, k, grid, ρw) / ℑxzᶠᵃᶠ(i, j, k, grid, ρ)
@inline momentum_flux_ρwv(i, j, k, grid, ρ, ρv, ρw) = ℑzᵃᵃᶠ(i, j, k, grid, Ay_ψᵃᵃᶠ, ρv) * ℑyᵃᶠᵃ(i, j, k, grid, ρw) / ℑyzᵃᶠᶠ(i, j, k, grid, ρ)
@inline div_ρuũ(i, j, k, grid, ρ, ρũ) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶠᵃᵃ(i, j, k, grid, momentum_flux_ρuu, ρ, ρũ.ρu)
+ δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρuv, ρ, ρũ.ρu, ρũ.ρv)
+ δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρuw, ρ, ρũ.ρu, ρũ.ρw)))
@inline div_ρvũ(i, j, k, grid, ρ, ρũ) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρvu, ρ, ρũ.ρu, ρũ.ρv)
+ δyᵃᶠᵃ(i, j, k, grid, momentum_flux_ρvv, ρ, ρũ.ρv)
+ δzᵃᵃᶜ(i, j, k, grid, momentum_flux_ρvw, ρ, ρũ.ρv, ρũ.ρw)))
@inline div_ρwũ(i, j, k, grid, ρ, ρũ) =
(1/Vᵃᵃᶠ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, momentum_flux_ρwu, ρ, ρũ.ρu, ρũ.ρw)
+ δyᵃᶜᵃ(i, j, k, grid, momentum_flux_ρwv, ρ, ρũ.ρv, ρũ.ρw)
+ δzᵃᵃᶠ(i, j, k, grid, momentum_flux_ρww, ρ, ρũ.ρw)))
#####
##### Viscous dissipation
#####
@inline strain_rate_tensor_ux(i, j, k, grid::AbstractGrid{FT}, ρ, ρu) where FT = FT(1/3) * ∂xᶜᵃᵃ(i, j, k, grid, U_over_ρ, ρ, ρu)
@inline strain_rate_tensor_vy(i, j, k, grid::AbstractGrid{FT}, ρ, ρv) where FT = FT(1/3) * ∂yᵃᶜᵃ(i, j, k, grid, V_over_ρ, ρ, ρv)
@inline strain_rate_tensor_wz(i, j, k, grid::AbstractGrid{FT}, ρ, ρw) where FT = FT(1/3) * ∂zᵃᵃᶜ(i, j, k, grid, W_over_ρ, ρ, ρw)
@inline strain_rate_tensor_uy(i, j, k, grid, ρ, ρu, ρv) = ∂yᵃᶠᵃ(i, j, k, grid, U_over_ρ, ρ, ρu) + ∂xᶠᵃᵃ(i, j, k, grid, V_over_ρ, ρ, ρv)
@inline strain_rate_tensor_uz(i, j, k, grid, ρ, ρu, ρw) = ∂zᵃᵃᶠ(i, j, k, grid, U_over_ρ, ρ, ρu) + ∂xᶠᵃᵃ(i, j, k, grid, W_over_ρ, ρ, ρw)
@inline strain_rate_tensor_vx(i, j, k, grid, ρ, ρv, ρu) = ∂xᶠᵃᵃ(i, j, k, grid, V_over_ρ, ρ, ρv) + ∂yᵃᶠᵃ(i, j, k, grid, U_over_ρ, ρ, ρu)
@inline strain_rate_tensor_vz(i, j, k, grid, ρ, ρv, ρw) = ∂zᵃᵃᶠ(i, j, k, grid, V_over_ρ, ρ, ρv) + ∂yᵃᶠᵃ(i, j, k, grid, W_over_ρ, ρ, ρw)
@inline strain_rate_tensor_wx(i, j, k, grid, ρ, ρw, ρu) = ∂xᶠᵃᵃ(i, j, k, grid, W_over_ρ, ρ, ρw) + ∂zᵃᵃᶠ(i, j, k, grid, U_over_ρ, ρ, ρu)
@inline strain_rate_tensor_wy(i, j, k, grid, ρ, ρw, ρv) = ∂yᵃᶠᵃ(i, j, k, grid, W_over_ρ, ρ, ρw) + ∂zᵃᵃᶠ(i, j, k, grid, V_over_ρ, ρ, ρv)
@inline viscous_flux_ux(i, j, k, grid, νᶜᶜᶜ, ρ, ρu) = @inbounds ρ[i, j, k] * νᶜᶜᶜ * Axᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_ux(i, j, k, grid, ρ, ρu)
@inline viscous_flux_vy(i, j, k, grid, νᶜᶜᶜ, ρ, ρv) = @inbounds ρ[i, j, k] * νᶜᶜᶜ * Ayᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_vy(i, j, k, grid, ρ, ρv)
@inline viscous_flux_wz(i, j, k, grid, νᶜᶜᶜ, ρ, ρw) = @inbounds ρ[i, j, k] * νᶜᶜᶜ * Azᵃᵃᵃ(i, j, k, grid) * strain_rate_tensor_wz(i, j, k, grid, ρ, ρw)
@inline viscous_flux_uy(i, j, k, grid, νᶠᶠᶜ, ρ, ρu, ρv) = ℑxyᶠᶠᵃ(i, j, k, grid, ρ) * νᶠᶠᶜ * Ayᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_uy(i, j, k, grid, ρ, ρu, ρv)
@inline viscous_flux_uz(i, j, k, grid, νᶠᶜᶠ, ρ, ρu, ρw) = ℑxzᶠᵃᶠ(i, j, k, grid, ρ) * νᶠᶜᶠ * Azᵃᵃᵃ(i, j, k, grid) * strain_rate_tensor_uz(i, j, k, grid, ρ, ρu, ρw)
@inline viscous_flux_vx(i, j, k, grid, νᶠᶠᶜ, ρ, ρv, ρu) = ℑxyᶠᶠᵃ(i, j, k, grid, ρ) * νᶠᶠᶜ * Axᵃᵃᶜ(i, j, k, grid) * strain_rate_tensor_vx(i, j, k, grid, ρ, ρv, ρu)
@inline viscous_flux_vz(i, j, k, grid, νᶜᶠᶠ, ρ, ρv, ρw) = ℑyzᵃᶠᶠ(i, j, k, grid, ρ) * νᶜᶠᶠ * Azᵃᵃᵃ(i, j, k, grid) * strain_rate_tensor_vz(i, j, k, grid, ρ, ρv, ρw)
@inline viscous_flux_wx(i, j, k, grid, νᶠᶜᶠ, ρ, ρw, ρu) = ℑxzᶠᵃᶠ(i, j, k, grid, ρ) * νᶠᶜᶠ * Axᵃᵃᶠ(i, j, k, grid) * strain_rate_tensor_wx(i, j, k, grid, ρ, ρw, ρu)
@inline viscous_flux_wy(i, j, k, grid, νᶜᶠᶠ, ρ, ρw, ρv) = ℑyzᵃᶠᶠ(i, j, k, grid, ρ) * νᶜᶠᶠ * Ayᵃᵃᶠ(i, j, k, grid) * strain_rate_tensor_wy(i, j, k, grid, ρ, ρw, ρv)
@inline ∂ⱼτ₁ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶠᵃᵃ(i, j, k, grid, viscous_flux_ux, closure.ν, ρ, ρũ.ρu)
+ δyᵃᶜᵃ(i, j, k, grid, viscous_flux_uy, closure.ν, ρ, ρũ.ρu, ρũ.ρv)
+ δzᵃᵃᶜ(i, j, k, grid, viscous_flux_uz, closure.ν, ρ, ρũ.ρu, ρũ.ρw)))
@inline ∂ⱼτ₂ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
(1/Vᵃᵃᶜ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, viscous_flux_vx, closure.ν, ρ, ρũ.ρv, ρũ.ρu)
+ δyᵃᶠᵃ(i, j, k, grid, viscous_flux_vy, closure.ν, ρ, ρũ.ρv)
+ δzᵃᵃᶜ(i, j, k, grid, viscous_flux_vz, closure.ν, ρ, ρũ.ρv, ρũ.ρw)))
@inline ∂ⱼτ₃ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
(1/Vᵃᵃᶠ(i, j, k, grid)
* ( δxᶜᵃᵃ(i, j, k, grid, viscous_flux_wx, closure.ν, ρ, ρũ.ρw, ρũ.ρu)
+ δyᵃᶜᵃ(i, j, k, grid, viscous_flux_wy, closure.ν, ρ, ρũ.ρw, ρũ.ρv)
+ δzᵃᵃᶠ(i, j, k, grid, viscous_flux_wz, closure.ν, ρ, ρũ.ρw)))
@inline u₁∂ⱼτ₁ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
U_over_ρ(i, j, k, grid, ρ, ρũ.ρu) * ∂ⱼτ₁ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)
@inline u₂∂ⱼτ₂ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
V_over_ρ(i, j, k, grid, ρ, ρũ.ρv) * ∂ⱼτ₂ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)
@inline u₃∂ⱼτ₃ⱼ(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
W_over_ρ(i, j, k, grid, ρ, ρũ.ρw) * ∂ⱼτ₃ⱼ(i, j, k, grid, closure, ρ, ρũ, args...)
@inline Q_dissipation(i, j, k, grid, closure::ConstantIsotropicDiffusivity, ρ, ρũ, args...) =
( ℑxᶜᵃᵃ(i, j, k, grid, u₁∂ⱼτ₁ⱼ, closure, ρ, ρũ, args...)
+ ℑyᵃᶜᵃ(i, j, k, grid, u₂∂ⱼτ₂ⱼ, closure, ρ, ρũ, args...)
+ ℑzᵃᵃᶜ(i, j, k, grid, u₃∂ⱼτ₃ⱼ, closure, ρ, ρũ, args...))