1 subroutine vortGLB (y, x, shp, shgl, ilwork, vortG) 2 3 use pointer_data 4 include "common.h" 5 include "mpif.h" 6c 7 real*8 y(nshg,ndof), x(numnp,nsd) 8c 9 real*8 shp(MAXTOP,maxsh,MAXQPT), 10 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 11 12 integer ilwork(nlwork) 13 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 14 real*8 vortG(nshg,4) ! The first three components are vorticity vector, the 4th one is Q criterion 15 real*8 vortIntG(nshg,4) 16 real*8 lpmassG(nshg) 17 real*8 tmp(nshg) 18 19 vortG=0 20 lpmassG=0 21 vortIntG=0 22 23 do iblk = 1, nelblk 24 nenl = lcblk(5,iblk) ! no. of vertices per element 25 iel = lcblk(1,iblk) 26 lelCat = lcblk(2,iblk) 27 lcsyst = lcblk(3,iblk) 28 iorder = lcblk(4,iblk) 29 nenl = lcblk(5,iblk) ! no. of vertices per element 30 nshl = lcblk(10,iblk) 31 mattyp = lcblk(7,iblk) 32 ndofl = lcblk(8,iblk) 33 nsymdl = lcblk(9,iblk) 34 npro = lcblk(1,iblk+1) - iel 35 inum = iel + npro - 1 36 ngauss = nint(lcsyst) 37c 38 allocate (tmpshp(nshl,MAXQPT)) 39 allocate (tmpshgl(nsd,nshl,MAXQPT)) 40 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 41 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 42 43 if(myrank.eq.0)then 44 if(mod(iblk,100).eq.0)then 45 write(*,*)'Compute vortIntG:', real(iblk)/nelblk 46 endif 47 endif 48 call vortELE(y,x,tmpshp,tmpshgl,mien(iblk)%p,vortIntG,lpmassG) 49 50 deallocate(tmpshp) 51 deallocate(tmpshgl) 52 enddo 53 54 if(myrank.eq.0)write(*,*)'Communicating...' 55 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 56 call commu(lpmassG,ilwork,1,'in ') 57 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 58 do i=1,4 59 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 60 tmp(:)=vortIntG(:,i) 61 call commu(tmp,ilwork,1,'in ') 62 vortIntG(:,i)=tmp(:) 63 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 64 enddo 65 call MPI_BARRIER(MPI_COMM_WORLD,ierr) 66 67 do i=1,4 68 vortG(:,i)=vortIntG(:,i)/lpmassG(:) 69 enddo 70 71! call write_field(myrank,'a'//char(0),'vortG'//char(0),5, 72! & vortG, 'd'//char(0), nshg,4,lstep) 73 74 return 75 end 76ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 77 78 subroutine vortELE(y,x,shp1,shgl1,ien,vortIntG,lpmassG) 79 80 81 include "common.h" 82 83 real*8 y(nshg,ndof), x(numnp,nsd) 84 real*8 shp1(nshl,MAXQPT), shgl1(nsd,nshl,MAXQPT) 85 real*8 shp2(nshl,ngauss), shgl2(nsd,nshl,ngauss) 86 real*8 shp(npro,nshl), shgl(npro,nsd,nshl) 87 real*8 sgn(npro,nshl) 88 integer ien(npro,nshl) 89 real*8 yl(npro,nshl,ndof) 90 real*8 xl(npro,nenl,nsd) 91 real*8 dxidx(npro,nsd,nsd) 92 real*8 WdetJ(npro) 93 real*8 vortIntL(npro,nshl,4) 94 real*8 lpmass(npro,nshl) 95 real*8 volum(npro) 96 real*8 alph2(npro) 97 real*8 vortIntG(nshg,4) 98 real*8 lpmassG(nshg) 99 real*8 shg(npro,nshl,nsd) 100 real*8 g1yi(npro,nflow) 101 real*8 g2yi(npro,nflow) 102 real*8 g3yi(npro,nflow) 103 104 call getsgn(ien,sgn) 105 shp2(:,1:ngauss)=shp1(:,1:ngauss) 106 shgl2(:,:,1:ngauss)=shgl1(:,:,1:ngauss) 107 108 call localy (y, yl, ien, ndof, 'gather ') 109 call localx (x, xl, ien, nsd, 'gather ') 110 111 112 vortIntL=0 113 lpmass=0 114 volum=0 115 alph2=0 116 117 118 do intp=1,ngauss 119 120 call getshp (shp2, shgl2, sgn, shp, shgl ) 121 122 call e3metric( xl, shgl, dxidx, 123 & shg, WdetJ) 124 125 g1yi = zero 126 g2yi = zero 127 g3yi = zero 128 129 130 do n=1,nshl 131 g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * yl(:,n,2) ! du/dx 132 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * yl(:,n,2) ! du/dy 133 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * yl(:,n,2) ! du/dz 134 135 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * yl(:,n,3) ! dv/dx 136 g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * yl(:,n,3) ! dv/dy 137 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * yl(:,n,3) ! dv/dz 138 139 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * yl(:,n,4) ! dw/dx 140 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * yl(:,n,4) ! dw/dy 141 g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * yl(:,n,4) ! dw/dz 142 143 alph2(:) = alph2(:) + shp(:,n)*shp(:,n)*WdetJ(:) 144 145 enddo 146 volum(:) = volum(:) + WdetJ(:) 147 148c--------- compute vorticity and Q criterion 149 do n=1,nshl 150 151 lpmass(:,n) = lpmass(:,n) + WdetJ(:)*shp(:,n)*shp(:,n) 152 153 vortIntL(:,n,1) = vortIntL(:,n,1) + 154 & (g2yi(:,4)-g3yi(:,3)) !dw/dy-dv/dz 155 & *WdetJ(:)*shp(:,n) 156 157 vortIntL(:,n,2) = vortIntL(:,n,2) + 158 & (g3yi(:,2)-g1yi(:,4)) !du/dz-dw/dx 159 & *WdetJ(:)*shp(:,n) 160 161 vortIntL(:,n,3) = vortIntL(:,n,3) + 162 & (g1yi(:,3)-g2yi(:,2)) !dv/dx-du/dy 163 & *WdetJ(:)*shp(:,n) 164 165 166 vortIntL(:,n,4) = vortIntL(:,n,4) + 167 & ( 0.5*(g2yi(:,4)-g3yi(:,3))**2 !dw/dy-dv/dz 168 & +0.5*(g3yi(:,2)-g1yi(:,4))**2 !du/dz-dw/dx 169 & +0.5*(g1yi(:,3)-g2yi(:,2))**2 !dv/dx-du/dy 170 & -g1yi(:,2)**2-g2yi(:,3)**2-g3yi(:,4)**2 171 & -0.5*(g2yi(:,4)+g3yi(:,3))**2 172 & -0.5*(g3yi(:,2)+g1yi(:,4))**2 173 & -0.5*(g1yi(:,3)+g2yi(:,2))**2 )*0.5 174 & *WdetJ(:)*shp(:,n) 175 176 enddo 177 178 enddo ! end of loop over Gaussian points 179 180 do n=1,nshl 181 lpmass(:,n)=lpmass(:,n)*volum(:)/alph2(:) 182 enddo 183 184 call local(vortIntG,vortIntL,ien,4,'scatter ') 185 call local(lpmassG,lpmass,ien,1,'scatter ') 186 187 return 188 end 189 190ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 191 192 subroutine PRcalcDuct (y) 193 194 include "common.h" 195 include "mpif.h" 196 197 real*8 y(nshg,ndof) 198 real*8 Mach(nshg),PresRec(nshg),PresCoeff(nshg),velocity(nshg,3) 199 real*8 Ptotal_contra, Pstatic_throat 200 open(911,file='PresConst.inp') 201 read(911,*) Ptotal_contra 202 read(911,*) Pstatic_throat 203 close(911) 204 write(*,*)'Ptotal_contra=',Ptotal_contra 205 write(*,*)'Pstatic_throat=',Pstatic_throat 206c in y the variables are u v w p T 207 Mach(:)=sqrt(y(:,1)**2+y(:,2)**2+y(:,3)**2)/sqrt(1.4*287*y(:,5)) 208 PresRec(:)=y(:,4)*(1+0.2*Mach(:)**2)**3.5/Ptotal_contra 209 PresCoeff(:)=(y(:,4)-Pstatic_throat)/(Ptotal_contra-Pstatic_throat) 210 velocity(:,1:3)=y(:,1:3) 211 write(*,*)'writing Mach, PR and Cp' 212 call write_field(myrank,'a'//char(0),'M'//char(0),1,Mach, 213 & 'd'//char(0),nshg,1,lstep) 214 call write_field(myrank,'a'//char(0),'PR'//char(0),2,PresRec, 215 & 'd'//char(0),nshg,1,lstep) 216 call write_field(myrank,'a'//char(0),'Cp'//char(0),2,PresCoeff, 217 & 'd'//char(0),nshg,1,lstep) 218 call write_field(myrank,'a'//char(0),'V'//char(0),1,velocity, 219 & 'd'//char(0),nshg,3,lstep) 220 end 221