[danhnam]

Dear Foamers!

I am developing real gas models such as EoS, thermodynamics and transport (including viscosity, thermal conductivity and Diffusion coefficients) in thermophysicalModels lib for Combustion solvers using OF-6.
In diffusion model implementation, I create 2 new functions to return mass diffusivity in hePsiThermo class as following:

template<class BasicPsiThermo, class MixtureType>
Foam::tmp<Foam::volScalarField>
Foam::hePsiThermo<BasicPsiThermo, MixtureType>::Dimix
(   
    label speciei
) const
{   
    Info << "******Check Dimix function is called from hePsiThermo or not*********" << endl;
    Info << "The value of Dimix[i] with ith =  " << speciei << " is = " << Dimix_[speciei] << endl;
    return Dimix_[speciei];
}

template<class BasicPsiThermo, class MixtureType>
Foam::tmp<Foam::scalarField>
Foam::hePsiThermo<BasicPsiThermo, MixtureType>::Dimix
(   
    const label speciei,
    const label patchi
) const
{
    Info << "******Check Dimix(patchi) function is called from hePsiThermo or not*********" << endl;
    Info << "******Check in Dimix(i, patchi), index is = " << speciei << " and patchi is = " << patchi << endl;
    Info << "******Check the number of patchi in Dimix is= " << Dimix_[speciei].boundaryField().size() << endl;

    return Dimix_[speciei].boundaryField()[patchi];
}

In this class, I made a new protected data member, say Dimix_:

        //- store list of mass diffusion coefficients 
        PtrList<volScalarField> Dimix_;

This data member Dimix_ is initialized in constructor as following:

    Dimix_(MixtureType::numberOfSpecies());
    forAll(Dimix_, i)
    {
        Dimix_.set
        (
            i,
            new volScalarField
            (
                IOobject
                (
                    this->phasePropertyName("thermo:Dimix"),
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh,
                dimensionSet(0, 2, -1, 0, 0)
            )
        );
    }

After that this Dimix_ is updated in the same way rho_ or mu_ can be updated in psiThermo class by calling a function in chungTransport class that I created before.

The compilation was fine. But I got Segmentation fault when running reactingFoam solver even though I did not call any new function Dimix(speciei) or Dimix(speciei, patchi) like this:

Starting time loop

Courant Number mean: 2.73659e-05 max: 0.0017574
deltaT = 1.19999e-06
Time = 1.19999e-06

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
******Check Dimix(patchi) function is called from hePsiThermo or not*********
******Check in Dimix(i, patchi), index is = 47721810817392 and patchi is = 21400536
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in "/lib64/libc.so.6"
#3  Foam::UPtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >::operator[](long) const at ??:?
#4  Foam::hePsiThermo<Foam::psiReactionThermo, Foam::SpecieMixture<Foam::mixReactingMixture<Foam::chungTransport<Foam::species::thermo<Foam::realJanafThermo<Foam::soaveRedlichKwong<Foam::mixSpecie> >, Foam::sensibleEnthalpy> > > > >::Dimix(long, long) const at ??:?
#5  Foam::combustionModels::laminar<Foam::psiReactionThermo>::R(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&) const at ??:?
#6  ? at ??:?
#7  __libc_start_main in "/lib64/libc.so.6"
#8  ? at ??:?
Segmentation fault

As far as I know this caused when a program tries to access memory outside its bounds.
But I don't understand why even though I did not call Dimix(speciei, patchi) function while the program seems to be passed by a random value of speciei into this function.

Does anyone know what the reason of this behavior could be?
Thanks in advance for your help.

[t.teschner]

well, let's look at the stack trace (helpful for future debugging as well):

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in "/lib64/libc.so.6"
#3  Foam::UPtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >::operator[](long) const at ??:?
#4  Foam::hePsiThermo<Foam::psiReactionThermo, Foam::SpecieMixture<Foam::mixReactingMixture<Foam::chungTransport<Foam::species::thermo<Foam::realJanafThermo<Foam::soaveRedlichKwong<Foam::mixSpecie> >, Foam::sensibleEnthalpy> > > > >::Dimix(long, long) const at ??:?
#5  Foam::combustionModels::laminar<Foam::psiReactionThermo>::R(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&) const at ??:?
#6  ? at ??:?
#7  __libc_start_main in "/lib64/libc.so.6"
#8  ? at ??:?
Segmentation fault

Whenever you call a function, you put the memory required in a new so called stack frame. That stack frame just gets inserted on top of the current stack. In you case you have 8 active stack frames. We start by looking at the top most (here stack frame 8) and work our way down, in this way we get an idea which function got called in which order which will help us understand the problem. The first meaningful stack frame (i.e. openfoam related, not from the system libraries) is stack frame 5, calling the function R within laminar class (pressumably within the combustionModel namespace). That function calls Dimix, so this function was called, even though you may not have called it (at this point you would use a debugger to step through your program to understand where the error is comming from). But the function was called, different from your claim, which is also supported by the print statement that you inserted into the class.

So then the next thing you want to look at is what is the error. In this case it is a segmentation fault which is pretty much always memory related (from the top of my head I can't come up with any example where that would not be the case but others may come up with examples here where that would not be memory related).

So, now we know where to look and what to look for. You're function that is called looks like this:

template<class BasicPsiThermo, class MixtureType>
Foam::tmp<Foam::scalarField>
Foam::hePsiThermo<BasicPsiThermo, MixtureType>::Dimix
(   
    const label speciei,
    const label patchi
) const
{
    Info << "******Check Dimix(patchi) function is called from hePsiThermo or not*********" << endl;
    Info << "******Check in Dimix(i, patchi), index is = " << speciei << " and patchi is = " << patchi << endl;
    Info << "******Check the number of patchi in Dimix is= " << Dimix_[speciei].boundaryField().size() << endl;

    return Dimix_[speciei].boundaryField()[patchi];
}

For speciei you get 47721810817392 and for patchi you get 21400536. The size of a unsigned long int is 4294967295 and for a signed long int 2147483647, which are both close to these values and so it seems that these values are not passed in correctly into your function. So you may want to go back to the function R which is calling the Dimixclass (we know that from the stack frame) and see how these parameters were passed.

There are two possibilities here:

1.) Within the R function, we already passed garbage into the Dimixfunction (which means that something went wrong there and we need to spend our debugging effort on that function) or

2.) We pass the values correctly into the hePsiThemro class but the compiler is doing a conversion in the back for us, for example from a signed to an unsigned int. At compile time there is no error because we only know at runtime the values so the compiler can not make a check for that and thus will not complain. This is likely the error here.

To illustrate point 2, consider the following example:

#include <iostream>

void print(unsigned short i) {
    std::cout << "i is: " << i << std::endl;
} 

int main() {
    print(27); // OK, prints 27 which is less then 65535
    print(65537); // also OK, but now prints 1, due to stack overflow
    return 0;
}

the print function will take an unsigned short integer, which can store any integer up to 65535. So if we pass 27 to the function, it prints 27. If we pass 65537 to the function, which is 2 numbers higher than what can be stored, it will start again from zero, so 65536 == 0, 65537 == 1, 65538 == 2 and so on (this is called a stack overflow).

In this case you need to make sure that the variables passed to Dimixare the same as in the R function (i.e. you define both here as a lable, are they also defined in R the same way?) Looking at stack frame 4, the values passed in are Dimix(long, long), so my best guess here is that label is not equal to long on your machine.

Quite lengthy, I admit (I might be wrong), well, but hope this is helpful anyways ...

[danhnam]

Hi Tom!

Thank you so much for your helpful explanation.
Well, I have checked the function R as your suggestion. But it might not be the possibilities 2 because I could not see the relation or the argument that could pass the parameter into my new Dimix function.

This is the code of function R in laminar class. Where R is a function returns the fuel consumption rate matrix based on the name of species. hePsiThermo class including new functions Dimix(speciei) and Dimix(speciei, patchi) is the derived class of ReactionThermo in this case.

template<class ReactionThermo>
Foam::tmp<Foam::fvScalarMatrix>
Foam::combustionModels::laminar<ReactionThermo>::R(volScalarField& Y) const
{
    Info << "*****Check R function is called from laminar class or not ****" << endl;//Nam

    tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));

    Info << "*****Check in R function - after make tSu ****" << endl; //Nam
    fvScalarMatrix& Su = tSu.ref();

    Info << "*****Check in R function - after make Su ****" << endl; //Nam

    if (this->active())
    {
    Info << "*****Check in R function - after active ****" << endl; //Nam
    Info << "*****Check in R function - Y.member() = " << Y.member() << endl; //Nam

        const label specieI =
            this->thermo().composition().species()[Y.member()];

            Info << "check specieI in laminar class is = " << specieI << endl; //Nam 

        Su += this->chemistryPtr_->RR(specieI);
    }
    return tSu;
}

I also put some lines to check where the error occurred. And this is the result when running reactingFoam solver.

Starting time loop

Courant Number mean: 2.73659e-05 max: 0.0017574
deltaT = 1.19999e-06
Time = 1.19999e-06

diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
*****Check R function is called from laminar class or not ****
*****Check in R function - after make tSu ****
*****Check in R function - after make Su ****
*****Check in R function - after active ****
*****Check in R function - Y.member() = O2
******Check Dimix(patchi) function is called from hePsiThermo or not*********
******Check in Dimix(i, patchi), index is = 47272830693744 and patchi is = 32541656
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in "/lib64/libc.so.6"
#3  Foam::UPtrList<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >::operator[](long) const at ??:?
#4  Foam::hePsiThermo<Foam::psiReactionThermo, Foam::SpecieMixture<Foam::mixReactingMixture<Foam::chungTransport<Foam::species::thermo<Foam::realJanafThermo<Foam::soaveRedlichKwong<Foam::mixSpecie> >, Foam::sensibleEnthalpy> > > > >::Dimix(long, long) const at ??:?
#5  Foam::combustionModels::laminar<Foam::psiReactionThermo>::R(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&) const at ??:?
#6  ? at ??:?
#7  __libc_start_main in "/lib64/libc.so.6"
#8  ? at ??:?
Segmentation fault

It looks like the program could not perform this line of code:

        const label specieI =
            this->thermo().composition().species()[Y.member()];

I still could not see the relation that led the program call Dimix(i, patchi) with i = 47272830693744 and patchi is = 32541656.
One more interesting thing here is the values of speicei and patchi are passed differ from time to time.

******Check in Dimix(i, patchi), index is = 47562861859184 and patchi is = 44112856

******Check in Dimix(i, patchi), index is = 47948489318768 and patchi is = 19065816

******Check in Dimix(i, patchi), index is = 46943863427440 and patchi is = 26770392

ANW, I think I have to use the debugger which I haven't used before in this case. And pls let's me know if you have any new idea or new suggestion for me.

Once again, Thank you so much for your time.

[t.teschner]

Getting different values for these variables is a clear sign that the memory for them has not been allocated and they are reading whatever variable (garbage) they find at that memory address when they are accidentally called.

This line

const label specieI =
            this->thermo().composition().species()[Y.member()];

Is causing the issue, as you rightly identified. Try the following

Info << "***** Temp 1 *****" << endl;
const auto temp1 = this->thermo();
Info << "***** Temp 2 *****" << endl;
const auto temp2 = temp1.composition();
Info << "***** Temp 3 *****" << endl;
const auto temp3 = temp2.species();
Info << "***** Temp 4 *****" << endl;
const auto temp 4 = temp3[Y.member()];

In this way you know which of these functions is causing the problem. You could also try to remove the const here, i.e.

Info << "***** Temp 1 *****" << endl;
auto temp1 = this->thermo();
Info << "***** Temp 2 *****" << endl;
auto temp2 = temp1.composition();
Info << "***** Temp 3 *****" << endl;
auto temp3 = temp2.species();
Info << "***** Temp 4 *****" << endl;
auto temp 4 = temp3[Y.member()];

Not that I think it would make a difference, but sometimes it is just good to cover all basis, just to make sure this is not causing the issue.

[danhnam]

Hi Tom!

Thank you for your help.
I am going to try to print more information about related classes and check your suggestions. I will let you know the results.

Thanks you so much!

[danhnam]

Dear Tom!

I did try to add some lines of code as your suggestion but it seems to be not in that way. Compiler complained like this (the reason was clearly shown in the message):

PaSR/../laminar/laminar.C: In member function ‘virtual Foam::tmp<Foam::fvMatrix<double> > Foam::combustionModels::laminar<ReactionThermo>::R(Foam::volScalarField&) const’:
PaSR/../laminar/laminar.C:148:21: error: expected initializer before numeric constant
     const auto temp 4 = temp3[Y.member()];
                     ^
PaSR/../laminar/laminar.C: In instantiation of ‘Foam::tmp<Foam::fvMatrix<double> > Foam::combustionModels::laminar<ReactionThermo>::R(Foam::volScalarField&) const [with ReactionThermo = Foam::rhoReactionThermo; Foam::volScalarField = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]’:
PaSR/PaSRs.C:40:1:   required from here
PaSR/../laminar/laminar.C:142:37: error: cannot allocate an object of abstract type ‘Foam::rhoReactionThermo’
     const auto temp1 = this->thermo();
                                     ^

I also tried some other ways to check what is the main reason.
First I commented out 2 new Dimix(speciei) and Dimix(speciei, patchi) functions while remaining declaration and updating part of Dimix_ data member. The error disappeared.

Then I made a new simple function to return just a scalar like this:

            virtual scalar Di(label speciei) const
            {
                  Info << "*******Check Di(i) is called from psiThermo or not*******" << endl;
                  Info << "Di[i] with ith = " << speciei << " is :" << endl;
                  return 1;
            }

The same error happened. Then I commented it out this new simple function, error disappeard.

This would means that whenever I put a any new function (I did tried with different returing Type) the program was crashed. Actually by this line of code:

        const label specieI =
            this->thermo().composition().species()[Y.member()];

And my temporary conclusion is there could be a incompatibility between some classes e.g., psiThermo, hePsiThermo ... and combustion model classes when I make a new function in psiThermo or hePsiThermo class.
It's not by the PtrList itself.

Now I need to track back some function like thermo() in combustionModel classes, composition() in heThermo and species() in MixtureType to see what going on.

I also draw an incomplete diagram to show how solver and classes that I am dealing with interface together for everyone can easier understand my case.
https://ibb.co/dGY3BQx

And pls let's me know if you have any idea about the possiblilities when you look at the above class diagram interface.

And for others, any comment about this problem would be appreciated.

Thank you so much!

[danhnam]

Dear Tom and everyone!

I am very happy to say that I have solved my problem.
The reason is very simple.
It can be seen in the class diagram posted above that combustionModel classes are based on Thermo Classes e.g., psiThermo, psiReactionThermo ...
After I develop these class, I should connect them with combustionModel lib and compile this lib again. But I did not. So That why when I run the solver, there was a conflict in these libraries. It's my careless.

Solution is modifying Make/options file of combustionModel as following to make sure it is updated based on modified Thermo classes:

EXE_INC = \
    -I$(WM_PROJECT_USER_DIR)/src/transportModels/compressible/lnInclude \
    -I$(WM_PROJECT_USER_DIR)/src/thermophysicalModels/basic/lnInclude \
    -I$(WM_PROJECT_USER_DIR)/src/thermophysicalModels/specie/lnInclude \
    -I$(WM_PROJECT_USER_DIR)/src/thermophysicalModels/reactionThermo/lnInclude \
    -I$(WM_PROJECT_USER_DIR)/src/thermophysicalModels/chemistryModel/lnInclude \
    -I$(WM_PROJECT_USER_DIR)/src/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(WM_PROJECT_USER_DIR)/src/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

LIB_LIBS = \
    -L$(FOAM_USER_LIBBIN) \
    -lcompressibleTransportModels \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -lchemistryModel \
    -lfiniteVolume \
    -lmeshTools

Now my new thermophysical library for real gas, which including SRK for EoS, Chung (1998) for visc and thermal conductivity and Fick's Law for mass diffusivity using binary diffusion coefficients based on Fuller and Takahashi, work well.
This library would be published soon.

Thank you all for your time.

[t.teschner]

Good to hear you managed to fix it. This sound like an annoying error (typically you would get an "undefined reference" error at runtime if you don't compile against the required library). But that explains why you seemingly never called the function but still got a segmentation fault. Thanks for sharing the solution, might be of value for someone else in the future.