!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! IMPROVE_OSOLVE Feb. 2007 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine improve_osolve (osolve,ov,params,threadinfo,istep,iter,total,step, &
inc,refine_level,boxes,cube_faces, &
increase_current_level)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine calculates the second invariant of the strain rate, and a value
! for the criterion used by the improve_octree subroutine, inside each element
! by using the velocity information at the nodes. Then, the routine improves the
! osolve on which the solution is calculated. Also, it refines the octree if the
! user has defined a box (or several boxes) in which the octree is artificially
! refined to a desired level. boxes and nboxes are read from the input file.
! The routine then refines user supplied rectangular areas of each of the six
! faces of the cubic simulation domain. The information are stored in cube_faces
! read from the input file. Finally, the octree is smoothed and super-smoothed
! if the user has set ismooth to 1
!
! ov is the velocity octree containing the velocity solution
! osolve is the octree to be improved
! refine_ratio is a parameter read in input.txt. When multiplied by the
! maximum of the criterion previously computed, one obtains the threshold
! used to determine whether a leaf is to be refined or not.
! refine_level is the level at which the osolve octree is to be refined
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
use gauss
use invariants
!use mpi
use threads
implicit none
include 'mpif.h'
type (octreesolve) osolve
type (octreev) ov
type (parameters) params
type (thread) threadinfo
integer istep,iter
double precision total,step,inc
integer refine_level
type (box) boxes(params%nboxes)
type (face), dimension(6) :: cube_faces
logical increase_current_level
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision jcb(3,3),jcbi(3,3)
double precision,dimension(:),allocatable :: x,y,z,vx,vy,vz
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
double precision,dimension(:),allocatable :: crit,crittemp,xc,yc,zc
double precision :: eps,r,s,t,volume,w
double precision :: exx,eyy,ezz,exy,eyz,ezx
double precision critmax,xxx,yyy,zzz,b,l
double precision :: x0,y0,z0,x1,y1,z1,E2d
integer err,ierr,nproc,iproc,k,i,iint,nint,mpe
integer :: j,nn,nb_of_boxes,level_box,nx,ny,nz,level,nv,nh
character*72 shift
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
eps=tiny(eps)
nint=8
mpe=8
shift=' '
!------------------------------------------------------------------
! first compute the refinement criterion measured on the ov octree
!------------------------------------------------------------------
call show_time (total,step,inc,1,'compute ref. criterion on ov$')
allocate (x(mpe),y(mpe),z(mpe),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'x/y/z','improve_osolve',size(x)+size(y)+size(z),'dp',+1)
allocate (vx(mpe),vy(mpe),vz(mpe),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'vx/y/z','improve_osolve',size(vx)+size(vy)+size(vz),'dp',+1)
allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',+1)
allocate (crit(ov%nleaves),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'crit','improve_osolve',size(crit),'dp',+1)
allocate (crittemp(ov%nleaves),stat=threadinfo%err)
if (params%debug.gt.1) call heap (threadinfo,'crittemp','improve_osolve',size(crittemp),'dp',+1)
crit=0.d0
crittemp=0.d0
r=0.d0
s=0.d0
t=0.d0
do i=1+iproc,ov%nleaves,nproc
do k=1,mpe
x(k)=ov%x(ov%icon(k,i))
y(k)=ov%y(ov%icon(k,i))
z(k)=ov%z(ov%icon(k,i))
vx(k)=ov%unode(ov%icon(k,i))
vy(k)=ov%vnode(ov%icon(k,i))
vz(k)=ov%wnode(ov%icon(k,i))
enddo
call compute_dhdx_dhdy_dhdz (mpe,r,s,t,x,y,z,dhdx,dhdy,dhdz,volume)
exx=0.d0
eyy=0.d0
ezz=0.d0
exy=0.d0
eyz=0.d0
ezx=0.d0
do k=1,mpe
exx=exx+dhdx(k)*vx(k)
eyy=eyy+dhdy(k)*vy(k)
ezz=ezz+dhdz(k)*vz(k)
exy=exy+(dhdx(k)*vy(k)+dhdy(k)*vx(k))/2.d0
eyz=eyz+(dhdy(k)*vz(k)+dhdz(k)*vy(k))/2.d0
ezx=ezx+(dhdz(k)*vx(k)+dhdx(k)*vz(k))/2.d0
enddo
call deviatoric (params%invariants_2d,exx,eyy,ezz,exy,eyz,ezx)
E2d=second_invariant(exx,eyy,ezz,exy,eyz,ezx)
select case (params%refine_criterion)
case(1)
crittemp(i)=sqrt(E2d)
! crittemp(i)=log(1.d0+sqrt(E2d))
case(2)
crittemp(i)=sqrt(exx**2+eyy**2+ezz**2)
case(3)
crittemp(i)=sqrt(E2d)*abs(x(1)-x(2))
case default
crittemp(i)=0.d0
end select
end do
call mpi_allreduce (crittemp,crit,ov%nleaves,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
if (params%debug.gt.1) call heap (threadinfo,'x/y/z','improve_osolve',size(x)+size(y)+size(z),'dp',-1)
deallocate (x,y,z)
if (params%debug.gt.1) call heap (threadinfo,'vx/vy/vz','improve_osolve',size(vx)+size(vy)+size(vz),'dp',-1)
deallocate (vx,vy,vz)
if (params%debug.gt.1) call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',-1)
deallocate (dhdx,dhdy,dhdz)
if (params%debug.gt.1) call heap (threadinfo,'crittemp','improve_osolve',size(crittemp),'dp',-1)
deallocate (crittemp)
!---------------------------------------------------
! we can now use crit to refine the osolve octree
!---------------------------------------------------
if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
call show_time (total,step,inc,1,'improve osolve for criterion$')
if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,es15.6)') shift//'max(crit) =',maxval(crit)
if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,es15.6)') shift//'max(crit,w_l_i_f)=',maxval(crit,ov%whole_leaf_in_fluid)
if(istep.eq.1 .and. iter.eq.1) then
critmax=maxval(crit)
else
critmax=maxval(crit,ov%whole_leaf_in_fluid)
end if
if (critmax.gt.tiny(critmax)) then
do i=1,ov%nleaves
!if(increase_current_level) then
! if (crit(i).gt. critmax*2.d0*params%refine_ratio ) then
! xxx=0.5d0*(ov%x(ov%icon(1,i))+ov%x(ov%icon(8,i)))
! yyy=0.5d0*(ov%y(ov%icon(1,i))+ov%y(ov%icon(8,i)))
! zzz=0.5d0*(ov%z(ov%icon(1,i))+ov%z(ov%icon(8,i)))
! if (ov%whole_leaf_in_fluid(i)) &
! call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
! end if
!else
if (crit(i).gt. critmax*params%refine_ratio ) then
xxx=0.5d0*(ov%x(ov%icon(1,i))+ov%x(ov%icon(8,i)))
yyy=0.5d0*(ov%y(ov%icon(1,i))+ov%y(ov%icon(8,i)))
zzz=0.5d0*(ov%z(ov%icon(1,i))+ov%z(ov%icon(8,i)))
if (ov%whole_leaf_in_fluid(i)) &
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
end if
!end if
enddo
end if
if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after crit refinement ',osolve%octree(2)
if (params%debug.gt.1) call heap (threadinfo,'crit','improve_osolve',size(crit),'dp',-1)
deallocate (crit)
!-----------------------------------------------
! osolve has yet to be refined in boxes, if any
!-----------------------------------------------
if (params%nboxes.gt.0) then
call show_time (total,step,inc,1,'Improve osolve for imposed boxes$')
do nn=1,params%nboxes
x0=boxes(nn)%x0
x1=boxes(nn)%x1
y0=boxes(nn)%y0
y1=boxes(nn)%y1
z0=boxes(nn)%z0
z1=boxes(nn)%z1
level_box=boxes(nn)%level
nx=int((x1-x0)/2.d0**(-level_box))! nb of box particles to insert on the x axis
ny=int((y1-y0)/2.d0**(-level_box))! nb of box particles to insert on the y axis
nz=int((z1-z0)/2.d0**(-level_box)) ! nb of box particles to insert on the z axis
do i=1,nx
do j=1,ny
do k=1,nz
xxx=x0+(real(i)-.5d0)*2.d0**(-level_box)
yyy=y0+(real(j)-.5d0)*2.d0**(-level_box)
zzz=z0+(real(k)-.5d0)*2.d0**(-level_box)
if (.not.((xxx>=x0).and.(xxx<=x1) .and. &
(yyy>=y0).and.(yyy<=y1) .and. &
(zzz>=z0).and.(zzz<=z1))) call stop_run ('pb with box cloud$')
! if (iproc.eq.0) print*,i,j,k
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level_box)
enddo
enddo
enddo
enddo
if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after box refinement ',osolve%octree(2)
endif
!-----------------------------------------------
! refine the octree on the cube faces
!----------------------------------------------
if (params%ref_on_faces) then
call show_time (total,step,inc,1,'Improve osolve on cube faces$')
do k=1,6
! if (k==5) then
! level=refine_level
! else
level=cube_faces(k)%level
! end if
l=cube_faces(k)%l
r=cube_faces(k)%r
t=cube_faces(k)%t
b=cube_faces(k)%b
nv=int((t-b)/2.d0**(-(level+1))) ! nb of box particles to insert on the vertical axis
nh=int((r-l)/2.d0**(-(level+1))) ! nb of box particles to insert on the horizontal axis
if (l > r .or. b > t) call stop_run ('pb in improve_octree_on_faces$')
select case (k)
case (1) ! x=0 face
xxx=0.00001
do i=1,nh
do j=1,nv
yyy=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
case (2) ! x=1 face
xxx=0.99999
do i=1,nh
do j=1,nv
yyy=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
case (3) ! y=0 face
yyy=0.00001
do i=1,nh
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
case (4) ! y=1 face
yyy=0.99999
do i=1,nh
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
case (5) ! z=0 face
zzz=0.00001
do i=1,nh
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
yyy=b+(real(j)-.5d0)*2.d0**(-(level+1))
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
case (6) ! z=1 face
zzz=0.99999
do i=1,nh
do j=1,nv
xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
yyy=b+(real(j)-.5d0)*2.d0**(-(level+1))
call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
enddo
enddo
end select
enddo
if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after faces refinement ',osolve%octree(2)
end if
!-----------------------------------------------
! finally smooth the octree
!-----------------------------------------------
call show_time (total,step,inc,1,'Smoothen octree after refinement$')
call octree_smoothen (osolve%octree,osolve%noctree)
if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after smoothing ',osolve%octree(2)
!-----------------------------------------------
! super smooth the octree if needed
!-----------------------------------------------
if (params%ismooth) then
call show_time (total,step,inc,1,'Super smoothing of osolve$')
call octree_super_smoothen (osolve%octree,osolve%noctree)
if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after super smoothing ',osolve%octree(2)
end if
if (params%debug.gt.1) call flush(threadinfo%Logunit)
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|