!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! DEFINE_BC_SEGMENTED_S_LINE_PARABOLA August 2015 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine define_bc_segmented_s_line_parabola(params,osolve,vo,bcdef,nest)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine )))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine assigns the velocity boundary conditions for the segmented s-line
! geometry
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
!use mpi
implicit none
include 'mpif.h'
type (parameters) params
type (octreesolve) osolve
type (void) vo
type (bc_definition) bcdef
type (nest_info) nest
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer i,iproc,nproc,ierr
double precision :: eps,lsf0,pi,lorig,h,x1,x2,phi,yend,cper,cscl,xstart,ystart
double precision :: theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth
double precision :: ywidth,xdisp,ydisp,nb
double precision,dimension(:),allocatable :: x0,ldisp
integer ie,ij,j,jp,nelemx,nelemz
double precision :: base,startp,endp,alpha,kink1,kink2,flytt,a,d
double precision :: wmax,uend,vmag,dipangle
double precision e(1),f(1),y(1)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
eps=1.d-10
osolve%kfix=0
osolve%kfixt=0
pi=atan(1.d0)*4.d0
l=bcdef%bc_parameters(1)
base=bcdef%bc_parameters(2)
startp=bcdef%bc_parameters(3)
endp=bcdef%bc_parameters(4)
vin=bcdef%bc_parameters(5)
alpha=bcdef%bc_parameters(6)
kink1=bcdef%bc_parameters(7)
kink2=bcdef%bc_parameters(8)
nelemx=idint(bcdef%bc_parameters(9))
nelemz=idint(bcdef%bc_parameters(10))
nb=2**params%levelmax_oct
dxy=1.d0/2**(params%levelmax_oct+1.d0)
alpha=alpha*pi/180
! allocate(x0(osolve%nnode))
! do i=1,osolve%nnode
! if (osolve%y(i).le.ystart) then
! x0(i)=xstart
! else if (osolve%y(i).le.yend) then
! x0(i)=(osolve%y(i)-ystart)*tan(theta)+xstart
! else
! x0(i)=xend
! endif
! enddo
flytt=(kink2-kink1)* tan(alpha)
a = (base-l) /((endp-startp)*(endp-startp))
d = startp
wmax=(2*vin*l)/(endp-startp)
uend=wmax/(2*a*(endp-startp))
vmag=sqrt(uend**2.d0+wmax**2.d0)
do i=1,osolve%nnode
if (osolve%x(i).lt.eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=1.d0
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
endif
if (osolve%x(i).gt.1.d0-eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
endif
if (osolve%y(i).lt.eps) then
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
endif
if (osolve%y(i).gt.1.d0-eps) then
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
endif
if (osolve%z(i).lt.eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
osolve%kfixt(i)=1 ; osolve%temp(i)=1.d0
if (osolve%y(i) .lt. kink1) then
if (osolve%x(i) .le. startp) then
osolve%u(i)=1.d0
elseif (osolve%x(i) .gt. startp .AND. osolve%x(i) .lt. endp) then
dipangle=atan(2*a*(osolve%x(i)-d))
osolve%u(i) = vmag*cos(dipangle)
osolve%w(i) = vmag*sin(dipangle)
else
osolve%u(i)=0.d0
osolve%v(i)=0.d0
osolve%w(i)=0.d0
endif
endif
if ( osolve%y(i) .ge. kink1 .AND. osolve%y(i) .lt. kink2 ) then
if (osolve%x(i).lt. (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)) + startp)) then
osolve%u(i)=1.d0
elseif (osolve%x(i) .gt. (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)) + (startp)) .and. osolve%x(i) .lt. (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1))) + endp) then
osolve%u(i) = -(a*sqrt(8*(osolve%x(i)- (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1))) )*d)+0.2)*vin/1
osolve%w(i) = 10*2*a*sqrt( (osolve%x(i)- (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)))) * (osolve%x(i)- (osolve%y(i) * ((flytt)/(kink2-kink1)) - ((flytt*kink1)/(kink2-kink1)))) - d*d )*vin/30
else
osolve%u(i)=0.d0
endif
endif
if ( osolve%y(i) .ge. kink2 .AND. osolve%y(i) .lt. (1.0-eps) ) then
if (osolve%x(i) .le. (startp + flytt)) then
osolve%u(i)=1.d0
! elseif (osolve%x(i) .gt. (startp+flytt) .AND. osolve%x(i) .lt. (startp+flytt+0.01)) then
!! osolve%w(i)= (-2*a*(base-d)*vin+(osolve%x(i)-flytt))
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (10*vin*(osolve%x(i)-flytt))/(2.9d0)
! elseif (osolve%x(i) .gt. (startp+flytt+0.01) .AND. osolve%x(i) .lt. (startp+flytt+0.04)) then
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (6*vin*(osolve%x(i)-flytt))/(2.9d0)
! elseif (osolve%x(i) .gt. (startp+flytt+0.04) .AND. osolve%x(i) .lt. (startp+flytt+0.05)) then
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (4*vin*(osolve%x(i)-flytt))/(2.9d0)
! elseif (osolve%x(i) .gt. (startp+flytt+0.05) .AND. osolve%x(i) .lt. (endp+flytt-0.05)) then
! osolve%w(i)= (-2*a*(base-d)*vin)/(2.9d0)
! osolve%u(i)= (2*vin*(osolve%x(i)-flytt))/(2.9d0)
elseif (osolve%x(i) .gt. (startp+flytt) .AND. osolve%x(i) .lt. (endp+flytt)) then
! osolve%w(i)= (2*a*(osolve%x(i)-(d+flytt)))*vin
! osolve%u(i)= (2*vin*(osolve%x(i)-flytt))
osolve%u(i) = -(a*sqrt(8*(osolve%x(i)-flytt)*d)+0.2)*vin/1
osolve%w(i) = 10*2*a*sqrt( (osolve%x(i)-flytt)*(osolve%x(i)-flytt) - d*d )*vin/30
else
osolve%u(i)=0.d0
osolve%v(i)=0.d0
osolve%w(i)=0.d0
endif
endif
endif
if (osolve%z(i).gt.1.d0-eps) then
osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
osolve%kfixt(i)=1 ; osolve%temp(i)=0.d0
endif
if (.not.vo%influid(i)) then
osolve%kfixt(i)=1
osolve%temp(i)=0.d0
endif
if (abs(bcdef%utrans).gt.eps) then
if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)+bcdef%utrans
endif
if (abs(bcdef%vtrans).gt.eps) then
if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)+bcdef%vtrans
endif
if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)*vin
if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)*vin
if (osolve%kfix((i-1)*3+3)==1) osolve%w(i)=osolve%w(i)*vin
enddo
if (params%isobc) then
call define_isostasy_bc(params,osolve,bcdef)
endif
end
!------------------------------------------------------------------------------|