all terms for CTM figures (bugfix: solved missing self allocations for primary fields)

added saving derivatives of averaged values
This commit is contained in:
Yeongdo Park 2021-07-15 05:50:29 +09:00
parent 2e20f720f2
commit 02a54649f7
3 changed files with 135 additions and 12 deletions

View file

@ -264,7 +264,7 @@ class ExpToCode(Transformer):
log = lambda self : "log"
exp = lambda self : "exp"
sqrt = lambda self : "sqrt"
abs = lambda self : "dabs"
abs = lambda self : "abs"
rxn_rate = lambda self : "rxn_rate"
udf = lambda self, a : a.value
@ -877,6 +877,7 @@ class Stage4():
''' analyze liveness and allocate array '''
def __init__ (self, src):
self.src = src
self.primary = src.primary
self.derived = src.derived
self.averaged = src.averaged
self.dependency = src.dependency
@ -922,7 +923,7 @@ class Stage4():
livesets = [set([])] + [set(np.asarray(l)[row]) for row in mask.T]
var2arr = {}
var2arr = { p : p for p in self.primary }
for i, (s0, s1) in enumerate(zip(livesets[:-1], livesets[1:])):
@ -952,23 +953,33 @@ class Stage4():
def write_avg_codes (self, avglist):
avg_array_write = '''
real*8, dimension(nxp) :: xbuffer
integer :: i
open (200, file="qEdge_X.dat")
write (200,*) output_header
do i=1,nxp
write (200,'({0}e20.10)') real(i)*hxp, {1}
end do
close (200)
open (200, file="d1.dat")
{2}
close (200)
open (200, file="d2.dat")
{3}
close (200)
'''
avgarr = "{}(i)"
deriv1_avgarr = """call ddx1d ( xbuffer, {} ) ; write (200,*) xbuffer"""
deriv2_avgarr = """call d2dx1d ( xbuffer, {} ) ; write (200,*) xbuffer"""
write_avg = avg_array_write.format(
len(self.averaged)+1,
", ".join(map(avgarr.format, avglist))
", ".join(map(avgarr.format, avglist)),
"\n".join(map(deriv1_avgarr.format, avglist)),
"\n".join(map(deriv2_avgarr.format, avglist)),
)
return write_avg

View file

@ -44,6 +44,19 @@ contains
end subroutine m_calculate_finalize
subroutine ddx1d(dst, src)
real*8, dimension(1,nxp), intent(in) :: src
real*8, dimension(1,nxp), intent(out) :: dst
integer :: i, j ,k
integer :: ju
call dfnonp(nxp, hxp, src, dst, 1, 1)
end subroutine ddx1d
subroutine ddx(dst, src)
real*8, dimension(nxp,nyp,nzp), intent(in) :: src
@ -111,6 +124,19 @@ contains
end subroutine ddz
subroutine d2dx1d(dst, src)
real*8, dimension(nxp), intent(in) :: src
real*8, dimension(nxp), intent(out) :: dst
integer :: i, j ,k
integer :: ju
call d2fnonp(nxp, hxp, src, dst, 1, 1)
end subroutine d2dx1d
subroutine d2dx(dst, src)
real*8, dimension(nxp,nyp,nzp), intent(in) :: src

View file

@ -35,57 +35,143 @@ wnz_fluc = w' * nz'
vn(latex="v_n") = unx + vny + wnz
divn(latex="(\nabla \cdot \vec{n})") = ddx(nx) + ddy(ny) + ddz(nz)
divn(latex="(\nabla \cdot \mathbf{n})") = ddx(nx) + ddy(ny) + ddz(nz)
divn2 = sqr(divn)
# grad(n) : grad(n)
gradn_gradn(latex="\nabla \vec{n} : \nabla \vec{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(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'
vn_divn (latex="\vec{v}'' \cdot \vec{n}'' (\nabla\cdot \vec{n})''") = (u' * nx' + v' * ny' + w' * nz') * 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
avg { u, v, w, c, f, lapc, ddx_c, ddy_c, ddz_c, d2dx_c, d2dy_c, d2dz_c, uc }
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,
sd,
unxsd,
vn, sd, vnsd,
vn_lvl, sd_lvl, vnsd_lvl,
absk, divn,
absk_lvl,
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,
sd_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
}