3
3
#include < cmath>
4
4
#include < iomanip>
5
5
#include < iostream>
6
- #include < vector>
7
6
#include < numeric>
7
+ #include < vector>
8
8
9
9
using std::sin;
10
10
using std::sqrt;
11
11
12
12
// Simplifies strided access on a 1d buffer
13
- template <class T >
13
+ template <class T >
14
14
class StridedAccess {
15
- public:
16
- StridedAccess (T*first, int stride) : first(first), stride(stride) {};
17
- // This accessor allows for row-major access
18
- T& operator ()(int i, int j) { return first[j*stride+i]; }
19
- private:
20
- T*first;
21
- int stride;
15
+ public:
16
+ StridedAccess (T *first, int stride)
17
+ : first(first), stride(stride){};
18
+ // This accessor allows for row-major access
19
+ T &operator ()(int i, int j)
20
+ {
21
+ return first[j * stride + i];
22
+ }
23
+
24
+ private:
25
+ T * first;
26
+ int stride;
22
27
};
23
28
24
29
extern " C" {
@@ -35,16 +40,16 @@ void dgesv_(
35
40
36
41
/* Function for fluid_nl i.e. non-linear */
37
42
int fluidComputeSolutionSerial (
38
- double const * const velocity_old,
39
- double const * const pressure_old,
40
- double const * const crossSectionLength_old,
41
- double const * const crossSectionLength,
42
- double t,
43
- int N,
44
- double kappa,
45
- double tau,
46
- double * velocity,
47
- double * pressure)
43
+ double const *const velocity_old,
44
+ double const *const pressure_old,
45
+ double const *const crossSectionLength_old,
46
+ double const *const crossSectionLength,
47
+ double t,
48
+ int N,
49
+ double kappa,
50
+ double tau,
51
+ double * velocity,
52
+ double * pressure)
48
53
{
49
54
const double PI = 3.141592653589793 ;
50
55
@@ -59,67 +64,52 @@ int fluidComputeSolutionSerial(
59
64
60
65
// Used as Ax = b
61
66
// i.e. LHS*x = Res
62
- std::vector<double > Res (2 * N + 2 );
63
- std::vector<double > LHS_buffer (std::pow (2 * N + 2 , 2 ));
64
- StridedAccess<double > LHS (LHS_buffer.data (), 2 * N + 2 );
67
+ std::vector<double > Res (2 * N + 2 );
68
+ std::vector<double > LHS_buffer (std::pow (2 * N + 2 , 2 ));
69
+ StridedAccess<double > LHS (LHS_buffer.data (), 2 * N + 2 );
65
70
66
71
/* LAPACK Variables */
67
- int nlhs = (2 * N + 2 );
68
- int nrhs = 1 ;
72
+ int nlhs = (2 * N + 2 );
73
+ int nrhs = 1 ;
69
74
std::vector<int > ipiv (nlhs);
70
75
71
76
/* Stabilization Intensity */
72
77
const double alpha = 0.0 ; // (N * kappa * tau) / (N * tau + 1);
73
- const double L = 10.0 ;
74
- const double dx = L / kappa; // 1.0 / (N * kappa);
78
+ const double L = 10.0 ;
79
+ const double dx = L / kappa; // 1.0 / (N * kappa);
75
80
76
81
// k is the iteration counter
77
- for (int k = 1 ; ; ++k) {
82
+ for (int k = 1 ;; ++k) {
78
83
std::fill (Res.begin (), Res.end (), 0.0 );
79
84
80
85
for (int i = 1 ; i < N; i++) {
81
- /* Momentum */
86
+ /* Momentum */
82
87
Res[i] = (velocity_old[i] * crossSectionLength_old[i] - velocity[i] * crossSectionLength[i]) * dx / tau;
83
88
84
- Res[i] += 0.25 * (- crossSectionLength[i + 1 ] * velocity[i] * velocity[i + 1 ]
85
- - crossSectionLength[i] * velocity[i] * velocity[i + 1 ]);
89
+ Res[i] += 0.25 * (-crossSectionLength[i + 1 ] * velocity[i] * velocity[i + 1 ] - crossSectionLength[i] * velocity[i] * velocity[i + 1 ]);
86
90
87
- Res[i] += 0.25 * (- crossSectionLength[i + 1 ] * velocity[i] * velocity[i]
88
- - crossSectionLength[i] * velocity[i] * velocity[i]
89
- + crossSectionLength[i] * velocity[i - 1 ] * velocity[i]
90
- + crossSectionLength[i - 1 ] * velocity[i - 1 ] * velocity[i]);
91
+ Res[i] += 0.25 * (-crossSectionLength[i + 1 ] * velocity[i] * velocity[i] - crossSectionLength[i] * velocity[i] * velocity[i] + crossSectionLength[i] * velocity[i - 1 ] * velocity[i] + crossSectionLength[i - 1 ] * velocity[i - 1 ] * velocity[i]);
91
92
92
- Res[i] += 0.25 * (+ crossSectionLength[i - 1 ] * velocity[i - 1 ] * velocity[i - 1 ]
93
- + crossSectionLength[i] * velocity[i - 1 ] * velocity[i - 1 ]);
93
+ Res[i] += 0.25 * (+crossSectionLength[i - 1 ] * velocity[i - 1 ] * velocity[i - 1 ] + crossSectionLength[i] * velocity[i - 1 ] * velocity[i - 1 ]);
94
94
95
- Res[i] += 0.25 *(+ crossSectionLength[i - 1 ] * pressure[i - 1 ]
96
- + crossSectionLength[i] * pressure[i - 1 ]
97
- - crossSectionLength[i - 1 ] * pressure[i]
98
- + crossSectionLength[i + 1 ] * pressure[i]
99
- - crossSectionLength[i] * pressure[i + 1 ]
100
- - crossSectionLength[i + 1 ] * pressure[i + 1 ]);
95
+ Res[i] += 0.25 * (+crossSectionLength[i - 1 ] * pressure[i - 1 ] + crossSectionLength[i] * pressure[i - 1 ] - crossSectionLength[i - 1 ] * pressure[i] + crossSectionLength[i + 1 ] * pressure[i] - crossSectionLength[i] * pressure[i + 1 ] - crossSectionLength[i + 1 ] * pressure[i + 1 ]);
101
96
102
97
/* Continuity */
103
98
Res[i + N + 1 ] = (crossSectionLength_old[i] - crossSectionLength[i]) * dx / tau;
104
- Res[i + N + 1 ] += 0.25 *(+ crossSectionLength[i - 1 ] * velocity[i - 1 ]
105
- + crossSectionLength[i] * velocity[i - 1 ]
106
- + crossSectionLength[i - 1 ] * velocity[i]
107
- - crossSectionLength[i + 1 ] * velocity[i]
108
- - crossSectionLength[i] * velocity[i + 1 ]
109
- - crossSectionLength[i + 1 ] * velocity[i + 1 ]);
99
+ Res[i + N + 1 ] += 0.25 * (+crossSectionLength[i - 1 ] * velocity[i - 1 ] + crossSectionLength[i] * velocity[i - 1 ] + crossSectionLength[i - 1 ] * velocity[i] - crossSectionLength[i + 1 ] * velocity[i] - crossSectionLength[i] * velocity[i + 1 ] - crossSectionLength[i + 1 ] * velocity[i + 1 ]);
110
100
111
101
Res[i + N + 1 ] += alpha * (pressure[i - 1 ] - 2 * pressure[i] + pressure[i + 1 ]);
112
102
}
113
103
114
104
/* Boundary */
115
105
116
106
/* Velocity Inlet is prescribed */
117
- const double u0 = 10.0 ;
118
- const double ampl = 3.0 ;
119
- const double frequency = 10.0 ;
120
- const double t_shift = 0.0 ;
107
+ const double u0 = 10.0 ;
108
+ const double ampl = 3.0 ;
109
+ const double frequency = 10.0 ;
110
+ const double t_shift = 0.0 ;
121
111
const double velocity_in = u0 + ampl * sin (frequency * (t + t_shift) * PI);
122
- Res[0 ] = velocity_in - velocity[0 ];
112
+ Res[0 ] = velocity_in - velocity[0 ];
123
113
124
114
/* Pressure Inlet is linearly interpolated */
125
115
Res[N + 1 ] = -pressure[0 ] + 2 * pressure[1 ] - pressure[2 ];
@@ -129,17 +119,15 @@ int fluidComputeSolutionSerial(
129
119
130
120
/* Pressure Outlet is "non-reflecting" */
131
121
const double tmp2 = sqrt (c_mk2 - pressure_old[N] / 2 ) - (velocity[N] - velocity_old[N]) / 4 ;
132
- Res[2 * N + 1 ] = -pressure[N] + 2 * (c_mk2 - std::pow (tmp2, 2 ));
122
+ Res[2 * N + 1 ] = -pressure[N] + 2 * (c_mk2 - std::pow (tmp2, 2 ));
133
123
134
124
// compute norm of residual
135
125
const double norm_1 = std::sqrt (
136
- std::inner_product (Res.begin (), Res.end (), Res.begin (), 0.0 )
137
- );
126
+ std::inner_product (Res.begin (), Res.end (), Res.begin (), 0.0 ));
138
127
139
128
const double norm_2 = std::sqrt (
140
129
std::inner_product (pressure, pressure + chunkLength, pressure, 0.0 ) +
141
- std::inner_product (velocity, velocity + chunkLength, velocity, 0.0 )
142
- );
130
+ std::inner_product (velocity, velocity + chunkLength, velocity, 0.0 ));
143
131
const double norm = norm_1 / norm_2;
144
132
145
133
// NOTE tolerance is 1e-10 and max iterations is 1000 in python
@@ -153,34 +141,25 @@ int fluidComputeSolutionSerial(
153
141
154
142
for (int i = 1 ; i < N; i++) {
155
143
// Momentum, Velocity
156
- LHS (i, i - 1 ) +=0.25 *(-2 * crossSectionLength[i - 1 ] * velocity[i - 1 ]
157
- -2 * crossSectionLength[i] * velocity[i - 1 ]
158
- - crossSectionLength[i] * velocity[i]
159
- - crossSectionLength[i - 1 ] * velocity[i]);
160
-
161
- LHS (i, i) += crossSectionLength[i] * dx / tau;
162
- LHS (i, i) += 0.25 * (+ crossSectionLength[i + 1 ] * velocity[i + 1 ]
163
- + crossSectionLength[i] * velocity[i + 1 ]
164
- +2 * crossSectionLength[i + 1 ] * velocity[i]
165
- +2 * crossSectionLength[i] * velocity[i]
166
- - crossSectionLength[i] * velocity[i - 1 ]
167
- - crossSectionLength[i - 1 ] * velocity[i - 1 ]);
168
- LHS (i, i + 1 ) += 0.25 * (+ crossSectionLength[i + 1 ] * velocity[i]
169
- + crossSectionLength[i] * velocity[i]);
144
+ LHS (i, i - 1 ) += 0.25 * (-2 * crossSectionLength[i - 1 ] * velocity[i - 1 ] - 2 * crossSectionLength[i] * velocity[i - 1 ] - crossSectionLength[i] * velocity[i] - crossSectionLength[i - 1 ] * velocity[i]);
145
+
146
+ LHS (i, i) += crossSectionLength[i] * dx / tau;
147
+ LHS (i, i) += 0.25 * (+crossSectionLength[i + 1 ] * velocity[i + 1 ] + crossSectionLength[i] * velocity[i + 1 ] + 2 * crossSectionLength[i + 1 ] * velocity[i] + 2 * crossSectionLength[i] * velocity[i] - crossSectionLength[i] * velocity[i - 1 ] - crossSectionLength[i - 1 ] * velocity[i - 1 ]);
148
+ LHS (i, i + 1 ) += 0.25 * (+crossSectionLength[i + 1 ] * velocity[i] + crossSectionLength[i] * velocity[i]);
170
149
171
150
// Momentum, Pressure
172
151
LHS (i, N + 1 + i - 1 ) += -0.25 * crossSectionLength[i - 1 ] - 0.25 * crossSectionLength[i];
173
- LHS (i, N + 1 + i) += 0.25 * crossSectionLength[i - 1 ] - 0.25 * crossSectionLength[i + 1 ];
174
- LHS (i, N + 1 + i + 1 ) += 0.25 * crossSectionLength[i] + 0.25 * crossSectionLength[i + 1 ];
152
+ LHS (i, N + 1 + i) += 0.25 * crossSectionLength[i - 1 ] - 0.25 * crossSectionLength[i + 1 ];
153
+ LHS (i, N + 1 + i + 1 ) += 0.25 * crossSectionLength[i] + 0.25 * crossSectionLength[i + 1 ];
175
154
176
155
// Continuity, Velocity
177
156
LHS (i + N + 1 , i - 1 ) += -0.25 * crossSectionLength[i - 1 ] - 0.25 * crossSectionLength[i];
178
- LHS (i + N + 1 , i) += -0.25 * crossSectionLength[i - 1 ] + 0.25 * crossSectionLength[i + 1 ];
179
- LHS (i + N + 1 , i + 1 ) += 0.25 * crossSectionLength[i] + 0.25 * crossSectionLength[i + 1 ];
157
+ LHS (i + N + 1 , i) += -0.25 * crossSectionLength[i - 1 ] + 0.25 * crossSectionLength[i + 1 ];
158
+ LHS (i + N + 1 , i + 1 ) += 0.25 * crossSectionLength[i] + 0.25 * crossSectionLength[i + 1 ];
180
159
181
160
// Continuity, Pressure
182
161
LHS (i + N + 1 , N + 1 + i - 1 ) -= alpha;
183
- LHS (i + N + 1 , N + 1 + i) += 2 * alpha;
162
+ LHS (i + N + 1 , N + 1 + i) += 2 * alpha;
184
163
LHS (i + N + 1 , N + 1 + i + 1 ) -= alpha;
185
164
}
186
165
0 commit comments