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 (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-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')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     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

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