Compare commits

...

4 commits

Author SHA1 Message Date
c07070ae31 Change element order in data structure to allow vectorization
Signed-off-by: Christoph Niethammer <niethammer@hlrs.de>
2024-03-13 13:49:45 +01:00
d6a0d64176 Change order of loops
Try to avoid branching in loops as it prevents vectorization

Signed-off-by: Christoph Niethammer <niethammer@hlrs.de>
2024-03-13 13:46:44 +01:00
39bee2e2e6 Reduce number of timesteps for testing to 1000
Signed-off-by: Christoph Niethammer <niethammer@hlrs.de>
2024-03-13 13:26:07 +01:00
3fbfd989b2 Output simulation parameters to console at start
Signed-off-by: Christoph Niethammer <niethammer@hlrs.de>
2024-03-13 12:49:41 +01:00

View file

@ -28,12 +28,23 @@ int main()
* *
* Changing of any of them might result in NaN error. * Changing of any of them might result in NaN error.
*/ */
int dt = 100000; int dt = 1000;
double h = 0.007; double h = 0.007;
std::cout << "--------------------------------------------------" << std::endl;
std::cout << "Simulation parameters" << std::endl;
std::cout << "--------------------------------------------------" << std::endl;
std::cout << "l: " << l << std::endl;
std::cout << "dx: " << dx << std::endl;
std::cout << "kmax: " << kmax << std::endl;
std::cout << "dt: " << dt << std::endl;
std::cout << "h: " << h << std::endl;
std::cout << "--------------------------------------------------" << std::endl;
// Stores the approximated values for each step // Stores the approximated values for each step
double Q[kmax][dt]; double Q[dt][kmax];
double P[kmax][dt]; double P[dt][kmax];
// Prepare initial condition // Prepare initial condition
for (int i = 0; i < dt; i++) for (int i = 0; i < dt; i++)
@ -44,12 +55,12 @@ int main()
// Boundary condition for the last spatial step // Boundary condition for the last spatial step
if (k == kmax - 1) if (k == kmax - 1)
{ {
Q[k][i] = 0; Q[i][k] = 0;
P[k][i] = -202.0; // P = -202 when Q = 10 P[i][k] = -202.0; // P = -202 when Q = 10
} else } else
{ {
Q[k][i] = 0; Q[i][k] = 0;
P[k][i] = 0; P[i][k] = 0;
} }
} }
} }
@ -61,14 +72,13 @@ int main()
*/ */
for (int i = 2; i < dt; i++) for (int i = 2; i < dt; i++)
{ {
// Loop over spatial steps // Loop over spatial steps for flow
for (int k = 1; k < kmax; k++) for (int k = 1; k < kmax; k++) {
{ Q[i + 1][k] = Q[i][k] + h * (alpha * (P[i][k - 1] - P[i][k]) - beta * Q[i][k] * abs(Q[i][k]));
Q[k][i + 1] = Q[k][i] + h * (alpha * (P[k - 1][i] - P[k][i]) - beta * Q[k][i] * abs(Q[k][i])); }
// Loop over steps for pressure, DO NOT calculate pressure for the last element as it is given as a boundary condition
// DO NOT calculate pressure for the last element such as it is given by default for (int k = 1; k < kmax - 1; k++) {
if (k != kmax - 1) P[i + 1][k] = P[k][i - 2] + h * (gamma * (Q[i][k] - Q[i][k + 1]));
P[k][i + 1] = P[k][i - 2] + h * (gamma * (Q[k][i] - Q[k + 1][i]));
} }
} }
@ -81,7 +91,7 @@ int main()
std::cout << i << ": "; std::cout << i << ": ";
for(int k = 0; k < 10; k++) for(int k = 0; k < 10; k++)
std::cout << std::setw(10) << Q[i][k]; std::cout << std::setw(10) << Q[k][i];
std::cout << std::endl; std::cout << std::endl;
} }
@ -94,7 +104,7 @@ int main()
std::cout << i << ": "; std::cout << i << ": ";
for(int k = 0; k < 10; k++) for(int k = 0; k < 10; k++)
std::cout << std::setw(10) << Q[i][dt - 10 + k]; std::cout << std::setw(10) << Q[dt - 10 + k][i];
std::cout << std::endl; std::cout << std::endl;
} }