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 hat.test;
26
27 import hat.Accelerator;
28 import hat.ComputeContext;
29 import hat.NDRange;
30
31 import static hat.KernelContext.*;
32 import hat.backend.Backend;
33 import optkl.ifacemapper.BoundSchema;
34 import optkl.ifacemapper.Buffer;
35 import optkl.ifacemapper.Schema;
36 import hat.test.exceptions.HATAssertionError;
37 import hat.test.exceptions.HATExpectedPrecisionError;
38 import jdk.incubator.code.Reflect;
39 import hat.test.annotation.HatTest;
40 import hat.test.exceptions.HATAsserts;
41
42 import java.lang.invoke.MethodHandles;
43 import java.util.Random;
44
45 import static hat.test.TestNbody.Universe.*;
46
47 public class TestNbody {
48
49 public interface Universe extends Buffer {
50 int length();
51
52 interface Body extends Struct {
53 float x();
54 float y();
55 float z();
56 float vx();
57 float vy();
58 float vz();
59 void x(float x);
60 void y(float y);
61 void z(float z);
62 void vx(float vx);
63 void vy(float vy);
64 void vz(float vz);
65 }
66
67 Body body(long idx);
68
69 Schema<Universe> schema = Schema.of(Universe.class, resultTable -> resultTable
70 .arrayLen("length").array("body", array -> array
71 .fields("x", "y", "z", "vx", "vy", "vz")
72 )
73 );
74 static Universe create(Accelerator accelerator, int length) {
75 return BoundSchema.of(accelerator ,schema, length).allocate();
76 }
77 }
78
79 @Reflect
80 static public void nbodyKernel( Universe universe, float mass, float delT, float espSqr) {
81 float accx = 0.0f;
82 float accy = 0.0f;
83 float accz = 0.0f;
84 Body body = universe.body(GIX());
85
86 for (int i = 0; i < universe.length(); i++) {
87 Body otherBody = universe.body(i);
88 float dx = otherBody.x() - body.x();
89 float dy = otherBody.y() - body.y();
90 float dz = otherBody.z() - body.z();
91 float invDist = (float) (1.0f / Math.sqrt(((dx * dx) + (dy * dy) + (dz * dz) + espSqr)));
92 float s = mass * invDist * invDist * invDist;
93 accx = accx + (s * dx);
94 accy = accy + (s * dy);
95 accz = accz + (s * dz);
96 }
97 accx = accx * delT;
98 accy = accy * delT;
99 accz = accz * delT;
100 body.x(body.x() + (body.vx() * delT) + accx * .5f * delT);
101 body.y(body.y() + (body.vy() * delT) + accy * .5f * delT);
102 body.z(body.z() + (body.vz() * delT) + accz * .5f * delT);
103 body.vx(body.vx() + accx);
104 body.vy(body.vy() + accy);
105 body.vz(body.vz() + accz);
106 }
107
108 @Reflect
109 public static void nbodyCompute(ComputeContext cc, Universe universe, final float mass, final float delT, final float espSqr) {
110 cc.dispatchKernel(NDRange.of1D(universe.length()), () -> nbodyKernel( universe, mass, delT, espSqr));
111 }
112
113 public static void computeSequential(Universe universe, float mass, float delT, float espSqr) {
114 float accx = 0.0f;
115 float accy = 0.0f;
116 float accz = 0.0f;
117 for (int j = 0; j < universe.length(); j++) {
118 Body body = universe.body(j);
119
120 for (int i = 0; i < universe.length(); i++) {
121 Body otherBody = universe.body(i);
122 float dx = otherBody.x() - body.x();
123 float dy = otherBody.y() - body.y();
124 float dz = otherBody.z() - body.z();
125 float invDist = (float) (1.0f / Math.sqrt(((dx * dx) + (dy * dy) + (dz * dz) + espSqr)));
126 float s = mass * invDist * invDist * invDist;
127 accx = accx + (s * dx);
128 accy = accy + (s * dy);
129 accz = accz + (s * dz);
130 }
131 accx = accx * delT;
132 accy = accy * delT;
133 accz = accz * delT;
134 body.x(body.x() + (body.vx() * delT) + accx * .5f * delT);
135 body.y(body.y() + (body.vy() * delT) + accy * .5f * delT);
136 body.z(body.z() + (body.vz() * delT) + accz * .5f * delT);
137 body.vx(body.vx() + accx);
138 body.vy(body.vy() + accy);
139 body.vz(body.vz() + accz);
140 }
141 }
142
143 @HatTest
144 @Reflect
145 public void testNbody() {
146 final int NUM_BODIES = 1024;
147 var accelerator = new Accelerator(MethodHandles.lookup(), Backend.FIRST);
148 Universe universe = create(accelerator, NUM_BODIES);
149 Universe universeSeq = create(accelerator, NUM_BODIES);
150 final float delT = .1f;
151 final float espSqr = 0.1f;
152 final float mass = .5f;
153
154 Random random = new Random(71);
155 for (int bodyIdx = 0; bodyIdx < NUM_BODIES; bodyIdx++) {
156 Body b = universe.body(bodyIdx);
157
158 final float theta = (float) (Math.random() * Math.PI * 2);
159 final float phi = (float) (Math.random() * Math.PI * 2);
160 final float radius = (float) (Math.random() * 100.f);
161
162 // get random 3D coordinates in sphere
163 b.x((float) (radius * Math.cos(theta) * Math.sin(phi)));
164 b.y((float) (radius * Math.sin(theta) * Math.sin(phi)));
165 b.z((float) (radius * Math.cos(phi)));
166 b.vx(random.nextFloat(1));
167 b.vy(random.nextFloat(1));
168 b.vz(random.nextFloat(1));
169
170 // Copy random values into the other body to check results
171 Body seqBody = universeSeq.body(bodyIdx);
172 seqBody.x(b.x());
173 seqBody.y(b.y());
174 seqBody.z(b.z());
175
176 seqBody.vx(b.vx());
177 seqBody.vy(b.vy());
178 seqBody.vz(b.vz());
179
180 }
181
182 accelerator.compute(computeContext -> {
183 TestNbody.nbodyCompute(computeContext, universe, mass, delT, espSqr);
184 });
185
186 computeSequential(universeSeq, espSqr, mass, espSqr);
187
188 // Check results
189 float delta = 0.1f;
190 for (int i = 0; i < NUM_BODIES; i++) {
191 Body hatBody = universe.body(i);
192 Body seqBody = universeSeq.body(i);
193 IO.println(i);
194 try {
195 HATAsserts.assertEquals(seqBody.x(), hatBody.x(), delta);
196 HATAsserts.assertEquals(seqBody.y(), hatBody.y(), delta);
197 HATAsserts.assertEquals(seqBody.z(), hatBody.z(), delta);
198 HATAsserts.assertEquals(seqBody.vx(), hatBody.vx(), delta);
199 HATAsserts.assertEquals(seqBody.vy(), hatBody.vy(), delta);
200 HATAsserts.assertEquals(seqBody.vz(), hatBody.vz(), delta);
201 } catch (HATAssertionError hatAssertionError) {
202 throw new HATExpectedPrecisionError(hatAssertionError.getMessage());
203 }
204 }
205 }
206 }