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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
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
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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 chekcing for cut cells
! in case the penalty parameter is not uniform in the fe grid
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
implicit none
type (parameters) params
type (octreesolve) osolve
type (octreev) ov
type (void) vo
type (material) mat(0:params%nmat)
!type(parameters) params
!integer nface
!integer nnode
!integer icon(params%mpe,nleaves)
!double precision x(nnode),y(nnode),z(nnode)
!double precision lsf(nnode,nlsf)
!integer nlsf
!double precision u(nnode),v(nnode),w(nnode)
!double precision temp(nnode)
!double precision pressure(nleaves)
!double precision strain(nnode)
!type (void) vo
!type (octreesolve) octree
!integer noctree
!double precision eviscosity(nleaves)
!------------------------------------------------------------------------------|
!(((((((((((((((( 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,volc
double precision,dimension(:),allocatable::press
integer,dimension(:),allocatable::levs
character*72 :: shift
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|