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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! COMPUTE_DIVERGENCE Oct. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine compute_divergence (params,osolve,ov,vo,istep,iter)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine calculates the divergence inside each element
! by using the velocity information at the nodes.
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
use gauss
implicit none
type(parameters) params
type (octreev) ov
type (octreesolve) osolve
type (void) vo
integer istep
integer iter
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision,dimension(:),allocatable :: x,y,z,vx,vy,vz
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
double precision,dimension(:),allocatable :: divergence,divergencetemp
double precision :: div,volume
integer err,ierr,nproc,iproc,k,i,iint,nint,mpe
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
nint=8
mpe=params%mpe
allocate (divergence(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc divergence in compute_divergence$')
allocate (divergencetemp(osolve%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc divergencetemp in compute_divergence$')
allocate (x(mpe),y(mpe),z(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc x,y,z in compute_divergence$')
allocate (vx(mpe),vy(mpe),vz(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vx,vy,vz in compute_divergence$')
allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=err) ; if (err.ne.0) call stop_run ('Error alloc dhdx,dhdy,dhdz in compute_divergence$')
divergence=0.d0
divergencetemp=0.d0
do i=1+iproc,osolve%nleaves,nproc
div=0.d0
if (vo%leaf(i).eq.0) then
do k=1,mpe
x(k) =osolve%x(osolve%icon(k,i))
y(k) =osolve%y(osolve%icon(k,i))
z(k) =osolve%z(osolve%icon(k,i))
vx(k)=ov%unode(osolve%icon(k,i))
vy(k)=ov%vnode(osolve%icon(k,i))
vz(k)=ov%wnode(osolve%icon(k,i))
enddo
do iint=1,nint
call compute_dhdx_dhdy_dhdz (mpe,rr(iint),ss(iint),tt(iint),x,y,z,dhdx,dhdy,dhdz,volume)
do k=1,mpe
div=div+(dhdx(k)*vx(k)+dhdy(k)*vy(k)+dhdz(k)*vz(k))*volume*ww(iint)
enddo
enddo
end if
divergencetemp(i)=div
end do
call mpi_reduce (divergencetemp, divergence, osolve%nleaves,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
if (iproc.eq.0) write (8,*) 'divergence ',minval(divergence),maxval(divergence)
call DoRuRe_div_stats (params%doDoRuRe,istep,iter,osolve%nleaves,divergence)
deallocate (divergence)
deallocate (divergencetemp)
deallocate (x,y,z)
deallocate (vx,vy,vz)
deallocate (dhdx,dhdy,dhdz)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|