Skip to content
Snippets Groups Projects
Select Git revision
  • feature/mumps
  • develop
  • feature/parabolic_subduction
  • feature/hdf5out
  • feature/optimisedMoveCloud
  • feature/thrust_isobc
  • feature/particleAdvectionFast
  • feature/orography_sinclair
  • feature/mixed_cloudmat
  • feature/improved_cloud_material_numbering
  • feature/multibox2
  • feature/EGU2020
  • feature/mumps_final
  • feature/oroprecip_tuebingen
  • MfN_Berlin
  • feature/cloudmat-sphere
  • feature/baselevel
  • feature/LK_manuscript_2018
  • feature/sb04_oroprecip
  • feature/oroprecip2
  • v0.5.1
  • v0.5
  • v0.4.3.0
  • v0.4.2.0
  • v0.4.1.0
  • v0.4.0.2
  • v0.4.0.1
  • v0.4
  • v0.4rc1
  • DOUAR-WSMP_v0.1.1
  • DOUAR-WSMP_v0.1
  • 2009-09-02
  • fromDavidWhipp
33 results
Blame Permalink
interpolate_ov_on_osolve.f90 12.30 KiB 
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              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

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)


!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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|