correct node fix and multi-block pastr

Author: fangjian

pastr/src/pastr_data_convert.F90

@@ -163,8 +163,7 @@ end function dat_out_cal_2d
     function dat_out_cal_3d(data_in,name_in,name_out,x,y,z) result(data_out)
 
       use pastr_commvar, only : im,jm,km,const2
-      use pastr_flowvis, only : schlieren_xy
-      use pastr_gradients, only : grad_xy
+      use pastr_gradients, only : grad_3d
       use pastr_thermo_phys, only : sos
 
       real(wp),intent(in) :: data_in(:,:,:,:)
@@ -191,8 +190,107 @@ function dat_out_cal_3d(data_in,name_in,name_out,x,y,z) result(data_out)
 
         enddo
 
+
+        if(name_out(i)=='lci') then
+  
+            allocate(is(3))
+            is=0
+            do j=1,nvar_in
+              if(name_in(j)=='u1') is(1)=j
+              if(name_in(j)=='u2') is(2)=j
+              if(name_in(j)=='u3') is(3)=j
+            enddo
+  
+            if(.not. allocated(du1)) then
+              allocate(du1(1:3,0:im,0:jm,0:km),du2(1:3,0:im,0:jm,0:km),du3(1:3,0:im,0:jm,0:km))
+              du1=grad_3d(data_in(:,:,:,is(1)),x,y,z)
+              du2=grad_3d(data_in(:,:,:,is(2)))
+              du3=grad_3d(data_in(:,:,:,is(3)))
+            endif
+            
+  
+            data_out(:,:,:,i)=lambdacical(du1,du2,du3,im,jm,km)
+  
+            deallocate(is)
+  
+        endif
+
       enddo
 
+
     end function dat_out_cal_3d
 
+    function lambdacical(du,dv,dw,im,jm,km)
+
+    use pastr_utility, only: cube_root
+
+    real(wp) :: lambdacical(0:im,0:jm,0:km)
+    !
+    integer,intent(in) :: im,jm,km
+    real(wp),intent(in) :: du(1:3,0:im,0:jm,0:km),                      &
+                              dv(1:3,0:im,0:jm,0:km),                      &
+                              dw(1:3,0:im,0:jm,0:km)
+    !
+    integer :: i,j,k
+    real(wp) :: del,a11,a12,a13,a21,a22,a23,a31,a32,a33,Q,R,var1,var2,  &
+               var3,delta,lambmax,duref
+    real(wp),save :: lambref=0._wp
+    
+    duref=1.d0
+
+    lambmax=0._wp
+    !$omp parallel default(shared) private(i,j,k,del,a11,a12,a13,a21,  &
+    !$omp    a22,a23,a31,a32,a33,q,r,var1,var2,var3,delta)
+    !
+    !$omp do
+    do k=0,km
+    do j=0,jm
+    do i=0,im
+      !
+      del=-(du(1,i,j,k)/duref+dv(2,i,j,k)/duref+dw(3,i,j,k)/duref)/3._wp
+      !
+      a11=du(1,i,j,k)/duref+del; a12=du(2,i,j,k)/duref;     a13=du(3,i,j,k)/duref
+      a21=dv(1,i,j,k)/duref;     a22=dv(2,i,j,k)/duref+del; a23=dv(3,i,j,k)/duref
+      a31=dw(1,i,j,k)/duref;     a32=dw(2,i,j,k)/duref;     a33=dw(3,i,j,k)/duref+del
+
+      Q=-0.5_wp*(a11*a11+a12*a21+a13*a31+ a21*a12+a22*a22+a23*a32+      &
+                    a31*a13+a32*a23+a33*a33 )
+      R=-1._wp/3._wp*(a11*a11*a11+a11*a12*a21+a11*a13*a31+               &
+                            a12*a21*a11+a12*a22*a21+a12*a23*a31+               &
+                            a13*a31*a11+a13*a32*a21+a13*a33*a31+               &
+                            a21*a11*a12+a21*a12*a22+a21*a13*a32+               &
+                            a22*a21*a12+a22*a22*a22+a22*a23*a32+               &
+                            a23*a31*a12+a23*a32*a22+a23*a33*a32+               &
+                            a31*a11*a13+a31*a12*a23+a31*a13*a33+               &
+                            a32*a21*a13+a32*a22*a23+a32*a23*a33+               &
+                            a33*a31*a13+a33*a32*a23+a33*a33*a33                )
+      !
+      delta=Q**3/27._wp+R**2/4._wp
+      !
+      ! To solve the eqation: lambda**3+Q*lambda+R=0
+      if(delta>=tiny(0._wp) ) then
+        var1=cube_root(-0.5_wp*R+sqrt(delta))
+        var2=cube_root(-0.5_wp*R-sqrt(delta))
+        var3=0.5_wp*sqrt(3._wp)*(var1-var2)
+        lambdacical(i,j,k)=var3*var3
+        
+      else
+        lambdacical(i,j,k)=0._wp
+      end if
+
+      lambmax=max(lambmax,lambdacical(i,j,k))
+
+
+    end do
+    end do
+    end do
+    !$omp end do
+    !$omp end parallel
+
+    print*,' ** max lambdaci=',lambmax
+    !
+    print*,' ** Swirling strength λ calculated'
+    !
+  end function lambdacical
+
 end module pastr_data_convert

pastr/src/pastr_field_view.F90

@@ -258,7 +258,8 @@ subroutine write_3d_field(filein,fileout,nfirst,nlast,format)
 
           call H5ReadArray(r8dat,im,jm,km,trim(nam_var_in(m)),trim(file2read))
 
-          dat_var_in(:,:,:,m)=real(r8dat)
+          dat_var_in(:,:,:,m)=r8dat
+
         enddo
 
         dat_var_out=dat_out_cal(dat_var_in,nam_var_in,nam_var_out,x,y,z)

pastr/src/pastr_gradients.F90

@@ -132,13 +132,13 @@ subroutine gridjacobian_3d(x,y,z,ddi,ddj,ddk)
       !$OMP DO
       do j=0,jm
       do i=0,im
-        dx(3,i,j,:)=dfdk(x(i,j,:))
-        dy(3,i,j,:)=dfdk(y(i,j,:))
-        dz(3,i,j,:)=dfdk(z(i,j,:))
+        dx(3,i,j,:)=dfdk(x(i,j,:),geom=.true.)
+        dy(3,i,j,:)=dfdk(y(i,j,:),geom=.true.)
+        dz(3,i,j,:)=dfdk(z(i,j,:),geom=.true.)
       end do
       end do
       !$OMP END DO
-      !
+
       !$OMP DO
       do k=0,km
       do j=0,jm
@@ -177,7 +177,7 @@ subroutine gridjacobian_3d(x,y,z,ddi,ddj,ddk)
       !$OMP END DO
       !
       !$OMP END PARALLEL
-      !
+
       deallocate(dx,dy,dz)
       !
       print*, ' ** Grid Jacobian matrix is calculated'
@@ -316,12 +316,25 @@ function dfdj(var)
 
     end function dfdj
     !
-    function dfdk(var)
+    function dfdk(var,geom)
 
       real(wp) ,allocatable :: dfdk(:)
       real(wp),intent(in) :: var(0:km)
+      logical,intent(in),optional :: geom
+
+      logical :: lgeom
+
+      if(present(geom)) then
+        lgeom=geom
+      else
+        lgeom=.false.
+      endif
 
-      dfdk =diff6ec(var,km,lkhomo)
+      if(lgeom) then
+        dfdk =diff6ec_geom(var,km,lkhomo)
+      else
+        dfdk =diff6ec(var,km,lkhomo)
+      endif
 
     end function dfdk
 
@@ -387,5 +400,59 @@ end function diff6ec
     !+-------------------------------------------------------------------+
     !| The end of the function diff6ec.                                  |
     !+-------------------------------------------------------------------+
+    function diff6ec_geom(vin,dim,homo) result(vout)
+      
+      integer,intent(in) :: dim
+      logical,intent(in) :: homo
+      real(8),intent(in) :: vin(0:dim)
+      real(8) :: vout(0:dim)
+      
+      ! local data
+      integer :: i
+
+      if(homo) then
+
+        vout(0)=0.75d0 *((vin(1)-vin(0))-(vin(dim-1)-vin(dim)))-           &
+                0.15d0 *((vin(2)-vin(0))-(vin(dim-2)-vin(dim)))+           &
+                num1d60*((vin(3)-vin(0))-(vin(dim-3)-vin(dim)))     
+        vout(1)=0.75d0 *((vin(2)-vin(0))-(vin(0)-vin(0)))    -           &
+                0.15d0 *((vin(3)-vin(0))-(vin(dim-1)-vin(dim)))+           &
+                num1d60*((vin(4)-vin(0))-(vin(dim-2)-vin(dim)))       
+        vout(2)=0.75d0 *((vin(3)-vin(0))-(vin(1)-vin(0)))    -           &
+                0.15d0 *((vin(4)-vin(0))-(vin(0)-vin(0)))    +           &
+                num1d60*((vin(5)-vin(0))-(vin(dim-1)-vin(dim)))
+
+        vout(dim-2) =0.75d0 *(vin(dim-1)-vin(dim-3))-  &
+                     0.15d0 *(vin(dim)  -vin(dim-4))+  &
+                     num1d60*(vin(1)-vin(0)    -(vin(dim-5)-vin(dim)))
+        vout(dim-1)=0.75d0 *(vin(dim)-vin(dim-2))-     &
+                    0.15d0 *(vin(1)-vin(0)  +vin(dim)-vin(dim-3))+     &
+                    num1d60*(vin(2)-vin(0)  +vin(dim)-vin(dim-4))
+        vout(dim)  =0.75d0 *(vin(1)-vin(0)  +vin(dim)-vin(dim-1))-       &
+                    0.15d0 *(vin(2)-vin(0)  +vin(dim)-vin(dim-2))+       &
+                    num1d60*(vin(3)-vin(0)  +vin(dim)-vin(dim-3))
+      else
+
+        vout(0)=-0.5d0*vin(2)+2.d0*vin(1)-1.5d0*vin(0)
+        vout(1)=0.5d0*(vin(2)-vin(0))
+        vout(2)=num2d3*(vin(3)-vin(1))-num1d12*(vin(4)-vin(0))
+
+        vout(dim-2) =num2d3*(vin(dim-1)-vin(dim-3))-      &
+                     num1d12*(vin(dim)  -vin(dim-4))
+        vout(dim-1)=0.5d0*(vin(dim)-vin(dim-2))
+        vout(dim)  =0.5d0*vin(dim-2)-2.d0*vin(dim-1)+1.5d0*vin(dim)
+
+      endif
+
+      do i=3,dim-3
+        vout(i)  =0.75d0 *(vin(i+1)-vin(i-1))-            &
+                  0.15d0 *(vin(i+2)-vin(i-2))+            &
+                  num1d60*(vin(i+3)-vin(i-3))
+      enddo
+
+    end function diff6ec_geom
+    !+-------------------------------------------------------------------+
+    !| The end of the function diff6ec.                                  |
+    !+-------------------------------------------------------------------+
 
 end module pastr_gradients