ppm_pse/0000755000175000017500000000000011455632530011203 5ustar omaromarppm_pse/pse_module_realloc.f0000644000175000017500000001222511445416205015207 0ustar omaromar !------------------------------------------------------------------------- ! Module : pse_module_realloc !------------------------------------------------------------------------- ! ! Purpose : This module contains subroutines to dynamically ! reallocate arrays. ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_module_realloc.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:26 ivos ! CBL version of the PPM library ! ! Revision 1.1 2005/04/27 01:19:11 ivos ! Initial commit. Tested. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum,Hirschengraben 84 ! CH-8092 Zurich,Switzerland !------------------------------------------------------------------------- MODULE pse_module_realloc IMPLICIT NONE !------------------------------------------------------------------------- ! Reallocate a 1d or 2d array to a new length,preserving its ! contents. !------------------------------------------------------------------------- INTERFACE reallocate MODULE PROCEDURE reallocate_rv,reallocate_rm,reallocate_iv, & reallocate_im,reallocate_hv END INTERFACE CONTAINS !------------------------------------------------------------------------- ! The source !------------------------------------------------------------------------- ! for real vectors function reallocate_rv(p,n) USE pse_global IMPLICIT NONE REAL(MK),DIMENSION(:),POINTER :: p,reallocate_rv INTEGER,INTENT(IN) :: n INTEGER :: nold,ierr ALLOCATE(reallocate_rv(n),STAT=ierr) IF (ierr .NE. 0) THEN CALL pwrite(rank,'pse_module_realloc','Not enough memory',ierr) RETURN ENDIF IF (.NOT. ASSOCIATED(p)) RETURN nold = SIZE(p) reallocate_rv(1:MIN(nold,n))=p(1:MIN(nold,n)) DEALLOCATE(p) END function reallocate_rv ! for integer vectors function reallocate_iv(p,n) USE pse_global IMPLICIT NONE INTEGER,DIMENSION(:),POINTER :: p,reallocate_iv INTEGER,INTENT(IN) :: n INTEGER :: nold,ierr ALLOCATE(reallocate_iv(n),STAT=ierr) IF (ierr .NE. 0) THEN CALL pwrite(rank,'pse_module_realloc','Not enough memory',ierr) RETURN ENDIF IF (.NOT. ASSOCIATED(p)) RETURN nold = SIZE(p) reallocate_iv(1:MIN(nold,n))=p(1:MIN(nold,n)) DEALLOCATE(p) END function reallocate_iv ! for character vectors function reallocate_hv(p,n) USE pse_global IMPLICIT NONE CHARACTER(1),DIMENSION(:),POINTER :: p,reallocate_hv INTEGER,INTENT(IN) :: n INTEGER :: nold,ierr ALLOCATE(reallocate_hv(n),STAT=ierr) IF (ierr .NE. 0) THEN CALL pwrite(rank,'pse_module_realloc','Not enough memory',ierr) RETURN ENDIF IF (.NOT. ASSOCIATED(p)) RETURN nold = SIZE(p) reallocate_hv(1:MIN(nold,n))=p(1:MIN(nold,n)) DEALLOCATE(p) END function reallocate_hv ! for real matrices function reallocate_rm(p,n1,n2) USE pse_global IMPLICIT NONE REAL(MK),DIMENSION(:,:),POINTER :: p,reallocate_rm INTEGER,INTENT(IN) :: n1,n2 INTEGER :: nold,mold,ierr ALLOCATE(reallocate_rm(n1,n2),STAT=ierr) IF (ierr .NE. 0) THEN CALL pwrite(rank,'pse_module_realloc','Not enough memory',ierr) RETURN ENDIF IF (.NOT. ASSOCIATED(p)) RETURN nold = SIZE(p,1) mold = SIZE(p,2) reallocate_rm(1:MIN(nold,n1),1:MIN(mold,n2)) = & & p(1:MIN(nold,n1),1:MIN(mold,n2)) DEALLOCATE(p) END function reallocate_rm ! for integer matrices function reallocate_im(p,n1,n2) USE pse_global IMPLICIT NONE INTEGER,DIMENSION(:,:),POINTER :: p,reallocate_im INTEGER,INTENT(IN) :: n1,n2 INTEGER :: nold,mold,ierr ALLOCATE(reallocate_im(n1,n2),STAT=ierr) IF (ierr .NE. 0) THEN CALL pwrite(rank,'pse_module_realloc','Not enough memory',ierr) RETURN ENDIF IF (.NOT. ASSOCIATED(p)) RETURN nold = SIZE(p,1) mold = SIZE(p,2) reallocate_im(1:MIN(nold,n1),1:MIN(mold,n2)) = & & p(1:MIN(nold,n1),1:MIN(mold,n2)) DEALLOCATE(p) END function reallocate_im END MODULE pse_module_realloc ppm_pse/pse_time.f0000644000175000017500000000613211455355442013165 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_time !------------------------------------------------------------------------- ! ! Purpose : Retreives the current time (cpu). Uses either ! MPI_Wtime, etime or the f90 CPU_TIME intrinsic, ! based on the DEFINES that are set (see ! ppm_define.h): ! __MPI uses MPI_Wtime ! __ETIME uses etime ! none uses CPU_TIME ! ! Input : ! ! Input/output : ! ! Output : timing (F) current CPU clock time ! ! Routines : etime (C intrinsic) ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_time(timing) !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- IMPLICIT NONE #ifdef __MPI INCLUDE 'mpif.h' #endif !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- REAL(8), INTENT( OUT) :: timing !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER :: it,itr #ifdef __ETIME REAL(8) :: array(2) !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- REAL(8), EXTERNAL :: etime #endif #ifdef __MPI !------------------------------------------------------------------------- ! Call the MPI function MPI_Wtime !------------------------------------------------------------------------- !timing = MPI_Wtime() CALL SYSTEM_CLOCK(it,itr) timing = REAL(it,8)/REAL(itr,8) #else !------------------------------------------------------------------------- ! Call the C routine: etime to get the cpu time !------------------------------------------------------------------------- #ifdef __ETIME timing = etime(array) #else CALL CPU_TIME(timing) #endif #endif !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE RETURN END SUBROUTINE pse_time ppm_pse/pse_param.h0000644000175000017500000000245611445416205013330 0ustar omaromar !------------------------------------------------------------------------- ! Header file : pse_param !------------------------------------------------------------------------- ! ! Purpose : Parameters for the ppm_pse code. ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_param.h,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.1 2004/06/02 15:32:00 ivos ! Initial checkin. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Precision !------------------------------------------------------------------------- ! INTEGER, PARAMETER :: MK = ppm_kind_single INTEGER, PARAMETER :: MK = ppm_kind_double INTEGER, PARAMETER :: MAXCHAR = ppm_char ppm_pse/Make.ifc0000644000175000017500000000172011455632507012547 0ustar omaromar #---------------------------------------------------------------------- # PATHS (need to be adjusted) #---------------------------------------------------------------------- # PPM Library PPMDIR = ../libppm/ PPMNUMDIR = ../libppmnumerics/ # METIS Library METISDIR = /home/user/sw/metis/gcc # MPI Library MPIDIR = /opt/openmpi/1.4.3/intel #---------------------------------------------------------------------- # Compiler settings #---------------------------------------------------------------------- F77 = ifort FC = ifort CC = icc INC = -I$(PPMDIR)/include -I$(PPMNUMDIR)/include #DEBUG = -C -g FFLAGS = -FR -O3 $(INC) $(DEBUG) # change FFLAGS for debugging purpose #FFLAGS = -FR $(INC) $(DEBUG) CFLAGS = -O3 LD = $(F77) LDFLAGS = $(FFLAGS) -static CPP = cpp LIBS = -L$(PPMDIR)/lib -L$(PPMNUMDIR)/lib -lppm -lppmnumerics -L$(METISDIR)/lib -lmetis LDFLAGS = $(LIBS) ppm_pse/pse_read_ctrl.f0000644000175000017500000004515611445431763014177 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_read_ctrl !------------------------------------------------------------------------- ! ! Purpose : Read Control input file ! ! Input : cfile (C) name of control file to be read ! ! Input/output : ! ! Output : info (I) return status ! ! Routines : substart ! substop ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_read_ctrl(cfile,info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- CHARACTER(LEN=*) , INTENT(IN ) :: cfile INTEGER , INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 INTEGER :: i,j,idx INTEGER :: iUnit,istat INTEGER :: ilen,ios,iline,ilenctrl CHARACTER(LEN=MAXCHAR) :: cbuf,cvalue,carg LOGICAL :: lExist !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_read_ctrl',t0,info) !------------------------------------------------------------------------- ! Define unit number !------------------------------------------------------------------------- iUnit = 20 !------------------------------------------------------------------------- ! Check that the ctrl file exists !------------------------------------------------------------------------- ilenctrl = LEN_TRIM(cfile) IF (ilenctrl.LT.1) THEN CALL pwrite(rank,'pse_read_ctrl','No Ctrl file given!',info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Check that the ctrl file exists !------------------------------------------------------------------------- INQUIRE(FILE=cfile,EXIST=lExist) IF (.NOT.lExist) THEN WRITE(cbuf,'(2A)')'No such file: ',cfile(1:ilenctrl) CALL pwrite(rank,'pse_read_ctrl',cbuf,info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Open the file !------------------------------------------------------------------------- OPEN(iUnit,FILE=cfile,IOSTAT=ios,ACTION='READ') IF (ios.NE.0) THEN WRITE(cbuf,'(2A)')'Failed to open file: ',cfile(1:ilenctrl) CALL pwrite(rank,'pse_read_ctrl',cbuf,info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Scan file !------------------------------------------------------------------------- iline = 0 DO !---------------------------------------------------------------------- ! Increment line counter !---------------------------------------------------------------------- iline = iline + 1 !---------------------------------------------------------------------- ! Read line !---------------------------------------------------------------------- READ(iUnit,'(A)',END=100,ERR=200) cbuf ilen = LEN_TRIM(cbuf) !---------------------------------------------------------------------- ! Debug print of Ctrl file !---------------------------------------------------------------------- IF (debug.GT.1) THEN CALL pwrite(rank,'pse_read_ctrl',cbuf(1:ilen),info) ENDIF !---------------------------------------------------------------------- ! Skip comment or empty lines !---------------------------------------------------------------------- IF (ilen.GT.0.AND.cbuf(1:1).NE.'#') THEN !------------------------------------------------------------------- ! Remove space !------------------------------------------------------------------- j = 0 DO i=1,ilen IF (cbuf(i:i).NE.' ') THEN j = j + 1 cbuf(j:j) = cbuf(i:i) ENDIF ENDDO !------------------------------------------------------------------- ! Update length of string !------------------------------------------------------------------- ilen = j !------------------------------------------------------------------- ! Find position of = !------------------------------------------------------------------- idx = INDEX(cbuf,'=') !------------------------------------------------------------------- ! Exit if missing !------------------------------------------------------------------- IF (idx.LT.0) THEN WRITE(cbuf,'(A,I5)')'Incorrect line: ',iline CALL pwrite(rank,'pse_read_ctrl',cbuf,info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------- ! Get argument and value !------------------------------------------------------------------- carg = ADJUSTL(cbuf(1:idx-1)) cvalue = ADJUSTL(cbuf(idx+1:ilen)) !------------------------------------------------------------------- ! Convert to upper case !------------------------------------------------------------------- CALL UpperCase(carg,idx-1,info) !------------------------------------------------------------------- ! Get input file format !------------------------------------------------------------------- IF (carg.EQ.'INPUTFMT') THEN IF (cvalue(1:1).EQ.'F'.OR.cvalue(1:1).EQ.'f') THEN inputfmt = 'FORMATTED' ELSE inputfmt = 'UNFORMATTED' ENDIF !------------------------------------------------------------------- ! Get output file format !------------------------------------------------------------------- ELSEIF (carg.EQ.'OUTPUTFMT') THEN IF (cvalue(1:1).EQ.'F'.OR.cvalue(1:1).EQ.'f') THEN outputfmt = 'FORMATTED' ELSE outputfmt = 'UNFORMATTED' ENDIF !------------------------------------------------------------------- ! Get diagnostics file format !------------------------------------------------------------------- ELSEIF (carg.EQ.'DIAGNOSTICSFMT') THEN IF (cvalue(1:1).EQ.'F'.OR.cvalue(1:1).EQ.'f') THEN diagfmt = 'FORMATTED' ELSE diagfmt = 'UNFORMATTED' ENDIF !------------------------------------------------------------------- ! Get neighbor list type !------------------------------------------------------------------- ELSEIF (carg.EQ.'NEIGHLIST') THEN IF (cvalue(1:1).EQ.'V'.OR.cvalue(1:1).EQ.'v') THEN neighlist = 1 ELSEIF (cvalue(1:1).EQ.'C'.OR.cvalue(1:1).EQ.'c') THEN neighlist = 2 ELSEIF (cvalue(1:1).EQ.'N'.OR.cvalue(1:1).EQ.'n') THEN neighlist = 3 ENDIF !------------------------------------------------------------------- ! Get restart flag (for continuation run) !------------------------------------------------------------------- ELSEIF (carg.EQ.'RESTART') THEN READ(cvalue,'(L)',IOSTAT=ios,ERR=200) restart !------------------------------------------------------------------- ! Use lookup table? !------------------------------------------------------------------- ELSEIF (carg.EQ.'LOOKUP') THEN READ(cvalue,'(L)',IOSTAT=ios,ERR=200) lookup !------------------------------------------------------------------- ! Rebalance load based on true particle costs? !------------------------------------------------------------------- ELSEIF (carg.EQ.'REBALANCE') THEN READ(cvalue,'(L)',IOSTAT=ios,ERR=200) rebalance !------------------------------------------------------------------- ! Get debug level !------------------------------------------------------------------- ELSEIF (carg.EQ.'DEBUG') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) debug !------------------------------------------------------------------- ! Get problem dimensionality !------------------------------------------------------------------- ELSEIF (carg.EQ.'DIMENSIONS') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) ndim !------------------------------------------------------------------- ! Get core size of particles !------------------------------------------------------------------- ELSEIF (carg.EQ.'EPSILON') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) eps !------------------------------------------------------------------- ! Get cutoff for PP interactions !------------------------------------------------------------------- ELSEIF (carg.EQ.'CUTOFF') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) cutoff !------------------------------------------------------------------- ! Get symmetry falg !------------------------------------------------------------------- ELSEIF (carg.EQ.'SYMMETRY') THEN READ(cvalue,'(L)',IOSTAT=ios,ERR=200) symmetry !------------------------------------------------------------------- ! Get kernel type !------------------------------------------------------------------- ELSEIF (carg.EQ.'KERNEL') THEN READ(cvalue,'(I4)',IOSTAT=ios,ERR=200) kernel !------------------------------------------------------------------- ! Get final time step !------------------------------------------------------------------- ELSEIF (carg.EQ.'NSTEP') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) nstep !------------------------------------------------------------------- ! Get number of bins for triangles !------------------------------------------------------------------- ELSEIF (carg.EQ.'NBINS') THEN READ(cvalue,*,IOSTAT=ios,ERR=10) NBINY,NBINZ GOTO 20 10 READ(cvalue,*,IOSTAT=ios,ERR=200) NBINY NBINZ = NBINY GOTO 20 20 CONTINUE !------------------------------------------------------------------- ! Get initialization grid size !------------------------------------------------------------------- ELSEIF (carg.EQ.'NGRID') THEN READ(cvalue,*,IOSTAT=ios,ERR=30) Ngrid(1),Ngrid(2),Ngrid(3) GOTO 50 30 READ(cvalue,*,IOSTAT=ios,ERR=40) Ngrid(1),Ngrid(2) Ngrid(3) = Ngrid(1) GOTO 50 40 READ(cvalue,*,IOSTAT=ios,ERR=200) Ngrid(1) Ngrid(2) = Ngrid(1) Ngrid(3) = Ngrid(1) 50 CONTINUE !------------------------------------------------------------------- ! Get time step !------------------------------------------------------------------- ELSEIF (carg.EQ.'DT') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) dt !------------------------------------------------------------------- ! Get frequency of post processing !------------------------------------------------------------------- ELSEIF (carg.EQ.'FREQOUTPUT') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) freqoutput !------------------------------------------------------------------- ! Get frequency of diagnogstics !------------------------------------------------------------------- ELSEIF (carg.EQ.'FREQDIAG') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) freqdiag !------------------------------------------------------------------- ! Get frequency of restart file output !------------------------------------------------------------------- ELSEIF (carg.EQ.'FREQRESTART') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) freqrestart !------------------------------------------------------------------- ! Get modelarea !------------------------------------------------------------------- ELSEIF (carg.EQ.'DOMAIN') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) min_phys(1:ndim), & & max_phys(1:ndim) !------------------------------------------------------------------- ! Get case name !------------------------------------------------------------------- ELSEIF (carg.EQ.'CASENAME') THEN casename = ADJUSTL(cvalue) !------------------------------------------------------------------- ! Get result file !------------------------------------------------------------------- ELSEIF (carg.EQ.'OUTPUTFILE') THEN outputfile = ADJUSTL(cvalue) !------------------------------------------------------------------- ! Get log file !------------------------------------------------------------------- ELSEIF (carg.EQ.'LOGFILE') THEN logfile = ADJUSTL(cvalue) !------------------------------------------------------------------- ! Get restart file !------------------------------------------------------------------- ELSEIF (carg.EQ.'RESTARTFILE') THEN rstfile = ADJUSTL(cvalue) !------------------------------------------------------------------- ! Get input file !------------------------------------------------------------------- ELSEIF (carg.EQ.'INPUTFILE') THEN inputfile = ADJUSTL(cvalue) !------------------------------------------------------------------- ! Get abort file !------------------------------------------------------------------- ELSEIF (carg.EQ.'ABORTFILE') THEN abortfile = ADJUSTL(cvalue) !------------------------------------------------------------------- ! Get diagnostics file !------------------------------------------------------------------- ELSEIF (carg.EQ.'DIAGNOSTICSFILE') THEN diagfile = ADJUSTL(cvalue) !------------------------------------------------------------------- ! Get number of species !------------------------------------------------------------------- ELSEIF (carg.EQ.'SPECIES') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) NSPEC !------------------------------------------------------------------- ! Get diffusion constant of species !------------------------------------------------------------------- ELSEIF (carg(1:7).EQ.'ISODIFF') THEN ALLOCATE(diff(NSPEC), STAT=istat) IF (istat .NE. 0) THEN WRITE(cbuf,'(2A)') & & 'Unable to allocate diffusion cst vector' CALL pwrite(rank,'pse_read_ctrl',cbuf,info) info = -2 GOTO 9999 ENDIF READ(cvalue,*,IOSTAT=ios,ERR=200) diff(1:NSPEC) !------------------------------------------------------------------- ! Get bleached box coordinates !------------------------------------------------------------------- ELSEIF (carg.EQ.'ROI') THEN READ(cvalue,*,IOSTAT=ios,ERR=200) min_roi,max_roi ENDIF ENDIF ENDDO !------------------------------------------------------------------------- ! Something went wrong !------------------------------------------------------------------------- 200 CONTINUE WRITE(cbuf,'(A,I5,2A)') 'Error reading line: ',iline, & & ' of file: ',cfile(1:ilenctrl) ilen = LEN_TRIM(cbuf) CALL pwrite(rank,'pse_read_ctrl',cbuf(1:ilen),info) info = -1 GOTO 9999 !------------------------------------------------------------------------- ! End of file !------------------------------------------------------------------------- 100 CONTINUE !------------------------------------------------------------------------- ! Close file !------------------------------------------------------------------------- CLOSE(iUnit) !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_read_ctrl',t0,info) RETURN END SUBROUTINE pse_read_ctrl ppm_pse/UpperCase.f0000644000175000017500000000603011445416206013237 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : UpperCase !------------------------------------------------------------------------- ! ! Purpose : This routine converts a character string to upper case ! letters. ! ! Input : ilen (I) length of input string ! ! Input/output : string (C) input: string to be converted. ! output: upper case version of it ! ! Output : into (I) return status ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: UpperCase.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.1 2004/06/02 15:28:40 ivos ! Initial implementation. ! !------------------------------------------------------------------------- ! Jens Honore Walther ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE UpperCase(string,ilen,info) IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- CHARACTER(LEN=*), INTENT(INOUT) :: string INTEGER , INTENT(IN ) :: ilen INTEGER , INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER :: i,j INTEGER :: i1,i2,i3,iadd !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- info = 0 !------------------------------------------------------------------------- ! Convert to upper case !------------------------------------------------------------------------- i1 = IACHAR('a') - 1 i2 = IACHAR('z') + 1 i3 = IACHAR('A') iadd = i3 - i1 - 1 DO i=1,ilen j = IACHAR(string(i:i)) IF (j.GT.i1.AND.j.LT.i2) THEN string(i:i) = CHAR(j+iadd) ENDIF ENDDO !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE RETURN END SUBROUTINE UpperCase ppm_pse/pse_create_particles.f0000644000175000017500000002057011455412773015543 0ustar omaromar !------------------------------------------------------------------------- ! Function : pse_create_particles !------------------------------------------------------------------------- ! ! Purpose : Creates the particles and allocates particle ! memory. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status. 0 on success. ! ! Return value : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_create_particles(info) !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE #ifdef __MPI INCLUDE "mpif.h" #endif !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER , INTENT(OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 INTEGER :: isub,jsub,i,j,k,Npn INTEGER, DIMENSION(3) :: ndatamax REAL(MK) :: x,y,z CHARACTER(LEN=MAXCHAR) :: cbuf INTEGER :: local_roinp !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_create_particles',t0,info) !------------------------------------------------------------------------- ! Sort triangles into bins !------------------------------------------------------------------------- CALL pse_sort_triangles(info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_create_particles','Sorting triangles failed!',& & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Determine max extent of any sub !------------------------------------------------------------------------- ndatamax = 0 DO isub=1,topo%nsublist jsub = topo%isublist(isub) IF (ndata(1,jsub) .GT. ndatamax(1)) ndatamax(1) = ndata(1,jsub) IF (ndata(2,jsub) .GT. ndatamax(2)) ndatamax(2) = ndata(2,jsub) IF (ndata(3,jsub) .GT. ndatamax(3)) ndatamax(3) = ndata(3,jsub) ENDDO !------------------------------------------------------------------------- ! Get indicator function !------------------------------------------------------------------------- ALLOCATE(flags(ndatamax(1),ndatamax(2),ndatamax(3),topo%nsublist),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_create_particles', & & 'Computing indicator function failed!',info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Get indicator function !------------------------------------------------------------------------- CALL pse_get_indicator(info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_create_particles', & & 'Computing indicator function failed!',info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Determine local number of particles. The -1 is needed to prevent ! overlapping particles on internal sub boundaries. For freespace ! b.c. this should in principle check if the bdr is external and ! then go to the full ndata. !------------------------------------------------------------------------- Npart = 0 DO isub=1,topo%nsublist jsub = topo%isublist(isub) DO k=1,ndata(3,jsub)-1 DO j=1,ndata(2,jsub)-1 DO i=1,ndata(1,jsub)-1 IF (flags(i,j,k,isub) .EQ. -1) THEN Npart = Npart + 1 ENDIF ENDDO ENDDO ENDDO ENDDO !------------------------------------------------------------------------- ! Allocate particle memory. !------------------------------------------------------------------------- ALLOCATE(xp(ndim,Npart),wp(NSPEC,Npart),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_create_particles', & & 'Failed allocate particle memory.',info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Store particle positions and initial values (based on ROI) !------------------------------------------------------------------------- wp = 1.0_MK Npn = 0 roinp = 0 local_roinp = 0 DO isub=1,topo%nsublist jsub = topo%isublist(isub) DO k=1,ndata(3,jsub)-1 z = topo%min_subd(3,jsub) + REAL(k-1,MK)*dhz DO j=1,ndata(2,jsub)-1 y = topo%min_subd(2,jsub) + REAL(j-1,MK)*dhy DO i=1,ndata(1,jsub)-1 x = topo%min_subd(1,jsub) + REAL(i-1,MK)*dhx IF (flags(i,j,k,isub) .EQ. -1) THEN Npn = Npn + 1 xp(1,Npn) = x xp(2,Npn) = y xp(3,Npn) = z IF ((x.GT.min_roi(1)).AND.(x.LT.max_roi(1)).AND. & & (y.GT.min_roi(2)).AND.(y.LT.max_roi(2)).AND. & & (z.GT.min_roi(3)).AND.(z.LT.max_roi(3))) THEN wp(1:NSPEC,Npn) = 0.0_MK local_roinp = local_roinp + 1 ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO IF (Npn .NE. Npart) THEN CALL pwrite(rank,'pse_create_particles', & & 'MAJOR HAVOC. LOST PARTICLES !!!.',info) info = -1 GOTO 9999 ENDIF #ifdef __MPI ! calculate number of all particles in ROI CALL MPI_Allreduce(local_roinp,roinp,1,MPI_INTEGER,MPI_SUM,comm,info) #else roinp = local_roinp #endif !------------------------------------------------------------------------- ! Output !------------------------------------------------------------------------- !WRITE(cbuf,'(A,I9)') 'Number of particles on local processor: ',Npart !CALL pwrite(rank,'pse_create_particles',cbuf,info) !------------------------------------------------------------------------- ! Deallocate indicator function !------------------------------------------------------------------------- IF (ALLOCATED(flags)) DEALLOCATE(flags,STAT=info) !------------------------------------------------------------------------- ! Dump if needed !------------------------------------------------------------------------- IF (dumpparticles) THEN WRITE(cbuf,'(A,I2.2,A)') 'initial_',rank,'.out' OPEN(40,FILE=cbuf,STATUS='REPLACE',ACTION='WRITE') DO i=1,Npart WRITE(40,'(4E16.4)') xp(1:3,i), wp(1,i) ENDDO CLOSE(40) ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_create_particles',t0,info) RETURN END SUBROUTINE pse_create_particles ppm_pse/pse_init.f0000644000175000017500000002061311455355740013173 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_init !------------------------------------------------------------------------- ! ! Purpose : Initialize the PSE solver and read all input files. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status. 0 on success. ! ! Routines : pse_read_input ! pse_defaults ! pse_read_ctrl ! ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_init(info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE ppm_module_user_util USE ppm_module_io USE ppm_module_comp_part USE pse_global USE pse_module_io IMPLICIT NONE !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- #ifdef __MPI INCLUDE 'mpif.h' #endif !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER :: ilen,ios,tolexp REAL(MK), DIMENSION(3) :: bsize CHARACTER(LEN = MAXCHAR) :: cbuf REAL(8) :: t0 !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- info = 0 #ifdef __MPI !------------------------------------------------------------------------- ! Initialise MPI !------------------------------------------------------------------------- comm = MPI_COMM_WORLD CALL MPI_Init(info) IF (info .NE. 0) THEN WRITE(*,*) 'FAILED TO INITIALIZE MPI. ABORTING!' GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Get rank of processor !------------------------------------------------------------------------- CALL MPI_Comm_Rank(comm,rank,info) !------------------------------------------------------------------------- ! Get total number of processors !------------------------------------------------------------------------- CALL MPI_Comm_Size(comm,Nproc,info) !------------------------------------------------------------------------- ! Get name of processor !------------------------------------------------------------------------- CALL MPI_Get_Processor_Name(cbuf,ilen,info) #else !------------------------------------------------------------------------- ! Initialise number of processors !------------------------------------------------------------------------- Nproc = 1 rank = 0 #endif !------------------------------------------------------------------------- ! Define defaults !------------------------------------------------------------------------- CALL pse_defaults(info) IF (info .NE. 0) GOTO 9999 !------------------------------------------------------------------------- ! Get input arguments !------------------------------------------------------------------------- CALL pse_arguments(info) IF (info.NE.0) GOTO 9999 !------------------------------------------------------------------------- ! Read control file !------------------------------------------------------------------------- CALL pse_read_ctrl(ctrlfile,info) IF (info.NE.0) GOTO 9999 !------------------------------------------------------------------------- ! Check parameter validity !------------------------------------------------------------------------- CALL pse_check_ctrl(info) IF (info .NE. 0) GOTO 9999 #ifdef __MPI !------------------------------------------------------------------------- ! Define MPI precision !------------------------------------------------------------------------- IF (MK .EQ. KIND(1.0E0)) THEN MPPREC = MPI_REAL ELSEIF (MK .EQ. KIND(1.0D0)) THEN MPPREC = MPI_DOUBLE_PRECISION ELSE info = -1 GOTO 9999 ENDIF #endif !------------------------------------------------------------------------- ! Open log file !------------------------------------------------------------------------- IF (Nproc .LT. 10) THEN WRITE(cbuf,'(2A,I1.1)') TRIM(logfile),'.',rank ELSEIF (Nproc .LT. 100) THEN WRITE(cbuf,'(2A,I2.2)') TRIM(logfile),'.',rank ELSEIF (Nproc .LT. 1000) THEN WRITE(cbuf,'(2A,I3.3)') TRIM(logfile),'.',rank ELSE WRITE(cbuf,'(2A,I5.5)') TRIM(logfile),'.',rank ENDIF OPEN(pse_log_unit,FILE=cbuf,ACTION='WRITE',IOSTAT=ios) CLOSE(pse_log_unit,STATUS='DELETE') OPEN(pse_log_unit,FILE=cbuf,ACTION='WRITE',IOSTAT=ios) !------------------------------------------------------------------------- ! Compute global constants !------------------------------------------------------------------------- CALL pse_constants(info) IF (info.NE.0) GOTO 9999 !------------------------------------------------------------------------- ! Initialise the ppm library !------------------------------------------------------------------------- tolexp = INT(LOG10(EPSILON(cutoff)))+10 CALL ppm_init(ndim,MK,tolexp,comm,debug,info,ppm_log_unit) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_init','Failed to initialize ppm library.',info) GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Set ppm log file unit !------------------------------------------------------------------------- CALL ppm_io_set_unit(6,0,ppm_log_unit,info) IF (info .NE. 0) GOTO 9999 !------------------------------------------------------------------------- ! Write some of the parameters !------------------------------------------------------------------------- CALL pse_dump_parameters(info) !------------------------------------------------------------------------- ! Compute correction factor !------------------------------------------------------------------------- bsize(1) = dhx bsize(2) = dhy bsize(3) = dhz CALL ppm_comp_pp_correct(bsize,cutoff,ppm_kernel,eps,kpar,2.0_MK, & & kappa,info) IF (info.NE.0) GOTO 9999 kpar(2) = kpar(2)*kappa !------------------------------------------------------------------------- ! Compute table look up tables !------------------------------------------------------------------------- IF (lookup) THEN CALL ppm_comp_pp_mk_table(ppm_kernel,kpar,ntable,cutoff2,ktab, & & drtableinv,info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_init', & & 'Failed to create kernel lookup table.',info) GOTO 9999 ENDIF ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_init',t0,info) RETURN END SUBROUTINE pse_init ppm_pse/._pse_rhs.f0000644000175000017500000000057311451413033013226 0ustar omaromarMac OS X  2I{ATTRû"í{¨Ó¨Ó%com.apple.metadata:kMDItemWhereFromsbplist00¢_Qhttps://ppm.inf.ethz.ch/trac/export/641/ppmclients/ppm_pse/branches/1.0/pse_rhs.f_Nhttps://ppm.inf.ethz.ch/trac/browser/ppmclients/ppm_pse/branches/1.0/pse_rhs.f _°ppm_pse/ppm_pse.f0000644000175000017500000006052711455632307013031 0ustar omaromar !------------------------------------------------------------------------- ! Program : ppm_pse !------------------------------------------------------------------------- ! ! Purpose : 2d and 3d parallel simulations of diffusion using ! the PSE method. ! ! Note : This is a simple test client for PPM 1.2 !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! ETH Zurich ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- PROGRAM ppm_pse !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE ppm_module_map USE ppm_module_topo USE ppm_module_loadbal USE ppm_module_decomp USE ppm_module_neighlist USE ppm_module_ode USE ppm_module_user_util USE ppm_module_data USE pse_global USE pse_module_io USE pse_module_comp IMPLICIT NONE !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- #ifdef __MPI INCLUDE 'mpif.h' #endif !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER :: isymm,info,i,j,k INTEGER :: nstart,istep,topo_id,mesh_id INTEGER :: ntime INTEGER, DIMENSION(6) :: bcdef REAL(8) :: t0,t1,t2,tsum,t2sum,tdiff,tdiffmax REAL(8) :: t12,t22,tsum2,t2sum2,tdiff2 REAL(8) :: tminmax,tmaxmax REAL(MK) :: stime REAL(MK), DIMENSION(3 ) :: bsize INTEGER :: decomp,assig CHARACTER(LEN=MAXCHAR) :: cbuf INTEGER :: ode_wp_id,ode_wp_scheme INTEGER :: ode_wp_nstage,ode_wp_bfrsize INTEGER :: istage LOGICAL :: ode_adaptive REAL(MK), DIMENSION(4) :: ode_wp_time INTEGER , DIMENSION(1,1) :: ode_wp_ipack REAL(MK), DIMENSION(1,1) :: ode_wp_rpack REAL(MK), DIMENSION(:,:), POINTER :: ode_wp_bfr INTEGER :: n, nsubs REAL(MK), DIMENSION(:,:), POINTER :: minsub,maxsub INTEGER, DIMENSION(:), POINTER :: sub2proc REAL(MK), DIMENSION(3) :: hx !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- info = 0 CALL pse_init(info) IF (info .LT. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to initialize.',info) GOTO 8000 ENDIF IF (info .GT. 0) GOTO 9999 CALL substart('ppm_pse',t0,info) !------------------------------------------------------------------------- ! Define the boundary conditions !------------------------------------------------------------------------- bcdef(1:6) = ppm_param_bcdef_freespace !------------------------------------------------------------------------- ! Nullify data pointers !------------------------------------------------------------------------- NULLIFY(xp) NULLIFY(wp) NULLIFY(dwp) !--------------------------------------------------------------------- ! Not a restart: begin from step 0 !--------------------------------------------------------------------- nstart = 0 stime = 0.0_MK !--------------------------------------------------------------------- ! Check if an input file was specified !--------------------------------------------------------------------- IF (LEN_TRIM(inputfile) .GT. 0) THEN !----------------------------------------------------------------- ! Create a temporary topology for processing the geometry !----------------------------------------------------------------- topo_id = 0 mesh_id = -1 assig = ppm_param_assign_internal decomp = ppm_param_decomp_cuboid CALL ppm_mktopo(topo_id,mesh_id,xp,0,decomp,assig,min_phys, & & max_phys,bcdef,(/1,1,1/),sub_cost,istart,ndata,Ngrid,info) IF (info .NE. 0) THEN CALL pwrite(rank,'isd_init', & & 'Failed to create temporary topology.',info) info = -1 GOTO 9999 ENDIF CALL ppm_topo_get(topo_id,topo,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to get topology.',info) GOTO 8000 ENDIF !----------------------------------------------------------------- ! If we have a geometry input file, read it. !----------------------------------------------------------------- CALL pse_read_geometry(inputfile,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to read geometry.',info) GOTO 8000 ENDIF !----------------------------------------------------------------- ! Initialize the particles inside the geometry !----------------------------------------------------------------- CALL pse_create_particles(info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to create particles.',info) GOTO 8000 ENDIF ELSE info = 1 CALL pwrite(rank,'ppm_pse','No geometry input file defined',info) GOTO 8000 ENDIF !------------------------------------------------------------------------- ! End of init !------------------------------------------------------------------------- CALL pwrite(rank,'ppm_pse','Initialization complete.',info) IF(rank .EQ. 0) WRITE(pse_log_unit,'(A)') 'Initialization complete.' 1000 CONTINUE !------------------------------------------------------------------------- ! Timing !------------------------------------------------------------------------- #ifdef __MPI CALL MPI_Barrier(comm,info) #endif CALL pse_time2(t1) !------------------------------------------------------------------------- ! Perform domain decomposition !------------------------------------------------------------------------- IF (cutoff .GT. 0.0_MK) THEN !--------------------------------------------------------------------- ! Subdomain decomposition !--------------------------------------------------------------------- assig = ppm_param_assign_internal ppm_decomp = ppm_param_decomp_cuboid IF (ASSOCIATED(pcost)) THEN CALL ppm_mktopo(topo_id,xp,Npart,ppm_decomp,assig,min_phys,& & max_phys,bcdef,cutoff,sub_cost,info,pcost) ELSE CALL ppm_mktopo(topo_id,xp,Npart,ppm_decomp,assig,min_phys,& & max_phys,bcdef,cutoff,sub_cost,info) ENDIF ELSE topo_id = ppm_param_topo_undefined !------------------------------------------------------------------------ ! The following lines are commented because starting from version 1.2 ! of PPM ring topologies are not created anymore, particles can be ! treated in a ring topology implicitely !------------------------------------------------------------------------ ENDIF IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to set up topology.',info) GOTO 8000 ENDIF !------------------------------------------------------------------------- ! Timing !------------------------------------------------------------------------- CALL pse_time2(t2) WRITE(cbuf,'(A,F12.6)') 'TIMING: topology setup took: ',(t2-t1) CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) !------------------------------------------------------------------------- ! Free pcost: no longer needed !------------------------------------------------------------------------- IF (ASSOCIATED(pcost)) DEALLOCATE(pcost) !------------------------------------------------------------------------- ! Map the particles onto the topology !------------------------------------------------------------------------- CALL ppm_map_part_global(topo_id,xp,Npart,info) ! positions CALL ppm_map_part_push(wp,NSPEC,Npart,info) ! strengths CALL ppm_map_part_send(Npart,Mpart,info) ! send CALL ppm_map_part_pop(wp,NSPEC,Npart,Mpart,info) ! strengths CALL ppm_map_part_pop(xp,ndim,Npart,Mpart,info) ! positions IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to perform global mapping.',info) GOTO 8000 ENDIF Npart = Mpart !---------------------------------------------------------------------- ! Get ghost particle positions. This is done outside of the loop ! since they will never move. Has to be done BEFORE building ! neighbor lists. !---------------------------------------------------------------------- isymm = 0 IF (symmetry) isymm = 1 IF (cutoff .GT. 0.0_MK) THEN CALL ppm_map_part_ghost_get(topo_id,xp,ndim,Npart,isymm,cutoff,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to get ghosts.',info) GOTO 8000 ENDIF CALL ppm_map_part_send(Npart,Mpart,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to send ghosts.',info) GOTO 8000 ENDIF CALL ppm_map_part_pop(xp,ndim,Npart,Mpart,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to pop ghost positions.',info) GOTO 8000 ENDIF ENDIF !------------------------------------------------------------------------- ! Build neighbor lists once as particles do not move !------------------------------------------------------------------------- IF (cutoff .GT. 0.0_MK) THEN IF (neighlist .EQ. 1) THEN !------------------------------------------------------------------ ! Generate Verlet lists !------------------------------------------------------------------ CALL ppm_neighlist_vlist(topo_id,xp,Mpart,cutoff,skin,symmetry,& & vlist,nvlist,info) ELSEIF (neighlist .EQ. 2) THEN !----------------------------------------------------------------- ! Boxes need to be cutoff in all directions !----------------------------------------------------------------- bsize(1:3) = cutoff !----------------------------------------------------------------- ! Generate cell lists !----------------------------------------------------------------- CALL ppm_neighlist_clist(topo_id,xp,Mpart,bsize,symmetry,clist,& & Nm,info) ENDIF IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to build neighbor lists.',info) GOTO 8000 ENDIF ENDIF #ifdef __MPI !------------------------------------------------------------------------- ! Rebalance !------------------------------------------------------------------------- IF (rebalance) THEN IF (Nproc .LT. 2) THEN CALL pwrite(rank,'ppm_pse', & & 'Rebalancing makes no sense on 1 proc. Skipping.',info) rebalance = .FALSE. ELSE CALL pwrite(rank,'ppm_pse','Rebalancing computational load...', & & info) IF(rank .EQ. 0) WRITE(pse_log_unit,'(A)') '**** REBALANCE STARTS' !----------------------------------------------------------------- ! Allocate particle cost !----------------------------------------------------------------- ALLOCATE(pcost(Npart),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to allocate pcost.',info) GOTO 8000 ENDIF !----------------------------------------------------------------- ! Determine true particle costs !----------------------------------------------------------------- IF (cutoff .GT. 0.0_MK .AND. neighlist .LT. 3) THEN IF (neighlist .NE. 1) THEN !----------------------------------------------------------- ! Count lengths of Verlet lists !----------------------------------------------------------- CALL ppm_neighlist_vlist(topo_id,xp,Mpart,cutoff,skin,& & symmetry,vlist,nvlist,info,lstore=.FALSE.) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse', & & 'Failed to determine particle costs',info) GOTO 8000 ENDIF ENDIF DO i=1,Npart pcost(i) = REAL(nvlist(i),MK) ENDDO ELSE !-------------------------------------------------------------- ! All particles are equally expensive in the ring topology !-------------------------------------------------------------- DO i=1,Npart pcost(i) = 1.0_MK ENDDO ENDIF !----------------------------------------------------------------- ! Redo domain decomposition and mapping !----------------------------------------------------------------- rebalance = .FALSE. GOTO 1000 ENDIF ! Nproc .GT. 1 ENDIF #endif !------------------------------------------------------------------------- ! Allocate dwp !------------------------------------------------------------------------- ALLOCATE(dwp(NSPEC,Mpart),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to allocate dwp.',info) GOTO 8000 ENDIF !------------------------------------------------------------------------- ! Write initial diagnostics !------------------------------------------------------------------------- IF (.NOT. restart) THEN CALL pse_write_diag(stime,20,.TRUE.,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse', & & 'Writing intitial diagnostics failed.',info) GOTO 8000 ENDIF ENDIF !------------------------------------------------------------------------- ! Write initial condition to output file !------------------------------------------------------------------------- IF (.NOT. restart) THEN CALL pse_write_output(stime,nstart,20,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Writing intitial condition failed.',& & info) GOTO 8000 ENDIF ENDIF !------------------------------------------------------------------------- ! Init ODE module !------------------------------------------------------------------------- CALL ppm_ode_init(topo_id,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to init ODE module',info) GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Create ODE mode !------------------------------------------------------------------------- ode_wp_id = -1 ! auto assign ode_adaptive = .FALSE. ode_wp_scheme = ppm_param_ode_scheme_eulerf CALL ppm_ode_create_ode(ode_wp_id,ode_wp_bfrsize,ode_wp_nstage, & & ode_wp_scheme,ode_wp_scheme,ode_adaptive,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to create ODE mode',info) GOTO 9999 ENDIF ode_wp_nstage = 1 ode_wp_rpack(1,1) = 0.0_MK ode_wp_ipack(1,1) = ode_wp_nstage !------------------------------------------------------------------------- ! Allocate stage buffer !------------------------------------------------------------------------- ALLOCATE(ode_wp_bfr(ode_wp_bfrsize*NSPEC,Npart),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Failed to allocate ODE buffer',info) GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Set time information !------------------------------------------------------------------------- ntime = 0 tsum = 0.0 t2sum = 0.0 tsum2 = 0.0 t2sum2 = 0.0 tmaxmax = 0.0 tminmax = HUGE(tminmax) ode_wp_time(1) = stime ! start ode_wp_time(2) = nstep*dt ! end ode_wp_time(3) = stime ! current ode_wp_time(4) = dt ! step CALL ppm_ode_start(info) !------------------------------------------------------------------------- ! Do time stepping !------------------------------------------------------------------------- istep = nstart DO WHILE (.NOT. ppm_ode_alldone(info)) !---------------------------------------------------------------------- ! Timing !---------------------------------------------------------------------- #ifdef __MPI CALL MPI_Barrier(comm,info) #endif CALL pse_time2(t1) CALL pse_time(t12) !---------------------------------------------------------------------- ! Loop over ODE stages !---------------------------------------------------------------------- DO istage=1,ode_wp_nstage !------------------------------------------------------------------ ! PP interactions !------------------------------------------------------------------ CALL ppm_ode_step(ode_wp_id,xp,wp,dwp,NSPEC,Npart,ode_wp_bfr, & & istage,ode_wp_time,pse_rhs,ipackdata=ode_wp_ipack, & & rpackdata=ode_wp_rpack,info=info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','ODE step failed.',info) GOTO 8000 ENDIF ENDDO ! ODE stage !---------------------------------------------------------------------- ! Sync with ODE solver !---------------------------------------------------------------------- istep = istep + 1 dt = ode_wp_time(3) - stime stime = ode_wp_time(3) !---------------------------------------------------------------------- ! Timing !---------------------------------------------------------------------- CALL pse_time2(t2) CALL pse_time(t22) tdiff = t2-t1 tsum = tsum + tdiff t2sum = t2sum + (tdiff*tdiff) tdiff2 = t22-t12 ntime = ntime + 1 tdiffmax = tdiff2 #ifdef __MPI CALL MPI_Allreduce(tdiff2,tdiffmax,1,MPI_DOUBLE_PRECISION, & & MPI_MAX,comm,info) #endif __MPI tsum2 = tsum2 + tdiffmax t2sum2 = t2sum2 + (tdiffmax*tdiffmax) IF (tdiffmax .GT. tmaxmax) tmaxmax = tdiffmax IF (tdiffmax .LT. tminmax) tminmax = tdiffmax !---------------------------------------------------------------------- ! Output time step !---------------------------------------------------------------------- IF (rank .EQ. 0) THEN WRITE(cbuf,'(A,I6,A,E18.4)') 'Time step completed: ',istep, & & ' dt=',dt CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) WRITE(cbuf,'(A,E18.4)') 'Current simulation time is: ',stime CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ENDIF !---------------------------------------------------------------------- ! Write diagnostics file if needed !---------------------------------------------------------------------- IF (MOD(istep,freqdiag) .EQ. 0) THEN CALL pse_write_diag(stime,20,.FALSE.,info) IF (info .NE. 0) THEN CALL pwrite(rank,'ppm_pse','Writing diagnostics failed.',info) GOTO 8000 ENDIF ENDIF ENDDO ! time loop !------------------------------------------------------------------------- ! Print timing results !------------------------------------------------------------------------- IF (rank .EQ. 0) THEN ! Average time of MAX times WRITE(cbuf,'(A,F12.6)') 'TIMING: Avg(Max) seconds per time step: ', & & tsum2/REAL(ntime,8) CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ! Maximum of MAX times WRITE(cbuf,'(A,F12.6)') 'TIMING: Max(Max) seconds per time step: ', & & tmaxmax CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ! Minimum of MAX times WRITE(cbuf,'(A,F12.6)') 'TIMING: Min(Max) seconds per time step: ', & & tminmax CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ! Variance of MAX times WRITE(cbuf,'(A,F12.6)') 'TIMING: Var(Max) seconds per time step: ', & & (t2sum2-((tsum2*tsum2)/REAL(ntime,8)))/REAL((ntime-1),8) CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ENDIF ! Average time per time step on this processor WRITE(cbuf,'(A,F12.6)') 'TIMING: local Avg seconds per time step: ', & & tsum/REAL(ntime,8) CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ! Variance of time over time steps on this processor WRITE(cbuf,'(A,F12.6)') 'TIMING: local Var seconds per time step: ', & & (t2sum-((tsum*tsum)/REAL(ntime,8)))/REAL((ntime-1),8) CALL pwrite(rank,'ppm_pse',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) !------------------------------------------------------------------------- ! Finalize !------------------------------------------------------------------------- IF (rank.EQ.0) THEN CALL pwrite(rank,'ppm_pse','Done. Goodbye...',info) ENDIF CALL pse_finalize(0,info) IF (ASSOCIATED(dwp)) DEALLOCATE(dwp,STAT=info) !------------------------------------------------------------------------- ! Skip error handler !------------------------------------------------------------------------- GOTO 9999 !------------------------------------------------------------------------- ! Error handler !------------------------------------------------------------------------- 8000 CONTINUE CALL pse_finalize(1,info) IF (ALLOCATED(diff)) DEALLOCATE(diff,STAT=info) IF (ASSOCIATED(ktab)) DEALLOCATE(ktab,STAT=info) IF (ALLOCATED(kpar)) DEALLOCATE(kpar,STAT=info) IF (ASSOCIATED(dwp)) DEALLOCATE(dwp,STAT=info) IF (ASSOCIATED(wp)) DEALLOCATE(wp,STAT=info) IF (ASSOCIATED(xp)) DEALLOCATE(xp,STAT=info) IF (ASSOCIATED(sub_cost)) DEALLOCATE(sub_cost,STAT=info) IF (ASSOCIATED(pcost)) DEALLOCATE(pcost,STAT=info) !------------------------------------------------------------------------- ! End !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('ppm_pse',t0,info) END PROGRAM ppm_pse ppm_pse/pse_constants.f0000644000175000017500000001403111451402500014220 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_constants !------------------------------------------------------------------------- ! ! Purpose : Compute values of some global constants and set the ! ppm kernel parameters. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_constants(info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 REAL(MK) :: eps2,eps2inv,eps4inv,eps5inv REAL(MK) :: dh2,dh3 !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_constants',t0,info) !------------------------------------------------------------------------- ! Mesh spacings !------------------------------------------------------------------------- dhx = (max_phys(1)-min_phys(1))/REAL(Ngrid(1)-1,MK) dhy = (max_phys(2)-min_phys(2))/REAL(Ngrid(2)-1,MK) dhz = (max_phys(3)-min_phys(3))/REAL(Ngrid(3)-1,MK) pvol = dhx*dhy*dhz pvolinv = 1.0_MK/pvol !------------------------------------------------------------------------- ! Extended domain with ghost layer (for triangle sorting) !------------------------------------------------------------------------- min_physg(1) = min_phys(1)-cutoff min_physg(2) = min_phys(2)-cutoff min_physg(3) = min_phys(3)-cutoff max_physg(1) = max_phys(1)+cutoff max_physg(2) = max_phys(2)+cutoff max_physg(3) = max_phys(3)+cutoff !------------------------------------------------------------------------- ! Compute PI !------------------------------------------------------------------------- IF (MK .EQ. KIND(1.0D0)) THEN pi = ACOS(-1.0_MK) ELSE pi = ACOS(-1.0_MK) ENDIF pi2 = pi*pi !------------------------------------------------------------------------- ! Kernel core size !------------------------------------------------------------------------- eps2 = eps*eps eps2inv = 1.0_MK/eps2 eps4inv = eps2inv*eps2inv eps5inv = eps4inv/eps !------------------------------------------------------------------------- ! Intern-particle distance !------------------------------------------------------------------------- dh2 = dhx*dhy dh3 = dh2*dhz !------------------------------------------------------------------------- ! Cutoff !------------------------------------------------------------------------- cutoff2 = cutoff*cutoff !------------------------------------------------------------------------- ! Kernel parameters !------------------------------------------------------------------------- ALLOCATE(kpar(2),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_constants','Failed to allocate kpar.',info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Set ppm decomposition types !------------------------------------------------------------------------- ppm_decomp = ppm_param_decomp_tree !------------------------------------------------------------------------- ! Set kernel parameters !------------------------------------------------------------------------- IF (ndim .EQ. 2) THEN IF (kernel .EQ. 1) THEN ppm_kernel = ppm_param_kernel_laplace2d_p2 kpar(1) = eps2inv kpar(2) = (15.0_MK*dh2*eps4inv)/pi2 ELSE CALL pwrite(rank,'pse_constants','Invalid kernel',info) info = -1 GOTO 9999 ENDIF ELSE IF (kernel .EQ. 1) THEN ppm_kernel = ppm_param_kernel_laplace3d_p2 kpar(1) = eps2inv kpar(2) = (15.0_MK*dh3*eps5inv)/pi2 ELSE CALL pwrite(rank,'pse_constants','Invalid kernel',info) info = -1 GOTO 9999 ENDIF ENDIF !------------------------------------------------------------------------- ! PSE needs no skin for the Verlet lists as P do not move !------------------------------------------------------------------------- skin = 0.0_MK !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_constants',t0,info) RETURN END SUBROUTINE pse_constants ppm_pse/pse_time2.f0000644000175000017500000000577511445416206013256 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_time !------------------------------------------------------------------------- ! ! Purpose : Retreives the current time (cpu). Uses either ! MPI_Wtime, etime or the f90 CPU_TIME intrinsic, ! based on the DEFINES that are set (see ! ppm_define.h): ! __MPI uses MPI_Wtime ! __ETIME uses etime ! none uses CPU_TIME ! ! Input : ! ! Input/output : ! ! Output : timing (F) current CPU clock time ! ! Routines : etime (C intrinsic) ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_time2(timing) !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- IMPLICIT NONE #ifdef __MPI INCLUDE 'mpif.h' #endif !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- REAL(8), INTENT( OUT) :: timing !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- #ifdef __ETIME REAL(8) :: array(2) !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- REAL(8), EXTERNAL :: etime #endif #ifdef __MPI !------------------------------------------------------------------------- ! Call the MPI function MPI_Wtime !------------------------------------------------------------------------- timing = MPI_Wtime() #else !------------------------------------------------------------------------- ! Call the C routine: etime to get the cpu time !------------------------------------------------------------------------- #ifdef __ETIME timing = etime(array) #else CALL CPU_TIME(timing) #endif #endif !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE RETURN END SUBROUTINE pse_time2 ppm_pse/pse_module_comp.f0000644000175000017500000000364311445416205014530 0ustar omaromar !------------------------------------------------------------------------- ! Module : pse_module_comp !------------------------------------------------------------------------- ! ! Purpose : This module contains all subroutines needed for the ! actual Particle-Particle computations. ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_module_comp.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.3 2004/07/27 07:57:45 ivos ! Updated to new ppm moduel structure and use of the ppm PP ! computation routines instead of internal ones. ! ! Revision 1.2 2004/06/03 16:11:48 ivos ! Added pse_comp_pp_cell. ! ! Revision 1.1 2004/06/02 15:28:33 ivos ! Initial implementation. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- MODULE pse_module_comp !---------------------------------------------------------------------- ! Define interface to pse_rhs !---------------------------------------------------------------------- INTERFACE pse_rhs MODULE PROCEDURE pse_rhs END INTERFACE !---------------------------------------------------------------------- ! Include the source !---------------------------------------------------------------------- CONTAINS #include "pse_rhs.f" END MODULE pse_module_comp ppm_pse/pse_get_indicator.f0000644000175000017500000001060411445416205015033 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_get_indicator !------------------------------------------------------------------------- ! ! Purpose : Determine which grid points are inside and which ! ones are outside of the closed surface and ! initialize the level function to a binary indicator ! function, -1 inside, 1 outside. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status. 0 on success. ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_get_indicator.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:26 ivos ! CBL version of the PPM library ! ! Revision 1.1 2005/04/27 01:19:09 ivos ! Initial commit. Tested. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_get_indicator(info) !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER , INTENT(OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 INTEGER :: i,j,k,isub,jsub,ff REAL(MK) :: x,y,z REAL(MK), DIMENSION(3) :: ptx !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_get_indicator',t0,info) !CALL pwrite(rank,'pse_get_indicator','Determining indicator...',info) !------------------------------------------------------------------------- ! Loop over all local subs !------------------------------------------------------------------------- DO isub=1,topo%nsublist jsub = topo%isublist(isub) DO k=1,ndata(3,jsub) z = topo%min_subd(3,jsub) + REAL(k-1,MK)*dhz DO j=1,ndata(2,jsub) y = topo%min_subd(2,jsub) + REAL(j-1,MK)*dhy DO i=1,ndata(1,jsub) x = topo%min_subd(1,jsub) + REAL(i-1,MK)*dhx ptx(1) = x ptx(2) = y ptx(3) = z !--------------------------------------------------------- ! Determine if point is inside domain !--------------------------------------------------------- CALL pse_point_inside(ptx,ff,info) IF (info .NE. 0) GOTO 9999 flags(i,j,k,isub) = ff ENDDO ENDDO ENDDO ENDDO !------------------------------------------------------------------------- ! Free memory and deallocate surface triangulation !------------------------------------------------------------------------- IF(ASSOCIATED(triang)) DEALLOCATE(triang,STAT=info) CALL pse_alloc_binlist(2,0,0,info) !CALL pwrite(rank,'pse_get_indicator','Done with indicator.',info) !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_get_indicator',t0,info) RETURN END SUBROUTINE pse_get_indicator ppm_pse/._Ctrl0000644000175000017500000000036011451406106012240 0ustar omaromarMac OS X  2¾ðATTRö ›ð˜X˜Xcom.apple.quarantineq/0000;4c1b6d69;Google\x20Chrome;8DB9488C-D362-4124-8532-7D538D9E983D|com.google.Chromeppm_pse/pse_write_diag.f0000644000175000017500000003173111445416206014343 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_write_diag !------------------------------------------------------------------------- ! ! Purpose : This routine writes diagnostics to a file. ! For the parallel version, all processors send ! their information to rank 0, who writes the file. ! ! Input : ts (F) simulation time ! iUnit (I) IO unit to be used for file output ! lfirst (L) .TRUE. if this is the first ! CALL to this routine in the current ! run. If this is not a restart run, ! it will cause old diag files to be ! overwritten. Otherwise, data is ! appended. ! ! The global variables xp,wp are read. ! ! Input/output : ! ! Output : info (I) return status (0 if no error) ! ! Routines : ppm_io_open ! ppm_io ! ppm_io_close ! substart ! substop ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_write_diag.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.4 2004/07/27 07:57:46 ivos ! Updated to new ppm moduel structure and use of the ppm PP ! computation routines instead of internal ones. ! ! Revision 1.3 2004/06/21 13:20:48 ivos ! output and log entry are now precise. ! ! Revision 1.2 2004/06/04 08:46:33 ivos ! Added argument lfirst to allow old diag files to be overwritten if ! needed. ! ! Revision 1.1 2004/06/03 16:11:07 ivos ! Initial implementation. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_write_diag(ts,iUnit,lfirst,info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global USE ppm_module_io IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- REAL(MK) , INTENT(IN ) :: ts INTEGER , INTENT(IN ) :: iUnit LOGICAL , INTENT(IN ) :: lfirst INTEGER , INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 REAL(MK), DIMENSION(:), ALLOCATABLE :: mass_tot,mass_roi,fvec REAL(MK) :: mass CHARACTER(LEN=MAXCHAR) :: cbuf INTEGER :: i,j,ispec,ifmt,jUnit,ilen,posn INTEGER :: n !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_write_diag',t0,info) !------------------------------------------------------------------------- ! Check if diagnostics are enabled !------------------------------------------------------------------------- IF (LEN_TRIM(diagfile) .LT. 1) THEN CALL pwrite(rank,'pse_write_diag', & & 'Diagnostics are disabled. Skipping.',info) GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Format !------------------------------------------------------------------------- IF (diagfmt(1:1) .EQ. 'f' .OR. diagfmt(1:1) .EQ. 'F') THEN ifmt = ppm_param_io_ascii ELSE ifmt = ppm_param_io_binary ENDIF !------------------------------------------------------------------------- ! Allocate memory !------------------------------------------------------------------------- ALLOCATE(mass_tot(NSPEC),mass_roi(NSPEC),fvec(NSPEC+1),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_diag', & & 'Failed to allocate memory.',info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Compute diagnostics !------------------------------------------------------------------------- mass = 0.0_MK mass_tot = 0.0_MK mass_roi = 0.0_MK IF (ndim .EQ. 3) THEN DO ispec=1,NSPEC ! all species j = 0 n = 0 DO i=1,Npart ! all REAL particles, not the ghosts! mass = wp(ispec,i) n = n + 1 mass_tot(ispec) = mass_tot(ispec) + mass !------------------------------------------------------------- ! Box-shaped ROI !------------------------------------------------------------- IF ((xp(1,i).GE.min_roi(1)).AND.(xp(1,i).LE.max_roi(1)).AND.& & (xp(2,i).GE.min_roi(2)).AND.(xp(2,i).LE.max_roi(2)).AND.& & (xp(3,i).GE.min_roi(3)).AND.(xp(3,i).LE.max_roi(3))) THEN !------------------------------------------------------------- ! Spherical ROI of radius min_roi(1) !------------------------------------------------------------- ! IF (xp(1,i)*xp(1,i) + xp(2,i)*xp(2,i) + xp(3,i)*xp(3,i) & ! & .LT. min_roi(1)*min_roi(1)) THEN mass_roi(ispec) = mass_roi(ispec) + mass j = j + 1 ENDIF ENDDO IF (n.GT.0) mass_tot(ispec) = mass_tot(ispec)/REAL(n,MK) ! mass_roi(ispec) = (3.0_MK*mass_roi(ispec))/ & ! & (4.0_MK*pi*(min_roi(1)**3)) IF (debug .GT. 0) THEN WRITE(cbuf,'(A,I3,A,I6)') & & 'Particles in ROI for species ',ispec,': ',j CALL pwrite(rank,'pse_write_diag',cbuf,info) WRITE(cbuf,'(A,I3,A,F12.3)') & & 'Mass in ROI for species ',ispec,': ',mass_roi(ispec) CALL pwrite(rank,'pse_write_diag',cbuf,info) ENDIF ENDDO ELSE DO ispec=1,NSPEC ! all species j = 0 n = 0 DO i=1,Npart ! all REAL particles, not the ghosts! mass = wp(ispec,i) n = n + 1 mass_tot(ispec) = mass_tot(ispec) + mass IF ((xp(1,i).GE.min_roi(1)).AND.(xp(1,i).LE.max_roi(1)).AND.& & (xp(2,i).GE.min_roi(2)).AND.(xp(2,i).LE.max_roi(2))) THEN mass_roi(ispec) = mass_roi(ispec) + mass j = j + 1 ENDIF ENDDO IF (n.GT.0) mass_tot(ispec) = mass_tot(ispec)/REAL(n,MK) IF (debug .GT. 0) THEN WRITE(cbuf,'(A,I3,A,I6)') & & 'Particles in ROI for species ',ispec,': ',j CALL pwrite(rank,'pse_write_diag',cbuf,info) WRITE(cbuf,'(A,I3,A,F12.3)') & & 'Mass in ROI for species ',ispec,': ',mass_roi(ispec) CALL pwrite(rank,'pse_write_diag',cbuf,info) ENDIF ENDDO ENDIF !------------------------------------------------------------------------- ! Normalize with number of processors !------------------------------------------------------------------------- DO ispec=1,NSPEC mass_tot(ispec) = mass_tot(ispec) / REAL(Nproc,MK) mass_roi(ispec) = mass_roi(ispec) / REAL(roinp,MK) ENDDO !------------------------------------------------------------------------- ! Append to file. If this is the first output of a new run, replace. !------------------------------------------------------------------------- posn = ppm_param_io_append IF (lfirst .AND. .NOT. restart) posn = ppm_param_io_replace !------------------------------------------------------------------------- ! Open ppm I/O unit for centralized I/O !------------------------------------------------------------------------- jUnit = iUnit WRITE(cbuf,'(2A)') TRIM(diagfile),'.ROI' CALL ppm_io_open(jUnit,cbuf,ppm_param_io_write,posn,ifmt, & & ppm_param_io_centralized,info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_diag','Failed to open unit', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Write mass in ROI to file !------------------------------------------------------------------------- fvec(1) = ts/REAL(Nproc,MK) ! since we do a sum reduce in ppm_io fvec(2:NSPEC+1) = mass_roi(1:NSPEC) WRITE(cbuf,'(A,I6,A)') '(F12.6,',NSPEC,'(E18.6))' ilen = LEN_TRIM(cbuf) CALL ppm_io(jUnit,fvec,ppm_param_io_write,ppm_param_io_sum, & & IOFMT=cbuf(1:ilen),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_diag','Writing fvec failed', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Close file !------------------------------------------------------------------------- CALL ppm_io_close(jUnit,info) !------------------------------------------------------------------------- ! Open ppm I/O unit for centralized appending I/O !------------------------------------------------------------------------- jUnit = iUnit WRITE(cbuf,'(2A)') TRIM(diagfile),'.ALL' CALL ppm_io_open(jUnit,cbuf,ppm_param_io_write,posn,ifmt, & & ppm_param_io_centralized,info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_diag','Failed to open unit', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Write total mass to file !------------------------------------------------------------------------- fvec(1) = ts/REAL(Nproc,MK) ! since we do a sum reduce in ppm_io fvec(2:NSPEC+1) = mass_tot(1:NSPEC) WRITE(cbuf,'(A,I6,A)') '(F12.6,',NSPEC,'(E18.6))' ilen = LEN_TRIM(cbuf) CALL ppm_io(jUnit,fvec,ppm_param_io_write,ppm_param_io_sum, & & IOFMT=cbuf(1:ilen),STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_diag','Writing fvec failed', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Close file !------------------------------------------------------------------------- CALL ppm_io_close(jUnit,info) !------------------------------------------------------------------------- ! Deallocate memory !------------------------------------------------------------------------- DEALLOCATE(mass_tot,mass_roi,fvec,STAT=info) !------------------------------------------------------------------------- ! Confirm !------------------------------------------------------------------------- IF (rank .EQ. 0) THEN IF (lfirst .AND. .NOT. restart) THEN WRITE(cbuf,'(2A)') 'Diagnostics written to file ',TRIM(diagfile) ELSE WRITE(cbuf,'(2A)') 'Diagnostics appended to file ',TRIM(diagfile) ENDIF CALL pwrite(rank,'pse_write_diag',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_write_diag',t0,info) RETURN END SUBROUTINE pse_write_diag ppm_pse/pse_dump_parameters.f0000644000175000017500000004060111451403104015377 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_dump_parameters !------------------------------------------------------------------------- ! ! Purpose : Prints some of the parameters to stdout as well as ! the log file. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status. 0 on success. ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_dump_parameters(info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global USE ppm_module_user IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER, DIMENSION(8) :: tval CHARACTER(LEN = MAXCHAR) :: cbuf CHARACTER(LEN = 20) :: cfmt REAL(8) :: t0 INTEGER :: ilen !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_dump_parameters',t0,info) !------------------------------------------------------------------------- ! Root only !------------------------------------------------------------------------- IF (rank .EQ. 0) THEN !---------------------------------------------------------------------- ! Get time and date !---------------------------------------------------------------------- CALL DATE_AND_TIME(values=tval) WRITE(cbuf,'(A,I4.4,A,5(I2.2,A))') ' Date and Time of Start: ', & & tval(1),'-',tval(2),'-',tval(3),' ',tval(5),':',tval(6),':', & & tval(7),' ' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I4)') ' Number of processors: ',nproc ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Control and debug (command line parameters) !---------------------------------------------------------------------- WRITE(cbuf,'(2A)') ' Control file read: ',TRIM(ctrlfile) ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I1.1)') ' Debug level for this run: ',debug ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(2A)') ' Log file name: ',TRIM(logfile) ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Run case name !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- CASE DETAILS ---------------------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(2A)') ' This is problem case ',TRIM(casename) ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Restart or input file !---------------------------------------------------------------------- IF (restart) THEN WRITE(cbuf,'(2A)') ' Restarted from file: ',& & TRIM(rstfile) ELSE IF (LEN_TRIM(inputfile) .GT. 0) THEN WRITE(cbuf,'(2A)') ' Geometry file: ',& & TRIM(inputfile) ELSE WRITE(cbuf,'(A)') ' Particles created internally' ENDIF ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (ndim .GT. 2) THEN WRITE(cbuf,'(A,3I6)') ' Grid resolution: ', & & Ngrid(1:3) ELSE WRITE(cbuf,'(A,2I6)') ' Grid resolution: ', & & Ngrid(1:2) ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Problem definition !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- PROBLEM DEFINITION ---------------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (.NOT. restart) THEN WRITE(cbuf,'(A,I1.1)') ' Space dimensionality: ',ndim ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I3)') ' Number of species: ',NSPEC ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (ndim .EQ. 2) THEN WRITE(cbuf,'(A,4(F12.3))') ' Computational domain: ', & & min_phys(1:ndim),max_phys(1:ndim) ELSE WRITE(cbuf,'(A,6(F12.3))') ' Computational domain: ', & & min_phys(1:ndim),max_phys(1:ndim) ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) ENDIF WRITE(cfmt,'(A,I6,A)') '(A,',NSPEC,'F8.3)' WRITE(cbuf,cfmt) ' Diffusion constants: ',diff(1:NSPEC) ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Method parameters !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- SOLVER METHOD PARAMETERS ---------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (kernel .EQ. 1) THEN WRITE(cbuf,'(2A)') ' PSE kernel type: ','Poly_2nd' ELSEIF (kernel .EQ. 2) THEN WRITE(cbuf,'(2A)') ' PSE kernel type: ','Gauss_2nd' ELSEIF (kernel .EQ. 3) THEN WRITE(cbuf,'(2A)') ' PSE kernel type: ','Gauss_4th' ELSEIF (kernel .EQ. 4) THEN WRITE(cbuf,'(2A)') ' PSE kernel type: ','Gauss_6th' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (.NOT. restart) THEN WRITE(cbuf,'(A,E12.4)') ' PSE kernel radius: ',eps ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (ndim .GT. 2) THEN WRITE(cbuf,'(A,3E12.3)') ' Inter-particle distances: ', & & dhx,dhy,dhz ELSE WRITE(cbuf,'(A,2E12.3)') ' Inter-particle distances: ', & & dhx,dhy ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) ENDIF IF (cutoff .GT. 0.0_MK) THEN WRITE(cbuf,'(A,E12.4)') ' P-P interaction cutoff: ',cutoff ELSE WRITE(cbuf,'(2A)') ' P-P interaction cutoff: ','none' ENDIF IF (neighlist .EQ. 1) THEN WRITE(cbuf,'(2A)') ' Neighbor list type: ','Verlet' ELSEIF (neighlist .EQ. 2) THEN WRITE(cbuf,'(2A)') ' Neighbor list type: ','Cell' ELSEIF (neighlist .EQ. 3) THEN WRITE(cbuf,'(2A)') ' Neighbor list type: ','None' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (symmetry) THEN WRITE(cbuf,'(2A)') ' Symmetric P-P interactn.: ','yes' ELSE WRITE(cbuf,'(2A)') ' Symmetric P-P interactn.: ','no' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (lookup) THEN WRITE(cbuf,'(2A)') ' Use kernel lookup table : ','yes' ELSE WRITE(cbuf,'(2A)') ' Use kernel lookup table : ','no' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Time integration !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- TIME INTEGRATION -----------------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I8)') ' Number of time steps: ',nstep ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,E12.4)') ' Base time step size: ',dt ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I3)') ' Number of STS stages: ',NSTS ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,E12.4)') ' STS kick-in time: ',STSstart ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Domain decomposition !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- DOMAIN DECOMPOSITION -------------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(2A)') ' Domain decomposition used: octtree' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Diagnostics !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- DIAGNOSTIC OUTPUT ----------------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (LEN_TRIM(diagfile) .LT. 1) THEN WRITE(cbuf,'(2A)') ' Diagnostics file: ','disabled' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) ELSE IF (ndim .EQ. 2) THEN WRITE(cbuf,'(A,4(F12.3))') ' Region of interest: ', & & min_roi(1:ndim),max_roi(1:ndim) ELSE WRITE(cbuf,'(A,6(F12.3))') ' Region of interest: ', & & min_roi(1:ndim),max_roi(1:ndim) ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I6)') ' Time steps between diag: ',freqdiag ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(2A)') ' Diagnostics file: ',TRIM(diagfile) ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) ENDIF !---------------------------------------------------------------------- ! Output !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- RESULT DATA OUTPUT ---------------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I6)') ' Time steps between output:',freqoutput ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(2A)') ' Output file name stub: ',TRIM(outputfile) ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Restart !---------------------------------------------------------------------- WRITE(cbuf,'(A)') '---- PROGRAM CONTROL FLAGS ------------------------' ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) WRITE(cbuf,'(A,I6)') ' Time steps between rst: ',freqrestart ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) !---------------------------------------------------------------------- ! Flags !---------------------------------------------------------------------- IF (checkmap) THEN WRITE(cbuf,'(2A)') ' Check mapping results: ','yes' ELSE WRITE(cbuf,'(2A)') ' Check mapping results: ','no' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (checkunique) THEN WRITE(cbuf,'(2A)') ' Check particle uniqueness:','yes' ELSE WRITE(cbuf,'(2A)') ' Check particle uniqueness:','no' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (probeproc) THEN WRITE(cbuf,'(2A)') ' Probe processor speeds: ','yes' ELSE WRITE(cbuf,'(2A)') ' Probe processor speeds: ','no' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (dumpparticles) THEN WRITE(cbuf,'(2A)') ' Dump particle positions: ','yes' ELSE WRITE(cbuf,'(2A)') ' Dump particle positions: ','no' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) IF (rebalance) THEN WRITE(cbuf,'(2A)') ' Rebalance comput load: ','yes' ELSE WRITE(cbuf,'(2A)') ' Rebalance comput load: ','no' ENDIF ilen = LEN_TRIM(cbuf) WRITE(pse_log_unit,'(A)') cbuf(1:ilen) WRITE(*,'(A)') cbuf(1:ilen) ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_dump_parameters',t0,info) RETURN END SUBROUTINE pse_dump_parameters ppm_pse/Ctrl0000644000175000017500000002353211455632152012037 0ustar omaromar#-------------------------------------------------------------------------- # Sample control file for ppm_pse # # Edit the settings below # # IFS, 8.6.2006 #-------------------------------------------------------------------------- #-------------------------------------------------------------------------- # Case: # casename: name of this case. "ppm_pse" is the default. # Must be between 1 and 30 characters long! #-------------------------------------------------------------------------- casename = BOXTEST #-------------------------------------------------------------------------- # Input: # inputfile: name of geometry input file to be read. Valid file # formats are OpenInventor (.iv) ASCII files or # raw (.surf) files containing a closed and connected # triangulated surface. # Ngrid : Particles will be initialized using this number # of grid points (incl. boundary points) in each # direction in the bounding box of the geometry. If a # geometry is specified, particles will only be created # inside the geometry. # nBins : number of bins in y and z direction for the fast # sorting of the triangles of the geometry file. #-------------------------------------------------------------------------- inputfile = box.surf # 1 CPU, h = 0.03125 #Ngrid = 129,129,129 # 8 CPUs, h = 0.015625 #Ngrid = 257,257,257 # 64 CPUs, h = 0.0078125 #Ngrid = 513,513,513 # 512 CPUs, h = 0.00390625 #Ngrid = 1025,1025,1025 # 4096 CPUs, h = 0.001953125 #Ngrid = 2049,2049,2049 Ngrid = 43,43,43 #Ngrid = 65,65,65 nBins = 2,2 #-------------------------------------------------------------------------- # Output: # freqoutput: frequency for writing output files (in time steps). # These files contain the full concentration field. # outputfile: name stub of ouput files. .out will be # appended # outputfmt : format of output file. One of: # FORMATTED (i.e. ASCII text) # UNFORMATTED (i.e. binary) #-------------------------------------------------------------------------- freqoutput = 100000 outputfile = pse_output outputfmt = UNFORMATTED #-------------------------------------------------------------------------- # Restart: # restart : set to .TRUE. if this is a continuation run. In # this case, the input file will not be read. # restartfile: name of the restart file which will be read if # restart = .TRUE. # freqrestart: frequency for writing restart files (in time steps) #-------------------------------------------------------------------------- restart = FALSE restartfile = pse_restart00000011.rst freqrestart = 1000 #-------------------------------------------------------------------------- # Problem definition: # dimensions: space dimension. Either 2 or 3. # domain : computational domain. First 2 (or 3) numbers are # lower x,y,(z) bounds, second 2 (or 3) upper ones. # species : number of particle species # isodiff : isotropic diffusion constants for all species. Comma # separated list of length [species] #-------------------------------------------------------------------------- dimensions = 3 domain = 0.0, 0.0, 0.0, 4.0, 4.0, 4.0 species = 1 isodiff = 0.03 #-------------------------------------------------------------------------- # Method parameters: # kernel : PSE kernel type. One of: # 1 2nd order polynomial kernel of Cottet (deflt) # 2 2nd order Gaussian kernel # epsilon : PSE kernel core radius. Must not be smaller than 1/h. # 1/h is a good choice where h is the grid width, i.e. # (domain length)/(Ngrid-1). # cutoff : cutoff for particle-particle interactions. If cutoff # .LT. 0 is specified, no cutoff will be used and the # full N**2 problem is solved. This will be SLOW. # For kernel type 1, a cutoff of about 2*eps is good, # kernel 2 needs 3 to 5*eps. # symmetry: set this to .TRUE. to use symmetry in # particle-particle interactions, to .FALSE. otherwise #-------------------------------------------------------------------------- kernel = 1 # 1 CPU #epsilon = 0.034 #cutoff = 0.07 # 8 CPU #epsilon = 0.017 #cutoff = 0.035 # 64 CPU #epsilon = 0.008 #cutoff = 0.017 # 512 CPU #epsilon = 0.004 #cutoff = 0.009 # 512 CPU #epsilon = 0.004 #cutoff = 0.009 # 4096 CPU #epsilon = 0.002 #cutoff = 0.003 epsilon = 0.1 cutoff = 0.21 symmetry = TRUE #-------------------------------------------------------------------------- # Time stepping: # NSTS : Number of stages used in the STS scheme. Set this to # 1 to recover explicit Euler. The STS acceleration scheme # is described in # # V. Alexiades, G. Amiez, and P.-A. Gremaud. Super-time-stepping # acceleration of explicit schemes for parabolic problems. Comm. # Numer. Meth. Eng., 12(1):31-42, 1996. # # STSstart: time at which to switch from Euler to STS # nstep: total number of time steps to be performed in the # simulation. # dt : time step size for Euler. STS will automatically # adjust this. For stability, this has to be smaller than # h^2/(2*D) #-------------------------------------------------------------------------- nSTS = 1 STSstart = 3.0 nstep = 20 # 1 CPU #dt = 0.005 # 8 CPU #dt = 0.001 # 64 CPU #dt = 0.0003 # 512 CPU #dt = 0.0001 # 4096 CPU #dt = 0.00002 dt = 0.01 #-------------------------------------------------------------------------- # Diagnostics: # freqdiag: frequency of diagnostics output in time steps # diagnosticsfile: name of file to which diagnostics will be # APPENDED. # diagnosticsfmt: format of diagnostics file. One of: # FORMATTED (i.e. ASCII text) # UNFORMATTED (i.e. binary) # roi : extent of region of interest for diagnostics. # First 2 (or 3) numbers are lower x,y,(z) bounds, # second 2 (or 3) upper ones. The average value of all # particles in this region is the diagnostic. #-------------------------------------------------------------------------- freqdiag = 1 diagnosticsfile = pse_diag diagnosticsfmt = FORMATTED roi = 2.0, 2.0, 0.0, 3.0, 3.0, 4.0 #-------------------------------------------------------------------------- # Program control parameters: # debug: ppm debug level. 0,1, or 2. Set to 0 for normal run. # abortfile: name of the abort file. At every time step, the # program will check if a file with this name is # present in the directory where it was started. If # such a file is found, the simulation terminates # gracefully, writing restart and output files. # logfile: name of the log file used by ppm_pse # neighlist: Type of neighbor list to be used. One of: # VERLET (Verlet lists) # CELL (Cell list) # NONE (Ring topology. Cutoff ignored) # Verlet lists are the fastest, but need more memory. # If you run out of memory, try CELL lists. NONE is # for testing purposes only and is extemely SLOW. # domainDecomp: Type of domain decomposition to be used (if # cutoff is .GT. 0.0) in a parallel run. One of: # OctTree (default) # BinTree_prune # pruned_cell # bisection # xPencil # yPencil # zPencil # xySlab # xzSlab # yzSlab # cuboid # lookup: Use kernel lookup table? Say .FALSE. if unsure. # checkmap: Check to result of the global map? This only makes # sense in a parallel run. # checkunique: Check that particles are unique? Say .FALSE. if # unsure. # probeproc: Probe relative processor speeds and use this info # for load balancing? This only makes sense in # a parallel run. # dumpparticles: Dump particle positions in ASCII files? # WARNING: large files will be written!! # rebalance: Rebalance computational load according to true # computational costs of all particles. This # requires the topology setup to be done twice. And # considerably slows down the initialization of the # simulation. It is only worth it in parallel runs with # more than 32 processors. #-------------------------------------------------------------------------- debug = 0 abortfile = ABORT logfile = ppm_pse.log neighlist = CELL domainDecomp = octtree lookup = FALSE checkmap = FALSE checkunique = FALSE probeproc = FALSE dumpparticles = FALSE ppm_pse/pse_alloc_binlist.f0000644000175000017500000001450211451403574015041 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_alloc_binlist !------------------------------------------------------------------------- ! ! Purpose : This routine allocates memory for the triangle bin ! lists. If called with flag=0,initial ! allocation is done. flag=1 enlarges the existing ! arrays and flag=2 deallocates all arrays. ! For the parallel version,all processors send ! their information to rank 0,who writes the file. ! ! Input : flag (I) Action to perform. One of: ! 0 initial allocation ! 1 reallocation (grow) ! 2 deallocate ! iy (I) y index of bin to allocate list ! for. ! iz (I) z index of bin to allocate list ! for. ! ! Input/output : ! ! Output : info (I) return status (0 if no error) ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_alloc_binlist.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:26 ivos ! CBL version of the PPM library ! ! Revision 1.1 2005/04/27 01:19:07 ivos ! Initial commit. Tested. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum,Hirschengraben 84 ! CH-8092 Zurich,Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_alloc_binlist(flag,iy,iz,info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global USE pse_module_realloc IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER ,INTENT(IN ) :: flag,iy,iz INTEGER ,INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- ! enlarge arrays by STEP entries everytime this is called INTEGER,PARAMETER :: STEP = 10 INTEGER :: i,j INTEGER :: llen CHARACTER(LEN=256) :: cbuf REAL(8) :: t0 !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_alloc_binlist',t0,info) !------------------------------------------------------------------------- ! flag=0: Initialize !------------------------------------------------------------------------- IF (flag .EQ. 0) THEN !---------------------------------------------------------------------- ! Bin lists (contain all triangles inside bin) !---------------------------------------------------------------------- ! allocate matrix to hold lengths of lists ! first index: biny,second index: binz ALLOCATE(nbin(NBINY,NBINZ),STAT=info) if (info .NE. 0) THEN CALL pwrite(rank,'pse_alloc_binlist','Not enough memory',info) info = -1 GOTO 9999 ENDIF nbin(1:NBINY,1:NBINZ) = 0 ! allocate array of pointers to lists ! first index: biny,second index: binz ALLOCATE(bin(NBINY,NBINZ),STAT=info) if (info .NE. 0) THEN CALL pwrite(rank,'pse_alloc_binlist','Not enough memory',info) info = -1 GOTO 9999 ENDIF ! nullify all pointers (associate no lists at the beginning) DO i=1,NBINY DO j=1,NBINZ NULLIFY(bin(i,j)%list) ENDDO ENDDO !------------------------------------------------------------------------ ! flag=1: Enlarge lists !------------------------------------------------------------------------ ELSEIF (flag .EQ. 1) THEN IF (ASSOCIATED(bin(iy,iz)%list)) THEN ! if list is already assigned,reallocate it llen = SIZE(bin(iy,iz)%list) bin(iy,iz)%list=>reallocate(bin(iy,iz)%list,llen+STEP) ELSE ! if there is no list yet,allocate one ALLOCATE(bin(iy,iz)%list(STEP),STAT=info) IF (info .NE. 0) THEN WRITE(cbuf,'(A,I5,A,I5,A)') 'Error reallocating bin list(' & & ,iy,',',iz,')' CALL pwrite(rank,'pse_alloc_binlist',cbuf,info) info = -1 GOTO 9999 ENDIF ENDIF !------------------------------------------------------------------------ ! flag=2: Deallocate all lists !------------------------------------------------------------------------ ELSEIF (flag .EQ. 2) THEN DO i=1,NBINY DO j=1,NBINZ IF (ASSOCIATED(bin(i,j)%list)) DEALLOCATE(bin(i,j)%list) ENDDO ENDDO IF (ALLOCATED(bin)) DEALLOCATE(bin) IF (ALLOCATED(nbin)) DEALLOCATE(nbin) ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_alloc_binlist',t0,info) RETURN END SUBROUTINE pse_alloc_binlist ppm_pse/pse_write_output.f0000644000175000017500000002014311451404256014771 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_write_output !------------------------------------------------------------------------- ! ! Purpose : This routine writes an output file for ! postprocessing. ! For the parallel version, all processors send ! their information to rank 0, who writes the file. ! ! Input : ts (F) current simulation time ! istep (I) time step number completed ! iUnit (I) IO unit to be used for file output ! ! The global variables xp,wp are read. ! ! Input/output : ! ! Output : info (I) return status (0 if no error) ! ! Routines : ppm_io_open ! ppm_io ! ppm_io_close ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_write_output.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.5 2004/09/01 09:05:40 ivos ! bugifx: do not write out ghosts! ! ! Revision 1.4 2004/07/27 07:57:46 ivos ! Updated to new ppm moduel structure and use of the ppm PP ! computation routines instead of internal ones. ! ! Revision 1.3 2004/06/04 12:04:05 ivos ! Now writes eps instead of eps2. ! ! Revision 1.2 2004/06/02 15:45:19 ivos ! Added log outputs. ! ! Revision 1.1 2004/06/02 15:28:37 ivos ! Initial implementation. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_write_output(ts,istep,iUnit,info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global USE ppm_module_io IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- REAL(MK) , INTENT(IN ) :: ts INTEGER , INTENT(IN ) :: istep,iUnit INTEGER , INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 INTEGER , DIMENSION(4 ) :: ivec REAL(MK), DIMENSION(2 ) :: fvec INTEGER :: ifmt,jUnit,i CHARACTER(LEN=MAXCHAR) :: filename,cbuf !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_write_output',t0,info) WRITE(filename,'(A,I8.8,A)') TRIM(outputfile),istep,'.out' !------------------------------------------------------------------------- ! Format !------------------------------------------------------------------------- IF (outputfmt(1:1) .EQ. 'f' .OR. outputfmt(1:1) .EQ. 'F') THEN ifmt = ppm_param_io_ascii ELSE ifmt = ppm_param_io_binary ENDIF OPEN(43,FILE=filename,STATUS='REPLACE',ACTION='WRITE') DO i=1,Npart IF (((xp(2,i).GT.-0.5_MK*dhy).AND.(xp(2,i).LT.0.5_MK*dhy)).AND. & & ((xp(3,i).GT.-0.5_MK*dhz).AND.(xp(3,i).LT.0.5_MK*dhz))) THEN WRITE(43,'(4F12.6)') xp(1:3,i),wp(1,i) ENDIF ENDDO CLOSE(43) GOTO 9999 !------------------------------------------------------------------------- ! Open ppm I/O unit for centralized I/O !------------------------------------------------------------------------- jUnit = iUnit CALL ppm_io_open(jUnit,filename,ppm_param_io_write,ppm_param_io_replace,& & ifmt,ppm_param_io_centralized,info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_output','Failed to open unit', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Write Ntot to file !------------------------------------------------------------------------- CALL ppm_io(jUnit,Npart,ppm_param_io_write,ppm_param_io_sum,STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_output','Writing Ntot failed', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Write MK,istep,NSPEC,ndim to file !------------------------------------------------------------------------- ivec(1) = MK ivec(2) = istep ivec(3) = NSPEC ivec(4) = ndim CALL ppm_io(jUnit,ivec,ppm_param_io_write,ppm_param_io_root,STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_output','Writing ivec failed', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Write ts,eps to file !------------------------------------------------------------------------- fvec(1) = ts fvec(2) = eps CALL ppm_io(jUnit,fvec(1:2),ppm_param_io_write, & & ppm_param_io_root,STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_output','Writing fvec failed', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Write xp !------------------------------------------------------------------------- CALL ppm_io(jUnit,xp(1:ndim,1:Npart),ppm_param_io_write, & & ppm_param_io_concat,STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_output','Writing xp failed', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Write wp !------------------------------------------------------------------------- CALL ppm_io(jUnit,wp(1:NSPEC,1:Npart),ppm_param_io_write, & & ppm_param_io_concat,STAT=info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_write_output','Writing xp failed', & & info) info = -1 GOTO 9999 ENDIF !------------------------------------------------------------------------- ! Close file !------------------------------------------------------------------------- CALL ppm_io_close(jUnit,info) !------------------------------------------------------------------------- ! Confirm !------------------------------------------------------------------------- IF (rank .EQ. 0) THEN WRITE(cbuf,'(2A)') 'Output written to file ',TRIM(filename) CALL pwrite(rank,'pse_write_output',cbuf,info) WRITE(pse_log_unit,'(A)') TRIM(cbuf) ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_write_output',t0,info) RETURN END SUBROUTINE pse_write_output ppm_pse/pse_finalize.f0000644000175000017500000001507211445416205014025 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_finalize !------------------------------------------------------------------------- ! ! Purpose : Clean up global memory and close open units. This ! also finalizes ppm and mpi. ! ! Input : errstat (I) error status. 0 if this is a ! normal termination. .NE. 0 if ! this is a program abort. ! ! Input/output : ! ! Output : info (I) return status. 0 on success. ! ! Routines : substart ! substop ! ppm_finalize ! MPI_Finalize ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_finalize.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.6 2005/04/27 01:21:44 ivos ! Included geometry processing directly into the client. No more ! preprocessor is needed. This program now directly reads OpenInventor ! files. Also, changed the time stepper to ppm_module_ode. ! ! Revision 1.5 2004/07/30 13:04:12 ivos ! bugfix: ktab needs to be POINTER (is allocated within ppm_comp_ ! pp_mk_table) and tested for ASSOCIATED and not ALLOCATED. ! ! Revision 1.4 2004/07/27 07:57:44 ivos ! Updated to new ppm moduel structure and use of the ppm PP ! computation routines instead of internal ones. ! ! Revision 1.3 2004/06/09 09:16:57 ivos ! Added errstat argument. ! ! Revision 1.2 2004/06/03 07:53:28 walther ! Bug fix: added ppm_module_user. ! ! Revision 1.1 2004/06/02 15:28:30 ivos ! Initial implementation. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_finalize(errstat,info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global USE ppm_module_neighlist USE ppm_module_user_util IMPLICIT NONE !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- #ifdef __MPI INCLUDE 'mpif.h' #endif !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER, INTENT(IN ) :: errstat INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_finalize',t0,info) !------------------------------------------------------------------------- ! Deallocate global memory !------------------------------------------------------------------------- IF (ALLOCATED(diff)) DEALLOCATE(diff,STAT=info) IF (ASSOCIATED(ktab)) DEALLOCATE(ktab,STAT=info) IF (ALLOCATED(kpar)) DEALLOCATE(kpar,STAT=info) IF (ASSOCIATED(wp)) DEALLOCATE(wp,STAT=info) IF (ASSOCIATED(xp)) DEALLOCATE(xp,STAT=info) IF (ASSOCIATED(dwp)) DEALLOCATE(dwp,STAT=info) ! IF (ASSOCIATED(min_sub)) DEALLOCATE(min_sub,STAT=info) ! IF (ASSOCIATED(max_sub)) DEALLOCATE(max_sub,STAT=info) ! IF (ASSOCIATED(sub2proc)) DEALLOCATE(sub2proc,STAT=info) ! IF (ASSOCIATED(isublist)) DEALLOCATE(isublist,STAT=info) ! IF (ASSOCIATED(proc_speed)) DEALLOCATE(proc_speed,STAT=info) ! IF (ASSOCIATED(sub_cost)) DEALLOCATE(sub_cost,STAT=info) ! IF (ASSOCIATED(pcost)) DEALLOCATE(pcost,STAT=info) !------------------------------------------------------------------------- ! Deallocate neighbor lists !------------------------------------------------------------------------- CALL ppm_clist_destroy(clist,info) IF (ASSOCIATED(Nm)) DEALLOCATE(Nm,STAT=info) IF (ASSOCIATED(vlist)) DEALLOCATE(vlist,STAT=info) IF (ASSOCIATED(nvlist)) DEALLOCATE(nvlist,STAT=info) !------------------------------------------------------------------------- ! Close log file !------------------------------------------------------------------------- CLOSE(pse_log_unit,IOSTAT=info) !------------------------------------------------------------------------- ! Finalize PPM !------------------------------------------------------------------------- CALL ppm_finalize(info) #ifdef __MPI IF (errstat .EQ. 0) THEN !--------------------------------------------------------------------- ! Finalize MPI !--------------------------------------------------------------------- CALL MPI_Finalize(info) IF (info .NE. 0) THEN WRITE(*,*) 'FAILED TO FINALIZE MPI. BAILING OUT!' GOTO 9999 ENDIF ELSE !--------------------------------------------------------------------- ! Abort MPI !--------------------------------------------------------------------- CALL MPI_Abort(comm,errstat,info) IF (info .NE. 0) THEN WRITE(*,*) 'FAILED TO ABORT MPI!' GOTO 9999 ENDIF ENDIF #endif !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_finalize',t0,info) RETURN END SUBROUTINE pse_finalize ppm_pse/pse_check_ctrl.f0000644000175000017500000002502311445416204014321 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_check_ctrl !------------------------------------------------------------------------- ! ! Purpose : Check validity of parameter settings. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status. 0 on success. ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_check_ctrl(info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 INTEGER :: i !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_check_ctrl',t0,info) !------------------------------------------------------------------------- ! Check case name !------------------------------------------------------------------------- i = LEN_TRIM(casename) IF (i .LT. 1) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL NOTICE: Case name too short. Reset to ppm_pse.',info) casename = 'ppm_pse' ENDIF IF (i .GT. 30) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL NOTICE: Case name too long. Truncated to 30 chars.',info) casename = casename(1:20) ENDIF !------------------------------------------------------------------------- ! Check dimensionality !------------------------------------------------------------------------- IF (ndim .NE. 2 .AND. ndim .NE. 3) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: Space dimensionality must be either 2 or 3.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check number of stages !------------------------------------------------------------------------- IF (NSTS .LT. 1 .OR. NSTS .GT. 999998) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: NSTS must be between 1 and 999998.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check epsilon !------------------------------------------------------------------------- IF (eps .LE. 0.0_MK) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: epsilon must be > 0.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check neighbor list type !------------------------------------------------------------------------- IF (neighlist .EQ. 3 .AND. cutoff .GT. 0.0_MK) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL NOTICE: Cutoff will be ignored without neighbor lists',info) cutoff = -1.0_MK info = 0 ENDIF !------------------------------------------------------------------------- ! Check kernel type !------------------------------------------------------------------------- IF (kernel .EQ. 3 .OR. kernel .EQ. 4) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL NOTICE: Kernel is not yet implemented. Reset to 1.',info) kernel = 1 ENDIF IF (kernel .NE. 1 .AND. kernel .NE. 2) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: Invalid kernel specified. Reset to default.',info) info = 0 kernel = 1 ENDIF !------------------------------------------------------------------------- ! Check overlap !------------------------------------------------------------------------- IF (dhx/eps .GT. 1.0_MK) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL WARNING: Overlap h/eps is > 1 in x direction.',info) info = 0 ENDIF IF (dhy/eps .GT. 1.0_MK) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL WARNING: Overlap h/eps is > 1 in y direction.',info) info = 0 ENDIF IF (dhz/eps .GT. 1.0_MK) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL WARNING: Overlap h/eps is > 1 in z direction.',info) info = 0 ENDIF !------------------------------------------------------------------------- ! Check number of bins !------------------------------------------------------------------------- IF ((NBINY .LT. 1) .OR. (NBINZ .LT. 1)) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: NBINS must be > 0 in all directions.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check grid size !------------------------------------------------------------------------- IF (Ngrid(1) .LE. 0 .OR. Ngrid(2) .LE. 0 .OR. Ngrid(3) .LE. 0) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: Grid size must be > 0 in all directions.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check time step !------------------------------------------------------------------------- IF (dt .LE. 0.0_MK) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: time step dt must be > 0.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check time step number !------------------------------------------------------------------------- IF (nstep .LE. 0) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL WARNING: Performing 0 time steps.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check output frequency !------------------------------------------------------------------------- IF (freqoutput .LE. 0) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: Output frequency must be > 0.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check diagnostics frequency !------------------------------------------------------------------------- IF (nstep .LE. 0) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: Diagnostics frequency must be > 0.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check restart frequency !------------------------------------------------------------------------- IF (nstep .LE. 0) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: Restart frequency must be > 0.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check number of species !------------------------------------------------------------------------- IF (NSPEC .LE. 0) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: Number of species must be > 0.',info) info = -1 ENDIF !------------------------------------------------------------------------- ! Check diffusion constants !------------------------------------------------------------------------- DO i=1,NSPEC IF (diff(i) .LT. 0.0_MK) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: All diffusion constants must be >= 0.',info) info = -1 ENDIF ENDDO !------------------------------------------------------------------------- ! Check domain size !------------------------------------------------------------------------- DO i=1,ndim IF (min_phys(i) .GT. max_phys(i)) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: min_phys must be > max_phys in all dims',info) info = -1 ENDIF ENDDO !------------------------------------------------------------------------- ! Check roi !------------------------------------------------------------------------- DO i=1,ndim IF (min_roi(i) .GT. max_roi(i)) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: min_roi must be > max_roi in all dims',info) info = -1 ENDIF IF (min_roi(i) .LT. min_phys(i)) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: min_roi must be >= min_phys in all dims',info) info = -1 ENDIF IF (max_roi(i) .GT. max_phys(i)) THEN CALL pwrite(rank,'pse_check_ctrl', & & 'CTRL ERROR: max_roi must be <= max_phys in all dims',info) info = -1 ENDIF ENDDO !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_check_ctrl',t0,info) RETURN END SUBROUTINE pse_check_ctrl ppm_pse/pse_global.f0000644000175000017500000003025511451402172013457 0ustar omaromar !------------------------------------------------------------------------- ! Module : pse_global !------------------------------------------------------------------------- ! ! Purpose : Declare global types and variables for ppm_pse. ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- MODULE pse_global !---------------------------------------------------------------------- ! Modules !---------------------------------------------------------------------- USE ppm_module_substart USE ppm_module_substop USE ppm_module_typedef USE ppm_module_neighlist, ONLY: ppm_type_ptr_to_clist PRIVATE :: ppm_type_ptr_to_clist !---------------------------------------------------------------------- ! Global TYPEs !---------------------------------------------------------------------- TYPE ptr_to_list INTEGER, DIMENSION(:), POINTER :: list END TYPE !---------------------------------------------------------------------- ! Header file !---------------------------------------------------------------------- INCLUDE 'ppm_numerics.h' INCLUDE 'pse_param.h' !---------------------------------------------------------------------- ! Positions, strengths, strength deltas, and particle costs !---------------------------------------------------------------------- REAL(MK), DIMENSION(:,:), POINTER :: xp,wp,dwp REAL(MK), DIMENSION(: ), POINTER :: pcost !---------------------------------------------------------------------- ! Topology !---------------------------------------------------------------------- REAL(MK), DIMENSION(3) :: min_phys,max_phys REAL(MK), DIMENSION(3) :: min_physg,max_physg REAL(MK), DIMENSION(: ), POINTER :: proc_speed REAL(MK), DIMENSION(: ), POINTER :: sub_cost INTEGER , DIMENSION(:,:), POINTER :: istart INTEGER , DIMENSION(:,:), POINTER :: ndata TYPE(ppm_t_topo), POINTER :: topo !---------------------------------------------------------------------- ! Mollifier kernel size !---------------------------------------------------------------------- REAL(MK) :: eps !---------------------------------------------------------------------- ! Domain decomposition !---------------------------------------------------------------------- INTEGER :: ppm_decomp !---------------------------------------------------------------------- ! Number of particles in each direction for internal initialization !---------------------------------------------------------------------- INTEGER, DIMENSION(3) :: Ngrid !---------------------------------------------------------------------- ! PSE kernel type and parameters !---------------------------------------------------------------------- INTEGER :: kernel INTEGER :: ppm_kernel REAL(MK), DIMENSION(:), ALLOCATABLE :: kpar !---------------------------------------------------------------------- ! Inter-particle distance and particle volume !---------------------------------------------------------------------- REAL(MK) :: dhx,dhy,dhz REAL(MK) :: kappa REAL(MK) :: pvol,pvolinv !---------------------------------------------------------------------- ! PP interaction cutoff !---------------------------------------------------------------------- REAL(MK) :: cutoff,cutoff2 !---------------------------------------------------------------------- ! Skin for Verlet lists !---------------------------------------------------------------------- REAL(MK) :: skin !---------------------------------------------------------------------- ! Neighbor list type: 1=Verlet, 2=Cell, 3=None !---------------------------------------------------------------------- INTEGER :: neighlist !---------------------------------------------------------------------- ! Number of particles on local processor !---------------------------------------------------------------------- ! without ghosts (i.e. real particles) INTEGER :: Npart ! total particles including the ghosts INTEGER :: Mpart ! total number of particles in ROI INTEGER :: roinp !---------------------------------------------------------------------- ! Surface triangulation !---------------------------------------------------------------------- INTEGER :: ntri REAL(MK), DIMENSION(:,:,:), POINTER :: triang !---------------------------------------------------------------------- ! Bin lists of triangles !---------------------------------------------------------------------- TYPE(ptr_to_list), DIMENSION(:,:),ALLOCATABLE :: bin INTEGER, DIMENSION(:,:), ALLOCATABLE :: nbin INTEGER :: NBINY,NBINZ REAL(MK) :: dby,dbz !---------------------------------------------------------------------- ! Geometry tolerance !---------------------------------------------------------------------- REAL(MK) :: lmyeps !---------------------------------------------------------------------- ! Inside/outside flags for all grid points/sub-domains !---------------------------------------------------------------------- INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: flags !---------------------------------------------------------------------- ! Region of interest for diagnostics output !---------------------------------------------------------------------- REAL(MK), DIMENSION(3) :: min_roi,max_roi !---------------------------------------------------------------------- ! Time stepping !---------------------------------------------------------------------- REAL(MK) :: dt INTEGER :: nstep !---------------------------------------------------------------------- ! Number of stages in the STS time stepper and kick-in time !---------------------------------------------------------------------- INTEGER :: NSTS REAL(MK) :: STSstart !---------------------------------------------------------------------- ! Diffusion constants for all species !---------------------------------------------------------------------- REAL(MK), DIMENSION(: ), ALLOCATABLE :: diff !---------------------------------------------------------------------- ! Math. constant !---------------------------------------------------------------------- REAL(MK) :: pi,pi2 !---------------------------------------------------------------------- ! Number of species !---------------------------------------------------------------------- INTEGER :: NSPEC !---------------------------------------------------------------------- ! Use of symmetry !---------------------------------------------------------------------- LOGICAL :: symmetry !---------------------------------------------------------------------- ! Number of dimensions !---------------------------------------------------------------------- INTEGER :: ndim !---------------------------------------------------------------------- ! Program control flags !---------------------------------------------------------------------- LOGICAL :: checkmap LOGICAL :: checkunique LOGICAL :: probeproc LOGICAL :: dumpparticles LOGICAL :: rebalance !---------------------------------------------------------------------- ! Kernel lookup table. lookup determines if it is used. !---------------------------------------------------------------------- LOGICAL :: lookup REAL(MK), DIMENSION(: ), POINTER :: ktab REAL(MK) :: drtableinv INTEGER :: ntable !---------------------------------------------------------------------- ! I/O and control !---------------------------------------------------------------------- INTEGER :: freqoutput,freqdiag INTEGER :: freqrestart CHARACTER(LEN=MAXCHAR) :: logfile,ctrlfile CHARACTER(LEN=MAXCHAR) :: inputfile,abortfile CHARACTER(LEN=MAXCHAR) :: outputtab,outputfile CHARACTER(LEN=MAXCHAR) :: inputfmt,outputfmt CHARACTER(LEN=MAXCHAR) :: rstfile,diagfile CHARACTER(LEN=MAXCHAR) :: prgname,diagfmt CHARACTER(LEN=MAXCHAR) :: casename LOGICAL :: restart INTEGER :: debug INTEGER :: ppm_log_unit INTEGER :: pse_log_unit !---------------------------------------------------------------------- ! MPI stuff !---------------------------------------------------------------------- INTEGER :: rank,Nproc,comm INTEGER :: MPPREC !---------------------------------------------------------------------- ! Verlet and cell list memory !---------------------------------------------------------------------- ! number of cells in all directions INTEGER, DIMENSION(:,:), POINTER :: Nm ! cell list TYPE(ppm_type_ptr_to_clist), DIMENSION(:), POINTER :: clist ! Verlet lists INTEGER, DIMENSION(:,:), POINTER :: vlist ! lengths of Verlet lists INTEGER, DIMENSION( :), POINTER :: nvlist END MODULE pse_global ppm_pse/box.surf0000644000175000017500000000140311445416204012667 0ustar omaromar#ASCII triangulated surface description Surface description for a cube of size 4x4x4 ============================================ M: total number of triangles ---------------------------- 12 Trianges oriented such that P1P2xP1P3 points outward ---------------------------------------------------- P1x P1y P1z P2x P2y P2z P3x P3y P3z 0.0 0.0 0.0 4.0 0.0 0.0 4.0 0.0 4.0 0.0 0.0 0.0 4.0 0.0 4.0 0.0 0.0 4.0 4.0 0.0 0.0 4.0 4.0 0.0 4.0 4.0 4.0 4.0 0.0 0.0 4.0 4.0 4.0 4.0 0.0 4.0 0.0 0.0 4.0 4.0 0.0 4.0 0.0 4.0 4.0 0.0 4.0 4.0 4.0 0.0 4.0 4.0 4.0 4.0 0.0 0.0 0.0 0.0 0.0 4.0 0.0 4.0 4.0 0.0 0.0 0.0 0.0 4.0 4.0 0.0 4.0 0.0 0.0 0.0 0.0 4.0 4.0 0.0 4.0 0.0 0.0 0.0 0.0 0.0 0.0 4.0 0.0 4.0 4.0 0.0 0.0 4.0 0.0 0.0 4.0 4.0 4.0 4.0 0.0 0.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 0.0 ppm_pse/pwrite.f0000644000175000017500000001017011445416206012662 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pwrite !------------------------------------------------------------------------- ! ! Purpose : Subroutine for parallel writing for debugging purposes. ! ! Input : rank (I) MPI rank of writing processor ! caller (C) name of calling subroutine ! cbuf (C) message ! ! Input/output : ! ! Output : info (I) error status ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pwrite.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.1 2004/06/02 15:28:38 ivos ! Initial implementation. ! !------------------------------------------------------------------------- ! Jens Honore Walther ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pwrite(rank,caller,cbuf,info) IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER :: rank CHARACTER(LEN=*) :: caller,cbuf INTEGER :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- CHARACTER(LEN=256) :: cformat INTEGER :: ios INTEGER :: icaller,ibuf !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Get length of messages !------------------------------------------------------------------------- icaller = LEN_TRIM(caller) ibuf = LEN_TRIM(cbuf) !------------------------------------------------------------------------- ! Define the print format !------------------------------------------------------------------------- IF (rank.LT.0) THEN cformat = '(4A)' ELSEIF (rank.LT.10) THEN cformat = '(A,I1,4A)' ELSEIF (rank.LT.100) THEN cformat = '(A,I2,4A)' ELSEIF (rank.LT.1000) THEN cformat = '(A,I3,4A)' ELSE cformat = '(A,I4,4A)' ENDIF !------------------------------------------------------------------------- ! Do the print !------------------------------------------------------------------------- IF (rank.LT.0) THEN WRITE(*,cformat,IOSTAT=ios) & & '(', & & caller(1:icaller), & & ') : ', & & cbuf(1:ibuf) ELSE WRITE(*,cformat,IOSTAT=ios) & & '[',rank,'] (', & & caller(1:icaller), & & ') : ', & & cbuf(1:ibuf) ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE RETURN END SUBROUTINE pwrite ppm_pse/pse_defaults.f0000644000175000017500000002434711451377704014047 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_defaults !------------------------------------------------------------------------- ! ! Purpose : Define default/initial values for some of the ! global variables. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_defaults(info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- #ifdef __MPI INCLUDE 'mpif.h' #endif !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- REAL(8) :: t0 !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_defaults',t0,info) !------------------------------------------------------------------------- ! Default run name !------------------------------------------------------------------------- casename = 'ppm_pse' !------------------------------------------------------------------------- ! Default Ctrl file !------------------------------------------------------------------------- ctrlfile = 'Ctrl' !------------------------------------------------------------------------- ! Log file units !------------------------------------------------------------------------- ppm_log_unit = 99 pse_log_unit = 88 !------------------------------------------------------------------------- ! No surface triangulation !------------------------------------------------------------------------- ntri = 0 NULLIFY(triang) !------------------------------------------------------------------------- ! Set geometry tolerance !------------------------------------------------------------------------- IF (MK .EQ. KIND(1.0D0)) THEN lmyeps = 1D-10 ELSE lmyeps = 1E-5 ENDIF !------------------------------------------------------------------------- ! Use 100 bins if not specified !------------------------------------------------------------------------- NBINY = 100 NBINZ = 100 !------------------------------------------------------------------------- ! Default restart file !------------------------------------------------------------------------- rstfile = 'restart.dat' !------------------------------------------------------------------------- ! Default input filename (actually empty, thus no default is offered) !------------------------------------------------------------------------- inputfile = ' ' inputfmt = 'FORMATTED' outputtab = ' ' outputfile= ' ' outputfmt = 'FORMATTED' !------------------------------------------------------------------------- ! Default particle grid size if no input file is given !------------------------------------------------------------------------- Ngrid(1) = 32 Ngrid(2) = 32 Ngrid(3) = 32 !------------------------------------------------------------------------- ! Default STS scheme uses 1 stage (Euler) and starts at time 0 !------------------------------------------------------------------------- NSTS = 1 STSstart = 0.0_MK !------------------------------------------------------------------------- ! Default decomposition type: binary tree !------------------------------------------------------------------------- ppm_decomp = -1 !------------------------------------------------------------------------- ! Default log filename !------------------------------------------------------------------------- logfile = 'pse.log' !------------------------------------------------------------------------- ! Default diagnostics filename !------------------------------------------------------------------------- diagfile = '' !------------------------------------------------------------------------- ! Default computational domain is unit cude !------------------------------------------------------------------------- min_phys(1:3) = 0.0_MK max_phys(1:3) = 1.0_MK !------------------------------------------------------------------------- ! Default region of interest is same as comput. domain !------------------------------------------------------------------------- min_phys(1:3) = 0.0_MK max_phys(1:3) = 1.0_MK !------------------------------------------------------------------------- ! Discretisation correction factor !------------------------------------------------------------------------- kappa = 1.0_MK !------------------------------------------------------------------------- ! Kernel type defaults to 2nd order polynomial (Cottet) !------------------------------------------------------------------------- kernel = 1 !------------------------------------------------------------------------- ! Kernel core size !------------------------------------------------------------------------- eps = 0.0_MK !------------------------------------------------------------------------- ! Inter-particle distances !------------------------------------------------------------------------- dhx = 0.0_MK dhy = 0.0_MK dhz = 0.0_MK !------------------------------------------------------------------------- ! Default debug value !------------------------------------------------------------------------- debug = 0 !------------------------------------------------------------------------- ! Default uses symmetry !------------------------------------------------------------------------- symmetry = .TRUE. !------------------------------------------------------------------------- ! Default uses no cutoff !------------------------------------------------------------------------- cutoff = -1 !------------------------------------------------------------------------- ! Default uses no neighbor lists !------------------------------------------------------------------------- neighlist = 3 !------------------------------------------------------------------------- ! Default is 3D simulation !------------------------------------------------------------------------- ndim = 3 !------------------------------------------------------------------------- ! Default: do not check mapping results, do not check particle ! uniqueness, do not probe processor speeds, do not dump part posns, do ! not rebalance computational load. !------------------------------------------------------------------------- rebalance = .FALSE. !------------------------------------------------------------------------- ! Default is no restart !------------------------------------------------------------------------- restart = .FALSE. !------------------------------------------------------------------------- ! Default time stepping values !------------------------------------------------------------------------- nstep = 0 !------------------------------------------------------------------------- ! Default output frequencies !------------------------------------------------------------------------- freqdiag = 1 freqoutput = 1 freqrestart = 100 !------------------------------------------------------------------------- ! Lookup table is not used by default !------------------------------------------------------------------------- lookup = .FALSE. !------------------------------------------------------------------------- ! Default number of entries in table loop-up (if used) !------------------------------------------------------------------------- ntable = 262145 !------------------------------------------------------------------------- ! Default number of species and diffusion constant !------------------------------------------------------------------------- NSPEC = 1 !------------------------------------------------------------------------- ! Abort file !------------------------------------------------------------------------- abortfile = 'ABORT' !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_defaults',t0,info) RETURN END SUBROUTINE pse_defaults ppm_pse/pse_point_inside.f0000644000175000017500000002250511451403662014707 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_point_inside !------------------------------------------------------------------------- ! ! Purpose : This routine checks if the point given by p(1:3) is ! inside the closes surface described by triang(...). ! flag returns -1 if the point is inside and 1 if it ! is outside. ! ! Input : p(3) (F) Point to test. ! ! Input/output : ! ! Output : flags (I) returns 1 if point is in domain, 0 ! otherwise. ! info (I) return status (0 if no error) ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_point_inside.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:27 ivos ! CBL version of the PPM library ! ! Revision 1.1 2005/04/27 01:19:12 ivos ! Initial commit. Tested. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum,Hirschengraben 84 ! CH-8092 Zurich,Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_point_inside(p,flag,info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- REAL(MK), DIMENSION(3) ,INTENT(IN ) :: p INTEGER ,INTENT( OUT) :: flag,info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER :: i, j, nip REAL(MK) :: phi REAL(MK), DIMENSION(3) :: a, b REAL(MK), DIMENSION(3,3) :: mat INTEGER :: idxy, idxz CHARACTER(LEN=256) :: cbuf !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Follow x direction and intersect with all triangles !------------------------------------------------------------------------- nip = 0 mat(1,1) = 0.0 mat(2,1) = 0.0 idxy = INT((p(2)-min_physg(2))/dby) + 1 IF (idxy .LT. 1) idxy = 1 ! account for boundary points IF (idxy .GT. NBINY) idxy = NBINY idxz = INT((p(3)-min_physg(3))/dbz) + 1 IF (idxz .LT. 1) idxz = 1 ! account for boundary points IF (idxz .GT. NBINZ) idxz = NBINZ IF (debug .GT. 5) THEN WRITE(cbuf,'(A,2I4)') 'bin indices: ',idxy,idxz CALL pwrite(rank,'pse_point_inside',cbuf,info) WRITE(cbuf,'(A,I8)') 'nbin=',nbin(idxy,idxz) CALL pwrite(rank,'pse_point_inside',cbuf,info) OPEN(70,FILE='bin.deb',STATUS='replace') WRITE(70,*) min_physg(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,*) min_physg(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(69,*) min_physg(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz)*dbz WRITE(70,*) min_physg(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz)*dbz WRITE(70,*) min_physg(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,'(A)') WRITE(70,'(A)') WRITE(70,*) max_phys(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,'(A)') WRITE(70,'(A)') WRITE(70,*) min_physg(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,'(A)') WRITE(70,'(A)') WRITE(70,*) min_physg(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz-1)*dbz WRITE(70,'(A)') WRITE(70,'(A)') WRITE(70,*) min_physg(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy)*dby, & & min_physg(3)+(idxz)*dbz WRITE(70,'(A)') WRITE(70,'(A)') WRITE(70,*) min_physg(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz)*dbz WRITE(70,*) max_phys(1), min_physg(2)+(idxy-1)*dby, & & min_physg(3)+(idxz)*dbz CLOSE(70) ENDIF !------------------------------------------------------------------------- ! Intersect with all triangles in the same bin !------------------------------------------------------------------------- DO j=1,nbin(idxy,idxz) ! get index of next triangle i = bin(idxy,idxz)%list(j) IF (debug .GT. 5) THEN OPEN(70,FILE='tlist.deb',POSITION='append') WRITE(70,*) triang(:,1,i) WRITE(70,*) triang(:,2,i) WRITE(70,*) triang(:,3,i) WRITE(70,*) triang(:,1,i) WRITE(70,'(A)') WRITE(70,'(A)') CLOSE(70) ENDIF ! edge vectors a(1) = triang(1,2,i)-triang(1,1,i) a(2) = triang(2,2,i)-triang(2,1,i) a(3) = triang(3,2,i)-triang(3,1,i) b(1) = triang(1,3,i)-triang(1,1,i) b(2) = triang(2,3,i)-triang(2,1,i) b(3) = triang(3,3,i)-triang(3,1,i) ! determinant of system matrix phi = -a(3)*b(2)+a(2)*b(3) ! for v=(1,0,0) IF (ABS(phi) .GT. lmyeps) THEN ! if intersection point exists phi = 1.0_MK/phi ! inverse of system matrix mat(1,2) = phi*b(3) mat(1,3) = -phi*b(2) mat(2,2) = -phi*a(3) mat(2,3) = phi*a(2) mat(3,1) = phi*(-a(3)*b(2)+a(2)*b(3)) mat(3,2) = phi*(a(3)*b(1)-a(1)*b(3)) mat(3,3) = phi*(-a(2)*b(1)+a(1)*b(2)) ! right hand side a = p-triang(:,1,i) ! solution of 3x3 system b = (mat(:,1)*a(1))+(mat(:,2)*a(2))+(mat(:,3)*a(3)) !------------------------------------------------------------------- ! Intersect inside triangle ? !------------------------------------------------------------------- IF (.NOT.(b(1)+b(2) .GT. 1.0_MK) .AND. .NOT.(b(1) .LT. 0.0_MK) & .AND. (b(2) .GT. lmyeps) .AND. (b(3) .GT. lmyeps)) THEN IF (debug .GT. 5) THEN WRITE(cbuf,'(A,I6)') 'intersected triangle ',i CALL pwrite(rank,'pse_point_inside',cbuf,info) OPEN(70,file='triang.deb',position='append') WRITE(70,*) triang(:,1,i) WRITE(70,*) triang(:,2,i) WRITE(70,*) triang(:,3,i) WRITE(70,*) triang(:,1,i) WRITE(70,'(A)') WRITE(70,'(A)') CLOSE(70) OPEN(70,FILE='points.deb',POSITION='append') WRITE(70,*) triang(:,1,i)+b(1)*(triang(:,2,i)-triang(:,1,i)) & +b(2)*(triang(:,3,i)-triang(:,1,i)) CLOSE(70) WRITE(70,*) b CLOSE(70) ENDIF nip=nip+1 ENDIF ENDIF ENDDO !------------------------------------------------------------------------- ! Odd or even number of intersections ? !------------------------------------------------------------------------- IF (MOD(nip,2) .EQ. 0) THEN ! outside flag = 1 ELSE ! inside flag = -1 ENDIF IF (debug .GT. 5) THEN WRITE(cbuf,'(A,I2)') 'point_inside=',flag CALL pwrite(rank,'pse_point_inside',cbuf,info) ENDIF !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE RETURN END SUBROUTINE pse_point_inside ppm_pse/Make.gcc0000644000175000017500000000203311455632453012540 0ustar omaromar #---------------------------------------------------------------------- # PATHS (need to be adjusted) #---------------------------------------------------------------------- # PPM Library PPMDIR = ../libppm/ PPMNUMDIR = ../libppmnumerics/ # METIS Library METISDIR = /home/user/sw/metis/gcc # MPI Library MPIDIR = /opt/openmpi1.4.3/ # C++ Compiler libraries CCDIR = /usr/lib #---------------------------------------------------------------------- # Compiler settings #---------------------------------------------------------------------- F77 = gfortran FC = gfortran CC = gcc INC = -I$(PPMDIR)/include -I$(PPMNUMDIR)/include #DEBUG = -C -g #FFLAGS = -ffree-form -O3 $(INC) # change FFLAGS for debugging purpose FFLAGS = -ffree-form -O3 $(INC) $(DEBUG) CFLAGS = -O LD = $(F77) LDFLAGS = $(FFLAGS) -static CPP = cpp LIBS = -L$(PPMDIR)/lib -L$(PPMNUMDIR)/lib -lppm -lppmnumerics -L$(METISDIR)/lib -lmetis -L$(MPIDIR)/lib LDFLAGS = $(LIBS) ppm_pse/Makefile0000644000175000017500000001413011455632520012641 0ustar omaromar#---------------------------------------------------------------------- # Makefile for PSE client #---------------------------------------------------------------------- SHELL = /bin/sh TARGET = ppm_pse BIN = ppm_pse #---------------------------------------------------------------------- # Compiler settings (Linux) #---------------------------------------------------------------------- include Make.ifc #DEFINE = -D__MPI -D__NEWMPI DEFINE = #---------------------------------------------------------------------- # Search path for RCS files #---------------------------------------------------------------------- VPATH = ./RCS #---------------------------------------------------------------------- # Additional suffix rules #---------------------------------------------------------------------- .SUFFIXES : .inc .inc,v .f,v .c,v .f,v.f : co $*.f .c,v.c : co $*.c .inc,v.inc : co $*.inc #---------------------------------------------------------------------- # Binary directory #---------------------------------------------------------------------- bindir = . #---------------------------------------------------------------------- # Default target #---------------------------------------------------------------------- all: $(TARGET) #---------------------------------------------------------------------- # Object files: #---------------------------------------------------------------------- OBJS = pse_global.o\ pse_constants.o\ pse_module_io.o\ pse_module_comp.o\ pse_defaults.o\ pwrite.o\ pse_check_abort.o\ pse_check_ctrl.o\ UpperCase.o\ pse_time.o\ pse_time2.o\ pse_module_realloc.o\ pse_alloc_binlist.o\ pse_sort_triangles.o\ pse_point_inside.o\ pse_get_indicator.o\ pse_arguments.o\ pse_dump_parameters.o\ pse_create_particles.o\ pse_init.o\ pse_finalize.o\ ppm_pse.o #---------------------------------------------------------------------- # Dependencies: #---------------------------------------------------------------------- pse_global.o: pse_global.f $(F77) $(FFLAGS) -c pse_global.f pse_module_realloc.o: pse_module_realloc.f pse_global.o $(CPP) $(DEFINE) -P pse_module_realloc.f __pse_module_realloc.f $(F77) $(FFLAGS) -c __pse_module_realloc.f -o pse_module_realloc.o pse_alloc_binlist.o: pse_alloc_binlist.f pse_module_realloc.o $(CPP) $(DEFINE) -P pse_alloc_binlist.f __pse_alloc_binlist.f $(F77) $(FFLAGS) -c __pse_alloc_binlist.f -o pse_alloc_binlist.o pse_sort_triangles.o: pse_sort_triangles.f pse_global.o pse_module_realloc.o $(CPP) $(DEFINE) -P pse_sort_triangles.f __pse_sort_triangles.f $(F77) $(FFLAGS) -c __pse_sort_triangles.f -o pse_sort_triangles.o pse_point_inside.o: pse_point_inside.f pse_global.o $(CPP) $(DEFINE) -P pse_point_inside.f __pse_point_inside.f $(F77) $(FFLAGS) -c __pse_point_inside.f -o pse_point_inside.o pse_get_indicator.o: pse_get_indicator.f pse_global.o $(CPP) $(DEFINE) -P pse_get_indicator.f __pse_get_indicator.f $(F77) $(FFLAGS) -c __pse_get_indicator.f -o pse_get_indicator.o pse_create_particles.o: pse_create_particles.f pse_global.o $(CPP) $(DEFINE) -P pse_create_particles.f __pse_create_particles.f $(F77) $(FFLAGS) -c __pse_create_particles.f -o pse_create_particles.o pse_module_io.o: pse_module_io.f pse_read_ctrl.f \ pse_write_output.f \ pse_write_diag.f pse_read_geometry.f $(CPP) $(DEFINE) -P pse_module_io.f __pse_module_io.f $(F77) $(FFLAGS) -c __pse_module_io.f -o pse_module_io.o pse_module_comp.o: pse_module_comp.f pse_rhs.f $(CPP) $(DEFINE) -P pse_module_comp.f __pse_module_comp.f $(F77) $(FFLAGS) -c __pse_module_comp.f -o pse_module_comp.o ppm_pse.o: ppm_pse.f pse_module_io.f pse_module_comp.f $(CPP) $(DEFINE) -P ppm_pse.f __ppm_pse.f $(F77) $(FFLAGS) -c __ppm_pse.f -o ppm_pse.o pse_time.o: pse_time.f $(CPP) $(DEFINE) -P pse_time.f __pse_time.f $(F77) $(FFLAGS) -c __pse_time.f -o pse_time.o pse_time2.o: pse_time2.f $(CPP) $(DEFINE) -P pse_time2.f __pse_time2.f $(F77) $(FFLAGS) -c __pse_time2.f -o pse_time2.o pse_init.o: pse_init.f pse_module_io.f $(CPP) $(DEFINE) -P pse_init.f __pse_init.f $(F77) $(FFLAGS) -c __pse_init.f -o pse_init.o pse_finalize.o: pse_finalize.f $(CPP) $(DEFINE) -P pse_finalize.f __pse_finalize.f $(F77) $(FFLAGS) -c __pse_finalize.f -o pse_finalize.o pse_arguments.o: pse_arguments.f $(CPP) $(DEFINE) -P pse_arguments.f __pse_arguments.f $(F77) $(FFLAGS) -c __pse_arguments.f -o pse_arguments.o pse_defaults.o: pse_defaults.f $(CPP) $(DEFINE) -P pse_defaults.f __pse_defaults.f $(F77) $(FFLAGS) -c __pse_defaults.f -o pse_defaults.o pwrite.o: pwrite.f $(F77) $(FFLAGS) -c pwrite.f pse_constants.o: pse_constants.f $(F77) $(FFLAGS) -c pse_constants.f pse_dump_parameters.o: pse_dump_parameters.f $(F77) $(FFLAGS) -c pse_dump_parameters.f pse_check_abort.o: pse_check_abort.f $(CPP) $(DEFINE) -P pse_check_abort.f __pse_check_abort.f $(F77) $(FFLAGS) -c __pse_check_abort.f -o pse_check_abort.o pse_check_ctrl.o: pse_check_ctrl.f $(CPP) $(DEFINE) -P pse_check_ctrl.f __pse_check_ctrl.f $(F77) $(FFLAGS) -c __pse_check_ctrl.f -o pse_check_ctrl.o UpperCase.o: UpperCase.f $(F77) $(FFLAGS) -c UpperCase.f #---------------------------------------------------------------------- # link #---------------------------------------------------------------------- $(TARGET): $(OBJS) $(LD) $(OBJS) -o $(BIN) $(LDFLAGS) $(LIBS) #---------------------------------------------------------------------- # Install #---------------------------------------------------------------------- install: $(BIN) (cp -f $(BIN) $(bindir)) #---------------------------------------------------------------------- # Clean #---------------------------------------------------------------------- new: clean ppm_pse clean: rm -f __*.f rm -f $(OBJS) rm -f *.lst rm -f *.l rm -f *.L ppm_pse/NOTICE0000644000175000017500000000114211445416204012102 0ustar omaromarLIABILITY NOTICE IN NO EVENT SHALL THE ETH BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE ETH HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. THE ETH SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE ETH HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. ppm_pse/README0000644000175000017500000000051511451412651012060 0ustar omaromarPPM PSE code ============ This is a stripped down version of the original ppm_pse client. For the full code please see: http://www.mosaic.ethz.ch/Downloads/ppm_pse You may use this code as an example or starting point to develop your own PPM client applications. For questions please contact Omar Awile (omar.awile@inf.ethz.ch). ppm_pse/pse_sort_triangles.f0000644000175000017500000003153311451403336015261 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_sort_triangles !------------------------------------------------------------------------- ! ! Purpose : This routine sets up the bin lists by assigning the ! triangles to bins of equal y and z coordinates. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status (0 if no error) ! ! Routines : ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! $Log: pse_sort_triangles.f,v $ ! Revision 1.1.1.1 2007/07/13 10:18:16 ivos ! CBL version of the PPM library ! ! Revision 1.1 2005/04/27 01:19:14 ivos ! Initial commit. Tested. ! !------------------------------------------------------------------------- ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum,Hirschengraben 84 ! CH-8092 Zurich,Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_sort_triangles(info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global USE pse_module_realloc IMPLICIT NONE !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER ,INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER :: i, j, k, p, ix1, ix2 INTEGER :: minbiny INTEGER :: maxbiny, minbinz, maxbinz INTEGER :: idxy, idxz REAL(MK), DIMENSION(4) :: bb REAL(MK), DIMENSION(2) :: a, b, x1, y1, v1, v2 REAL(MK) :: det, det2 REAL(MK) :: alpha, beta LOGICAL :: NotIn REAL(MK), DIMENSION(2,3) :: edge REAL(MK), DIMENSION(2,2) :: side REAL(8) :: t0 !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_sort_triangles',t0,info) !------------------------------------------------------------------------- ! Initialise the bin lists !------------------------------------------------------------------------- CALL pse_alloc_binlist(0,0,0,info) IF (info .NE. 0) THEN CALL pwrite(rank,'pse_sort_triangles','Failed to allocate bin lists',& & info) info = -1 GOTO 9999 ENDIF dby = (max_physg(2)-min_physg(2))/REAL(NBINY,MK) dbz = (max_physg(3)-min_physg(3))/REAL(NBINZ,MK) !------------------------------------------------------------------------- ! Assign triangles to linear bins !------------------------------------------------------------------------- !CALL pwrite(rank,'pse_sort_triangles','Sorting triangles into bins...',& !& info) DO i=1,ntri ! loop over all triangles minbiny = NBINY+1 maxbiny = 0 minbinz = NBINZ+1 maxbinz = 0 x1 = triang(2:3,1,i) edge(:,1) = triang(2:3,2,i)-x1 edge(:,2) = triang(2:3,3,i)-x1 edge(:,3) = triang(2:3,3,i)-triang(2:3,2,i) a = edge(:,1) b = edge(:,2) det = a(1)*b(2)-a(2)*b(1) det = 1.0_MK/det !---------------------------------------------------------------------- ! First look at the triangle vertices !---------------------------------------------------------------------- DO j=1,3 ! loop over all vertices !------------------------------------------------------------------- ! Determine index of bin of this vertex !------------------------------------------------------------------- idxy = INT((triang(2,j,i)-min_physg(2))/dby) + 1 IF (idxy .EQ. 0) idxy = 1 ! include boundaries at walls IF (idxy .EQ. NBINY+1) idxy = NBINY ! skip triangles outside the extended domain IF (idxy .LT. 1) CYCLE IF (idxy .GT. NBINY) CYCLE idxz = INT((triang(3,j,i)-min_physg(3))/dbz) + 1 IF (idxz .EQ. 0) idxz = 1 ! include boundaries at walls IF (idxz .EQ. NBINZ+1) idxz = NBINZ ! skip triangles outside the extended domain IF (idxz .LT. 1) CYCLE IF (idxz .GT. NBINZ) CYCLE !------------------------------------------------------------------- ! Add triangle to this bin if not already there !------------------------------------------------------------------- IF (.NOT. ASSOCIATED(bin(idxy,idxz)%list)) THEN CALL pse_alloc_binlist(1,idxy,idxz,info) IF (info .NE. 0) GOTO 9999 ENDIF NotIn = .FALSE. IF (nbin(idxy,idxz) .EQ. 0) THEN NotIn = .TRUE. ELSEIF (bin(idxy,idxz)%list(nbin(idxy,idxz)) .NE. i) THEN NotIn = .TRUE. ENDIF IF (NotIn) THEN p = nbin(idxy,idxz)+1 nbin(idxy,idxz) = p ! enlarge this list if nedded IF (p .GT. SIZE(bin(idxy,idxz)%list)) THEN CALL pse_alloc_binlist(1,idxy,idxz,info) IF (info .NE. 0) GOTO 9999 ENDIF bin(idxy,idxz)%list(p) = i ! update bounding box of cells IF (idxy .LT. minbiny) minbiny = idxy IF (idxy .GT. maxbiny) maxbiny = idxy IF (idxz .LT. minbinz) minbinz = idxz IF (idxz .GT. maxbinz) maxbinz = idxz ENDIF ENDDO !---------------------------------------------------------------------- ! Loop over all bins in the bounding box of the vertices !---------------------------------------------------------------------- DO j=minbiny,maxbiny bb(1) = min_physg(2)+(j-1)*dby ! lower y bb(2) = min_physg(2)+j*dby ! upper y DO k=minbinz,maxbinz !---------------------------------------------------------------- ! If no list for this bin is associated, create one. !---------------------------------------------------------------- IF (.NOT. ASSOCIATED(bin(j,k)%list)) THEN CALL pse_alloc_binlist(1,j,k,info) IF (info .NE. 0) GOTO 9999 ENDIF NotIn = .FALSE. !---------------------------------------------------------------- ! If this triangle does not already belong to this bin, ! check it. !---------------------------------------------------------------- IF (nbin(j,k) .EQ. 0) THEN NotIn = .TRUE. ELSEIF (bin(j,k)%list(nbin(j,k)) .NE. i) THEN NotIn = .TRUE. ENDIF IF (NotIn) THEN bb(3) = min_physg(3)+(k-1)*dbz ! lower z bb(4) = min_physg(3)+k*dbz ! upper z ! edge vectors of this bin (projected onto y/z plane) ! only store two since the other two are identical side(1,1) = bb(2)-bb(1) side(2,1) = 0.0_MK side(1,2) = 0.0_MK side(2,2) = bb(4)-bb(3) DO ix1=1,2 y1(1) = bb(ix1)-x1(1) DO ix2=3,4 y1(2) = bb(ix2)-x1(2) alpha = det*(b(2)*y1(1)-b(1)*y1(2)) beta = det*(-a(2)*y1(1)+a(1)*y1(2)) IF ((alpha .GE. 0.0_MK) .AND. (beta .GE. 0.0_MK) & .AND. (alpha+beta .LE. 1.0_MK)) THEN ! point is inside triangle: add triangle to this bin p = nbin(j,k)+1 nbin(j,k) = p ! enlarge this list if nedded IF (p .GT. SIZE(bin(j,k)%list)) THEN CALL pse_alloc_binlist(1,j,k,info) IF (info .NE. 0) GOTO 9999 ENDIF bin(j,k)%list(p) = i GOTO 250 ! next bin ENDIF ENDDO ENDDO !------------------------------------------------------------- ! If no vertex of the bin is inside the triangle, check ! if any edge of the triangle intersects any edge of the ! bin. !------------------------------------------------------------- DO ix1=1,3 ! loop over triangles edges v1 = edge(:,ix1) DO ix2=1,4 ! loop over bins edges ! choose adge vector MOD 2 since bin is axis parallel v2 = -side(:,MOD(ix2,2)+1) ! system determinant (DO not overwrite det !!) det2 = v1(1)*v2(2)-v1(2)*v2(1) det2 = 1.0_MK/det2 ! determine rhs (starting points of lines depending ! on what edge we are at IF (ix2 .EQ. 1) THEN y1 = (/bb(1), bb(3)/) ELSEIF (ix2 .EQ. 2) THEN y1 = (/bb(1), bb(3)/) ELSEIF (ix2 .EQ. 3) THEN y1 = (/bb(2), bb(3)/) ELSEIF (ix2 .EQ. 4) THEN y1 = (/bb(1), bb(4)/) ENDIF IF (ix1 .EQ. 1 .OR. ix1 .EQ. 2) THEN y1 = y1 - x1 ELSE y1 = y1 - triang(2:3,2,i) ENDIF ! solve for intersection point alpha = det2*(v2(2)*y1(1)-v2(1)*y1(2)) beta = det2*(-v1(2)*y1(1)+v1(1)*y1(2)) !------------------------------------------------------- ! If intersected, add triangle to this bin !------------------------------------------------------- IF ((alpha .GE. 0.0_MK) .AND. (beta .GE. 0.0_MK) .AND. & (alpha .LE. 1.0_MK) .AND. (beta .LE. 1.0_MK)) THEN ! point is inside triangle: add triangle to this bin p = nbin(j,k)+1 nbin(j,k) = p ! enlarge this list if nedded IF (p .GT. SIZE(bin(j,k)%list)) THEN CALL pse_alloc_binlist(1,j,k,info) IF (info .NE. 0) GOTO 9999 ENDIF bin(j,k)%list(p) = i GOTO 250 ! next bin ENDIF ENDDO ENDDO 250 ENDIF ENDDO ENDDO ENDDO !------------------------------------------------------------------------- ! prune bin lists to actual SIZE to save memory !------------------------------------------------------------------------- DO i=1,NBINY DO j=1,NBINZ bin(i,j)%list=>reallocate(bin(i,j)%list,nbin(i,j)) ENDDO ENDDO !CALL pwrite(rank,'pse_sort_triangles','Done sorting triangles.',info) !------------------------------------------------------------------------- ! Return !------------------------------------------------------------------------- 9999 CONTINUE CALL substop('pse_sort_triangles',t0,info) RETURN END SUBROUTINE pse_sort_triangles ppm_pse/pse_arguments.f0000644000175000017500000002171411455355615014241 0ustar omaromar !------------------------------------------------------------------------- ! Subroutine : pse_arguments !------------------------------------------------------------------------- ! ! Purpose : Read and process command line arguments and print ! defines. ! ! Input : ! ! Input/output : ! ! Output : info (I) return status. 0 on success. ! ! Routines : iargc ! ! Remarks : ! ! References : ! ! Revisions : !------------------------------------------------------------------------- ! !------------------------------------------------------------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! CH-8092 Zurich, Switzerland !------------------------------------------------------------------------- SUBROUTINE pse_arguments(info) !------------------------------------------------------------------------- ! Modules !------------------------------------------------------------------------- USE pse_global IMPLICIT NONE !------------------------------------------------------------------------- ! Includes !------------------------------------------------------------------------- #ifdef __MPI INCLUDE 'mpif.h' #endif !------------------------------------------------------------------------- ! Arguments !------------------------------------------------------------------------- INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !------------------------------------------------------------------------- INTEGER :: i INTEGER :: info2 INTEGER :: nargc INTEGER :: ios INTEGER :: ilen,ilenprg CHARACTER(LEN = MAXCHAR) :: cbuf REAL(8) :: t0 !------------------------------------------------------------------------- ! Externals !------------------------------------------------------------------------- EXTERNAL iargc INTEGER iargc !------------------------------------------------------------------------- ! Initialise !------------------------------------------------------------------------- CALL substart('pse_arguments',t0,info) !------------------------------------------------------------------------- ! Root only !------------------------------------------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0,prgname) ilenprg = LEN_TRIM(prgname) !---------------------------------------------------------------------- ! Read input arguments !---------------------------------------------------------------------- #ifdef __MPI #ifdef __SUPER_UX nargc = iargc() i = 0 #elif defined __SunOS nargc = iargc() - 6 i = 0 #elif defined __HP_UX nargc = iargc() + 1 i = 1 #elif defined __IRIX64 nargc = iargc() - 4 i = 1 #elif defined __Linux !nargc = iargc() - 4 nargc = COMMAND_ARGUMENT_COUNT() i = 0 #elif defined __MacOS nargc = iargc() - 4 i = 0 #else !nargc = iargc() nargc = COMMAND_ARGUMENT_COUNT() i = 0 #endif #else !nargc = iargc() nargc = COMMAND_ARGUMENT_COUNT() i = 0 #endif DO WHILE (i.LT.nargc) i = i + 1 CALL getarg(i,cbuf) ilen = LEN_TRIM(cbuf) !------------------------------------------------------------------- ! Handle flags !------------------------------------------------------------------- IF (cbuf(1:1).EQ.'-') THEN !---------------------------------------------------------------- ! help !---------------------------------------------------------------- IF (cbuf(1:2).EQ.'-h') THEN WRITE(*,'(3A)')'Usage: ',prgname(1:ilenprg), & & ' [-h] [-v] [-d debuglevel] ' info = 1 GOTO 100 !---------------------------------------------------------------- ! Version !---------------------------------------------------------------- ELSEIF (cbuf(1:2).EQ.'-v') THEN WRITE(*,'(A)') 'This is ppm_pse, compiled with:' #ifdef __MPI WRITE(*,'(A)') ' -D__MPI Parallel version with MPI support' #endif #ifdef __NEWMPI WRITE(*,'(A)') ' -D__NEWMPI Support for MPI version 2.x' #endif #ifdef __SUPER_UX WRITE(*,'(A)') ' -D__SUPER_UX SUPPER UNIX fast array operations' #endif #ifdef __SunOS WRITE(*,'(A)') ' -D__SunOS SUN OS specific support' #endif #ifdef __IRIX64 WRITE(*,'(A)') ' -D__IRIX64 64bit IRIX support' #endif #ifdef __PGF90 WRITE(*,'(A)') ' -D__PGF90 Portland Group Fortran90 extensions' #endif #ifdef __XLF WRITE(*,'(A)') ' -D__XLF IBM XLF Fortran compiler support' #endif #ifdef __SXF90 WRITE(*,'(A)') ' -D__SXF90 NEC SX Fortran compiler support' #endif #ifdef __VECTOR WRITE(*,'(A)') ' -D__VECTOR Vector processor optimizations' #endif #ifdef __Linux WRITE(*,'(A)') ' -D__Linux Linux specific support' #endif #ifdef __MacOS WRITE(*,'(A)') ' -D__MacOS Mac OS X specific support' #endif #ifdef __ETIME WRITE(*,'(A)') ' -D__ETIME etime timing support' #endif info = 1 GOTO 100 !---------------------------------------------------------------- ! Debug flag !---------------------------------------------------------------- ELSEIF (cbuf(1:2).EQ.'-d') THEN !------------------------------------------------------------- ! Check that an argument has been passed !------------------------------------------------------------- i = i + 1 IF (i.GT.nargc) THEN WRITE(*,'(A)')'No argument passed' GOTO 9999 ENDIF !------------------------------------------------------------- ! Get the argument !------------------------------------------------------------- CALL getarg(i,cbuf) !------------------------------------------------------------- ! Convert the argument !------------------------------------------------------------- ilen = LEN_TRIM(cbuf) READ(cbuf,*,IOSTAT=ios) debug !------------------------------------------------------------- ! Check status of conversion !------------------------------------------------------------- IF (ios.NE.0) THEN WRITE(*,'(A)')'Invalid argument: ',cbuf(1:ilen) GOTO 9999 ENDIF !---------------------------------------------------------------- ! Unknown flag !---------------------------------------------------------------- ELSE WRITE(*,'(2A)')'Unknown flag: ',cbuf(1:ilen) GOTO 9999 ENDIF i = i + 1 CALL getarg(i,cbuf) ilen = LEN_TRIM(cbuf) i = i + 1 CALL getarg(i,cbuf) ilen = LEN_TRIM(cbuf) !------------------------------------------------------------------- ! Handle flags !--------------------------------------------------------or procD-------c hel# !------ A------------------------------------------------------ ! Convert the argument r !------ gument g (I) ret------------------ ! Convert the argument r !--------- WRITE(*,'(A)') ' ument,'(_BCIX6 ,1,,'(_ !--0, ----- IF 50 i + 1 ument,'(_BCIX6 ,1,,'(_ !--0, -----) ument,'(_BCIX6 cD------,! Exte,,'(_----------0, -----) ument,'(_BCIX6 ,! Exte,,'(_----------0, -----) ---------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0,',t0,info) RETURN END SUBROUTINE pse_sort_triangles --------------------rgc() Subroutine : pse_arguments 35561501424_comp.o: pse_mo------------ ! Subroutine : pse_arguments ustar omar11730204504omarLIABILITY NOTICE IN NO EVENT SHALL THE --------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0 ! Purpose : Read and pr1424_comp.o: pse_---------------------------------------- IF (rank .EQ. 0) THEN CALL getargL getarg(0 ! ! Input _coms surface .o: pse_ script WRITE(*,turn status. 0 on from ei ASCII .surf pse Open--- orRITE(*,turn status. 0 on fpse. Ine psetype utide(2lyHEN d from RITE(*,turn status. 0 on tne pseL getargL getarg(0Routines : pse (C) N ------- pseto _comL getargL getarg(0Routiemarks : ! ! References : ! ! Revisions : !-------------------------------------------L getarg(0,'---------------------L getarg(0,'------------------------------------------------------------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! C !------------------------------------------------------------------------- SUBROUTINE pse_sort_triangles(info) !------------------------------------------------------------------------- ! Modules !--------------_comp.o: pse_( pse ----------------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0 IMPLICIT NONE !------------------------------------------------------------------------- ! Argument---------------------------------------------------------- INTEGER ,INTENT( OUT) :: info !------------------------------------------------------------------------- ! Lo--------- =256------------------IN rnals pse rgument- :: a, b, x1, ------------------------------------------------------------------- INTEGER :: i, j, k, p, ix1, ix2 INTEGER :: minbiny INTEGER :: maxbiny, minb ! Lo--------- =80 !-----------------als psefm (I) rernals !------------------------------------------------ =256--------------------------- !------------------------------------ EXTERNAL iargc INTEGER iargc !------------------------------------------------------------------------- ! Initiali ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! C !------------------------------------------------------------------------- ! Root only !-------------_comp.o: pse_--------------------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0,nargulk)%l .o: pse_---------------------------------------- IF (rank .EQ. 0) THEN CALL getarALL geus ofrg(i,cb ',nar------rface script Ws pse im( pse ) ument----------- _comp.o: pse_--us of---------------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0((triang(2 pseforma (I) ret------------------ ! Convert the argument r !----------OPEN(20,FILE= pse -versUS='OLD'--AD(20,------- psefm IF 50( psefm ------'#ASCII ulk)%l --rface script W passed -------------------------------------------- #ifdef __MPI #ifdef __SUPER_UX narg.surf orma (I) re -------------------------------------------- #ifdef __MPI #ifdef __SUPER_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX AD(20,*) R_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX AD(20,------R_UX ALL geus ofrgA,I6,i,cb 'S-rface ulk) Ws ' -----------------'---------- ument----------- _comp.o: pse_--us of------ R_UX f ! lo,j,i): x,y,z; sed: number;,k,ird: e this linumberR_UX menOCA ge. NBINZ+13 ),versioi GOTO 9999 ENDIF ument----------- _comp.o: pse_--'Not enough bin(i,'= i NBINZ,MK) = i + 1 A------- = 0 =AD(20,*) . NBINZ: ,. NBINZ: 2 ,. NBINZ: 3(1)*y1(2)) sts to actual SCLOSE(20)*y1(2))---------------------------------------------------------- #ifdef __MPI #ifdef __SUPER_UX f) orma (I) re -------------------------------------------- #ifdef __MPI #ifdef __SUPER_UX ument----------- _comp.o: pse_--'f) pseforma !'= i NBINZ,MK) !----------------------- uLOSE(20)*y1(2))ument----------- _comp.o: pse_--' !---rnar------rface .o: pse_---------------------------------------------------- 9999 CONTINUE CALL substop('pse_sort_triangles',t0,info) RETURN END SUBROUTINE pse_sort_triangles _comp.o: pse_----------- =_comp.o: pse_- 35561501424_hustar omar452041201302omarPPM PSE code ============ This is a strrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr------------------------------------- 9999 CONTINUE CALL substop('pse_sort_triangleFun : Read and procepr1424_hin lists by assigning the ! triangles to bins of equal y and z coordinates. ! ! In--righ nd ---------ENTSA -----RITE(*,turn status. 0 on IF ac 2) d!---THENL getargL getarg(0Routines : xpl(:----(F) ticse RITE(*,turn status. 0 on wpl(:----(F) ticsestrengthRITE(*,turn status. 0 on lda 0 on lda ---wpRITE(*,turn status. 0 on Np 0 on l ::number --- ticsesL getargL getarg(0Routiemarks : ! ! References dwpl(:,:)(F) C ng ticsestrength duetoRITE(*,turn status. 0 on IF ac.RITE(*,turn status. 0 on ! ! Revisions : !--------------------------------!-------RITE(*,turn status. 0 on RITE(*,turn status. 0 on 3556def _pp_r(3)/) turn status. 0 on 3556def _pp_verle-RITE(*,turn status. 0 on 3556def _pp_minb--------L getarg(0,'-------------lda nd Np 2) need(NBINYformSUPER_ ! ! Output 3556 !-z.chrequire ipack, rpack, lpack ! ! Output 2) (currently unused) dummyt : iNOT. Ae ODERITE(*,turn status. 0 on v1(1r !------------------'------------------------------------------------------------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini --------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini C !------------------------les(info) ---------- !------------------------------------------------------------------------- ! Modules FUNCTIONr1424_hi(o_id,xpl,wpl,dwpl,lda,Np,ipack,lpack,rpack,---------------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0 IMP--------------------------------------------------------------- ! Includes !----556 !-user !----556 !-user_nume s !----556 !-da,2) s !-io) s ! Argument--------------s !-------------"------"---------------------------------------------------- INTEGER, INTENT( OUT) :: info !------------------------------------------------------------------------- ! Local variables !-urn status. 0 on ---------INon ----------------o_id (I) rernals---------------:,:), PO-----urn status. 0 on :: xpl (I) rernals---------------:,:), PO-----urn status. 0 on :: wpl (I) rernals---------------:,:), PO-----urn status. 0 on :: dwpl (I) re !-urn status. 0 on ---------INon ---------------lda,Np (I) re !-urn status. 0 on ---------------------------------------- :: maxbiny, ------------1424_hin lists ------------:,:), -------INon -, OPTIONAL-----pack ! ------------:,:), -------INon -, OPTIONAL----lpack ! rnals---------------:,:), -------INon -, OPTIONAL----rpack ! ---------------------------------------------------- ! Local variables :: i, j, k, p, ix1, ix2 INTEGER :: minbiny INTEGER rnals !------------------------------------ :: maxbin----p,iinfo,iiymm ! ---------------------------------------------------- ! Local variables !------------------------------------------------------------------------- ! Initiali ! Institute of Computational Science ! ETH Zentrum, Hirschengraben 84 ! C !----------------------------------------------------------------rum, Hirschengraben 84 ! !-------------_h----------------_h-g(2:3,3,i)-iiymmg(2:3,3,i)- iymmpse_)-iiymmg(21-------------------------------------------------------- CALL pse_alloc_binlist(0,0,0,info) ze dwp---------------------------------------------------------------------- ! Initiali ! I= 0 ps',' WRITE(dwpl(1,ip DO ix2=3,4 e the triangl------------------------------------------------ CALL pse_alloc_binlist(0,0,0,U strength ---ghost- ticses.triangl------------------------------------------------ CALL pse_alloc_binlist(0,0an arutoff0) GOTiangle: ENDIF iymmpse_)- !---------------!hreside(1,Ae mapp(3), : iymmgis used, stil det2 = v1! current det2 = v1ument-556 a _p t_ghost_get(o_id,xpl,---m,' t,iiymm,-----------------------rutoff= i GOTO 9999 ENDIF 2 = v1ument----------- 55615'------------------------------REAL(NBINYget-ghost---------------------------- i + 1 CALL etarg(i,cbuf) -------------= v1ument-556 a _p t_push(wpl,NSPEC,' t,i--------------- GOTO 9999 ENDIF 2 ument----------- 55615'--------------------------REAL(NBINYpush-ghost-strength----------------------- i + 1 C-------------= v1ument-556 a _p t_s---(' t,M t,i--------------- GOTO 9999 ENDIF 2ument----------- 55615'--REAL(NBINYs----ghost----------------------- i + 1 C-------------= v1ument-556 a _p t_pop(wpl,NSPEC,' t,M t,i--------------- GOTO 9999 ENDIF 2 ument----------- 55615'--------------------------REAL(NBINYpop-ghost-strength----------------------- i + 1 C-------------= v1 iymmpse_)- !---------------!hINYempty ufferR_UX -= v1ument-556 a _p t_pop(xpl,---m,' t,M t,i------------------- GOTO 9999 ENDIF 2 2 ument----------- 55615'------------------------------REAL(NBINYpop-ghost-posi --------------------------- i + 1 CALL etarg(i,cbuf) ------ C------- triangl------------------------------------------------ CALL pse_alloc_binlist(0,0,0,PP IF acINYevalu Laplaciantriangl------------------------------------------------ CALL pse_alloc_binlist(0,0an arutoff0)LGOTiangle,i) neigho) = (/bb(1), bb(4)/) ument-556def _pp_r(3)(xpl,Np,wpl,lda,iymmpse_,-556k el--------------------k ,dwpl,--------------------------o iymmpse_)- !---------------!hU ghost- as wel det2 = v1tIn =eigho) = (/bb(1), bb(3)/) tIn op oup ENDIF ument-556def _pp_verle-(xpl,M t,wpl,lda,iymmpse_,----------------------------ktab,drtj, kinv,nv-----v-----dwpl,------------------------t aument-556def _pp_verle-(xpl,M t,wpl,lda,iymmpse_,-----------------------------556k el-k ,nv-----v-----dwpl,------------------------t-------------------hat an a=eigho) = (/b ENDIF tIn op oup ENDIF ument-556def _pp_minb(o_id,xpl,M t,wpl,lda,iymmpse_,-----------------------------ktab,drtj, kinv,Nm,c-----rutoff2-dwpl,------------------------t aument-556def _pp_minb(o_id,xpl,M t,wpl,lda,iymmpse_,-----------------------------556k el-k ,Nm,c-----rutoff2-dwpl,------------------------t-------------------htarg(i,cbuf) - ! O - ticsesL getarrrrrrrrrtIn =eigho) = (/bb(1), bb(3)/) tIn op oup ENDIF ument-556def _pp_verle-(xpl,Np,wpl,lda,iymmpse_,----------------------------ktab,drtj, kinv,nv-----v-----dwpl,------------------------t aument-556def _pp_verle-(xpl,Np,wpl,lda,iymmpse_,-----------------------------556k el-k ,nv-----v-----dwpl,------------------------t-------------------hat an a=eigho) = (/b ENDIF tIn op oup ENDIF ument-556def _pp_minb(o_id,xpl,Np,wpl,lda,iymmpse_,----------------------------ktab,drtj, kinv,Nm,c-----rutoff2-dwpl,------------------------t aument-556def _pp_minb(o_id,xpl,Np,wpl,lda,iymmpse_,-----------------------------556k el-k ,Nm,c-----rutoff2-dwpl,------------------------t-------------------htarg(i,cbuf) ------ C------ C GOTO 9999 ENDIF ument----------- 55615'--PP IF acfEAL(N------------------- !-----------------triangl------------------------------------------------ CALL pse_alloc_binlist(0,0,0,S----ghostsribuba f.EQ. 2) us-----ymmpse_triangl------------------------------------------------ CALL pse_alloc_binlist(0,0an a-ymmpse_= nbin(rutoff0) GOTiangle: ENDIF ! W) d!- !need xp = 0.0_MK a stde(d IF n(2lyH----------F ! MK posi 2) no bQ. ggrega sour.0_ __M-------F ! anyway...-------F ! Inputmesses IF n(2 D, so DO -----EQ.w toRITE(*,F ! get-ghost--EQ.h INYdo a ghost_get. gaangles edges ument-556 a _p t_ghost_put(o_id,i--------------- GOTO 9999 ENDIF 2 ument----------- 55615'--------------------------REAL(NBINYput-ghost------------------------ i + 1 C-------------= v1ument-556 a _p t_push(dwpl,NSPEC,' t,i--------------- GOTO 9999 ENDIF 2ument----------- 55615'--REAL(NBINYpush-ghost----------------------- i + 1 C-------------= v1ument-556 a _p t_s---(' t,M t,i--------------- GOTO 9999 ENDIF 2ument----------- 55615'--REAL(NBINYs----ghost----------------------- i + 1 C-------------= v1ument-556 a _p t_ghost_pop(dwpl,NSPEC,' t,M t,i--------------- GOTO 9999 ENDIF 2 ument----------- 55615'--------------------------REAL(NBINYpop-ghost-strength----------------------- i + 1 C------ C---------------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0 ultiplyHUPER_ diffus Wssta (I) ret------------------ ! Convert the argument r !---------- lda = (/bb(1), bb(3)/) = 0 ps',' WRITE( E(dwpl(1,ip dwpl(1,ip *diff ! on wha-------------at an alda = (/b ENDIF= 0 ps',' WRITE( E(dwpl(1,ip dwpl(1,ip *diff ! on wha( E(dwpl(2,ip dwpl(2,ip *diff 2 ! on wha-------------at an alda = (/bb(1), bb(4)/) = 0 ps',' WRITE( E(dwpl(1,ip dwpl(1,ip *diff ! on wha( E(dwpl(2,ip dwpl(2,ip *diff 2 ! on wha( E(dwpl(3,ip dwpl(3,ip *diff 3 ! on wha-------------at b(4)/) = 0 infos',lda ! on wha( E(= 0 ps',' WRITE( E(( E(dwpl(iinfo,ip dwpl(iinfo,ip *diff iinfo ! on wha( E(------------------------------------------------------------ ! Convert the argument r !------ guments',t0,info) RETURN END SUBROUTINE pse_sort_triangles _h------------FUNCTIONr1424_hippm_pse/pse_arguments.f00006440001750000175000000021714114555555555555535561501556nume s.h000064400017500001750000000217141145555555555555ustar omar072334541620014057omarPPM PSE code ============ This is a strrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr------------------------- ! Convert the argument r !------ gument IM------- : Read and procep155 nume s header pse-------------------------- ! Convert the argument r !------ gumetriangl------------------------------------------------ CALL pse_alloc_binlist(0,0,0,Pendif #iP ticseMesh-Librae_=(PPM)ist(0,0,0,--------------------------------------------------- SUBROUTINE pse_arguments(info) !------------------------------------------------------------------------- ! Modules !------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0(( quadra,t0e r IMP-NOT. Ae unit e this li(BEM)ist(0,0,-------------------------------- IF (rank .EQ. 0) THEN CALL getar !---PARAME---u---15561endm_bem_quadr !-ce IF !---PARAME---u---15561endm_bem_quadr !-no---- 2 !---PARAME---u---15561endm_bem_quadr !- - 3L getar !---PARAME---u---15561endm_bem_quadr !-cpose 4L getar !---PARAME---u---15561endm_bem_quadr !-stroud 5L getar !---PARAME---u---15561endm_bem_quadr !-hamme 3= - 4 !---PARAME---u---15561endm_bem_quadr !-hamme 4= 7 4 !---PARAME---u---15561endm_bem_quadr !-hamme 7= 8 4 !---PARAME---u---15561endm_bem_quadr !-hamme 12= 9 !------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0(( basiP-Nun P-NOT. Ae unit e this li(BEM)ist(0,0,-------------------------------- IF (rank .EQ. 0) THEN CALL getar !---PARAME---u---15561endm_bem_basiP_st se !---PARAME---u---15561endm_bem_basiP_utpu- 2 !---PARAME---u---15561endm_bem_basiP_quad se 3 !------------------------------------- IF (rank .EQ. 0) THEN CALL getarg(0(( field v1(1rin lists by assigning the ! triangles to bins of equal y and z coor !---PARAME---u---15561endm_v1(1r_mg (2:3,3,i)- !---PARAME---u---15561endm_v1(1r_OS (2 !---PARAME---u---15561endm_v1(1r_fft (22 n lists by assigning the ! triangles to bins of equal y anL getarg(0(( equa 55 v1(1ri n lists by assigning the ! triangles to bins of equal y anL getar !---PARAME---u---15561endm_eq_poiss (2 !------------------------------------- IF (rank .EQ. 0) THEN CAL getarg(0(( order or pnit diff------ !------------------------------------- IF (rank .EQ. 0) THEN CA L getar !---PARAME---u---15561endm_order_1 (2 !---PARAME---u---15561endm_order_2 (22 !---PARAME---u---15561endm_order_3 (23L getar !---PARAME---u---15561endm_order_4 (24 !------------------------------------- IF (rank .EQ. 0) THEN CAL getarg(0(( smooSwitzerland -------------------------------- IF (rank .EQ. 0) THEN CAL getar !---PARAME---u---15561endm_mooSw_rbOS (2 !------------------------------------- IF (rank .EQ. 0) THEN CAAL getarg(0ODE 1endm(triin lists by assigning the ! triangles to bins of equal y and z coor !---PARAME---u---15561endm_o--_UBROm!- !rf (2 !---PARAME---u---15561endm_o--_UBROm!-tvdrk2 (22 !---PARAME---u---15561endm_o--_UBROm!-midrk2 (23L getar !---PARAME---u---15561endm_o--_UBROm!-rk4 (24 !---PARAME---u---15561endm_o--_UBROm!-trapez (25L getar !---PARAME---u---15561endm_o--_UBROm!-tvdrk3 (2- 4 !---PARAME---u---15561endm_o--_UBROm!-s------- = 7 n============ This is a strrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr35561501424 !-ioustar omar0577137751214omarPPM PSE code ============ This is a strrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr !------------------------------------- IF (rank .EQ. 0) THEN CAAL getarg(0 IM------- : Read and proc !-io) --------------------------------- IF (rank .EQ. 0) THEN CALL getargL getarg(0 ! ! Input !s allr--------need(NB or L getarg(0000000000000000I/O.RITE(*,turn status. 0 on--- guments'---------------------L getarg(0,'------------------------------------------------------------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini $Log: 1424 !-io,v $------------------ 1.1.1.1 2007/07/13 10:18:16 ivoRITE(*,turCBL __SunOS --------PM librae_-------------------------- 1.3 2005/04/27 01:21:45 ivoRITE(*,tur-------l .o: pse_f __MacO----directly into client. Not r-----------pre __MacOS i--need(N. Input __gndm now-directly _coms Open--- orRITE(*,tur psesTSA so, c ng d -----stepper INYp556 !-z.c !------------------'------ 1.2 2004/06/03 16:11:25 ivoRITE(*,turAd--l 1424-----_diag !------------------'------ 1.1 2004/06/02 15:28:33 ivoRITE(*,tur--fo) implea.RITE(*,t----------------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini !------------------------------------------------------------------------- SUBROUTINE pse_arguments(info) !------------------------------------------------------------------------- ! Modules !----MODUL !-io) ------------------------------------------------------- ! Modules ) ----(0(( IF face o_compcD--) ------------------------------------------------------- ! Modules ) --------FACE_compcD--) -------MODUL PROCEDURE_compcD--) --------FACE ) ------------------------------------------------------- ! Modules ) ----(0(( IF face o_comp.o: pse_-------------------------------------------------------------- ! Modules ) --------FACE_comp.o: pse_--------------MODUL PROCEDURE_comp.o: pse_---------------FACE ) ------------------------------------------------------- ! Modules ) ----(0(( IF face o-----_marks -------------------------------------------------------------- ! Modules ) --------FACE-----_marks --------------MODUL PROCEDURE-----_marks ---------------FACE ------------------------------------------------------- ! Modules ) ----(0(( IF face o-----_diag-------------------------------------------------------------- ! Modules ) --------FACE-----_diag--------------MODUL PROCEDURE-----_diag---------------FACE ------------------------------------------------------- ! Modules ) ----(0------- sour.0 ------------------------------------------------------- ! Modules ) ---AINS #i------ "_compcD--.f" #i------ "_comp.o: pse_mo" #i------ "-----_marks mo" #i------ "-----_diag o" MODUL !-io)rrrr35561501424cEN _abortustar omar077454541620014477omarPPM PSE code ============ This is a strrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr ------------------------------------------------------- ! Modules -zarini ! Purpose : Read and pr 1424cEN _abort) --------------------------------- IF (rank .EQ. 0) THEN CALL getargL getarg(0 ! ! R----0 cEN s f.a- ticular pses ex D WRITE(*,turn status. 0 on l ::directoe_m If so, all __MacOS s haltL getargL getarg(0Routines :L getargL getarg(0Routiemarks : ! ! References lEx ences(L) .TRUE. f.---- pseex D, .FALSE.RITE(*,turn status. 0 on oSwitwis ! ! Output : ! ! Revisions : !----------------------------------!------------------'-------------Could cEN -----ity---- pse (length etc). Ine ! ! Output Abort flict--EPER_ Ct :: p L getarg(0000000000000000systemD, thu--EQ.re it eo ChkAbort !------------------'------------------------------------------------------------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini $Log: 1424cEN _abort,v $------------------ 1.1.1.1 2007/07/13 10:18:16 ivoRITE(*,turCBL __SunOS --------PM librae_-------------------------- 1.1 2004/06/02 15:28:27 ivoRITE(*,tur--fo) implea.RITE(*,t----------------------------- ! Jens Honore Walther ! Ivo F. Sbalzarini !------------------------------------------------------------------------- SUBROUTINE pse_arguments(info)