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 
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 }
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;
119     Rectangle(int startX, int startY, int width, int height)
120             : xRange(startX, width),
121               yRange(startY, height) {}
122     int StartX() const {
123         return xRange.first;
124     }
125     int StartY() const {
126         return yRange.first;
127     }
128     int Width() const {
129         return xRange.second;
130     }
131     int Height() const {
132         return yRange.second;
133     }
134     int EndX() const {
135         return xRange.first + xRange.second;
136     }
137     int EndY() const {
138         return yRange.first + yRange.second;
139     }
140 };
141 
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 
161 void Universe::SerialUpdateStress() {
162     Rectangle area(0, 0, UniverseWidth - 1, UniverseHeight - 1);
163     UpdateStress(area);
164 }
165 
166 struct UpdateStressBody {
167     Universe& u_;
168     UpdateStressBody(Universe& u) : u_(u) {}
169     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 
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 
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 
190 void Universe::SerialUpdateVelocity() {
191     UpdateVelocity(Rectangle(1, 1, UniverseWidth - 1, UniverseHeight - 1));
192 }
193 
194 struct UpdateVelocityBody {
195     Universe& u_;
196     UpdateVelocityBody(Universe& u) : u_(u) {}
197     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 
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 
210 void Universe::SerialUpdateUniverse() {
211     UpdatePulse();
212     SerialUpdateStress();
213     SerialUpdateVelocity();
214 }
215 
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 
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 
237 void Universe::SetDrawingMemory(const drawing_memory& dmem) {
238     drawingMemory = dmem;
239 }
240