compute_pressure.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              COMPUTE_PRESSURE    Oct. 2008                                   |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine compute_pressure (params,osolve,ov,vo,mat)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine is used to calculate the pressure
! It is an equivalent to build_system but in a stripped down version

! nleaves is the number of leaves/fe
! nface is the number of bad faces
! nnode is the number of nodes
! icon is the element-node connectivity matrix
! mpe is the number of nodes per elements
! x,y,z are the coordinates of the nodes of the fe grid
! lsf contains the nodal values of the nlsf level set functions
! mat is the material matrix for the nmat materials
! materialn contains the material number associated to each lsf
! u,v,w are the nodal velocities from the previous time step/iteration
! temp and pressure are the nodal temperature and pressure interpolated from the
!   previous time step/iteration
! vo is the structure containing the information on where the void is
! levelcut is the maximum number of levels used for dividing cut cells
! levelapprox is used to improve the estimate of the positive volume
!   in cut cells by another subdivision

! Like in build_system, we go through a routine called pressure_cut
! that deals with the integration of cut cells with precision
! It is necessary to go through this lengthy process of checking for cut cells
! in case the penalty parameter is not uniform in the fe grid

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions
!use mpi

implicit none

include 'mpif.h'

type (parameters) params
type (octreesolve) osolve
type (octreev) ov
type (void) vo
type (material) mat(0:params%nmat)

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer err,iproc,nproc,ierr,ileaves,level,i,k,levelmax,icut
double precision,dimension(:,:),allocatable::lsf_el
double precision r0,s0,t0,rst
double precision,dimension(:),allocatable::press
!integer,dimension(:),allocatable::levs
character*72  :: shift

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

shift=' '

allocate (press(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc press in compute_pressure$')

press=0.d0

do ileaves=1+iproc,osolve%nleaves,nproc
   ! first test if in the void
   if (vo%leaf(ileaves).eq.0) then
      level=0
      levelmax=params%levelcut
      allocate (lsf_el(params%mpe,osolve%nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsf_el in compute_pressure$')
         do k=1,params%mpe
            do i=1,osolve%nlsf
               lsf_el(k,i)=osolve%lsf(osolve%icon(k,ileaves),i)
            enddo
         enddo
      r0=-1.d0
      s0=-1.d0
      t0=-1.d0
      rst=2.d0
      call pressure_cut (params,level,levelmax,osolve%icon(1,ileaves),         &
                         osolve%x,osolve%y,osolve%z,mat,params%materialn,      &
                         ov%unode,ov%vnode,ov%wnode,osolve%temp,press(ileaves),&
                         osolve%nnode,lsf_el,osolve%nlsf,r0,s0,t0,rst,icut,    &
                         ileaves,osolve%eviscosity(ileaves))
      deallocate (lsf_el)
   endif
enddo

osolve%pressure=0.d0
call mpi_allreduce (press,osolve%pressure,osolve%nleaves,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)

deallocate (press)

if (iproc.eq.0 .and. params%debug>=1) then
   write(*,'(a,2e13.5)') shift//'raw pressures :',minval(osolve%pressure),maxval(osolve%pressure)
end if


!12/02/2008
! modification by Jean (no need to rescale pressures by element size
! as it is not done at element level (in make_pressure) any longer
!allocate (levs(nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc levs in compute_pressure$')
!call octree_find_element_level (octree,noctree,levs,nleaves)
!pressure=pressure*(8.d0**levs)
!deallocate (levs)

return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|