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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|