U

Author: xu-23

src/dfGenMatrixIncompressible/GenFvMatrix.H

@@ -202,4 +202,6 @@ void check_fvmatrix_equal(const fvScalarMatrix& answer, const fvScalarMatrix& ch
 } // End namespace Foam
 
 
+
+
 // ************************************************************************* //

src/dfGenMatrixIncompressible/GenMatrix_U.C

@@ -19,6 +19,10 @@
 
 #include "gaussGrad.H"
 #include "extrapolatedCalculatedFvPatchField.H"
+#include <mpi.h>
+
+double UEqn_build_pre_time = 0., UEqn_build_item1_intern_time = 0., UEqn_build_item1_bound_time = 0., UEqn_build_item2_intern_time = 0., UEqn_build_item2_bound_time = 0., UEqn_build_item3_time = 0., UEqn_build_item3_last_time = 0;
+double startU, endU;
 
 namespace Foam{
 
@@ -252,6 +256,9 @@ GenMatrix_U(
     const volScalarField& p, 
     incompressible::turbulenceModel& turbulence
 ){
+    UEqn_build_pre_time = 0., UEqn_build_item1_intern_time = 0., UEqn_build_item1_bound_time = 0., UEqn_build_item2_intern_time = 0., UEqn_build_item2_bound_time = 0., UEqn_build_item3_time = 0., UEqn_build_item3_last_time = 0;
+    startU = MPI_Wtime();
+
     word name("div("+rhoPhi.name()+','+U.name()+')');
 
     const fvMesh& mesh = U.mesh();
@@ -276,8 +283,29 @@ GenMatrix_U(
     );
     fvVectorMatrix& fvm = tfvm_DDT.ref();
 
+    // -------------
     scalar rDeltaT = 1.0/mesh.time().deltaTValue();
 
+    Info << "UEqn gaussLaplacianSchemeFvmLaplacian" << endl;
+    
+    const volScalarField alphaEff = turbulence.alpha()*rho*turbulence.nuEff();
+    // tmp<fv::snGradScheme<vector>> tsnGradScheme_(new fv::orthogonalSnGrad<vector>(mesh));  // 正交修正方案
+    tmp<fv::snGradScheme<vector>> tsnGradScheme_(new fv::correctedSnGrad<vector>(mesh));  // 对应snGradSchemes{default         corrected;}
+
+    surfaceScalarField alphaEfff = fvc::interpolate(alphaEff);    ///  
+    
+    // gammaMagSf = alphaEfff * magSf
+    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf
+    (
+        alphaEfff * mesh.magSf()
+    );
+    
+    const surfaceScalarField& deltaCoeffs = tsnGradScheme_().deltaCoeffs(U);
+    // -------------
+
+    endU = MPI_Wtime();
+    UEqn_build_pre_time += endU - startU;
+
     // interField
     // fvm.lower() = -weights.primitiveField()*rhoPhi.primitiveField();
     // // fvm.upper() = fvm.lower() + rhoPhi.primitiveField();
@@ -289,6 +317,13 @@ GenMatrix_U(
     //     *U.oldTime().primitiveField()*mesh.Vsc();
 
 
+    // interField
+    // fvm_laplace.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
+    // fvm_laplace.negSumDiag();
+    // fvm -=fvm_laplace;
+
+    startU = MPI_Wtime();
+
     vector* __restrict__ sourcePtr = fvm.source().begin();
     scalar* __restrict__ diagPtr = fvm.diag().begin();  
     scalar* __restrict__ lowerPtr = fvm.lower().begin();
@@ -297,6 +332,9 @@ GenMatrix_U(
     const labelUList& l = fvm.lduAddr().lowerAddr();
     const labelUList& u = fvm.lduAddr().upperAddr();
 
+    const scalar* const __restrict__ gammaMagSfPtr = gammaMagSf.primitiveField().begin();
+    const scalar* const __restrict__ deltaCoeffsPtr = deltaCoeffs.primitiveField().begin();
+
     const scalar* const __restrict__ weightsPtr = weights.primitiveField().begin();
     const scalar* const __restrict__ rhoPhiPtr = rhoPhi.primitiveField().begin();
     const scalar* const __restrict__ rhoPtr = rho.primitiveField().begin();
@@ -308,31 +346,43 @@ GenMatrix_U(
     const label nFaces = fvm.lower().size();
     const label nCells = fvm.diag().size();
 
+    for (label celli = 0; celli < nCells; celli++)
+    {
+        diagPtr[celli] = 0.0;
+    }
+
     for (label facei = 0; facei < nFaces; facei++)
     {
         scalar flux = weightsPtr[facei] * rhoPhiPtr[facei];
         lowerPtr[facei] = -flux;
         upperPtr[facei] = -flux + rhoPhiPtr[facei];
+        
+        diagPtr[l[facei]] += flux;         
+        diagPtr[u[facei]] += flux - rhoPhiPtr[facei];  
     }
 
     for (label celli = 0; celli < nCells; celli++)
     {
-        diagPtr[celli] = 0.0; // vector::zero;
+        diagPtr[celli] += rDeltaT * rhoPtr[celli] * meshVscPtr[celli];
+        sourcePtr[celli] = rDeltaT * rhoOldTimePtr[celli] * UOldTimePtr[celli] * meshVscPtr[celli];
     }
 
-    // 
     for (label facei = 0; facei < nFaces; facei++)
     {
-        diagPtr[l[facei]] -= lowerPtr[facei];  
-        diagPtr[u[facei]] -= upperPtr[facei]; 
-    }
+        scalar coeffs = deltaCoeffsPtr[facei] * gammaMagSfPtr[facei];
 
-    for (label celli = 0; celli < nCells; celli++)
-    {
-        diagPtr[celli] += rDeltaT * rhoPtr[celli] * meshVscPtr[celli];
-        sourcePtr[celli] = rDeltaT * rhoOldTimePtr[celli] * UOldTimePtr[celli] * meshVscPtr[celli];
+        upperPtr[facei] -= coeffs;
+        lowerPtr[facei] -= coeffs;
+        
+        diagPtr[l[facei]] += coeffs; 
+        diagPtr[u[facei]] += coeffs;  
     }
 
+    endU = MPI_Wtime();
+    UEqn_build_item1_intern_time += endU - startU;
+
+    startU = MPI_Wtime();
+
     // boundaryField
     forAll(U.boundaryField(), patchi)
     {
@@ -350,77 +400,11 @@ GenMatrix_U(
         fvm += fvc::surfaceIntegrate(rhoPhi*gcs.interpScheme().correction(U));
     }
 
+    endU = MPI_Wtime();
+    UEqn_build_item1_bound_time += endU - startU;
     // --------------------------------------------------------------------------------------------------------------
+    startU = MPI_Wtime();
 
-    Info << "UEqn gaussLaplacianSchemeFvmLaplacian" << endl;
-     
-    const volScalarField alphaEff = turbulence.alpha()*rho*turbulence.nuEff();
-    // tmp<fv::snGradScheme<vector>> tsnGradScheme_(new fv::orthogonalSnGrad<vector>(mesh));  // 正交修正方案
-    tmp<fv::snGradScheme<vector>> tsnGradScheme_(new fv::correctedSnGrad<vector>(mesh));  // 对应snGradSchemes{default         corrected;}
-
-    surfaceScalarField alphaEfff = fvc::interpolate(alphaEff);    ///  
-    
-    // gammaMagSf = alphaEfff * magSf
-    GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf
-    (
-        alphaEfff * mesh.magSf()
-    );
-    
-    const surfaceScalarField& deltaCoeffs = tsnGradScheme_().deltaCoeffs(U);
-    // 
-    tmp<fvVectorMatrix> tfvm_Laplacian
-    (
-        new fvVectorMatrix
-        (
-            U,
-            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*U.dimensions()
-        )
-    );
-    fvVectorMatrix& fvm_laplace = tfvm_Laplacian.ref();
-    
-    // interField
-    // fvm_laplace.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
-    // fvm_laplace.negSumDiag();
-
-    vector* __restrict__ sourcePtrs = fvm_laplace.source().begin();
-    scalar* __restrict__ diagPtrs = fvm_laplace.diag().begin();  
-    scalar* __restrict__ lowerPtrs = fvm_laplace.lower().begin();
-    scalar* __restrict__ upperPtrs = fvm_laplace.upper().begin();
-
-    const labelUList& ls = fvm_laplace.lduAddr().lowerAddr();
-    const labelUList& us = fvm_laplace.lduAddr().upperAddr();
-
-    // 获取场数据指针
-    const scalar* const __restrict__ gammaMagSfPtrs = gammaMagSf.primitiveField().begin();
-    const scalar* const __restrict__ deltaCoeffsPtrs = deltaCoeffs.primitiveField().begin();
-    const scalar* const __restrict__ meshVPtrs = mesh.V().begin();
-
-    const label nFacess = fvm_laplace.lower().size();
-    const label nCellss = fvm_laplace.diag().size();
-
-    // 初始化矩阵
-    for (label celli = 0; celli < nCellss; celli++)
-    {
-        diagPtrs[celli] = 0.0;
-        sourcePtrs[celli] = vector::zero;
-    }
-
-    // 设置内部场矩阵系数
-    for (label facei = 0; facei < nFacess; facei++)
-    {
-        scalar coeffs = deltaCoeffsPtrs[facei] * gammaMagSfPtrs[facei];
-        upperPtrs[facei] = +coeffs;
-        lowerPtrs[facei] = +coeffs;
-    }
-
-    // 对角线求和
-    for (label facei = 0; facei < nFacess; facei++)
-    {
-        diagPtrs[ls[facei]] -= lowerPtrs[facei];  
-        diagPtrs[us[facei]] -= upperPtrs[facei]; 
-    }
-
-    // ------------------------------------
     forAll(U.boundaryField(), patchi)
     {
         const fvPatchVectorField& pvf = U.boundaryField()[patchi];
@@ -429,47 +413,70 @@ GenMatrix_U(
 
         if (pvf.coupled())
         {
-            fvm_laplace.internalCoeffs()[patchi] =
+            fvm.internalCoeffs()[patchi] -=
                 pGamma * pvf.gradientInternalCoeffs(pDeltaCoeffs);
-            fvm_laplace.boundaryCoeffs()[patchi] =
+            fvm.boundaryCoeffs()[patchi] -=
                -pGamma * pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
         }
         else
         {
-            fvm_laplace.internalCoeffs()[patchi] = pGamma * pvf.gradientInternalCoeffs();
-            fvm_laplace.boundaryCoeffs()[patchi] = -pGamma * pvf.gradientBoundaryCoeffs();
+            fvm.internalCoeffs()[patchi] -= pGamma * pvf.gradientInternalCoeffs();
+            fvm.boundaryCoeffs()[patchi] -= -pGamma * pvf.gradientBoundaryCoeffs();
         }
     }
 
-    // 非正交修正,否则不需要
     if (mesh.fluxRequired(U.name()))
     {
-        fvm_laplace.faceFluxCorrectionPtr() = new
-        GeometricField<vector, fvsPatchField, surfaceMesh>
+        tmp<GeometricField<vector, fvsPatchField, surfaceMesh>> tCorrection
         (
             gammaMagSf * tsnGradScheme_().correction(U)
         );
 
-        fvm_laplace.source() -=
-            mesh.V() *
-            fvc::div
-            (
-                *fvm_laplace.faceFluxCorrectionPtr()
-            )().primitiveField();
+        tmp<GeometricField<vector, fvPatchField, volMesh>> tDivCorrection
+        (
+            fvc::div(tCorrection())
+        );
+        
+        const vector* const __restrict__ divCorrectionPtrs = 
+            tDivCorrection().primitiveField().begin();
+        
+        for (label celli = 0; celli < nCells; celli++)
+        {
+            sourcePtr[celli] += meshVscPtr[celli] * divCorrectionPtrs[celli];
+        }
+        
+        fvm.faceFluxCorrectionPtr() = new
+        GeometricField<vector, fvsPatchField, surfaceMesh>
+        (
+            -tCorrection
+        );
     }
     else
     {
-        fvm_laplace.source() -=
-            mesh.V() *
-            fvc::div
-            (
-                gammaMagSf * tsnGradScheme_().correction(U)
-            )().primitiveField();
+        tmp<GeometricField<vector, fvsPatchField, surfaceMesh>> tCorrection
+        (
+            gammaMagSf * tsnGradScheme_().correction(U)
+        );
+        
+        tmp<GeometricField<vector, fvPatchField, volMesh>> tDivCorrection
+        (
+            fvc::div(tCorrection())
+        );
+        
+        const vector* const __restrict__ divCorrectionPtrs = 
+            tDivCorrection().primitiveField().begin();
+        
+        for (label celli = 0; celli < nCells; celli++)
+        {
+            sourcePtr[celli] += meshVscPtr[celli] * divCorrectionPtrs[celli];
+        }
     }
 
-    fvm -= fvm_laplace;
+    endU = MPI_Wtime();
+    UEqn_build_item2_bound_time += endU - startU;
 
     // --------------------------------------------------------------------------------------------------------------
+    startU = MPI_Wtime();
 
     Info << "UEqn Calculating grad(U) using extracted functions" << endl;
     tmp<surfaceInterpolationScheme<vector>> tinterpScheme =
@@ -485,9 +492,14 @@ GenMatrix_U(
     const volTensorField& gradU = tgGradU();
     // const volScalarField alphaEff = turbulence.alpha()*rho*turbulence.nuEff();
 
+    endU = MPI_Wtime();
+    UEqn_build_item3_time += endU - startU;
 
+    startU = MPI_Wtime();
 
     // -------------------------------------------------------
+
+    // -------------------------------------------------------  
     // interFoam    
     // const alphaField& alpha;
 
@@ -500,6 +512,7 @@ GenMatrix_U(
         // + turbulence.divDevRhoReff(rho, U)  
         // - fvc::div((turbulence.alpha()*rho*turbulence.nuEff())*dev2(T(fvc::grad(U))))
         - fvc::div((turbulence.alpha()*rho*turbulence.nuEff())*dev2(T(gradU)))
+        // - fvc::div((turbulence.alpha()*rho*turbulence.nuEff())*dev2T_gradU)
 
         // - fvm::laplacian(turbulence.alpha()*rho*turbulence.nuEff(), U)
         // - fvOptions(rho, U)
@@ -514,6 +527,20 @@ GenMatrix_U(
     //     -fvc::grad(p)
     // );     
 
+    endU = MPI_Wtime();
+    UEqn_build_item3_last_time += endU - startU;
+
+    Info<< "========Time Spent in UEqn build once=================="<< endl;
+    Info << "U pre Time :   " << UEqn_build_pre_time << endl;
+    Info << "U item1 intern Time : " << UEqn_build_item1_intern_time << endl;
+    Info << "U item1 bound Time : " << UEqn_build_item1_bound_time << endl;
+    // Info << "U item2 intern Time : " << UEqn_build_item2_intern_time << endl;
+    Info << "U item2 bound Time : " << UEqn_build_item2_bound_time << endl;
+    Info << "U item3 Time : " << UEqn_build_item3_time << endl;
+    Info << "U item3 last Time : " << UEqn_build_item3_last_time << endl;
+    Info<< "============================================"<<nl<< endl;
+
+    // return tfvm_DDT;
     return tfvm;
 }