410-bwaves Report

**410-bwaves**
Exec Time (%)	Codelet Name	Error (%)
99.8	shell_lam_shell__159	0.1
51.4	block_solver_mat_times_vec__167	23.0
1.4	flux_lam_flux__58	0.3
29.2	jacobian_lam_jacobian__30	0.6
56.4	block_solver_bi_cgstab_block__52	4.4
51.4	block_solver_mat_times_vec__167	23.0

Exec Time : 1.657044e+10 cycles

Matching : 99.8%

                
        subroutine bi_cgstab_block(x,b,a,axp,ayp,azp,axm,aym,azm,
     $  epsilon,nb,nx,ny,nz,maxit)
c       /* Bi-CGSTAB */
        implicit none
        integer nb,nx,ny,nz,kit, i,j,k,l,m
        integer maxit

        real*8 x(nb,nx,ny,nz),b(nb,nx,ny,nz),r(nb,nx,ny,nz),
     1  rhat(nb,nx,ny,nz),p(nb,nx,ny,nz),v(nb,nx,ny,nz),t(nb,nx,ny,nz)

        real*8 a(nb,nb,nx,ny,nz),
     1  axp(nb,nb,nx,ny,nz),ayp(nb,nb,nx,ny,nz),azp(nb,nb,nx,ny,nz),
     2  axm(nb,nb,nx,ny,nz),aym(nb,nb,nx,ny,nz),azm(nb,nb,nx,ny,nz)

        real*8 epsilon,r2,alpha,beta,rho,oldrho,omega,oldomega,tmp,tmp1


        kit = 0

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 do l=1,nb
                    r(l,i,j,k)=0.0d0
                 enddo
              enddo
           enddo
        enddo

        call mat_times_vec(r,x,a,axp,ayp,azp,axm,aym,azm,
     $  nb,nx,ny,nz)
        r2=0.0d0

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 do l=1,nb
                    r(l,i,j,k) = b(l,i,j,k) - r(l,i,j,k)
                    r2 =r2+r(l,i,j,k)**2
                    p(l,i,j,k) = 0.0d0
                    v(l,i,j,k) = 0.0d0
                    t(l,i,j,k)=0.0d0
                    rhat(l,i,j,k) = r(l,i,j,k)
                 enddo
              enddo
           enddo
        enddo

        oldrho = 1
        oldomega = 1
        alpha = 1

c       write(10,*) 'r2==', r2
        DO WHILE ( (r2 - epsilon). GT. 0.0d0 )
                kit = kit + 1
                rho=0.0d0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            rho =rho+rhat(l,i,j,k)*r(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                beta = (rho / oldrho)*(alpha/oldomega)

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            p(l,i,j,k) = r(l,i,j,k) + beta * 
     1                      ( p(l,i,j,k) - oldomega * v(l,i,j,k) )
                         enddo
                      enddo
                   enddo
                enddo
                call mat_times_vec(v,p,a,axp,ayp,azp,axm,aym,azm,
     $          nb,nx,ny,nz)
                tmp=0.0d0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            tmp=tmp+rhat(l,i,j,k)*v(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                alpha = rho / tmp

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            r(l,i,j,k)=r(l,i,j,k)-alpha * v(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                call mat_times_vec(t,r,a,axp,ayp,azp,axm,aym,azm,
     $          nb,nx,ny,nz)
                tmp=0.0d0
                tmp1=0.0d0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            tmp=tmp+t(l,i,j,k)*r(l,i,j,k)
                            tmp1=tmp1+t(l,i,j,k)*t(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                omega = tmp/tmp1

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            x(l,i,j,k)= x(l,i,j,k)+alpha*p(l,i,j,k)+
     1                      omega*r(l,i,j,k)
                            r(l,i,j,k)=r(l,i,j,k)-omega*t(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                oldrho = rho
                oldomega = omega
                r2=0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            r2=r2+r(l,i,j,k)*r(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo

c               print *,'iteration ',kit,'  |residual|^2 =',r2
        END DO
c       write(30,*) 'convergence after ',kit,' iterations.'
       write(10,*) '|residual|^2 =', r2
        maxit = maxit + kit

        RETURN
        END

        subroutine mat_times_vec(y,x,a,axp,ayp,azp,axm,aym,azm,
     $  nb,nx,ny,nz)
        implicit none
        integer nb,nx,ny,nz,i,j,k,m,l,kit,im1,ip1,jm1,jp1,km1,kp1

        real*8 y(nb,nx,ny,nz),x(nb,nx,ny,nz)

        real*8 a(nb,nb,nx,ny,nz),
     1  axp(nb,nb,nx,ny,nz),ayp(nb,nb,nx,ny,nz),azp(nb,nb,nx,ny,nz),
     2  axm(nb,nb,nx,ny,nz),aym(nb,nb,nx,ny,nz),azm(nb,nb,nx,ny,nz)


      do k=1,nz
         km1=mod(k+nz-2,nz)+1
         kp1=mod(k,nz)+1
         do j=1,ny
            jm1=mod(j+ny-2,ny)+1
            jp1=mod(j,ny)+1
            do i=1,nx
               im1=mod(i+nx-2,nx)+1
               ip1=mod(i,nx)+1
               do l=1,nb
                  y(l,i,j,k)=0.0d0
                  do m=1,nb
                     y(l,i,j,k)=y(l,i,j,k)+
     1               a(l,m,i,j,k)*x(m,i,j,k)+
     2               axp(l,m,i,j,k)*x(m,ip1,j,k)+
     3               ayp(l,m,i,j,k)*x(m,i,jp1,k)+
     4               azp(l,m,i,j,k)*x(m,i,j,kp1)+
     5               axm(l,m,i,j,k)*x(m,im1,j,k)+
     6               aym(l,m,i,j,k)*x(m,i,jm1,k)+
     7               azm(l,m,i,j,k)*x(m,i,j,km1)
                  enddo
               enddo
            enddo
         enddo
        enddo          

 

c        y=x
c        where (mask) y=tmp
        return
        end

invivoblock_solver_mat_times_vec__167

Real time: 8.521994e+09 -- Predicted time: 1.107973e+10

Callcount: 851

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
393	### 1	404.635172490719	1.309768e+07	1.022026e+07	21.9
793	### 2	399.496157122142	1.281759e+07	9.627160e+06	24.8
564	### 3	47.7306514716025	1.381433e+07	1.160517e+07	15.9

                      subroutine jacobian(q,je,jv,kt,kx,ky,kz,Re,Pr,gm,nx,ny,nz,
     $step,ax)

      implicit none
      integer nx,ny,nz,ax,ish,jsh,ksh,i,j,k,ip1,jp1,kp1
      real*8 kt,kx,ky,kz,gm,Re,Pr,al0,al1,al2,al3,al4,al5,al6,step

      real*8 q(5,nx,ny,nz)

      real*8 je(5,5,nx,ny,nz),jv(5,5,nx,ny,nz)

      real*8 u,v,w,fi2,alf,tht,a1,mu,ro,us,vs,ws,ros
C     Compute shift distance
      ish=0
      jsh=0
      ksh=0
      if (ax.eq.1) ish=1
      if (ax.eq.2) jsh=1
      if (ax.eq.3) ksh=1

      al0=(kx**2+ky**2+kz**2)/Pr
      al1=4.0d0/3.0d0*kx**2+ky**2+kz**2
      al2=kx*ky
      al3=1./3.0d0*kx*kz
      al4=kx**2+4.0d0/3.0d0*ky**2+kz**2
      al5=1.0d0/3.0d0*ky*kz
      al6=kx**2+ky**2+4.0d0/3.0d0*kz**2


      do k=1,nz
         kp1=mod(k,nz+1-ksh)+ksh
         do j=1,ny
            jp1=mod(j,ny+1-jsh)+jsh
            do i=1,nx
               ip1=mod(i,nx+1-ish)+ish

C     Initialize support variables
               ro=q(1,i,j,k)
               u=q(2,i,j,k)/ro
               v=q(3,i,j,k)/ro
               w=q(4,i,j,k)/ro
               fi2=0.5d0*(gm-1.0d0)*(u*u+v*v+w*w)
               alf=gm*q(5,i,j,k)/ro
               a1=alf-fi2
               tht=kx*u+ky*v+kz*w
               mu=((gm-1.0d0)*(q(5,i,j,k)/ro-0.5d0*
     1              (u*u+v*v+w*w)))**0.75d0

C     Initialize an jacobian for euiler part

               je(1,1,i,j,k)= kt
               je(2,1,i,j,k)= -u*tht+kx*fi2
               je(3,1,i,j,k)= -v*tht+ky*fi2
               je(4,1,i,j,k)= -w*tht+kz*fi2
               je(5,1,i,j,k)= tht*(2.0d0*fi2-alf)

               je(1,2,i,j,k)= kx
               je(2,2,i,j,k)= kt+tht-(gm-2.0d0)*kx*u
               je(3,2,i,j,k)= kx*v-(gm-1.0d0)*ky*u
               je(4,2,i,j,k)= kx*w-(gm-1.0d0)*kz*u
               je(5,2,i,j,k)= kx*a1-(gm-1.0d0)*u*tht

               je(1,3,i,j,k)= ky
               je(2,3,i,j,k)= ky*u-(gm-1.0d0)*kx*v
               je(3,3,i,j,k)= kt+tht-(gm-2.0d0)*ky*v
               je(4,3,i,j,k)= ky*w-(gm-1.0d0)*kz*v
               je(5,3,i,j,k)= ky*a1-(gm-1.0d0)*v*tht

               je(1,4,i,j,k)= kz
               je(2,4,i,j,k)= kz*u-(gm-1.0d0)*kx*w
               je(3,4,i,j,k)= kz*v-(gm-1.0d0)*ky*w
               je(4,4,i,j,k)= kt+tht-(gm-2.0d0)*kz*w
               je(5,4,i,j,k)= kz*a1-(gm-1.0d0)*w*tht

               je(1,5,i,j,k)= 0.0d0
               je(2,5,i,j,k)= (gm-1.0d0)*kx     
               je(3,5,i,j,k)= (gm-1.0d0)*ky
               je(4,5,i,j,k)= (gm-1.0d0)*kz
               je(5,5,i,j,k)= gm*tht+kt

C     Initialize an jacobian for viscous part


               ros=q(1,ip1,jp1,kp1)
               us=q(2,ip1,jp1,kp1)/ros
               vs=q(3,ip1,jp1,kp1)/ros
               ws=q(4,ip1,jp1,kp1)/ros
               a1=(1.0d0/ros-1.0d0/ro)/step               
               mu=(mu +((gm-1.0d0)*(q(5,ip1,jp1,kp1)/ros-
     1         0.5d0*(us*us+vs*vs+ws*ws)))**0.75d0)/2.0d0


               jv(1,1,i,j,k)=0.0d0
               jv(2,1,i,j,k)=mu/step*(
     $         al1*(u/ro-us/ros)+al2*(v/ro-vs/ros)+al3*(w/ro-ws/ros))
               jv(3,1,i,j,k)=mu/step*(
     $         al2*(u/ro-us/ros)+al4*(v/ro-vs/ros)+al5*(w/ro-ws/ros))
               jv(4,1,i,j,k)=mu/step*(
     $         al3*(u/ro-us/ros)+al5*(v/ro-vs/ros)+al6*(w/ro-ws/ros))
               jv(5,1,i,j,k)=mu/step*(
     $         al1*(u**2/ro-us**2/ros)+2.0d0*al2*(u*v/ro-us*vs/ros)+
     $         2.0d0*al3*(u*w/ro-us*ws/ros)+al4*(v**2/ro-vs**2/ros)+
     $         al6*(w**2/ro-ws**2/ros)+2.0d0*al5*(v*w/ro-vs*ws/ros)+
     $         al0*(q(5,i,j,k)/ro**2-q(5,ip1,jp1,kp1)/ros**2)+
     $         al0*((u*u+v*v+w*w)/ro-(us**2+vs**2+ws**2)/ros))

               jv(1,2,i,j,k)=0.0d0
               jv(2,2,i,j,k)=mu*al1*a1
               jv(3,2,i,j,k)=mu*al2*a1
               jv(4,2,i,j,k)=mu*al3*a1
               jv(5,2,i,j,k)=-jv(2,1,i,j,k)-mu*al0*(us/ros-u/ro)/step

               jv(1,3,i,j,k)=0.0d0
               jv(2,3,i,j,k)=mu*al2*a1
               jv(3,3,i,j,k)=mu*al4*a1
               jv(4,3,i,j,k)=mu*al5*a1
               jv(5,3,i,j,k)=-jv(3,1,i,j,k)-mu*al0*(vs/ros-v/ro)/step

               jv(1,4,i,j,k)=0.0d0
               jv(2,4,i,j,k)=mu*al3*a1
               jv(3,4,i,j,k)=mu*al5*a1
               jv(4,4,i,j,k)=mu*al6*a1
               jv(5,4,i,j,k)=-jv(4,1,i,j,k)-mu*al0*(ws/ros-w/ro)/step

               jv(1,5,i,j,k)=0.0d0
               jv(2,5,i,j,k)=0.0d0
               jv(3,5,i,j,k)=0.0d0
               jv(4,5,i,j,k)=0.0d0
               jv(5,5,i,j,k)=mu*al0*a1
            enddo
         enddo
      enddo
      return
      end

invivojacobian_lam_jacobian__30

Real time: 4.851621e+09 -- Predicted time: 4.818960e+09

Callcount: 75

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
67	### 1	58.5640267331168	1.667783e+07	1.707306e+07	2.3
8	### 2	3.00225848911494	3.147402e+08	3.153064e+08	0.1
11	### 3	3.00331428049316	3.441076e+08	3.443595e+08	0.0
14	### 4	3.00234592734919	2.707702e+08	2.710894e+08	0.1
17	### 5	3.0018995683394	2.113549e+08	2.121556e+08	0.3
20	### 6	3.00521274160711	9.259428e+07	9.296912e+07	0.4
23	### 7	3.0141263061733	4.583992e+07	4.619410e+07	0.7

                
        subroutine bi_cgstab_block(x,b,a,axp,ayp,azp,axm,aym,azm,
     $  epsilon,nb,nx,ny,nz,maxit)
c       /* Bi-CGSTAB */
        implicit none
        integer nb,nx,ny,nz,kit, i,j,k,l,m
        integer maxit

        real*8 x(nb,nx,ny,nz),b(nb,nx,ny,nz),r(nb,nx,ny,nz),
     1  rhat(nb,nx,ny,nz),p(nb,nx,ny,nz),v(nb,nx,ny,nz),t(nb,nx,ny,nz)

        real*8 a(nb,nb,nx,ny,nz),
     1  axp(nb,nb,nx,ny,nz),ayp(nb,nb,nx,ny,nz),azp(nb,nb,nx,ny,nz),
     2  axm(nb,nb,nx,ny,nz),aym(nb,nb,nx,ny,nz),azm(nb,nb,nx,ny,nz)

        real*8 epsilon,r2,alpha,beta,rho,oldrho,omega,oldomega,tmp,tmp1


        kit = 0

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 do l=1,nb
                    r(l,i,j,k)=0.0d0
                 enddo
              enddo
           enddo
        enddo

        call mat_times_vec(r,x,a,axp,ayp,azp,axm,aym,azm,
     $  nb,nx,ny,nz)
        r2=0.0d0

        do k=1,nz
           do j=1,ny
              do i=1,nx
                 do l=1,nb
                    r(l,i,j,k) = b(l,i,j,k) - r(l,i,j,k)
                    r2 =r2+r(l,i,j,k)**2
                    p(l,i,j,k) = 0.0d0
                    v(l,i,j,k) = 0.0d0
                    t(l,i,j,k)=0.0d0
                    rhat(l,i,j,k) = r(l,i,j,k)
                 enddo
              enddo
           enddo
        enddo

        oldrho = 1
        oldomega = 1
        alpha = 1

c       write(10,*) 'r2==', r2
        DO WHILE ( (r2 - epsilon). GT. 0.0d0 )
                kit = kit + 1
                rho=0.0d0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            rho =rho+rhat(l,i,j,k)*r(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                beta = (rho / oldrho)*(alpha/oldomega)

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            p(l,i,j,k) = r(l,i,j,k) + beta * 
     1                      ( p(l,i,j,k) - oldomega * v(l,i,j,k) )
                         enddo
                      enddo
                   enddo
                enddo
                call mat_times_vec(v,p,a,axp,ayp,azp,axm,aym,azm,
     $          nb,nx,ny,nz)
                tmp=0.0d0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            tmp=tmp+rhat(l,i,j,k)*v(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                alpha = rho / tmp

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            r(l,i,j,k)=r(l,i,j,k)-alpha * v(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                call mat_times_vec(t,r,a,axp,ayp,azp,axm,aym,azm,
     $          nb,nx,ny,nz)
                tmp=0.0d0
                tmp1=0.0d0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            tmp=tmp+t(l,i,j,k)*r(l,i,j,k)
                            tmp1=tmp1+t(l,i,j,k)*t(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                omega = tmp/tmp1

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            x(l,i,j,k)= x(l,i,j,k)+alpha*p(l,i,j,k)+
     1                      omega*r(l,i,j,k)
                            r(l,i,j,k)=r(l,i,j,k)-omega*t(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo
                oldrho = rho
                oldomega = omega
                r2=0

                do k=1,nz
                   do j=1,ny
                      do i=1,nx
                         do l=1,nb
                            r2=r2+r(l,i,j,k)*r(l,i,j,k)
                         enddo
                      enddo
                   enddo
                enddo

c               print *,'iteration ',kit,'  |residual|^2 =',r2
        END DO
c       write(30,*) 'convergence after ',kit,' iterations.'
       write(10,*) '|residual|^2 =', r2
        maxit = maxit + kit

        RETURN
        END

        subroutine mat_times_vec(y,x,a,axp,ayp,azp,axm,aym,azm,
     $  nb,nx,ny,nz)
        implicit none
        integer nb,nx,ny,nz,i,j,k,m,l,kit,im1,ip1,jm1,jp1,km1,kp1

        real*8 y(nb,nx,ny,nz),x(nb,nx,ny,nz)

        real*8 a(nb,nb,nx,ny,nz),
     1  axp(nb,nb,nx,ny,nz),ayp(nb,nb,nx,ny,nz),azp(nb,nb,nx,ny,nz),
     2  axm(nb,nb,nx,ny,nz),aym(nb,nb,nx,ny,nz),azm(nb,nb,nx,ny,nz)


      do k=1,nz
         km1=mod(k+nz-2,nz)+1
         kp1=mod(k,nz)+1
         do j=1,ny
            jm1=mod(j+ny-2,ny)+1
            jp1=mod(j,ny)+1
            do i=1,nx
               im1=mod(i+nx-2,nx)+1
               ip1=mod(i,nx)+1
               do l=1,nb
                  y(l,i,j,k)=0.0d0
                  do m=1,nb
                     y(l,i,j,k)=y(l,i,j,k)+
     1               a(l,m,i,j,k)*x(m,i,j,k)+
     2               axp(l,m,i,j,k)*x(m,ip1,j,k)+
     3               ayp(l,m,i,j,k)*x(m,i,jp1,k)+
     4               azp(l,m,i,j,k)*x(m,i,j,kp1)+
     5               axm(l,m,i,j,k)*x(m,im1,j,k)+
     6               aym(l,m,i,j,k)*x(m,i,jm1,k)+
     7               azm(l,m,i,j,k)*x(m,i,j,km1)
                  enddo
               enddo
            enddo
         enddo
        enddo          

 

c        y=x
c        where (mask) y=tmp
        return
        end

invivoblock_solver_bi_cgstab_block__52

Real time: 9.346048e+09 -- Predicted time: 9.781475e+09

Callcount: 25

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	3.09020560490978	8.348555e+07	9.862624e+07	15.3
5	### 2	4.288745255644	1.926579e+08	1.601882e+08	16.8
9	### 3	2.99999967970211	3.244180e+08	2.747442e+08	15.3
12	### 4	4.04350372155428	3.677495e+08	3.687185e+08	0.2
15	### 5	4.14890823579954	4.689913e+08	4.584233e+08	2.2
19	### 6	2.99843673850849	5.739560e+08	5.731351e+08	0.1
20	### 7	1.99802514917976	6.656611e+08	6.657009e+08	0.0
22	### 8	1.99711849508498	6.209930e+08	6.197234e+08	0.2

                      subroutine flux(q,e,f,g,ev,fv,gv,Re,Pr,gm,nx,ny,nz,
     $     dx,dy,dz)
      implicit none
      integer nx,ny,nz
      real*8 gm,Re,Pr,dx,dy,dz

      real*8 q(5,nx,ny,nz),e(5,nx,ny,nz),f(5,nx,ny,nz),g(5,nx,ny,nz),
     1     ev(5,nx,ny,nz),fv(5,nx,ny,nz),gv(5,nx,ny,nz)

      real*8 u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz),p(nx,ny,nz),
     1     ro(nx,ny,nz),
     2     mu(nx,ny,nz)

      real*8 t0,t1,t2,t3
      integer im1,im2,ip1,ip2,jm1,jm2,jp1,jp2,km1,km2,kp1,kp2,i,j,k,l
      real*8 dx2,dy2,dz2

      dx2=2.0d0*dx
      dy2=2.0d0*dy
      dz2=2.0d0*dz

      do k=1,nz
         do j=1,ny
            do i=1,nx

               ro(i,j,k)=q(1,i,j,k)
               u(i,j,k)=q(2,i,j,k)/ro(i,j,k)
               v(i,j,k)=q(3,i,j,k)/ro(i,j,k)
               w(i,j,k)=q(4,i,j,k)/ro(i,j,k)
               p(i,j,k)=(gm-1.0d0)*(q(5,i,j,k)-0.5d0*ro(i,j,k)*
     1         (u(i,j,k)**2+v(i,j,k)**2+w(i,j,k)**2))
               mu(i,j,k)=(gm*p(i,j,k)/ro(i,j,k))**0.75d0

C     Euler's fluxes
               e(1,i,j,k)=ro(i,j,k)*u(i,j,k)
               e(2,i,j,k)=ro(i,j,k)*u(i,j,k)*u(i,j,k)+p(i,j,k)
               e(3,i,j,k)=ro(i,j,k)*u(i,j,k)*v(i,j,k)
               e(4,i,j,k)=ro(i,j,k)*u(i,j,k)*w(i,j,k)
               e(5,i,j,k)=u(i,j,k)*(q(5,i,j,k)+p(i,j,k))


               f(1,i,j,k)=ro(i,j,k)*v(i,j,k)
               f(2,i,j,k)=ro(i,j,k)*v(i,j,k)*u(i,j,k)
               f(3,i,j,k)=ro(i,j,k)*v(i,j,k)*v(i,j,k)+p(i,j,k)
               f(4,i,j,k)=ro(i,j,k)*v(i,j,k)*w(i,j,k)
               f(5,i,j,k)=v(i,j,k)*(q(5,i,j,k)+p(i,j,k))

               g(1,i,j,k)=ro(i,j,k)*w(i,j,k)
               g(2,i,j,k)=ro(i,j,k)*w(i,j,k)*u(i,j,k)
               g(3,i,j,k)=ro(i,j,k)*w(i,j,k)*v(i,j,k)
               g(4,i,j,k)=ro(i,j,k)*w(i,j,k)*w(i,j,k)+p(i,j,k)
               g(5,i,j,k)=w(i,j,k)*(q(5,i,j,k)+p(i,j,k))
            enddo
         enddo
      enddo


      do k=1,nz
         km2=mod(k+nz-3,nz)+1
         km1=mod(k+nz-2,nz)+1
         kp1=mod(k,nz)+1
         kp2=mod(k+1,nz)+1
         do j=1,ny
            jm2=mod(j+ny-3,ny)+1
            jm1=mod(j+ny-2,ny)+1
            jp1=mod(j,ny)+1
            jp2=mod(j+1,ny)+1
            do i=1,nx
               im2=mod(i+nx-3,nx)+1
               im1=mod(i+nx-2,nx)+1
               ip1=mod(i,nx)+1
               ip2=mod(i+1,nx)+1      

         

C     Viscous fluxes
               ev(1,i,j,k)=0.0d0
               t0=0.5*(mu(i,j,k)+mu(ip1,j,k))
               t3=gm*p(i,j,k)/ro(i,j,k)

               t1=(v(i,jp1,k)-v(i,jm1,k))/dy2
               t2=(w(i,j,kp1)-w(i,j,km1))/dz2
               ev(2,i,j,k)=t0/3.0d0*(4.0d0*(u(ip1,j,k)-u(i,j,k))/dx-
     1              (t1+(v(ip1,jp1,k)-v(ip1,jm1,k))/dy2+
     2              t2+(w(ip1,j,kp1)-w(ip1,j,km1))/dz2))

               t1=(u(i,jp1,k)-u(i,jm1,k))/dy2
               ev(3,i,j,k)=t0*((t1+(u(ip1,jp1,k)
     1                        -u(ip1,jm1,k))/dy2)/2.0d0+
     2              (v(ip1,j,k)-v(i,j,k))/dx)

               t2=(u(i,j,kp1)-u(i,j,km1))/dz2
               ev(4,i,j,k)=t0*((t2+(u(ip1,j,kp1)
     1              -u(ip1,j,km1))/dz2)/2.0d0+
     2              (w(ip1,j,k)-w(i,j,k))/dx)

               ev(5,i,j,k)=0.5d0*((u(ip1,j,k)+u(i,j,k))*ev(2,i,j,k)+
     1              (v(ip1,j,k)+v(i,j,k))*ev(3,i,j,k)+
     2              (w(ip1,j,k)+w(i,j,k))*ev(4,i,j,k))+
     3              t0/Pr/(gm-1.0d0)*(gm*p(ip1,j,k)/ro(ip1,j,k)-t3)/dx

c     ************************************************************   
               fv(1,i,j,k)=0.0d0
               t0=(mu(i,j,k)+mu(i,jp1,k))/2.0d0

               t1=(v(ip1,j,k)-v(im1,j,k))/dx2
               fv(2,i,j,k)=t0*(((v(ip1,jp1,k)
     2              -v(im1,jp1,k))/dx2+t1)/2.0d0+
     1              (u(i,jp1,k)-u(i,j,k))/dy)

               t1=(u(ip1,j,k)-u(im1,j,k))/dx2
               t2=(w(i,j,kp1)-w(i,j,km1))/dz2
               fv(3,i,j,k)=t0/3.0d0*(4.0d0*(v(i,jp1,k)-v(i,j,k))/dy-
     1         ((u(ip1,jp1,k)-u(im1,jp1,k))/dx2+t1+
     2         (w(i,jp1,kp1)-w(i,jp1,km1))/dz2+t2))

               fv(4,i,j,k)=t0*(
     1         0.5d0*((u(ip1,jp1,k)-u(im1,jp1,k))/dx2+t1)+
     2         (w(i,jp1,k)-w(i,j,k))/dy)

               fv(5,i,j,k)=0.5*(
     1         (u(i,jp1,k)+u(i,j,k))*fv(2,i,j,k)+
     2         (v(i,jp1,k)+v(i,j,k))*fv(3,i,j,k)+
     3         (w(i,jp1,k)+w(i,j,k))*fv(4,i,j,k))+
     4         t0/Pr/(gm-1.0d0)*(gm*p(i,jp1,k)/ro(i,jp1,k)-t3)/dy

C     *************************************************************

               gv(1,i,j,k)=0.0d0
               t0=(mu(i,j,k)+mu(i,j,kp1))/2.0d0
               t1=(w(ip1,j,k)-w(im1,j,k))/dx2  
               gv(2,i,j,k)=t0*(
     1         ((w(ip1,j,kp1)-w(im1,j,kp1))/dx2+t1)/2.0d0+
     2         (u(i,j,kp1)-u(i,j,k))/dz)

               t1=(w(i,jp1,k)-w(i,jm1,k))/dy2
               gv(3,i,j,k)=t0*(
     1         ((w(i,jp1,kp1)-w(i,jm1,kp1))/dy2+t1)/2.0d0+
     2         (v(i,j,kp1)-v(i,j,k))/dz)

               t1=(u(ip1,j,k)-u(im1,j,k))/dx2
               t2=(v(i,jp1,k)-v(i,jm1,k))/dy2
               gv(4,i,j,k)=t0/3.0d0*(4.0d0*(w(i,j,kp1)-w(i,j,k))/dz-
     1         ((u(ip1,j,kp1)-u(im1,j,kp1))/dx2+t1+
     2          (v(i,jp1,kp1)-v(i,jm1,kp1))/dy2+t2))

                gv(5,i,j,k)=0.5d0*(
     1          (u(i,j,kp1)+u(i,j,k))*gv(2,i,j,k)+
     2          (v(i,j,kp1)+v(i,j,k))*gv(3,i,j,k)+
     4          (w(i,j,kp1)+w(i,j,k))*gv(4,i,j,k))+
     5          t0/Pr/(gm-1.0d0)*(gm*p(i,j,kp1)/ro(i,j,kp1)-t3)/dz

                

             enddo
          enddo
       enddo 


       
      return
      end

invivoflux_lam_flux__58

Real time: 2.336126e+08 -- Predicted time: 2.326995e+08

Callcount: 25

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	8.372797e+06	1.017800e+07	17.7
13	### 2	23.9173139456854	9.379260e+06	9.352484e+06	0.2

                      subroutine shell(Re,Pr,nx,ny,nz,
     $nuim,nuex2,nuex4,cfl,scheme,conf,ni,maxit) 
      implicit none
      COMMON /scale/ ro1,p1,ro2,p2
      integer nx,ny,nz,ni,n,i,j,k,scheme,conf,si,sj,sk,l,m
      integer im2,im1,ip1,ip2,jm2,jm1,jp1,jp2,km2,km1,kp1,kp2
      integer maxit
      real*8 nuim,nuex2,nuex4,cfl,epsilon,dx,dy,dz,dt,rad
      real*8 gm,u1,v1,w1,p1,ro1,u2,v2,w2,p2,ro2,Re,Pr,cfll
     $,time,t1,t2,t3,t4,t5,t6,t7,t8,mu
      real*8 dqnorm

      real*8 q(5,nx,ny,nz),dq(5,nx,ny,nz),rhs(5,nx,ny,nz),e(5,nx,ny,nz),
     1     f(5,nx,ny,nz),g(5,nx,ny,nz),ev(5,nx,ny,nz),fv(5,nx,ny,nz),
     2     gv(5,nx,ny,nz),diss(5,nx,ny,nz)

      real*8 a(5,5,nx,ny,nz),be(5,5,nx,ny,nz),ce(5,5,nx,ny,nz),
     1     av(5,5,nx,ny,nz),bv(5,5,nx,ny,nz),cv(5,5,nx,ny,nz)

      real*8 ae(5,5,nx,ny,nz),axp(5,5,nx,ny,nz),ayp(5,5,nx,ny,nz),
     1     azp(5,5,nx,ny,nz),axm(5,5,nx,ny,nz),aym(5,5,nx,ny,nz),
     2     azm(5,5,nx,ny,nz),ident(5,5,nx,ny,nz)


       real*8 u(nx,ny,nz),v(nx,ny,nz),w(nx,ny,nz),p(nx,ny,nz),
     1 ro(nx,ny,nz),at(nx,ny,nz)

      real*8 xyz(3,nx,ny,nz)


C       This particular problem is periodic only




C     Data initialization
      epsilon=1.d-16

      si=nx/8
      sj=ny/8
      sk=nz/8
      dx=1.0d0/(nx-1)
      dy=1.0d0/(ny-1)
      dz=1.0d0/(nz-1)
      time=0.0d0



        do k = 1,nz
           t3=(k-1)*dz
           do j=1,ny
              t2=(j-1)*dy
              do i=1,nx
                 t1=(i-1)*dx
                 xyz(1,i,j,k)=t1
                 xyz(2,i,j,k)=t2
                 xyz(3,i,j,k)=t3
                 do l=1,5
                    rhs(l,i,j,k)=0.0d0
                    do m=1,5
                       a(m,l,i,j,k)=0.0d0
                       axp(m,l,i,j,k)=0.0d0
                       ayp(m,l,i,j,k)=0.0d0
                       azp(m,l,i,j,k)=0.0d0
                       axm(m,l,i,j,k)=0.0d0
                       aym(m,l,i,j,k)=0.0d0
                       azm(m,l,i,j,k)=0.0d0
                       ident(m,l,i,j,k)=0.0d0
                    enddo
                    ident(l,l,i,j,k)=1.0d0
                    rhs(l,i,j,k)=0.0d0
                    dq(l,i,j,k)=0.0d0
                 enddo
              enddo
           enddo
        enddo              

C Flow Initialization

      gm=1.4d0
      u1=0.0d0
      v1=0.0d0
      w1=0.0d0
      p1=1.0d0
      ro1=1.0d0
      u2=0.0d0
      v2=0.0d0
      w2=0.0d0
      p2=0.1d0
      ro2=0.1d0

      t1=ro2
      t2=ro2*u2
      t3=ro2*v2
      t4=ro2*w2
      t5=p2/(gm-1.0d0)+0.5d0*ro2*(u2**2+v2**2+w2**2)


        do k = 1,nz
           do j=1,ny
              do i=1,nx
                 q(1,i,j,k)=t1
                 q(2,i,j,k)=t2
                 q(3,i,j,k)=t3
                 q(4,i,j,k)=t4
                 q(5,i,j,k)=t5
              enddo
           enddo
        enddo

      t1=ro1
      t2=ro1*u1
      t3=ro1*v1
      t4=ro1*w1
      t5=p1/(gm-1.0d0)+0.5d0*ro1*(u1**2+v1**2+w1**2)

      if (conf.EQ.0) then

         do k=nz/2-sk,nz/2+sk
            do j=ny/2-sj,ny/2+sj
               do i=nx/2-si,nx/2+si
                 q(1,i,j,k)=t1
                 q(2,i,j,k)=t2
                 q(3,i,j,k)=t3
                 q(4,i,j,k)=t4
                 q(5,i,j,k)=t5
              enddo
           enddo
        enddo
      else
         rad=DMIN1((si*dx)**2,(sj*dy)**2,(sk*dz)**2)

         do k=1,nz
            t8=((k-nz/2)*dz)**2
            do j=1,ny
               t7=((j-ny/2)*dy)**2
               do i=1,nx
                  t6=((i-nx/2)*dx)**2 +t7 + t8
                  if (t6.LE.rad) then
                     q(1,i,j,k)=t1
                     q(2,i,j,k)=t2
                     q(3,i,j,k)=t3
                     q(4,i,j,k)=t4
                     q(5,i,j,k)=t5
                  endif
               enddo
            enddo
         enddo
      endif

      

C Propagation in time



      time = 0.0d0

      do n=1,ni
C     Time step definition

      cfll=0.1d0+(n-1.0d0)*cfl/20.0d0
      if (cfll.ge.cfl) cfll=cfl
      t8=0.0d0

      do k=1,nz
         do j=1,ny
            do i=1,nx
               t1=q(1,i,j,k)
               t2=q(2,i,j,k)/t1
               t3=q(3,i,j,k)/t1
               t4=q(4,i,j,k)/t1
               t5=(gm-1.0d0)*(q(5,i,j,k)-0.5d0*t1*(t2*t2+t3*t3+t4*t4))
               t6=dSQRT(gm*t5/t1)
               mu=gm*Pr*(gm*t5/t1)**0.75d0*2.0d0/Re/t1
               t7=((dabs(t2)+t6)/dx+mu/dx**2)**2 +
     1            ((dabs(t3)+t6)/dy+mu/dy**2)**2 +
     2            ((dabs(t4)+t6)/dz+mu/dz**2)**2
               t7=DSQRT(t7)
               t8=max(t8,t7)
            enddo
         enddo
      enddo
      dt=cfll / t8


C Left hand side
      call jacobian(q,ae,av,0.0d0,1.0d0,0.0d0,0.0d0,Re,Pr,gm,
     1     nx,ny,nz,dx,1)
      call jacobian(q,be,bv,0.0d0,0.0d0,1.0d0,0.0d0,Re,Pr,gm,
     1     nx,ny,nz,dy,2)
      call jacobian(q,ce,cv,0.0d0,0.0d0,0.0d0,1.0d0,Re,Pr,gm,
     1     nx,ny,nz,dz,3)


      do k=1,nz
         km1=mod(k-2+nz,nz)+1
         kp1=mod(k,nz)+1
         do j=1,ny
            jm1=mod(j-2+ny,ny)+1
            jp1=mod(j,ny)+1
            do i=1,nx
               im1=mod(i-2+nx,nx)+1
               ip1=mod(i,nx)+1
               do l=1,5
                  do m=1,5
                     t1=ident(m,l,i,j,k)
                     a(m,l,i,j,k)=t1 - 0.5d0*dt*(
     1               (av(m,l,i,j,k)-av(m,l,im1,j,k))/dx+
     2               (bv(m,l,i,j,k)-bv(m,l,i,jm1,k))/dy+
     3               (cv(m,l,i,j,k)-cv(m,l,i,j,km1))/dz)/Re
     4               + 2.0d0*nuim*dt*(1.0d0/dx+1.0d0/dy+1.0d0/dz)*t1     

                     axp(m,l,i,j,k)=0.5d0*dt/dx*
     1               (ae(m,l,ip1,j,k)-av(m,l,i,j,k)/Re) 
     2               - nuim*dt/dx * t1
                     axm(m,l,i,j,k)=-0.5d0*dt/dx*
     1                    (ae(m,l,im1,j,k)-av(m,l,im1,j,k)/Re)
     2               - nuim*dt/dx * t1 

                     ayp(m,l,i,j,k)=0.5d0*dt/dy*
     1               (be(m,l,i,jp1,k)-bv(m,l,i,j,k)/Re)
     2               - nuim*dt/dy * t1
                     aym(m,l,i,j,k)=-0.5d0*dt/dy*
     1               (be(m,l,i,jm1,k)-bv(m,l,i,jm1,k)/Re)
     2               - nuim*dt/dy * t1

                     azp(m,l,i,j,k)=0.5d0*dt/dz*
     1               (ce(m,l,i,j,kp1)-cv(m,l,i,j,k)/Re)
     2               - nuim*dt/dz * t1
                     azm(m,l,i,j,k)=-0.5d0*dt/dz*
     1               (ce(m,l,i,j,km1)-cv(m,l,i,j,km1)/Re)
     2               - nuim*dt/dz * t1

                  enddo
               enddo
            enddo
         enddo
      enddo

      

C Right hand side
      call flux(q,e,f,g,ev,fv,gv,Re,Pr,gm,nx,ny,nz,dx,dy,dz)


      do k=1,nz
         km2=mod(k+nz-3,nz)+1
         km1=mod(k+nz-2,nz)+1
         kp1=mod(k,nz)+1
         kp2=mod(k+1,nz)+1
         do j=1,ny
            jm2=mod(j+ny-3,ny)+1
            jm1=mod(j+ny-2,ny)+1
            jp1=mod(j,ny)+1
            jp2=mod(j+1,ny)+1
            do i=1,nx
               im2=mod(i+nx-3,nx)+1
               im1=mod(i+nx-2,nx)+1
               ip1=mod(i,nx)+1
               ip2=mod(i+1,nx)+1
               do l=1,5
                  t1= -0.5d0*dt*(
     1            (e(l,ip1,j,k)-e(l,im1,j,k))/dx +
     2            (f(l,i,jp1,k)-f(l,i,jm1,k))/dy +
     3            (g(l,i,j,kp1)-g(l,i,j,km1))/dz) +
     4            dt/Re*((ev(l,i,j,k)-ev(l,im1,j,k))/dx +
     5                  (fv(l,i,j,k)-fv(l,i,jm1,k))/dy +
     6                  (gv(l,i,j,k)-gv(l,i,j,km1))/dz)

C     Artificial Viscosity
C     Explicit dissipation - second order
                  t2=dt*nuex2*(
     1            (q(l,ip1,j,k)-2.0d0*q(l,i,j,k)+q(l,im1,j,k))/dx +
     2            (q(l,i,jp1,k)-2.0d0*q(l,i,j,k)+q(l,i,jm1,k))/dy +
     3            (q(l,i,j,kp1)-2.0d0*q(l,i,j,k)+q(l,i,j,km1))/dz)

C     Explicit dissipation - fourth order
                  t2=t2-dt*nuex4*(
     1            (q(l,ip2,j,k)-4.0d0*q(l,ip1,j,k)+6.0d0*q(l,i,j,k)-
     2            4.0d0*q(l,im1,j,k)+q(l,im2,j,k))/dx +
     3            (q(l,i,jp2,k)-4.0d0*q(l,i,jp1,k)+6.0d0*q(l,i,j,k)-
     4            4.0d0*q(l,i,jm1,k)+q(l,i,jm2,k))/dy +
     5            (q(l,i,j,kp2)-4.0d0*q(l,i,j,kp1)+6.0d0*q(l,i,j,k)-
     6            4.0d0*q(l,i,j,km1)+q(l,i,j,km2))/dx )

                  rhs(l,i,j,k)=t1+t2
C     initial guess is rhs
                  dq(l,i,j,k)=rhs(l,i,j,k)
               enddo
            enddo
         enddo
      enddo
     
c      write(10,*) 'Time step: ',n,'  dt: ',dt 

      if (scheme.eq.1) call bi_cgstab_block(dq,rhs,a,axp,
     $ ayp,azp,axm,aym,azm,epsilon,5,nx,ny,nz,maxit)


      do k=1,nz
         do j=1,ny
            do i=1,nx
               do l=1,5
                  q(l,i,j,k)=q(l,i,j,k)+dq(l,i,j,k)
               enddo
            enddo
         enddo
      enddo

c      write(10,*) n, q(1,2,2,2)
c      write(10,*) q(2,2,2,2),q(3,2,2,2),q(4,2,2,2)
c      write(10,*) q(5,2,2,2)



      time=time+dt
      enddo

      dqnorm = 0.0d0
      do k=1,nz
         do j=1,ny
            do i=1,nx
               do l=1,5
                  dqnorm = dqnorm + dq(l,i,j,k)*dq(l,i,j,k)
               enddo
            enddo
         enddo
      enddo

      write(30,*) 'dqnorm ==', dqnorm 

      return
      end

invivoshell_lam_shell__159

Real time: 1.654916e+10 -- Predicted time: 1.651822e+10

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	1.651822e+10	1.659105e+10	0.4

410-bwaves

Exec Time : 1.657044e+10 cycles

Matching : 99.8%

__invivo__block_solver_mat_times_vec__167

Real time: 8.521994e+09 -- Predicted time: 1.107973e+10

Callcount: 851

__invivo__jacobian_lam_jacobian__30

Real time: 4.851621e+09 -- Predicted time: 4.818960e+09

Callcount: 75

__invivo__block_solver_bi_cgstab_block__52

Real time: 9.346048e+09 -- Predicted time: 9.781475e+09

Callcount: 25

__invivo__flux_lam_flux__58

Real time: 2.336126e+08 -- Predicted time: 2.326995e+08

Callcount: 25

__invivo__shell_lam_shell__159

Real time: 1.654916e+10 -- Predicted time: 1.651822e+10

Callcount: 1

invivoblock_solver_mat_times_vec__167

invivojacobian_lam_jacobian__30

invivoblock_solver_bi_cgstab_block__52

invivoflux_lam_flux__58

invivoshell_lam_shell__159