+/**\r
+A D implementation of the _bh_ Olden benchmark.\r
+The Olden benchmark implements the Barnes-Hut benchmark\r
+that is decribed in:\r
+\r
+J. Barnes and P. Hut, "A hierarchical o(N log N) force-calculation algorithm",\r
+Nature, 324:446-449, Dec. 1986\r
+\r
+The original code in the Olden benchmark suite is derived from the\r
+ftp://hubble.ifa.hawaii.edu/pub/barnes/treecode\r
+source distributed by Barnes.\r
+\r
+Java code converted to D by leonardo maffi, V.1.0, Dec 25 2009.\r
+\r
+Removed output unless an option is passed by Leandro Lucarella, 2010-08-04.\r
+Downloaded from http://www.fantascienza.net/leonardo/js/\r
+ http://www.fantascienza.net/leonardo/js/dolden_bh.zip\r
+*/\r
+\r
+\r
+version (Tango) {\r
+ import tango.stdc.stdio: printf, sprintf, fprintf, stderr;\r
+ import tango.stdc.stdlib: exit, malloc, free;\r
+ import tango.stdc.time: CLOCKS_PER_SEC, clock;\r
+ import tango.math.Math: sqrt, floor, PI, pow;\r
+\r
+ import Integer = tango.text.convert.Integer;\r
+ alias Integer.parse toInt;\r
+} else {\r
+ import std.c.stdio: printf, sprintf, fprintf, stderr;\r
+ import std.c.stdlib: exit, malloc, free;\r
+ import std.c.time: CLOCKS_PER_SEC, clock;\r
+ import std.math: sqrt, floor, PI, pow;\r
+\r
+ version (D_Version2) {\r
+ import std.conv: to;\r
+ alias to!(int, char[]) toInt;\r
+ } else {\r
+ import std.conv: toInt;\r
+ }\r
+}\r
+\r
+\r
+double myclock() {\r
+ return clock() / cast(double)CLOCKS_PER_SEC;\r
+}\r
+\r
+\r
+/*\r
+Basic uniform random generator: Minimal Standard in Park and\r
+Miller (1988): "Random Number Generators: Good Ones Are Hard to\r
+Find", Comm. of the ACM, 31, 1192-1201.\r
+Parameters: m = 2^31-1, a=48271.\r
+\r
+Adapted from Pascal code by Jesper Lund:\r
+http://www.gnu-pascal.de/crystal/gpc/en/mail1390.html\r
+*/\r
+struct Random {\r
+ const int m = int.max;\r
+ const int a = 48_271;\r
+ const int q = m / a;\r
+ const int r = m % a;\r
+ int seed;\r
+\r
+ public double uniform(double min, double max) {\r
+ int k = seed / q;\r
+ seed = a * (seed - k * q) - r * k;\r
+ if (seed < 1)\r
+ seed += m;\r
+ double r = cast(double)seed / m;\r
+ return r * (max - min) + min;\r
+ }\r
+}\r
+\r
+\r
+interface Enumeration {\r
+ bool hasMoreElements();\r
+ Object nextElement();\r
+}\r
+\r
+\r
+/**\r
+A class representing a three dimensional vector that implements\r
+several math operations. To improve speed we implement the\r
+vector as an array of doubles rather than use the exising\r
+code in the java.util.Vec3 class.\r
+*/\r
+final class Vec3 {\r
+ /// The number of dimensions in the vector\r
+ // Can't be changed because of crossProduct()\r
+ const int NDIM = 3;\r
+\r
+ /// An array containing the values in the vector.\r
+ private double data[];\r
+\r
+ /// Construct an empty 3 dimensional vector for use in Barnes-Hut algorithm.\r
+ this() {\r
+ data.length = NDIM;\r
+ data[] = 0.0;\r
+ }\r
+\r
+ /// Create a copy of the vector. Return a clone of the math vector\r
+ public Vec3 clone() {\r
+ Vec3 v = new Vec3;\r
+ v.data.length = NDIM;\r
+ v.data[] = this.data[];\r
+ return v;\r
+ }\r
+\r
+ /**\r
+ Return the value at the i'th index of the vector.\r
+ @param i the vector index\r
+ @return the value at the i'th index of the vector.\r
+ */\r
+ final double value(int i) {\r
+ return data[i];\r
+ }\r
+\r
+ /**\r
+ Set the value of the i'th index of the vector.\r
+ @param i the vector index\r
+ @param v the value to store\r
+ */\r
+ final void value(int i, double v) {\r
+ data[i] = v;\r
+ }\r
+\r
+ /**\r
+ Set one of the dimensions of the vector to 1.0\r
+ param j the dimension to set.\r
+ */\r
+ final void unit(int j) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = (i == j) ? 1.0 : 0.0;\r
+ }\r
+\r
+ /**\r
+ Add two vectors and the result is placed in this vector.\r
+ @param u the other operand of the addition\r
+ */\r
+ final void addition(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] += u.data[i];\r
+ }\r
+\r
+ /**\r
+ * Subtract two vectors and the result is placed in this vector.\r
+ * This vector contain the first operand.\r
+ * @param u the other operand of the subtraction.\r
+ **/\r
+ final void subtraction(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] -= u.data[i];\r
+ }\r
+\r
+ /**\r
+ Subtract two vectors and the result is placed in this vector.\r
+ @param u the first operand of the subtraction.\r
+ @param v the second opernd of the subtraction\r
+ */\r
+ final void subtraction(Vec3 u, Vec3 v) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = u.data[i] - v.data[i];\r
+ }\r
+\r
+ /**\r
+ Multiply the vector times a scalar.\r
+ @param s the scalar value\r
+ **/\r
+ final void multScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] *= s;\r
+ }\r
+\r
+ /**\r
+ * Multiply the vector times a scalar and place the result in this vector.\r
+ * @param u the vector\r
+ * @param s the scalar value\r
+ **/\r
+ final void multScalar(Vec3 u, double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = u.data[i] * s;\r
+ }\r
+\r
+ /**\r
+ * Divide each element of the vector by a scalar value.\r
+ * @param s the scalar value.\r
+ **/\r
+ final void divScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] /= s;\r
+ }\r
+\r
+ /**\r
+ * Return the dot product of a vector.\r
+ * @return the dot product of a vector.\r
+ **/\r
+ final double dotProduct() {\r
+ double s = 0.0;\r
+ for (int i = 0; i < NDIM; i++)\r
+ s += data[i] * data[i];\r
+ return s;\r
+ }\r
+\r
+ final double absolute() {\r
+ double tmp = 0.0;\r
+ for (int i = 0; i < NDIM; i++)\r
+ tmp += data[i] * data[i];\r
+ return sqrt(tmp);\r
+ }\r
+\r
+ final double distance(Vec3 v) {\r
+ double tmp = 0.0;\r
+ for (int i = 0; i < NDIM; i++)\r
+ tmp += ((data[i] - v.data[i]) * (data[i] - v.data[i]));\r
+ return sqrt(tmp);\r
+ }\r
+\r
+ final void crossProduct(Vec3 u, Vec3 w) {\r
+ data[0] = u.data[1] * w.data[2] - u.data[2] * w.data[1];\r
+ data[1] = u.data[2] * w.data[0] - u.data[0] * w.data[2];\r
+ data[2] = u.data[0] * w.data[1] - u.data[1] * w.data[0];\r
+ }\r
+\r
+ final void incrementalAdd(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] += u.data[i];\r
+ }\r
+\r
+ final void incrementalSub(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] -= u.data[i];\r
+ }\r
+\r
+ final void incrementalMultScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] *= s;\r
+ }\r
+\r
+ final void incrementalDivScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] /= s;\r
+ }\r
+\r
+ final void addScalar(Vec3 u, double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = u.data[i] + s;\r
+ }\r
+\r
+ public char* toCString() {\r
+ // this is not generic code at all\r
+ char* result = cast(char*)malloc(100);\r
+ if (result == null) exit(1);\r
+ sprintf(result, "%.17f %.17f %.17f ", data[0], data[1], data[2]);\r
+ return result;\r
+ }\r
+}\r
+\r
+\r
+/// A class that represents the common fields of a cell or body data structure.\r
+abstract class Node {\r
+ // mass of the node\r
+ double mass;\r
+\r
+ // Position of the node\r
+ Vec3 pos;\r
+\r
+ // highest bit of int coord\r
+ const int IMAX = 1073741824;\r
+\r
+ // potential softening parameter\r
+ const double EPS = 0.05;\r
+\r
+ /// Construct an empty node\r
+ protected this() {\r
+ mass = 0.0;\r
+ pos = new Vec3();\r
+ }\r
+\r
+ abstract Node loadTree(Body p, Vec3 xpic, int l, Tree root);\r
+\r
+ abstract double hackcofm();\r
+\r
+ abstract HG walkSubTree(double dsq, HG hg);\r
+\r
+ static final int oldSubindex(Vec3 ic, int l) {\r
+ int i = 0;\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ if ((cast(int)ic.value(k) & l) != 0)\r
+ i += Cell.NSUB >> (k + 1);\r
+ return i;\r
+ }\r
+\r
+ public char* toCString() {\r
+ char* result = cast(char*)malloc(130);\r
+ if (result == null) exit(1);\r
+ char* pos_str = pos.toCString();\r
+ sprintf(result, "%f : %s", mass, pos_str);\r
+ free(pos_str);\r
+ return result;\r
+ }\r
+\r
+ /// Compute a single body-body or body-cell interaction\r
+ final HG gravSub(HG hg) {\r
+ Vec3 dr = new Vec3();\r
+ dr.subtraction(pos, hg.pos0);\r
+\r
+ double drsq = dr.dotProduct() + (EPS * EPS);\r
+ double drabs = sqrt(drsq);\r
+\r
+ double phii = mass / drabs;\r
+ hg.phi0 -= phii;\r
+ double mor3 = phii / drsq;\r
+ dr.multScalar(mor3);\r
+ hg.acc0.addition(dr);\r
+ return hg;\r
+ }\r
+\r
+ /**\r
+ * A sub class which is used to compute and save information during the\r
+ * gravity computation phase.\r
+ **/\r
+ static protected final class HG {\r
+ /// Body to skip in force evaluation\r
+ Body pskip;\r
+\r
+ /// Point at which to evaluate field\r
+ Vec3 pos0;\r
+\r
+ /// Computed potential at pos0\r
+\r
+ double phi0;\r
+\r
+ /// computed acceleration at pos0\r
+ Vec3 acc0;\r
+\r
+ /**\r
+ * Create a HG object.\r
+ * @param b the body object\r
+ * @param p a vector that represents the body\r
+ **/\r
+ this(Body b, Vec3 p) {\r
+ pskip = b;\r
+ pos0 = p.clone();\r
+ phi0 = 0.0;\r
+ acc0 = new Vec3();\r
+ }\r
+ }\r
+}\r
+\r
+\r
+/// A class used to representing particles in the N-body simulation.\r
+final class Body : Node {\r
+ Vec3 vel, acc, newAcc;\r
+ double phi = 0.0;\r
+ private Body next, procNext;\r
+\r
+ /// Create an empty body.\r
+ this() {\r
+ vel = new Vec3();\r
+ acc = new Vec3();\r
+ newAcc = new Vec3();\r
+ }\r
+\r
+ /**\r
+ * Set the next body in the list.\r
+ * @param n the body\r
+ **/\r
+ final void setNext(Body n) {\r
+ next = n;\r
+ }\r
+\r
+ /**\r
+ * Get the next body in the list.\r
+ * @return the next body\r
+ **/\r
+ final Body getNext() {\r
+ return next;\r
+ }\r
+\r
+ /**\r
+ * Set the next body in the list.\r
+ * @param n the body\r
+ **/\r
+ final void setProcNext(Body n) {\r
+ procNext = n;\r
+ }\r
+\r
+ /**\r
+ * Get the next body in the list.\r
+ * @return the next body\r
+ **/\r
+ final Body getProcNext() {\r
+ return procNext;\r
+ }\r
+\r
+ /**\r
+ * Enlarge cubical "box", salvaging existing tree structure.\r
+ * @param tree the root of the tree.\r
+ * @param nsteps the current time step\r
+ **/\r
+ final void expandBox(Tree tree, int nsteps) {\r
+ Vec3 rmid = new Vec3();\r
+\r
+ bool inbox = icTest(tree);\r
+ while (!inbox) {\r
+ double rsize = tree.rsize;\r
+ rmid.addScalar(tree.rmin, 0.5 * rsize);\r
+\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ if (pos.value(k) < rmid.value(k)) {\r
+ double rmin = tree.rmin.value(k);\r
+ tree.rmin.value(k, rmin - rsize);\r
+ }\r
+\r
+ tree.rsize = 2.0 * rsize;\r
+ if (tree.root !is null) {\r
+ Vec3 ic = tree.intcoord(rmid);\r
+ if (ic is null)\r
+ throw new Exception("Value is out of bounds");\r
+ int k = oldSubindex(ic, IMAX >> 1);\r
+ Cell newt = new Cell();\r
+ newt.subp[k] = tree.root;\r
+ tree.root = newt;\r
+ inbox = icTest(tree);\r
+ }\r
+ }\r
+ }\r
+\r
+ /**\r
+ * Check the bounds of the body and return true if it isn't in the\r
+ * correct bounds.\r
+ **/\r
+ final bool icTest(Tree tree) {\r
+ double pos0 = pos.value(0);\r
+ double pos1 = pos.value(1);\r
+ double pos2 = pos.value(2);\r
+\r
+ // by default, it is in bounds\r
+ bool result = true;\r
+\r
+ double xsc = (pos0 - tree.rmin.value(0)) / tree.rsize;\r
+ if (!(0.0 < xsc && xsc < 1.0))\r
+ result = false;\r
+\r
+ xsc = (pos1 - tree.rmin.value(1)) / tree.rsize;\r
+ if (!(0.0 < xsc && xsc < 1.0))\r
+ result = false;\r
+\r
+ xsc = (pos2 - tree.rmin.value(2)) / tree.rsize;\r
+ if (!(0.0 < xsc && xsc < 1.0))\r
+ result = false;\r
+\r
+ return result;\r
+ }\r
+\r
+ /**\r
+ * Descend Tree and insert particle. We're at a body so we need to\r
+ * create a cell and attach this body to the cell.\r
+ * @param p the body to insert\r
+ * @param xpic\r
+ * @param l\r
+ * @param tree the root of the data structure\r
+ * @return the subtree with the new body inserted\r
+ **/\r
+ final Node loadTree(Body p, Vec3 xpic, int l, Tree tree) {\r
+ // create a Cell\r
+ Cell retval = new Cell();\r
+ int si = subindex(tree, l);\r
+ // attach this Body node to the cell\r
+ retval.subp[si] = this;\r
+\r
+ // move down one level\r
+ si = oldSubindex(xpic, l);\r
+ Node rt = retval.subp[si];\r
+ if (rt !is null)\r
+ retval.subp[si] = rt.loadTree(p, xpic, l >> 1, tree);\r
+ else\r
+ retval.subp[si] = p;\r
+ return retval;\r
+ }\r
+\r
+ /**\r
+ * Descend tree finding center of mass coordinates\r
+ * @return the mass of this node\r
+ **/\r
+ final double hackcofm() {\r
+ return mass;\r
+ }\r
+\r
+ /// iteration of the bodies\r
+ int opApply(int delegate(ref Body) dg) {\r
+ int result;\r
+ Body current = this;\r
+ while (current !is null) {\r
+ result = dg(current);\r
+ current = current.next;\r
+ if (result)\r
+ break;\r
+ }\r
+ return result;\r
+ }\r
+\r
+ /// inverse iteration of the bodies\r
+ int opApplyReverse(int delegate(ref Body) dg) {\r
+ int result;\r
+ Body current = this;\r
+ while (current !is null) {\r
+ result = dg(current);\r
+ current = current.procNext;\r
+ if (result)\r
+ break;\r
+ }\r
+ return result;\r
+ }\r
+\r
+\r
+ /**\r
+ * Return an enumeration of the bodies\r
+ * @return an enumeration of the bodies\r
+ **/\r
+ final Enumeration elements() {\r
+ // a local class that implements the enumerator\r
+ static final class Enumerate : Enumeration {\r
+ private Body current;\r
+ public this(Body outer) {\r
+ //this.current = Body.this;\r
+ this.current = outer;\r
+ }\r
+ public bool hasMoreElements() {\r
+ return current !is null;\r
+ }\r
+ public Object nextElement() {\r
+ Object retval = current;\r
+ current = current.next;\r
+ return retval;\r
+ }\r
+ }\r
+\r
+ return new Enumerate(this);\r
+ }\r
+\r
+ final Enumeration elementsRev() {\r
+ // a local class that implements the enumerator\r
+ static class Enumerate : Enumeration {\r
+ private Body current;\r
+ public this(Body outer) {\r
+ //this.current = Body.this;\r
+ this.current = outer;\r
+ }\r
+ public bool hasMoreElements() {\r
+ return current !is null;\r
+ }\r
+ public Object nextElement() {\r
+ Object retval = current;\r
+ current = current.procNext;\r
+ return retval;\r
+ }\r
+ }\r
+\r
+ return new Enumerate(this);\r
+ }\r
+\r
+ /**\r
+ * Determine which subcell to select.\r
+ * Combination of intcoord and oldSubindex.\r
+ * @param t the root of the tree\r
+ **/\r
+ final int subindex(Tree tree, int l) {\r
+ Vec3 xp = new Vec3();\r
+\r
+ double xsc = (pos.value(0) - tree.rmin.value(0)) / tree.rsize;\r
+ xp.value(0, floor(IMAX * xsc));\r
+\r
+ xsc = (pos.value(1) - tree.rmin.value(1)) / tree.rsize;\r
+ xp.value(1, floor(IMAX * xsc));\r
+\r
+ xsc = (pos.value(2) - tree.rmin.value(2)) / tree.rsize;\r
+ xp.value(2, floor(IMAX * xsc));\r
+\r
+ int i = 0;\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ if ((cast(int)xp.value(k) & l) != 0)\r
+ i += Cell.NSUB >> (k + 1);\r
+ return i;\r
+ }\r
+\r
+ /**\r
+ * Evaluate gravitational field on the body.\r
+ * The original olden version calls a routine named "walkscan",\r
+ * but we use the same name that is in the Barnes code.\r
+ **/\r
+ final void hackGravity(double rsize, Node root) {\r
+ Vec3 pos0 = pos.clone();\r
+ HG hg = new HG(this, pos);\r
+ hg = root.walkSubTree(rsize * rsize, hg);\r
+ phi = hg.phi0;\r
+ newAcc = hg.acc0;\r
+ }\r
+\r
+ /// Recursively walk the tree to do hackwalk calculation\r
+ final HG walkSubTree(double dsq, HG hg) {\r
+ if (this != hg.pskip)\r
+ hg = gravSub(hg);\r
+ return hg;\r
+ }\r
+\r
+ /**\r
+ * Return a string represenation of a body.\r
+ * @return a string represenation of a body.\r
+ **/\r
+ public char* toCString() {\r
+ char* result = cast(char*)malloc(140);\r
+ if (result == null) exit(1);\r
+ char* super_str = super.toCString();\r
+ sprintf(result, "Body %s", super_str);\r
+ free(super_str);\r
+ return result;\r
+ }\r
+}\r
+\r
+\r
+\r
+/// A class used to represent internal nodes in the tree\r
+final class Cell : Node {\r
+ // subcells per cell\r
+ const int NSUB = 1 << Vec3.NDIM;\r
+\r
+ /**\r
+ * The children of this cell node. Each entry may contain either\r
+ * another cell or a body.\r
+ **/\r
+ Node[] subp;\r
+ Cell next;\r
+\r
+ this() {\r
+ subp.length = NSUB;\r
+ }\r
+\r
+ /**\r
+ * Descend Tree and insert particle. We're at a cell so\r
+ * we need to move down the tree.\r
+ * @param p the body to insert into the tree\r
+ * @param xpic\r
+ * @param l\r
+ * @param tree the root of the tree\r
+ * @return the subtree with the new body inserted\r
+ **/\r
+ Node loadTree(Body p, Vec3 xpic, int l, Tree tree) {\r
+ // move down one level\r
+ int si = oldSubindex(xpic, l);\r
+ Node rt = subp[si];\r
+ if (rt !is null)\r
+ subp[si] = rt.loadTree(p, xpic, l >> 1, tree);\r
+ else\r
+ subp[si] = p;\r
+ return this;\r
+ }\r
+\r
+ /**\r
+ * Descend tree finding center of mass coordinates\r
+ * @return the mass of this node\r
+ **/\r
+ double hackcofm() {\r
+ double mq = 0.0;\r
+ Vec3 tmpPos = new Vec3();\r
+ Vec3 tmpv = new Vec3();\r
+ for (int i = 0; i < NSUB; i++) {\r
+ Node r = subp[i];\r
+ if (r !is null) {\r
+ double mr = r.hackcofm();\r
+ mq = mr + mq;\r
+ tmpv.multScalar(r.pos, mr);\r
+ tmpPos.addition(tmpv);\r
+ }\r
+ }\r
+ mass = mq;\r
+ pos = tmpPos;\r
+ pos.divScalar(mass);\r
+\r
+ return mq;\r
+ }\r
+\r
+ /// Recursively walk the tree to do hackwalk calculation\r
+ HG walkSubTree(double dsq, HG hg) {\r
+ if (subdivp(dsq, hg)) {\r
+ for (int k = 0; k < Cell.NSUB; k++) {\r
+ Node r = subp[k];\r
+ if (r !is null)\r
+ hg = r.walkSubTree(dsq / 4.0, hg);\r
+ }\r
+ } else {\r
+ hg = gravSub(hg);\r
+ }\r
+ return hg;\r
+ }\r
+\r
+ /**\r
+ * Decide if the cell is too close to accept as a single term.\r
+ * @return true if the cell is too close.\r
+ **/\r
+ bool subdivp(double dsq, HG hg) {\r
+ Vec3 dr = new Vec3();\r
+ dr.subtraction(pos, hg.pos0);\r
+ double drsq = dr.dotProduct();\r
+\r
+ // in the original olden version drsp is multiplied by 1.0\r
+ return drsq < dsq;\r
+ }\r
+\r
+ /**\r
+ * Return a string represenation of a cell.\r
+ * @return a string represenation of a cell.\r
+ **/\r
+ public char* toCString() {\r
+ char* result = cast(char*)malloc(140);\r
+ if (result == null) exit(1);\r
+ char* super_str = super.toCString();\r
+ sprintf(result, "Cell %s", super_str);\r
+ free(super_str);\r
+ return result;\r
+ }\r
+}\r
+\r
+\r
+/**\r
+* A class that represents the root of the data structure used\r
+* to represent the N-bodies in the Barnes-Hut algorithm.\r
+**/\r
+final class Tree {\r
+ Vec3 rmin;\r
+ double rsize;\r
+\r
+ /// A reference to the root node.\r
+ Node root;\r
+\r
+ /// The complete list of bodies that have been created.\r
+ private Body bodyTab;\r
+\r
+ /// The complete list of bodies that have been created - in reverse.\r
+ private Body bodyTabRev;\r
+\r
+ /// Construct the root of the data structure that represents the N-bodies.\r
+ this() {\r
+ rmin = new Vec3();\r
+ rsize = -2.0 * -2.0;\r
+ root = null;\r
+ bodyTab = null;\r
+ bodyTabRev = null;\r
+ rmin.value(0, -2.0);\r
+ rmin.value(1, -2.0);\r
+ rmin.value(2, -2.0);\r
+ }\r
+\r
+ /**\r
+ * Return an enumeration of the bodies.\r
+ * @return an enumeration of the bodies.\r
+ **/\r
+ final Enumeration bodies() {\r
+ return bodyTab.elements();\r
+ }\r
+\r
+ /**\r
+ * Return an enumeration of the bodies - in reverse.\r
+ * @return an enumeration of the bodies - in reverse.\r
+ **/\r
+ final Enumeration bodiesRev() {\r
+ return bodyTabRev.elementsRev();\r
+ }\r
+\r
+ /**\r
+ * Create the testdata used in the benchmark.\r
+ * @param nbody the number of bodies to create\r
+ **/\r
+ final void createTestData(int nbody) {\r
+ Vec3 cmr = new Vec3();\r
+ Vec3 cmv = new Vec3();\r
+ Body head = new Body();\r
+ Body prev = head;\r
+\r
+ double rsc = 3.0 * PI / 16.0;\r
+ double vsc = sqrt(1.0 / rsc);\r
+ int seed = 123;\r
+ Random rnd = Random(seed);\r
+\r
+ for (int i = 0; i < nbody; i++) {\r
+ Body p = new Body();\r
+\r
+ prev.setNext(p);\r
+ prev = p;\r
+ p.mass = 1.0 / cast(double)nbody;\r
+\r
+ double t1 = rnd.uniform(0.0, 0.999);\r
+ t1 = pow(t1, (-2.0 / 3.0)) - 1.0;\r
+ double r = 1.0 / sqrt(t1);\r
+\r
+ double coeff = 4.0;\r
+ for (int k = 0; k < Vec3.NDIM; k++) {\r
+ r = rnd.uniform(0.0, 0.999);\r
+ p.pos.value(k, coeff * r);\r
+ }\r
+\r
+ cmr.addition(p.pos);\r
+\r
+ double x, y;\r
+ do {\r
+ x = rnd.uniform(0.0, 1.0);\r
+ y = rnd.uniform(0.0, 0.1);\r
+ } while (y > (x * x * pow(1.0 - x * x, 3.5)));\r
+\r
+ double v = sqrt(2.0) * x / pow(1 + r * r, 0.25);\r
+\r
+ double rad = vsc * v;\r
+ double rsq;\r
+ do {\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ p.vel.value(k, rnd.uniform(-1.0, 1.0));\r
+ rsq = p.vel.dotProduct();\r
+ } while (rsq > 1.0);\r
+ double rsc1 = rad / sqrt(rsq);\r
+ p.vel.multScalar(rsc1);\r
+ cmv.addition(p.vel);\r
+ }\r
+\r
+ // mark end of list\r
+ prev.setNext(null);\r
+\r
+ // toss the dummy node at the beginning and set a reference to the first element\r
+ bodyTab = head.getNext();\r
+\r
+ cmr.divScalar(cast(double)nbody);\r
+ cmv.divScalar(cast(double)nbody);\r
+\r
+ prev = null;\r
+\r
+ for (Enumeration e = bodyTab.elements(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ b.pos.subtraction(cmr);\r
+ b.vel.subtraction(cmv);\r
+ b.setProcNext(prev);\r
+ prev = b;\r
+ }\r
+\r
+ // set the reference to the last element\r
+ bodyTabRev = prev;\r
+ }\r
+\r
+\r
+ /**\r
+ * Advance the N-body system one time-step.\r
+ * @param nstep the current time step\r
+ **/\r
+ void stepSystem(int nstep) {\r
+ // free the tree\r
+ root = null;\r
+\r
+ makeTree(nstep);\r
+\r
+ // compute the gravity for all the particles\r
+ for (Enumeration e = bodyTabRev.elementsRev(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ b.hackGravity(rsize, root);\r
+ }\r
+ vp(bodyTabRev, nstep);\r
+ }\r
+\r
+ /**\r
+ * Initialize the tree structure for hack force calculation.\r
+ * @param nsteps the current time step\r
+ **/\r
+ private void makeTree(int nstep) {\r
+ for (Enumeration e = bodiesRev(); e.hasMoreElements();) {\r
+ Body q = cast(Body)e.nextElement();\r
+ if (q.mass != 0.0) {\r
+ q.expandBox(this, nstep);\r
+ Vec3 xqic = intcoord(q.pos);\r
+ if (root is null) {\r
+ root = q;\r
+ } else {\r
+ root = root.loadTree(q, xqic, Node.IMAX >> 1, this);\r
+ }\r
+ }\r
+ }\r
+ root.hackcofm();\r
+ }\r
+\r
+ /**\r
+ * Compute integerized coordinates.\r
+ * @return the coordinates or null if rp is out of bounds\r
+ **/\r
+ final Vec3 intcoord(Vec3 vp) {\r
+ Vec3 xp = new Vec3();\r
+\r
+ double xsc = (vp.value(0) - rmin.value(0)) / rsize;\r
+ if (0.0 <= xsc && xsc < 1.0)\r
+ xp.value(0, floor(Node.IMAX * xsc));\r
+ else\r
+ return null;\r
+\r
+ xsc = (vp.value(1) - rmin.value(1)) / rsize;\r
+ if (0.0 <= xsc && xsc < 1.0)\r
+ xp.value(1, floor(Node.IMAX * xsc));\r
+ else\r
+ return null;\r
+\r
+ xsc = (vp.value(2) - rmin.value(2)) / rsize;\r
+ if (0.0 <= xsc && xsc < 1.0)\r
+ xp.value(2, floor(Node.IMAX * xsc));\r
+ else\r
+ return null;\r
+\r
+ return xp;\r
+ }\r
+\r
+ static final private void vp(Body p, int nstep) {\r
+ Vec3 dacc = new Vec3();\r
+ Vec3 dvel = new Vec3();\r
+ double dthf = 0.5 * BH.DTIME;\r
+\r
+ for (Enumeration e = p.elementsRev(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ Vec3 acc1 = b.newAcc.clone();\r
+ if (nstep > 0) {\r
+ dacc.subtraction(acc1, b.acc);\r
+ dvel.multScalar(dacc, dthf);\r
+ dvel.addition(b.vel);\r
+ b.vel = dvel.clone();\r
+ }\r
+\r
+ b.acc = acc1.clone();\r
+ dvel.multScalar(b.acc, dthf);\r
+\r
+ Vec3 vel1 = b.vel.clone();\r
+ vel1.addition(dvel);\r
+ Vec3 dpos = vel1.clone();\r
+ dpos.multScalar(BH.DTIME);\r
+ dpos.addition(b.pos);\r
+ b.pos = dpos.clone();\r
+ vel1.addition(dvel);\r
+ b.vel = vel1.clone();\r
+ }\r
+ }\r
+}\r
+\r
+\r
+final class BH {\r
+ /// The user specified number of bodies to create.\r
+ private static int nbody = 0;\r
+\r
+ /// The maximum number of time steps to take in the simulation\r
+ private static int nsteps = 10;\r
+\r
+ /// Should we print information messsages\r
+ private static bool printMsgs = false;\r
+\r
+ /// Should we print detailed results\r
+ private static bool printResults = false;\r
+\r
+ const double DTIME = 0.0125;\r
+ private const double TSTOP = 2.0;\r
+\r
+ public static final void main(char[][] args) {\r
+ parseCmdLine(args);\r
+\r
+ if (printMsgs)\r
+ printf("nbody = %d\n", nbody);\r
+\r
+ auto start0 = myclock();\r
+ Tree root = new Tree();\r
+ root.createTestData(nbody);\r
+ auto end0 = myclock();\r
+ if (printMsgs)\r
+ printf("Bodies created\n");\r
+\r
+ auto start1 = myclock();\r
+ double tnow = 0.0;\r
+ int i = 0;\r
+ while ((tnow < TSTOP + 0.1 * DTIME) && (i < nsteps)) {\r
+ root.stepSystem(i++);\r
+ tnow += DTIME;\r
+ }\r
+ auto end1 = myclock();\r
+\r
+ if (printResults) {\r
+ int j = 0;\r
+ for (Enumeration e = root.bodies(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ char* str_ptr = b.pos.toCString();\r
+ printf("body %d: %s\n", j++, str_ptr);\r
+ free(str_ptr);\r
+ }\r
+ }\r
+\r
+ if (printMsgs) {\r
+ printf("Build time %.2f\n", end0 - start0);\r
+ printf("Compute Time %.2f\n", end1 - start1);\r
+ printf("Total Time %.2f\n", end1 - start0);\r
+ printf("Done!\n");\r
+ }\r
+ }\r
+\r
+ private static final void parseCmdLine(char[][] args) {\r
+ int i = 1;\r
+ while (i < args.length && args[i][0] == '-') {\r
+ char[] arg = args[i++];\r
+\r
+ // check for options that require arguments\r
+ if (arg == "-b") {\r
+ if (i < args.length)\r
+ nbody = toInt(args[i++]);\r
+ else\r
+ throw new Exception("-l requires the number of levels");\r
+ } else if (arg == "-s") {\r
+ if (i < args.length)\r
+ nsteps = toInt(args[i++]);\r
+ else\r
+ throw new Exception("-l requires the number of levels");\r
+ } else if (arg == "-m") {\r
+ printMsgs = true;\r
+ } else if (arg == "-p") {\r
+ printResults = true;\r
+ } else if (arg == "-h") {\r
+ usage();\r
+ }\r
+ }\r
+\r
+ if (nbody == 0)\r
+ usage();\r
+ }\r
+\r
+ /// The usage routine which describes the program options.\r
+ private static final void usage() {\r
+ fprintf(stderr, "usage: BH -b <size> [-s <steps>] [-p] [-m] [-h]\n");\r
+ fprintf(stderr, " -b the number of bodies\n");\r
+ fprintf(stderr, " -s the max. number of time steps (default=10)\n");\r
+ fprintf(stderr, " -p (print detailed results)\n");\r
+ fprintf(stderr, " -m (print information messages\n");\r
+ exit(0);\r
+ }\r
+}\r
+\r
+\r
+void main(char[][] args) {\r
+ BH.main(args);\r
+}\r
projects / software / dgc / dgcbench.git / commitdiff