1 /*
2 * Copyright (c) 2025, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
23 * questions.
24 */
25 package experiments;
26
27 import hat.Accelerator;
28 import hat.ComputeContext;
29 import hat.KernelContext;
30 import static hat.KernelContext.*;
31 import hat.NDRange;
32 import hat.backend.Backend;
33 import jdk.incubator.code.Reflect;
34 import optkl.ifacemapper.BoundSchema;
35 import optkl.ifacemapper.Buffer;
36 import optkl.ifacemapper.MappableIface;
37 import optkl.ifacemapper.Schema;
38
39 import java.lang.foreign.MemorySegment;
40 import java.lang.invoke.MethodHandles;
41 import java.util.Random;
42
43 import static optkl.ifacemapper.MappableIface.RO;
44 import static optkl.ifacemapper.MappableIface.RW;
45
46 public class NBodyF32x4 {
47 public interface Universe extends Buffer {
48 long length();
49
50 interface Body extends Struct {
51 float x();
52
53 float y();
54
55 float z();
56
57 float w();
58
59 float vx();
60
61 float vy();
62
63 float vz();
64
65 float vw();
66
67 void x(float x);
68
69 void y(float y);
70
71 void z(float z);
72
73 void w(float z);
74
75 void vx(float vx);
76
77 void vy(float vy);
78
79 void vz(float vz);
80
81 void vw(float vw);
82 }
83
84 Body body(long idx);
85
86 Schema<Universe> schema = Schema.of(Universe.class, resultTable -> resultTable
87 .arrayLen("length")
88 .pad(8)
89 .array("body", array -> array
90 .fields("x", "y", "z", "w", "vx", "vy", "vz", "vw")
91 )
92 );
93
94 static Universe create(Accelerator accelerator, int length) {
95 return BoundSchema.of(accelerator, schema, length).allocate();
96 }
97 }
98
99 @Reflect
100 static public void nbodyKernel(@RO KernelContext kc, @RW Universe universe, float mass, float delT, float espSqr) {
101 float accx = 0.0f;
102 float accy = 0.0f;
103 float accz = 0.0f;
104 Universe.Body body = universe.body(GIX());
105
106 for (int i = 0; i < universe.length(); i++) {
107 Universe.Body otherBody = universe.body(i);
108 float dx = otherBody.x() - body.x();
109 float dy = otherBody.y() - body.y();
110 float dz = otherBody.z() - body.z();
111 float invDist = (float) (1.0f / Math.sqrt(((dx * dx) + (dy * dy) + (dz * dz) + espSqr)));
112 float s = mass * invDist * invDist * invDist;
113 accx = accx + (s * dx);
114 accy = accy + (s * dy);
115 accz = accz + (s * dz);
116 }
117 accx = accx * delT;
118 accy = accy * delT;
119 accz = accz * delT;
120 body.x(body.x() + (body.vx() * delT) + accx * .5f * delT);
121 body.y(body.y() + (body.vy() * delT) + accy * .5f * delT);
122 body.z(body.z() + (body.vz() * delT) + accz * .5f * delT);
123 body.vx(body.vx() + accx);
124 body.vy(body.vy() + accy);
125 body.vz(body.vz() + accz);
126 }
127
128 @Reflect
129 public static void nbodyCompute(@RO ComputeContext cc, @RW Universe universe, final float mass, final float delT, final float espSqr) {
130 var ndrange = NDRange.of1D((int)universe.length());
131 cc.dispatchKernel(ndrange, kernelContext -> nbodyKernel(kernelContext, universe, mass, delT, espSqr));
132 }
133
134 public static void computeSequential(Universe universe, float mass, float delT, float espSqr) {
135
136 var ndrange = NDRange.of1D((int)universe.length());
137 KernelContext kernelContext = new KernelContext(ndrange);
138 //We can't do this once we refactor to static KerneContext
139 throw new RuntimeException("We need NDRANGE for this");
140 // for (GIX() = 0; GIX() < GSX(); GIX()++) {
141 // nbodyKernel(kernelContext,universe,mass,delT,espSqr);
142 // }
143 }
144
145 @Reflect
146 public static void main(String[] args) {
147 final int NUM_BODIES = 1024;
148 var accelerator = new Accelerator(MethodHandles.lookup(), Backend.FIRST);
149 Universe universe = Universe.create(accelerator, NUM_BODIES);
150
151 final float delT = .1f;
152 final float espSqr = 0.1f;
153 final float mass = .5f;
154
155 Random random = new Random(71);
156 for (int bodyIdx = 0; bodyIdx < NUM_BODIES; bodyIdx++) {
157 Universe.Body b = universe.body(bodyIdx);
158
159 final float theta = (float) (Math.random() * Math.PI * 2);
160 final float phi = (float) (Math.random() * Math.PI * 2);
161 final float radius = (float) (Math.random() * 100.f);
162
163 // get random 3D coordinates in sphere
164 b.x((float) (radius * Math.cos(theta) * Math.sin(phi)));
165 b.y((float) (radius * Math.sin(theta) * Math.sin(phi)));
166 b.z((float) (radius * Math.cos(phi)));
167 b.vx(random.nextFloat(1));
168 b.vy(random.nextFloat(1));
169 b.vz(random.nextFloat(1));
170 }
171 Universe universeSeq = Universe.create(accelerator, NUM_BODIES);
172 MemorySegment from = MappableIface.getMemorySegment(universe);
173 MemorySegment toSeq = MappableIface.getMemorySegment(universeSeq);
174 toSeq.copyFrom(from);
175
176 accelerator.compute(computeContext -> nbodyCompute(computeContext, universe, mass, delT, espSqr));
177
178 computeSequential(universeSeq, espSqr, mass, espSqr);
179
180 System.out.println("Delta = "+averageDisplacementError(universe,universeSeq));
181 }
182
183
184
185 /**
186 * Compares two sets of positions and returns the average Euclidean error.
187 * @return The Average Displacement Error (ADE)
188 */
189 public static double averageDisplacementError(Universe lhs, Universe rhs) {
190 double totalError = 0;
191 for (int i = 0; i < lhs.length(); i++) {
192 var rightBody = lhs.body(i);
193 var leftBody = rhs.body(i);
194 double dx = rightBody.x() - leftBody.x();
195 double dy = rightBody.y() - leftBody.y();
196 double dz = rightBody.z() - leftBody.z();
197 totalError += Math.sqrt(dx * dx + dy * dy + dz * dz);
198 }
199 return totalError / lhs.length();
200 }
201
202 }