[mturcios777]

I was looking at the thread discussing alternate chemistry models and Cantera Implementation:

http://www.cfd-online.com/Forums/ope...1-cantera.html

There was a message dated October last year about making alternateChemistryModel work with OF >= 1.6. I was wondering two things:

1) Has there been any further work with having alternateChemistryModel interface with later versions of OpenFOAM?

2) I currently have a solver where I call an external utility to calculate reaction rate arrays for use in YEqn.H where the transport equation of select species are solved. I was thinking of doing something more general, but it seems with all the difficulty it may be easier to simply create a second member function like chemistry::calculate2 that takes the names of the select species and calls the external program. Its a pain to have multiple copies of reaction rates and species, so this would simplify and make things consistent.

3) On a more general note, I'm looking for more documentation on the way chemistryModel is templatized. The biggest question I have is why there are both an ODEChemistryModel and an ode solver. In particular, if I wanted to implement something similar to what is described in 2, should I do this in ODEChemistryModel or ode?

Thanks in advance Foam Friends!

[marupio]

I'll tackle 3)

ode should really be named conventionalOdeChemistrySolver.

There are several chemistry models, only one of which uses ODE: ODEChemistryModel. I'm sure you figured that out.

But, the ODEChemistryModel can be solved in several ways: even though it has an ODE set of equations, it doesn't mean it has to use a conventional ODE solver (e.g. RK, SIBS, etc..). So there are several solvers for ODEChemistryModel. These are templated on chemistrySolver: EulerImplicit, sequential, and ode. ode is the one that uses the conventional (and well established) ode solver algorithms.

[mturcios777]

What you say now makes total sense, it would be nice to implement those changes.

This is going to be a long entry due to the code. Sorry in advance

I'm currently using ODEChemistry as a base (as I still need some of the functionality provided by ODEChemistry) by adding some private members and public access function. The public access is for both external access without changes and to write the fields in the time directories. For example, I usually need rho written out for post-processing, so in createFields I use:

    volScalarField rho
    (
     IOobject
     (
      "rho",
      runTime.timeName(),
      mesh,
      IOobject::NO_READ,
      IOobject::AUTO_WRITE
      ),
     thermo.rho()
     );

So, in ODEChemistryModel I have protected entities such as:

//ODEChemistryModel.H
volScalarField dummy_;
...
//ODEChemistryModel.C
dummy_
(
                 IOobject
             (
                 "dummy",
                 mesh.time().timeName(),
                 mesh,
                 IOobject::MUST_READ,
                 IOobject::NO_WRITE
                 ),
             mesh,
             dimless
         )
)

I also have publc const access functions (dummy()), and the code seems to compile correctly. In my solver I am using psiChemistry, and I'm trying to gain access to that function to do

    volScalarField dummy
    (
     IOobject
     (
      "rho",
      runTime.timeName(),
      mesh,
      IOobject::NO_READ,
      IOobject::AUTO_WRITE
      ),
     chemistry.dummy()
     );

When I first tried this I was told that there psiChemistry didn't have a dummy() member, so I tried adding it as a virtual function to psiChemistry to no avail. I'm still working out the inheritance (its been a while since my last OO programming course). Any hints on what I;m doing wrong or if there is a better way of implementing this?

[marupio]

"Any hints on what I;m doing wrong or if there is a better way of implementing this?"

I'm not sure what you are trying to do in the big picture, so I can't say if there's a better way of implementing it. I look at changing the core code as a last resort, but I have had to do it from time to time. Before I do, I look to see if I can pull the objects out of OpenFOAM and reimplement them in my own library with my custom features. Sometimes it's too big an endeavour.

What are you trying to do? These dummy fields look a little awkward - why do you think you need them instead of a reference? You realize the dummy fields are copies, and will not track any changes to the original fields.

If you are committed to modifying the core by adding these dummy fields and the cascade of virtual access functions, I think it is failing because in psiChemistry you need to define a default implementation so the other objects to inherit it don't fail. Use the notImplemented idiom.

[mturcios777]

Originally Posted by marupio

"Any hints on what I;m doing wrong or if there is a better way of implementing this?"

I'm not sure what you are trying to do in the big picture, so I can't say if there's a better way of implementing it. I look at changing the core code as a last resort, but I have had to do it from time to time. Before I do, I look to see if I can pull the objects out of OpenFOAM and reimplement them in my own library with my custom features. Sometimes it's too big an endeavour.

What are you trying to do? These dummy fields look a little awkward - why do you think you need them instead of a reference? You realize the dummy fields are copies, and will not track any changes to the original fields.

If you are committed to modifying the core by adding these dummy fields and the cascade of virtual access functions, I think it is failing because in psiChemistry you need to define a default implementation so the other objects to inherit it don't fail. Use the notImplemented idiom.

I see what you mean. I've made a copy of the ODEChemistry library for modification, and have changed the name for compilation.

I was trying to follow what is done in other classes. For example, in hsCombustionThermo, there is a protected member called hs_ that is created in much the same way as dummy_, and two access functions (one constant, one non-constant) that just return references to hs_. In createFields.H, there is a

volScalarField& hs = thermo.hs();

Which is the non-constant access that is used in transport equations. I'll be updating the dummy_ fields inside the chemistry library, but may need access to them externally for debugging and post-processing. So I'll need to do

dummy = chemistry.dummy();

to get that version updated (something similar is done with rho).

The extra fields I need are related to my external chemistry utility. I actually only really need the chemical source terms for select species to be taken from my external utility, but the way ODEChemistry works as-is doesn't allow that (I don't think).

Thanks for the pointers, I'll keep working on this and post a solution in case anyone needs to do something like this.

[mturcios777]

Alright, I'm still getting the hang of the inheritance (I think I need to have pure virtual versions of the access functions in basicChemistryModel.H), but I have a more pressing question about constructors. In ODEChemistryModel.H, the source term fields are created by the following:

RR_(nSpecie_)
 {
     // create the fields for the chemistry sources
     forAll(RR_, fieldI)
     {
         RR_.set
         (
             fieldI,
             new scalarField(mesh.nCells(), 0.0)
         );
     }
 
     Info<< "ODEChemistryModel: Number of species = " << nSpecie_
          << " and reactions = " << nReaction_ << endl;
 }

where RR_ is a PtrList<scalarField>. I want to do something similar with my new entity where it is a volumeScalarField, something like

RR_(nSpecie_)
 {
     // create the fields for the chemistry sources
     forAll(RR_, fieldI)
     {
         RR_.set
         (
             fieldI,
             new scalarField(mesh.nCells(), 0.0)
         );
     }
 
     Info<< "ODEChemistryModel: Number of species = " << nSpecie_
          << " and reactions = " << nReaction_ << endl;
 }

dummy_(nDummy_) //Error here!
 {
     // create the fields
     forAll(dummy_, fieldI)
     {
         dummy_.set
         (
             fieldI,
             new volScalarField(
                                            IOobject
                                            (
                                                "dummy"+fieldI,
                                                 mesh.time().timeName(),
                                                  mesh,
                                                  IOobject::NO_READ,
                                                  IOobject::NO_WRITE
                                             ),
                                             mesh,
                                             dimless
                                             )

         );
     }
}

I can't compile this, the error I get is:

error: expected constructor, destructor, or type conversion before ‘(’ token

I'm not very sure of the syntax at all, but its seems like it should work from what I read of the PtList::set should be able to take a correctly constructed volScalarField. I think I'm missing something between the two constructors, but I don't know what I can be. Still reviewing my comp-sci, so apologies if the solution is very basic.

[marupio]

You are taking the constructor out of context. The constructor you are using is for the ODEChemistryModel class, not for RR_.

RR_(nSpecie_) // Constructor initialization list for ODEChemistryModel
 { // These braces are the constructor body for ODEChemistryModel
     // create the fields for the chemistry sources
     forAll(RR_, fieldI)
     {
         RR_.set
         (
             fieldI,
             new scalarField(mesh.nCells(), 0.0)
         );
     }
}

Maybe read up a little on constructor initialization lists, and constructor bodies.

[mturcios777]

Not having done intense class implementation since second year undergrad I am truly sorry for my ignorance. Lots of Google and review later, I have been able add the members and have been making loads of modifications on top of a copy of ODEChemistryModel. I haven't changed the souce file names yet or the arguments of the Makefile.

Now I want to make a separate chemistryModel (lets call it myChemistryModel) such that I can call it through the constant/chemistryProperties file:

psiChemistryModel myChemistryModel<gasThermoPhysics>

or

rhoChemistryModel myChemistryModel<gasThermoPhysics>

For consistency’s sake I'll be renaming the folder myChemistryModel, as well as the source files and the class to myChemistryModel. I believe I also need to modify the psiChemistryModels.H and rhoChemistryModels.H to include the new model. What I am really unsure of is what modifications (if any) need to be done on the chemistrySolver files. In particular it looks like all the solvers are built onto ODEChemistryType templates, and I suspect I need to duplicate this code and adapt it to myChemistryModel, but I'm not sure how. This is really big I know, but any pointers at all (even literature) would be beneficial.