diff --git a/Code/.gitignore b/Code/.gitignore new file mode 100644 index 0000000..63a7748 --- /dev/null +++ b/Code/.gitignore @@ -0,0 +1 @@ +*.mod diff --git a/Code/Makefile b/Code/Makefile index 3082aed..ed8c06c 100644 --- a/Code/Makefile +++ b/Code/Makefile @@ -60,7 +60,7 @@ clean-exec: OBJS = params.o grids.o levels.o fourier.o forces.o $(PARALLEL_SKY) \ inout.o coulomb.o trivial.o densities.o fragments.o twobody.o \ energies.o static.o meanfield.o dynamic.o pairs.o moment.o \ - main3d.o user.o external.o ylm.o + main3d.o user.o external.o ylm.o abso_bc.o # LINKERFLAGS_SKY = $(COMPILERFLAGS_SKY) # @@ -124,7 +124,9 @@ static.o : static.f90 params.o densities.o grids.o \ dynamic.o : dynamic.f90 params.o densities.o grids.o trivial.o \ levels.o moment.o energies.o inout.o meanfield.o $(PARALLEL_SKY) \ - twobody.o external.o Makefile + twobody.o external.o abso_bc.o Makefile + +abso_bc.o : abso_bc.f90 params.o grids.o levels.o densities.o meanfield.o $(PARALLEL_SKY) Makefile main3d.o : main3d.f90 params.o fourier.o forces.o densities.o meanfield.o \ levels.o grids.o fragments.o $(PARALLEL_SKY) dynamic.o static.o coulomb.o user.o \ diff --git a/Code/abso_bc.f90 b/Code/abso_bc.f90 new file mode 100644 index 0000000..f52baed --- /dev/null +++ b/Code/abso_bc.f90 @@ -0,0 +1,535 @@ +!------------------------------------------------------------------------------ +! MODULE: abso_bc +!------------------------------------------------------------------------------ +! DESCRIPTION: +!> @brief +!!Package containing absorbing boundary conditions and +!!subsequent analysis of observables from electron emission. +!!One may also consider to perform the cumulation of total +!!absorbed density within the absorbing loop. +!------------------------------------------------------------------------------ +MODULE abso_bc + USE Params + USE Grids, ONLY: dx,dy,dz,x,y,z,wxyz + USE Levels, ONLY: nstmax,psi,npsi,npmin,wocc,isospin,nprot,nneut + USE Meanfield + USE Parallel + IMPLICIT NONE + SAVE + !*********************************************************************** + ! * + ! General variables for the absorbing bounds * + ! * + !*********************************************************************** + INTEGER,PARAMETER :: nocc=1 !< degeneracy of states (outdated) + INTEGER :: ispherabso=0 !< switch to spherical absorbing bounds + INTEGER :: iangabso=0 + INTEGER :: ipes=0 + INTEGER :: ifabsorbital=0 !< switch to orbit-wise accumulation + INTEGER :: nangtheta=0 !< number of reference points for outgoing spectra + INTEGER :: nangphi=0 + INTEGER :: jescmaskorb=0 + INTEGER :: jescmask=10 + REAL(db) :: powabs=0.0175D0 !< power of absorbing mask + + REAL(db),PRIVATE :: xango + REAL(db),PRIVATE :: yango + REAL(db),PRIVATE :: zango + INTEGER,PRIVATE :: ix + INTEGER,PRIVATE :: iy + INTEGER,PRIVATE :: iz + INTEGER,PRIVATE :: nst + INTEGER,PRIVATE :: iq + INTEGER,PRIVATE :: nta + INTEGER,PRIVATE :: ita + INTEGER,PRIVATE :: nabso=0 ! @brief +!!Apply absorbing bounds, optionally accumulate absorbed density. +!--------------------------------------------------------------------------- + SUBROUTINE absbc(nabsorb,it,nt,time) + INTEGER,INTENT(in) :: it,nt,nabsorb + REAL(db),INTENT(in) :: time + REAL(db) :: weight + INTEGER :: iqa + + LOGICAL,SAVE :: firstcall = .TRUE. + REAL(db),ALLOCATABLE,DIMENSION(:,:,:,:) :: rhoabso_all + + !------------------------------------------------------------ + + nabso = nabsorb + IF(nabsorb <= 0) RETURN + + ! set calling parameters to internal common variables + ita = it + nta = nt + timea = time + + ! initializations + + IF (firstcall) THEN + IF(nabsorb>=NX/2.OR.nabsorb>=NY/2 .OR. nabsorb>=NZ/2) & + STOP 'ABSBC: problem with nabsorb' + CALL init_absbc() + firstcall = .FALSE. + END IF + + ! applying absorbing masks and accumulating absorption + + + + + ! apply mask function (and accumulate absorption per state) + ! and compute new density + + DO nst = 1,nstloc + iqa = isospin(globalindex(nst)) + weight = nocc*wocc(nst) + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + IF(tgridabso(ix,iy,iz)) THEN + rhoabso(ix,iy,iz,iqa) = rhoabso(ix,iy,iz,iqa) + & + weight*spherloss(ix,iy,iz)*( & + abs2(psi(ix,iy,iz,1,nst))+abs2(psi(ix,iy,iz,2,nst)) ) + psi(ix,iy,iz,1,nst)=absomask(ix,iy,iz)*psi(ix,iy,iz,1,nst) + psi(ix,iy,iz,2,nst)=absomask(ix,iy,iz)*psi(ix,iy,iz,2,nst) + END IF + END DO;END DO;END DO + END DO + + IF(ifabsorbital==1) THEN + DO nst = 1,nstloc + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + IF(tgridabso(ix,iy,iz)) THEN + rhoescmaskorb(ix,iy,iz,nst) = rhoescmaskorb(ix,iy,iz,nst) & + + spherloss(ix,iy,iz)*( & + abs2(psi(ix,iy,iz,1,nst))+abs2(psi(ix,iy,iz,2,nst)) ) + END IF + END DO;END DO;END DO + END DO + END IF + + + + ! print escaped nucleons + + IF(tmpi .AND. (MOD(ita,mprint)==1 .OR. MOD(ita,mrest)==0 .OR. & + MOD(ita,jescmask) == 0 .OR. ita == nta) ) THEN + ALLOCATE(rhoabso_all(NX,NY,NZ,2)) + CALL mpi_barrier (mpi_comm_world, mpi_ierror) + CALL mpi_allreduce(rhoabso,rhoabso_all,2*NX*NY*NZ, & + mpi_double_precision,mpi_sum, & + mpi_comm_world,mpi_ierror) + IF(MOD(ita,mprint)==1) CALL nescape(rhoabso_all) + CALL escmask(rhoabso_all) + ELSE + IF(MOD(ita,mprint)==1) CALL nescape(rhoabso) + CALL escmask(rhoabso) + END IF + + ! prints wavefunction on measuring points on file + + ! IF (ipes /= 0) CALL evalMP() ! not yet implemented + + IF(MOD(ita,mrest)==0) THEN + IF(tmpi) THEN + CALL wrtabso(ita,rhoabso_all) + ELSE + CALL wrtabso(ita,rhoabso) + END IF + WRITE(6,*) ' Augment restart file at it=',ita + ENDIF + IF(ALLOCATED(rhoabso_all)) DEALLOCATE(rhoabso_all) + + RETURN + END SUBROUTINE absbc +!--------------------------------------------------------------------------- +! DESCRIPTION: init_absbc +!> @brief +!!Initializes absorbing boundary conditions. +!--------------------------------------------------------------------------- + + SUBROUTINE init_absbc() + + NAMELIST /absobc/ ispherabso,iangabso,ipes,ifabsorbital, & + nangtheta,nangphi,jescmaskorb,jescmask,powabs + + !READ(5,absobc) + + ALLOCATE(tgridabso(NX,NY,NZ)) + ALLOCATE(absomask(NX,NY,NZ)) + ALLOCATE(spherloss(NX,NY,NZ)) + IF(ispherabso /= 0) THEN + CALL init_spherabso() + ELSE + CALL init_abso() + END IF + IF(ifabsorbital==1) THEN + ALLOCATE(rhoescmaskorb(NX,NY,NZ,nstloc)) + IF(.NOT.trestart) rhoescmaskorb = 0D0 + END IF + ALLOCATE(rhoabso(NX,NY,NZ,2)) + IF(.NOT.trestart) rhoabso = 0D0 + ! IF (ipes /= 0) CALL init_MP() ! not yet implemented + IF(trestart) THEN + CALL readabso() + WRITE(6,*) ' Reading abso-info from restart file' + ENDIF + + END SUBROUTINE init_absbc + + + + + REAL(db) FUNCTION dist_min() + + ! determines minimum radius for last non-absorbing point + + REAL(db) :: rmin + + xango=0D0 + yango=0D0 + zango=0D0 + + rmin = abs(x(NX-nabso)-xango) + rmin = min(rmin,abs(x(nabso+1)-xango)) + rmin = min(rmin,abs(y(NY-nabso)-yango)) + rmin = min(rmin,abs(y(nabso+1)-yango)) + rmin = min(rmin,abs(z(NZ-nabso)-zango)) + rmin = min(rmin,abs(z(nabso+1)-zango)) + + dist_min = rmin + + + RETURN + END FUNCTION dist_min +!--------------------------------------------------------------------------- +! DESCRIPTION: init_sphereabso +!> @brief +!!Initializes sphericalabsorbing boundary conditions. +!--------------------------------------------------------------------------- + SUBROUTINE init_spherabso() + + ! Initializes mask function for spherical boundary conditions + + REAL(db) :: dmin1,dmin2,bcrad,dmin12,dmin22,dist,dist2,cosact + + !------------------------------------------------------------ + + bcrad = nabso*dx + dmin1 = dist_min() + IF (DMIN1 < 0.) STOP 'Error in abso: dmin1<0' + dmin2 = dmin1+bcrad + dmin22=dmin2**2 + dmin12=DMIN1**2 + + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + dist2 = x(ix)**2+y(iy)**2+z(iz)**2 + IF (dist2 <= dmin12) THEN + absomask(ix,iy,iz)=1D0 + tgridabso(ix,iy,iz)=.false. + ELSE IF (dist2 > dmin22) THEN + absomask(ix,iy,iz)=0D0 + tgridabso(ix,iy,iz)=.true. + ELSE + dist = MAX(1D-20,SQRT(dist2)) + cosact = COS((dist-dmin1)*0.5D0*pi/bcrad) + IF(cosact > 0D0) THEN + absomask(ix,iy,iz)= cosact**powabs + ELSE + absomask(ix,iy,iz)= 0D0 + END IF + tgridabso(ix,iy,iz)=.true. + END IF + spherloss(ix,iy,iz) = 1D0-absomask(ix,iy,iz)**2 + END DO;END DO;END DO + + RETURN + END SUBROUTINE init_spherabso +!--------------------------------------------------------------------------- +! DESCRIPTION: init_abso +!> @brief +!!Initializes normal absorbing boundary conditions. +!--------------------------------------------------------------------------- + SUBROUTINE init_abso() + USE Parallel + + ! Initializes mask for rectangular absorbing boundaries conditions + + REAL(db),ALLOCATABLE,DIMENSION(:) :: xmask,ymask,zmask + + LOGICAL :: wflagabs=.true. + + ALLOCATE(xmask(NX)) + ALLOCATE(ymask(NY)) + ALLOCATE(zmask(NZ)) + + ! prepare mask functions in each direction separately + + zmask = 1D0 + DO iz = 1,nabso + zmask(iz) = COS(pi*0.5D0*(nabso+1.0-iz)/nabso)**powabs + zmask(NZ+1-iz) = zmask(iz) + END DO + + ymask = 1D0 + DO iy = 1,nabso + ymask(iy) = COS(pi*0.5D0*(nabso+1.0-iy)/nabso)**powabs + ymask(NY+1-iy) = ymask(iy) + END DO + + xmask = 1D0 + DO ix = 1,nabso + xmask(ix) = COS(pi*0.5D0*(nabso+1.0D0-ix)/nabso)**powabs + xmask(NX+1-ix) = xmask(ix) + END DO + + IF(mpi_myproc == 0 .AND. wflagabs) THEN + WRITE(6,'(a)') ' ZMASK:' + WRITE(6,'(1x,5(1pg12.4))') zmask + WRITE(6,'(a)') ' YMASK:' + WRITE(6,'(1x,5(1pg12.4))') ymask + WRITE(6,'(a)') ' XMASK:' + WRITE(6,'(1x,5(1pg12.4))') xmask + END IF + + + ! compose to one mask function on all grid points + + FORALL(ix=1:NX,iy=1:NY,iz=1:NZ) + absomask(ix,iy,iz)=xmask(ix)*ymask(iy)*zmask(iz) + tgridabso(ix,iy,iz)= absomask(ix,iy,iz) < 0.999999999999D0 + spherloss(ix,iy,iz) = 1D0-absomask(ix,iy,iz)**2 + END FORALL + + RETURN + END SUBROUTINE init_abso + + + + !-----escmask--------------------------------------------------------- + + SUBROUTINE escmask(rhoabsoprint) + USE Parallel + + ! print collected information on escaping nucleons + + REAL(db),INTENT(IN),DIMENSION(:,:,:,:) :: rhoabsoprint + + CHARACTER(LEN=5) :: scratchline + CHARACTER(LEN=5) :: readstring + CHARACTER(LEN=4) :: str + + !-------------------------------------------------------------------- + + + + + IF(ifabsorbital==1) THEN + IF (MOD(ita,jescmaskorb) == 0 .OR. ita==nta) THEN + DO nst=1,nstmax + IF(nstmax > 9999) STOP 'ERROR: Too many states for escmaskOrb' + IF(mpi_myproc == node(nst)) THEN + WRITE(scratchline,'(i5)') 10000+nst + READ(scratchline,'(a)') readstring + str = readstring(2:5) + OPEN(588,STATUS='unknown', FILE='escmaskOrb.'//str) + WRITE(588,'(a,f12.4,a,i5,/a,3i5/a)') & + '# distribution of escaped nucleon at time=',timea,& + ' for state=',nst, & + '#nx,ny,nz:',NX,NY,NZ, & + '# x y z loss of density' + + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + WRITE(588,'(3f12.4,1e17.7)') & + x(ix),y(iy),z(iz),rhoescmaskorb(ix,iy,iz,node(nst)) + END DO; END DO; END DO + + CLOSE(588) + END IF + END DO + END IF + END IF + + + IF(mpi_myproc == 0) THEN + IF (MOD(ita,jescmask) == 0 .OR. ita == nta) THEN + DO iq=1,2 + IF(iq == 1)THEN + OPEN(589,STATUS='unknown',FILE='escmask.1') + WRITE(589,'(a,f12.4/a,3i5/a)') & + '# total distribution of escaped protons at time=',timea, & + '#nx,ny,nz:',NX,NY,NZ, & + '# x y z loss of density' + ELSE + OPEN(589,STATUS='unknown',FILE='escmask.2') + WRITE(589,'(a,f12.4/a,3i5/a)') & + '# total distribution of escaped neutrons at time=',timea, & + '#nx,ny,nz:',NX,NY,NZ, & + '# x y z loss of density' + END IF + + DO iz=1,NZ; DO iy=1,NY; DO ix=1,NX + WRITE(589,'(3f12.4,1e17.7)') & + x(ix),y(iy),z(iz),rhoabsoprint(ix,iy,iz,iq) + END DO; END DO; END DO + + CLOSE(589) + END DO + END IF + END IF + + + RETURN + END SUBROUTINE escmask + + + + + SUBROUTINE nescape(rhoabsoprint) + + ! Compute and print total number of escaped nucleons + + REAL(db),INTENT(IN),DIMENSION(:,:,:,:) :: rhoabsoprint + + LOGICAL,SAVE :: firstcall = .TRUE. + LOGICAL :: texist + REAL(db) :: aprot,aneut + + !------------------------------------------------------------ + + IF(firstcall) THEN + INQUIRE(file='nescape.res',exist=texist) + IF(.NOT.texist) THEN + OPEN(unit=scratch,file='nescape.res') + WRITE(scratch,'(a)') & + '# number of escaped nucleons ', & + '# time lost protons lost neutrons ' + CLOSE(unit=scratch) + END IF + firstcall = .FALSE. + END IF + + aprot = 1D0*nprot !SUM(fcharg(1:2)) + aneut = 1D0*nneut !SUM(fmass(1:2))-aprot + IF(mpi_myproc == 0) THEN + OPEN(unit=scratch,file='nescape.res',POSITION='APPEND') + WRITE(scratch,'(1x,f10.4,4(1pg20.10))') timea, & + aprot-wxyz*sum(rho(:,:,:,1)),aneut-wxyz*sum(rho(:,:,:,2)), & + wxyz*sum(rhoabsoprint(:,:,:,1)),wxyz*sum(rhoabsoprint(:,:,:,2)) + CLOSE(unit=scratch) + END IF + + RETURN + END SUBROUTINE nescape + + !------------------------------------------------------------ + + + SUBROUTINE wrtabso(iteration,rhoabsoprint) + ! + USE Parallel, ONLY: mpi_myproc,mpi_nprocs,nstloc + IMPLICIT NONE + INTEGER,INTENT(IN) :: iteration + INTEGER :: nst,lenstr + CHARACTER(120) :: rsfp + REAL(db),INTENT(IN),DIMENSION(:,:,:,:) :: rhoabsoprint + ! + !*********************************************************************** + ! * + ! complement wffile by info on accumulated absorption * + ! * + !*********************************************************************** + ! + ! first write general information only with processor 0 + ! + IF(wffile=='none'.OR.wffile=='NONE') THEN + WRITE(6,*) " wffile='NONE' --> no file written " + RETURN + ENDIF + lenstr = LEN_TRIM(wffile) + IF(mpi_myproc==0) THEN + OPEN(UNIT=scratch2,FILE=wffile(1:lenstr)//'.a',STATUS='REPLACE', & + FORM='UNFORMATTED') + WRITE(scratch2) rhoabsoprint + ENDIF + ! now if we are running with one processor, also write the wave functions + ! onto that file; otherwise open new individualized file + ! + IF(ifabsorbital==1) THEN + IF(mpi_nprocs>1) THEN + IF(mpi_myproc/=0) CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(rsfp,'(I3.3,''.'',A)') mpi_myproc,wffile(1:lenstr)//'.a' + OPEN(UNIT=scratch2,FILE=rsfp,STATUS='REPLACE', & + FORM='UNFORMATTED') + ENDIF + ! + DO nst = 1,nstloc + WRITE(scratch2) rhoescmaskOrb(:,:,:,nst) + ENDDO + END IF + CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(6,*) 'Restart file augmented at iteration ',iteration + END SUBROUTINE wrtabso + + SUBROUTINE readabso() + ! + USE Parallel, ONLY: mpi_myproc,mpi_nprocs,nstloc + IMPLICIT NONE + INTEGER :: nst,lenstr + CHARACTER(120) :: rsfp + ! + !*********************************************************************** + ! * + ! complement wffile by info on accumulated absorption * + ! * + !*********************************************************************** + ! + ! first write general information only with processor 0 + ! + IF(wffile=='none'.OR.wffile=='NONE') THEN + WRITE(6,*) " wffile='NONE' --> no file read " + RETURN + ENDIF + lenstr = LEN_TRIM(wffile) + IF(mpi_myproc==0) THEN + OPEN(UNIT=scratch2,FILE=wffile(1:lenstr)//'.a',STATUS='OLD', & + FORM='UNFORMATTED') + READ(scratch2) rhoabso + ENDIF + ! now if we are running with one processor, also write the wave functions + ! onto that file; otherwise open new individualized file + ! + IF(ifabsorbital==1) THEN + IF(mpi_nprocs>1) THEN + IF(mpi_myproc/=0) CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(rsfp,'(I3.3,''.'',A)') mpi_myproc,wffile(1:lenstr)//'.a' + OPEN(UNIT=scratch2,FILE=rsfp,STATUS='OLD', & + FORM='UNFORMATTED') + ENDIF + ! + DO nst = 1,nstloc + READ(scratch2) rhoescmaskOrb(:,:,:,nst) + ENDDO + END IF + CLOSE(UNIT=scratch2,STATUS='KEEP') + WRITE(6,*) 'Absorption data read from restart file' + END SUBROUTINE readabso + + PURE REAL(db) FUNCTION abs2(c) + COMPLEX(db),INTENT(IN) :: c + abs2 = REAL(c)**2 + AIMAG(c)**2 + RETURN + END FUNCTION abs2 + + +END MODULE abso_bc diff --git a/Code/dynamic.f90 b/Code/dynamic.f90 index 80c3d80..1c714c9 100644 --- a/Code/dynamic.f90 +++ b/Code/dynamic.f90 @@ -22,8 +22,10 @@ MODULE DYNAMIC USE Inout, ONLY: write_wavefunctions,write_densities, plot_density, & sp_properties,start_protocol USE External + USE abso_bc IMPLICIT NONE SAVE + INTEGER :: nabsorb=0 INTEGER :: nt !< the number of the final time step to be !! calculated. In case of a restart this is smaller than the total number of time steps. REAL(db) :: dt !< the physical time increment in units of fm/c. @@ -49,7 +51,7 @@ MODULE DYNAMIC !!parameters. !--------------------------------------------------------------------------- SUBROUTINE getin_dynamic - NAMELIST /dynamic/ nt,dt,mxpact,mrescm,rsep,texternal + NAMELIST /dynamic/ nt,dt,mxpact,mrescm,rsep,texternal,nabsorb READ(5,dynamic) IF(wflag) THEN WRITE(*,*) '***** Parameters for the dynamic calculation *****' @@ -58,6 +60,7 @@ SUBROUTINE getin_dynamic WRITE(*,'(A,F7.2,A)') ' The calculation stops at ',rsep, & ' fm fragment separation' WRITE(*,'(A,I3)') ' Power limit in operator expansion:',mxpact + WRITE(*,'(A,I3)') ' Number of absorbing points:',nabsorb ENDIF IF(texternal) CALL getin_external END SUBROUTINE getin_dynamic @@ -265,19 +268,31 @@ SUBROUTINE dynamichf current=0.0D0 sdens=0.0D0 sodens=0.0D0 - ! propagate to end of step, accumulate densities + ! propagate to end of step; accumulate densities only if no absorbing boundaries !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst,ps4) SCHEDULE(STATIC) & !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) DO nst=1,nstloc ps4=psi(:,:,:,:,nst) CALL tstep(isospin(globalindex(nst)),mxpact,ps4) - CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & + IF(nabsorb == 0) CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & ps4,rho,tau,current,sdens,sodens) psi(:,:,:,:,nst)=ps4 ENDDO !$OMP END PARALLEL DO ! sum up over nodes - IF(tmpi) CALL collect_densities + IF(nabsorb == 0 .AND. tmpi) CALL collect_densities + ! Apply absorbing boundary conditions and compute densities from masked wavefunctions + IF(nabsorb > 0) THEN + CALL absbc(nabsorb,iter,nt,time) + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(nst) SCHEDULE(STATIC) & + !$OMP REDUCTION(+:rho,tau,current,sdens,sodens) + DO nst=1,nstloc + CALL add_density(isospin(globalindex(nst)),wocc(globalindex(nst)), & + psi(:,:,:,:,nst),rho,tau,current,sdens,sodens) + ENDDO + !$OMP END PARALLEL DO + IF(tmpi) CALL collect_densities + END IF ! Step 4: eliminate center-of-mass motion if desired IF(mrescm/=0) THEN IF(MOD(iter,mrescm)==0) THEN