1
2 /*
3 * Copyright (c) 2024, Oracle and/or its affiliates. All rights reserved.
4 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
5 *
6 * This code is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU General Public License version 2 only, as
8 * published by the Free Software Foundation. Oracle designates this
9 * particular file as subject to the "Classpath" exception as provided
10 * by Oracle in the LICENSE file that accompanied this code.
11 *
12 * This code is distributed in the hope that it will be useful, but WITHOUT
13 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 * version 2 for more details (a copy is included in the LICENSE file that
16 * accompanied this code).
17 *
18 * You should have received a copy of the GNU General Public License version
19 * 2 along with this work; if not, write to the Free Software Foundation,
20 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
21 *
22 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
23 * or visit www.oracle.com if you need additional information or have any
24 * questions.
25 */
26 package nbody;
27
28 import hat.Accelerator;
29 import hat.Accelerator.Compute;
30 import hat.ComputeContext;
31
32 import static hat.KernelContext.*;
33 import hat.NDRange;
34 import hat.buffer.F32Array;
35 import hat.buffer.S32RGBAImage;
36
37 import jdk.incubator.code.Reflect;
38
39 import javax.swing.JFrame;
40 import javax.swing.JPanel;
41 import java.awt.Dimension;
42 import java.awt.Graphics;
43 import java.awt.Graphics2D;
44 import java.awt.RenderingHints;
45 import java.awt.image.BufferedImage;
46 import java.awt.image.DataBufferInt;
47 import java.awt.image.WritableRaster;
48 import java.lang.invoke.MethodHandles;
49 import java.util.stream.IntStream;
50
51 import hat.types.vec3;
52 import optkl.ifacemapper.MappableIface;
53
54 import static hat.types.F32.inversesqrt;
55 import static hat.types.F32.pow;
56 import static hat.types.vec3.dot;
57 import static hat.types.vec3.mul;
58 import static hat.types.vec3.sub;
59 import static hat.types.vec3.vec3;
60
61 public class Main extends JFrame implements Runnable {
62
63 public enum Mode {HAT, JavaMt, JavaSeq}
64
65 public static class DirectRasterPanel extends JPanel {
66 private BufferedImage bufferedImage;
67 private WritableRaster writableRaster;
68 private DataBufferInt dataBuffer;
69 private S32RGBAImage image;
70
71 public DirectRasterPanel(Accelerator acc, int width, int height) {
72 setPreferredSize(new Dimension(width, height));
73 this.image = S32RGBAImage.create(acc, width, height);
74 this.bufferedImage = new BufferedImage(width, height, BufferedImage.TYPE_INT_RGB);
75 this.writableRaster = bufferedImage.getRaster();
76 this.dataBuffer = ((DataBufferInt) writableRaster.getDataBuffer());
77 }
78
79 @Override
80 protected void paintComponent(Graphics g) {
81 Graphics2D g2d = (Graphics2D) g;
82 g2d.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);
83 g2d.drawImage(bufferedImage, 0, 0, null);
84 }
85 }
86
87 final Accelerator acc;
88 final Mode mode;
89
90 final float delT = .9f;
91 final int width;
92 final int height;
93 final float espSqr = 0.9f;
94
95 final float mass = .9f;
96 final int bodies;
97
98 final F32Array xyzPosFloatArr;
99 final F32Array xyzVelFloatArr;
100 DirectRasterPanel panel;
101
102 public Main(Accelerator acc, Mode mode, int bodyCount, int width, int height) {
103 super("NBody Sans OpenGL");
104 this.mode = mode;
105 this.acc = acc;
106 this.bodies = bodyCount;
107 this.width = width;
108 this.height = height;
109 this.xyzPosFloatArr = F32Array.create(acc, bodies * 4);
110 this.xyzVelFloatArr = F32Array.create(acc, bodies * 4);
111 panel = new DirectRasterPanel(acc, width, height);
112 final float maxDist = width / 2;
113
114 //System.out.println(bodies + " particles");
115
116 for (int body = 0; body < bodies; body++) {
117 final float theta = (float) (Math.random() * Math.PI * 2);
118 final float phi = (float) (Math.random() * Math.PI * 2);
119 final float radius = (float) (Math.random() * maxDist);
120
121 var radialx =
122 (float) (radius * Math.cos(theta) * Math.sin(phi) + width / 2);
123 var radialy =
124 (float) (radius * Math.sin(theta) * Math.sin(phi) + height / 2);
125 var radialz =
126 (float) (radius * Math.cos(phi) + Math.min(width, height) / 2);
127
128 xyzPosFloatArr.array(body * 4 + 0, radialx);
129 xyzPosFloatArr.array(body * 4 + 1, radialy);
130 xyzPosFloatArr.array(body * 4 + 2, radialz);
131 }
132
133 add(panel);
134 pack();
135 setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
136 setVisible(true);
137 new Thread(this).start();
138 }
139
140 @Reflect
141 public static void runSansVec(int bodyIdx, int bodies, F32Array xyzPos, F32Array xyzVel, S32RGBAImage image, int imageWidth, float mass, float delT, float espSqr) {
142 final int FAR = 500;
143 final int MID = 300;
144 final int NEAR = 100;
145 float accx = 0.f;
146 float accy = 0.f;
147 float accz = 0.f;
148
149 final float myPosx = xyzPos.array(bodyIdx * 4 + 0);
150 final float myPosy = xyzPos.array(bodyIdx * 4 + 1);
151 final float myPosz = xyzPos.array(bodyIdx * 4 + 2);
152 final float myVelx = xyzVel.array(bodyIdx * 4 + 0);
153 final float myVely = xyzVel.array(bodyIdx * 4 + 1);
154 final float myVelz = xyzVel.array(bodyIdx * 4 + 2);
155
156 for (int body = 0; body < bodies; body++) {
157 final float dx = xyzPos.array(body * 4 + 0) - myPosx;
158 final float dy = xyzPos.array(body * 4 + 1) - myPosy;
159 final float dz = xyzPos.array(body * 4 + 2) - myPosz;
160 final float invDist = 1 / (float) Math.sqrt((dx * dx) + (dy * dy) + (dz * dz) + espSqr);
161 final float s = mass * invDist * invDist * invDist;
162 accx = accx + (s * dx);
163 accy = accy + (s * dy);
164 accz = accz + (s * dz);
165 }
166 accx = accx * delT;
167 accy = accy * delT;
168 accz = accz * delT;
169
170 float fx = myPosx + (myVelx + accx * .5f) * delT;
171 float fy = myPosy + (myVely + accx * .5f) * delT;
172 float fz = myPosz + (myVelz + accx * .5f) * delT;
173 xyzPos.array(bodyIdx * 4 + 0, fx);
174 xyzPos.array(bodyIdx * 4 + 1, fy);
175 xyzPos.array(bodyIdx * 4 + 2, fz);
176
177 xyzVel.array(bodyIdx * 4 + 0, myVelx + accx);
178 xyzVel.array(bodyIdx * 4 + 1, myVely + accy);
179 xyzVel.array(bodyIdx * 4 + 2, myVelz + accz);
180
181 int x = (int) fx;
182 int y = (int) fy;
183 int z = (int) fz;
184 if (x > 1 && x < imageWidth - 2 && y > 1 && y < imageWidth - 2) {
185 // Calculate brightness based on depth (Z)
186 int brightness = (255 - (z / imageWidth * 255));
187 int color = (brightness << 16) | (brightness << 8) | brightness;
188 int pos = ((y * imageWidth) + x);
189 image.data(pos, color);
190 if (z < FAR) {
191 image.data(pos + 1, color);
192 image.data(pos - 1, color);
193 image.data(pos + imageWidth, color);
194 image.data(pos - imageWidth, color);
195 if (z < MID) {
196 image.data(pos + imageWidth + 1, color);
197 image.data(pos + imageWidth - 1, color);
198 image.data(pos - imageWidth + 1, color);
199 image.data(pos - imageWidth - 1, color);
200 if (z < NEAR) {
201 image.data(pos + imageWidth * 2 + 2, color);
202 image.data(pos + imageWidth * 2 - 2, color);
203 image.data(pos - imageWidth * 2 + 2, color);
204 image.data(pos - imageWidth * 2 - 2, color);
205 }
206 }
207 }
208 }
209 }
210
211 @Reflect
212 public static void run(int bodyIdx, int bodies, F32Array xyzPos, F32Array xyzVel, S32RGBAImage image, int imageWidth, float mass, float delT, float espSqr) {
213 final int FAR = 500;
214 final int MID = 300;
215 final int NEAR = 100;
216 var acc = vec3(0f);
217 var idx = bodyIdx * 4;
218 var myPos = vec3(xyzPos.array(idx),xyzPos.array(idx+1),xyzPos.array(idx+2));
219 var myVel = vec3(xyzVel.array(idx),xyzVel.array(idx+1),xyzVel.array(idx+2));
220 for (int i = 0; i < bodies*4; i+=4) {
221 var delta = sub(vec3(xyzPos.array(i),xyzPos.array(i+1),xyzPos.array(i+2)),myPos);
222 float invDist = inversesqrt(dot(delta,delta)+espSqr);
223 acc = vec3.add(acc,mul(delta,pow(invDist,3f)*mass));
224 }
225 acc = mul(acc,delT);
226 var f = vec3.add(myPos,mul(vec3.add(myVel,mul(acc,.5f)),delT));
227
228 xyzPos.array(idx, f.x());
229 xyzPos.array(idx+1, f.y());
230 xyzPos.array(idx+2, f.z());
231
232 xyzVel.array(idx, myVel.x() + acc.x());
233 xyzVel.array(idx+1, myVel.y() + acc.y());
234 xyzVel.array(idx+2, myVel.z() + acc.z());
235
236 int x = (int) f.x();
237 int y = (int) f.y();
238 int z = (int) f.z();
239 if (x > 1 && x < imageWidth - 2 && y > 1 && y < imageWidth - 2) {
240 // Calculate brightness based on depth (Z)
241 int brightness = (255 - (z / imageWidth * 255));
242 int color = (brightness << 16) | (brightness << 8) | brightness;
243 int pos = ((y * imageWidth) + x);
244 image.data(pos, color);
245 if (z < FAR) {
246 image.data(pos + 1, color);
247 image.data(pos - 1, color);
248 image.data(pos + imageWidth, color);
249 image.data(pos - imageWidth, color);
250 if (z < MID) {
251 image.data(pos + imageWidth + 1, color);
252 image.data(pos + imageWidth - 1, color);
253 image.data(pos - imageWidth + 1, color);
254 image.data(pos - imageWidth - 1, color);
255 if (z < NEAR) {
256 image.data(pos + imageWidth * 2 + 2, color);
257 image.data(pos + imageWidth * 2 - 2, color);
258 image.data(pos - imageWidth * 2 + 2, color);
259 image.data(pos - imageWidth * 2 - 2, color);
260 }
261 }
262 }
263 }
264 }
265
266 @Reflect
267
268 static public void nbodyKernel(
269 int bodies,
270 F32Array xyzPos,
271 F32Array xyzVel,
272 S32RGBAImage image,
273 int imageWidth,
274 float mass,
275 float delT,
276 float espSqr
277 ) {
278 run(GIX(), bodies, xyzPos, xyzVel, image, imageWidth, mass, delT, espSqr);
279 }
280
281 @Reflect
282 public static void clearImage(S32RGBAImage image) {
283 image.data(GIX(), 0);
284 }
285
286 @Reflect
287 public static void nbodyCompute(
288 ComputeContext cc,
289 int bodies,
290 F32Array xyzPos,
291 F32Array xyzVel,
292 S32RGBAImage image,
293 int imageWidth,
294 float mass,
295 float delT,
296 float espSqr) {
297 float cmass = mass;
298 float cdelT = delT;
299 float cespSqr = espSqr;
300 int cbodies = bodies;
301 int cimageWidth = imageWidth;
302
303 cc.dispatchKernel(NDRange.of1D(imageWidth * image.height()), () -> clearImage( image));
304
305 cc.dispatchKernel(NDRange.of1D(bodies), () -> nbodyKernel( cbodies, xyzPos, xyzVel, image, cimageWidth, cmass, cdelT, cespSqr));
306 }
307
308
309 public void run() {
310 while (true) {
311 // long startNs = System.nanoTime();
312 float cmass = mass;
313 float cdelT = delT;
314 float cespSqr = espSqr;
315 int cbodies = bodies;
316 int cimageWidth = width;
317 F32Array cxyzPosFloatArr = xyzPosFloatArr;
318 F32Array cxyzVelFloatArr = xyzVelFloatArr;
319 S32RGBAImage cimage = panel.image;
320
321 switch (mode) {
322 case HAT -> {
323 acc.compute((@Reflect Compute)
324 cc -> nbodyCompute(cc, cbodies, cxyzPosFloatArr, cxyzVelFloatArr, cimage, cimageWidth, cmass, cdelT, cespSqr));
325 }
326 case JavaMt -> {
327 MappableIface.getMemorySegment(panel.image).fill((byte) 0x00); // Dont do this if using HAT! ;)
328 IntStream.range(0, bodies).parallel().forEach(
329 i -> run(i, cbodies, cxyzPosFloatArr, cxyzVelFloatArr, cimage, cimageWidth, cmass, cdelT, cespSqr));
330 }
331 case JavaSeq -> {
332 MappableIface.getMemorySegment(panel.image).fill((byte) 0x00); // Dont do this if using HAT! ;)
333 IntStream.range(0, bodies).forEach(
334 i -> run(i, cbodies, cxyzPosFloatArr, cxyzVelFloatArr, cimage, cimageWidth, cmass, cdelT, cespSqr));
335 }
336 }
337 panel.image.syncToRasterDataBuffer(panel.dataBuffer);
338
339 // long endNs = System.nanoTime();
340 // System.out.println((endNs - startNs) / 1000000 + "ms");
341 repaint();
342 try {
343 Thread.sleep(1);
344 } catch (Exception e) {
345 }
346 }
347 }
348
349 static void main(String[] args) {
350 new Main(new Accelerator(MethodHandles.lookup()),
351 args.length > 0 ?
352 switch (args[0]) {
353 case "HAT" -> Mode.HAT;
354 case "JavaMT" -> Mode.JavaMt;
355 case "JavaSeq" -> Mode.JavaSeq;
356 default -> throw new IllegalStateException("No such mode as " + args[0]);
357 } : Mode.HAT
358 , 4096 * 4, 1024, 1024);
359
360 }
361 }
362