! ================================================================ ! TREELOAD: routines to handle tree construction and verification. ! ================================================================ ! ---------------------------------------------------------------- ! MKTREE: initialize the tree structure for the force calculation. ! ---------------------------------------------------------------- SUBROUTINE MKTREE INCLUDE 'treedefs.f90' DOUBLE PRECISION CPUT1, CPUT2 ! ---------------------------------- ! Start timing of tree construction. ! ---------------------------------- ! CALL SECOND(CPUT1) ! CALL SECNDS(CPUT1) CPUT1=MPI_WTIME() ! -------------------------------------- ! Expand root volume to hold all bodies. ! -------------------------------------- CALL EXPBOX ! -------------------------- ! Load bodies into the tree. ! -------------------------- CALL LDTREE ! ------------------------------------------------ ! Compute masses, center of mass coordinates, etc. ! ------------------------------------------------ CALL HACKCM ! ------------------------------------- ! Compute timing for tree construction. ! ------------------------------------- ! CALL SECOND(CPUT2) ! CALL SECNDS(CPUT2) CPUT2=MPI_WTIME() CPUMKT = CPUT2 - CPUT1 END ! ----------------------------------------------------- ! EXPBOX: enlarge system cube to include all particles. ! ----------------------------------------------------- SUBROUTINE EXPBOX INTRINSIC MAX INCLUDE 'treedefs.f90' DOUBLE PRECISION :: XYZMAX INTEGER J, I ! ------------------------------ ! Find maximum coordinate value. ! ------------------------------ XYZMAX = 0 DO I = 1, NBODY DO J = 1, NDIM XYZMAX = MAX(XYZMAX, ABS(POS(I,J))) END DO END DO ! ----------------------------------------- ! Expand box by factors of 2 until it fits. ! ----------------------------------------- DO WHILE (XYZMAX .GE. RSIZE/2.0) RSIZE = 2.0 * RSIZE END DO END ! --------------------------------------------------------------------- ! LDTREE: construct tree body-by-body. This phase initializes the SUBP ! array and loads the geometric midpoint into the MID of each cell. ! --------------------------------------------------------------------- SUBROUTINE LDTREE INCLUDE 'treedefs.f90' INTEGER MKCELL, K, P ! --------------------------------------- ! Deallocate current tree, begin new one. ! --------------------------------------- NCELL = 0 ROOT = MKCELL() ! ------------------------------------------ ! Initialize midpoint and size of root cell. ! ------------------------------------------ DO K = 1, NDIM MID(ROOT,K) = 0.0 END DO CLSIZE(ROOT) = RSIZE ! --------------------------------------------- ! Load bodies into the new tree, one at a time. ! --------------------------------------------- DO P = 1, NBODY CALL LDBODY(P) END DO END ! ---------------------------------------- ! LDBODY: load body P into tree structure. ! ---------------------------------------- SUBROUTINE LDBODY(P) INTEGER P INCLUDE 'treedefs.f90' INTEGER Q, QIND, SBINDX, MKCELL, C, K, P0 ! --------------------------------------------- ! Start Q,QIND pair in correct subcell of root. ! --------------------------------------------- Q = ROOT QIND = SBINDX(P, Q) ! ----------------------------------------------------- ! Loop descending tree until an empty subcell is found. ! ----------------------------------------------------- DO WHILE (SUBP(Q, QIND) .NE. NULL) ! -------------------------------------- ! On reaching another body, extend tree. ! -------------------------------------- IF (SUBP(Q, QIND) .LT. INCELL) THEN ! ------------------------------------------- ! Allocate an empty cell to hold both bodies. ! ------------------------------------------- C = MKCELL() ! ------------------------------------------------------ ! Locate midpoint of new cell wrt. parent, and set size. ! ------------------------------------------------------ DO K = 1, NDIM IF (POS(P,K) .GE. MID(Q,K)) THEN MID(C,K) = MID(Q,K) + CLSIZE(Q)/4.0 ELSE MID(C,K) = MID(Q,K) - CLSIZE(Q)/4.0 ENDIF END DO CLSIZE(C) = CLSIZE(Q) / 2.0 ! ------------------------------------------------------ ! Store old body in appropriate subcell within new cell. ! ------------------------------------------------------ P0 = SUBP(Q, QIND) SUBP(C, SBINDX(P0, C)) = P0 ! --------------------------------------------- ! Link new cell into tree in place of old body. ! --------------------------------------------- SUBP(Q, QIND) = C ENDIF ! -------------------------------------------------------- ! At this point, the node indexed by Q,QIND is known to be ! a cell, so advance to the next level of tree, and loop. ! -------------------------------------------------------- Q = SUBP(Q, QIND) QIND = SBINDX(P, Q) END DO ! --------------------------------------------- ! Found place in tree for P, so store it there. ! --------------------------------------------- SUBP(Q, QIND) = P END ! ------------------------------------------------------- ! SBINDX: compute subcell index for node P within cell Q. ! ------------------------------------------------------- INTEGER FUNCTION SBINDX(P, Q) INTEGER P, Q INCLUDE 'treedefs.f90' INTEGER K ! --------------------------------------------------- ! Initialize subindex to point to lower left subcell. ! --------------------------------------------------- SBINDX = 1 ! --------------------------------- ! Loop over all spatial dimensions. ! --------------------------------- DO K = 1, NDIM IF (POS(P,K) .GE. MID(Q,K)) & SBINDX = SBINDX + 2 ** (3 - K) END DO END ! --------------------------------------------------------- ! MKCELL: function to allocate a cell, returning its index. ! --------------------------------------------------------- INTEGER FUNCTION MKCELL() INCLUDE 'treedefs.f90' INTEGER I ! ---------------------------------------------------------- ! Terminate simulation if no remaining space for a new cell. ! ---------------------------------------------------------- IF (NCELL .GE. MXCELL) & CALL TERROR(' MKCELL: NO MORE MEMORY') ! ---------------------------------------------------- ! Increment cell counter, initialize new cell pointer. ! ---------------------------------------------------- NCELL = NCELL + 1 MKCELL = NCELL + MXBODY ! -------------------------------------- ! Zero pointers to subcells of new cell. ! -------------------------------------- DO I = 1, NSUBC SUBP(MKCELL,I) = NULL END DO END ! ---------------------------------------------------------------- ! HACKCM: compute cell masses, c.m. positions, check tree structure, ! assign critical radii, and compute quadrupole moments. Whew! ! ---------------------------------------------------------------- SUBROUTINE HACKCM INCLUDE 'treedefs.f90' INTEGER IND(MXCELL), P, Q, I, J, K, L, M, N DOUBLE PRECISION :: POS0(NDIM), DIST2 ! --------------------------------------- ! List cells in order of decreasing size. ! --------------------------------------- CALL BFLIST(IND) ! -------------------------------------------- ! Loop processing cells from smallest to root. ! -------------------------------------------- DO I = NCELL, 1, -1 P = IND(I) ! -------------------------------------------------------------- ! Zero accumulators for this cell. A temporary variable is used ! for the c. of m. so as to preserve the stored midpoints. ! -------------------------------------------------------------- MASS(P) = 0.0 DO K = 1, NDIM POS0(K) = 0.0 END DO ! ------------------------------------------------------------- ! Compute cell properties as sum of properties of its subcells. ! ------------------------------------------------------------- DO J = 1, NSUBC Q = SUBP(P,J) ! ------------------------------ ! Only access cells which exist. ! ------------------------------ IF (Q .NE. NULL) THEN ! ------------------------------------------------------- ! Sum properties of subcells to obtain values for cell P. ! ------------------------------------------------------- MASS(P) = MASS(P) + MASS(Q) DO K = 1, NDIM POS0(K) = POS0(K) + MASS(Q) * POS(Q,K) END DO ENDIF END DO ! -------------------------------------------------------- ! Normalize center of mass coordinates by total cell mass. ! -------------------------------------------------------- DO K = 1, NDIM POS0(K) = POS0(K) / MASS(P) END DO ! ----------------------------------------------------------------- ! Check tree, compute cm-to-mid distance, and assign cell position. ! ----------------------------------------------------------------- DIST2 = 0.0 DO K = 1, NDIM IF (POS0(K) .LT. MID(P,K) - CLSIZE(P)/2.0 .OR. & POS0(K) .GE. MID(P,K) + CLSIZE(P)/2.0) THEN WRITE(ULOG, '(/,1X,''TREE ERROR'',2I6,3E14.6)') & P, K, POS0(K), MID(P,K), CLSIZE(P) CALL TERROR(' HACKCM: TREE STRUCTURE ERROR') ENDIF DIST2 = DIST2 + (POS0(K) - MID(P,K))**2 ! -------------------------------------------------------- ! Copy cm position to cell. This overwrites the midpoint. ! -------------------------------------------------------- POS(P,K) = POS0(K) END DO ! ------------------------------------------------------------ ! Assign critical radius for cell, adding offset from midpoint ! for more accurate forces. This overwrites the cell size. ! ------------------------------------------------------------ RCRIT2(P) = (CLSIZE(P) / THETA + SQRT(DIST2))**2 END DO ! ---------------------------------------- ! Compute quadrupole moments, if required. ! ---------------------------------------- IF (USQUAD) THEN ! ------------------------------- ! Loop processing cells as above. ! ------------------------------- DO I = NCELL, 1, -1 P = IND(I) ! -------------------------------------------- ! Zero accumulator for quad moments of cell P. ! -------------------------------------------- DO K = 1, NQUAD QUAD(P,K) = 0.0 END DO ! -------------------------------- ! Loop over descendents of cell P. ! -------------------------------- DO J = 1, NSUBC Q = SUBP(P,J) IF (Q .NE. NULL) THEN ! -------------------------------------------------------- ! Sum properties of subcell Q to obtain values for cell P. ! -------------------------------------------------------- DO M = 1, MIN(2,NDIM) DO N = M, NDIM L = (M-1) * (NDIM-1) + N QUAD(P,L) = QUAD(P,L) + 3.0 * MASS(Q) * & (POS(Q,M) - POS(P,M)) * (POS(Q,N) - POS(P,N)) IF (M .EQ. N) THEN DO K = 1, NDIM QUAD(P,L) = QUAD(P,L) - MASS(Q) * & (POS(Q,K) - POS(P,K))**2 END DO ENDIF ! ------------------------------------------- ! If Q itself is a cell, add its moments too. ! ------------------------------------------- IF (Q .GE. INCELL) & QUAD(P,L) = QUAD(P,L) + QUAD(Q,L) END DO END DO ENDIF END DO END DO ENDIF END ! ----------------------------------------------------------------- ! BFLIST: list cells in breadth-first order, from largest (root) to ! smallest. Thanks to Jun Makino for this elegant routine. ! ----------------------------------------------------------------- SUBROUTINE BFLIST(IND) INTEGER IND(*) INCLUDE 'treedefs.f90' INTEGER FACELL, LACELL, NACELL, K, I ! ----------------------------------------- ! Start scan with root as only active cell. ! ----------------------------------------- IND(1) = ROOT FACELL = 1 LACELL = 1 ! ----------------------------------- ! Loop while active cells to process. ! ----------------------------------- DO WHILE (FACELL .LE. LACELL) ! ---------------------------------------------- ! Start counting active cells in next iteration. ! ---------------------------------------------- NACELL = LACELL ! --------------------------------------- ! Loop over subcells of each active cell. ! --------------------------------------- DO K = 1, NSUBC DO I = FACELL, LACELL ! ------------------------------------------- ! Add all cells on next level to active list. ! ------------------------------------------- IF (SUBP(IND(I),K) .GE. INCELL) THEN NACELL = NACELL + 1 IND(NACELL) = SUBP(IND(I),K) ENDIF END DO END DO ! ------------------------------------------------------ ! Advance first and last active cell indicies, and loop. ! ------------------------------------------------------ FACELL = LACELL + 1 LACELL = NACELL END DO ! -------------------------------------------------- ! Above loop should list all cells; check the count. ! -------------------------------------------------- IF (NACELL .NE. NCELL) & CALL TERROR(' BFLIST: INCONSISTENT CELL COUNT') END