UEqn Laplacian

Author: xu-23

src/dfGenMatrixIncompressible/GenMatrix_U.C

@@ -58,7 +58,7 @@ GenMatrix_U(
     // // fvm.upper() = fvm.lower() + rhoPhi.primitiveField();
     // fvm.upper() = -weights.primitiveField()*rhoPhi.primitiveField() + rhoPhi.primitiveField();
     // fvm.negSumDiag();   
-    // fvm.diag() += rDeltaT*rho.primitiveField()*mesh.Vsc(); // 再加上瞬态项
+    // fvm.diag() += rDeltaT*rho.primitiveField()*mesh.Vsc(); 
     // fvm.source() = rDeltaT
     //     *rho.oldTime().primitiveField()
     //     *U.oldTime().primitiveField()*mesh.Vsc();
@@ -95,11 +95,11 @@ GenMatrix_U(
         diagPtr[celli] = 0.0; // vector::zero;
     }
 
-    // 然后从 off-diagonal 项累积到对角线
+    // 
     for (label facei = 0; facei < nFaces; facei++)
     {
-        diagPtr[l[facei]] -= lowerPtr[facei];  // 从下三角累积
-        diagPtr[u[facei]] -= upperPtr[facei];  // 从上三角累积
+        diagPtr[l[facei]] -= lowerPtr[facei];  
+        diagPtr[u[facei]] -= upperPtr[facei]; 
     }
 
     for (label celli = 0; celli < nCells; celli++)
@@ -125,10 +125,142 @@ GenMatrix_U(
         fvm += fvc::surfaceIntegrate(rhoPhi*gcs.interpScheme().correction(U));
     }
 
-    // -------------------------------------------------------
+    // --------------------------------------------------------------------------------------------------------------
+
+    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];
+        const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
+        const fvsPatchScalarField& pDeltaCoeffs = deltaCoeffs.boundaryField()[patchi];
+        
+        if (pvf.coupled())
+        {
+            fvm_laplace.internalCoeffs()[patchi] =
+                pGamma * pvf.gradientInternalCoeffs(pDeltaCoeffs);
+            fvm_laplace.boundaryCoeffs()[patchi] =
+               -pGamma * pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
+        }
+        else
+        {
+            fvm_laplace.internalCoeffs()[patchi] = pGamma * pvf.gradientInternalCoeffs();
+            fvm_laplace.boundaryCoeffs()[patchi] = -pGamma * pvf.gradientBoundaryCoeffs();
+        }
+    }
+
+    // 非正交修正,否则不需要
+    if (mesh.fluxRequired(U.name()))
+    {
+        fvm_laplace.faceFluxCorrectionPtr() = new
+        GeometricField<vector, fvsPatchField, surfaceMesh>
+        (
+            gammaMagSf * tsnGradScheme_().correction(U)
+        );
+        
+        fvm_laplace.source() -=
+            mesh.V() *
+            fvc::div
+            (
+                *fvm_laplace.faceFluxCorrectionPtr()
+            )().primitiveField();
+    }
+    else
+    {
+        fvm_laplace.source() -=
+            mesh.V() *
+            fvc::div
+            (
+                gammaMagSf * tsnGradScheme_().correction(U)
+            )().primitiveField();
+    }
+
+    fvm -= fvm_laplace;
+
+    // --------------------------------------------------------------------------------------------------------------
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+    
+    
     // -------------------------------------------------------
     // interFoam    
     // const alphaField& alpha;
@@ -141,7 +273,8 @@ GenMatrix_U(
         // + MRF.DDt(rho, U) 
         // + turbulence.divDevRhoReff(rho, U)  
         - fvc::div((turbulence.alpha()*rho*turbulence.nuEff())*dev2(T(fvc::grad(U))))
-        - fvm::laplacian(turbulence.alpha()*rho*turbulence.nuEff(), U)
+        
+        // - fvm::laplacian(turbulence.alpha()*rho*turbulence.nuEff(), U)
         // - fvOptions(rho, U)
     );
 

src/dfGenMatrixIncompressible/gaussGrad.C

@@ -323,4 +323,41 @@ gaussGradSchemeGrad
 
 } // End namespace Foam
 
-// ************************************************************************* //
+// ************************************************************************* //
+
+// template<class Type>
+// tmp
+// <
+//     GeometricField
+//     <
+//         typename outerProduct<vector,Type>::type, fvPatchField, volMesh
+//     >
+// >
+// grad
+// (
+//     const GeometricField<Type, fvPatchField, volMesh>& vf
+// )
+// {
+//     return fvc::grad(vf, "grad(" + vf.name() + ')');
+// }
+
+// template<class Type>
+// tmp
+// <
+//     GeometricField
+//     <
+//         typename outerProduct<vector,Type>::type, fvPatchField, volMesh
+//     >
+// >
+// grad
+// (
+//     const GeometricField<Type, fvPatchField, volMesh>& vf,
+//     const word& name
+// )
+// {
+//     return fv::gradScheme<Type>::New
+//     (
+//         vf.mesh(),
+//         vf.mesh().gradScheme(name)
+//     )().grad(vf, name);
+// }