@@ -32,20 +32,22 @@ parse_impl(int argc, char *argv[])
32
32
return argv [1 ];
33
33
}
34
34
else {
35
- fprintf (stderr , "USAGE: %s [scipy_ode_dopri5 | sundials_cvode | jl_diffeq]\n" , argv [0 ]);
35
+ fprintf (stderr , "USAGE: %s [scipy_ode_dopri5 | sundials_cvode | jl_diffeq]\n" ,
36
+ argv [0 ]);
36
37
exit (EXIT_FAILURE );
37
38
}
38
39
}
39
40
}
40
41
41
42
static int
42
- compute_initial_condition_ (size_t N , OIFArrayF64 * u0 , OIFArrayF64 * grid , double * dx , double * dt_max )
43
+ compute_initial_condition_ (size_t N , OIFArrayF64 * u0 , OIFArrayF64 * grid , double * dx ,
44
+ double * dt_max )
43
45
{
44
46
double a = 0.0 ;
45
47
double b = 2.0 ;
46
48
double * x = grid -> data ;
47
49
* dx = (b - a ) / N ;
48
-
50
+
49
51
for (int i = 0 ; i < N ; ++ i ) {
50
52
x [i ] = a + i * (* dx );
51
53
}
@@ -72,7 +74,7 @@ rhs(double t, OIFArrayF64 *y, OIFArrayF64 *rhs_out, void *user_data)
72
74
double * u = y -> data ;
73
75
double * udot = rhs_out -> data ;
74
76
75
- double dx = * ((double * ) user_data );
77
+ double dx = * ((double * )user_data );
76
78
77
79
double * flux = malloc (N * sizeof (double ));
78
80
if (flux == NULL ) {
@@ -99,18 +101,17 @@ rhs(double t, OIFArrayF64 *y, OIFArrayF64 *rhs_out, void *user_data)
99
101
}
100
102
101
103
for (int i = 0 ; i < N - 1 ; ++ i ) {
102
- flux_hat [i ] = 0.5 * ( flux [ i ] + flux [ i + 1 ]) -
103
- 0.5 * local_sound_speed * (u [i + 1 ] - u [i ]);
104
+ flux_hat [i ] =
105
+ 0.5 * ( flux [ i ] + flux [ i + 1 ]) - 0.5 * local_sound_speed * (u [i + 1 ] - u [i ]);
104
106
}
105
107
106
108
for (int i = 1 ; i < N - 1 ; ++ i ) {
107
109
udot [i ] = -1.0 / dx * (flux_hat [i ] - flux_hat [i - 1 ]);
108
110
}
109
- double f_rb = 0.5 * (flux [0 ] + flux [N - 1 ]) -
110
- 0.5 * local_sound_speed * (u [0 ] - u [N - 1 ]);
111
+ double f_rb = 0.5 * (flux [0 ] + flux [N - 1 ]) - 0.5 * local_sound_speed * (u [0 ] - u [N - 1 ]);
111
112
double f_lb = f_rb ;
112
113
udot [0 ] = -1.0 / dx * (flux_hat [0 ] - f_lb );
113
- udot [N - 1 ] = -1.0 / dx * (f_rb - flux_hat [N - 2 ]);
114
+ udot [N - 1 ] = -1.0 / dx * (f_rb - flux_hat [N - 2 ]);
114
115
115
116
retval = 0 ;
116
117
@@ -146,7 +147,6 @@ main(int argc, char *argv[])
146
147
double dt_max ;
147
148
int status = 1 ; // Aux variable to check for errors.
148
149
149
-
150
150
status = compute_initial_condition_ (N , y0 , grid , & dx , & dt_max );
151
151
assert (status == 0 );
152
152
@@ -179,7 +179,7 @@ main(int argc, char *argv[])
179
179
180
180
// Time step.
181
181
double dt = dt_max ;
182
- int n_time_steps = (int ) (t_final / dt + 1 );
182
+ int n_time_steps = (int )(t_final / dt + 1 );
183
183
for (int i = 0 ; i < n_time_steps ; ++ i ) {
184
184
double t = t0 + (i + 1 ) * dt ;
185
185
if (t > t_final ) {
@@ -212,4 +212,3 @@ main(int argc, char *argv[])
212
212
213
213
return retval ;
214
214
}
215
-
0 commit comments