# problems creating a sink term with scalarCodedSource

 Register Blogs Members List Search Today's Posts Mark Forums Read

June 6, 2018, 22:32
problems creating a sink term with scalarCodedSource
#1
New Member

Stephen Waite
Join Date: May 2013
Location: Auckland, New Zealand
Posts: 29
Rep Power: 13
Hi guys,

I am tyring to create a sink term that will remove some amount of a passive scalar depending on the amount of said scalar within the cell.

So I thought I would start of with the scalarCodedSource function within fvOptions. I have an example from online of this working (attached as buoyantPimpleFoam).

So I added a simple passive transport equation to a Piso solver for a volScalarField psi. This works fine for a scalarSemiImplicitSource, but when I change this to a scalarCodedSource (both options are in passiveScalarPisoFOAM attachment, along with the make files for the solver, works for OF 4.1 and 5), I start getting the following error.

Code:
```Courant Number mean: 0 max: 0
smoothSolver:  Solving for Ux, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
GAMG:  Solving for p, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
GAMG:  Solving for p, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0, global = 0, cumulative = 0
Using dynamicCode for fvOption:: sourceTime at line 51 in "/hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/constant/fvOptions.psiSource.scalarCodedSourceCoeffs"
Invoking "wmake -s libso /hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/dynamicCode/sourceTime"
wmake libso /hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/dynamicCode/sourceTime
ln: ./lnInclude
wmkdep: codedFvOptionTemplate.C
Ctoo: codedFvOptionTemplate.C
/hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/constant/fvOptions.psiSource.scalarCodedSourceCoeffs: In member function ‘virtual void Foam::fv::sourceTimeFvOptionscalarSource::addSup(Foam::fvMatrix<double>&, Foam::label)’:
/hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/constant/fvOptions.psiSource.scalarCodedSourceCoeffs:73:18: error: passing ‘const scalarField {aka const Foam::Field<double>}’ as ‘this’ argument of ‘void Foam::Field<Type>::operator-=(const Type&) [with Type = double]’ discards qualifiers [-fpermissive]
/hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/constant/fvOptions.psiSource.scalarCodedSourceCoeffs: In member function ‘virtual void Foam::fv::sourceTimeFvOptionscalarSource::addSup(const volScalarField&, Foam::fvMatrix<double>&, Foam::label)’:
/hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/constant/fvOptions.psiSource.scalarCodedSourceCoeffs:73:18: error: passing ‘const scalarField {aka const Foam::Field<double>}’ as ‘this’ argument of ‘void Foam::Field<Type>::operator-=(const Type&) [with Type = double]’ discards qualifiers [-fpermissive]
make: *** [Make/linux64GccDPInt32Opt/codedFvOptionTemplate.o] Error 1

--> FOAM FATAL IO ERROR:

file: /hpc_ntot/swai013/PostDoc/SimulationWork/passiveScalarPisoFOAM/Example/square/constant/fvOptions.psiSource.scalarCodedSourceCoeffs from line 51 to line 82.

From function void Foam::codedBase::createLibrary(Foam::dynamicCode&, const Foam::dynamicCodeContext&) const
in file db/dynamicLibrary/codedBase/codedBase.C at line 206.

FOAM exiting```
So the problem seems to be that it doesn’t like to take the volScalarField psi, which is strange because again scalarSemiImplicitSource is happy to use it.

Going back to the working case, the field used (h) is maybe a scalarField, not a volScalarField (I found it pretty tough going understanding the solver and where it actually creates h, but it looks like it is created in heThermo.C), and so that's the only difference I can think of.

Hopefully this is just a simple fix because of my poor coding skills, if anyone knows how to fix this, or knows of a better way to add in a sink term, that would be great!

Cheers
Stephen
Attached Files
 buoyantPimpleFoam.zip (127.4 KB, 27 views) passiveScalarPisoFOAM.zip (12.4 KB, 27 views)

 June 7, 2018, 03:40 #2 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 38 Hi, const specifier is, in general, for values you do not want to change. So this Code: ```const scalarField& psiS = eqn.source(); psiS -= sqr(time.value());``` is not very reasonable. In fact, compiler error message tells you about your logic flaw: Code: `error: passing ‘const scalarField {aka const Foam::Field}’ as ‘this’ argument of ‘void Foam::Field::operator-=(const Type&)` Just remove const before scalarField. Stephen Waite likes this.

 June 7, 2018, 05:35 #3 New Member   Stephen Waite Join Date: May 2013 Location: Auckland, New Zealand Posts: 29 Rep Power: 13 Fantastic that fixed that issue, although I don't understand why it could be a const field for the example case, but not for mine (although I'm happy to accept that it is). I dont know if you are familiar with this fvOptions...option but it looks like eqn.source() is grabbing the rhs vector for the whole domain, not just the selectionMode. Which is fine, but when I try and just apply source values to the cellZone Code: ```codeAddSup #{ //const Time& time = mesh().time(); //const scalarField& V = mesh_.V(); label cellZoneID = mesh_.cellZones().findZoneID("scalarBlob"); const labelList& cells = mesh_.cellZones()[cellZoneID]; scalarField& psiS = eqn.source(); forAll(cells, cellI) { psiS[cellI] -= 0.0001; } #};``` It sort of does something, it applies to cells at the bottom of the domain, not in the cellZone Any thoughts? Thanks for your help

 June 7, 2018, 06:16 #4 Senior Member   Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,938 Rep Power: 38 Hi, In example case const scalarField is volume of the mesh: Code: `const scalarField& V = mesh_.V();` equation source is non-constant reference: Code: `scalarField& heSource = eqn.source();` Your second problem has almost nothing to do with fvOptions. In this cycle: Code: ```forAll(cells, cellI) { psiS[cellI] -= 0.0001; }``` cellI is just counter, which goes from 0 to a size of cells list (so, it seems, numbering of cells in your mesh starts from bottom). You need real cell IDs, if you would like to add source in cell zone, to get them use cells list: Code: ```forAll(cells, i) { const label cellIdx = cells[i]; psiS[cellIdx] -= 0.0001; }``` Stephen Waite likes this.

 June 7, 2018, 06:26 #5 New Member   Stephen Waite Join Date: May 2013 Location: Auckland, New Zealand Posts: 29 Rep Power: 13 Haha dam I realised the cellI mistake when I was on the bus and was hoping to get home and change the post before you saw it!! And you are absolutely right about the scalarField I cant even copy and paste properly I cant believe I didnt notice that. I have spent the better part of the day trying to figure that out. Thanks again for your help your a life saver! Cheers Stephen

 Tags passivescalar, scalarcodedsource