Skip to content

Commit c82d61c

Browse files
author
nianhuawong
committed
二阶steger-warming还是算不了
1 parent 847921b commit c82d61c

File tree

4 files changed

+35
-56
lines changed

4 files changed

+35
-56
lines changed

‎src/Flux_Solver.cpp‎

Lines changed: 19 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -707,8 +707,8 @@ void Flux_Solver::Steger_Warming_Scheme_Interp_X()
707707
double eps = 1e-4;
708708
VInt2D& marker = mesh->Get_Marker();
709709

710-
VDouble fluxVector1(num_of_prim_vars);
711-
VDouble fluxVector2(num_of_prim_vars);
710+
VDouble fluxVector11(num_of_prim_vars);
711+
VDouble fluxVector22(num_of_prim_vars);
712712
for (int j = jst; j <= jed; j++)
713713
{
714714
for (int i = 0; i <= ied + 1; i++)//每个通量点的值都有了
@@ -722,7 +722,7 @@ void Flux_Solver::Steger_Warming_Scheme_Interp_X()
722722
double u = qField1[i][j][IU];
723723
double v = qField1[i][j][IV];
724724
double p = qField1[i][j][IP];
725-
double a = sqrt(gama * p / rho);//声速取绝对值
725+
double a = sqrt(fabs(gama * p / rho));//声速取绝对值
726726

727727
double lmd1 = u;
728728
double lmd3 = u - a;
@@ -738,13 +738,13 @@ void Flux_Solver::Steger_Warming_Scheme_Interp_X()
738738
int kkk = 1;
739739
}
740740

741-
Steger_Flux_F(fluxVector1, rho, u, v, p, lmd_m);
741+
Steger_Flux_F(fluxVector11, rho, u, v, p, lmd_m);
742742

743743
rho = qField2[i][j][IR];
744744
u = qField2[i][j][IU];
745745
v = qField2[i][j][IV];
746746
p = qField2[i][j][IP];
747-
a = sqrt(gama * p / rho);//声速取绝对值
747+
a = sqrt(fabs(gama * p / rho));//声速取绝对值
748748

749749
lmd1 = u;
750750
lmd3 = u - a;
@@ -760,12 +760,12 @@ void Flux_Solver::Steger_Warming_Scheme_Interp_X()
760760
int kkk = 1;
761761
}
762762

763-
Steger_Flux_F(fluxVector2, rho, u, v, p, lmd_p);
763+
Steger_Flux_F(fluxVector22, rho, u, v, p, lmd_p);
764764

765-
fluxVector[i][j][IR] = fluxVector1[IR] + fluxVector2[IR];
766-
fluxVector[i][j][IU] = fluxVector1[IU] + fluxVector2[IU];
767-
fluxVector[i][j][IV] = fluxVector1[IV] + fluxVector2[IV];
768-
fluxVector[i][j][IP] = fluxVector1[IP] + fluxVector2[IP];
765+
fluxVector[i][j][IR] = fluxVector11[IR] + fluxVector22[IR];
766+
fluxVector[i][j][IU] = fluxVector11[IU] + fluxVector22[IU];
767+
fluxVector[i][j][IV] = fluxVector11[IV] + fluxVector22[IV];
768+
fluxVector[i][j][IP] = fluxVector11[IP] + fluxVector22[IP];
769769
}
770770
}
771771
}
@@ -775,8 +775,8 @@ void Flux_Solver::Steger_Warming_Scheme_Interp_Y()
775775
double eps = 1e-4;
776776
VInt2D& marker = mesh->Get_Marker();
777777

778-
VDouble fluxVector1(num_of_prim_vars);
779-
VDouble fluxVector2(num_of_prim_vars);
778+
VDouble fluxVector11(num_of_prim_vars);
779+
VDouble fluxVector22(num_of_prim_vars);
780780

781781
for (int j = 0; j <= jed + 1; j++) //每个通量点的值都有了
782782
{
@@ -807,7 +807,7 @@ void Flux_Solver::Steger_Warming_Scheme_Interp_Y()
807807
int kkk = 1;
808808
}
809809

810-
Steger_Flux_G(fluxVector1, rho, u, v, p, mu_m);
810+
Steger_Flux_G(fluxVector11, rho, u, v, p, mu_m);
811811

812812
rho = qField2[i][j][IR];
813813
u = qField2[i][j][IU];
@@ -829,12 +829,12 @@ void Flux_Solver::Steger_Warming_Scheme_Interp_Y()
829829
int kkk = 1;
830830
}
831831

832-
Steger_Flux_G(fluxVector2, rho, u, v, p, mu_p);
832+
Steger_Flux_G(fluxVector22, rho, u, v, p, mu_p);
833833

834-
fluxVector[i][j][IR] = fluxVector1[IR] + fluxVector2[IR];
835-
fluxVector[i][j][IU] = fluxVector1[IU] + fluxVector2[IU];
836-
fluxVector[i][j][IV] = fluxVector1[IV] + fluxVector2[IV];
837-
fluxVector[i][j][IP] = fluxVector1[IP] + fluxVector2[IP];
834+
fluxVector[i][j][IR] = fluxVector11[IR] + fluxVector22[IR];
835+
fluxVector[i][j][IU] = fluxVector11[IU] + fluxVector22[IU];
836+
fluxVector[i][j][IV] = fluxVector11[IV] + fluxVector22[IV];
837+
fluxVector[i][j][IP] = fluxVector11[IP] + fluxVector22[IP];
838838
}
839839
}
840840
}
@@ -852,7 +852,7 @@ void Flux_Solver::Steger_Flux_F(VDouble& fluxVector, double rho, double u, doubl
852852

853853
void Flux_Solver::Steger_Flux_G(VDouble& fluxVector, double rho, double u, double v, double p, VDouble mu)
854854
{
855-
double a = sqrt(fabs(gama * p / rho)); //声速取绝对值
855+
double a = sqrt(fabs(gama * p / rho));//声速取绝对值
856856
double h = Enthalpy(rho, u, v, p, gama);
857857

858858
fluxVector[IR] = rho / (2 * gama) * (2 * (gama - 1) * mu[0] + mu[1] + mu[2]);

‎src/QlQr_Solver.cpp‎

Lines changed: 8 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -74,11 +74,11 @@ void QlQr_Solver::Boundary_QlQr_MUSCL_X()
7474
{
7575
qField1[0][j][iVar] = qField[0][j][iVar]; //虚拟点i=0, ied + 1的值还没有
7676
qField2[0][j][iVar] = qField[1][j][iVar];
77-
qField1[1][j][iVar] = qField[1][j][iVar];
78-
qField2[1][j][iVar] = qField[2][j][iVar];
77+
//qField1[1][j][iVar] = qField[1][j][iVar];
78+
//qField2[1][j][iVar] = qField[2][j][iVar];
7979

80-
qField1[ied ][j][iVar] = qField[ied ][j][iVar];
81-
qField2[ied ][j][iVar] = qField[ied + 1][j][iVar];
80+
//qField1[ied ][j][iVar] = qField[ied ][j][iVar];
81+
//qField2[ied ][j][iVar] = qField[ied + 1][j][iVar];
8282
qField1[ied + 1][j][iVar] = qField[ied + 1][j][iVar];
8383
qField2[ied + 1][j][iVar] = qField[ied + 2][j][iVar];
8484
}
@@ -88,30 +88,18 @@ void QlQr_Solver::Boundary_QlQr_MUSCL_X()
8888
void QlQr_Solver::Boundary_QlQr_MUSCL_Y()
8989
{
9090
//VInt2D& marker = mesh->Get_Marker();
91-
//for (int i = ist; i < ied; i++)
92-
//{
93-
// for (int j = 0; j < num_half_point_y; j++)
94-
// {
95-
// if (marker[i][j] == 0) continue;
96-
// for (int iVar = 0; iVar < num_of_prim_vars; iVar++)
97-
// {
98-
// qField1[i][j][iVar] = qField[i][j ][iVar];
99-
// qField2[i][j][iVar] = qField[i][j + 1][iVar];
100-
// }
101-
// }
102-
//}
10391

10492
for (int i = ist; i <= ied; i++)
10593
{
10694
for (int iVar = 0; iVar < num_of_prim_vars; iVar++)//虚拟点j=0, jed + 1的值还没有
10795
{
10896
qField1[i][0][iVar] = qField[i][0][iVar];
10997
qField2[i][0][iVar] = qField[i][1][iVar];
110-
qField1[i][1][iVar] = qField[i][1][iVar];
111-
qField2[i][1][iVar] = qField[i][2][iVar];
98+
//qField1[i][1][iVar] = qField[i][1][iVar];
99+
//qField2[i][1][iVar] = qField[i][2][iVar];
112100

113-
qField1[i][jed ][iVar] = qField[i][jed ][iVar];
114-
qField2[i][jed ][iVar] = qField[i][jed + 1][iVar];
101+
//qField1[i][jed ][iVar] = qField[i][jed ][iVar];
102+
//qField2[i][jed ][iVar] = qField[i][jed + 1][iVar];
115103
qField1[i][jed + 1][iVar] = qField[i][jed + 1][iVar];
116104
qField2[i][jed + 1][iVar] = qField[i][jed + 2][iVar];
117105
}

‎src/Spatial_Derivative.cpp‎

Lines changed: 4 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -56,16 +56,12 @@ void Spatial_Derivative::Spatial_Derivative_X()
5656
{
5757
VInt2D& marker = mesh->Get_Marker();
5858

59-
//for (int j = jst; j < jed - 1; j++)
60-
//{
61-
// for (int i = ist; i < ied - 1; i++)
62-
// {
6359
#ifdef _OPENMP
6460
#pragma omp parallel for
6561
#endif
66-
for (int i = 1; i < num_half_point_x; i++)
62+
for (int i = 1; i <= ied + 1; i++)//i=0和i=ied+2的空间导数项没有值
6763
{
68-
for (int j = jst; j < jed; j++)
64+
for (int j = jst; j <= jed; j++)
6965
{
7066
if (marker[i][j] == 0) continue;
7167
if (i == 61 && j == 62)
@@ -84,23 +80,18 @@ void Spatial_Derivative::Spatial_Derivative_X()
8480
}
8581
}
8682
}
87-
//i=0和i=total_num_pointx的空间导数项没有值
8883
}
8984

9085
void Spatial_Derivative::Spatial_Derivative_Y()
9186
{
9287
VInt2D& marker = mesh->Get_Marker();
9388

94-
//for (int i = ist; i < ied - 1; i++)
95-
//{
96-
// for (int j = jst; j < jed - 1; j++)
97-
// {
9889
#ifdef _OPENMP
9990
#pragma omp parallel for
10091
#endif
101-
for (int i = ist; i < ied; i++)
92+
for (int i = ist; i <= ied; i++)
10293
{
103-
for (int j = 1; j < num_half_point_y; j++)
94+
for (int j = 1; j <= jed + 1; j++)//j=0和j=jed+2的空间导数项没有值
10495
{
10596
if (marker[i][j] == 0) continue;
10697
if (i == 16 && j == 2)//((i==60||i == 61) && j == 62)

‎src/Time_Integral.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -80,7 +80,7 @@ void Update_Flowfield_X(int iStage)
8080
VDouble qConservative1(num_of_prim_vars);
8181
#pragma omp for
8282
#endif
83-
for (int i = 0; i <= ied; i++)
83+
for (int i = 1; i <= ied + 1; i++)
8484
{
8585
for (int j = jst; j <= jed; j++)
8686
{
@@ -108,7 +108,7 @@ void Update_Flowfield_X(int iStage)
108108
if (qPrimitive1[IR] < minimumDensityLimit || qPrimitive1[IP] < minimumPressureLimit)
109109
{
110110
int kkk = 1;
111-
//SolutionFix(qPrimitive1, i, j);
111+
SolutionFix(qPrimitive1, i, j);
112112
}
113113
//RK公式里左端项，q1、q2、q3，即下一stage的q值，还要继续用该值计算rhs(q1)、rhs(q2)
114114
qField_N1[i][j] = qPrimitive1;
@@ -153,7 +153,7 @@ void Update_Flowfield_Y(int iStage)
153153
#endif
154154
for (int i = ist; i <= ied; i++)
155155
{
156-
for (int j = 0; j <= jed; j++)
156+
for (int j = 1; j <= jed + 1; j++)
157157
{
158158
if (i == 16 && j == 2)
159159
{
@@ -179,7 +179,7 @@ void Update_Flowfield_Y(int iStage)
179179
if (qPrimitive1[IR] < minimumDensityLimit || qPrimitive1[IP] < minimumPressureLimit)
180180
{
181181
int kkk = 1;
182-
//SolutionFix(qPrimitive1, i, j);
182+
SolutionFix(qPrimitive1, i, j);
183183
}
184184
//RK公式里左端项，q1、q2、q3，即下一stage的q值，还要继续用该值计算rhs(q1)、rhs(q2)
185185
qField_N1[i][j] = qPrimitive1;

0 commit comments

Comments
 (0)