Peano 4
Loading...
Searching...
No Matches
AnalyticalSolution.cpp
Go to the documentation of this file.
1// solves the exact riemann problem
2// original code from bruce fryxell: https://cococubed.com/code_pages/exact_riemann.shtml
3
4#include <iostream>
5#include <vector>
6#include <string>
7#include <cmath>
8#include <fstream>
9#include <limits>
10
11// solves shock tube equation
12double f(double p4, double p1, double p5, double rho1, double rho5, double gamma)
13{
14 // local variables
15 double z, c1, c5, gm1, gp1, g2, fact;
16
17 z = (p4 / p5 - 1.0);
18 c1 = std::sqrt(gamma * p1 / rho1);
19 c5 = std::sqrt(gamma * p5 / rho5);
20
21 gm1 = gamma - 1.0;
22 gp1 = gamma + 1.0;
23 g2 = 2.0 * gamma;
24
25 fact = gm1 / g2 * (c5 / c1) * z / std::sqrt(1.0 + gp1 / g2 * z);
26 fact = std::pow(1.0 - fact, g2 / gm1);
27
28 return p1 * fact - p4;
29}
30
31int main(int argn, char **argc)
32{
33 if (argn > 2)
34 {
35 int npts = std::stoi(argc[1]); // number of cells per row/column
36 double t = std::stod(argc[2]); // time at which solution is desired
37
38 // declare variables
39 int itmax, iter, i;
40 std::vector<double> x(npts), rho(npts), u(npts), p(npts);
41 double rhol, pl, ul, rhor, pr, ur, gamma, xi, xl, xr,
42 rho1, p1, u1, rho5, p5, u5, p40, p41, f0, eps, f1, p4,
43 error, z, c5, gm1, gp1, gmfac1, gmfac2, fact, u4, rho4, w,
44 p3, u3, rho3, c1, c3, xsh, xcd, xft, xhd, dx;
45
46 // set initial conditions
47
48 // state at left of discontinuity
49 rhol = 1.0;
50 pl = 1.0;
51 ul = 0.0;
52
53 // state at right of discontinuity
54 rhor = 0.125;
55 pr = 0.1;
56 ur = 0.0;
57
58 // equation of state
59 gamma = 1.4;
60
61 // location of discontinuity at t = 0
62 xi = 0.5;
63
64 // spatial interval over which to compute solution
65 xl = 0.0;
66 xr = 1.0;
67
68 // begin solution
69 rho1 = rhol;
70 p1 = pl;
71 u1 = ul;
72 rho5 = rhor;
73 p5 = pr;
74 u5 = ur;
75
76 // solve for post-shock pressure by secant method
77 // initial guesses
78
79 p40 = p1;
80 p41 = p5;
81 f0 = f(p40, p1, p5, rho1, rho5, gamma);
82
83 // maximum number of iterations and maxium allowable relative error
84 itmax = 200;
85 eps = std::numeric_limits<double>::min();
86 iter = 0;
87 while (iter < itmax)
88 {
89 f1 = f(p41, p1, p5, rho1, rho5, gamma);
90 if (f1 == f0)
91 break;
92 p4 = p41 - (p41 - p40) * f1 / (f1 - f0);
93 error = std::abs(p4 - p41) / p41;
94 if (error < eps)
95 break;
96 p40 = p41;
97 p41 = p4;
98 f0 = f1;
99 iter++;
100 }
101
102 // compute post-shock density and velocity
103 z = (p4 / p5 - 1.0);
104 c5 = std::sqrt(gamma * p5 / rho5);
105
106 gm1 = gamma - 1.0;
107 gp1 = gamma + 1.0;
108 gmfac1 = 0.5 * gm1 / gamma;
109 gmfac2 = 0.5 * gp1 / gamma;
110
111 fact = std::sqrt(1.0 + gmfac2 * z);
112
113 u4 = c5 * z / (gamma * fact);
114 rho4 = rho5 * (1.0 + gmfac2 * z) / (1.0 + gmfac1 * z);
115
116 // shock speed
117 w = c5 * fact;
118
119 // compute values at foot of rarefaction
120 p3 = p4;
121 u3 = u4;
122 rho3 = rho1 * std::pow(p3 / p1, 1.0 / gamma);
123
124 // compute positions of waves
125 c1 = sqrt(gamma * p1 / rho1);
126 c3 = sqrt(gamma * p3 / rho3);
127
128 xsh = xi + w * t;
129 xcd = xi + u3 * t;
130 xft = xi + (u3 - c3) * t;
131 xhd = xi - c1 * t;
132
133 // compute solution as a function of position
134 dx = (xr - xl) / (npts - 1);
135 for (int i = 0; i < npts; i++)
136 x[i] = xl + dx * (i + 0.5);
137
138 for (int i = 0; i < npts; i++)
139 {
140 if (x[i] < xhd)
141 {
142 rho[i] = rho1;
143 p[i] = p1;
144 u[i] = u1;
145 }
146 else if (x[i] < xft)
147 {
148 u[i] = 2.0 / gp1 * (c1 + (x[i] - xi) / t);
149 fact = 1.0 - 0.5 * gm1 * u[i] / c1;
150 rho[i] = rho1 * std::pow(fact, (2.0 / gm1));
151 p[i] = p1 * std::pow(fact, 2.0 * gamma / gm1);
152 }
153 else if (x[i] < xcd)
154 {
155 rho[i] = rho3;
156 p[i] = p3;
157 u[i] = u3;
158 }
159 else if (x[i] < xsh)
160 {
161 rho[i] = rho4;
162 p[i] = p4;
163 u[i] = u4;
164 }
165 else
166 {
167 rho[i] = rho5;
168 p[i] = p5;
169 u[i] = u5;
170 }
171 }
172
173 std::ofstream output("sod_analytical.csv");
174 output << "x-momentum, density, energy\n";
175 for (int i = 0; i < npts; i++)
176 {
177 double e = p[i] / (gamma - 1) + 0.5 * u[i] * u[i];
178 output << rho[i] * u[i] << "," << rho[i] << "," << e << "\n";
179 }
180 output.close();
181 }
182 else
183 {
184 std::cout << "Error provide following details:\n1)number of cells per row/column \n2)time of solution\n";
185 }
186}
int main(int argn, char **argc)
Main routine of the SPH code.
double f(double p4, double p1, double p5, double rho1, double rho5, double gamma)
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx