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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! INTERPOLATE_OV_ON_OSOLVE Apr. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine interpolate_ov_on_osolve (osolve,ov,iter,params,threadinfo)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine transfers the velocity solution from the velocity octree onto
! the osolve octree. It also performs a transfer of the pressure from the nodes
! of the ov octree to the leaves of the solution octree.
! At the end, the velocity octree is redimensioned to fit the dimension of the
! solve octree and is thus ready for the next solution step
! osolve is the solution octree
! ov is the velocity octree
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use threads
use definitions
implicit none
type (parameters) params
type (octreesolve) osolve
type (octreev) ov
type (thread) threadinfo
integer iter
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer err,ie_ov,ie_level,ie_loc,nproc,ierr,iproc,k,ie,ie_osolve,i,j
double precision :: x,y,z,x0_leaf,y0_leaf,z0_leaf,dxyz_leaf
double precision,dimension(:),allocatable::pressure,pressurep
double precision,dimension(:),allocatable::u,v,w,temp
double precision,dimension(:),allocatable::countnode
double precision,dimension(:),allocatable::ov_nodepressure
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=' '
!=====C -> N=====================================
!allocate(countnode(ov%nnode))
!allocate(ov_nodepressure(ov%nnode))
!
!ov_nodepressure =0.d0
!countnode =0.d0
!
!do i=1,ov%nleaves
! do j=1,8
! k=ov%icon(j,i)
! ov_nodepressure(k) = ov_nodepressure(k) + ov%pressure(i)
! countnode(k) = countnode(k) + 1
! end do
!end do
!
!do i=1,ov%nnode
! if (countnode(i).gt.0.d0) ov_nodepressure(i) =ov_nodepressure(i) / countnode(i)
!end do
!
!deallocate(countnode)
!=====interpolate ov on osolve===================
if (iproc.eq.0) write(8,*) 'opla1' ; call flush(8)
allocate (pressure(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc pressure in main$')
allocate (pressurep(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc pressurep in main$')
allocate (u(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc u in main$')
allocate (v(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc v in main$')
allocate (w(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc w in main$')
allocate (temp(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc temp in main$')
u=0.d0
v=0.d0
w=0.d0
temp=0.d0
pressure=0.d0
! interpolate velocity solution onto solve octree
do k=1+iproc,osolve%nnode,nproc
call octree_interpolate_many &
(5,ov%octree,ov%noctree,ov%icon, &
ov%nleaves,ov%nnode, &
osolve%x(k),osolve%y(k),osolve%z(k), &
ov%unode,u(k), &
ov%vnode,v(k), &
ov%wnode,w(k), &
ov%temp,temp(k), &
ov%temporary_nodal_pressure,pressure(k))
enddo
!do k=1+iproc,osolve%nnode,nproc
! call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%unode ,ov%x,ov%y,ov%z,ov%nnode,osolve%x(k),osolve%y(k),osolve%z(k),u (k))
! call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%vnode ,ov%x,ov%y,ov%z,ov%nnode,osolve%x(k),osolve%y(k),osolve%z(k),v (k))
! call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%wnode ,ov%x,ov%y,ov%z,ov%nnode,osolve%x(k),osolve%y(k),osolve%z(k),w (k))
! call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%temp ,ov%x,ov%y,ov%z,ov%nnode,osolve%x(k),osolve%y(k),osolve%z(k),temp (k))
! call octree_interpolate3 (ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%temporary_nodal_pressure,ov%x,ov%y,ov%z,ov%nnode,osolve%x(k),osolve%y(k),osolve%z(k),pressure(k))
!enddo
if (params%debug>=1 .and. iproc.eq.0) then
write(*,'(a,2e11.4)') shift//'ov%pressure: ',minval(ov%temporary_nodal_pressure),maxval(ov%temporary_nodal_pressure)
endif
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
osolve%u=0.d0
osolve%v=0.d0
osolve%w=0.d0
osolve%temp=0.d0
pressurep=0.d0
!do i=1,osolve%nnode
! if (ieee_is_nan(u(i))) print *,'field -> NaN in u$'
! if (ieee_is_nan(v(i))) print *,'field -> NaN in v$'
! if (ieee_is_nan(w(i))) print *,'field -> NaN in w$'
! if (ieee_is_nan(temp(i))) print *,'field -> NaN in temp$'
! if (ieee_is_nan(pressure(i))) print *,'field -> NaN in pressure$'
!end do
call mpi_allreduce (u, osolve%u, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (v, osolve%v, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (w, osolve%w, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (temp, osolve%temp, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (pressure, pressurep, osolve%nnode, mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
!if (iproc.eq.0) print *,'int.ov.on.osolve: min/max pressurep',minval(pressurep),maxval(pressurep)
!=====N -> C=====================================
! further interpolation from nodes to leaves for pressure
do ie=1,osolve%nleaves
osolve%pressure(ie)=sum(pressurep(osolve%icon(1:8,ie)))/8.d0
enddo
osolve%spressure=osolve%pressure
if (params%debug>=1 .and. iproc.eq.0) then
write(*,'(a,2e11.4)') shift//'osolve%pressure: ',minval(osolve%pressure),maxval(osolve%pressure)
endif
deallocate (pressure)
deallocate (pressurep)
deallocate (u)
deallocate (v)
deallocate (w)
deallocate (temp)
! prepare the velo octree to receive the solution and carry
! it to the next step/iteration
call heap (threadinfo,'ov%octree','interpolate_ov...',size(ov%octree),'int',-1) ; deallocate (ov%octree)
call heap (threadinfo,'ov%x','interpolate_ov...',size(ov%x),'dp',-1) ; deallocate (ov%x)
call heap (threadinfo,'ov%y','interpolate_ov...',size(ov%y),'dp',-1) ; deallocate (ov%y)
call heap (threadinfo,'ov%z','interpolate_ov...',size(ov%z),'dp',-1) ; deallocate (ov%z)
call heap (threadinfo,'ov%unode','interpolate_ov...',size(ov%unode),'dp',-1) ; deallocate (ov%unode)
call heap (threadinfo,'ov%vnode','interpolate_ov...',size(ov%vnode),'dp',-1) ; deallocate (ov%vnode)
call heap (threadinfo,'ov%wnode','interpolate_ov...',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode)
call heap (threadinfo,'ov%temp','interpolate_ov...',size(ov%temp),'dp',-1) ; deallocate (ov%temp)
call heap (threadinfo,'ov%icon','interpolate_ov...',size(ov%icon),'int',-1) ; deallocate (ov%icon)
call heap (threadinfo,'ov%temporary...','interpolate_ov...',size(ov%temporary_nodal_pressure),'dp',-1)
deallocate (ov%temporary_nodal_pressure)
call heap (threadinfo,'ov%whole_leaf...','interpolate_ov...',size(ov%whole_leaf_in_fluid),'bool',-1)
deallocate (ov%whole_leaf_in_fluid)
ov%noctree=osolve%noctree
ov%nleaves=osolve%nleaves
ov%nnode=osolve%nnode
allocate (ov%octree(ov%noctree),stat=threadinfo%err)
call heap (threadinfo,'ov%octree','interpolate_ov...',size(ov%octree),'int',+1)
allocate (ov%icon(8,ov%nleaves),stat=threadinfo%err)
call heap (threadinfo,'ov%icon','interpolate_ov...',size(ov%icon),'int',+1)
allocate (ov%x(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%x','interpolate_ov...',size(ov%x),'dp',+1)
allocate (ov%y(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%y','interpolate_ov...',size(ov%y),'dp',+1)
allocate (ov%z(ov%nnode),stat=threadinfo%err) ; call heap (threadinfo,'ov%z','interpolate_ov...',size(ov%z),'dp',+1)
allocate (ov%unode(ov%nnode),stat=threadinfo%err)
call heap (threadinfo,'ov%unode','interpolate_ov...',size(ov%unode),'dp',+1)
allocate (ov%vnode(ov%nnode),stat=threadinfo%err)
call heap (threadinfo,'ov%vnode','interpolate_ov...',size(ov%vnode),'dp',+1)
allocate (ov%wnode(ov%nnode),stat=threadinfo%err)
call heap (threadinfo,'ov%wnode','interpolate_ov...',size(ov%wnode),'dp',+1)
allocate (ov%temp(ov%nnode),stat=threadinfo%err)
call heap (threadinfo,'ov%temp','interpolate_ov...',size(ov%temp),'dp',+1)
allocate (ov%temporary_nodal_pressure(ov%nnode),stat=threadinfo%err)
call heap (threadinfo,'ov%temporary_nodal_pressure','interpolate_ov...',size(ov%temporary_nodal_pressure),'dp',+1)
allocate (ov%whole_leaf_in_fluid(ov%nleaves),stat=threadinfo%err)
call heap (threadinfo,'ov%whole_leaf_in_fluid','interpolate_ov...',size(ov%whole_leaf_in_fluid),'bool',+1)
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
!allocate (ov%octree(ov%noctree),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%octree in main$')
!allocate (ov%icon(8,ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%icon in main$')
!allocate (ov%x(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%x in main$')
!allocate (ov%y(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%y in main$')
!allocate (ov%z(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%z in main$')
!allocate (ov%unode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%unode in main$')
!allocate (ov%vnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%vnode in main$')
!allocate (ov%wnode(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wnode in main$')
!allocate (ov%temp(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%temp in main$')
!allocate (ov%temporary_nodal_pressure(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%temp_nodal_pressure in main$')
!allocate (ov%whole_leaf_in_fluid(ov%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%whole_leaf_in_fluid in main$')
if (iproc.eq.0) write(8,*) 'opla9' ; call flush(8)
ov%icon=osolve%icon
ov%octree=osolve%octree
! beware that the sizes of osolve%x,y and z are not equal to osolve%nnode
! they are in fact overdimensioned because we do not know, a priori, the
! number of nodes in an octree; hence the explicit index length
ov%x(1:ov%nnode)=osolve%x(1:osolve%nnode)
ov%y(1:ov%nnode)=osolve%y(1:osolve%nnode)
ov%z(1:ov%nnode)=osolve%z(1:osolve%nnode)
ov%temp=osolve%temp
ov%unode=0.d0
ov%vnode=0.d0
ov%wnode=0.d0
!ov%pressure=0.d0
ov%temporary_nodal_pressure=0.d0
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|