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 hat.NDRange;
31 import hat.backend.Backend;
32 import jdk.incubator.code.Reflect;
33 import optkl.ifacemapper.BoundSchema;
34 import optkl.ifacemapper.Buffer;
35 import optkl.ifacemapper.MappableIface;
36 import optkl.ifacemapper.Schema;
37
38 import java.lang.foreign.MemorySegment;
39 import java.lang.invoke.MethodHandles;
40 import java.util.Random;
41
42 import static optkl.ifacemapper.MappableIface.RO;
43 import static optkl.ifacemapper.MappableIface.RW;
44
45 public class NBodyF32x4 {
46 public interface Universe extends Buffer {
47 long length();
48
49 interface Body extends Struct {
50 float x();
51
52 float y();
53
54 float z();
55
56 float w();
57
58 float vx();
59
60 float vy();
61
62 float vz();
63
64 float vw();
65
66 void x(float x);
67
68 void y(float y);
69
70 void z(float z);
71
72 void w(float z);
73
74 void vx(float vx);
75
76 void vy(float vy);
77
78 void vz(float vz);
79
80 void vw(float vw);
81 }
82
83 Body body(long idx);
84
85 Schema<Universe> schema = Schema.of(Universe.class, resultTable -> resultTable
86 .arrayLen("length")
87 .pad(8)
88 .array("body", array -> array
89 .fields("x", "y", "z", "w", "vx", "vy", "vz", "vw")
90 )
91 );
92
93 static Universe create(Accelerator accelerator, int length) {
94 return BoundSchema.of(accelerator, schema, length).allocate();
95 }
96 }
97
98 @Reflect
99 static public void nbodyKernel(@RO KernelContext kc, @RW Universe universe, float mass, float delT, float espSqr) {
100 float accx = 0.0f;
101 float accy = 0.0f;
102 float accz = 0.0f;
103 Universe.Body body = universe.body(kc.gix);
104
105 for (int i = 0; i < universe.length(); i++) {
106 Universe.Body otherBody = universe.body(i);
107 float dx = otherBody.x() - body.x();
108 float dy = otherBody.y() - body.y();
109 float dz = otherBody.z() - body.z();
110 float invDist = (float) (1.0f / Math.sqrt(((dx * dx) + (dy * dy) + (dz * dz) + espSqr)));
111 float s = mass * invDist * invDist * invDist;
112 accx = accx + (s * dx);
113 accy = accy + (s * dy);
114 accz = accz + (s * dz);
115 }
116 accx = accx * delT;
117 accy = accy * delT;
118 accz = accz * delT;
119 body.x(body.x() + (body.vx() * delT) + accx * .5f * delT);
120 body.y(body.y() + (body.vy() * delT) + accy * .5f * delT);
121 body.z(body.z() + (body.vz() * delT) + accz * .5f * delT);
122 body.vx(body.vx() + accx);
123 body.vy(body.vy() + accy);
124 body.vz(body.vz() + accz);
125 }
126
127 @Reflect
128 public static void nbodyCompute(@RO ComputeContext cc, @RW Universe universe, final float mass, final float delT, final float espSqr) {
129 var ndrange = NDRange.of1D((int)universe.length());
130 cc.dispatchKernel(ndrange, kernelContext -> nbodyKernel(kernelContext, universe, mass, delT, espSqr));
131 }
132
133 public static void computeSequential(Universe universe, float mass, float delT, float espSqr) {
134 var ndrange = NDRange.of1D((int)universe.length());
135 KernelContext kernelContext = new KernelContext(ndrange);
136 for (kernelContext.gix = 0; kernelContext.gix < kernelContext.gsx; kernelContext.gix++) {
137 nbodyKernel(kernelContext,universe,mass,delT,espSqr);
138 }
139 }
140
141 @Reflect
142 public static void main(String[] args) {
143 final int NUM_BODIES = 1024;
144 var accelerator = new Accelerator(MethodHandles.lookup(), Backend.FIRST);
145 Universe universe = Universe.create(accelerator, NUM_BODIES);
146
147 final float delT = .1f;
148 final float espSqr = 0.1f;
149 final float mass = .5f;
150
151 Random random = new Random(71);
152 for (int bodyIdx = 0; bodyIdx < NUM_BODIES; bodyIdx++) {
153 Universe.Body b = universe.body(bodyIdx);
154
155 final float theta = (float) (Math.random() * Math.PI * 2);
156 final float phi = (float) (Math.random() * Math.PI * 2);
157 final float radius = (float) (Math.random() * 100.f);
158
159 // get random 3D coordinates in sphere
160 b.x((float) (radius * Math.cos(theta) * Math.sin(phi)));
161 b.y((float) (radius * Math.sin(theta) * Math.sin(phi)));
162 b.z((float) (radius * Math.cos(phi)));
163 b.vx(random.nextFloat(1));
164 b.vy(random.nextFloat(1));
165 b.vz(random.nextFloat(1));
166 }
167 Universe universeSeq = Universe.create(accelerator, NUM_BODIES);
168 MemorySegment from = MappableIface.getMemorySegment(universe);
169 MemorySegment toSeq = MappableIface.getMemorySegment(universeSeq);
170 toSeq.copyFrom(from);
171
172 accelerator.compute(computeContext -> nbodyCompute(computeContext, universe, mass, delT, espSqr));
173
174 computeSequential(universeSeq, espSqr, mass, espSqr);
175
176 System.out.println("Delta = "+averageDisplacementError(universe,universeSeq));
177 }
178
179
180
181 /**
182 * Compares two sets of positions and returns the average Euclidean error.
183 * @return The Average Displacement Error (ADE)
184 */
185 public static double averageDisplacementError(Universe lhs, Universe rhs) {
186 double totalError = 0;
187 for (int i = 0; i < lhs.length(); i++) {
188 var rightBody = lhs.body(i);
189 var leftBody = rhs.body(i);
190 double dx = rightBody.x() - leftBody.x();
191 double dy = rightBody.y() - leftBody.y();
192 double dz = rightBody.z() - leftBody.z();
193 totalError += Math.sqrt(dx * dx + dy * dy + dz * dz);
194 }
195 return totalError / lhs.length();
196 }
197
198 }