# fsd weighted (vn+sd) correlations

[u, v, w, y]

c(latex="c") = 1.0 - y
wrate (latex="\dot{\omega}") = rxn_rate(c)
f (latex="\Sigma_f'")= sqrt (sqr(ddx(c)) + sqr(ddy(c)) + sqr(ddz(c)))

uc = u' * c'

lapc(latex="\Delta c") = d2dx(c) + d2dy(c) + d2dz(c)
sd(latex="S_d") = ($rod * lapc + wrate) / f

nx(latex="n_x") = - ddx(c) / f
ny(latex="n_y") = - ddy(c) / f
nz(latex="n_z") = - ddz(c) / f

nx_lvl = sqr(nx')
ny_lvl = sqr(ny')
nz_lvl = sqr(nz')

nx2 = sqr(nx)
ny2 = sqr(ny)
nz2 = sqr(nz)

nn = nx2 + ny2 + nz2

unx = u * nx
vny = v * ny
wnz = w * nz

unx_fluc = u' * nx'
vny_fluc = v' * ny'
wnz_fluc = w' * nz'

vn(latex="v_n") = unx + vny + wnz

divn(latex="(\nabla \cdot \mathbf{n})") = ddx(nx) + ddy(ny) + ddz(nz)

divn2 = sqr(divn)

# grad(n) : grad(n)
gradn_gradn(latex="\nabla \mathbf{n} : \nabla \mathbf{n}") = sqr(ddx(nx)) + sqr(ddy(ny)) + sqr(ddz(nz)) + 2*(ddx(ny)*ddy(nx)) + 2*(ddx(nz)*ddz(nx)) + 2*(ddy(nz)*ddz(ny))

gradn_gradn_lvl = gradn_gradn' * gradn_gradn'

absk = abs(divn)

divn_lvl = sqr(divn')
absk_lvl = sqr(absk')

unxsd = u + nx * sd

vnsd = vn+sd

vn_lvl = sqr(vn')
sd_lvl = sqr(sd')
vnsd_lvl = sqr(vnsd')

u_lvl = sqr(u')
v_lvl = sqr(v')
w_lvl = sqr(w')

vnsd_nx = vnsd * nx
vnsd_nx_fluc = vnsd' * nx'
vnsd_nx_lvl = sqr(vnsd_nx')
# u * nx_lvl
# nx * unx_fluc
vn_nx = vn' * nx'
sd_nx = sd' * nx'
u_nx_nx = u' * sqr(nx')
v_ny_nx = v' * ny' * nx'
w_nz_nx = w' * nz' * nx'

vnsd_divn = vnsd * divn
vnsd_divn_fluc = vnsd' * divn'
nx_divn = nx' * divn'
u_divn = u' * divn'
v_n_divn (latex="\mathbf{v}'' \cdot \mathbf{n}'' (\nabla\cdot \mathbf{n})''") = (u' * nx' + v' * ny' + w' * nz') * divn'
sd_divn = sd' * divn'

vn_divn_ = vn * divn
vn_divn = vn' * divn'
sd_divn_ = sd * divn

divt_v (latex="\nabla_T \cdot \mathbf{v}") = (1.-nx) * ddx(u) + (1.-ny) * ddy(v) + (1.-nz) * ddz(w)

vnsd_nx_nx_flux = vnsd_nx' * nx'
vnsd_divn_nx_flux = vnsd_divn' * nx'

vnsd_nx_divn_flux = vnsd_nx' * divn'
vnsd_divn_divn_flux = vnsd_divn' * divn'

gradtx_vnsd (latex="\nabla_{Tx} \left(v_n + S_d \right)") = (1.-nx) * ddx(vnsd)
gradty_vnsd (latex="\nabla_{Ty} \left(v_n + S_d \right)") = (1.-ny) * ddy(vnsd)
gradtz_vnsd (latex="\nabla_{Tz} \left(v_n + S_d \right)") = (1.-nz) * ddz(vnsd)

lapt_vnsd (latex="\nabla^2_T \left(v_n + S_d\right)") = (1.-nx) * ddx(gradtx_vnsd) + (1.-ny) * ddy(gradty_vnsd) + (1.-nz) * ddz(gradtz_vnsd)

vnsd_gradn_gradn = vnsd * gradn_gradn

vnsd_nx_absk = vnsd_nx' * absk'
vnsd_divn_absk = vnsd_divn' * absk'

sign_lapt_vnsd = $sign(1.0, divn) * lapt_vnsd
sign_vnsd_gradn_gradn = $sign(1.0, divn) * vnsd_gradn_gradn

dcdn = nx * ddx(c) + ny * ddy(c) + nz * ddz(c)
d2cdn2 = nx * ddx(dcdn) + ny * ddy(dcdn) + nz * ddz(dcdn)





avg {
u, v, w,
u_lvl, v_lvl, w_lvl,
y, c, f,
ddx_c, ddy_c, ddz_c,
lapc, d2dx_c, d2dy_c, d2dz_c,
uc,
dcdn, d2cdn2
}

# <>u
avg y { u, v, w, c, f, lapc, ddx_c, ddy_c, ddz_c, d2dx_c, d2dy_c, d2dz_c }

# <>b
avg c { u, v, w, c, f, lapc, ddx_c, ddy_c, ddz_c, d2dx_c, d2dy_c, d2dz_c }

# <>K
avg d2cdn2 { nx, ny, nz }

# <>f
avg f { 
u, v, w, 
u_lvl, v_lvl, w_lvl,
nx, ny, nz, 
nx2, ny2, nz2, 
nn,
unxsd,
vn, sd, vnsd,
vn_lvl, sd_lvl, vnsd_lvl,
absk, divn,
absk_lvl, divn_lvl,
nx_lvl, ny_lvl, nz_lvl,
unx, vny, wnz,
unx_fluc, vny_fluc, wnz_fluc,
divt_v,
vnsd_nx,
vnsd_nx_fluc,
vnsd_nx_lvl,
vn_nx,
sd_nx,
u_nx_nx ,
v_ny_nx ,
w_nz_nx ,
vnsd_divn,
vnsd_divn_fluc,
nx_divn,
u_divn,
v_n_divn,
sd_divn,
sd_divn_,
vn_divn,
vn_divn_,
vnsd_nx_nx_flux,
vnsd_nx_divn_flux,
vnsd_divn_nx_flux,
vnsd_divn_divn_flux,
gradtx_vnsd,
lapt_vnsd,
vnsd_gradn_gradn,
vnsd_divn_absk,
sign_lapt_vnsd,
sign_vnsd_gradn_gradn
}