Compare commits

..

No commits in common. "c07070ae31736eb811463a0874c830b95b609378" and "86a23af45e42960db4ff68bc25e811923a2e8be9" have entirely different histories.

View file

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