-
Douglas Guptill authoredDouglas Guptill authored
improve_osolve.f90 14.88 KiB
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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 threads
use definitions
use gauss
use invariants
implicit none
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
INCLUDE 'mpif.h'
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) ; 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) ; 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) ; call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',+1)
allocate (crit(ov%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'crit','improve_osolve',size(crit),'dp',+1)
allocate (crittemp(ov%nleaves),stat=threadinfo%err) ; 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
! 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 (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)
call heap (threadinfo,'x/y/z','improve_osolve',size(x)+size(y)+size(z),'dp',-1) ; deallocate (x,y,z)
call heap (threadinfo,'vx/vy/vz','improve_osolve',size(vx)+size(vy)+size(vz),'dp',-1) ; deallocate (vx,vy,vz)
call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',-1) ; deallocate (dhdx,dhdy,dhdz)
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_particles(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_particles(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)
write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after crit refinement ',osolve%octree(2)
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$')
call octree_create_from_particles(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)
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_particles(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_particles(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_particles(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_particles(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_particles(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_particles(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)
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)
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)
write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after super smoothing ',osolve%octree(2)
end if
call flush(threadinfo%Logunit)
return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|