# gaussLaplacianScheme(3/3): Adding source terms – fvOptions

Here, we continue the modification of the solver myThermalConductionSolver using fvOptions. Please see the previous post

su =1, sp = 0 case

• If we keep su=1 in constant/tansportProperties, the source is modified from 9(-150 -80 -80 -70 0 0 -130 -60 -60) to 9(-150.01 -80.01 -80.01 -70.01 -0.01 -0.01 -130.01 -60.01 -60.01).
• su is multiplied by volume 0.01(0.1m*0.1m*1m here) and added to the source term.upper, lower and diag remain same.
• Temp is 9(371.854 380.969 382.795 350.053 350.096 350.11 328.217 319.151 317.34);

su = 0, sp =1 case

• Only diag value before adding bc effect changes from 9(-0.2 -0.3 -0.2 -0.3 -0.4 -0.3 -0.2 -0.3 -0.2) to 9(-0.19 -0.29 -0.19 -0.29 -0.39 -0.29 -0.19 -0.29 -0.19)
• This is because sp is +ve and it is added to the fvMatrix after multiplying with volume 0.01
• upper, lower and source remain same
• T = 9(385.473 404.352 409.393 369.939 386.458 392.282 340.771 340.612 341.768);

su = 0, sp = -1 case

• Only diag value before adding bc effect changes from 9(-0.2 -0.3 -0.2 -0.3 -0.4 -0.3 -0.2 -0.3 -0.2) to 9(-0.19 -0.29 -0.19 -0.29 -0.39 -0.29 -0.19 -0.29 -0.19)
• This is because sp is +ve and it is added to the fvMatrix after multiplying with volume 0.01
• upper, lower and source remain same
• T = 9(385.473 404.352 409.393 369.939 386.458 392.282 340.771 340.612 341.768);

Explanation in cfd-online by Alexey Matveichev1 fvOptions is an object of Foam::fv::optionList class. It has operator ():

            //- Return source for equation
template<class Type>
tmp<fvMatrix<Type> > operator()
(
GeometricField<Type, fvPatchField, volMesh>& fld
);

which in turn is implemented like this:

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
(
GeometricField<Type, fvPatchField, volMesh>& fld
)
{
return this->operator()(fld, fld.name());
}

template<class Type>
Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
(
GeometricField<Type, fvPatchField, volMesh>& fld,
const word& fieldName
)
{
checkApplied();
const dimensionSet ds = fld.dimensions()/dimTime*dimVolume;
tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(fld, ds));
fvMatrix<Type>& mtx = tmtx();

forAll(*this, i)
{
option& source = this->operator[](i);
label fieldI = source.applyToField(fieldName);
if (fieldI != -1)
{
source.setApplied(fieldI);
if (source.isActive())
{
if (debug)
{
Info<< "Applying source " << source.name() << " to field "
<< fieldName << endl;
}
}
}
}

return tmtx;
}

The awesome explanation

Explanation in cfd-onlie by Zeppo(Sergei) 2

I would like to reproduce what is told in the above post !!

Change the myThermalConductionSolver to accept fvOptions. A new solver is made myTCSfvOptions

Used the same case 03su0sp1 case but removed su ad sp entry from constant/transportProperties and added fvOptions in system. See the case here.

fvOptions with su = 0, sp =1

• diag from 9(-0.19 -0.29 -0.19 -0.29 -0.39 -0.29 -0.19 -0.29 -0.19) in non-fvOptions case to 9(-0.21 -0.31 -0.21 -0.31 -0.41 -0.31 -0.21 -0.31 -0.21) in fvOptions case
• This is because sp is added in the rhs which become -ve when taken into lhs. Hence it is same as the case su0sp-1

Now we can try with sp =-1

fvOptions with su = 0, sp =-1

• diagonal elements not affected. diag = 9(-0.2 -0.3 -0.2 -0.3 -0.4 -0.3 -0.2 -0.3 -0.2)
• source changed to 9(-153 -83 -83 -73 -3 -3 -133 -63 -63)
• i.e., since sp is -ve, explicit evaluation is carried out but it is not same as su=1 case, which had source = 9(-150.01 -80.01 -80.01 -70.01 -0.01 -0.01 -130.01 -60.01 -60.01). What is the difference? It is sp*vol*T = 1*0.01*300 = 3 is added to source.
• So, in the earlier case, su is supposed to be the constant value in source linearisation (which already have T in it. Say Su =AV(Tinf^4+3Tc^4) for radiation)
• This sort of answers this question