Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
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
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
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$')
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! extract material information to be passed to solver
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
increase_current_level=.false.
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 $')
Dave Whipp
committed
if (iproc.eq.0) write(8,*) 'Doing step ',istep
if (istep.gt.1) params%init_e2d=.false.
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
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
call show_time (total,step,inc,1,'refine surface$')
call refine_surface(params,surface(is),surface0(is),threadinfo,nadd,nrem,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+nrem)
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)
Loading
Loading full blame...