Particle Density calculation only within domain

find_leaf crashes if particle lies outside domain too far.

src/update_tracking_cloud.f90

@@ -47,11 +47,12 @@ double precision dt,current_time
 !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
 !------------------------------------------------------------------------------|
 
+logical,dimension(:),allocatable::indomain,indomainp
 integer iproc,nproc,ierr,err,iter
 integer i,j,k,trgridshape,gridstep,tracksi,ie
 integer ntrpx,ntrpy,ntrpz
 integer nsurfd,nreset,leaf,level,locc
-integer,dimension(:),allocatable::trdensloss
+integer,dimension(:),allocatable::trdensloss,trdenslossp
 integer,dimension(:,:),allocatable::trackswap
 double precision u,v,w,rat
 double precision x0,y0,z0,xi,yi,zi,x,y,z,dxyz
@@ -68,6 +69,8 @@ character(4) cistep
 call mpi_comm_size (mpi_comm_world,nproc,ierr)
 call mpi_comm_rank (mpi_comm_world,iproc,ierr)
 
+!write (*,'(i2,a,i6)') iproc,': Starting update_trcloud, np = ',trcl%np
+
 !------------------------------------------------------------------------------|
 ! Creating Tracking Cloud at first timestep
 !------------------------------------------------------------------------------|
@@ -174,6 +177,8 @@ end if
 ! -- taken form update_cloud_fields.f90                                        |
 !------------------------------------------------------------------------------|
 
+!write (*,'(i2,a,i6)') iproc,': start interpolation, np = ',trcl%np
+
 !Get upper surface level set function
 allocate (lsf(osolve%nnode),stat=err) 
 if (err.ne.0) call stop_run ('Error alloc lsf in update_cloud_fields')
@@ -221,6 +226,8 @@ do i = 1,trcl%np
     if (trcl%lsf0(i) > 0) nsurfd = nsurfd+1
 end do
 
+!write (*,'(i2,a,i4)') iproc,': interpolation done, nsurfd = ',nsurfd
+
 savetracks: if (nsurfd>0) then
     !increase tr acks array for newly surfaced particles
     tracksize: if (trcl%ntracks > 0) then
@@ -264,6 +271,8 @@ end if savetracks
 ! The log of all surfaced tracks is saved.                                     |
 !------------------------------------------------------------------------------|
 
+!write (*,'(i2,a6)') iproc,': savetracks done, write to disk'
+
 if (iproc.eq.0 .AND. mod(istep-params%tr_activation,params%tr_savstep)==0) then
     call int_to_char (cistep,4,istep)
     open (20,file='TROUT/tracking_t'//cistep//'.bin',status='replace',form='unformatted')
@@ -295,6 +304,7 @@ end if
 ! taken from move_cloud.f90                                                    |
 !------------------------------------------------------------------------------|
 
+!write (*,'(i2,a6)') iproc,': Start moving particles'
 allocate (xbuf(trcl%np),stat=err) 
 if (err.ne.0) call stop_run ('Error alloc xbuf in update_trcloud')
 allocate (ybuf(trcl%np),stat=err) 
@@ -361,6 +371,8 @@ deallocate (zbuf)
 ! Check whether density referene array is still valid. If not, recalculate from 
 ! initial particle positions.
 
+!write (*,'(i2,a)') iproc,': start reset particles'
+
 if (trcl%ne_init .ne. osolve%nleaves) then
     if (iproc == 0) write (*,*) 'Start calculation reference array'
     !get rid of old array
@@ -380,23 +392,51 @@ if (trcl%ne_init .ne. osolve%nleaves) then
     enddo
 end if
 
+
 !calculate the decrease in each element
+allocate (trdenslossp(trcl%ne_init),stat=err) 
+if (err.ne.0) call stop_run ('error in alloc trdenslossp in update_trcloud')
 allocate (trdensloss(trcl%ne_init),stat=err) 
 if (err.ne.0) call stop_run ('error in alloc trdensloss in update_trcloud')
+
+allocate (indomain(trcl%ne_init),stat=err) 
+if (err.ne.0) call stop_run ('error in alloc indomain in update_trcloud')
+allocate (indomainp(trcl%ne_init),stat=err) 
+if (err.ne.0) call stop_run ('error in alloc indomainp in update_trcloud')
+
 trdensloss=0
-do i=1,trcl%np
-   call octree_find_leaf (osolve%octree,osolve%noctree, & 
-                          trcl%x(i),trcl%y(i),trcl%z(i), &
-                          leaf,level,locc,x0,y0,z0,dxyz)
-   trdensloss(leaf)=trdensloss(leaf)+1
+trdenslossp=0
+
+indomain=.false.
+indomainp=.false.
+
+do i=1+iproc,trcl%np,nproc
+   !write (*,'(a,3F10.6)') 'position', trcl%x(i), trcl%y(i), trcl%z(i)
+   indomainp(i) = 0.0d0<trcl%x(i) .AND. trcl%x(i)<1.d0 .AND. &
+       0.0d0<trcl%y(i) .AND. trcl%y(i)<1.d0 .AND. &
+       0.0d0<trcl%z(i) .AND. trcl%lsf0(i)<0.d0
+   if (indomainp(i)) then
+       call octree_find_leaf (osolve%octree,osolve%noctree, & 
+                              trcl%x(i),trcl%y(i),trcl%z(i), &
+                              leaf,level,locc,x0,y0,z0,dxyz)
+       !write (*,*) 'np, leaf, trdensloss', i, leaf, trdensloss(leaf)+1
+       trdensloss(leaf)=trdensloss(leaf)+1
+   end if
 enddo
-trdensloss = trcl%trdens_init-trdensloss
+
+call mpi_allreduce (trdensloss,trdenslossp,trcl%np,mpi_integer,mpi_sum,mpi_comm_world,ierr)
+call mpi_allreduce (indomainp,indomain,trcl%np,mpi_logical,mpi_lor,mpi_comm_world,ierr)
+
+trdensloss = trcl%trdens_init-trdenslossp
+
+deallocate (indomainp)
+deallocate (trdenslossp)
+
+!write (*,'(i2,a)') iproc,': trdensloss calculated, start reset routine'
 
 nreset=0
 do i=1,trcl%np
-    if (trcl%x(i)>1.d0 .OR. trcl%x(i)<0.d0 .OR. &
-        trcl%y(i)>1.d0 .OR. trcl%y(i)<0.d0 .OR. &
-        trcl%lsf0(i)>0.d0 .OR. trcl%z(i)<0.d0) then
+    if (.NOT. indomain(i)) then
 
         !find leaf with hightes loss of particles
         ie = maxloc(trdensloss,DIM=1)
@@ -426,6 +466,7 @@ if (iproc==0) write (*,'(a,i4,a,i4,a)') 'update_trcloud:',nsurfd, &
     ' particles surfaced and in total ',nreset,' particles reset'
 
 deallocate (trdensloss)
+deallocate (indomain)
 
 return
 end