! ============================================================= ! TREEGRAV: routines to perform hierarchical force calculation. ! ============================================================= ! -------------------------------------------------------- ! ACCEL: routine to compute the acceleration and potential ! for all bodies in the system. ! -------------------------------------------------------- SUBROUTINE ACCEL INCLUDE 'treedefs.f90' DOUBLE PRECISION CPUT1, CPUT2 INTEGER P INTEGER :: NC_temp,NB_temp,NT_temp ! ---------------------------------------- ! Initialize timing for force calculation. ! ---------------------------------------- ! CALL SECOND(CPUT1) ! CALL SECNDS(CPUT1) CPUT1=MPI_WTIME() ! -------------------------------------------- ! Initialize the force evaluation diagnostics. ! -------------------------------------------- NBTOT = 0 NCTOT = 0 NTMAX = 0 ! ---------------------------------------------------------------- ! Distribute body data to nodes ! ---------------------------------------------------------------- CALL SEND_RAW_DATA !MPI_INDEX comes from this subroutine ! --------------------------------------------------- ! Loop over bodies, computing force for each in turn. ! --------------------------------------------------- DO P = mpi_index, NBODY CALL TRWALK(P) END DO ! ---------------------------------------------------------------- ! Recieve Processed Data and Diagnostics from Nodes ! ---------------------------------------------------------------- CALL RECV_PROCESSED_DATA ! ----------------------------------------------- ! Compute average number of force terms per body. ! ----------------------------------------------- NTAVG = (NBTOT + NCTOT) / NBODY ! ------------------------------------- ! Compute timing for force calculation. ! ------------------------------------- ! CALL SECOND(CPUT2) ! CALL SECNDS(CPUT2) CPUT2=MPI_WTIME() CPUACC = CPUT2 - CPUT1 END ! -------------------------------------------------------------- ! TRWALK: recursive routine to walk the tree computing forces on ! particle P. ! -------------------------------------------------------------- SUBROUTINE TRWALK(P) INTEGER P INCLUDE 'treedefs.f90' INTEGER MXSPTR PARAMETER(MXSPTR = 256) DOUBLE PRECISION :: EPS2, PHI0, ACC0(NDIM), POS0(NDIM) DOUBLE PRECISION :: DX, DY, DZ, DR2, DR2INV, DRINV, PHIM, DR5INV, PHIQ DOUBLE PRECISION :: N_sq,ACC_Mond,AX,AY,AZ,A_not,Newton,Numerator,Denominator INTEGER Q, NBTERM, NCTERM, SPTR, STACK(MXSPTR), K LOGICAL SKPSLF EPS2 = EPS*EPS ! ---------------------------------------------------------- ! Zero potential and acceleration for subsequent summations. ! ---------------------------------------------------------- PHI0 = 0.0 ACC0(1) = 0.0 ACC0(2) = 0.0 ACC0(3) = 0.0 ! ----------------------------------------- ! Copy position of P for quicker reference. ! ----------------------------------------- POS0(1) = POS(P,1) POS0(2) = POS(P,2) POS0(3) = POS(P,3) ! ---------------------------------------------- ! Initialize counters for number of force terms. ! ---------------------------------------------- NBTERM = 0 NCTERM = 0 SKPSLF = .FALSE. ! ---------------------------------- ! Push the root cell onto the stack. ! ---------------------------------- SPTR = 1 STACK(SPTR) = ROOT ! ------------------------------------- ! Loop while nodes on stack to process. ! ------------------------------------- DO WHILE (SPTR .GT. 0) ! -------------------------- ! Pop node off top of stack. ! -------------------------- Q = STACK(SPTR) SPTR = SPTR - 1 ! --------------------------------------------- ! Compute distance to center-of-mass of node Q. ! --------------------------------------------- DX = POS0(1) - POS(Q,1) DY = POS0(2) - POS(Q,2) DZ = POS0(3) - POS(Q,3) DR2 = DX*DX + DY*DY + DZ*DZ ! ------------------------------- ! Classify Q as a body or a cell. ! ------------------------------- IF (Q .LT. INCELL) THEN ! ----------------------------------- ! A body: check for self-interaction. ! ----------------------------------- IF (Q .NE. P) THEN ! ------------------------------ ! Compute body-body interaction. ! ------------------------------ DR2INV = 1.0 / (DR2 + EPS2) PHIM = MASS(Q) * SQRT(DR2INV) PHI0 = PHI0 - PHIM PHIM = PHIM * DR2INV ACC0(1) = ACC0(1) - PHIM * DX ACC0(2) = ACC0(2) - PHIM * DY ACC0(3) = ACC0(3) - PHIM * DZ NBTERM = NBTERM + 1 ELSE ! ------------------------------------------- ! Remember that self-interaction was skipped. ! ------------------------------------------- SKPSLF = .TRUE. ENDIF ELSE ! -------------------------------------------- ! A cell: test if interaction can be accepted. ! -------------------------------------------- IF (DR2 .GE. RCRIT2(Q)) THEN ! ---------------------------------------- ! Accepted: compute body-cell interaction. ! ---------------------------------------- DR2INV = 1.0 / (DR2 + EPS2) DRINV = SQRT(DR2INV) PHIM = MASS(Q) * DRINV PHI0 = PHI0 - PHIM PHIM = PHIM * DR2INV ACC0(1) = ACC0(1) - PHIM * DX ACC0(2) = ACC0(2) - PHIM * DY ACC0(3) = ACC0(3) - PHIM * DZ ! ------------------------------------ ! Optionally include quadrupole terms. ! ------------------------------------ IF (USQUAD) THEN DR5INV = DR2INV * DR2INV * DRINV PHIQ = DR5INV * & (0.5 * ((DX*DX - DZ*DZ) * QUAD(Q,1) + & (DY*DY - DZ*DZ) * QUAD(Q,4)) + & DX*DY * QUAD(Q,2) + DX*DZ * QUAD(Q,3) + & DY*DZ * QUAD(Q,5)) PHI0 = PHI0 - PHIQ PHIQ = 5.0 * PHIQ * DR2INV ACC0(1) = ACC0(1) - PHIQ*DX + DR5INV * & (DX*QUAD(Q,1) + DY*QUAD(Q,2) + DZ*QUAD(Q,3)) ACC0(2) = ACC0(2) - PHIQ*DY + DR5INV * & (DX*QUAD(Q,2) + DY*QUAD(Q,4) + DZ*QUAD(Q,5)) ACC0(3) = ACC0(3) - PHIQ*DZ + DR5INV * & (DX*QUAD(Q,3) + DY*QUAD(Q,5) - & DZ*(QUAD(Q,1) + QUAD(Q,4))) ENDIF NCTERM = NCTERM + 1 ELSE ! ----------------------------------- ! Rejected: examine children of cell. ! ----------------------------------- DO K = 1, NSUBC ! -------------------------------------- ! Push existing children onto the stack. ! -------------------------------------- IF (SUBP(Q,K) .NE. NULL) THEN IF (SPTR .GE. MXSPTR) & CALL TERROR(' TRWALK: STACK OVERFLOW') SPTR = SPTR + 1 STACK(SPTR) = SUBP(Q,K) ENDIF END DO ENDIF ENDIF END DO ! --------------------------------------------------- ! Check that self-interaction was explicitly skipped. ! --------------------------------------------------- IF (.NOT. SKPSLF) & CALL TERROR(' TRWALK: MISSED SELF-INTERACTION') ! ---------------------------------------------- ! Calculate Mond Accelerations ! ---------------------------------------------- ! N_sq = ACC0(1)*ACCO(1) + ACC0(2)*ACC0(2) + ACC0(3)*ACC0(3) ! Newton= SQRT(N_sq) ! ! AX = ACC0(1) / Newton ! AY = ACC0(2) / Newton ! AZ = ACC0(3) / Newton ! ! ACC_MOND = Newton * (0.5 + 0.5 * SQRT(1.0 + 4.0 * a_not**2 / (N_sq)))**0.5 ! ! ACCO(1) = AX * ACC_MOND ! ACC0(2) = AY * ACC_MOND ! ACC0(3) = AZ * ACC_MOND ! ---------------------------------------------- ! Copy total potential and acceleration to body. ! ---------------------------------------------- PHI(P) = PHI0 ACC(P,1) = ACC0(1) ACC(P,2) = ACC0(2) ACC(P,3) = ACC0(3) ! ------------------------------------- ! Update force calculation diagnostics. ! ------------------------------------- NBTOT = NBTOT + NBTERM NCTOT = NCTOT + NCTERM NTMAX = MAX(NTMAX, NCTERM + NBTERM) END