Skip to content

Commit 79c8e39

Browse files
committed
GenMatrix_p
1 parent 74123a4 commit 79c8e39

File tree

7 files changed

+559
-103
lines changed

7 files changed

+559
-103
lines changed

src/dfGenMatrixIncompressible/GenFvMatrix.H

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,43 @@ gaussGradSchemeGrad
7575
);
7676

7777

78+
// fvc::div
79+
template<class Type>
80+
tmp<GeometricField<Type, fvPatchField, volMesh>>
81+
gaussConvectionSchemeFvcDiv
82+
(
83+
const surfaceScalarField& faceFlux,
84+
const GeometricField<Type, fvPatchField, volMesh>& vf
85+
);
86+
87+
template<class Type>
88+
tmp<GeometricField<Type, fvPatchField, volMesh>>
89+
gaussConvectionSchemeFvcDiv
90+
(
91+
const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
92+
);
93+
94+
template<class Type>
95+
tmp<GeometricField<Type, fvPatchField, volMesh>>
96+
gaussConvectionSchemeFvcDiv
97+
(
98+
const tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>& tssf
99+
);
100+
101+
template<class Type>
102+
tmp
103+
<
104+
GeometricField
105+
<
106+
typename innerProduct<vector, Type>::type, fvPatchField, volMesh
107+
>
108+
>
109+
gaussDivFvcdiv
110+
(
111+
const GeometricField<Type, fvPatchField, volMesh>& vf
112+
);
113+
114+
78115
// fvm::laplacian
79116

80117
template<class Type>
@@ -143,6 +180,14 @@ GenMatrix_U(
143180
incompressible::turbulenceModel& turbulence
144181
);
145182

183+
template<class Type>
184+
tmp<fvMatrix<Type>>
185+
GenMatrix_p(
186+
const GeometricField<scalar, fvsPatchField, surfaceMesh>& rAUf,
187+
const GeometricField<Type, fvPatchField, volMesh>& p_rgh,
188+
const tmp<GeometricField<double, fvsPatchField, surfaceMesh>>& phiHbyA
189+
);
190+
146191
void check_fvmatrix_equal(const fvVectorMatrix& answer, const fvVectorMatrix& check, const word& name);
147192
void check_fvmatrix_equal(const fvScalarMatrix& answer, const fvScalarMatrix& check, const word& name);
148193

src/dfGenMatrixIncompressible/GenMatrix_U.C

Lines changed: 16 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -45,18 +45,18 @@ GenMatrix_U(
4545
rho.dimensions()*U.dimensions()*dimVol/dimTime
4646
)
4747
);
48-
fvVectorMatrix& fvm_DDT = tfvm_DDT.ref();
48+
fvVectorMatrix& fvm = tfvm_DDT.ref();
4949

5050
scalar rDeltaT = 1.0/mesh.time().deltaTValue();
5151

5252
// interField
53-
fvm_DDT.lower() = -weights.primitiveField()*rhoPhi.primitiveField();
54-
// fvm_DDT.upper() = fvm_DDT.lower() + rhoPhi.primitiveField();
55-
fvm_DDT.upper() = -weights.primitiveField()*rhoPhi.primitiveField() + rhoPhi.primitiveField();
56-
fvm_DDT.negSumDiag(); // diag[i] = - (sum_of_lower_coeffs[i] + sum_of_upper_coeffs[i]) 先设置对流项
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]) 先设置对流项
5757

58-
fvm_DDT.diag() += rDeltaT*rho.primitiveField()*mesh.Vsc(); // 再加上瞬态项
59-
fvm_DDT.source() = rDeltaT
58+
fvm.diag() += rDeltaT*rho.primitiveField()*mesh.Vsc(); // 再加上瞬态项
59+
fvm.source() = rDeltaT
6060
*rho.oldTime().primitiveField()
6161
*U.oldTime().primitiveField()*mesh.Vsc();
6262

@@ -67,25 +67,25 @@ GenMatrix_U(
6767
const fvsPatchScalarField& patchFlux = rhoPhi.boundaryField()[patchi];
6868
const fvsPatchScalarField& pw = weights.boundaryField()[patchi];
6969

70-
fvm_DDT.internalCoeffs()[patchi] = patchFlux*psf.valueInternalCoeffs(pw);
71-
fvm_DDT.boundaryCoeffs()[patchi] = -patchFlux*psf.valueBoundaryCoeffs(pw);
70+
fvm.internalCoeffs()[patchi] = patchFlux*psf.valueInternalCoeffs(pw);
71+
fvm.boundaryCoeffs()[patchi] = -patchFlux*psf.valueBoundaryCoeffs(pw);
7272
}
7373

7474
// correct
7575
if (gcs.interpScheme().corrected())
7676
{
77-
fvm_DDT += fvc::surfaceIntegrate(rhoPhi*gcs.interpScheme().correction(U));
77+
fvm += fvc::surfaceIntegrate(rhoPhi*gcs.interpScheme().correction(U));
7878
}
7979

8080

8181

82-
// scalar* __restrict__ diagPtr_ddt = fvm_DDT.diag().begin();
83-
// scalar* __restrict__ sourcePtr_ddt = fvm_DDT.source().begin();
84-
// scalar* __restrict__ lowerPtr_ddt = fvm_DDT.lower().begin();
85-
// scalar* __restrict__ upperPtr_ddt = fvm_DDT.upper().begin();
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();
8686

87-
// const labelUList& l = fvm_DDT.lduAddr().lowerAddr();
88-
// const labelUList& u = fvm_DDT.lduAddr().upperAddr();
87+
// const labelUList& l = fvm.lduAddr().lowerAddr();
88+
// const labelUList& u = fvm.lduAddr().upperAddr();
8989

9090
// interFoam
9191

Lines changed: 100 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,10 @@
11
#include "GenFvMatrix.H"
2-
#include "multivariateGaussConvectionScheme.H"
3-
#include "gaussConvectionScheme.H"
2+
3+
#include "gaussLaplacianScheme.H"
4+
#include "surfaceInterpolate.H"
5+
#include "fvcDiv.H"
6+
#include "fvcGrad.H"
7+
#include "fvMatrices.H"
48
#include "snGradScheme.H"
59
#include "linear.H"
610
#include "orthogonalSnGrad.H"
@@ -15,69 +19,122 @@
1519

1620
namespace Foam{
1721

18-
1922
template<class Type>
20-
tmp
21-
<
22-
GeometricField
23-
<
24-
typename outerProduct<vector, Type>::type,
25-
fvPatchField,
26-
volMesh
27-
>
28-
>
23+
tmp<fvMatrix<Type>>
2924
GenMatrix_p(
30-
const GeometricField<Type, fvPatchField, volMesh>& vsf
25+
const GeometricField<scalar, fvsPatchField, surfaceMesh>& rAUf,
26+
const GeometricField<Type, fvPatchField, volMesh>& p_rgh,
27+
const tmp<GeometricField<double, fvsPatchField, surfaceMesh>>& phiHbyA
3128
){
32-
word name("grad(" + vsf.name() + ')');
33-
34-
const fvMesh& mesh = U.mesh();
35-
36-
37-
38-
39-
40-
41-
42-
43-
44-
45-
46-
47-
48-
49-
50-
51-
52-
53-
54-
55-
56-
57-
58-
59-
60-
29+
Info << "pEqn gaussLaplacianSchemeFvmLaplacian" << endl;
30+
31+
const fvMesh& mesh = p_rgh.mesh();
32+
tmp<fv::snGradScheme<Type>> tsnGradScheme_(new fv::orthogonalSnGrad<Type>(mesh));
6133

34+
GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf
35+
(
36+
rAUf*mesh.magSf()
37+
);
6238

39+
const surfaceScalarField& deltaCoeffs = tsnGradScheme_().deltaCoeffs(p_rgh);
6340

41+
// -----
42+
tmp<fvMatrix<Type>> tfvm_Lap
43+
(
44+
new fvScalarMatrix
45+
(
46+
p_rgh,
47+
deltaCoeffs.dimensions()*gammaMagSf.dimensions()*p_rgh.dimensions()
48+
)
49+
);
6450

51+
fvScalarMatrix& fvm = tfvm_Lap.ref();
52+
53+
fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
54+
fvm.negSumDiag();
55+
56+
forAll(p_rgh.boundaryField(), patchi)
57+
{
58+
const fvPatchField<Type>& pvf = p_rgh.boundaryField()[patchi];
59+
const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
60+
const fvsPatchScalarField& pDeltaCoeffs =
61+
deltaCoeffs.boundaryField()[patchi];
62+
63+
if (pvf.coupled())
64+
{
65+
fvm.internalCoeffs()[patchi] =
66+
pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
67+
fvm.boundaryCoeffs()[patchi] =
68+
-pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
69+
}
70+
else
71+
{
72+
fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
73+
fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
74+
}
75+
}
76+
77+
// -----
78+
79+
if (tsnGradScheme_().corrected())
80+
{
81+
if (mesh.fluxRequired(p_rgh.name()))
82+
{
83+
fvm.faceFluxCorrectionPtr() = new
84+
GeometricField<Type, fvsPatchField, surfaceMesh>
85+
(
86+
gammaMagSf*tsnGradScheme_().correction(p_rgh)
87+
);
88+
89+
fvm.source() -=
90+
mesh.V()*
91+
fvc::div
92+
(
93+
*fvm.faceFluxCorrectionPtr()
94+
)().primitiveField();
95+
}
96+
else
97+
{
98+
fvm.source() -=
99+
mesh.V()*
100+
fvc::div
101+
(
102+
gammaMagSf*tsnGradScheme_().correction(p_rgh)
103+
)().primitiveField();
104+
}
105+
}
106+
107+
// -------------------------------------------------------
108+
109+
110+
111+
// interFoam pEqn.H
65112

66113
tmp<fvScalarMatrix> tfvm
67114
(
68115
new fvScalarMatrix
69116
(
70-
(fvm::laplacian(rAUf, p_rgh))
117+
(
118+
tfvm_Lap
119+
// + fvm::laplacian(rAUf, p_rgh)
120+
)
71121
==
72122
(fvc::div(phiHbyA))
73123
)
74124
);
75125

76126
return tfvm;
77127

128+
}
78129

79130

80131
}
81132

82133

134+
namespace Foam {
135+
template tmp<fvMatrix<scalar>> GenMatrix_p<scalar>(
136+
const GeometricField<scalar, fvsPatchField, surfaceMesh>&,
137+
const GeometricField<scalar, fvPatchField, volMesh>&,
138+
const tmp<GeometricField<double, fvsPatchField, surfaceMesh>>&
139+
);
83140
}

src/dfGenMatrixIncompressible/Utils.C

Lines changed: 22 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -124,31 +124,31 @@ void check_field_boundary_equal(const volVectorField& answer, const volVectorFie
124124
}
125125
}
126126

127-
// void check_field_boundary_equal(const surfaceScalarField& answer, const surfaceScalarField& check, const word& name){
128-
// check_field_error(answer, check, "field");
129-
// forAll(answer.boundaryField(), patchi)
130-
// {
131-
// check_field_error(answer.boundaryField()[patchi], check.boundaryField()[patchi], name + " boundaryField_" + std::to_string(patchi));
132-
// }
133-
// }
127+
void check_field_boundary_equal(const surfaceScalarField& answer, const surfaceScalarField& check, const word& name){
128+
check_field_error(answer, check, "field");
129+
forAll(answer.boundaryField(), patchi)
130+
{
131+
check_field_error(answer.boundaryField()[patchi], check.boundaryField()[patchi], name + " boundaryField_" + std::to_string(patchi));
132+
}
133+
}
134134

135-
// void check_fvmatrix_equal(const fvScalarMatrix& answer,const fvScalarMatrix& check, const word& name){
136-
// if(answer.source().begin() == check.source().begin()){
137-
// SeriousError << "answer and check are the same matrix." << endl;
138-
// MPI_Abort(PstreamGlobals::MPI_COMM_FOAM, -1);
139-
// }
135+
void check_fvmatrix_equal(const fvScalarMatrix& answer,const fvScalarMatrix& check, const word& name){
136+
// if(answer.source().begin() == check.source().begin()){
137+
// SeriousError << "answer and check are the same matrix." << endl;
138+
// MPI_Abort(PstreamGlobals::MPI_COMM_FOAM, -1);
139+
// }
140140

141-
// check_field_error(answer.source(),check.source(), name + " source");
142-
// check_field_error(answer.diag(),check.diag(), name + " diag");
143-
// check_field_error(answer.lower(),check.lower(), name + " lower");
144-
// check_field_error(answer.upper(),check.upper(), name + " upper");
141+
check_field_error(answer.source(),check.source(), name + " source");
142+
check_field_error(answer.diag(),check.diag(), name + " diag");
143+
check_field_error(answer.lower(),check.lower(), name + " lower");
144+
check_field_error(answer.upper(),check.upper(), name + " upper");
145145

146-
// for(label patchi = 0; patchi < const_cast<fvScalarMatrix&>(answer).internalCoeffs().size(); ++patchi)
147-
// {
148-
// check_field_error(const_cast<fvScalarMatrix&>(answer).internalCoeffs()[patchi],const_cast<fvScalarMatrix&>(check).internalCoeffs()[patchi], name + " internalCoeffs_" + std::to_string(patchi));
149-
// check_field_error(const_cast<fvScalarMatrix&>(answer).boundaryCoeffs()[patchi],const_cast<fvScalarMatrix&>(check).boundaryCoeffs()[patchi], name + " boundaryCoeffs_" + std::to_string(patchi));
150-
// }
151-
// }
146+
for(label patchi = 0; patchi < const_cast<fvScalarMatrix&>(answer).internalCoeffs().size(); ++patchi)
147+
{
148+
check_field_error(const_cast<fvScalarMatrix&>(answer).internalCoeffs()[patchi],const_cast<fvScalarMatrix&>(check).internalCoeffs()[patchi], name + " internalCoeffs_" + std::to_string(patchi));
149+
check_field_error(const_cast<fvScalarMatrix&>(answer).boundaryCoeffs()[patchi],const_cast<fvScalarMatrix&>(check).boundaryCoeffs()[patchi], name + " boundaryCoeffs_" + std::to_string(patchi));
150+
}
151+
}
152152

153153
// void check_fvmatrix_equal(const fvVectorMatrix& answer,const fvVectorMatrix& check, const word& name){
154154
// if(answer.source().begin() == check.source().begin()){

0 commit comments

Comments
 (0)