[megacrout]

Hi Foamers,

I ran into the error message shown in the title while trying to modify porousExplicitSourceReactingParcelFoam. This solver works with icoPoly8ThermoPhysics for the computation of the thermophysical properties but I´d like it to use gasThermoPhysics (Janaf/Sutherland) instead (I use only gaseous species). In other words, I´d like to use the data I used with reactingFoam.

What I did is the following: replace in the .C file of the solver and in its createFields file:
- h(s)ReactionThermo by h(s)CombustionThermo
- rhoChemistryModel by psiChemistryModel

All other functions are already the same. The (more complete) error message I got while compiling the new solver is:

pEqn.H:6: error: no match for ´operator-=´ in ´Foam::basicPsiThermo::rho() -= Foam:perator*(const Foam::GeometricField<double, PatchField, GeoMesh>&, const Foam::GeometricField<double, PatchField, GeoMesh>&) [with PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&)((const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>*)p)))´

pEqn.H:54: error: no match for ´operator+=´ in ´Foam::basicPsiThermo::rho() += Foam:perator*(const Foam::GeometricField<double, PatchField, GeoMesh>&, const Foam::GeometricField<double, PatchField, GeoMesh>&) [with PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&)((const Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>*)p)))´

I read that this kind of error message pops up when members of an equation are not of the same kind (many threads on this). Lines 6 and 54 of pEqn.H correspond to:

thermo.rho() -= psi*p;
and
thermo.rho() += psi*p;

So, I am not sure how why psi and/or p are now of a different type, simply because I modified the way they get computed and, more importantly, I do not know how to solve this problem (I´ve used reactingFoam quite a lot but never modified any solver). Other threads did not help me much.

Can someone help on this?
Thanks in advance.

Tibo

[adhiraj]

I read that this kind of error message pops up when members of an equation are not of the same kind (many threads on this).

I thought if both objects are volScalarFields but with different units, when you try to assign one to the other, you will get an error while running. The compiler does not catch errors in units.
I don't fully understand why these two lines are not working, But it may
have something to do with fact that the function in basicPsiThermo (thermo.rho()) returns a reference to density with the keyword 'tmp'.
Try defining a volScalarField& 'RhoDefined' and assign the thermo.rho() to RhoDefined.

volScalarField& RhoDefined=thermo.rho();

Then use 'RhoDefined' for these two lines.

[megacrout]

Hi Adhiraj,

Thanks for the reply. I first thought your explanation was not right because the pEqn.H file starts with:

rho = thermo.rho();

so I expected an error message to appear for that line too. I gave the solution you suggested a shot anyway and got following message:

pEqn.H:3: error: invalid initialization of non-const reference of type ´Foam::volScalarField&´from a temporary of type ´Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >´

So you were right about the tmp nature of thermo.rho() but apparently not about the way to call it into pEqn.H (were you?). Can you help on this second error message?

It might be straightforward stuff, but it is pretty new to me. Thanks in advance for your help and patience.

[adhiraj]

I see now. The variable 'rho' is the new variable defined in createFields.H, and it gets the value from 'thermo.rho()'. Why don't you just use 'rho' for your computations? If you need a different variable, follow what is being done for 'rho'.

[megacrout]

I really need two parameters with the value of the computed density, as you suggested first. Why? The pEqn.H file looks like this:

 rho = thermo.rho();
thermo.rho() -= psi*p;
[recomputing of p and psi using rho]
thermo.rho() += psi*p;

I think this talks for itself =D But as you pointed out, thermo.rho() being a tmp, it needs to be called and stored in a new variable.

I got one last error message while compiling, quite a strange one actually:

 [: 41: -a: unexpected operator

No more details. Turned out I wrote:

 EXE = $(USER_APPBIN)/blablaFoam

instead of

 EXE = $(FOAM_USER_APPBIN)/blablaFoam

in my new solver´s Make/files file. Fixed that, solver compiled.
I still have to check that I can run my case as desired with that new solver, of course, but so far so good.

Thanks for all your help so far.

Tibo

[megacrout]

Update.

The solver compiled but I could not run any simulation with it. As far as a understand the corresponding error message, my implementation was not compatible with some other functionality of the original solver that I had turned off. More specifically, I got

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

where

Selecting CompositionModel SingleMixtureFraction

used to appear. From this I located the problem somehwere in constant/reactingCloud1Properties but I didn´t get much more into this.
Instead - because I do not use most functionalities of porousBlablaReactingParcelFoam - I modified the solver further and deleted all functions that I do not use. So what I am left with is a very simplified porousBlaReactingParcelFoam based on Janaf/Sutherland which is the same as to say a porous reactingFoam.

Good news is, I don´t get any bad_cast anymore. Problem is, I get:

--> FOAM FATAL ERROR: 
attempted to assign to a const reference to constant object
From function tmp<T>::operator=(const tmp<T>& t)
in file /home/tibo/OpenFOAM/OpenFOAM-1.7.1/src/OpenFOAM/lnInclude/tmpI.H at line 270.
FOAM aborting

in the log file, on the line on which information about the first computation of rho should appear.

This error message can be found in only three threads in this forum and those do not help me much. I suppose the mistake comes from the way I implemented 'rhoDefined' but I am not sure what exactly. Any idea?

Thx.

[adhiraj]

Can you see the code line which causes this error?

[megacrout]

"No". Entire log file:

Create time
Create mesh for time = 0
 
Reading chemistry properties
　
Reading g
 
Reading thermophysical properties
Selecting psiChemistryModel ODEChemistryModel<gasThermoPhysics>
Selecting thermodynamics package hsPsiMixtureThermo<reactingMixture<gasThermoPhysics>>
Selecting chemistryReader foamChemistryReader
Selecting chemistrySolver ode
Selecting ODE solver SIBS
ODEChemistryModel: Number of species = 7 and reactions = 5
 
Reading field U
 
Reading/calculating face flux field phi
Creating turbulence model
Selecting turbulence model type RASModel
Selecting RAS turbulence model laminar
Creating multi-variate interpolation scheme
Creating porous zones
Creating porous zone: filter
 
Courant Number mean: 0.004158804042 max: 0.03048215486
Starting time loop
Courant Number mean: 0.004158804042 max: 0.03048215486
 
deltaT = 1.199904008e-05
Time = 1.1999e-05
 
Solving chemistry
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 1.03157411e-05, Final residual = 3.322125927e-10, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 5.89709041e-05, Final residual = 1.737820003e-09, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.00572078929, Final residual = 2.038681855e-07, No Iterations 1
DILUPBiCG: Solving for CH4, Initial residual = 0.006611153486, Final residual = 9.122130143e-12, No Iterations 2
DILUPBiCG: Solving for O2, Initial residual = 5.812337977e-05, Final residual = 2.667974859e-14, No Iterations 2
DILUPBiCG: Solving for CO, Initial residual = 5.812390347e-05, Final residual = 3.671877178e-09, No Iterations 1
DILUPBiCG: Solving for CO2, Initial residual = 9.529344853e-07, Final residual = 3.084257859e-11, No Iterations 1
DILUPBiCG: Solving for H2, Initial residual = 5.812281717e-05, Final residual = 3.668396256e-09, No Iterations 1
DILUPBiCG: Solving for H2O, Initial residual = 7.498871945e-05, Final residual = 1.330169042e-09, No Iterations 1
DILUPBiCG: Solving for hs, Initial residual = 0.0006554718076, Final residual = 3.97606595e-08, No Iterations 1
T gas min/max = 800, 817.6749767
GAMG: Solving for p, Initial residual = 0.2773997563, Final residual = 4.469223421e-07, No Iterations 15
　
--> FOAM FATAL ERROR: 
attempted to assign to a const reference to constant object
From function tmp<T>::operator=(const tmp<T>& t)
in file /home/tibo/OpenFOAM/OpenFOAM-1.7.1/src/OpenFOAM/lnInclude/tmpI.H at line 270.
FOAM aborting
 
#0 Foam::error::printStack(Foam::Ostream&) in "/home/tibo/OpenFOAM/OpenFOAM-1.7.1/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/tibo/OpenFOAM/OpenFOAM-1.7.1/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 
in "/home/tibo/OpenFOAM/tibo-1.7.1/applications/bin/linuxGccDPOpt/myJanafPorousReactingFoam"
#3 __libc_start_main in "/lib/i686/cmov/libc.so.6"
#4 
in "/home/tibo/OpenFOAM/tibo-1.7.1/applications/bin/linuxGccDPOpt/myJanafPorousReactingFoam"
Aborted

A similar, working reactingFoam case leads to:

Solving chemistry
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 4.75983e-05, Final residual = 8.53939e-11, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 0.000258606, Final residual = 5.58587e-10, No Iterations 1
DILUPBiCG: Solving for Uz, Initial residual = 0.000340047, Final residual = 3.04858e-10, No Iterations 1
DILUPBiCG: Solving for CH4, Initial residual = 2.62297e-05, Final residual = 7.19764e-11, No Iterations 1
DILUPBiCG: Solving for O2, Initial residual = 2.61903e-05, Final residual = 9.2364e-11, No Iterations 1
DILUPBiCG: Solving for CO, Initial residual = 2.62444e-05, Final residual = 4.18228e-11, No Iterations 1
DILUPBiCG: Solving for CO2, Initial residual = 2.0143e-05, Final residual = 3.10614e-11, No Iterations 1
DILUPBiCG: Solving for H2, Initial residual = 2.62166e-05, Final residual = 4.17723e-11, No Iterations 1
DILUPBiCG: Solving for H2O, Initial residual = 7.05618e-05, Final residual = 6.07149e-16, No Iterations 2
DILUPBiCG: Solving for hs, Initial residual = 3.079e-05, Final residual = 1.02736e-10, No Iterations 1
T gas min/max = 800, 921.852
DICPCG: Solving for p, Initial residual = 0.710627, Final residual = 5.66667e-07, No Iterations 9
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 5.55144e-12, global = -1.13361e-12, cumulative = -1.13361e-12
DICPCG: Solving for p, Initial residual = 0.000554997, Final residual = 2.10293e-07, No Iterations 5
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 4.28187e-12, global = -1.1607e-12, cumulative = -2.29431e-12
ExecutionTime = 2.86 s ClockTime = 3 s

so the problem is indeed in pEqn.H (s. e.g. reactingFoam.C or porousExplicitSourceReactingParcelFoam.C), where I implemented the rhoDefined parameter we talked about in previous posts...

[adhiraj]

Originally Posted by megacrout

"No".

I believe if you run in debug mode with gdb you can see which line and which file the error is.
It appears that your method tried to write something to a read only variable. Can you post the pEqn.H file as well as the creatFields.H file where you defined rhoDefined?

[megacrout]

I just read about debug mode, I had never used it before.
I have not compiled OF in debug mode yet, so I won´t have a complete report soon but I´ll give it a shot in top-level debugging mode and post the result later.

I attached createFields and pEqn. You´ll see I did not define rhoDefined in createFields because I only need it in pEqn, it can thus be created in pEqn (can´t it?).

Update:
I compiled the solver in fullDebug mode but I did not get anything special. Should I be looking for an output in a particular file? If not, I´d say it ran fine. Still have to run it in debug-OF though.

[megacrout]

All right. I should not need to recompile OF. I simply hadn´t used the top-level debug mode properly.

It says the error is in pEqn.H, line 68. That is:

thermo.rho() = rhoDefined;

I can see why it bugs there and I could actually tell from the beginning it was no pretty code but I do not know how to write it properly. I have the feeling I have to redefine the class of thermo.rho, i.e. thermo, but I do not know where and how this has to be done. Somewhere in createFields, I guess?

[adhiraj]

In one line of your pEqn.H you are trying to assign a value to 'thermo.rho()'. I am not sure that is allowed.
I dug a bit deep into the thermo files , but I don't have a real answer to some things that are there.

[adhiraj]

OK, you beat me to it.
One suggestion is that do something like this:

rho = thermo.rho();
rhoDefined=rho-psi*p;
<do whatever>
rhoDefined+=psi*rho;

[megacrout]

Originally Posted by adhiraj

assign a value to 'thermo.rho()'. I am not sure that is allowed.

No, it is not, indeed, and THAT is precisely my problem. And that is also why your last suggestion does not work. It is pretty much the same as what you suggested before but it does not solve the dirty step

thermo.rho() = rhoDefined;

Anyway, I´ll have a look at the thermo files. There probably is no better/easier way.

Thx again.

[megacrout]

Hey Adhiraj,

Modify the thermo package(s) was not the right way, as I pretty quickly expected. I found an alternative solution, which is very close to what you suggested.
My main problem was a basic misunderstanding of how the value of the density field gets computed/saved, i.e. of what "thermo.rho()" and "rho" are/do.

So, my pEqn file now looks like this:

volScalarField rhoDefined ((...), rho);
rhoDefined -= psi*p;
(recompute p and psi)
rhoDefined += psi*p;
rho = rhoDefined;
#include "rhoEqn.H"

i.e I removed the first line

rho = thermo.rho();

and included it in the solver.C file which now reads:

rho = thermo.rho();
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
#include "pEqn.H"
}

So, the right value of rho is used for the first iteration of the PISO loop thanks to the modified solver.C file and the right value of rho is used in the following iterations and eventually saved as the density field thanks to the modifications in the pEqn.H file (the value of the density field actually gets recomputed using again rho = thermo.rho(); and saved later in the solver.C file but I did not show it here...)

The solver compiled properly, simulation ran and results looks good to me. Definitely have to have a second look at it though.

Thx a lot for your help!