CFD Online Discussion Forums - How the boundary conditions are disposed into the matrix coefficient?

Hi the BC contributions (internalCoffs and boundaryCoeffs) to the matrix and RHS are calculated where the discrete version of the differential operators is assembled. These contributions are stored to be inserted into the matrix prior to the solution. In the case of a case of scalar transport and the div operator the lines are 98 in gaussConvectionScheme.C

Code:

00066 template<class Type>

00067 tmp<fvMatrix<Type> >

00068 gaussConvectionScheme<Type>::fvmDiv

00069 (

00070     const surfaceScalarField& faceFlux,

00071     const GeometricField<Type, fvPatchField, volMesh>& vf

00072 ) const

00073 {

00074     tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);

00075     const surfaceScalarField& weights = tweights();

00076 

00077     tmp<fvMatrix<Type> > tfvm

00078     (

00079         new fvMatrix<Type>

00080         (

00081             vf,

00082             faceFlux.dimensions()*vf.dimensions()

00083         )

00084     );

00085     fvMatrix<Type>& fvm = tfvm();

00086 

00087     fvm.lower() = -weights.internalField()*faceFlux.internalField();

00088     fvm.upper() = fvm.lower() + faceFlux.internalField();

00089     fvm.negSumDiag();

00090 

00091     forAll(vf.boundaryField(), patchI)

00092     {

00093         const fvPatchField<Type>& psf = vf.boundaryField()[patchI];

00094         const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchI];

00095         const fvsPatchScalarField& pw = weights.boundaryField()[patchI];

00096 

00097         fvm.internalCoeffs()[patchI] = patchFlux*psf.valueInternalCoeffs(pw);

00098         fvm.boundaryCoeffs()[patchI] = -patchFlux*psf.valueBoundaryCoeffs(pw);

00099     }

00100 

00101     if (tinterpScheme_().corrected())

00102     {

00103         fvm += fvc::surfaceIntegrate(faceFlux*tinterpScheme_().correction(vf));

00104     }

00105 

00106     return tfvm;

00107 }

and 106 and 109 in fvScalarMatrix.C

Code:

00095 template<>

00096 Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::fvSolver::solve

00097 (

00098     const dictionary& solverControls

00099 )

00100 {

00101     GeometricField<scalar, fvPatchField, volMesh>& psi =

00102         const_cast<GeometricField<scalar, fvPatchField, volMesh>&>

00103         (fvMat_.psi());

00104 

00105     scalarField saveDiag(fvMat_.diag());

00106     fvMat_.addBoundaryDiag(fvMat_.diag(), 0);

00107 

00108     scalarField totalSource(fvMat_.source());

00109     fvMat_.addBoundarySource(totalSource, false);

00110 

00111     // assign new solver controls

00112     solver_->read(solverControls);

00113 

00114     lduMatrix::solverPerformance solverPerf = solver_->solve

00115     (

00116         psi.internalField(),

00117         totalSource

00118     );

00119 

00120     solverPerf.print();

00121 

00122     fvMat_.diag() = saveDiag;

00123 

00124     psi.correctBoundaryConditions();

00125 

00126     psi.mesh().setSolverPerformance(psi.name(), solverPerf);

00127 

00128     return solverPerf;

00129 }

Regards.