1313#include "immiscibleIncompressibleTwoPhaseMixture.H"
1414#include "turbulentTransportModel.H"
1515
16+ #include "turbulenceModel.H"
17+ #include "autoPtr.H"
18+ #include "runTimeSelectionTables.H"
19+
1620namespace Foam {
1721
1822tmp < fvVectorMatrix >
@@ -50,15 +54,59 @@ GenMatrix_U(
5054 scalar rDeltaT = 1.0 /mesh .time ().deltaTValue ();
5155
5256 // interField
53- fvm .lower () = - weights .primitiveField ()* rhoPhi .primitiveField ();
54- // fvm.upper() = fvm.lower() + rhoPhi.primitiveField();
55- fvm .upper () = - weights .primitiveField ()* rhoPhi .primitiveField () + rhoPhi .primitiveField ();
56- fvm .negSumDiag (); // diag[i] = - (sum_of_lower_coeffs[i] + sum_of_upper_coeffs[i]) 先设置对流项
57-
58- fvm .diag () += rDeltaT * rho .primitiveField ()* mesh .Vsc (); // 再加上瞬态项
59- fvm .source () = rDeltaT
60- * rho .oldTime ().primitiveField ()
61- * U .oldTime ().primitiveField ()* mesh .Vsc ();
57+ // fvm.lower() = -weights.primitiveField()*rhoPhi.primitiveField();
58+ // // fvm.upper() = fvm.lower() + rhoPhi.primitiveField();
59+ // fvm.upper() = -weights.primitiveField()*rhoPhi.primitiveField() + rhoPhi.primitiveField();
60+ // fvm.negSumDiag();
61+ // fvm.diag() += rDeltaT*rho.primitiveField()*mesh.Vsc(); // 再加上瞬态项
62+ // fvm.source() = rDeltaT
63+ // *rho.oldTime().primitiveField()
64+ // *U.oldTime().primitiveField()*mesh.Vsc();
65+
66+
67+ vector * __restrict__ sourcePtr = fvm .source ().begin ();
68+ scalar * __restrict__ diagPtr = fvm .diag ().begin ();
69+ scalar * __restrict__ lowerPtr = fvm .lower ().begin ();
70+ scalar * __restrict__ upperPtr = fvm .upper ().begin ();
71+
72+ const labelUList & l = fvm .lduAddr ().lowerAddr ();
73+ const labelUList & u = fvm .lduAddr ().upperAddr ();
74+
75+ const scalar * const __restrict__ weightsPtr = weights .primitiveField ().begin ();
76+ const scalar * const __restrict__ rhoPhiPtr = rhoPhi .primitiveField ().begin ();
77+ const scalar * const __restrict__ rhoPtr = rho .primitiveField ().begin ();
78+ const scalar * const __restrict__ meshVscPtr = mesh .Vsc ()().begin ();
79+ const scalar * const __restrict__ rhoOldTimePtr = rho .oldTime ().primitiveField ().begin ();
80+
81+ const vector * const __restrict__ UOldTimePtr = U .oldTime ().primitiveField ().begin ();
82+
83+ const label nFaces = fvm .lower ().size ();
84+ const label nCells = fvm .diag ().size ();
85+
86+ for (label facei = 0 ; facei < nFaces ; facei ++ )
87+ {
88+ scalar flux = weightsPtr [facei ] * rhoPhiPtr [facei ];
89+ lowerPtr [facei ] = - flux ;
90+ upperPtr [facei ] = - flux + rhoPhiPtr [facei ];
91+ }
92+
93+ for (label celli = 0 ; celli < nCells ; celli ++ )
94+ {
95+ diagPtr [celli ] = 0.0 ; // vector::zero;
96+ }
97+
98+ // 然后从 off-diagonal 项累积到对角线
99+ for (label facei = 0 ; facei < nFaces ; facei ++ )
100+ {
101+ diagPtr [l [facei ]] -= lowerPtr [facei ]; // 从下三角累积
102+ diagPtr [u [facei ]] -= upperPtr [facei ]; // 从上三角累积
103+ }
104+
105+ for (label celli = 0 ; celli < nCells ; celli ++ )
106+ {
107+ diagPtr [celli ] += rDeltaT * rhoPtr [celli ] * meshVscPtr [celli ];
108+ sourcePtr [celli ] = rDeltaT * rhoOldTimePtr [celli ] * UOldTimePtr [celli ] * meshVscPtr [celli ];
109+ }
62110
63111 // boundaryField
64112 forAll (U .boundaryField (), patchi )
@@ -77,35 +125,26 @@ GenMatrix_U(
77125 fvm += fvc ::surfaceIntegrate (rhoPhi * gcs .interpScheme ().correction (U ));
78126 }
79127
128+ // -------------------------------------------------------
80129
130+
81131
82- // scalar* __restrict__ diagPtr_ddt = fvm.diag().begin();
83- // scalar* __restrict__ sourcePtr_ddt = fvm.source().begin();
84- // scalar* __restrict__ lowerPtr_ddt = fvm.lower().begin();
85- // scalar* __restrict__ upperPtr_ddt = fvm.upper().begin();
86-
87- // const labelUList& l = fvm.lduAddr().lowerAddr();
88- // const labelUList& u = fvm.lduAddr().upperAddr();
89-
132+ // -------------------------------------------------------
90133 // interFoam
134+ // const alphaField& alpha;
91135
92136 tmp < fvVectorMatrix > tfvm
93137 (
94- new fvVectorMatrix
95- (
96- (
97- tfvm_DDT
98- // + fvm::ddt(rho, U)
99- // + fvm::div(rhoPhi, U)
100- )
101- // + MRF.DDt(rho, U)
102-
103- ==
104- (- turbulence .divDevRhoReff (rho , U )) // fvOptions(rho, U)
105- )
138+ tfvm_DDT
139+ // + fvm::ddt(rho, U)
140+ // + fvm::div(rhoPhi, U)
141+ // + MRF.DDt(rho, U)
142+ // + turbulence.divDevRhoReff(rho, U)
143+ - fvc ::div ((turbulence .alpha ()* rho * turbulence .nuEff ())* dev2 (T (fvc ::grad (U ))))
144+ - fvm ::laplacian (turbulence .alpha ()* rho * turbulence .nuEff (), U )
145+ // - fvOptions(rho, U)
106146 );
107147
108-
109148 // @dfLowMachFoam
110149 // tmp<fvVectorMatrix> tUEqn
111150 // (
0 commit comments