HydroCODE_2D 0.1
This is a implementation of fully explict forward Euler scheme for 2-D Euler equations of motion on Eulerian coordinate
flux_solver.c
浏览该文件的文档.
1
6#include <stdio.h>
7#include <math.h>
8
9#include "../include/var_struc.h"
10#include "../include/inter_process.h"
11#include "../include/riemann_solver.h"
12
13
19void Roe_flux(struct i_f_var * ifv, struct i_f_var * ifv_R)
20{
21 const int dim = (int)config[0];
22 const double delta = 0.2;
23
24 double F[4];
25 double lambda_max;
26 if (dim == 1)
27 {
28 Roe_solver(F, &lambda_max, ifv, ifv_R, delta);
29 ifv->F_rho = F[0];
30 ifv->F_u = F[1];
31 ifv->F_e = F[2];
32 }
33 else if (dim == 2)
34 {
35 Roe_2D_solver(F, &lambda_max, ifv, ifv_R, delta);
36 ifv->F_rho = F[0];
37 ifv->F_u = F[1];
38 ifv->F_v = F[2];
39 ifv->F_e = F[3];
40 }
41}
42
43
49void HLL_flux(struct i_f_var * ifv, struct i_f_var * ifv_R)
50{
51 const int dim = (int)config[0];
52
53 double F[4];
54 double lambda_max;
55 if (dim == 2)
56 {
57 HLL_2D_solver(F, &lambda_max, ifv, ifv_R);
58 ifv->F_rho = F[0];
59 ifv->F_u = F[1];
60 ifv->F_v = F[2];
61 ifv->F_e = F[3];
62 }
63}
64
65
75int Riemann_exact_flux(struct i_f_var * ifv, struct i_f_var * ifv_R)
76{
77 const int dim = (int)config[0];
78 const double eps = config[4];
79 const double n_x = ifv->n_x, n_y = ifv->n_y;
80 double gamma_mid = ifv->gamma;
81 ifv->lambda_u = 0.0; ifv->lambda_v = 0.0;
82
83 int retval;
84
85 if (dim == 2)
86 {
87 double u, u_R;
88 u = ifv->U *n_x + ifv->V *n_y;
89 u_R = ifv_R->U*n_x + ifv_R->V*n_y;
90 ifv->V = -ifv->U *n_y + ifv->V *n_x;
91 ifv_R->V = -ifv_R->U*n_y + ifv_R->V*n_x;
92 ifv->U = u;
93 ifv_R->U = u_R;
94 }
95
96 double wave_speed[2], dire[6], mid[6], star[6];
97
98 linear_GRP_solver_Edir_Q1D(wave_speed, dire, mid, star, ifv, ifv_R, eps, INFINITY);
99
100 if((retval = star_dire_check(mid, dire, 2)))
101 return retval;
102
103 double rho_mid = mid[0], p_mid = mid[3], u_mid = mid[1], v_mid = mid[2];
104#ifdef MULTIFLUID_BASICS
105 double phi_mid = mid[4], z_a_mid = mid[5];
106 gamma_mid = mid[1] > 0.0 ? ifv->gamma : ifv_R->gamma;
107#endif
108 if (dim == 1)
109 {
110 u_mid = mid[1];
111 ifv->F_rho = rho_mid*u_mid;
112 ifv->F_u = ifv->F_rho*u_mid + p_mid;
113 }
114 if (dim == 2)
115 {
116 u_mid = mid[1]*n_x - mid[2]*n_y;
117 v_mid = mid[1]*n_y + mid[2]*n_x;
118 ifv->F_rho = rho_mid*(u_mid*n_x + v_mid*n_y);
119 ifv->F_u = ifv->F_rho*u_mid + p_mid*n_x;
120 ifv->F_v = ifv->F_rho*v_mid + p_mid*n_y;
121 }
122 ifv->F_e = (gamma_mid/(gamma_mid-1.0))*p_mid/rho_mid + 0.5*u_mid*u_mid;
123 if (dim >= 2)
124 ifv->F_e += 0.5*v_mid*v_mid;
125 ifv->F_e = ifv->F_rho*ifv->F_e;
126
127#ifdef MULTIFLUID_BASICS
128 ifv->F_phi = ifv->F_rho * phi_mid;
129 if ((_Bool)config[60])
130 ifv->F_gamma = ifv->F_rho*gamma_mid;
131 ifv->F_e_a = z_a_mid/(config[6]-1.0)*p_mid/rho_mid + 0.5*phi_mid*u_mid*u_mid;
132 if (dim >= 2)
133 ifv->F_e_a += 0.5*phi_mid*v_mid*v_mid;
134 ifv->F_e_a = ifv->F_rho*ifv->F_e_a;
135
136 ifv->U_qt_add_c = ifv->F_rho*u_mid*phi_mid;
137 if (dim >= 2)
138 ifv->V_qt_add_c = ifv->F_rho*v_mid*phi_mid;
139 ifv->U_qt_star = p_mid*n_x;
140 ifv->V_qt_star = p_mid*n_y;
141 ifv->P_star = p_mid/rho_mid*ifv->F_rho;
142#endif
143 return retval;
144}
145
146
158int GRP_2D_flux(struct i_f_var * ifv, struct i_f_var * ifv_R, const double tau)
159{
160 const double eps = config[4];
161 const double n_x = ifv->n_x, n_y = ifv->n_y;
162 double gamma_mid = ifv->gamma;
163 ifv->lambda_u = 0.0; ifv->lambda_v = 0.0;
164
165 int retval;
166
167 double u, u_R, d_u, d_u_R, t_u, t_u_R;
168 u = ifv->U *n_x + ifv->V *n_y;
169 u_R = ifv_R->U *n_x + ifv_R->V *n_y;
170 d_u = ifv->d_u *n_x + ifv->d_v *n_y;
171 d_u_R = ifv_R->d_u*n_x + ifv_R->d_v*n_y;
172 t_u = ifv->t_u *n_x + ifv->t_v *n_y;
173 t_u_R = ifv_R->t_u*n_x + ifv_R->t_v*n_y;
174 ifv->V = -ifv->U *n_y + ifv->V *n_x;
175 ifv_R->V = -ifv_R->U *n_y + ifv_R->V *n_x;
176 ifv->d_v = -ifv->d_u *n_y + ifv->d_v *n_x;
177 ifv_R->d_v = -ifv_R->d_u*n_y + ifv_R->d_v*n_x;
178 ifv->t_v = -ifv->t_u *n_y + ifv->t_v *n_x;
179 ifv_R->t_v = -ifv_R->t_u*n_y + ifv_R->t_v*n_x;
180 ifv->U = u;
181 ifv_R->U = u_R;
182 ifv->d_u = d_u;
183 ifv_R->d_u = d_u_R;
184 ifv->t_u = t_u;
185 ifv_R->t_u = t_u_R;
186
187 double wave_speed[2], dire[6], mid[6], star[6];
188
189 // linear_GRP_solver_Edir_G2D(wave_speed, dire, mid, star, ifv, ifv_R, eps, eps);
190 // linear_GRP_solver_Edir_G2D(wave_speed, dire, mid, star, ifv, ifv_R, eps, INFINITY);
191 linear_GRP_solver_Edir_Q1D(wave_speed, dire, mid, star, ifv, ifv_R, eps, eps);
192 // linear_GRP_solver_Edir_Q1D(wave_speed, dire, mid, star, ifv, ifv_R, eps, -0.0);
193 // linear_GRP_solver_Edir_Q1D(wave_speed, dire, mid, star, ifv, ifv_R, eps, INFINITY);
194
195 if((retval = star_dire_check(mid, dire, 2)))
196 return retval;
197
198 double rho_mid, p_mid, u_mid, v_mid;
199 rho_mid = mid[0] + 0.5*tau*dire[0];
200 u_mid = (mid[1] + 0.5*tau*dire[1])*n_x - (mid[2] + 0.5*tau*dire[2])*n_y;
201 v_mid = (mid[1] + 0.5*tau*dire[1])*n_y + (mid[2] + 0.5*tau*dire[2])*n_x;
202 p_mid = mid[3] + 0.5*tau*dire[3];
203#ifdef MULTIFLUID_BASICS
204 double phi_mid, z_a_mid;
205 phi_mid = mid[5] + 0.5*tau*dire[5];
206 z_a_mid = mid[4] + 0.5*tau*dire[4];
207 gamma_mid = 1.0/(z_a_mid/(config[6]-1.0)+(1.0-z_a_mid)/(config[106]-1.0))+1.0;
208#endif
209
210 ifv->F_rho = rho_mid*(u_mid*n_x + v_mid*n_y);
211 ifv->F_u = ifv->F_rho*u_mid + p_mid*n_x;
212 ifv->F_v = ifv->F_rho*v_mid + p_mid*n_y;
213 ifv->F_e = (gamma_mid/(gamma_mid-1.0))*p_mid/rho_mid + 0.5*(u_mid*u_mid + v_mid*v_mid);
214 ifv->F_e = ifv->F_rho*ifv->F_e;
215
216 ifv->U_int = (mid[1] + tau*dire[1])*n_x - (mid[2] + tau*dire[2])*n_y;
217 ifv->V_int = (mid[1] + tau*dire[1])*n_y + (mid[2] + tau*dire[2])*n_x;
218 ifv->RHO_int = mid[0] + tau*dire[0];
219 ifv->P_int = mid[3] + tau*dire[3];
220
221#ifdef MULTIFLUID_BASICS
222 ifv->F_phi = ifv->F_rho*phi_mid;
223 if ((_Bool)config[60])
224 ifv->F_gamma = ifv->F_rho*gamma_mid;
225 ifv->F_e_a = z_a_mid/(config[6]-1.0)*p_mid/rho_mid + 0.5*phi_mid*(u_mid*u_mid + v_mid*v_mid);
226 ifv->F_e_a = ifv->F_rho*ifv->F_e_a;
227 ifv->PHI = mid[5] + tau*dire[5];
228 ifv->Z_a = mid[4] + tau*dire[4];
229
230 ifv->U_qt_add_c = ifv->F_rho*u_mid*phi_mid;
231 ifv->V_qt_add_c = ifv->F_rho*v_mid*phi_mid;
232 ifv->U_qt_star = p_mid*n_x;
233 ifv->V_qt_star = p_mid*n_y;
234 ifv->P_star = p_mid/rho_mid*ifv->F_rho;
235#endif
236 return retval;
237}
GRP_2D_flux
int GRP_2D_flux(struct i_f_var *ifv, struct i_f_var *ifv_R, const double tau)
This function calculate Eulerian fluxes of 2-D Euler equations by 2-D GRP solver.
Definition: flux_solver.c:158
HLL_flux
void HLL_flux(struct i_f_var *ifv, struct i_f_var *ifv_R)
This function calculate Eulerian fluxes of 2-D Euler equations by HLL solver.
Definition: flux_solver.c:49
Riemann_exact_flux
int Riemann_exact_flux(struct i_f_var *ifv, struct i_f_var *ifv_R)
This function calculate Eulerian fluxes of 2-D Euler equations by Riemann solver.
Definition: flux_solver.c:75
Roe_flux
void Roe_flux(struct i_f_var *ifv, struct i_f_var *ifv_R)
This function calculate Eulerian fluxes of Euler equations by Roe solver.
Definition: flux_solver.c:19
config
double config[N_CONF]
Initial configuration data array.
Definition: hydrocode.c:123
star_dire_check
int star_dire_check(double *mid, double *dire, const int dim)
This function checks whether fluid variables of mid[] and dire[] are within the value range.
Definition: fluid_var_check.c:65
HLL_2D_solver
void HLL_2D_solver(double *F, double *lambda_max, const struct i_f_var *ifv_L, const struct i_f_var *ifv_R)
A HLL approxiamate Riemann solver for unsteady compressible inviscid single-component flow in two spa...
Definition: hll_2D_solver.c:25
Roe_2D_solver
void Roe_2D_solver(double *F, double *lambda_max, const struct i_f_var *ifv_L, const struct i_f_var *ifv_R, const double delta)
An approxiamate Riemann solver of Roe for unsteady compressible inviscid single-component flow in two...
Definition: roe_2D_solver.c:26
Roe_solver
void Roe_solver(double *F, double *lambda_max, const struct i_f_var *ifv_L, const struct i_f_var *ifv_R, const double delta)
An approxiamate Riemann solver of Roe for unsteady compressible inviscid single-component flow in one...
Definition: roe_solver.c:25
linear_GRP_solver_Edir_Q1D
void linear_GRP_solver_Edir_Q1D(double *wave_speed, double *D, double *U, double *U_star, const struct i_f_var *ifv_L, const struct i_f_var *ifv_R, const double eps, const double atc)
A Quasi-1D direct Eulerian GRP solver for unsteady compressible inviscid two-component flow in two sp...
Definition: linear_grp_solver_Edir_Q1D.c:39
i_f_var
Interfacial Fluid VARiables.
Definition: var_struc.h:116
i_f_var::V
double V
primitive variable values at t_{n}.
Definition: var_struc.h:122
i_f_var::V_int
double V_int
interfacial primitive variables at t_{n+1}.
Definition: var_struc.h:123
i_f_var::lambda_v
double lambda_v
grid moving velocity components in direction x and y.
Definition: var_struc.h:120
i_f_var::t_v
double t_v
Definition: var_struc.h:127
i_f_var::U_int
double U_int
Definition: var_struc.h:123
i_f_var::U
double U
Definition: var_struc.h:122
i_f_var::n_y
double n_y
x- and y-coordinates of the interfacial unit normal vector.
Definition: var_struc.h:117
i_f_var::gamma
double gamma
specific heat ratio.
Definition: var_struc.h:121
i_f_var::F_rho
double F_rho
Definition: var_struc.h:124
i_f_var::RHO_int
double RHO_int
Definition: var_struc.h:123
i_f_var::n_x
double n_x
Definition: var_struc.h:117
i_f_var::F_e
double F_e
Definition: var_struc.h:124
i_f_var::d_u
double d_u
Definition: var_struc.h:126
i_f_var::t_u
double t_u
Definition: var_struc.h:127
i_f_var::F_u
double F_u
Definition: var_struc.h:124
i_f_var::F_v
double F_v
interfacial fluxes at t_{n+1/2}.
Definition: var_struc.h:124
i_f_var::lambda_u
double lambda_u
Definition: var_struc.h:120
i_f_var::d_v
double d_v
Definition: var_struc.h:126
i_f_var::P_int
double P_int
Definition: var_struc.h:123