哦，
diag() = saveDiag;//恢复diag
是在括号之前的，所以每次的修改不会叠加。
之前理解错误的一点是所有的边界有关的系数都是放在internalCoeffs和boundaryCoeffs中的，不仅仅是coupled BC。从初始化时候的size可以看出来。
template<class Type>
Foam::fvMatrix<Type>::fvMatrix
(
const GeometricField<Type, fvPatchField, volMesh>& psi,
const dimensionSet& ds
)
:
lduMatrix(psi.mesh()),
psi_(psi),
dimensions_(ds),
source_(psi.size(), Zero),
internalCoeffs_(psi.mesh().boundary().size()), //和边界的face数大小一样。
boundaryCoeffs_(psi.mesh().boundary().size()),
faceFluxCorrectionPtr_(nullptr)
{
//...
不过非耦合的应该是internalCoeffs放比例系数，boundaryCoeffs放源项，而耦合的BC是internalCoeffs放owner的系数，boundaryCoeffs放的另一侧单元的系数。从addBoundarySource()的源代码可以看出。
template<class Type>
void Foam::fvMatrix<Type>::addBoundarySource
(
Field<Type>& source,
const bool couples
) const
{
forAll(psi_.boundaryField(), patchi)
{
const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
const Field<Type>& pbc = boundaryCoeffs_[patchi];
if (!ptf.coupled())//非耦合边界，只有owner，另一侧没有单元，pbc存的是边界源项。
{
addToInternalField(lduAddr().patchAddr(patchi), pbc, source);
}
else if (couples)//耦合边界，如果couples==true, 另一侧有单元，pbc存的是另一侧单元的系数。
{
const tmp<Field<Type>> tpnf = ptf.patchNeighbourField();
const Field<Type>& pnf = tpnf();
const labelUList& addr = lduAddr().patchAddr(patchi);
forAll(addr, facei)
{
source[addr[facei]] += cmptMultiply(pbc[facei], pnf[facei]);
}
}
}
}
