!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
program DOUAR
use threads
use definitions
implicit none
INCLUDE 'mpif.h'
integer i,j,is,iter,istep,niter,nedge,nedgen
integer ntriangle,ndof,ndoft,material0,nsurfacen
integer numarg,iargc
integer ierr,nproc,iproc,k,nrem
integer current_level
integer ie_ov,ie_loc,ie_level,iter_nl
integer err,iunit,cltemp
integer nremove,ninject,nnmax,nadd,naddp,ref_count
integer nleaves_new_mem1
integer nleaves_new_mem2
integer nleaves_old_mem1
integer nleaves_old_mem2
integer, external :: ioctree_number_of_elements
double precision x,y,z,z0,u,v,w
double precision exx,eyy,ezz,exy,eyz,ezx,e2dmax
double precision duvw,uvw,dtcourant,current_time
double precision umax,xxx,yyy,zzz,x0_leaf,y0_leaf,z0_leaf
double precision total,step,inc
real time1,time2,time3
character :: ic*7
character*3 :: ciproc
character*4 :: cs,cs2
character*40 :: string
character*72 :: shift
logical converge,increase_current_level,velocity_converged,usecourant
type (sheet),dimension(:),allocatable::surface,surface0
type (octreev) ov
type (octreelsf) olsf
type (octreesolve) osolve
type (material),dimension(:),allocatable::mat
type (void) vo
type (cloud) cl
type (topology),dimension(:),allocatable::tpl
type (box),dimension(:),allocatable::boxes
type (cross_section),dimension(:),allocatable::sections
type (face),dimension(6) :: cube_faces
type (edge),dimension(:),allocatable::ed,edswap
type (parameters) params
type (thread) threadinfo
type (ziso) zi
integer n_iproc_st,n_iproc_end,n_iproc
integer ldb,nrhs,n,nz_loc
integer, dimension(:), allocatable :: ia,ja
logical, dimension(:), allocatable :: iproc_col
double precision, dimension(:), allocatable :: avals
double precision, dimension(:,:), allocatable :: b
double precision, dimension(:), allocatable :: weightel
double precision, dimension(:), allocatable :: xyz_t
shift=' '
ndof=3
ndoft=1
params%mpe=8
nrhs = 1
nleaves_old_mem1=0
nleaves_old_mem2=0
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_init(ierr)
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
call write_streepjes(6,2)
!=====[allocate threadinfo and open mpi log and memory heap files]=========
params%nmpi=nproc
call int_to_char(ciproc,3,iproc)
threadinfo%Logunit=1000+iproc
open (unit=threadinfo%Logunit,file='./DEBUG/mpilogs/Log_'//ciproc//'.dat',status='replace')
write(threadinfo%Logunit,*) 'This is douar, last modifed 2009-12-21'
if (iproc.eq.0) write(*,*) 'This is douar, last modifed 2009-12-21, using ', nproc, ' processors.'
write(threadinfo%Logunit,'(a,i3)') 'Log file of mpi process',iproc
call heap_init(threadinfo,2000+iproc,'./DEBUG/mpilogs/mem_heap_'//ciproc//'.dat')
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! is there any input file to douar ?
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
numarg=iargc()
! 2009-09-03: Douglas Guptill
! mpich doesn't allow us to read the command line
! so we simulate an empty one
numarg = 0
if (numarg==0) then
params%infile='input.txt'
if (iproc.eq.0) write(*,*) 'program called with no argument'
else
call getarg (1,params%infile)
if (iproc.eq.0) then
write(*,*) 'program called with input file ',params%infile
end if
end if
call write_streepjes(6,2)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! read controlling parameters
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call read_controlling_parameters (params,threadinfo)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! open and prepare output files
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call io_DoRuRe2 (params,'init')
call show_time (total,step,inc,0,'Start of Computations$')
call show_time (total,step,inc,1,'Reading Input$')
allocate (surface ( params%ns ),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface in main$')
allocate (surface0 ( params%ns ),stat=err) ; if (err.ne.0) call stop_run ('Error alloc surface0 in main$')
allocate (mat (0:params%nmat),stat=err) ; if (err.ne.0) call stop_run ('Error alloc mat in main$')
allocate (params%materialn(0:params%ns ),stat=err); if (err.ne.0) call stop_run ('Error alloc materialn in main$')
if (params%nboxes.gt.0) then
allocate (boxes(params%nboxes),stat=err) ; if (err.ne.0) call stop_run ('Error alloc boxes arrays$')
else
allocate (boxes(1),stat=err) ! necessary to avoid nil size argument
endif
if (params%nsections.gt.0) then
allocate (sections(params%nsections),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cross_section arrays$')
else
allocate (sections(1),stat=err) ! necessary to avoid nil size argument
end if
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! read input file for all parameters of the run
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call read_input_file (params,threadinfo,material0,mat,surface,boxes,sections,cube_faces)
if (params%irestart==0) then
current_level=params%initial_refine_level
else
current_level=params%levelmax_oct
end if
!if (params%isobc) allocate(zisodisp(2**params%levelmax_oct+1,2**params%levelmax_oct+1),stat=err) ; if (err.ne.0) call stop_run ('Error alloc zisodisp in main$')
call show_time (total,step,inc,1,'Problem Setup$')
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! either reads or creates the surface(s)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call show_time (total,step,inc,1,'define surfaces$')
call define_surface(params,surface,threadinfo,total,step,inc,current_time)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! either reads or creates the cloud
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call show_time (total,step,inc,1,'define cloud$')
call define_cloud(cl,params,zi)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! extract material information to be passed to solver
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
params%materialn(0)=material0
do i=1,params%ns
params%materialn(i)=surface(i)%material
enddo
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! compute plastic parameters
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call compute_plastic_params (params,mat)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! either reads or creates the velocity octree
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call show_time (total,step,inc,1,'define velocity octree$')
call define_ov (ov,params,threadinfo)
istep=1+params%irestart
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! is the CFL condition used for the timestep ?
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
usecourant = (params%dt .lt. 0.d0)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! beginning of time stepping
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
do while (istep.le.params%nstep)
call write_streepjes(6,2)
call write_streepjes(8,2)
call show_time (total,step,inc,-istep,'Start of Step $')
if (iproc.eq.0) write(8,*) 'Doing step ',istep
call write_streepjes(6,2)
call write_streepjes(8,2)
call heap_hop1(threadinfo,istep)
!---------------------------------------------------------------------------|
!---------------------------------------------------------------------------|
! refining surfaces
!---------------------------------------------------------------------------|
!---------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Loop over surfaces$')
do is=1,params%ns
call int_to_char(ic,2,is)
call show_time (total,step,inc,1,'surface '//ic(1:2)//'$')
if (iproc.eq.0 ) then
write (8,*) 'Surface ',is,' has ',surface(is)%nsurface,'particles before refinement'
if (params%debug>=1) write(*,'(a,i2,a,i8,a)') shift//'S.',is,': ',surface(is)%nsurface,'points'
end if
if (surface(is)%activation_time.ge.current_time+tiny(current_time)) then
!------------------------------------------------------------------------!
! if surface is slave to top surface, no refinement
!------------------------------------------------------------------------!
deallocate(surface(is)%x,surface(is)%y,surface(is)%z)
deallocate(surface(is)%xn,surface(is)%yn,surface(is)%zn)
deallocate(surface(is)%r,surface(is)%s)
deallocate(surface(is)%icon)
surface(is)%nsurface=surface(1)%nsurface
surface(is)%nt=surface(1)%nt
allocate (surface(is)%x(surface(is)%nsurface),surface(is)%y(surface(is)%nsurface),surface(is)%z(surface(is)%nsurface))
allocate (surface(is)%xn(surface(is)%nsurface),surface(is)%yn(surface(is)%nsurface),surface(is)%zn(surface(is)%nsurface))
allocate (surface(is)%r(surface(is)%nsurface),surface(is)%s(surface(is)%nsurface))
allocate (surface(is)%icon(3,surface(is)%nt))
surface(is)%x=surface(1)%x
surface(is)%y=surface(1)%y
surface(is)%z=surface(1)%z
surface(is)%xn=surface(1)%xn
surface(is)%yn=surface(1)%yn
surface(is)%zn=surface(1)%zn
surface(is)%r=surface(1)%r
surface(is)%s=surface(1)%s
surface(is)%icon=surface(1)%icon(:,1:surface(is)%nt)
else
!------------------------------------------------------------------------!
! if surface is no slave to top surface, do refine
!------------------------------------------------------------------------!
ref_count=1
nadd=1
do while (nadd>0)
call show_time (total,step,inc,1,'refine surface$')
call refine_surface(params,surface(is),threadinfo,nadd,is,istep,ref_count)
call show_time (total,step,inc,1,'check delaunay$')
call check_delaunay (params,surface(is),is,istep,'delaunay_after_refine',ref_count)
call DoRuRe_surf_stats (params%doDoRuRe,istep,ref_count,is,surface(is)%nt,surface(is)%nsurface,nedge,nadd)
ref_count=ref_count+1
end do
endif
if (iproc.eq.0 ) then
write (8,*) 'Surface ',is,' has ',surface(is)%nsurface,'particles after refinement'
if (params%debug>=1) write(*,'(a,i2,a,i8,a)') shift//'S.',is,': ',surface(is)%nsurface,'points'
end if
end do
!---------------------------------------------------------------------------|
!---------------------------------------------------------------------------|
! Stores the geometry of the surfaces for the mid-point iterative scheme
!---------------------------------------------------------------------------|
!---------------------------------------------------------------------------|
do is=1,params%ns
surface0(is)%nsurface=surface(is)%nsurface
allocate (surface0(is)%r(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',+1)
allocate (surface0(is)%s(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',+1)
allocate (surface0(is)%x(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',+1)
allocate (surface0(is)%y(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',+1)
allocate (surface0(is)%z(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',+1)
allocate (surface0(is)%xn(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',+1)
allocate (surface0(is)%yn(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',+1)
allocate (surface0(is)%zn(surface0(is)%nsurface),stat=threadinfo%err)
call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',+1)
surface0(is)%r=surface(is)%r
surface0(is)%s=surface(is)%s
surface0(is)%x=surface(is)%x
surface0(is)%y=surface(is)%y
surface0(is)%z=surface(is)%z
surface0(is)%xn=surface(is)%xn
surface0(is)%yn=surface(is)%yn
surface0(is)%zn=surface(is)%zn
enddo
!---------------------------------------------------------------------------|
!---------------------------------------------------------------------------|
! beginning of grid iterations
!---------------------------------------------------------------------------|
!---------------------------------------------------------------------------|
iter =1
do while (iter.le.abs(params%griditer).or.iter.eq.1)
call write_streepjes(6,1)
call write_streepjes(8,1)
call show_time (total,step,inc,-iter,' Start of iteration $')
if(iproc.eq.0) write(8,*) 'Doing iteration ',iter
call write_streepjes(6,1)
call write_streepjes(8,1)
call heap_hop2(threadinfo,iter)
call show_time (total,step,inc,1,'Create octree solve$')
converge=.false.
! if (params%griditer .lt. 0 .and. current_level==params%levelmax_oct) converge=.true.
if (iproc.eq.0 .and. params%debug>=1) then
write(*,'(a,L1)') shift//'(1) params%griditer < 0 ',(params%griditer.lt.0)
write(*,'(a,L1)') shift//'(2) current_level==params%levelmax_oct',(current_level==params%levelmax_oct)
write(*,'(a,L1)') shift//'(3) increase_current_level ',increase_current_level
write(*,'(a,L1)') shift//'converge = (1) & (2) & (3) -> ',(params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level)
end if
if (params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level) converge=.true.
if (iproc.eq.0 .and. converge .and. params%debug>=1) then
write(*,'(a)') shift//'+++++++++++++++++'
write(*,'(a)') shift//'grid conv reached'
write(*,'(a)') shift//'+++++++++++++++++'
end if
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! creation of osolve
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! we create the large "solve" octree, ie the octree on which the finite
! elements will be built.
osolve%noctree=params%noctreemax
allocate (osolve%octree(osolve%noctree),stat=threadinfo%err)
call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',+1)
call octree_init (osolve%octree,osolve%noctree)
call octree_create_uniform (osolve%octree,osolve%noctree,params%leveluniform_oct)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! refinement/improvement of osolve
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call improve_osolve (osolve,ov,params,threadinfo,istep,iter,total,step,inc, &
current_level,boxes,cube_faces,increase_current_level )
nleaves_new_mem1=osolve%octree(2)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! find connectivity for osolve
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Build osolve icon$')
osolve%nleaves = ioctree_number_of_elements (osolve%octree,osolve%noctree)
osolve%nnode = osolve%nleaves*3
osolve%nlsf = params%ns
if (iproc.eq.0) write (8,*) osolve%nleaves,' leaves in solve octree '
allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',+1)
allocate (osolve%x(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%x', 'main',size(osolve%x),'dp',+1)
allocate (osolve%y(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%y', 'main',size(osolve%y),'dp',+1)
allocate (osolve%z(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%z', 'main',size(osolve%z),'dp',+1)
call octree_find_node_connectivity (osolve%octree,osolve%noctree, &
osolve%icon,osolve%nleaves,osolve%x,osolve%y,osolve%z,osolve%nnode)
! osolve%nnode has been changed by octree_find_node_connectivity, so re-size x,y,z
call octreesolve_shrink_xyz(osolve, threadinfo)
! now that osolve%nnode is known we can allocate osolve%lsf
allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err)
call heap (threadinfo,'osolve%lsf', 'main',size(osolve%lsf),'dp',+1)
osolve%lsf=0.d0
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! embed the surfaces in osolve
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Embed surface in osolve$')
do is=1,params%ns
write (ic(1:2),'(i2)') is
call show_time (total,step,inc,1,'embedding surface '//ic(1:2)//'$')
call embed_surface_in_octree (osolve,params,surface(is),is,istep,iter,threadinfo)
enddo
nleaves_new_mem2=osolve%octree(2)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! assess grid convergence
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Assessing octree stability$')
increase_current_level= (nleaves_old_mem2 .ge. nleaves_new_mem2 .and. &
(dble(nleaves_old_mem2-nleaves_new_mem2)) / (dble(nleaves_old_mem2+nleaves_new_mem2))*2.d0 &
.lt. params%octree_refine_ratio )
! increase_current_level= ((dble(nleaves_old_mem2-nleaves_new_mem2)) / (dble(nleaves_old_mem2+nleaves_new_mem2))*2.d0 &
! .lt. params%octree_refine_ratio )
if(iproc.eq.0 .and. params%debug>=1) then
write(*,'(a,i5,a)') shift//'current refine level ',current_level
write(*,'(a)') shift//'After criterion based refinement:'
write(*,'(a,i8,a)') shift//'nleaves_old= ',nleaves_old_mem1,' leaves'
write(*,'(a,i8,a)') shift//'nleaves_new= ',nleaves_new_mem1,' leaves'
write(*,'(a)') shift//'After surfaces embedding:'
write(*,'(a,i8,a)') shift//'nleaves_old= ',nleaves_old_mem2,' leaves'
write(*,'(a,i8,a)') shift//'nleaves_new= ',nleaves_new_mem2,' leaves'
write(*,'(a,l1)') shift//'C2: authorise increase of refine level = ',increase_current_level
end if
nleaves_old_mem1=nleaves_new_mem1
nleaves_old_mem2=nleaves_new_mem2
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! allocate fields of osolve
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
allocate (osolve%u(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',+1)
allocate (osolve%v(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',+1)
allocate (osolve%w(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',+1)
allocate (osolve%wiso(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%wiso','main',size(osolve%wiso),'dp',+1)
allocate (osolve%temp(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',+1)
allocate (osolve%pressure(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',+1)
allocate (osolve%spressure(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',+1)
allocate (osolve%strain(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',+1)
osolve%u = 0.d0
osolve%v = 0.d0
osolve%w = 0.d0
osolve%wiso = 0.d0
osolve%temp = 0.d0
osolve%pressure = 0.d0
osolve%spressure = 0.d0
osolve%strain = 0.d0
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! improve and update cloud fields
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Change 3D cloud$')
call update_cloud_structure (cl,osolve,params,ninject,nremove,istep)
! if (istep==1) call strain_history (params,osolve,cl)
! call strain_history (params,osolve,cl)
call DoRuRe_cloud_stats (params%doDoRuRe,istep,iter,cl%np,nremove,ninject,cl%strain,cl%temp,cl%press)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! interpolate velocity from ov to solve, prepare ov receive the new solution
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Interpolate velo onto osolve$')
call interpolate_ov_on_osolve (osolve,ov,iter,params,threadinfo)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! find bad faces of solve octree
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Find bad faces$')
osolve%nface=osolve%nleaves
allocate (osolve%iface(9,osolve%nface),stat=threadinfo%err)
call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',+1)
call octree_find_bad_faces (osolve%octree,osolve%noctree, &
osolve%iface,osolve%nface, &
osolve%icon,osolve%nleaves)
if (osolve%nface.gt.osolve%nleaves) call stop_run ('nface larger than nleaves$')
if (iproc.eq.0) write (8,*) osolve%nface,' bad faces in solve octree'
if (iproc.eq.0) write (8,'(a,i9)') shift//'nb of bad faces',osolve%nface
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! find void (which nodes are in the void) and which elements are entirely in
! the void and should not be included in the fe calculations
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Find void$')
allocate (vo%node(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'vo%node','main',size(vo%node),'int',+1)
allocate (vo%leaf(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',+1)
allocate (vo%face(osolve%nface),stat=threadinfo%err)
call heap (threadinfo,'vo%face','main',size(vo%face),'int',+1)
allocate (vo%rtf(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',+1)
allocate (vo%ftr(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',+1)
allocate (vo%influid(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',+1)
call find_void_nodes (params,osolve,vo,istep,ov,threadinfo)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! definition of the bc
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'Define boundary conditions$')
allocate (osolve%kfix(osolve%nnode*3),stat=threadinfo%err)
call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',+1)
allocate (osolve%kfixt(osolve%nnode),stat=threadinfo%err)
call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',+1)
call define_bc (params,osolve,vo,zi)
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! build matrix arrays, allocate memory for wsmp
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'wsmp setup1$')
n = vo%nnode * ndof
!---[topology]-----
allocate (tpl(n))
tpl(:)%nheightmax=17*ndof
do i=1,n
allocate (tpl(i)%icol(tpl(i)%nheightmax))
enddo
call heap (threadinfo,'tpl(:)%icol(:)','main',n*tpl(1)%nheightmax,'int',+1)
!---[topology]-----
call wsmp_setup1 (n,n_iproc_st,n_iproc_end,n_iproc,nz_loc,ldb,ndof,vo,osolve,tpl,params,threadinfo,istep,iter)
allocate(iproc_col(n),stat=threadinfo%err) ; call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',+1)
allocate(ja(nz_loc),stat=threadinfo%err) ; call heap (threadinfo,'ja','main',size(ja),'int',+1)
allocate(ia(n_iproc+1),stat=threadinfo%err) ; call heap (threadinfo,'ia','main',size(ia),'int',+1)
allocate(avals(nz_loc),stat=threadinfo%err) ; call heap (threadinfo,'avals','main',size(avals),'dp',+1)
allocate(b(ldb,nrhs),stat=threadinfo%err) ; call heap (threadinfo,'b','main',size(b),'dp',+1)
call show_time (total,step,inc,1,'wsmp setup2$')
call wsmp_setup2 (n,n_iproc,n_iproc_st,n_iproc_end,nz_loc,iproc_col, &
ia,ja,ndof,vo,osolve,tpl,params,threadinfo,istep,iter)
!---[topology]-----
call heap (threadinfo,'tpl(:)%icol(:)','main',n*tpl(1)%nheightmax,'int',-1)
do i=1,n
deallocate (tpl(i)%icol)
enddo
deallocate (tpl)
!---[topology]-----
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! allocate measurement arrays
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
allocate (osolve%eviscosity(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',+1)
allocate (osolve%q(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',+1)
allocate (osolve%crit(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',+1)
allocate (osolve%e2d(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',+1)
allocate (osolve%e3d(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',+1)
allocate (osolve%lode(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',+1)
allocate (osolve%is_plastic(osolve%nleaves),stat=threadinfo%err)
call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',+1)
allocate (weightel(osolve%nleaves))
call heap (threadinfo,'weightel','main',size(weightel),'bool',+1)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! beginning of non linear iterations
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
do iter_nl=1,abs(params%nonlinear_iterations)
if (iproc.eq.0) then
write(*,*) '-----------------------------------'
write(8,*) '-----------------------------------'
call show_time(total,step,inc,-iter_nl,' start of non-linear iteration $')
write(8,*) 'Doing inner iteration ',iter_nl
write(*,*) '-----------------------------------'
write(8,*) '-----------------------------------'
end if
call heap_hop3(threadinfo,iter_nl)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! transfering previous solution from ov to osolve
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (iter_nl.gt.1) then
do i=1,ov%nnode
osolve%u(i)=ov%unode(i)
osolve%v(i)=ov%vnode(i)
osolve%w(i)=ov%wnode(i)
osolve%wiso(i)=ov%wnodeiso(i)
enddo
end if
! call define_bc (params,osolve,vo)
!if (iproc.eq.0) then
!write(*,*) minval(osolve%u),maxval(osolve%u)
!write(*,*) minval(osolve%v),maxval(osolve%v)
!write(*,*) minval(osolve%w),maxval(osolve%w)
!end if
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! building the FEM matrix and rhs
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'build system$')
call build_system_wsmp (n,n_iproc,n_iproc_st,nz_loc,iproc_col,ia,ja,ldb,nrhs,avals,b, &
params,osolve,ndof,mat,vo,istep,iter,iter_nl,threadinfo,weightel)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! solve system with wsmp solver
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'wsmp solve$')
call solve_with_pwssmp (n,nz_loc,n_iproc,n_iproc_st,n_iproc_end,ia,ja, &
ldb,nrhs,avals,b,params,osolve,ov,vo,threadinfo,&
ndof,istep,iter,iter_nl)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! calculate pressures
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Calculate pressures$')
call compute_pressure (params,osolve,ov,vo,mat)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! smoothing pressures
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Smoothing pressures$')
call smooth_pressures (osolve,ov,params)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! leaf measurements
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'do leaf measurements$')
call do_leaf_measurements (osolve,ov,params,istep,iter,iter_nl)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! write all informations in a text file
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (params%debug>=2) then
call show_time (total,step,inc,1,'write global output$')
call write_global_output(params,istep,iter,current_time,osolve,ov,vo,surface,cl,zi,'debug')
call mpi_barrier (mpi_comm_world,ierr)
! call slices (params,osolve,ov,vo,sections,istep,iter,iter_nl)
end if
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! convergence criterion computation
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'compute convergence criterion$')
call compute_convergence_criterion(osolve,ov,vo,params,istep,iter,iter_nl, &
current_level,velocity_converged)
cltemp=current_level !!++!!
if (velocity_converged .and. increase_current_level) current_level=min(current_level+1,params%levelmax_oct)
if (iproc.eq.0 .and. params%debug>=1) then
write(*,'(a,L1)') shift//'increase_current_level ->',increase_current_level
write(*,'(a,i2)') shift//'current_level ->',current_level
end if
if (cltemp+1 == params%levelmax_oct) increase_current_level=.false. !!++!!
if (params%nonlinear_iterations.lt.0 .and. velocity_converged) exit
end do
call heap_hop3_f (threadinfo)
call DoRuRe_nonlin_stats(params%doDoRuRe,istep,iter,iter_nl)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! end of non linear iterations
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'slicing the cube$')
! call slices(params,osolve,ov,vo,sections,istep,iter,iter_nl)
!call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',-1) ; deallocate (osolve%q)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! estimating the divergence (to check if incompressibility has been respected)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Compute divergence$')
call compute_divergence (params,osolve,ov,vo,istep,iter)
!-------------------------------------------------------------------------!
!-------------------------------------------------------------------------!
! deallocate memory used by the solver and terminates the solver's job
!-------------------------------------------------------------------------!
!-------------------------------------------------------------------------!
call show_time (total,step,inc,1,'wsmp_cleanup$')
call heap (threadinfo,'ia','main',size(ia),'int',-1) ; deallocate(ia)
call heap (threadinfo,'ja','main',size(ja),'int',-1) ; deallocate(ja)
call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',-1) ; deallocate(iproc_col)
call heap (threadinfo,'avals','main',size(avals),'dp',-1) ; deallocate(avals)
call heap (threadinfo,'b','main',size(b),'dp',-1) ; deallocate(b)
call heap (threadinfo,'weightel','main',size(weightel),'dp',-1) ; deallocate(weightel)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! check for convergence
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (iter.eq.abs(params%griditer)) converge=.true.
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! deallocate osolve
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! note that osolve will be needed in the temperature calculations
! so if this is the last iteration (converge is true), we will not deallocate osolve
if (.not.converge) then
call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',-1) ; deallocate (osolve%octree)
call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
call heap (threadinfo,'osolve%x','main',size(osolve%x),'dp',-1) ; deallocate (osolve%x)
call heap (threadinfo,'osolve%y','main',size(osolve%y),'dp',-1) ; deallocate (osolve%y)
call heap (threadinfo,'osolve%z','main',size(osolve%z),'dp',-1) ; deallocate (osolve%z)
call heap (threadinfo,'osolve%lsf','main',size(osolve%lsf),'dp',-1) ; deallocate (osolve%lsf)
call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',-1) ; deallocate (osolve%kfix)
call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',-1) ; deallocate (osolve%iface)
call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',-1) ; deallocate (osolve%strain)
call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',-1) ; deallocate (osolve%pressure)
call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',-1) ; deallocate (osolve%spressure)
call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',-1) ; deallocate (osolve%kfixt)
call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1) ; deallocate (osolve%u)
call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1) ; deallocate (osolve%v)
call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1) ; deallocate (osolve%w)
call heap (threadinfo,'osolve%wiso','main',size(osolve%wiso),'dp',-1) ; deallocate (osolve%wiso)
call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1) ; deallocate (osolve%temp)
call heap (threadinfo,'osolve%crip','main',size(osolve%crit),'dp',-1) ; deallocate (osolve%crit)
call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1) ; deallocate (osolve%e2d)
call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',-1) ; deallocate (osolve%e3d)
call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',-1) ; deallocate (osolve%lode)
call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',-1) ; deallocate (osolve%is_plastic)
!call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1) ; deallocate (ov%wpreiso)
end if
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! finds courant condition time step in case this is the first iteration
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (iter.eq.1) then
umax=0.d0
do i=1,ov%nnode
if (vo%influid(i)) umax=max(umax,sqrt(ov%unode(i)**2+ov%vnode(i)**2+ov%wnode(i)**2))
enddo
dtcourant=.5d0**params%levelmax_oct/umax*params%courant
if (usecourant) then
params%dt=dtcourant
if(iproc.eq.0 .and. params%debug .ge. 1) write(*,'(a,es15.6)') shift//'dt (CFL) =',params%dt
end if
endif
if (converge) iter=abs(params%griditer)
call DoRuRe_dt_stats (params%doDoRuRe,istep,params%dt)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! reset surface geometry to original geometry
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time(total,step,inc,1,'Reset surface geometry$')
do is=1,params%ns
surface(is)%r=surface0(is)%r
surface(is)%s=surface0(is)%s
surface(is)%x=surface0(is)%x
surface(is)%y=surface0(is)%y
surface(is)%z=surface0(is)%z
surface(is)%xn=surface0(is)%xn
surface(is)%yn=surface0(is)%yn
surface(is)%zn=surface0(is)%zn
enddo
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! update surface geometry by midpoint rule if not converge
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (.not.converge) then
olsf%nnode=ov%nnode
olsf%nleaves=ov%nleaves
olsf%noctree=ov%noctree
allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)
call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
olsf%y(1:olsf%nnode)=ov%y(1:ov%nnode)
olsf%z(1:olsf%nnode)=ov%z(1:ov%nnode)
olsf%icon(:,1:olsf%nleaves)=ov%icon(:,1:ov%nleaves)
do is=1,params%ns
if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
call move_surface (surface(is),surface0(is),1,0,ov,params%dt/2.d0,params,istep,is)
else
surface(is)%x=surface(1)%x
surface(is)%y=surface(1)%y
surface(is)%z=surface(1)%z
surface(is)%xn=surface(1)%xn
surface(is)%yn=surface(1)%yn
surface(is)%zn=surface(1)%zn
surface(is)%r=surface(1)%r
surface(is)%s=surface(1)%s
endif
! erosion added by Jean on Dec 12 2007
if (material0.eq.0.and.params%erosion) then
call int_to_char(ic,2,is)
call show_time (total,step,inc,1,'Erode surface '//ic(1:2)//'$')
if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
call erosion (surface(is),olsf,is,params)
else
if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
surface(is)%z=surface(1)%z
endif
end if
enddo
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! compute isostasy and move surfaces again
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (params%isostasy) then
call show_time (total,step,inc,1,'Compute isostasy and adjust vertical velocity$')
!allocate(ov%wpreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wpreiso in main$')
call isostasy (params,weightel,ov,surface,mat,1,zi)
do is=1,params%ns
if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
call move_surface (surface(is),surface0(is),1,1,ov,params%dt/2.d0,params,istep,is)
else
surface(is)%x=surface(1)%x
surface(is)%y=surface(1)%y
surface(is)%z=surface(1)%z
surface(is)%xn=surface(1)%xn
surface(is)%yn=surface(1)%yn
surface(is)%zn=surface(1)%zn
surface(is)%r=surface(1)%r
surface(is)%s=surface(1)%s
endif
enddo
endif
call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',-1) ; deallocate (olsf%x)
call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',-1) ; deallocate (olsf%y)
call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',-1) ; deallocate (olsf%z)
call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',-1) ; deallocate (olsf%lsf)
call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',-1) ; deallocate (olsf%icon)
call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',-1) ; deallocate (olsf%octree)
endif
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! deallocate void
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! note that vo will be needed in the temperature calculations
! so if this is the last iteration (converge is true), we will not deallocate vo
if (.not.converge) then
call heap (threadinfo,'vo%node','main',size(vo%node),'int',-1) ; deallocate (vo%node)
call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',-1) ; deallocate (vo%leaf)
call heap (threadinfo,'vo%face','main',size(vo%face),'int',-1) ; deallocate (vo%face)
call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1) ; deallocate (vo%rtf)
call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1) ; deallocate (vo%ftr)
call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1) ; deallocate (vo%influid)
end if
if (converge) then
do is=1,params%ns
call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',-1) ; deallocate (surface0(is)%r)
call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',-1) ; deallocate (surface0(is)%s)
call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',-1) ; deallocate (surface0(is)%x)
call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',-1) ; deallocate (surface0(is)%y)
call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',-1) ; deallocate (surface0(is)%z)
call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',-1) ; deallocate (surface0(is)%xn)
call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',-1) ; deallocate (surface0(is)%yn)
call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',-1) ; deallocate (surface0(is)%zn)
enddo
endif
iter=iter+1
enddo
call heap_hop2_f (threadinfo)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! end of grid iterations
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! we update current_time and can only do it now because this is where the time step is known
current_time=current_time+params%dt
if (iproc.eq.0) write(8,*) 'current time =', current_time
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! solve for temperature
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (iproc.eq.0) then
write(*,*) '-----------------------------------------------------------------------'
write(8,*) '-----------------------------------------------------------------------'
call show_time (total,step,inc,1,'Temperature calculations $')
write (8,*) 'Start of temperature calculations'
write(*,*) '-----------------------------------------------------------------------'
write(8,*) '-----------------------------------------------------------------------'
end if
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! transfers velocity and temperature solution from ov to osolve
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
do i=1,ov%nnode
osolve%u(i)=ov%unode(i)
osolve%v(i)=ov%vnode(i)
osolve%w(i)=ov%wnode(i)
!osolve%wpreiso(i)=ov%wnodepreiso(i)
enddo
if (params%calculate_temp) then
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
! build matrix arrays, allocate memory for wsmp
!------------------------------------------------------------------------|
!------------------------------------------------------------------------|
call show_time (total,step,inc,1,'wsmp setup1$')
n = vo%nnode * ndoft
!---[topology]-----
allocate (tpl(n))
tpl%nheightmax=27*ndoft
do i=1,n
allocate (tpl(i)%icol(tpl(i)%nheightmax),stat=err)
enddo
call heap (threadinfo,'tpl(:)%icol(:)','main',n*tpl(1)%nheightmax,'int',+1)
!---[topology]-----
call wsmp_setup1 (n,n_iproc_st,n_iproc_end,n_iproc,nz_loc,ldb,ndoft,vo,osolve,tpl,params,threadinfo,istep,iter)
allocate(iproc_col(n),stat=threadinfo%err) ; call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',+1)
allocate(ja(nz_loc),stat=threadinfo%err) ; call heap (threadinfo,'ja','main',size(ja),'int',+1)
allocate(ia(n_iproc+1),stat=threadinfo%err) ; call heap (threadinfo,'ia','main',size(ia),'int',+1)
allocate(avals(nz_loc),stat=threadinfo%err) ; call heap (threadinfo,'avals','main',size(avals),'dp',+1)
allocate(b(ldb,nrhs),stat=threadinfo%err) ; call heap (threadinfo,'b','main',size(b),'dp',+1)
allocate(weightel(osolve%nleaves),stat=threadinfo%err) ; call heap (threadinfo,'weightel','main',size(weightel),'dp',+1)
call show_time (total,step,inc,1,'wsmp setup2$')
call wsmp_setup2 (n,n_iproc,n_iproc_st,n_iproc_end,nz_loc,iproc_col, &
ia,ja,ndoft,vo,osolve,tpl,params,threadinfo,istep,iter)
!---[topology]-----
call heap (threadinfo,'tpl(:)%icol(:)','main',n*tpl(1)%nheightmax,'int',-1)
do i=1,n
deallocate (tpl(i)%icol)
enddo
deallocate (tpl)
!---[topology]-----
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
! building the FEM matrix and rhs
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
call show_time (total,step,inc,1,'build system$')
call build_system_wsmp (n,n_iproc,n_iproc_st,nz_loc,iproc_col,ia,ja,ldb,nrhs,avals,b, &
params,osolve,ndoft,mat,vo,istep,iter,iter_nl,threadinfo, &
weightel)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! solve system with wsmp solver
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call solve_with_pwgsmp (n,nz_loc,n_iproc,n_iproc_st,n_iproc_end,ia,ja,ldb,&
nrhs,avals,b,params,osolve,ov,vo,threadinfo,ndoft,&
istep,iter,iter_nl)
!-------------------------------------------------------------------------!
!-------------------------------------------------------------------------!
! deallocate memory used by the solver and terminates the solver's job
!-------------------------------------------------------------------------!
!-------------------------------------------------------------------------!
call show_time (total,step,inc,1,'wsmp_cleanup$')
call heap (threadinfo,'ia','main',size(ia),'int',-1) ; deallocate(ia)
call heap (threadinfo,'ja','main',size(ja),'int',-1) ; deallocate(ja)
call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',-1) ; deallocate(iproc_col)
call heap (threadinfo,'avals','main',size(avals),'dp',-1) ; deallocate(avals)
call heap (threadinfo,'b','main',size(b),'dp',-1) ; deallocate(b)
call heap (threadinfo,'weightel','main',size(weightel),'dp',-1) ; deallocate(weightel)
else
if(iproc.eq.0) then
write(*,'(a,f15.7,a)') shift//'================================='
write(*,'(a)') shift//'skip temperature calculation'
write(*,'(a,f15.7,a)') shift//'================================='
end if
end if
!-------------------------------------------------------------------------!
!-------------------------------------------------------------------------!
call DoRuRe_temp_stats(params%doDoRuRe,istep,ov%nnode,ov%temp)
olsf%nnode=ov%nnode
olsf%nleaves=ov%nleaves
olsf%noctree=ov%noctree
allocate (olsf%octree(olsf%noctree),stat=threadinfo%err) ; call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
allocate (olsf%x(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
allocate (olsf%y(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
allocate (olsf%z(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err) ; call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)
call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)
olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
olsf%y(1:olsf%nnode)=ov%y(1:ov%nnode)
olsf%z(1:olsf%nnode)=ov%z(1:ov%nnode)
olsf%icon(:,1:olsf%nleaves)=ov%icon(:,1:ov%nleaves)
! write time step to Log file
if (iproc.eq.0) then
write (8,*) 'Current time step ',params%dt
write (8,*) 'Courant time step ',dtcourant
end if
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! move surfaces and apply erosion
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Advect the surfaces$')
do is=1,params%ns
write (ic(1:2),'(i2)') is
call show_time (total,step,inc,1,'Advect surface '//ic(1:2)//'$')
if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
call move_surface (surface(is),surface0(is),0,0,ov,params%dt,params,istep,is)
else
surface(is)%x=surface(1)%x
surface(is)%y=surface(1)%y
surface(is)%z=surface(1)%z
surface(is)%xn=surface(1)%xn
surface(is)%yn=surface(1)%yn
surface(is)%zn=surface(1)%zn
surface(is)%r=surface(1)%r
surface(is)%s=surface(1)%s
endif
! We apply erosion
if (material0.eq.0.and.params%erosion) then
call show_time (total,step,inc,1,'Erode surface '//ic(1:2)//'$')
if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
call erosion (surface(is),olsf,is,params)
else
if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
surface(is)%z=surface(1)%z
endif
end if
if (iproc.eq.0) write (8,*) 'Min-Max z surf ',is,':',minval(surface(is)%z),maxval(surface(is)%z)
allocate(surface(is)%u(surface(is)%nsurface))
allocate(surface(is)%v(surface(is)%nsurface))
allocate(surface(is)%w(surface(is)%nsurface))
enddo
!call interpolate_velocity_on_surface(params,surface,ov)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! compute isostasy, move surfaces again, update Moho displacement
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (params%isostasy) then
call show_time (total,step,inc,1,'Compute isostasy and adjust vertical velocity$')
!allocate(ov%wpreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wpreiso in main$')
call isostasy (params,weightel,ov,surface,mat,0,zi)
if (params%isobc) zi%zisodisp=zi%zisodisp-zi%zisoinc
do is=1,params%ns
if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
call move_surface (surface(is),surface0(is),0,1,ov,params%dt,params,istep,is)
else
surface(is)%x=surface(1)%x
surface(is)%y=surface(1)%y
surface(is)%z=surface(1)%z
surface(is)%xn=surface(1)%xn
surface(is)%yn=surface(1)%yn
surface(is)%zn=surface(1)%zn
surface(is)%r=surface(1)%r
surface(is)%s=surface(1)%s
endif
enddo
endif
call interpolate_velocity_on_surface(params,surface,ov)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! update strain, pressure, temperature on 3D cloud
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Update cloud fields$')
call update_cloud_fields (cl,ov,osolve,vo,params)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! advect cloud
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Advect the cloud$')
! call move_cloud (params,cl,ov%octree,ov%noctree,ov%unode,ov%vnode,ov%wnode,ov%nnode,ov%icon,ov%nleaves)
call move_cloud (params,cl,ov)
call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',-1) ; deallocate (olsf%x)
call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',-1) ; deallocate (olsf%y)
call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',-1) ; deallocate (olsf%z)
call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',-1) ; deallocate (olsf%lsf)
call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',-1) ; deallocate (olsf%icon)
call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',-1) ; deallocate (olsf%octree)
osolve%u=ov%unode
osolve%v=ov%vnode
osolve%w=ov%wnode
osolve%wiso=ov%wnodeiso
osolve%temp=ov%temp
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! write global output before deallocating arrays
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call show_time (total,step,inc,1,'Write output$')
call write_global_output (params,istep,iter-1,current_time,osolve,ov,vo,surface,cl,zi,'final')
call mpi_barrier (mpi_comm_world,ierr)
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! deallocate osolve and vo
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
do is=1,params%ns
deallocate(surface(is)%u)
deallocate(surface(is)%v)
deallocate(surface(is)%w)
end do
call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',-1) ; deallocate (osolve%octree)
call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1) ; deallocate (osolve%icon)
call heap (threadinfo,'osolve%x','main',size(osolve%x),'dp',-1) ; deallocate (osolve%x)
call heap (threadinfo,'osolve%y','main',size(osolve%y),'dp',-1) ; deallocate (osolve%y);
call heap (threadinfo,'osolve%z','main',size(osolve%z),'dp',-1) ; deallocate (osolve%z)
call heap (threadinfo,'osolve%lsf','main',size(osolve%lsf),'dp',-1) ; deallocate (osolve%lsf)
call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',-1) ; deallocate (osolve%kfix)
call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',-1) ; deallocate (osolve%iface)
call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',-1) ; deallocate (osolve%pressure)
call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',-1) ; deallocate (osolve%spressure)
call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',-1) ; deallocate (osolve%strain)
call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',-1) ; deallocate (osolve%kfixt)
call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1) ; deallocate (osolve%u)
call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1) ; deallocate (osolve%v)
call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1) ; deallocate (osolve%w)
call heap (threadinfo,'osolve%wiso','main',size(osolve%wiso),'dp',-1) ; deallocate (osolve%wiso)
call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1) ; deallocate (osolve%temp)
call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',-1) ; deallocate (osolve%crit)
call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1) ; deallocate (osolve%e2d)
call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',-1) ; deallocate (osolve%e3d)
call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',-1) ; deallocate (osolve%lode)
call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1) ; deallocate (osolve%eviscosity)
call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',-1) ; deallocate (osolve%is_plastic)
call heap (threadinfo,'vo%node','main',size(vo%node),'int',-1) ; deallocate (vo%node)
call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',-1) ; deallocate (vo%leaf)
call heap (threadinfo,'vo%face','main',size(vo%face),'int',-1) ; deallocate (vo%face)
call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1) ; deallocate (vo%rtf)
call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1) ; deallocate (vo%ftr)
call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1) ; deallocate (vo%influid)
!call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1) ; deallocate (ov%wpreiso)
call show_time (total,step,inc,1,'End of time step$')
! end of big time loop
istep=istep+1
enddo
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
! end of time stepping
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
call heap_hop1_f (threadinfo)
call show_time (total,step,inc,1,'End of run$')
call io_DoRuRe2 (params,'close')
deallocate (cl%x)
deallocate (cl%y)
deallocate (cl%z)
deallocate (cl%x0)
deallocate (cl%y0)
deallocate (cl%z0)
deallocate (cl%strain)
deallocate (cl%lsf0)
deallocate (cl%temp)
deallocate (cl%press)
deallocate (cl%tag)
call heap (threadinfo,'ov%octree','main',size(ov%octree),'int',-1) ; deallocate (ov%octree)
call heap (threadinfo,'ov%x','main',size(ov%x),'dp',-1) ; deallocate (ov%x)
call heap (threadinfo,'ov%y','main',size(ov%y),'dp',-1) ; deallocate (ov%y)
call heap (threadinfo,'ov%z','main',size(ov%z),'dp',-1) ; deallocate (ov%z)
call heap (threadinfo,'ov%unode','main',size(ov%unode),'dp',-1) ; deallocate (ov%unode)
call heap (threadinfo,'ov%vnode','main',size(ov%vnode),'dp',-1) ; deallocate (ov%vnode)
call heap (threadinfo,'ov%wnode','main',size(ov%wnode),'dp',-1) ; deallocate (ov%wnode)
call heap (threadinfo,'ov%wnodeiso','main',size(ov%wnodeiso),'dp',-1); deallocate (ov%wnodeiso)
call heap (threadinfo,'ov%icon','main',size(ov%icon),'int',-1) ; deallocate (ov%icon)
call heap (threadinfo,'ov%temp','main',size(ov%temp),'dp',-1) ; deallocate (ov%temp)
call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1)
deallocate (ov%temporary_nodal_pressure)
call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1)
deallocate (ov%whole_leaf_in_fluid)
do is=1,params%ns
call heap (threadinfo,'surface(is)%x','main',size(surface(is)%x),'dp',-1) ; deallocate (surface(is)%x)
call heap (threadinfo,'surface(is)%y','main',size(surface(is)%y),'dp',-1) ; deallocate (surface(is)%y)
call heap (threadinfo,'surface(is)%z','main',size(surface(is)%z),'dp',-1) ; deallocate (surface(is)%z)
call heap (threadinfo,'surface(is)%xn','main',size(surface(is)%xn),'dp',-1) ; deallocate (surface(is)%xn)
call heap (threadinfo,'surface(is)%yn','main',size(surface(is)%yn),'dp',-1) ; deallocate (surface(is)%yn)
call heap (threadinfo,'surface(is)%zn','main',size(surface(is)%zn),'dp',-1) ; deallocate (surface(is)%zn)
call heap (threadinfo,'surface(is)%r','main',size(surface(is)%r),'dp',-1) ; deallocate (surface(is)%r)
call heap (threadinfo,'surface(is)%s','main',size(surface(is)%s),'dp',-1) ; deallocate (surface(is)%s)
call heap (threadinfo,'surface(is)%icon','main',size(surface(is)%icon),'int',-1) ; deallocate (surface(is)%icon)
end do
deallocate (surface)
deallocate (surface0)
deallocate (mat)
!if (params%isobc) deallocate (zisodisp)
if (params%isobc) then
deallocate (zi%zisodisp)
deallocate (zi%zisoinc)
endif
deallocate (params%materialn)
call heap_final (threadinfo)
close (threadinfo%Logunit)
close (threadinfo%mem_heap_unit)
call mpi_finalize (ierr)
end