Skip to content

Commit 02be3f6

Browse files
pcarruscagpcarruscag
authored
Merge pull request #2486 from su2code/symmetry_jacobian
Improve the Euler wall Jacobian
2 parents c9f098d + 551d7c5 commit 02be3f6

File tree

27 files changed

+597
-571
lines changed

27 files changed

+597
-571
lines changed

‎SU2_CFD/include/solvers/CFVMFlowSolverBase.inl‎

Lines changed: 63 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -1141,41 +1141,54 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
11411141
for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area;
11421142
}
11431143

1144-
su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex);
1144+
/*--- Energy terms due to grid movement (aka work of pressure forces). ---*/
1145+
if (dynamic_grid) {
1146+
su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex);
11451147

1146-
/*--- Grid movement ---*/
1147-
if (dynamic_grid)
1148-
conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
1148+
conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint),
1149+
geometry->nodes->GetGridVel(iPoint));
11491150

1150-
/*--- Normal vector for this vertex (negate for outward convention). ---*/
1151-
for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim];
1152-
conv_numerics->SetNormal(Normal);
1151+
/*--- Normal vector for this vertex (negate for outward convention). ---*/
1152+
for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim];
1153+
conv_numerics->SetNormal(Normal);
11531154

1154-
for (auto iVar = 0u; iVar < nPrimVar; iVar++)
1155-
V_reflected[iVar] = nodes->GetPrimitive(iPoint, iVar);
1155+
for (auto iVar = 0u; iVar < nPrimVar; iVar++)
1156+
V_reflected[iVar] = nodes->GetPrimitive(iPoint, iVar);
11561157

1157-
su2double ProjVelocity_i = nodes->GetProjVel(iPoint, UnitNormal);
1158-
/*--- Adjustment to v.n due to grid movement. ---*/
1159-
if (dynamic_grid)
1158+
su2double ProjVelocity_i = nodes->GetProjVel(iPoint, UnitNormal);
1159+
/*--- Adjustment to v.n due to grid movement. ---*/
11601160
ProjVelocity_i -= GeometryToolbox::DotProduct(nDim, geometry->nodes->GetGridVel(iPoint), UnitNormal);
11611161

1162-
for (auto iDim = 0u; iDim < nDim; iDim++)
1163-
V_reflected[iDim + iVel] = nodes->GetVelocity(iPoint, iDim) - ProjVelocity_i * UnitNormal[iDim];
1162+
for (auto iDim = 0u; iDim < nDim; iDim++)
1163+
V_reflected[iDim + iVel] = nodes->GetVelocity(iPoint, iDim) - ProjVelocity_i * UnitNormal[iDim];
11641164

1165-
/*--- Get current solution at this boundary node ---*/
1166-
const su2double* V_domain = nodes->GetPrimitive(iPoint);
1165+
/*--- Get current solution at this boundary node. ---*/
1166+
const su2double* V_domain = nodes->GetPrimitive(iPoint);
11671167

1168-
/*--- Set Primitive and Secondary for numerics class. ---*/
1169-
conv_numerics->SetPrimitive(V_domain, V_reflected);
1170-
conv_numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint));
1168+
/*--- Set Primitive and Secondary for numerics class. ---*/
1169+
conv_numerics->SetPrimitive(V_domain, V_reflected);
1170+
conv_numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint));
11711171

1172-
/*--- Compute the residual using an upwind scheme. ---*/
1173-
auto residual = conv_numerics->ComputeResidual(config);
1172+
/*--- Compute the residual using an upwind scheme. ---*/
1173+
auto residual = conv_numerics->ComputeResidual(config);
11741174

1175-
/*--- We include an update of the continuity and energy here, this is important for stability since
1176-
* these fluxes include numerical diffusion. ---*/
1177-
for (auto iVar = 0u; iVar < nVar; iVar++) {
1178-
if (iVar < iVel || iVar >= iVel + nDim) LinSysRes(iPoint, iVar) += residual.residual[iVar];
1175+
/*--- Use just the energy fluxes to update the residual, adding the others would
1176+
* increase numerical diffusion which we wish to avoid if possible. ---*/
1177+
for (auto iVar = iVel + nDim; iVar < nVar; iVar++) {
1178+
LinSysRes(iPoint, iVar) += residual.residual[iVar];
1179+
}
1180+
if (implicit) {
1181+
auto* block = Jacobian.GetBlock(iPoint, iPoint);
1182+
/*--- But in the Jacobian we also include the mass flux, this allows some cases with
1183+
* motion to use larger CFL, for example pywrapper_translating_naca0012. ---*/
1184+
for (auto iVar = 0u; iVar < nVar; iVar++) {
1185+
if (iVar < iVel || iVar >= iVel + nDim) {
1186+
for (auto jVar = 0u; jVar < nVar; jVar++) {
1187+
block[iVar * nVar + jVar] += SU2_TYPE::GetValue(residual.jacobian_i[iVar][jVar]);
1188+
}
1189+
}
1190+
}
1191+
}
11791192
}
11801193

11811194
/*--- Explicitly set the velocity components normal to the symmetry plane to zero.
@@ -1184,7 +1197,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
11841197

11851198
su2double* solutionOld = nodes->GetSolution_Old(iPoint);
11861199

1187-
su2double gridVel[MAXNVAR] = {};
1200+
su2double gridVel[MAXNDIM] = {};
11881201
if (dynamic_grid) {
11891202
for (auto iDim = 0u; iDim < nDim; iDim++) {
11901203
gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim];
@@ -1215,7 +1228,30 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
12151228

12161229
/*--- Jacobian contribution for implicit integration. ---*/
12171230
if (implicit) {
1218-
Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
1231+
/*--- Modify the Jacobians according to the modification of the residual
1232+
* J_new = (I - n * n^T) * J where n = {0, nx, ny, nz, 0, ...} ---*/
1233+
su2double mat[MAXNVAR * MAXNVAR] = {};
1234+
1235+
for (auto iVar = 0u; iVar < nVar; iVar++)
1236+
mat[iVar * nVar + iVar] = 1;
1237+
for (auto iDim = 0u; iDim < nDim; iDim++)
1238+
for (auto jDim = 0u; jDim < nDim; jDim++)
1239+
mat[(iDim + iVel) * nVar + jDim + iVel] -= UnitNormal[iDim] * UnitNormal[jDim];
1240+
1241+
auto ModifyJacobian = [&](const unsigned long jPoint) {
1242+
su2double jac[MAXNVAR * MAXNVAR], newJac[MAXNVAR * MAXNVAR];
1243+
auto* block = Jacobian.GetBlock(iPoint, jPoint);
1244+
for (auto iVar = 0u; iVar < nVar * nVar; iVar++) jac[iVar] = block[iVar];
1245+
1246+
CBlasStructure().gemm(nVar, nVar, nVar, mat, jac, newJac, config);
1247+
1248+
for (auto iVar = 0u; iVar < nVar * nVar; iVar++)
1249+
block[iVar] = SU2_TYPE::GetValue(newJac[iVar]);
1250+
};
1251+
ModifyJacobian(iPoint);
1252+
for (size_t iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
1253+
ModifyJacobian(geometry->nodes->GetPoint(iPoint, iNeigh));
1254+
}
12191255
}
12201256

12211257
/*--- Correction for multigrid. ---*/

‎TestCases/Tutorials/design/Inc_Turbulent_Bend_Wallfunctions/optim.csv.ref‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,2 @@
11
ITER, avg_dp
2-
1, 63.82289754521027
2+
1, 64.90301114659223

‎TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref‎

Lines changed: 38 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1,39 +1,39 @@
11
VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP"
2-
0 , -2656.0 , 0.001
3-
1 , -7290.69 , 0.001
4-
2 , -11252.7 , 0.001
5-
3 , -12531.4 , 0.001
6-
4 , -10923.7 , 0.001
7-
5 , -6826.86 , 0.001
8-
6 , -843.71 , 0.001
9-
7 , 6462.74 , 0.001
10-
8 , 14656.3 , 0.001
11-
9 , 23308.0 , 0.001
12-
10 , 31688.2 , 0.001
13-
11 , 38407.2 , 0.001
14-
12 , 41274.9 , 0.001
15-
13 , 37812.2 , 0.001
16-
14 , 27093.1 , 0.001
17-
15 , 14049.9 , 0.001
18-
16 , 17583.8 , 0.001
19-
17 , 78620.8 , 0.001
20-
18 , 140614.0 , 0.001
21-
19 , -18263.2 , 0.001
22-
20 , -20877.9 , 0.001
23-
21 , -22859.6 , 0.001
24-
22 , -24751.5 , 0.001
25-
23 , -28392.5 , 0.001
26-
24 , -35514.1 , 0.001
27-
25 , -46855.4 , 0.001
28-
26 , -61697.7 , 0.001
29-
27 , -77631.1 , 0.001
30-
28 , -90537.8 , 0.001
31-
29 , -95036.8 , 0.001
32-
30 , -85899.0 , 0.001
33-
31 , -60885.3 , 0.001
34-
32 , -24333.4 , 0.001
35-
33 , 12334.3 , 0.001
36-
34 , 39948.1 , 0.001
37-
35 , 65452.2 , 0.001
38-
36 , 82215.8 , 0.001
39-
37 , -44066.2 , 0.001
2+
0 , 292.459 , 0.001
3+
1 , -8318.5 , 0.001
4+
2 , -16158.4 , 0.001
5+
3 , -21277.2 , 0.001
6+
4 , -23366.9 , 0.001
7+
5 , -22727.6 , 0.001
8+
6 , -19872.3 , 0.001
9+
7 , -15278.1 , 0.001
10+
8 , -9283.99 , 0.001
11+
9 , -2179.87 , 0.001
12+
10 , 5520.63 , 0.001
13+
11 , 12726.0 , 0.001
14+
12 , 17538.0 , 0.001
15+
13 , 17441.1 , 0.001
16+
14 , 10335.6 , 0.001
17+
15 , -2686.37 , 0.001
18+
16 , -10522.5 , 0.001
19+
17 , 24712.2 , 0.001
20+
18 , 166438.0 , 0.001
21+
19 , -15618.6 , 0.001
22+
20 , -14178.7 , 0.001
23+
21 , -12765.5 , 0.001
24+
22 , -12007.2 , 0.001
25+
23 , -13597.5 , 0.001
26+
24 , -19002.9 , 0.001
27+
25 , -28729.6 , 0.001
28+
26 , -41946.6 , 0.001
29+
27 , -56289.3 , 0.001
30+
28 , -67832.1 , 0.001
31+
29 , -71484.0 , 0.001
32+
30 , -62334.5 , 0.001
33+
31 , -38478.2 , 0.001
34+
32 , -4757.34 , 0.001
35+
33 , 26448.5 , 0.001
36+
34 , 45049.5 , 0.001
37+
35 , 60960.9 , 0.001
38+
36 , 83515.9 , 0.001
39+
37 , 8837.4 , 0.001

‎TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
11
VARIABLES="VARIABLE" , "DRAG" , "EFFICIENCY" , "FORCE_X" , "FORCE_Y" , "FORCE_Z" , "LIFT" , "MOMENT_X" , "MOMENT_Y" , "MOMENT_Z" , "SIDEFORCE"
2-
0 , 0.2405003164 , -108.3397202 , 0.2752804375 , -1.591326829 , 0.0 , -1.596953347 , 0.0 , 0.0 , 1.18511582 , 0.0
3-
1 , 0.4163088272 , -176.5661566 , 0.4692829784 , -2.423229574 , 0.0 , -2.432890265 , 0.0 , 0.0 , 1.020271816 , 0.0
4-
2 , 0.5558397416 , -231.7769399 , 0.6239540061 , -3.115569035 , 0.0 , -3.128439099 , 0.0 , 0.0 , 0.6093229584 , 0.0
2+
0 , 0.24008251 , -117.3444057 , 0.2742430499 , -1.56293638 , 0.0 , -1.568547024 , 0.0 , 0.0 , 1.189018284 , 0.0
3+
1 , 0.4064005433 , -189.77779 , 0.4586830911 , -2.391641926 , 0.0 , -2.401078899 , 0.0 , 0.0 , 1.030793484 , 0.0
4+
2 , 0.5421052294 , -249.5397676 , 0.6095878335 , -3.086770277 , 0.0 , -3.099333798 , 0.0 , 0.0 , 0.6218682473 , 0.0

‎TestCases/cont_adj_euler/wedge/of_grad_combo.dat.ref‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
VARIABLES="VARIABLE" , "GRADIENT" , "FINDIFF_STEP"
2-
0 , 0.00770904 , 0.0001
3-
1 , 0.00500468 , 0.0001
4-
2 , 0.00247269 , 0.0001
5-
3 , 0.000899035 , 0.0001
2+
0 , 0.00765473 , 0.0001
3+
1 , 0.00497838 , 0.0001
4+
2 , 0.0024697 , 0.0001
5+
3 , 0.00090216 , 0.0001

‎TestCases/disc_adj_fsi/Airfoil_2d/configFlow.cfg‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -51,10 +51,9 @@ TIME_DISCRE_FLOW= EULER_IMPLICIT
5151
% Linear solvers ------------------------------------------------------- %
5252
LINEAR_SOLVER= FGMRES
5353
LINEAR_SOLVER_PREC= ILU
54-
LINEAR_SOLVER_ERROR= 1E-3
55-
LINEAR_SOLVER_ITER= 20
56-
LINEAR_SOLVER_SMOOTHER_RELAXATION= 0.7
57-
DISCADJ_LIN_SOLVER= SMOOTHER
54+
LINEAR_SOLVER_ERROR= 1E-4
55+
LINEAR_SOLVER_ITER= 50
56+
DISCADJ_LIN_SOLVER= FGMRES
5857
DISCADJ_LIN_PREC= ILU
5958
NEWTON_KRYLOV= YES
6059
QUASI_NEWTON_NUM_SAMPLES= 999
@@ -91,4 +90,4 @@ VOLUME_ADJ_FILENAME= adjoint_fluid
9190
SURFACE_FILENAME= surface_fluid
9291
SURFACE_ADJ_FILENAME= adjoint_surface_fluid
9392
SCREEN_WRT_FREQ_INNER= 10
94-
CONV_FILENAME= history
93+
CONV_FILENAME= history

‎TestCases/disc_adj_fsi/dyn_fsi/config.cfg‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -11,16 +11,15 @@ MESH_FILENAME= mesh.su2
1111

1212
TIME_DOMAIN = YES
1313
TIME_ITER = 3
14-
TIME_STEP = 0.01
14+
TIME_STEP = 0.002
1515
UNST_ADJOINT_ITER= 10
1616
ITER_AVERAGE_OBJ= 3
1717

1818
OBJECTIVE_FUNCTION= REFERENCE_NODE
1919

20-
SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0]
20+
SCREEN_OUTPUT= TIME_ITER, OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], SENSITIVITY[1]
2121

2222
WRT_ZONE_HIST=YES
2323
%WRT_ZONE_CONV=YES
2424

25-
% Debug: Force Zero Grid Velocity = NO
2625
WRT_PERFORMANCE= NO

‎TestCases/disc_adj_fsi/dyn_fsi/configFEA.cfg‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ DESIGN_VARIABLE_FEA= YOUNG_MODULUS
1717

1818
% Solid properties ----------------------------------------------------- %
1919
MATERIAL_MODEL= NEO_HOOKEAN
20-
ELASTICITY_MODULUS= 7E8
20+
ELASTICITY_MODULUS= 2E9
2121
POISSON_RATIO= 0.35
2222
MATERIAL_DENSITY= 2700.0
2323

@@ -48,7 +48,6 @@ INCREMENTAL_LOAD= NO
4848
MESH_FILENAME= mesh.su2
4949
MESH_FORMAT= SU2
5050

51-
RESTART_SOL= NO
5251
SOLUTION_FILENAME= solution_solid.dat
5352
SOLUTION_ADJ_FILENAME= adjoint_solid.dat
5453
OUTPUT_WRT_FREQ= 1

‎TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -46,12 +46,15 @@ NUM_METHOD_GRAD= GREEN_GAUSS
4646
CFL_NUMBER= 1000.0
4747

4848
% Flow numerics -------------------------------------------------------- %
49-
CONV_NUM_METHOD_FLOW= JST
49+
CONV_NUM_METHOD_FLOW= ROE
50+
MUSCL_FLOW= YES
51+
SLOPE_LIMITER_FLOW= NONE
5052
JST_SENSOR_COEFF= ( 0.5, 0.02 )
53+
ENTROPY_FIX_COEFF= 0.01
5154
TIME_DISCRE_FLOW= EULER_IMPLICIT
5255

5356
% Linear solvers ------------------------------------------------------- %
54-
LINEAR_SOLVER= BCGSTAB
57+
LINEAR_SOLVER= FGMRES
5558
LINEAR_SOLVER_PREC= ILU
5659
LINEAR_SOLVER_ERROR= 1E-6
5760
LINEAR_SOLVER_ITER= 25
@@ -82,7 +85,6 @@ CONV_RESIDUAL_MINVAL= -7
8285
MESH_FILENAME= mesh.su2
8386
MESH_FORMAT= SU2
8487
%
85-
RESTART_SOL= NO
8688
SOLUTION_FILENAME= solution_fluid.dat
8789
SOLUTION_ADJ_FILENAME= adjoint_fluid.dat
8890
OUTPUT_WRT_FREQ= 1
@@ -95,5 +97,6 @@ SURFACE_FILENAME= surface_fluid
9597
SURFACE_ADJ_FILENAME= adjoint_surface_fluid
9698
SCREEN_WRT_FREQ_INNER= 10
9799
CONV_FILENAME= history
100+
SCREEN_OUTPUT= ITER, RMS_RES, LINSOL
98101
HISTORY_OUTPUT= ITER, RMS_RES, AERO_COEFF, TAVG_AERO_COEFF
99102
OUTPUT_FILES= NONE

‎TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
INDEX GRAD
2-
0 -3.533163525398154e-03
3-
1 -1.879941699367452e-03
4-
2 -8.078436962531523e-04
5-
3 -2.794845540524826e-04
6-
4 -2.790760674146328e-04
7-
5 -8.045490756184532e-04
8-
6 -1.867166054373919e-03
9-
7 -3.498473613679959e-03
2+
0 -4.570869215186209e-04
3+
1 -2.401466750293265e-04
4+
2 -9.134698389661222e-05
5+
3 -1.628087055796442e-05
6+
4 -1.741052457131011e-05
7+
5 -9.462489787924037e-05
8+
6 -2.452466259888652e-04
9+
7 -4.635483632874935e-04

0 commit comments

Comments
 (0)