Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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, pressure and strain are the nodal temperature, pressure and strain
! interpolated from the previous time step/iteration or, in the case of
! the strain, from the 3D cloud
! 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
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,dimension(:),allocatable::press
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%strain,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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|