Home Esplora Aiuto
Registrati Accedi
jarjonam
/
tidal_tails
1
0
Forka 0
File
Albero (Tree): 06320b2395
Rami (Branch) Tag
tidal_tails
/
docs
/
acceleration.html

acceleration.html 23 KB
Cronologia Originale

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432 
  1. <!doctype html>
    2. 

  2. <html lang="en">
    3. 

  3. <head>
    4. 

  4. <meta charset="utf-8">
    5. 

  5. <meta name="viewport" content="width=device-width, initial-scale=1, minimum-scale=1" />
    6. 

  6. <meta name="generator" content="pdoc 0.5.2" />
    7. 

  7. <title>acceleration API documentation</title>
    8. 

  8. <meta name="description" content="Defines the possible routines for computing the gravitational forces in the
    9. 

  9. simulation …" />
    10. 

  10. <link href='https://cdnjs.cloudflare.com/ajax/libs/normalize/8.0.0/normalize.min.css' rel='stylesheet'>
    11. 

  11. <link href='https://cdnjs.cloudflare.com/ajax/libs/10up-sanitize.css/8.0.0/sanitize.min.css' rel='stylesheet'>
    12. 

  12. <link href="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.12.0/styles/github.min.css" rel="stylesheet">
    13. 

  13. <style>.flex{display:flex !important}body{line-height:1.5em}#content{padding:20px}#sidebar{padding:30px;overflow:hidden}.http-server-breadcrumbs{font-size:130%;margin:0 0 15px 0}#footer{font-size:.75em;padding:5px 30px;border-top:1px solid #ddd;text-align:right}#footer p{margin:0 0 0 1em;display:inline-block}#footer p:last-child{margin-right:30px}h1,h2,h3,h4,h5{font-weight:300}h1{font-size:2.5em;line-height:1.1em}h2{font-size:1.75em;margin:1em 0 .50em 0}h3{font-size:1.4em;margin:25px 0 10px 0}h4{margin:0;font-size:105%}a{color:#058;text-decoration:none;transition:color .3s ease-in-out}a:hover{color:#e82}.title code{font-weight:bold}h2[id^="header-"]{margin-top:2em}.ident{color:#900}pre code{background:#f8f8f8;font-size:.8em;line-height:1.4em}code{background:#f2f2f1;padding:1px 4px;overflow-wrap:break-word}h1 code{background:transparent}pre{background:#f8f8f8;border:0;border-top:1px solid #ccc;border-bottom:1px solid #ccc;margin:1em 0;padding:1ex}#http-server-module-list{display:flex;flex-flow:column}#http-server-module-list div{display:flex}#http-server-module-list dt{min-width:10%}#http-server-module-list p{margin-top:0}.toc ul,#index{list-style-type:none;margin:0;padding:0}#index code{background:transparent}#index h3{border-bottom:1px solid #ddd}#index ul{padding:0}#index h4{font-weight:bold}#index h4 + ul{margin-bottom:.6em}#index .two-column{column-count:2}dl{margin-bottom:2em}dl dl:last-child{margin-bottom:4em}dd{margin:0 0 1em 3em}#header-classes + dl > dd{margin-bottom:3em}dd dd{margin-left:2em}dd p{margin:10px 0}.name{background:#eee;font-weight:bold;font-size:.85em;padding:5px 10px;display:inline-block;min-width:40%}.name:hover{background:#e0e0e0}.name > span:first-child{white-space:nowrap}.name.class > span:nth-child(2){margin-left:.4em}.name small{font-weight:normal}.inherited{color:#999;border-left:5px solid #eee;padding-left:1em}.inheritance em{font-style:normal;font-weight:bold}.desc h2{font-weight:400;font-size:1.25em}.desc h3{font-size:1em}.desc dt code{background:inherit}.source summary{color:#666;text-align:right;font-weight:400;font-size:.8em;text-transform:uppercase;cursor:pointer}.source pre{max-height:500px;overflow:auto;margin:0}.source pre code{overflow:visible}.hlist{list-style:none}.hlist li{display:inline}.hlist li:after{content:',\2002'}.hlist li:last-child:after{content:none}.hlist .hlist{display:inline;padding-left:1em}img{max-width:100%}.admonition{padding:.1em .5em}.admonition-title{font-weight:bold}.admonition.note,.admonition.info,.admonition.important{background:#aef}.admonition.todo,.admonition.versionadded,.admonition.tip,.admonition.hint{background:#dfd}.admonition.warning,.admonition.versionchanged,.admonition.deprecated{background:#fd4}.admonition.error,.admonition.danger,.admonition.caution{background:lightpink}</style>
    14. 

  14. <style media="screen and (min-width: 700px)">@media screen and (min-width:700px){#sidebar{width:30%}#content{width:70%;max-width:100ch;padding:3em 4em;border-left:1px solid #ddd}pre code{font-size:1em}.item .name{font-size:1em}main{display:flex;flex-direction:row-reverse;justify-content:flex-end}.toc ul ul,#index ul{padding-left:1.5em}.toc > ul > li{margin-top:.5em}}</style>
    15. 

  15. <style media="print">@media print{#sidebar h1{page-break-before:always}.source{display:none}}@media print{*{background:transparent !important;color:#000 !important;box-shadow:none !important;text-shadow:none !important}a[href]:after{content:" (" attr(href) ")";font-size:90%}a[href][title]:after{content:none}abbr[title]:after{content:" (" attr(title) ")"}.ir a:after,a[href^="javascript:"]:after,a[href^="#"]:after{content:""}pre,blockquote{border:1px solid #999;page-break-inside:avoid}thead{display:table-header-group}tr,img{page-break-inside:avoid}img{max-width:100% !important}@page{margin:0.5cm}p,h2,h3{orphans:3;widows:3}h1,h2,h3,h4,h5,h6{page-break-after:avoid}}</style>
    16. 

  16. </head>
    17. 

  17. <body>
    18. 

  18. <main>
    19. 

  19. <article id="content">
    20. 

  20. <header>
    21. 

  21. <h1 class="title"><code>acceleration</code> module</h1>
    22. 

  22. </header>
    23. 

  23. <section id="section-intro">
    24. 

  24. <p>Defines the possible routines for computing the gravitational forces in the
    25. 

  25. simulation.</p>
    26. 

  26. <p>All the methods in these file require a position (n, 3) vector, a mass (n, )
    27. 

  27. vector and an optional softening scale float.</p>
    28. 

  28. <details class="source">
    29. 

  29. <summary>Source code</summary>
    30. 

  30. <pre><code class="python">&#34;&#34;&#34;Defines the possible routines for computing the gravitational forces in the
    31. 

  31. simulation.
    32. 

  32. 

  33. 

  34. All the methods in these file require a position (n, 3) vector, a mass (n, )
    35. 

  35. vector and an optional softening scale float.&#34;&#34;&#34;
    36. 

  36. 

  37. import ctypes
    38. 

  38. import numpy.ctypeslib as ctl
    39. 

  39. import numpy as np
    40. 

  40. from numba import jit
    41. 

  41. 

  42. 

  43. def bruteForce(r_vec, mass, soft=0.):
    44. 

  44.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    45. 

  45.        Complexity O(n*m) where n is the total number of masses and m is the
    46. 

  46.        number of massive particles.
    47. 

  47. 

  48.     Parameters:
    49. 

  49.         r_vec (array): list of particles positions.
    50. 

  50.             Shape (n, 3) where n is the number of particles
    51. 

  51.         mass (array): list of particles masses.
    52. 

  52.             Shape (n,)
    53. 

  53.         soft (float): characteristic plummer softening length scale
    54. 

  54.     Returns:
    55. 

  55.         forces (array): list of forces acting on each particle.
    56. 

  56.             Shape (n, 3)
    57. 

  57.     &#34;&#34;&#34;
    58. 

  58.     # Only calculate forces from massive particles
    59. 

  59.     mask = mass!=0
    60. 

  60.     massMassive = mass[mask]
    61. 

  61.     rMassive_vec = r_vec[mask]
    62. 

  62.     #  x m x 1 matrix (m = number of massive particles) for broadcasting
    63. 

  63.     mass_mat = massMassive.reshape(1, -1, 1)
    64. 

  64.     # Calculate displacements
    65. 

  65.     # r_ten is the direction of the pairwise displacements. Shape (n, m, 3)
    66. 

  66.     # r_mat is the absolute distance of the pairwise displacements. (n, m, 1)
    67. 

  67.     r_ten = rMassive_vec.reshape(1, -1, 3) - r_vec.reshape(-1, 1, 3)
    68. 

  68.     r_mat = np.linalg.norm(r_ten, axis=-1, keepdims=True)
    69. 

  69.     # Avoid division by zeros
    70. 

  70.     # $a = M / (r + \epsilon)^2$, where $\epsilon$ is the softening scale
    71. 

  71.     # r_ten/r_mat gives the direction unit vector
    72. 

  72.     accel = np.divide(r_ten * mass_mat/(r_mat+soft)**2, r_mat, 
    73. 

  73.         where=r_ten.astype(bool), out=r_ten) # Reuse memory from r_ten
    74. 

  74.     return accel.sum(axis=1) # Add all forces on each particle
    75. 

  75. 

  76. @jit(nopython=True) # Numba annotation
    77. 

  77. def bruteForceNumba(r_vec, mass, soft=0.):
    78. 

  78.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    79. 

  79.     It is done in the same way as in bruteForce, but this
    80. 

  80.     method is ran through Numba&#34;&#34;&#34;
    81. 

  81.     mask = mass!=0
    82. 

  82.     massMassive = mass[mask]
    83. 

  83.     rMassive_vec = r_vec[mask]
    84. 

  84.     mass_mat = massMassive.reshape(1, -1, 1)
    85. 

  85.     r_ten = rMassive_vec.reshape(1, -1, 3) - r_vec.reshape(-1, 1, 3)
    86. 

  86.     # Avoid np.linalg.norm to allow Numba optimizations
    87. 

  87.     r_mat = np.sqrt(r_ten[:,:,0:1]**2 + r_ten[:,:,1:2]**2 + r_ten[:,:,2:3]**2)
    88. 

  88.     r_mat = np.where(r_mat == 0, np.ones_like(r_mat), r_mat)
    89. 

  89.     accel = r_ten/r_mat * mass_mat/(r_mat+soft)**2
    90. 

  90.     return accel.sum(axis=1) # Add all forces in each particle
    91. 

  91. 

  92. @jit(nopython=True) # Numba annotation
    93. 

  93. def bruteForceNumbaOptimized(r_vec, mass, soft=0.):
    94. 

  94.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    95. 

  95.     This is optimized for high performance with Numba. All massive particles
    96. 

  96.     must appear first.&#34;&#34;&#34;
    97. 

  97.     accel = np.zeros_like(r_vec) 
    98. 

  98.     # Use superposition to add all the contributions
    99. 

  99.     n = r_vec.shape[0] # Number of particles
    100. 

  100.     delta = np.zeros((3,)) # Only allocate this once
    101. 

  101.     for i in range(n):
    102. 

  102.         # Only consider pairs with at least one massive particle i
    103. 

  103.         if mass[i] == 0: break
    104. 

  104.         for j in range(i+1, n):
    105. 

  105.             # Explicitely separate components for high performance
    106. 

  106.             # i.e. do not do delta = r_vec[j] - r_vec[i]
    107. 

  107.             # (The effect of this is VERY relevant (x10) and has to do with
    108. 

  108.             # memory reallocation) Numba will vectorize the loops.
    109. 

  109.             for k in range(3): delta[k] = r_vec[j,k] - r_vec[i,k]
    110. 

  110.             r = np.sqrt(delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2])
    111. 

  111.             tripler = (r+soft)**2 * r
    112. 

  112.             
    113. 

  113.             # Compute acceleration on first particle
    114. 

  114.             mr3inv = mass[i]/(tripler)
    115. 

  115.             # Again, do NOT do accel[j] -= mr3inv * delta
    116. 

  116.             for k in range(3): accel[j,k] -= mr3inv * delta[k]
    117. 

  117.             
    118. 

  118.             # Compute acceleration on second particle
    119. 

  119.             # For pairs with one massless particle, no reaction force
    120. 

  120.             if mass[j] == 0: break
    121. 

  121.             # Otherwise, opposite direction (+)
    122. 

  122.             mr3inv = mass[j]/(tripler)
    123. 

  123.             for k in range(3): accel[i,k] += mr3inv * delta[k]
    124. 

  124.     return accel
    125. 

  125. 

  126. # C++ interface, load library
    127. 

  127. ACCLIB = None
    128. 

  128. def loadCPPLib():
    129. 

  129.     &#34;&#34;&#34;Loads the C++ shared library to the global variable ACCLIB. Must be
    130. 

  130.     called before using the library.&#34;&#34;&#34;
    131. 

  131.     global ACCLIB
    132. 

  132.     ACCLIB = ctypes.CDLL(&#39;cpp/acclib.so&#39;)
    133. 

  133.     # Define appropiate types for library functions
    134. 

  134.     doublepp = np.ctypeslib.ndpointer(dtype=np.uintp) # double**
    135. 

  135.     doublep = ctl.ndpointer(np.float64, flags=&#39;aligned, c_contiguous&#39;)#double*
    136. 

  136.     # Check cpp/acclib.cpp for function signatures
    137. 

  137.     ACCLIB.bruteForceCPP.argtypes = [doublepp, doublep, 
    138. 

  138.         ctypes.c_int, ctypes.c_double]
    139. 

  139.     ACCLIB.barnesHutCPP.argtypes = [doublepp, doublep, 
    140. 

  140.         ctypes.c_int, ctypes.c_double, ctypes.c_double, 
    141. 

  141.         ctypes.c_double, ctypes.c_double, ctypes.c_double]
    142. 

  142. 

  143. def bruteForceCPP(r_vec, m_vec, soft=0.):
    144. 

  144.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    145. 

  145.     This is ran in a shared C++ library through Brute Force (pairwise sums)
    146. 

  146.     Massive particles must appear first.&#34;&#34;&#34;
    147. 

  147.     # Convert array to data required by C++ library
    148. 

  148.     if ACCLIB is None: loadCPPLib() # Singleton pattern
    149. 

  149.     # Change type to be appropiate for calling library
    150. 

  150.     r_vec_c = (r_vec.ctypes.data + np.arange(r_vec.shape[0]) 
    151. 

  151.         * r_vec.strides[0]).astype(np.uintp)
    152. 

  152.     # Set return type as double*
    153. 

  153.     ACCLIB.bruteForceCPP.restype = np.ctypeslib.ndpointer(dtype=np.float64, 
    154. 

  154.         shape=(r_vec.shape[0]*3,))
    155. 

  155.     # Call the C++ function: double* bruteForceCPP
    156. 

  156.     accel =  ACCLIB.bruteForceCPP(r_vec_c, m_vec, r_vec.shape[0], soft)
    157. 

  157.     # Change shape to get the expected Numpy array (n, 3)
    158. 

  158.     accel.shape = (-1, 3)
    159. 

  159.     return accel
    160. 

  160. 

  161. def barnesHutCPP(r_vec, m_vec, soft=0.):
    162. 

  162.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    163. 

  163.     This is ran in a shared C++ library using a BarnesHut tree&#34;&#34;&#34;
    164. 

  164.     # Convert array to data required by C++ library
    165. 

  165.     if ACCLIB is None: loadCPPLib() # Singleton pattern
    166. 

  166.     # Change type to be appropiate for calling library
    167. 

  167.     r_vec_c = (r_vec.ctypes.data + np.arange(r_vec.shape[0]) 
    168. 

  168.         * r_vec.strides[0]).astype(np.uintp)
    169. 

  169.     # Set return type as double*
    170. 

  170.     ACCLIB.barnesHutCPP.restype = np.ctypeslib.ndpointer(dtype=np.float64, 
    171. 

  171.         shape=(r_vec.shape[0]*3,))
    172. 

  172.     # Explicitely pass the corner and size of the box for the top node
    173. 

  173.     px, py, pz = np.min(r_vec, axis=0)
    174. 

  174.     size = np.max(np.max(r_vec, axis=0) - np.min(r_vec, axis=0))
    175. 

  175.     # Call the C++ function: double* barnesHutCPP
    176. 

  176.     accel =  ACCLIB.barnesHutCPP(r_vec_c, m_vec, r_vec.shape[0], 
    177. 

  177.         size, px, py, pz, soft)
    178. 

  178.     # Change shape to get the expected Numpy array (n, 3)
    179. 

  179.     accel.shape = (-1, 3)
    180. 

  180.     return accel</code></pre>
    181. 

  181. </details>
    182. 

  182. </section>
    183. 

  183. <section>
    184. 

  184. </section>
    185. 

  185. <section>
    186. 

  186. </section>
    187. 

  187. <section>
    188. 

  188. <h2 class="section-title" id="header-functions">Functions</h2>
    189. 

  189. <dl>
    190. 

  190. <dt id="acceleration.barnesHutCPP"><code class="name flex">
    191. 

  191. <span>def <span class="ident">barnesHutCPP</span></span>(<span>r_vec, m_vec, soft=0.0)</span>
    192. 

  192. </code></dt>
    193. 

  193. <dd>
    194. 

  194. <section class="desc"><p>Calculates the acceleration generated by a set of masses on themselves.
    195. 

  195. This is ran in a shared C++ library using a BarnesHut tree</p></section>
    196. 

  196. <details class="source">
    197. 

  197. <summary>Source code</summary>
    198. 

  198. <pre><code class="python">def barnesHutCPP(r_vec, m_vec, soft=0.):
    199. 

  199.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    200. 

  200.     This is ran in a shared C++ library using a BarnesHut tree&#34;&#34;&#34;
    201. 

  201.     # Convert array to data required by C++ library
    202. 

  202.     if ACCLIB is None: loadCPPLib() # Singleton pattern
    203. 

  203.     # Change type to be appropiate for calling library
    204. 

  204.     r_vec_c = (r_vec.ctypes.data + np.arange(r_vec.shape[0]) 
    205. 

  205.         * r_vec.strides[0]).astype(np.uintp)
    206. 

  206.     # Set return type as double*
    207. 

  207.     ACCLIB.barnesHutCPP.restype = np.ctypeslib.ndpointer(dtype=np.float64, 
    208. 

  208.         shape=(r_vec.shape[0]*3,))
    209. 

  209.     # Explicitely pass the corner and size of the box for the top node
    210. 

  210.     px, py, pz = np.min(r_vec, axis=0)
    211. 

  211.     size = np.max(np.max(r_vec, axis=0) - np.min(r_vec, axis=0))
    212. 

  212.     # Call the C++ function: double* barnesHutCPP
    213. 

  213.     accel =  ACCLIB.barnesHutCPP(r_vec_c, m_vec, r_vec.shape[0], 
    214. 

  214.         size, px, py, pz, soft)
    215. 

  215.     # Change shape to get the expected Numpy array (n, 3)
    216. 

  216.     accel.shape = (-1, 3)
    217. 

  217.     return accel</code></pre>
    218. 

  218. </details>
    219. 

  219. </dd>
    220. 

  220. <dt id="acceleration.bruteForce"><code class="name flex">
    221. 

  221. <span>def <span class="ident">bruteForce</span></span>(<span>r_vec, mass, soft=0.0)</span>
    222. 

  222. </code></dt>
    223. 

  223. <dd>
    224. 

  224. <section class="desc"><p>Calculates the acceleration generated by a set of masses on themselves.
    225. 

  225. Complexity O(n*m) where n is the total number of masses and m is the
    226. 

  226. number of massive particles.</p>
    227. 

  227. <h2 id="parameters">Parameters</h2>
    228. 

  228. <dl>
    229. 

  229. <dt><strong><code>r_vec</code></strong> :&ensp;<code>array</code></dt>
    230. 

  230. <dd>list of particles positions.
    231. 

  231. Shape (n, 3) where n is the number of particles</dd>
    232. 

  232. <dt><strong><code>mass</code></strong> :&ensp;<code>array</code></dt>
    233. 

  233. <dd>list of particles masses.
    234. 

  234. Shape (n,)</dd>
    235. 

  235. <dt><strong><code>soft</code></strong> :&ensp;<code>float</code></dt>
    236. 

  236. <dd>characteristic plummer softening length scale</dd>
    237. 

  237. </dl>
    238. 

  238. <h2 id="returns">Returns</h2>
    239. 

  239. <dl>
    240. 

  240. <dt><strong><code>forces</code></strong> :&ensp;<code>array</code></dt>
    241. 

  241. <dd>list of forces acting on each particle.
    242. 

  242. Shape (n, 3)</dd>
    243. 

  243. </dl></section>
    244. 

  244. <details class="source">
    245. 

  245. <summary>Source code</summary>
    246. 

  246. <pre><code class="python">def bruteForce(r_vec, mass, soft=0.):
    247. 

  247.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    248. 

  248.        Complexity O(n*m) where n is the total number of masses and m is the
    249. 

  249.        number of massive particles.
    250. 

  250. 

  251.     Parameters:
    252. 

  252.         r_vec (array): list of particles positions.
    253. 

  253.             Shape (n, 3) where n is the number of particles
    254. 

  254.         mass (array): list of particles masses.
    255. 

  255.             Shape (n,)
    256. 

  256.         soft (float): characteristic plummer softening length scale
    257. 

  257.     Returns:
    258. 

  258.         forces (array): list of forces acting on each particle.
    259. 

  259.             Shape (n, 3)
    260. 

  260.     &#34;&#34;&#34;
    261. 

  261.     # Only calculate forces from massive particles
    262. 

  262.     mask = mass!=0
    263. 

  263.     massMassive = mass[mask]
    264. 

  264.     rMassive_vec = r_vec[mask]
    265. 

  265.     #  x m x 1 matrix (m = number of massive particles) for broadcasting
    266. 

  266.     mass_mat = massMassive.reshape(1, -1, 1)
    267. 

  267.     # Calculate displacements
    268. 

  268.     # r_ten is the direction of the pairwise displacements. Shape (n, m, 3)
    269. 

  269.     # r_mat is the absolute distance of the pairwise displacements. (n, m, 1)
    270. 

  270.     r_ten = rMassive_vec.reshape(1, -1, 3) - r_vec.reshape(-1, 1, 3)
    271. 

  271.     r_mat = np.linalg.norm(r_ten, axis=-1, keepdims=True)
    272. 

  272.     # Avoid division by zeros
    273. 

  273.     # $a = M / (r + \epsilon)^2$, where $\epsilon$ is the softening scale
    274. 

  274.     # r_ten/r_mat gives the direction unit vector
    275. 

  275.     accel = np.divide(r_ten * mass_mat/(r_mat+soft)**2, r_mat, 
    276. 

  276.         where=r_ten.astype(bool), out=r_ten) # Reuse memory from r_ten
    277. 

  277.     return accel.sum(axis=1) # Add all forces on each particle</code></pre>
    278. 

  278. </details>
    279. 

  279. </dd>
    280. 

  280. <dt id="acceleration.bruteForceCPP"><code class="name flex">
    281. 

  281. <span>def <span class="ident">bruteForceCPP</span></span>(<span>r_vec, m_vec, soft=0.0)</span>
    282. 

  282. </code></dt>
    283. 

  283. <dd>
    284. 

  284. <section class="desc"><p>Calculates the acceleration generated by a set of masses on themselves.
    285. 

  285. This is ran in a shared C++ library through Brute Force (pairwise sums)
    286. 

  286. Massive particles must appear first.</p></section>
    287. 

  287. <details class="source">
    288. 

  288. <summary>Source code</summary>
    289. 

  289. <pre><code class="python">def bruteForceCPP(r_vec, m_vec, soft=0.):
    290. 

  290.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    291. 

  291.     This is ran in a shared C++ library through Brute Force (pairwise sums)
    292. 

  292.     Massive particles must appear first.&#34;&#34;&#34;
    293. 

  293.     # Convert array to data required by C++ library
    294. 

  294.     if ACCLIB is None: loadCPPLib() # Singleton pattern
    295. 

  295.     # Change type to be appropiate for calling library
    296. 

  296.     r_vec_c = (r_vec.ctypes.data + np.arange(r_vec.shape[0]) 
    297. 

  297.         * r_vec.strides[0]).astype(np.uintp)
    298. 

  298.     # Set return type as double*
    299. 

  299.     ACCLIB.bruteForceCPP.restype = np.ctypeslib.ndpointer(dtype=np.float64, 
    300. 

  300.         shape=(r_vec.shape[0]*3,))
    301. 

  301.     # Call the C++ function: double* bruteForceCPP
    302. 

  302.     accel =  ACCLIB.bruteForceCPP(r_vec_c, m_vec, r_vec.shape[0], soft)
    303. 

  303.     # Change shape to get the expected Numpy array (n, 3)
    304. 

  304.     accel.shape = (-1, 3)
    305. 

  305.     return accel</code></pre>
    306. 

  306. </details>
    307. 

  307. </dd>
    308. 

  308. <dt id="acceleration.bruteForceNumba"><code class="name flex">
    309. 

  309. <span>def <span class="ident">bruteForceNumba</span></span>(<span>r_vec, mass, soft=0.0)</span>
    310. 

  310. </code></dt>
    311. 

  311. <dd>
    312. 

  312. <section class="desc"><p>Calculates the acceleration generated by a set of masses on themselves.
    313. 

  313. It is done in the same way as in bruteForce, but this
    314. 

  314. method is ran through Numba</p></section>
    315. 

  315. <details class="source">
    316. 

  316. <summary>Source code</summary>
    317. 

  317. <pre><code class="python">@jit(nopython=True) # Numba annotation
    318. 

  318. def bruteForceNumba(r_vec, mass, soft=0.):
    319. 

  319.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    320. 

  320.     It is done in the same way as in bruteForce, but this
    321. 

  321.     method is ran through Numba&#34;&#34;&#34;
    322. 

  322.     mask = mass!=0
    323. 

  323.     massMassive = mass[mask]
    324. 

  324.     rMassive_vec = r_vec[mask]
    325. 

  325.     mass_mat = massMassive.reshape(1, -1, 1)
    326. 

  326.     r_ten = rMassive_vec.reshape(1, -1, 3) - r_vec.reshape(-1, 1, 3)
    327. 

  327.     # Avoid np.linalg.norm to allow Numba optimizations
    328. 

  328.     r_mat = np.sqrt(r_ten[:,:,0:1]**2 + r_ten[:,:,1:2]**2 + r_ten[:,:,2:3]**2)
    329. 

  329.     r_mat = np.where(r_mat == 0, np.ones_like(r_mat), r_mat)
    330. 

  330.     accel = r_ten/r_mat * mass_mat/(r_mat+soft)**2
    331. 

  331.     return accel.sum(axis=1) # Add all forces in each particle</code></pre>
    332. 

  332. </details>
    333. 

  333. </dd>
    334. 

  334. <dt id="acceleration.bruteForceNumbaOptimized"><code class="name flex">
    335. 

  335. <span>def <span class="ident">bruteForceNumbaOptimized</span></span>(<span>r_vec, mass, soft=0.0)</span>
    336. 

  336. </code></dt>
    337. 

  337. <dd>
    338. 

  338. <section class="desc"><p>Calculates the acceleration generated by a set of masses on themselves.
    339. 

  339. This is optimized for high performance with Numba. All massive particles
    340. 

  340. must appear first.</p></section>
    341. 

  341. <details class="source">
    342. 

  342. <summary>Source code</summary>
    343. 

  343. <pre><code class="python">@jit(nopython=True) # Numba annotation
    344. 

  344. def bruteForceNumbaOptimized(r_vec, mass, soft=0.):
    345. 

  345.     &#34;&#34;&#34;Calculates the acceleration generated by a set of masses on themselves.
    346. 

  346.     This is optimized for high performance with Numba. All massive particles
    347. 

  347.     must appear first.&#34;&#34;&#34;
    348. 

  348.     accel = np.zeros_like(r_vec) 
    349. 

  349.     # Use superposition to add all the contributions
    350. 

  350.     n = r_vec.shape[0] # Number of particles
    351. 

  351.     delta = np.zeros((3,)) # Only allocate this once
    352. 

  352.     for i in range(n):
    353. 

  353.         # Only consider pairs with at least one massive particle i
    354. 

  354.         if mass[i] == 0: break
    355. 

  355.         for j in range(i+1, n):
    356. 

  356.             # Explicitely separate components for high performance
    357. 

  357.             # i.e. do not do delta = r_vec[j] - r_vec[i]
    358. 

  358.             # (The effect of this is VERY relevant (x10) and has to do with
    359. 

  359.             # memory reallocation) Numba will vectorize the loops.
    360. 

  360.             for k in range(3): delta[k] = r_vec[j,k] - r_vec[i,k]
    361. 

  361.             r = np.sqrt(delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2])
    362. 

  362.             tripler = (r+soft)**2 * r
    363. 

  363.             
    364. 

  364.             # Compute acceleration on first particle
    365. 

  365.             mr3inv = mass[i]/(tripler)
    366. 

  366.             # Again, do NOT do accel[j] -= mr3inv * delta
    367. 

  367.             for k in range(3): accel[j,k] -= mr3inv * delta[k]
    368. 

  368.             
    369. 

  369.             # Compute acceleration on second particle
    370. 

  370.             # For pairs with one massless particle, no reaction force
    371. 

  371.             if mass[j] == 0: break
    372. 

  372.             # Otherwise, opposite direction (+)
    373. 

  373.             mr3inv = mass[j]/(tripler)
    374. 

  374.             for k in range(3): accel[i,k] += mr3inv * delta[k]
    375. 

  375.     return accel</code></pre>
    376. 

  376. </details>
    377. 

  377. </dd>
    378. 

  378. <dt id="acceleration.loadCPPLib"><code class="name flex">
    379. 

  379. <span>def <span class="ident">loadCPPLib</span></span>(<span>)</span>
    380. 

  380. </code></dt>
    381. 

  381. <dd>
    382. 

  382. <section class="desc"><p>Loads the C++ shared library to the global variable ACCLIB. Must be
    383. 

  383. called before using the library.</p></section>
    384. 

  384. <details class="source">
    385. 

  385. <summary>Source code</summary>
    386. 

  386. <pre><code class="python">def loadCPPLib():
    387. 

  387.     &#34;&#34;&#34;Loads the C++ shared library to the global variable ACCLIB. Must be
    388. 

  388.     called before using the library.&#34;&#34;&#34;
    389. 

  389.     global ACCLIB
    390. 

  390.     ACCLIB = ctypes.CDLL(&#39;cpp/acclib.so&#39;)
    391. 

  391.     # Define appropiate types for library functions
    392. 

  392.     doublepp = np.ctypeslib.ndpointer(dtype=np.uintp) # double**
    393. 

  393.     doublep = ctl.ndpointer(np.float64, flags=&#39;aligned, c_contiguous&#39;)#double*
    394. 

  394.     # Check cpp/acclib.cpp for function signatures
    395. 

  395.     ACCLIB.bruteForceCPP.argtypes = [doublepp, doublep, 
    396. 

  396.         ctypes.c_int, ctypes.c_double]
    397. 

  397.     ACCLIB.barnesHutCPP.argtypes = [doublepp, doublep, 
    398. 

  398.         ctypes.c_int, ctypes.c_double, ctypes.c_double, 
    399. 

  399.         ctypes.c_double, ctypes.c_double, ctypes.c_double]</code></pre>
    400. 

  400. </details>
    401. 

  401. </dd>
    402. 

  402. </dl>
    403. 

  403. </section>
    404. 

  404. <section>
    405. 

  405. </section>
    406. 

  406. </article>
    407. 

  407. <nav id="sidebar">
    408. 

  408. <h1>Index</h1>
    409. 

  409. <div class="toc">
    410. 

  410. <ul></ul>
    411. 

  411. </div>
    412. 

  412. <ul id="index">
    413. 

  413. <li><h3><a href="#header-functions">Functions</a></h3>
    414. 

  414. <ul class="">
    415. 

  415. <li><code><a title="acceleration.barnesHutCPP" href="#acceleration.barnesHutCPP">barnesHutCPP</a></code></li>
    416. 

  416. <li><code><a title="acceleration.bruteForce" href="#acceleration.bruteForce">bruteForce</a></code></li>
    417. 

  417. <li><code><a title="acceleration.bruteForceCPP" href="#acceleration.bruteForceCPP">bruteForceCPP</a></code></li>
    418. 

  418. <li><code><a title="acceleration.bruteForceNumba" href="#acceleration.bruteForceNumba">bruteForceNumba</a></code></li>
    419. 

  419. <li><code><a title="acceleration.bruteForceNumbaOptimized" href="#acceleration.bruteForceNumbaOptimized">bruteForceNumbaOptimized</a></code></li>
    420. 

  420. <li><code><a title="acceleration.loadCPPLib" href="#acceleration.loadCPPLib">loadCPPLib</a></code></li>
    421. 

  421. </ul>
    422. 

  422. </li>
    423. 

  423. </ul>
    424. 

  424. </nav>
    425. 

  425. </main>
    426. 

  426. <footer id="footer">
    427. 

  427. <p>Generated by <a href="https://pdoc3.github.io/pdoc"><cite>pdoc</cite> 0.5.2</a>.</p>
    428. 

  428. </footer>
    429. 

  429. <script src="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.12.0/highlight.min.js"></script>
    430. 

  430. <script>hljs.initHighlightingOnLoad()</script>
    431. 

  431. </body>
    432. 

  432. </html>
    433.