DOUAR.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

program DOUAR

use threads
use definitions
!use mpi

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 (bc_definition) bcdef
type (nest_info) nest

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-WSMP v0.1 (2011-03-23)'
!if (iproc.eq.0) write(*,*)  'This is DOUAR-WSMP v0.1 (2011-03-23), 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')

params%nmpi=nproc
if (iproc.eq.0) write(*,'(a,i3,a)')  'This is DOUAR-WSMP v0.1 (2011-03-23), using ', nproc, ' processors.'

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     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(*,'(a)') 'Program called with no argument'
else
   call getarg (1,params%infile)
   if (iproc.eq.0) then
      write(*,'(a,a)') 'program called with input file ',params%infile
   end if
end if

call write_streepjes(6,2)
call show_time (total,step,inc,0,'Start of Computations$')
call show_time (total,step,inc,1,'Reading Input$')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     read debug level
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call read_controlling_parameters (params,threadinfo,'debug')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     allocate threadinfo, open MPI log and memory heap files
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call int_to_char(ciproc,3,iproc)
if (params%debug.gt.1) then
  threadinfo%Logunit=1000+iproc
  open (unit=threadinfo%Logunit,file='./DEBUG/mpilogs/Log_'//ciproc//'.dat',status='replace')
  write(threadinfo%Logunit,*) 'This is DOUAR-WSMP v0.1 (2011-03-23)'
  write(threadinfo%Logunit,'(a,i3)') 'Log file of mpi process',iproc
  write(threadinfo%Logunit,'(a16,i3)') 'debug',params%debug
  call heap_init(threadinfo,2000+iproc,'./DEBUG/mpilogs/mem_heap_'//ciproc//'.dat')
endif

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     read controlling parameters
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call read_controlling_parameters (params,threadinfo,'main')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     open and prepare output files
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call io_DoRuRe2 (params,'init')

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,nest,bcdef)

if (params%irestart==0) then