!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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
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)
!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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|