From 02a54649f705c00666d2c778dc65790cf7f3b0a0 Mon Sep 17 00:00:00 2001 From: Yeongdo Park Date: Thu, 15 Jul 2021 05:50:29 +0900 Subject: [PATCH] all terms for CTM figures (bugfix: solved missing self allocations for primary fields) added saving derivatives of averaged values --- code/code_gen/post.py | 21 +++++-- code/m_calculate.f90 | 26 ++++++++ code/termspecs/termspec.20210708 | 100 ++++++++++++++++++++++++++++--- 3 files changed, 135 insertions(+), 12 deletions(-) diff --git a/code/code_gen/post.py b/code/code_gen/post.py index 8914860..2fe3817 100644 --- a/code/code_gen/post.py +++ b/code/code_gen/post.py @@ -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 diff --git a/code/m_calculate.f90 b/code/m_calculate.f90 index b1ef3f9..8ef09f9 100644 --- a/code/m_calculate.f90 +++ b/code/m_calculate.f90 @@ -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 diff --git a/code/termspecs/termspec.20210708 b/code/termspecs/termspec.20210708 index 2322675..2adf479 100644 --- a/code/termspecs/termspec.20210708 +++ b/code/termspecs/termspec.20210708 @@ -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 }