Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! COMPACTION MARCH 2013 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine compaction(params,density_str,osolve,ov,mat,dtcur,write_compp)
use definitions
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine compares density profiles, for each possible element centre |
! location at the base of the model, between the start of the time step (past |
! density) and the current mechanical solution (current density). Compaction |
! is called after erosion and sedimentation have been applied. This subroutine|
! determines the velocity field necessary to account for increased density in |
! compactible materials, compared to the start of time step density |
! configuration. Velocities are stores on ____, to be added to the nodal |
! velocities used to move surfaces and advect the cloud. Current density is |
! updated with each call to compaction.
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
implicit none
Dave Whipp
committed
include 'mpif.h'
type (parameters) params
type (string) :: density_str(2**params%levelmax_oct,2**params%levelmax_oct)
! In how much detail do we need to specify density_str? We will not be
! allocating/deallocating and of the density_str entires in compaction or
! calculate_density_profile.
! Note we are passing in the full structure here.
! We also need material number for each element, and compaction parameters
! for compactible material types.
! We also need LSFs, to check if we're above/below the free surface.
! Suggest for now we consider only the free surface. This could be
! adapted to look for other LSFs, to better define transitions between
! material types.
type(octreesolve) osolve
type(octreev) ov
type (material) mat(0:params%nmat)
double precision :: dtcur
integer :: write_compp
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
Dave Whipp
committed
!integer :: i,mat_maxi!,
integer :: iproc,nproc,ierr
!logical :: mat_max_tie
double precision,dimension(:),allocatable :: str_velo,countnode,countednode
double precision :: xstr,ystr,zstr,x0,y0,z0,dxyz,zmax_elem
double precision :: str_length_inc,eps
integer :: i,j,k,size_str(2),leaf,level,loc,nstr_per_elem,nstr_pts_per_elem
integer :: inode
logical :: next_elem
! 100414 JCA
logical,dimension(:), allocatable :: calc_comp_velo
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
Dave Whipp
committed
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
eps=1d-10
Dave Whipp
committed
! Loop over all instances of density_str (i,j 1:2**params%levelmax_oct), for each
! string:
! 1. Call calculate_density_profile(density_str(i,j), ??)
! - this will also need material number for each element, and compaction
! parameters for all materials (some of which may not be compactible)
! 2. Compare density_str(i,j)%density and %densityp
! - look at each point on density_str(i,j), moving from base of model
! to intersection with free surface (found when density==0+/-smallnumber)
! 3. If density>densityp, calculate velocity required to account for compaction
! 4. Interpolate velocities from density_str locations onto element nodes.
! - easier to first interpolate onto element centers, looking over all
! elements?
! - will need to account for the number of velocity contributions to each
! node
! Loop through all instances of density_str (all possible element center
! locations at the base of the model
! Zero out existing osolve%compaction_density
osolve%compaction_density=0.d0
! Could calculate this before calc_density_profile and pass in
size_str=size(density_str(1,1)%density)
! Allocate str_velo
allocate(str_velo(size_str(1)))
! 100414 JCA
allocate(calc_comp_velo(size_str(1)))
! Allocate countnode
! Is there a smarter way to do this????????
allocate (countnode(osolve%nnode),countednode(osolve%nnode))
countnode=0
countednode=1
do i=1,osolve%nleaves
do j=5,8
! write (*,*) 'osolve%nleaves: ',osolve%nleaves
! write (*,*) 'i: ',i
! write (*,*) 'j: ',j
! write (*,*) 'ov%icon(j,i): ',ov%icon(j,i)
k = osolve%icon(j,i)
! write (*,*) 'k: ',k
countnode(k) = countnode(k) + 1
countednode(k)=0
enddo
enddo
countnode=countnode+countednode
deallocate(countednode)
! String length increment
Dave Whipp
committed
str_length_inc = (1.d0/(2.d0**dble(params%levelmax_oct))*params%vex)/dble(params%dstrnz)
do j=1,2**params%levelmax_oct
do i=1,2**params%levelmax_oct
! 100414 JCA
do k=1,size_str(1)
calc_comp_velo(k)=.false.
enddo
!Update current density, density_str(i,j)%density
call calculate_density_profile(params,density_str(i,j),osolve,mat,0,i,j, &
calc_comp_velo,write_compp)
Dave Whipp
committed
! if (j.eq.2**(params%levelmax_oct-1)) then
! if (i.eq.3 .or. i.eq.2**(params%levelmax_oct-1) .or. i.eq.2**(params%levelmax_oct)-3) then
! if (iproc.eq.0) then
! write (*,*) 'density string output (',i,',',j,')'
! write (*,*) 'density, densityp'
! do k=1,size_str(1)
! write (*,*) density_str(i,j)%density(k),', ',density_str(i,j)%densityp(k)
! enddo
! endif
! endif
! endif
Dave Whipp
committed
! Could calculate this before calc_density_profile and pass in
Dave Whipp
committed
xstr = (1.d0/(2.d0**dble(params%levelmax_oct)))*(0.5d0+dble(i-1))
ystr = (1.d0/(2.d0**dble(params%levelmax_oct)))*(0.5d0+dble(j-1))
next_elem = .true.
str_velo = 0.d0
do k=1,size_str(1)
Dave Whipp
committed
! Note we're now accounting for squishing
zstr = (1.d0/(2.d0**dble(params%levelmax_oct))*params%vex)*(0.5d0+dble(k-1))/dble(params%dstrnz)
if (next_elem) then
Dave Whipp
committed
! Here we assume octree_find_leaf expects and returns unsquished z values
call octree_find_leaf (osolve%octree,osolve%noctree,xstr,ystr,zstr/params%vex, &
leaf,level,loc,x0,y0,z0,dxyz)
nstr_per_elem = int(2.d0**(2.d0*dble(params%levelmax_oct - level)))
nstr_pts_per_elem = int(2.d0**(dble(params%levelmax_oct - level)) * &
dble(nstr_per_elem) * dble(params%dstrnz))
Dave Whipp
committed
! This takes unsquished values from octree_find_leaf and squishes
zmax_elem = (z0 + dxyz / 2.d0)*params%vex
next_elem = .false.
endif
osolve%compaction_density(leaf) = osolve%compaction_density(leaf) + density_str(i,j)%density(k)/dble(nstr_pts_per_elem)
! 100414 JCA
if (calc_comp_velo(k)) then
if (abs(density_str(i,j)%density(k)) < eps .or. abs(density_str(i,j)%densityp(k)) < eps) then
str_velo(k) = 0.d0
else
str_velo(k) = str_length_inc * ((density_str(i,j)%densityp(k)/density_str(i,j)%density(k))-1.d0)/dtcur
if (iproc==0) then
write (*,*) 'str_velo(k): ',str_velo(k)
write (*,*) 'str_length_inc: ',str_length_inc
write (*,*) 'density_str(i,j)%densityp(k): ',density_str(i,j)%densityp(k)
write (*,*) 'density_str(i,j)%density(k): ',density_str(i,j)%density(k)
write (*,*) 'dtcur: ',dtcur
endif
endif
else
str_velo(k) = 0.d0
endif
! JANICE2: Commented out debug output below
! write (*,*) 'str_velo(k): ',str_velo(k)
! write (*,*) 'str_length_inc: ',str_length_inc
! write (*,*) 'density_str(i,j)%densityp(k): ',density_str(i,j)%densityp(k)
! write (*,*) 'density_str(i,j)%density(k): ',density_str(i,j)%density(k)
! write (*,*) 'dtcur: ',dtcur
if (k > 1) str_velo(k) = str_velo(k) + str_velo(k-1)
if (zmax_elem - zstr .le. str_length_inc) then
do inode=5,8
! This is OK for a uniform octree at octreelvl_max, but the velocites
! should be scaled by distance to the nodes for elements with multiple
! strings!
if (iproc == 0) then
write (*,*) 'ov%wnodecompact(ov%icon(inode,leaf)) 1: ',ov%wnodecompact(ov%icon(inode,leaf))
endif
ov%wnodecompact(ov%icon(inode,leaf)) = ov%wnodecompact(ov%icon(inode,leaf)) + str_velo(k) / dble(nstr_per_elem)
if (iproc == 0) then
write (*,*) 'ov%wnodecompact(ov%icon(inode,leaf)) 2: ',ov%wnodecompact(ov%icon(inode,leaf))
write (*,*) 'str_velo(k): ',str_velo(k)
write (*,*) 'dble(nstr_per_elem): ',dble(nstr_per_elem)
endif
enddo
next_elem = .true.
endif
enddo
enddo
enddo
do i = 1,osolve%nnode
if (iproc == 0) then
write (*,*) 'ov%wnodecompact(i) 1: ',ov%wnodecompact(i)
endif
ov%wnodecompact(i) = ov%wnodecompact(i) / dble(countnode(i))
if (iproc == 0) then
write (*,*) 'ov%wnodecompact(i) 2: ',ov%wnodecompact(i)
write (*,*) 'dble(countnode(i)): ',dble(countnode(i))
endif
! JANICE2: Commented out debug output below
! write (*,*) 'ov%wnodecompact(i): ', ov%wnodecompact(i)
! write (*,*) 'dble(countnode(i)): ',dble(countnode(i))
! write (*,*) 'ov%wnode(i): ',ov%wnode(i)
ov%wnode(i)=ov%wnode(i)+ov%wnodecompact(i)
! JANICE2: Commented out debug output below
! write (*,*) 'ov%wnode(i) + wcompact: ',ov%wnode(i)
enddo
!Note: if storing velocity at a differently sized array, could compare
! density, densityp outside the above loop over all instances of density_str
! 100414 JCA
! deallocate calc_comp_velo? deallocate str_velo?
deallocate(countnode)
end subroutine compaction
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------