1 /*
2 Copyright (c) 2005-2021 Intel Corporation
3
4 Licensed under the Apache License, Version 2.0 (the "License");
5 you may not use this file except in compliance with the License.
6 You may obtain a copy of the License at
7
8 http://www.apache.org/licenses/LICENSE-2.0
9
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
15 */
16
17 #include <cmath>
18 #include <cstring>
19
20 #include "oneapi/tbb/blocked_range.h"
21 #include "oneapi/tbb/parallel_for.h"
22
23 #include "common/gui/video.hpp"
24
25 #ifdef _MSC_VER
26 // warning C4068: unknown pragma
27 #pragma warning(disable : 4068)
28 // warning C4351: new behavior: elements of array 'array' will be default initialized
29 #pragma warning(disable : 4351)
30 #endif
31
32 #include "universe.hpp"
33
34 const colorcomp_t MaterialColor[4][3] = {
35 // BGR
36 { 96, 0, 0 }, // WATER
37 { 0, 48, 48 }, // SANDSTONE
38 { 32, 32, 23 } // SHALE
39 };
40
InitializeUniverse(video const & colorizer)41 void Universe::InitializeUniverse(video const& colorizer) {
42 pulseCounter = pulseTime = 100;
43 pulseX = UniverseWidth / 3;
44 pulseY = UniverseHeight / 4;
45 // Initialize V, S, and T to slightly non-zero values, in order to avoid denormal waves.
46 for (int i = 0; i < UniverseHeight; ++i)
47 #pragma ivdep
48 for (int j = 0; j < UniverseWidth; ++j) {
49 T[i][j] = S[i][j] = V[i][j] = ValueType(1.0E-6);
50 }
51 for (int i = 1; i < UniverseHeight - 1; ++i) {
52 for (int j = 1; j < UniverseWidth - 1; ++j) {
53 float x = float(j - UniverseWidth / 2) / (UniverseWidth / 2);
54 ValueType t = (ValueType)i / (ValueType)UniverseHeight;
55 MaterialType m;
56 D[i][j] = 1.0;
57 // Coefficient values are fictitious, and chosen to visually exaggerate
58 // physical effects such as Rayleigh waves. The fabs/exp line generates
59 // a shale layer with a gentle upwards slope and an anticline.
60 if (t < 0.3f) {
61 m = WATER;
62 M[i][j] = 0.125;
63 L[i][j] = 0.125;
64 }
65 else if (fabs(t - 0.7 + 0.2 * exp(-8 * x * x) + 0.025 * x) <= 0.1) {
66 m = SHALE;
67 M[i][j] = 0.5;
68 L[i][j] = 0.6;
69 }
70 else {
71 m = SANDSTONE;
72 M[i][j] = 0.3;
73 L[i][j] = 0.4;
74 }
75 material[i][j] = m;
76 }
77 }
78 ValueType scale = 2.0f / (ValueType)ColorMapSize;
79 for (int k = 0; k < 4; ++k) {
80 for (int i = 0; i < ColorMapSize; ++i) {
81 colorcomp_t c[3];
82 ValueType t = (i - ColorMapSize / 2) * scale;
83 ValueType r = t > 0 ? t : 0;
84 ValueType b = t < 0 ? -t : 0;
85 ValueType g = 0.5f * fabs(t);
86 memcpy(c, MaterialColor[k], sizeof(c));
87 c[2] = colorcomp_t(r * (255 - c[2]) + c[2]);
88 c[1] = colorcomp_t(g * (255 - c[1]) + c[1]);
89 c[0] = colorcomp_t(b * (255 - c[0]) + c[0]);
90 ColorMap[k][i] = colorizer.get_color(c[2], c[1], c[0]);
91 }
92 }
93 // Set damping coefficients around border to reduce reflections from boundaries.
94 ValueType d = 1.0;
95 for (int k = DamperSize - 1; k > 0; --k) {
96 d *= 1 - 1.0f / (DamperSize * DamperSize);
97 for (int j = 1; j < UniverseWidth - 1; ++j) {
98 D[k][j] *= d;
99 D[UniverseHeight - 1 - k][j] *= d;
100 }
101 for (int i = 1; i < UniverseHeight - 1; ++i) {
102 D[i][k] *= d;
103 D[i][UniverseWidth - 1 - k] *= d;
104 }
105 }
106 drawingMemory = colorizer.get_drawing_memory();
107 }
UpdatePulse()108 void Universe::UpdatePulse() {
109 if (pulseCounter > 0) {
110 ValueType t = (pulseCounter - pulseTime / 2) * 0.05f;
111 V[pulseY][pulseX] += 64 * sqrt(M[pulseY][pulseX]) * exp(-t * t);
112 --pulseCounter;
113 }
114 }
115
116 struct Universe::Rectangle {
117 struct std::pair<int, int> xRange;
118 struct std::pair<int, int> yRange;
RectangleUniverse::Rectangle119 Rectangle(int startX, int startY, int width, int height)
120 : xRange(startX, width),
121 yRange(startY, height) {}
StartXUniverse::Rectangle122 int StartX() const {
123 return xRange.first;
124 }
StartYUniverse::Rectangle125 int StartY() const {
126 return yRange.first;
127 }
WidthUniverse::Rectangle128 int Width() const {
129 return xRange.second;
130 }
HeightUniverse::Rectangle131 int Height() const {
132 return yRange.second;
133 }
EndXUniverse::Rectangle134 int EndX() const {
135 return xRange.first + xRange.second;
136 }
EndYUniverse::Rectangle137 int EndY() const {
138 return yRange.first + yRange.second;
139 }
140 };
141
UpdateStress(Rectangle const & r)142 void Universe::UpdateStress(Rectangle const& r) {
143 drawing_area drawing(r.StartX(), r.StartY(), r.Width(), r.Height(), drawingMemory);
144 for (int i = r.StartY(); i < r.EndY(); ++i) {
145 drawing.set_pos(1, i - r.StartY());
146 #pragma ivdep
147 for (int j = r.StartX(); j < r.EndX(); ++j) {
148 S[i][j] += M[i][j] * (V[i][j + 1] - V[i][j]);
149 T[i][j] += M[i][j] * (V[i + 1][j] - V[i][j]);
150 int index = (int)(V[i][j] * (ColorMapSize / 2)) + ColorMapSize / 2;
151 if (index < 0)
152 index = 0;
153 if (index >= ColorMapSize)
154 index = ColorMapSize - 1;
155 color_t* c = ColorMap[material[i][j]];
156 drawing.put_pixel(c[index]);
157 }
158 }
159 }
160
SerialUpdateStress()161 void Universe::SerialUpdateStress() {
162 Rectangle area(0, 0, UniverseWidth - 1, UniverseHeight - 1);
163 UpdateStress(area);
164 }
165
166 struct UpdateStressBody {
167 Universe& u_;
UpdateStressBodyUpdateStressBody168 UpdateStressBody(Universe& u) : u_(u) {}
operator ()UpdateStressBody169 void operator()(const oneapi::tbb::blocked_range<int>& range) const {
170 Universe::Rectangle area(0, range.begin(), u_.UniverseWidth - 1, range.size());
171 u_.UpdateStress(area);
172 }
173 };
174
ParallelUpdateStress(oneapi::tbb::affinity_partitioner & affinity)175 void Universe::ParallelUpdateStress(oneapi::tbb::affinity_partitioner& affinity) {
176 oneapi::tbb::parallel_for(
177 oneapi::tbb::blocked_range<int>(0, UniverseHeight - 1), // Index space for loop
178 UpdateStressBody(*this), // Body of loop
179 affinity); // Affinity hint
180 }
181
UpdateVelocity(Rectangle const & r)182 void Universe::UpdateVelocity(Rectangle const& r) {
183 for (int i = r.StartY(); i < r.EndY(); ++i)
184 #pragma ivdep
185 for (int j = r.StartX(); j < r.EndX(); ++j)
186 V[i][j] =
187 D[i][j] * (V[i][j] + L[i][j] * (S[i][j] - S[i][j - 1] + T[i][j] - T[i - 1][j]));
188 }
189
SerialUpdateVelocity()190 void Universe::SerialUpdateVelocity() {
191 UpdateVelocity(Rectangle(1, 1, UniverseWidth - 1, UniverseHeight - 1));
192 }
193
194 struct UpdateVelocityBody {
195 Universe& u_;
UpdateVelocityBodyUpdateVelocityBody196 UpdateVelocityBody(Universe& u) : u_(u) {}
operator ()UpdateVelocityBody197 void operator()(const oneapi::tbb::blocked_range<int>& y_range) const {
198 u_.UpdateVelocity(
199 Universe::Rectangle(1, y_range.begin(), u_.UniverseWidth - 1, y_range.size()));
200 }
201 };
202
ParallelUpdateVelocity(oneapi::tbb::affinity_partitioner & affinity)203 void Universe::ParallelUpdateVelocity(oneapi::tbb::affinity_partitioner& affinity) {
204 oneapi::tbb::parallel_for(
205 oneapi::tbb::blocked_range<int>(1, UniverseHeight), // Index space for loop
206 UpdateVelocityBody(*this), // Body of loop
207 affinity); // Affinity hint
208 }
209
SerialUpdateUniverse()210 void Universe::SerialUpdateUniverse() {
211 UpdatePulse();
212 SerialUpdateStress();
213 SerialUpdateVelocity();
214 }
215
ParallelUpdateUniverse()216 void Universe::ParallelUpdateUniverse() {
217 /** Affinity is an argument to parallel_for to hint that an iteration of a loop
218 is best replayed on the same processor for each execution of the loop.
219 It is a static object because it must remember where the iterations happened
220 in previous executions. */
221 static oneapi::tbb::affinity_partitioner affinity;
222 UpdatePulse();
223 ParallelUpdateStress(affinity);
224 ParallelUpdateVelocity(affinity);
225 }
226
TryPutNewPulseSource(int x,int y)227 bool Universe::TryPutNewPulseSource(int x, int y) {
228 if (pulseCounter == 0) {
229 pulseCounter = pulseTime;
230 pulseX = x;
231 pulseY = y;
232 return true;
233 }
234 return false;
235 }
236
SetDrawingMemory(const drawing_memory & dmem)237 void Universe::SetDrawingMemory(const drawing_memory& dmem) {
238 drawingMemory = dmem;
239 }
240