Skip to content
Snippets Groups Projects
post.f90 44.3 KiB
Newer Older

MODULE definitions

   type structure
      integer,dimension(:),pointer::octree
      integer noctree,nnode,nleaves
      integer nelemr,nlsf 
      double precision,dimension(:),pointer::x,y,z
      double precision,dimension(:),pointer::u,v,w,wiso,temp,pressure,strain
      double precision,dimension(:),pointer::nodal_pressure,spressure
      double precision,dimension(:),pointer::e2d,eviscosity
      double precision,dimension(:),pointer::crit,dilatr
      integer,dimension(:,:),pointer::icon,iconr
      double precision,dimension(:,:),pointer:: lsf
      integer,dimension(:),pointer::on(:),is_plastic,matnum
      logical, dimension(:), pointer :: whole_leaf_in_fluid
   end type structure

   type sheet
      integer nsurface,nt
      double precision,dimension(:),pointer::r,s,x,y,z,xn,yn,zn,u,v,w
      integer,dimension(:,:),pointer::icon
      integer nhull,levelt,itype,material
   end type sheet

END MODULE definitions


!==============================================================================!
!==============================================================================!
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!==============================================================================!
!==============================================================================!
!                                                                              |
!              POST PROCESSING    Apr. 2008                                    |
!                                                                              |
!==============================================================================!
!==============================================================================!

program data_postprocessing

use definitions

implicit none

type (sheet),dimension(:),allocatable::surface
type(structure),SAVE :: ov
integer :: ie,mer,nn,ner,i,i1,i2,i3,i4,i5,kfx,kfy,kfz,kft
integer :: iunit,ilsf,nne,is,k,npcl,j,np,nl
integer :: output_u_field
integer :: output_v_field
integer :: output_w_field
integer :: output_velo_vect 
integer :: output_preiso_velo_vect 
Dave Whipp's avatar
Dave Whipp committed
integer :: output_iso_velo_vect 
integer :: output_disp_field
integer :: output_press_fieldn 
integer :: output_press_fielde 
integer :: output_smooth_press_fielde
integer :: output_countnode_field
integer :: output_temp_field
integer :: output_elem_matnum
integer :: output_e2d_fielde 
integer :: output_e2d_fieldn
integer :: output_eviscosity_fielde 
integer :: output_is_plastic_fielde
integer :: output_dilatr_fielde
integer :: output_crit_field
integer :: output_strain_field
integer :: output_lsf_field
integer :: output_surface_icon
integer :: output_strain_tensor
integer :: output_ps
integer :: output_cubes
integer :: output_rivers
integer :: output_regular
integer :: output_zisodisp
integer iproc,nproc,ierr,nnn,nnne,nnnep,levelmax,ix,iy,iz,nlinks,iordermin
integer iconmin,iconmax
integer,dimension(:),allocatable::nstrain,donor,order
integer,dimension(:),allocatable::levs,li
integer nz,ni,ijk,ptcnt,nb,dumpi
integer,dimension(:),allocatable::invoid,elvoid,rtf,ftr
integer myicon(100),nstep,ndir,iter,ii,lsf,levmax

logical,dimension(:),allocatable::influid,do_it,subset_leaves,instrain

character clsf*3,c4*4,cc4*4,dir*128

double precision :: eps,dil,current_time,activation_time,zmin,xx,yy,zz,maxe2d,dist
double precision, allocatable, dimension (:) :: ov_nodee2d, ov_nodecrit
double precision, allocatable, dimension (:) :: countnode 
double precision,dimension(:,:),allocatable::s,n1,n2,n3,zisodisp
double precision,dimension(:,:,:),allocatable::strain,strainn
double precision,dimension(:),allocatable::xi,yi,zi,ui,vi,wi,si,ei
double precision,dimension(:),allocatable::s11,s12,s13,s22,s23,s33,str11,str12,str13,str22,str23,str33
double precision,dimension(:),allocatable::azimuth1,azimuth3,dip1,dip3
double precision zmax,dz,l1,l2,l3,n11,n12,n13,n21,n22,n23,n31,n32,n33,con,dxy
double precision :: sselemx,sselemy,sselemz,xminls,yminls,zminls
double precision dxyz,x,y,z,xcut,dumpdp
integer icut 


!==============================================================================!
!==============================================================================!

call system('clear')

nne=8

print*,'Enter directory (default is OUT)'
read (*,'(a)') dir
ndir=len_trim(dir)
if (ndir.eq.0) then
  ndir=3
  dir(1:ndir)='OUT'
endif

print*,'Enter time step >'
read*,nstep
write (c4,'(i4)') nstep
if (nstep.lt.1000) c4(1:1)='0'
if (nstep.lt.100) c4(1:2)='00'
if (nstep.lt.10) c4(1:3)='000'

do i=1,1
  write (*,*) 'Is this a nested model? (y/n)'
  read (*,'(a)') nestin
  select case (trim(nestin))
  case ('y','Y','yes','Yes')
    nest=.true.
    do j=1,1
      write (*,*) 'Should the values be scaled to fit the large(st) scale model? (y/n)'
      read (*,'(a)') nestsclin
      select case (trim(nestsclin))
      case ('y','Y','yes','Yes')
        nestscl=.true.
      case ('n','N','no','No')
        nestscl=.false.
      case default 
        write (*,*) 'Error: Please input either "y" or "n"'
        j=0
      end select
    enddo
  case ('n','N','no','No')
    nest=.false.
  case default 
    write (*,*) 'Error: Please input either "y" or "n"'
    i=0
  end select
enddo

print*,'Is it a debug file ? (y or n)'
read (*,'(a)') cs
select case (cs)
case ('y')
   print*,'Which debug file: please enter grid iteration number'
   read*,iter
   write (cc4,'(i4)') iter
   if (iter.lt.1000) cc4(1:1)='0'
   if (iter.lt.100) cc4(1:2)='00'
   if (iter.lt.10) cc4(1:3)='000'
case ('n')
   print*,'it is not a debug file'
case default 
   print *,'error: please input y or n'
   stop
end select

!==============================================================================
!==============================================================================

write(*,*) '**************************************************************************'
write(*,*) '--------------------------------------------------------------------------'
write(*,*) '--------------------------------------------------------------------------'
write(*,*) '------------------- VTK file production in progress ----------------------'
write(*,*) '--------------------------------------------------------------------------'
write(*,*) '--------------------------------------------------------------------------'
if (cs=='y') then
   write(*,*) 'output file to be processed: ', '../'//dir(1:ndir)//'/time_'//c4//'_'//cc4//'.bin'
   open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'_'//cc4//'.bin',status='old',form='unformatted')
else
   write(*,*) 'output file to be processed: ', '../'//dir(1:ndir)//'/time_'//c4//'.bin'
   open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'.bin',status='old',form='unformatted')
end if
write(*,*) '--------------------------------------------------------------------------'
open(unit=77,file='input_of_outputs.txt',status='old')
read(77,*) output_u_field 
read(77,*) output_v_field 
read(77,*) output_w_field
read(77,*) output_velo_vect 
read(77,*) output_preiso_velo_vect 
Dave Whipp's avatar
Dave Whipp committed
read(77,*) output_iso_velo_vect 
read(77,*) output_disp_field 
read(77,*) output_press_fieldn 
read(77,*) output_press_fielde 
read(77,*) output_smooth_press_fielde
read(77,*) output_countnode_field
read(77,*) output_temp_field
read(77,*) output_elem_matnum
read(77,*) output_e2d_fielde 
read(77,*) output_e2d_fieldn
read(77,*) output_eviscosity_fielde 
read(77,*) output_is_plastic_fielde
read(77,*) output_dilatr_fielde
read(77,*) output_crit_field 
read(77,*) output_strain_field 
read(77,*) output_lsf_field 
read(77,*) output_surface_icon
read(77,*) output_strain_tensor
read(77,*) output_ps
read(77,*) output_cubes
read(77,*) output_rivers
read(77,*) output_regular
read(77,*) output_zisodisp

write(*,*) 'output u      field       ->',(output_u_field==1)
write(*,*) 'output v      field       ->',(output_v_field==1)
write(*,*) 'output w      field       ->',(output_w_field==1)
write(*,*) 'output velo   vectors     ->',(output_velo_vect==1) 
write(*,*) 'output preiso velo vectors->',(output_preiso_velo_vect==1) 
Dave Whipp's avatar
Dave Whipp committed
write(*,*) 'output iso velo vectors->',(output_preiso_velo_vect==1) 
write(*,*) 'output disp   field       ->',(output_disp_field==1)
write(*,*) 'output node press field   ->',(output_press_fieldn==1) 
write(*,*) 'output raw press field    ->',(output_press_fielde==1) 
write(*,*) 'output smooth press field ->',(output_smooth_press_fielde==1)
write(*,*) 'output countnode field    ->',(output_countnode_field==1)
write(*,*) 'output temp   field       ->',(output_temp_field==1)
write(*,*) 'output elem material #    ->',(output_elem_matnum==1)
write(*,*) 'output elemental e2d field->',(output_e2d_fielde==1) 
write(*,*) 'output nodal e2d field    ->',(output_e2d_fieldn==1)
write(*,*) 'output elem eff viscosity ->',(output_eviscosity_fielde==1) 
write(*,*) 'output element is plastic ->',(output_is_plastic_fielde==1)
write(*,*) 'output element dilat rate ->',(output_dilatr_fielde==1) 
write(*,*) 'output crit   field       ->',(output_crit_field==1)
write(*,*) 'output strain field       ->',(output_strain_field==1)
write(*,*) 'output lsf    field       ->',(output_lsf_field==1)
write(*,*) 'output surfaces icon      ->',(output_surface_icon==1)
write(*,*) 'output strain tensor      ->',(output_strain_tensor==1)
write(*,*) 'output principal stresses ->',(output_ps==1)
write(*,*) 'output elemental fields   ->',(output_cubes==1)
write(*,*) 'output rivers             ->',(output_rivers==1)
write(*,*) 'output regular            ->',(output_regular==1)
write(*,*) 'output zisodisp           ->',(output_zisodisp==1)
write(*,*) '--------------------------------------------------------------------------'

!==============================================================================!
!==============================================================================!

read (7) ov%noctree,ov%nnode,ov%nleaves,ner,ov%nlsf,npcl,current_time
nn=ov%nnode
allocate(ov%on(nn))
allocate(ov%x(nn),ov%y(nn),ov%z(nn))
allocate(ov%u(nn),ov%v(nn),ov%w(nn),ov%wiso(nn))
allocate(ov%pressure(ov%nleaves),ov%strain(nn))
allocate(ov%e2d(ov%nleaves),ov%eviscosity(ov%nleaves),ov%is_plastic(ov%nleaves))
allocate(ov%dilatr(ov%nleaves),ov%whole_leaf_in_fluid(ov%nleaves))
Dave Whipp's avatar
Dave Whipp committed
allocate(ov%crit(ov%nleaves),ov%matnum(ov%nleaves))
allocate(ov%temp(nn))
allocate(ov%lsf(nn,ov%nlsf))
allocate(influid(nn))
allocate(invoid(nn))
allocate(rtf(nn),ftr(nn),elvoid(nn))
allocate(surface(ov%nlsf))
allocate(ov%nodal_pressure(nn),ov%spressure(ov%nleaves))
allocate(ov%icon(8,ov%nleaves))
allocate(ov%octree(ov%noctree))

!========================
!=====[nodal values]=====
!========================

read(7)(ov%x(i),ov%y(i),ov%z(i),   &
        ov%u(i),ov%v(i),ov%w(i),   &
        ov%wiso(i),                &
        ov%lsf(i,1:ov%nlsf),       &
        ov%temp(i),                &
        ov%nodal_pressure(i),      &
        ov%strain(i),              &
        kfx,kfy,kfz,kft,i=1,nn)

!======================
!=====[icon array]=====
!======================

read(7) (ov%icon(1:8,ie),            &
         ov%pressure(ie),            & 
         ov%spressure(ie),           &
         ov%crit(ie),                &
         ov%e2d(ie),                 &
         ov%eviscosity(ie),          &
         ov%is_plastic(ie),          &
         ov%dilatr(ie),              &
         ov%whole_leaf_in_fluid(ie), & 
         ie=1,ov%nleaves)

!==============================
!=====[octree information]=====
!==============================

read(7) (ov%octree(i),i=1,ov%noctree)

!================================
!=====[bad face information]=====
!================================
mer=ner   ! mer=1
if(ner.gt.mer)mer=ner
allocate(ov%iconr(9,mer))
read(7) (ov%iconr(1:9,ie),ie=1,ner)

ov%nelemr=ner

!============================
!=====[void information]=====
!============================

read (7) (invoid(i),elvoid(i),ftr(i),rtf(i),influid(i),i=1,nn)
ov%on=0
do i=1,nn
   if (influid(i)) ov%on(i)=1
enddo

read(7) (i5,i=1,ner)

!==============================================================
!=====[measure work, strain, compute principal directions]=====
!==============================================================

allocate (s (3,ov%nleaves))
allocate (n1(3,ov%nleaves))
allocate (n2(3,ov%nleaves))
allocate (n3(3,ov%nleaves))
allocate (strain(3,3,ov%nleaves))

call find_volumes (ov%icon,ov%octree,ov%noctree,ov%nleaves,ov%lsf,ov%nlsf,ov%nnode,influid, &
                   ov%x,ov%y,ov%z,ov%u,ov%v,ov%w,s,n1,n2,n3,strain)

! rescaling of the values for plotting purposes

s=s/maxval(abs(s))*.01

strain=strain/maxval(abs(strain))*.01

!===============================
!=====[surface information]=====
!===============================

do is=1,ov%nlsf
   read(7) surface(is)%nsurface,activation_time,surface(is)%nt
   allocate (surface(is)%x(surface(is)%nsurface),  &
             surface(is)%y(surface(is)%nsurface),  &
             surface(is)%z(surface(is)%nsurface),  &
             surface(is)%xn(surface(is)%nsurface), &
             surface(is)%yn(surface(is)%nsurface), &
             surface(is)%zn(surface(is)%nsurface), &
             surface(is)%r(surface(is)%nsurface),  &
             surface(is)%s(surface(is)%nsurface),  &
             surface(is)%u(surface(is)%nsurface),  &
             surface(is)%v(surface(is)%nsurface),  &
             surface(is)%w(surface(is)%nsurface),  &
             surface(is)%icon(3,surface(is)%nt))
  
if (cs=='y') then
   read (7) (surface(is)%r(i),surface(is)%s(i), &
               surface(is)%x(i),surface(is)%y(i),surface(is)%z(i), &
               surface(is)%xn(i),surface(is)%yn(i),surface(is)%zn(i), &
               i=1,surface(is)%nsurface)
else
   read (7) (surface(is)%r(i),surface(is)%s(i), &
               surface(is)%x(i),surface(is)%y(i),surface(is)%z(i), &
               surface(is)%xn(i),surface(is)%yn(i),surface(is)%zn(i), &
               surface(is)%u(i),surface(is)%v(i),surface(is)%w(i), &
               i=1,surface(is)%nsurface)
end if

   read (7) (surface(is)%icon(1:3,i),i=1,surface(is)%nt)

end do

!=============================
!=====[cloud information]=====
!=============================

read (7) (dumpdp,      & 
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpdp,    & 
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpdp,    &
            dumpi,     &
            i=1,npcl)
!================================
!=====[isostasy information]=====
!================================

read (7) nb
dxy=1.d0/real(nb)
allocate(zisodisp(nb+1,nb+1))
read (7) ((zisodisp(i,j),j=1,nb+1),i=1,nb+1)
!============================
!=====[nest information]=====
!============================

read(7) sselemx,sselemy,sselemz,xminls,yminls,zminls

close(7)

!==============================================================================
!==============================================================================
!=====[C -> N]=================================================================
!==============================================================================
!==============================================================================

allocate(countnode(ov%nnode))
allocate(ov_nodee2d(ov%nnode))
allocate(ov_nodecrit(ov%nnode))

ov_nodee2d       = 0.d0
ov_nodecrit      = 0.d0
countnode        = 0.d0
 
do i=1,ov%nleaves
   do j=1,8
      k=ov%icon(j,i)
      ov_nodee2d(k)       = ov_nodee2d(k)      + ov%e2d(i)    
      ov_nodecrit(k)      = ov_nodecrit(k)     + ov%crit(i)    
      countnode(k)        = countnode(k)       + 1                     
   end do
end do

do i=1,ov%nnode
   if (countnode(i).gt.0.d0) then
      ov_nodee2d(i)      = ov_nodee2d(i)      / countnode(i)
      ov_nodecrit(i)     = ov_nodecrit(i)     / countnode(i)
   end if
end do

! Convert coordinates if using a nest
if (nest) then
  ov%x=ov%x*sselemx+xminls
  ov%y=ov%y*sselemy+yminls
  ov%z=ov%z*sselemz+zminls
endif

!==============================================================================
!==============================================================================

write(*,'(a,i10)') 'octree length              ',ov%noctree
write(*,'(a,i10)') 'nb of nodes                ',ov%nnode
write(*,'(a,i10)') 'nb of leaves               ',ov%nleaves
write(*,'(a,i10)') 'nb of bad faces            ',ner
write(*,'(a,i10)') 'nb of lsf                  ',ov%nlsf
write(*,*) '--------------------------------------------------------------------------'
write(*,'(a,2f30.20)') 'x                 range : ',minval(ov%x(1:nn)),       maxval(ov%x(1:nn))
write(*,'(a,2f30.20)') 'y                 range : ',minval(ov%y(1:nn)),       maxval(ov%y(1:nn))
write(*,'(a,2f30.20)') 'z                 range : ',minval(ov%z(1:nn)),       maxval(ov%z(1:nn))
write(*,'(a,2f30.20)') 'u                 range : ',minval(ov%u(1:nn)),       maxval(ov%u(1:nn))
write(*,'(a,2f30.20)') 'v                 range : ',minval(ov%v(1:nn)),       maxval(ov%v(1:nn))
write(*,'(a,2f30.20)') 'w                 range : ',minval(ov%w(1:nn)),       maxval(ov%w(1:nn))
write(*,'(a,2f30.20)') 'w (iso)           range : ',minval(ov%wiso(1:nn)),    maxval(ov%wiso(1:nn))
write(*,'(a,2f30.20)') 'pressure          range : ',minval(ov%pressure),      maxval(ov%pressure) 
write(*,'(a,2f30.20)') 'smooth pressure   range : ',minval(ov%spressure),     maxval(ov%spressure)
write(*,'(a,2f30.20)') 'pressure (nodes)  range : ',minval(ov%nodal_pressure),maxval(ov%nodal_pressure)
write(*,'(a,2f30.20)') 'e2d               range : ',minval(ov%e2d),           maxval(ov%e2d)
write(*,'(a,2f30.20)') 'e2d (nodes)       range : ',minval(ov_nodee2d),       maxval(ov_nodee2d)
write(*,'(a,2f30.20)') 'dilatation rate   range : ',minval(ov%dilatr),        maxval(ov%dilatr)
write(*,'(a,2f30.20)') 'crit              range : ',minval(ov%crit),          maxval(ov%crit)
write(*,'(a,2f30.20)') 'crit (nodes)      range : ',minval(ov_nodecrit),      maxval(ov_nodecrit)
write(*,'(a,2f30.20)') 'strain            range : ',minval(ov%strain(1:nn)),  maxval(ov%strain(1:nn))
Loading
Loading full blame...