[mkhm]

Dear Foamers,



I would like to write an utility where I could extract the member functions of a known class. For instance, I would like to extract the forward and backward reaction rates. I know that there is a class called ReactionProxy :

https://github.com/OpenFOAM/OpenFOAM...eactionProxy.C


where these two variables are virtual public member functions of it. So in createFields, I have to use a pointer to this class. But I have difficulties to figure out how I should write this pointer. I tried many variants such as :

autoPtr<ReactionProxy> pReaction(ReactionProxy::New(mesh));

but it did not work.



In chemkinReader.C, they use this class like :

case irreversible:
180         {
181             reactions_.append
182             (
183                 new IrreversibleReaction
184                 <Reaction, gasHThermoPhysics, ReactionRateType>
185                 (
186                     ReactionProxy<gasHThermoPhysics>
187                     (
188                         speciesTable_,
189                         lhs.shrink(),
190                         rhs.shrink(),
191                         speciesThermo_
192                     ),
193                     rr
194                 )
195             );

I tried to have inspiration from that part of the code, but I had errors.



Could someone explain how we can use public members of a class by using pointers ? I would be very grateful of your replies. I should write some postProcessing tools but as I am lacking C++ knowledge, I have lots of difficulty.



Thanks in advance for our replies.

[mkhm]

No one can help me ?

[wyldckat]

Greetings mkhm,

1. I've moved your thread to the sub-forum for Programming: https://www.cfd-online.com/Forums/op...ng-development - given that your question is mostly regarding programming
2. If you notice that you've posted a new post or thread in the wrong place, please contact the forum moderators, as indicated in the forum rules page, and I quote:
   

   If you have posted your question in the wrong place please contact a Forum Moderator to have it moved.

   (Sorry for nitpicking, but since you posted that no one was helping you, I felt compelled to point out why possibly you were not getting any answers... )
3. It would be useful to have more details on what you want to do exactly, because it's possible that you are trying to get things to work with the wrong approach.
   • What I mean is that it depends if you are trying to do new reaction calculations or if you are trying to load existing results and process them accordingly, because each strategy requires different approaches.
4. Knowing which solver you've used to simulate the case also helps, given that each solver can use a different mechanism for creating transport and reaction structures.
5. And the last question is if this is something that you only want to calculate after the simulation is complete or if it's something you want to do while the simulation is running.

Best regards,
Bruno

[mkhm]

Originally Posted by wyldckat

Greetings mkhm,

1. I've moved your thread to the sub-forum for Programming: https://www.cfd-online.com/Forums/op...ng-development - given that your question is mostly regarding programming
2. If you notice that you've posted a new post or thread in the wrong place, please contact the forum moderators, as indicated in the forum rules page, and I quote:
   
   (Sorry for nitpicking, but since you posted that no one was helping you, I felt compelled to point out why possibly you were not getting any answers... )
3. It would be useful to have more details on what you want to do exactly, because it's possible that you are trying to get things to work with the wrong approach.
   • What I mean is that it depends if you are trying to do new reaction calculations or if you are trying to load existing results and process them accordingly, because each strategy requires different approaches.
4. Knowing which solver you've used to simulate the case also helps, given that each solver can use a different mechanism for creating transport and reaction structures.
5. And the last question is if this is something that you only want to calculate after the simulation is complete or if it's something you want to do while the simulation is running.

Best regards,
Bruno

Hi Bruno

Thanks for your reply. I am so excited that finally I got an answer to my question. In below, my replies to your questions/hints:

1. Thanks for moving the thread. Indeed, my question is more about C++ programming. But as the context was the openFoam, I thought I might get a reply faster by posting it in post processings forum of OpenFoam. I was indeed wrong.

2. Thanks for mentionning that. It would help me to probably have an answer to my other thereads which did not get any attention.

3. I have to develop some post processing tools. As I am running the steady state simulations, I can as well think to write a tool which would be executed at the final step of my simulations. I tried to take the sensitivity analysis case of openFoam 1806 and compile it for OpenFoam 4.x where I have my modified solvers (My thread about that function object: How to create a function object in OpenFoam that runs properly ? ). However, this seems to be very complicated as this function object depends on many other classes not available in OpenFoam 4.x. So, I started to think to simplify the work and go step by step. For instance, writing an post processing tool where only some useful variables available by virtual functions are extracted. For example: there is this ProxyReaction ( https://github.com/OpenFOAM/OpenFOAM...eactionProxy.C ) or even Reaction (https://cpp.openfoam.org/v6/classFoam_1_1Reaction.html) where the reaction rate coefficients kf (forward) and kr (reverse) could be accessed by a pointer. However, I can not find a right way of writing an autoPointer to these classes to extact data. I explain you my strategy:

The constructers of reaction proxy are:

         //- Construct and return a clone
        virtual autoPtr<Reaction<ReactionThermo>> clone() const;

        //- Construct and return a clone with new speciesTable
        virtual autoPtr<Reaction<ReactionThermo>> clone
        (
            const speciesTable& species
        ) const;

I'll write something like:

    autoPtr<ReactionProxy<gasHThermoPhysics > >  pReaction
    (
            ReactionProxy<gasHThermoPhysics>::clone(species_) 

    );
and for species_:

Info<< "Reading Fake dictionary\n" << endl;
   IOdictionary Fake_
    (
        IOobject
        (
            "Fake",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );
 
            wordList speciesList = Fake_.lookup("species");
            speciesTable species_ = speciesList;
OR 

    autoPtr<ReactionProxy<gasHThermoPhysics > >  pReaction
    (
            ReactionProxy<gasHThermoPhysics>::clone() 

    );

In fisrt case, I have the following error:

cannot call member function ‘Foam::autoPtr<Foam::Reaction<ReactionThermo> > Foam::ReactionProxy<ReactionThermo>::clone(const speciesTable&) const [with ReactionThermo = Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> >; Foam::speciesTable = Foam::hashedWordList]’ without object
          ReactionProxy<gasHThermoPhysics>::clone(species_)

For me, species_ is an object to work on it.

In the second case :

cannot call member function ‘Foam::autoPtr<Foam::Reaction<ReactionThermo> > Foam::ReactionProxy<ReactionThermo>::clone() const [with ReactionThermo = Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> >]’ without object
          ReactionProxy<gasHThermoPhysics>::clone()

4. I took chemFoam solver's make folder, its create fields and the Utility_name.C does nothing for the moment.

5. Ideally, I want it to calculate after simulation. However, my purpose now, it is more learning the C++ by writing some simple tools to elaborate them afterwards.


Thanks a lot for your time. Don't hesitate to ask me if there is a need to clarify more some points.

Best regards,
Mary

[wyldckat]

Quick'ish answer:

Originally Posted by mkhm

4. I took chemFoam solver's make folder, its create fields and the Utility_name.C does nothing for the moment.

I think this wasn't the answer I was looking for... what I meant by my original 4th question, is that I need to know which solver you first run the simulation... because chemFoam is mostly an initialization/testing solver, if I remember correctly...


But if you really use chemFoam to solve your case... mmm... so, here is my train of thought:

1. No cloning should be necessary. You should be able to access the data directly, after getting access to the appropriate object...
2. In chemFoam, the appropriate object is "thermo", defined in the file "createFields.H", this line:
   

   psiReactionThermo& thermo = chemistry.thermo();

   It's just a guess, because it has Reaction in the name...
3. So in concept, within the chemFoam structure, you should do your new code inside the file "output.H", because you said you wanted to only do the math after the time step is completed.
4. First example, what you are trying to get access to is this method:

   Reaction<ReactionThermo>::kf()

5. But what we have access to is the object from the class "psiReactionThermo"... so how are the two connected...?
6. OK, so after being a bit lost in the code, I took a quick look into the tutorial cases for this solver and ended up searching for "kf" in the main "src" folder:
   

   grep -r kf $FOAM_SRC/

7. And the relevant file seems to be:
   

   [...]/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

8. "kf" is being used by the method:
   

   Foam::chemistryModel<CompType, ThermoType>::omega

   But it already has direct access to the "Reaction" object we want to access...
9. So from what I can see in this file, we need to access the "reactions_" list, in order to then retrieve the "kf" for each one on that list...
10. So if we look into the file "chemistryModel.H"... there is a method for getting access to them:
    

    inline const PtrList<Reaction<ThermoType>>& reactions() const;

11. So what we need access to is not the object "thermo", it's the "chemistry" one... therefore, in theory, it should be as easy as using something like this:
    

    const PtrList<Reaction<psiThermo>>& myReactionsAccess =  chemistry.reactions();
    forAll(myReactionsAccess, reactioni)
    {
        const Reaction<psiThermo>& reaction = myReactionsAccess[reactioni];
        
        Info << "For reaction index '" << reactioni << "', "
             << "here is the kf I'm looking for: " << reaction.kf()
             << endl:
    }

I have not tested building/compiling this, but this is the rough idea of what you should be aiming for, namely, a fairly simple implementation, instead of messing around with cloning and proxies... OpenFOAM is not coded as if it was Java

[mkhm]

Dear Bruno,



Thanks for your thorough reply.

I am using a modified version of rhoCentralFoam. The solver is changed in order to simulate a reactive flow.



I tried to implement your idea. I have the following error:

‘class Foam::psiChemistryModel’ has no member named ‘reactions’
      const PtrList<Reaction<psiThermo>>& myReactionsAccess =  chemistry.reactions();

psiChemistryModel is inherited from basicChemistryModel and not the chemistryModel. So, in order to use this function:

inline const PtrList<Reaction<ThermoType>>& reactions() const;

we should make a pointer to either the chemistryModel or a class which is inherited from it. But how, I don't know.

In Reaction proxy and TDAC (https://cpp.openfoam.org/v5/classFoa...3d3b916874f5ca) , I realised that they copied somehow parts of the chemistryModel in order to use some parameters like kf. But as my C++ knowledge is limited, I get lost by this whole information.

Do you have any ideas how to solve this problem ?

Best regards

[hyFoam]

Hi Mary,


If you're willing to modify the source code, you could add a public member function to the chemistryModel class in

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.H

to return kf.


For that, you may want to create a list of volScalarField storing kf in that class and update it for each cell and each species when kf is being calculated (Ctrl+f kf in this file).


Your new method would return kf for a given cell and a given species.



In your solver, I assume that you know how to update the species concentrations, typically:

reaction->correct();

With that in mind, you could well have something similar to:

const scalar kf = reaction->get_kf(speciei, celli);

Of course, you can arrange the syntax and the object returned (scalar, volScalarField).



For more details, here's a solver that is also based on rhoCentralFoam and reactingFoam:
https://github.com/vincentcasseau/hy...sible/hy2Foam/
You could check the definition of reaction and see how it's used in eqns/YEqn.H.



Please let us know if that answers your question.



Thanks,
Vince

[mkhm]

Thanks Vince for your reply. My objective was more to learn how to extract data from existing classes of OpenFoam. Anyway, thanks for proposing this alternative

[wyldckat]

Greetings to all,

OK... that's what happens when I don't test things myself... alright, continuing with the thought+search process:

1. We need to understand how "chemistryModel" is used by "psiChemistryModel". To do this, I ran the following commands:
   

   src
   cd thermophysicalModels/
   grep -r chemistryModel
   kate chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModels.C

   The last command was for opening the file in the text editor Kate.
2. Inside that file has blocks like this one:
   

       makeChemistryModel
       (
           chemistryModel,
           psiChemistryModel,
           constGasHThermoPhysics
       );

   Which means that it constructs a combo class from a template. So we need to look for "makeChemistryModel".
3. I ran the command:
   

   grep -r makeChemistryModel

   and got that it's the file "chemistryModel/chemistryModel/makeChemistryModel.H".
4. So the macro is this:
   

   #define makeChemistryModel(SS, Comp, Thermo)                                   \
                                                                                  \
       typedef SS<Comp, Thermo> SS##Comp##Thermo;                                 \
                                                                                  \
       defineTemplateTypeNameAndDebugWithName                                     \
       (                                                                          \
           SS##Comp##Thermo,                                                      \
           (#SS"<"#Comp"," + Thermo::typeName() + ">").c_str(),                   \
           0                                                                      \
       );

   From which what I care about is this:
   

   typedef SS<Comp, Thermo> SS##Comp##Thermo;

   which will construct:
   

   chemistryModel<psiChemistryModel, constGasHThermoPhysics>

   uh... ok, so we either need to access the parent class or check how the inheritance is carried...
5. So looking at "src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H", the "psiChemistryModel" is used as "CompType" in this definition:
   

   template<class CompType, class ThermoType>
   class chemistryModel
   :
       public CompType,
       public ODESystem

6. OK, this is one of those crazy situations where I refer to this as C++ voodoo... because this means that we need to typecast our access:
   

   psiChemistryModel& chemistry

   as the actual combo from "chemistryModel/chemistryModel/makeChemistryModel.H" for the model combo being used in our specific case.
   • Side note: Now I get the reason why you were trying to use "ReactionProxy"... but that does not exist in OpenFOAM 4, but it exist in OpenFOAM 6, but that was after some considerable reworking of the code.
7. So if I look at the tutorial case "combustion/chemFoam/gri", files "constant/thermophysicalProperties" and "constant/chemistryProperties", to look for ideas on which combo we need... OK, this isn't going to be easy... I guess it should be this combo:
   

       makeChemistryModel
       (
           chemistryModel,
           psiChemistryModel,
           gasHThermoPhysics
       );

   because in "constant/thermophysicalProperties" we have "perfectGas" being used:
   

   thermoType
   {
       type            hePsiThermo;
       mixture         reactingMixture;
       transport       sutherland;
       thermo          janaf;
       energy          sensibleEnthalpy;
       equationOfState perfectGas;
       specie          specie;
   }

8. So what we have to do is some dynamic casting magic and using the correct "thermoType":
   

       chemistryModel<psiChemistryModel, gasHThermoPhysics> &myChemistry =
           dynamic_cast<chemistryModel<psiChemistryModel, gasHThermoPhysics>&>(chemistry);
       const PtrList<Reaction<gasHThermoPhysics>>& myReactionsAccess =  myChemistry.reactions();
       forAll(myReactionsAccess, reactioni)
       {
           const Reaction<gasHThermoPhysics>& reaction = myReactionsAccess[reactioni];
   
           Info << "For reaction index '" << reactioni << "', "
               << "here is the kf I'm looking for: " << reaction.kf()
               << endl;
       }

   However, this still gets me some compilation issues:
   

   output.H:20:66: error: no matching function for call to ‘Foam::Reaction<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > >::kf() const’
                << "here is the kf I'm looking for: " << reaction.kf()
   [...]
   [...]/lnInclude/Reaction.C:444:14: note: candidate: Foam::scalar Foam::Reaction<ReactionThermo>::kf(Foam::scalar, Foam::scalar, const scalarField&) const [with ReactionThermo = Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> >; Foam::scalar = double; Foam::scalarField = Foam::Field<double>]
    Foam::scalar Foam::Reaction<ReactionThermo>::kf
                 ^
   [...]/lnInclude/Reaction.C:444:14: note:   candidate expects 3 arguments, 0 provided

   Oh... because it needs 3 arguments...
9. Which by looking for it here: https://cpp.openfoam.org/v4/classFoam_1_1Reaction.html
   

   kf (const scalar p, const scalar T, const scalarField &c) const

   means I'll give it some values just for kicks and re-adapt the code:
   

       chemistryModel<psiChemistryModel, gasHThermoPhysics> &myChemistry =
           dynamic_cast<chemistryModel<psiChemistryModel, gasHThermoPhysics>&>(chemistry);
       const PtrList<Reaction<gasHThermoPhysics>>& myReactionsAccess =  myChemistry.reactions();
       forAll(myReactionsAccess, reactioni)
       {
           const Reaction<gasHThermoPhysics>& reaction = myReactionsAccess[reactioni];
   
           Info << "For reaction index '" << reactioni << "', "
               << "here is the kf I'm looking for: " << reaction.kf(101325, 300, Y[0])
               << endl;
       }

10. Run wmake... and it finally compiled! Now I wonder if it crashes or not...
11. I ran my modified solver in the tutorial case I mentioned before and got this:
    

    Starting time loop
    
    deltaT = 9.9999e-08
    Time = 9.9999e-08
    
    Sh = -1005.62, T = 1000, p = 1.36789e+06, CH4 = 0.0551873
    For reaction index '0', here is the kf I'm looking for: 5198.52
    For reaction index '1', here is the kf I'm looking for: 2e+10
    For reaction index '2', here is the kf I'm looking for: 1.05704e+06
    For reaction index '3', here is the kf I'm looking for: 5.7e+10
    For reaction index '4', here is the kf I'm looking for: 8e+10
    For reaction index '5', here is the kf I'm looking for: 1.5e+10
    For reaction index '6', here is the kf I'm looking for: 1.5e+10
    For reaction index '7', here is the kf I'm looking for: 5.06e+10
    For reaction index '8', here is the kf I'm looking for: 2881.8
    For reaction index '9', here is the kf I'm looking for: 3e+10
    [...]
    For reaction index '283', here is the kf I'm looking for: 3.10601e+09
    For reaction index '284', here is the kf I'm looking for: -nan
    For reaction index '285', here is the kf I'm looking for: -nan
    For reaction index '286', here is the kf I'm looking for: -nan
    For reaction index '287', here is the kf I'm looking for: -nan
    For reaction index '288', here is the kf I'm looking for: -nan
    For reaction index '289', here is the kf I'm looking for: -nan
    For reaction index '290', here is the kf I'm looking for: -nan
    For reaction index '291', here is the kf I'm looking for: -nan
    For reaction index '292', here is the kf I'm looking for: -nan
    For reaction index '293', here is the kf I'm looking for: -nan
    For reaction index '294', here is the kf I'm looking for: -nan
    For reaction index '295', here is the kf I'm looking for: -nan
    For reaction index '296', here is the kf I'm looking for: -nan
    #0  Foam::error::printStack(Foam::Ostream&) at ??:?
    #1  Foam::sigFpe::sigHandler(int) at ??:?
    #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
    #3  Foam::ReversibleReaction<Foam::Reaction, Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> >, Foam::FallOffReactionRate<Foam::ArrheniusReactionRate, Foam::TroeFallOffFunction> >::kf(double, double, Foam::Field<double> const&) const at ??:?
    #4  ? at ??:?
    #5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
    #6  ? at ??:?
    Floating point exception (core dumped)

12. I'm guessing the problem is really because I simply threw in some pressure and temperature values, without caring for the case setup. But it seems to work just fine

Beyond this, yesterday I found another thread discussing accesses to new entries on boundary conditions and dictionaries in general, which might come in handy for you as well for something else: Accessing boundary patch dictionary within thermodynamic library



Best regards,
Bruno

[mkhm]

Thanks a lot Bruno.



It works fine for me. Without any crashes
Special thanks for sharing your thought+search process. It is a very useful teaching way. I learned a lot. And I am sure that I will read again and again this thread for other implementation processes.



Thanks for the time you put for this case.