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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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,wpreiso,temp
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
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 (wpreiso(osolve%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc wpreiso 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
wpreiso=0.d0
temp=0.d0
pressure=0.d0
! interpolate velocity solution onto solve octree
do k=1+iproc,osolve%nnode,nproc
! call octree_interpolate_five &
! (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))
call octree_interpolate_six &
(6,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%wnodepreiso,wpreiso(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
osolve%u=0.d0
osolve%v=0.d0
osolve%w=0.d0
osolve%wpreiso=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 (wpreiso, osolve%wpreiso, 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 (wpreiso)
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%wnodepreiso','interpolate_ov...',size(ov%wnodepreiso),'dp',-1) ; deallocate (ov%wnodepreiso)
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%wnodepreiso(ov%nnodepreiso),stat=threadinfo%err)
call heap (threadinfo,'ov%wnodepreiso','interpolate_ov...',size(ov%wnodepreiso),'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)
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
!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%wnodepreiso=0.d0
!ov%pressure=0.d0
ov%temporary_nodal_pressure=0.d0
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|