[Stephen Waite]

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.

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"
Creating new library in "dynamicCode/sourceTime/platforms/linux64GccDPInt32Opt/lib/libsourceTime_f2075c0c6387b361c5451a93d1f438c9885adef6.so"
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: 
Failed wmake "dynamicCode/sourceTime/platforms/linux64GccDPInt32Opt/lib/libsourceTime_f2075c0c6387b361c5451a93d1f438c9885adef6.so"


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

[alexeym]

Hi,

const specifier is, in general, for values you do not want to change. So this

const scalarField& psiS = eqn.source();
psiS -= sqr(time.value());

is not very reasonable. In fact, compiler error message tells you about your logic flaw:

error: passing ‘const scalarField {aka const Foam::Field<double>}’ as ‘this’ argument of ‘void Foam::Field<Type>::operator-=(const Type&)

Just remove const before scalarField.

[Stephen Waite]

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

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

[alexeym]

Hi,

In example case const scalarField is volume of the mesh:

const scalarField& V = mesh_.V();

equation source is non-constant reference:

scalarField& heSource = eqn.source();

Your second problem has almost nothing to do with fvOptions. In this cycle:

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:

forAll(cells, i)
{
  const label cellIdx = cells[i];
  psiS[cellIdx] -= 0.0001;
}

[Stephen Waite]

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