extend syntax, implement ddx, ddy, ddz

code/code_gen/code_gen.py

@@ -29,16 +29,26 @@ calc_grammar = """
     ?atom: NUMBER           -> number
          | "-" atom         -> neg
          | NAME             -> var
+         | "$" NAME         -> env
          | "(" sum ")"
+         | inlinefunc "(" sum ")" -> icall
          | mathfunc "(" sum ")" -> fcall
          | derivative "(" NAME ")" -> dnx
 
+    ?inlinefunc: "sqr" -> sqr
+             | "pow3" -> pow3
+
     ?mathfunc: "log" -> log
              | "exp" -> exp
              | "sqrt" -> sqrt
+             | "rxn_rate" -> rxn_rate
 
     ?derivative: "dx"  -> dx
                | "d2x" -> d2x
+               | "dy" -> dy
+               | "d2y" -> d2y
+               | "dz" -> dz
+               | "d2z" -> d2z
 
     %import common.CNAME -> NAME
     %import common.NUMBER
@@ -85,12 +95,21 @@ class CalculateTree(Transformer):
     def var(self, name):
         return name + "(i,j,k)"
 
+    def env(self, name):
+        return name
+
     def dnx (self, a, b):
-        signature = "{}_{}".format(a, b)
+        signature = "{}_{}".format(a.data, b)
 
-        self.derivatives[signature] = (a, b)
+        self.derivatives[signature] = (a.data, b)
 
-        return "{}_{}".format(a, b) + "(i,j,k)"
+        return "{}_{}".format(a.data, b) + "(i,j,k)"
+
+    def icall (self, a, b):
+        if a.data == "sqr":
+            return "(({0})*({0}))".format(b)
+        elif a.data == "pow3":
+            return "(({0})*({0})*({0}))".format(b)
 
     def fcall (self, a, b):
         return "( {} ( {} ) )".format(a, b)
@@ -112,9 +131,8 @@ class CalculateTree(Transformer):
 
     log = lambda self : "log"
     exp = lambda self : "exp"
-    dx  = lambda self : "dx"
-    d2x = lambda self : "d2x"
     sqrt = lambda self : "sqrt"
+    rxn_rate = lambda self : "rxn_rate"
 
     def array_decl (self):
 

code/code_gen/resources/m_template.f90

@@ -2,6 +2,7 @@ module m_{0[module_name]}
 
 use m_parameters
 use m_arrays
+use m_calculate
 
 implicit none
 

code/code_gen/terms.input

@@ -1,5 +1,5 @@
 [u, v, w, y]
 
 c_auto = 1 - y
-wrate_auto = 21000 * (1 - c_auto) * exp (- 26.7 / (1 + 3 * c_auto))
-dcdx_auto = dx (c_auto)
+wrate_auto = rxn_rate(c_auto)
+fsd_auto = sqrt (sqr(dx(c_auto)) + sqr(dy(c_auto)) + sqr(dz(c_auto)))

code/m_calculate.f90

@@ -7,30 +7,30 @@ module m_calculate
 
     real*8, allocatable :: work(:,:,:,:)
 
-    real*8, allocatable, dimension(:,:) :: xsrc, ysrc, zsrc
-    real*8, allocatable, dimension(:,:) :: xdst, ydst, zdst
+    real*8, allocatable, dimension(:,:) :: xsrc
+    real*8, allocatable, dimension(:,:) :: xdst
+
+    real*8, allocatable, dimension(:,:,:) :: zsrc
+    real*8, allocatable, dimension(:,:,:) :: zdst
 
     integer, parameter :: nb = 16
 
     private :: work, nb
-    private :: xsrc, ysrc, zsrc
-    private :: xdst, ydst, zdst
+    private :: xsrc, xdst, zsrc, zdst
 
 contains
 
      subroutine m_calculate_init
-     
+
      integer :: ierr
 
      allocate(work(nxp, nyp, nzp, 3), stat=ierr)
 
-     allocate(xsrc(nb, nxp), stat=ierr)
-     allocate(ysrc(nb, nyp), stat=ierr)
-     allocate(zsrc(nb, nzp), stat=ierr)
+     allocate(xsrc(nyp, nxp), stat=ierr)
+     allocate(xdst(nyp, nxp), stat=ierr)
 
-     allocate(xdst(nb, nxp), stat=ierr)
-     allocate(ydst(nb, nyp), stat=ierr)
-     allocate(zdst(nb, nzp), stat=ierr)
+     allocate(zsrc(nyp, nzp, nxp), stat=ierr)
+     allocate(zdst(nyp, nzp, nxp), stat=ierr)
 
      end subroutine m_calculate_init
 
@@ -40,40 +40,72 @@ contains
      deallocate(work)
 
      deallocate(xsrc)
-     deallocate(ysrc)
-     deallocate(zsrc)
-
      deallocate(xdst)
-     deallocate(ydst)
+
+     deallocate(zsrc)
      deallocate(zdst)
 
      end subroutine m_calculate_finalize
 
 
-     subroutine ddx(src, dst)
-        
+     subroutine ddx(dst, src)
+
         real*8, dimension(nxp,nyp,nzp) :: src, dst
 
         integer :: i, j ,k
 
 
         do k = 1,nzp
-        do j = 1,nyp,nb
 
-        call tp2(xsrc, src(:,j:j+nb-1,k), nb, nxp)
+        call tp2(xsrc, src(:,:,k), nyp, nxp)
 
-        call dfnonp(nxp, hxp, xsrc, xdst, nb, 1)
+        call dfnonp(nxp, hxp, xsrc, xdst, nyp, 1)
 
-        call tp2(dst(:,j:j+nb-1,k), xdst, nxp, nb)
+        call tp2(dst(:,:,k), xdst, nxp, nyp)
 
         end do
-        end do
 
      end subroutine ddx
 
+
+     subroutine ddy(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp) :: src, dst
+
+        integer :: i, j ,k
+
+
+        do k = 1,nzp
+
+        call dfp(nyp, hyp, src(:,:,k), dst(:,:,k), nxp, 2)
+
+        end do
+
+     end subroutine ddy
+
+
+     subroutine ddz(dst, src)
+
+        real*8, dimension(nxp,nyp,nzp) :: src, dst
+
+        integer :: i, j ,k
+
+
+        call tp2(zsrc, src, nyp, nzp*nxp)
+
+        do k = 1,nzp
+
+        call dfp(nzp, hzp, zsrc, zdst, nyp, 3)
+
+        end do
+
+        call tp2(dst, zdst, nxp, nyp*nzp)
+
+     end subroutine ddz
+
      subroutine tp(a, b, nx)
        ! a(nb,nx) = transpose(b(nx,nb))
-        
+
         integer,intent(in) :: nx
 
         real*8,intent(out) :: a(nb,nx)
@@ -85,26 +117,26 @@ contains
 
      subroutine tp2 (a, b, n1, n2)
        ! a = transpose(b)
-   
+
          implicit none
-   
+
          integer,intent(in) :: n1, n2
          real*8,intent(out) :: a(n1,n2)
          real*8,intent(in) :: b(n2,n1)
          integer :: i,j,ii,jj
-   
+
             DO jj=1,n2,nb
              DO ii=1,n1,nb
-   
+
             DO j=jj,jj+nb-1
              DO i=ii,ii+nb-1
                a(i,j) = b(j,i)
              ENDDO
             ENDDO
-   
+
              ENDDO
             ENDDO
-   
+
      end subroutine tp2
 
      real function rxn_rate (c)

code/m_terms.f90

@@ -2,13 +2,16 @@ module m_terms
 
 use m_parameters
 use m_arrays
+use m_calculate
 
 implicit none
 
 real*8, allocatable, dimension(:,:,:) :: c_auto
-real*8, allocatable, dimension(:,:,:) :: dcdx_auto
+real*8, allocatable, dimension(:,:,:) :: fsd_auto
 real*8, allocatable, dimension(:,:,:) :: wrate_auto
 real*8, allocatable, dimension(:,:,:) :: dx_c_auto
+real*8, allocatable, dimension(:,:,:) :: dy_c_auto
+real*8, allocatable, dimension(:,:,:) :: dz_c_auto
 
 
 contains
@@ -18,9 +21,11 @@ subroutine m_terms_init
 integer :: ierr
 
 allocate(c_auto(nxp,nyp,nzp), stat=ierr) ; c_auto = 0.
-allocate(dcdx_auto(nxp,nyp,nzp), stat=ierr) ; dcdx_auto = 0.
+allocate(fsd_auto(nxp,nyp,nzp), stat=ierr) ; fsd_auto = 0.
 allocate(wrate_auto(nxp,nyp,nzp), stat=ierr) ; wrate_auto = 0.
 allocate(dx_c_auto(nxp,nyp,nzp), stat=ierr) ; dx_c_auto = 0.
+allocate(dy_c_auto(nxp,nyp,nzp), stat=ierr) ; dy_c_auto = 0.
+allocate(dz_c_auto(nxp,nyp,nzp), stat=ierr) ; dz_c_auto = 0.
 
 
 end subroutine m_terms_init
@@ -29,9 +34,11 @@ end subroutine m_terms_init
 subroutine m_terms_finalize
 
 deallocate(c_auto)
-deallocate(dcdx_auto)
+deallocate(fsd_auto)
 deallocate(wrate_auto)
 deallocate(dx_c_auto)
+deallocate(dy_c_auto)
+deallocate(dz_c_auto)
 
 
 end subroutine m_terms_finalize
@@ -54,7 +61,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-dcdx_auto(i,j,k) = dx_c_auto(i,j,k)
+fsd_auto(i,j,k) = ( sqrt ( ( ( ((dx_c_auto(i,j,k))*(dx_c_auto(i,j,k))) + ((dy_c_auto(i,j,k))*(dy_c_auto(i,j,k))) ) + ((dz_c_auto(i,j,k))*(dz_c_auto(i,j,k))) ) ) )
 end do
 end do
 end do
@@ -63,7 +70,7 @@ end do
 do k = 1, nzp
 do j = 1, nyp
 do i = 1, nxp
-wrate_auto(i,j,k) = ( ( 21000.0 * ( 1.0 - c_auto(i,j,k) ) ) * ( exp ( ( ( - 26.7 ) / ( 1.0 + ( 3.0 * c_auto(i,j,k) ) ) ) ) ) )
+wrate_auto(i,j,k) = ( rxn_rate ( c_auto(i,j,k) ) )
 end do
 end do
 end do

code/makefile

@@ -2,7 +2,7 @@
 #Neuman flags
 #flags = -O3 -r8 -i8 -openmp -msse3 -axSSSE3 -m64 -fPIC -i_dynamic
 #Cluster_2(16.161) flags 
-flags = -Wall -O3 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native
+flags = -Wall -O3 -ffree-line-length-0 -fdefault-integer-8 -fdefault-double-8 -fdefault-real-8 -march=native
 
 ifdef CHECK
 flags += -fcheck=all