[Benben]

Hello,

I am coding a new boundary condition in OpenFoam 6. This boundary condition is to be used with the "reactingFoam" solver.
I want to acces the list of chemical species defined from this boundary condition.

I can access the mass fraction of a single specie using :

const Foam::fvPatchField<scalar>& W_N2 = this->patch().lookupPatchField<volScalarField, scalar>("N2");

However, i would like to have access to the complete list of species mass fractions. Just like in the reactingFoam solver :

const PtrList<volScalarField>& Y = thermo.composition().Y();

This doesnt work as thermo is only defined in the solver's "createFields.H" file.
Is there anyway to get access to it from the boundary condition ?

[zhangyan]

Maybe you can get thermo by lookupObject?

[Benben]

Hello, thank you for your interest in my problem.
It seems like it worked.

To create a "psiReactionThermo" object (which is the type of thermo), i have to include "psiReactionThermo.H" in my boundary condition header file, and add the lines :

-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude

to the "Make/options" file.

Then in the code of my boundary condition, i add :

const psiReactionThermo& thermo = this->db().lookupObject<psiReactionThermo>("thermophysicalProperties");
const PtrList<volScalarField>& Y = thermo.composition().Y();

It compiles and i yet have to try to use it.

However i have no clue on why i have access to it from the "ObjectRegistry". I am not very familiar with this class. In the declaration of "thermo", nothing is done to add it to a registry. What other variables do we have access to from there ? every variable created ?

[zhangyan]

https://openfoam.top/en/thermodynamicLIB/
You can have a look at this UML, and you'll find that IOdictionary is the topmost base class of psiReactionThermo.
Since IOdictionary is often registried to the mesh, so will psiReactionThermo.

[Benben]

I have another issue :

when calling

thermo.composition()

the program crashes and prints to the console

pure virtual method called

It seems that "psiReactionThermo" doesnt implement "composition()" (see documentation here).
How should i proceed to get the composition object ?

[Benben]

It's probable that i should define it as an "hePsiThermo" (documentation here) as this is what i am asking openfoam to use in my dictionnary in "constant/thermophysicalProperties":

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

I will try it.

[Benben]

"hePsiThermo" takes two templates class, and i can't manage to declare a instance of the class :/

i Tried :

const hePsiThermo<psiReactionThermo, reactingMixture<psiReactionThermo>>& thermo = this->db().lookupObject<hePsiThermo<psiReactionThermo, reactingMixture<psiReactionThermo>>("t");

And it gave me :

/opt/openfoam6/src/thermophysicalModels/reactionThermo/lnInclude/multiComponentMixture.H:61:28: error: cannot declare field 'Foam::multiComponentMixture<Foam::psiReactionThermo>::mixture_' to be of abstract type 'Foam::psiReactionThermo'
         mutable ThermoType mixture_;
                            ^
In file included from evaporatingZinc_RDMP/evaporatingZinc_RDMPFvPatchScalarField.H:81:0,
                 from evaporatingZinc_RDMP/evaporatingZinc_RDMPFvPatchScalarField.C:26:
/opt/openfoam6/src/thermophysicalModels/reactionThermo/lnInclude/psiReactionThermo.H:52:7: note:   because the following virtual functions are pure within 'Foam::psiReactionThermo':
 class psiReactionThermo

I must be doing things wrong, i will try to look for examples in the code of "thermo = lookupObject("

[zhangyan]

Try this?

const PtrList<volScalarField>& Y = 
dynamic_cast<hePsiThermo<psiReactionThermo,SpecieMixture<reactingMixture<gasHThermoPhysics>>>&>
(this->db().lookupObject<psiReactionThermo>("thermo")).composition().Y();

[Benben]

I tried as you proposed :

const psiReactionThermo& hepsi_thermo = this->db().lookupObject<psiReactionThermo>("thermophysicalProperties");
const hePsiThermo<psiReactionThermo,SpecieMixture<reactingMixture<gasHThermoPhysics>>>& thermo 
    = 
    dynamic_cast<const hePsiThermo<psiReactionThermo,SpecieMixture<reactingMixture<gasHThermoPhysics>>>&>(hepsi_thermo);

but it raised me

terminate called after throwing an instance of 'std::bad_cast'
  what():  std::bad_cast

If I understand correctly, "hePsiThermo<psiReactionThermo,SpecieMixture<react ingMixture<gasHThermoPhysics>>>" is not derived from "psiReactionThermo".

How did you know what types should be used for the cast ? In the code I only found the definition of "thermo" in "createFields.H". And it is defined as a psiReactionThermo :

autoPtr<psiReactionThermo> pThermo(psiReactionThermo::New(mesh));
psiReactionThermo& thermo = pThermo();

Is it "psiReactionThermo::New(mesh)" that would actually return an object of a type derived from "psiReactionThermo" ?

[Benben]

I tried by cheating with pointers :

const psiReactionThermo& hepsi_thermo = this->db().lookupObject<psiReactionThermo>("thermophysicalProperties");
    hePsiThermo<psiReactionThermo,SpecieMixture<reactingMixture<gasHThermoPhysics>>> 
* thermo = & hepsi_thermo;

It fails at compile, and displays :

error: cannot convert 
'const Foam::psiReactionThermo*'
to
'Foam::hePsiThermo<Foam::psiReactionThermo, Foam::SpecieMixture<Foam::reactingMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > > > >*'
in initialization
    hePsiThermo<psiReactionThermo,SpecieMixture<reactingMixture<gasHThermoPhysics>>> * thermo = & hepsi_thermo;

[Benben]

I managed to solve my problem.

The right way to do it was to define

#include "psiReactionThermo.H"
const psiReactionThermo& thermo = this->db().lookupObject<psiReactionThermo>("thermophysicalProperties");
const PtrList<volScalarField>& Y = thermo.composition().Y();

and add

    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude

To the "Make/options" file.

However, the first call of the boundary conditions is made when thermo is beeing built (as mass fractions, rho, and T are initialized by thermo). So, when the boundary condition is called for the first time, thermo is not completly initialized yet. So i had to implement a bypass so that the first call to the boundary condition initializes to some garbage, and doesnt call "thermo.composition()".

The boundary condition is still updated before solving the first time step, so no problem for the resolution (i verified this by outputing some text in the console). The only thing is that paraFoam will show some garbage for the boundary condition at time t=0.