From 856cbfbe49c403210cac3c38098e64ed63526be5 Mon Sep 17 00:00:00 2001 From: ignis Date: Mon, 22 Apr 2019 15:41:40 +0900 Subject: [PATCH] average and weighted average --- code/code_gen/code_gen.py | 55 +++++++----- code/code_gen/terms.input | 6 +- code/m_terms.f90 | 172 ++++++++++++++++++++++++++++++++++++-- 3 files changed, 202 insertions(+), 31 deletions(-) diff --git a/code/code_gen/code_gen.py b/code/code_gen/code_gen.py index 69072c7..5042eaa 100644 --- a/code/code_gen/code_gen.py +++ b/code/code_gen/code_gen.py @@ -17,9 +17,9 @@ calc_grammar = """ ?start: statement* - ?statement: NAME "=" sum -> assign_var - | NAME "=" avg -> assign_avg_var - | varlist + ?statement: NAME "=" sum -> assign_var + | avg "{" [sum ("," sum)*] "}" -> assign_avg_var + | varlist ?sum: product | sum "+" product -> add @@ -38,9 +38,7 @@ calc_grammar = """ | mathfunc "(" sum ")" -> fcall | derivative "(" NAME ")" -> dnx - ?avg: "<" sum ">" -> ucavg - | "<" sum ">u" -> cavgu - | "<" sum ">b" -> cavgb + avg: "avg" [NAME] ?inlinefunc: "sqr" -> sqr | "pow3" -> pow3 @@ -86,13 +84,13 @@ avg_array_sum = """ do k = 1, nzp do j = 1, nyp do i = 1, nxp -{0[0]}(i) = {0[0]}(i) + {0[1]} +{0}(i) = {0}(i) + {1} {2} end do end do end do """ -avg_array_divide = "{0[0]} = {0[0]} / denum" +avg_array_divide = "{0} = {0} / denum {1}" real_array_diff = "call {0[0]} ( {0[0]}_{0[1]}, {0[1]} )" @@ -116,19 +114,26 @@ class CalculateTree(Transformer): self.dependency[name.value] = dep return "{} = {}".format(name, value) - def assign_avg_var(self, name, (value, dep)): - self.averaged[name.value] = value - self.dependency[name.value] = dep - return "{} = {}".format(name, value) + def assign_avg_var(self, weight, *args): #name, (value, dep)): - def ucavg(self, (exp, dep)): - return (exp, dep) + fmt = "avg_{}" + if weight is not None: + wvalue, wdep = self.var(weight) + self.averaged[fmt.format(weight)] = wvalue, None + self.dependency[fmt.format(weight)] = wdep + fmt = weight + "_" + fmt - def cavgu(self, (exp, dep)): - return (exp, dep) + for i, (value, dep) in enumerate(args): + self.averaged[fmt.format(i)] = value, weight + self.dependency[fmt.format(i)] = dep + (wdep if weight is not None else []) - def cavgb(self, (exp, dep)): - return (exp, dep) + return "" + + def avg(self, *args): + try: + return args[0] + except IndexError: + return None def number(self, numeral): return (numeral, []) @@ -250,8 +255,9 @@ class CalculateTree(Transformer): for tup in self.derivatives.iteritems(): code_dict[tup[0]] = real_array_diff.format(tup[1]) + "\n" - for tup in self.averaged.iteritems(): - code_dict[tup[0]] = avg_array_sum.format(tup) + "\n" + wfmt = "* {}(i,j,k)" + for k, (v,w) in self.averaged.iteritems(): + code_dict[k] = avg_array_sum.format(k, v, wfmt.format(w) if w is not None else "") + "\n" for var in self.sort_vars(): f_code = f_code + code_dict[var] @@ -263,8 +269,13 @@ class CalculateTree(Transformer): f_code = "" - for tup in self.averaged.iteritems(): - f_code = f_code + avg_array_divide.format(tup) + "\n" + meanw = "/ avg_{}" + + for k, (v,w) in self.averaged.iteritems(): + if k.startswith("avg"): + f_code = avg_array_divide.format(k, meanw.format(w) if w is not None else "") + "\n" + f_code + else: + f_code = f_code + avg_array_divide.format(k, meanw.format(w) if w is not None else "") + "\n" return f_code diff --git a/code/code_gen/terms.input b/code/code_gen/terms.input index d3f1f32..5c7062b 100644 --- a/code/code_gen/terms.input +++ b/code/code_gen/terms.input @@ -3,4 +3,8 @@ c_auto = 1.0 - y wrate_auto = rxn_rate(c_auto) fsd_auto = sqrt (sqr(ddx(c_auto)) + sqr(ddy(c_auto)) + sqr(ddz(c_auto))) -cbar = + +avg { u, v, w} +avg c_auto { u, v, w} +avg y { u, v, w} +avg fsd_auto { u } diff --git a/code/m_terms.f90 b/code/m_terms.f90 index 881852c..1237c1e 100644 --- a/code/m_terms.f90 +++ b/code/m_terms.f90 @@ -12,7 +12,19 @@ real*8, allocatable, dimension(:,:,:) :: wrate_auto real*8, allocatable, dimension(:,:,:) :: ddz_c_auto real*8, allocatable, dimension(:,:,:) :: ddy_c_auto real*8, allocatable, dimension(:,:,:) :: ddx_c_auto -real*8, allocatable, dimension(:) :: cbar +real*8, allocatable, dimension(:) :: avg_fsd_auto +real*8, allocatable, dimension(:) :: fsd_auto_avg_0 +real*8, allocatable, dimension(:) :: avg_c_auto +real*8, allocatable, dimension(:) :: c_auto_avg_2 +real*8, allocatable, dimension(:) :: y_avg_2 +real*8, allocatable, dimension(:) :: c_auto_avg_1 +real*8, allocatable, dimension(:) :: y_avg_0 +real*8, allocatable, dimension(:) :: y_avg_1 +real*8, allocatable, dimension(:) :: avg_0 +real*8, allocatable, dimension(:) :: avg_1 +real*8, allocatable, dimension(:) :: avg_2 +real*8, allocatable, dimension(:) :: c_auto_avg_0 +real*8, allocatable, dimension(:) :: avg_y contains @@ -27,7 +39,19 @@ allocate(wrate_auto(nxp,nyp,nzp), stat=ierr) ; wrate_auto = 0. allocate(ddz_c_auto(nxp,nyp,nzp), stat=ierr) ; ddz_c_auto = 0. allocate(ddy_c_auto(nxp,nyp,nzp), stat=ierr) ; ddy_c_auto = 0. allocate(ddx_c_auto(nxp,nyp,nzp), stat=ierr) ; ddx_c_auto = 0. -allocate(cbar(nxp), stat=ierr) ; cbar = 0. +allocate(avg_fsd_auto(nxp), stat=ierr) ; avg_fsd_auto = 0. +allocate(fsd_auto_avg_0(nxp), stat=ierr) ; fsd_auto_avg_0 = 0. +allocate(avg_c_auto(nxp), stat=ierr) ; avg_c_auto = 0. +allocate(c_auto_avg_2(nxp), stat=ierr) ; c_auto_avg_2 = 0. +allocate(y_avg_2(nxp), stat=ierr) ; y_avg_2 = 0. +allocate(c_auto_avg_1(nxp), stat=ierr) ; c_auto_avg_1 = 0. +allocate(y_avg_0(nxp), stat=ierr) ; y_avg_0 = 0. +allocate(y_avg_1(nxp), stat=ierr) ; y_avg_1 = 0. +allocate(avg_0(nxp), stat=ierr) ; avg_0 = 0. +allocate(avg_1(nxp), stat=ierr) ; avg_1 = 0. +allocate(avg_2(nxp), stat=ierr) ; avg_2 = 0. +allocate(c_auto_avg_0(nxp), stat=ierr) ; c_auto_avg_0 = 0. +allocate(avg_y(nxp), stat=ierr) ; avg_y = 0. end subroutine m_terms_init @@ -41,7 +65,19 @@ deallocate(wrate_auto) deallocate(ddz_c_auto) deallocate(ddy_c_auto) deallocate(ddx_c_auto) -deallocate(cbar) +deallocate(avg_fsd_auto) +deallocate(fsd_auto_avg_0) +deallocate(avg_c_auto) +deallocate(c_auto_avg_2) +deallocate(y_avg_2) +deallocate(c_auto_avg_1) +deallocate(y_avg_0) +deallocate(y_avg_1) +deallocate(avg_0) +deallocate(avg_1) +deallocate(avg_2) +deallocate(c_auto_avg_0) +deallocate(avg_y) end subroutine m_terms_finalize @@ -52,6 +88,15 @@ subroutine m_terms_calculate_pass1 integer :: i, j, k +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +avg_1(i) = avg_1(i) + v(i,j,k) +end do +end do +end do + + do k = 1, nzp do j = 1, nyp do i = 1, nxp @@ -60,17 +105,98 @@ end do end do end do + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +avg_c_auto(i) = avg_c_auto(i) + c_auto(i,j,k) +end do +end do +end do + +call ddz ( ddz_c_auto, c_auto ) + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +c_auto_avg_2(i) = c_auto_avg_2(i) + w(i,j,k) * c_auto(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +c_auto_avg_1(i) = c_auto_avg_1(i) + v(i,j,k) * c_auto(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +c_auto_avg_0(i) = c_auto_avg_0(i) + u(i,j,k) * c_auto(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +y_avg_2(i) = y_avg_2(i) + w(i,j,k) * y(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +y_avg_0(i) = y_avg_0(i) + u(i,j,k) * y(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +y_avg_1(i) = y_avg_1(i) + v(i,j,k) * y(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +avg_0(i) = avg_0(i) + u(i,j,k) +end do +end do +end do + call ddy ( ddy_c_auto, c_auto ) do k = 1, nzp do j = 1, nyp do i = 1, nxp -cbar(i) = cbar(i) + c_auto(i,j,k) +avg_2(i) = avg_2(i) + w(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +avg_y(i) = avg_y(i) + y(i,j,k) end do end do end do -call ddz ( ddz_c_auto, c_auto ) call ddx ( ddx_c_auto, c_auto ) do k = 1, nzp @@ -82,6 +208,24 @@ end do end do +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +fsd_auto_avg_0(i) = fsd_auto_avg_0(i) + u(i,j,k) * fsd_auto(i,j,k) +end do +end do +end do + + +do k = 1, nzp +do j = 1, nyp +do i = 1, nxp +avg_fsd_auto(i) = avg_fsd_auto(i) + fsd_auto(i,j,k) +end do +end do +end do + + do k = 1, nzp do j = 1, nyp do i = 1, nxp @@ -102,7 +246,19 @@ real*8 :: denum denum=real(nfiles*nyp*nzp) -cbar = cbar / denum +avg_y = avg_y / denum +avg_2 = avg_2 / denum +avg_1 = avg_1 / denum +avg_0 = avg_0 / denum +avg_c_auto = avg_c_auto / denum +avg_fsd_auto = avg_fsd_auto / denum +fsd_auto_avg_0 = fsd_auto_avg_0 / denum / avg_fsd_auto +c_auto_avg_2 = c_auto_avg_2 / denum / avg_c_auto +y_avg_2 = y_avg_2 / denum / avg_y +c_auto_avg_1 = c_auto_avg_1 / denum / avg_c_auto +y_avg_0 = y_avg_0 / denum / avg_y +y_avg_1 = y_avg_1 / denum / avg_y +c_auto_avg_0 = c_auto_avg_0 / denum / avg_c_auto end subroutine m_terms_average_pass1 @@ -112,5 +268,5 @@ end module m_terms ! {'c_auto': '( 1.0 - y(i,j,k) )', 'fsd_auto': '( sqrt ( ( ( ((ddx_c_auto(i,j,k))*(ddx_c_auto(i,j,k))) + ((ddy_c_auto(i,j,k))*(ddy_c_auto(i,j,k))) ) + ((ddz_c_auto(i,j,k))*(ddz_c_auto(i,j,k))) ) ) )', 'wrate_auto': '( rxn_rate ( c_auto(i,j,k) ) )'} ! {'ddz_c_auto': ('ddz', 'c_auto'), 'ddy_c_auto': ('ddy', 'c_auto'), 'ddx_c_auto': ('ddx', 'c_auto')} -! {'fsd_auto': ['ddx_c_auto', 'ddy_c_auto', 'ddz_c_auto'], 'ddz_c_auto': ['c_auto'], 'c_auto': [], 'ddy_c_auto': ['c_auto'], 'cbar': ['c_auto'], 'ddx_c_auto': ['c_auto'], 'wrate_auto': ['c_auto']} -! ['c_auto', 'ddy_c_auto', 'cbar', 'ddz_c_auto', 'ddx_c_auto', 'fsd_auto', 'wrate_auto'] +! {'fsd_auto': ['ddx_c_auto', 'ddy_c_auto', 'ddz_c_auto'], 'avg_fsd_auto': ['fsd_auto'], u'c_auto_avg_2': ['c_auto'], 'avg_c_auto': ['c_auto'], 'avg_1': [], 'ddz_c_auto': ['c_auto'], 'c_auto': [], u'y_avg_2': [], u'c_auto_avg_1': ['c_auto'], u'y_avg_0': [], u'y_avg_1': [], 'avg_0': [], 'ddy_c_auto': ['c_auto'], 'avg_2': [], u'c_auto_avg_0': ['c_auto'], 'avg_y': [], u'fsd_auto_avg_0': ['fsd_auto'], 'ddx_c_auto': ['c_auto'], 'wrate_auto': ['c_auto']} +! ['avg_1', 'c_auto', 'avg_c_auto', 'ddz_c_auto', u'c_auto_avg_2', u'c_auto_avg_1', u'c_auto_avg_0', u'y_avg_2', u'y_avg_0', u'y_avg_1', 'avg_0', 'ddy_c_auto', 'avg_2', 'avg_y', 'ddx_c_auto', 'fsd_auto', u'fsd_auto_avg_0', 'avg_fsd_auto', 'wrate_auto']