[olwi]

Hello,

I have an equation for which one will often use only Neumann boundary conditions. If you're lucky an iterative solver will cough up a solution anyway, although at some arbitrary offset level that can be anything.

To deal with this in a reproducible way one would either have to give the variable a fix value at some boundary, or fix it to a given value in one single cell somewhere in the flow.

So: What is the most elegant openFOAM-way to force a variable to given value in one single cell? In a code other than OpenFOAM I would make some violence on the matrix coefficient for the cell in question, or design an implicit source term only in that cell. In OpenFOAM I'm note sure how to do that. But I'm sure it's very elegant...

Thanks i advance,

Ola

[mprinkey]

One way to do it, though it is not at all elegant, is to build a linearized source term just for that cell. Make the Ap contribution some large number, 1.e30 say. Make the Source contribution the desired value times that same large number. The result is that all of the neighbor coefficients end up in round off error and the linear solver will produce the specified value you intended (1.e30*value/1.e30 = value).

I've done this many times in Fluent. There shouldn't be any reason you couldn't do the same thing in OpenFOAM.

[mattijs]

Check e.g. icoFoam.C (around line 86):

pEqn.setReference(pRefCell, pRefValue);

[olwi]

Thanks for the quicke response, both of you.

Michael's idea: Yep, that's what I'd do in Fluent, too... I guess I would then construct a zero volScalarField in which I would set some cell with number "celli", say, as you propose? Works for me. How do I set/adress that single cell value?

Mattjis: I actually tried that, but it doesn't help. I don't understand why, though... Is that really what that line does? Should it be "activated" in some way? I have this in the createFields.H:
...
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

label epotRefCell = 0;
scalar epotRefValue = 0.0;
setRefCell(epot, mesh.solutionDict().subDict("MHD"), epotRefCell, epotRefValue);
...

and this in the main:
...
epotEqn.setReference(epotRefCell, epotRefValue);
epotEqn.solve();
...

It's just what you propose, Mattjis, so there's a mystery for you. If that would work, it clearly qualifies as the most elegant openFOAM-way of doing it!

/Ola

[hjasak]

Mattijs's comment is actually better: messing about with the central coefficient is so 90's :-) It will cause you trouble with residual normalisation and convergence, because too much of the residual is associated with one cell.

The setReference does the equivalent, but "properly" (sorry, Mike): it goes into the matrix, finds the equation, uses the fixed value to multiply it with the neighbouring coefficients to account for the effect in the rest of the domain and then sets off-diagonal coefficients to zero for the row where the solution is fixed.

This is in effect the same as making the diagonal really really big, but does not change the condition number of the matrix, achieves diagonal dominance in neighours in the proper manner and does not require any filtering in the solver or special treatment of Dirichlet points in AMG coarsening.

(after the advert, I think you can guess where it comes from...)

Hrv

[hjasak]

Forgot to say: this will work - check your code.

Also, for vector variables, where we only wish to fix some of the components (e.g. 2-D with symmetry planes), there is a component-wise equivalent, but it works using boundary coefficients.

Search for something like:

UEqn.setComponentReference(1, 0, vector::X, 0);

Hrv

[mprinkey]

"Mattijs's comment is actually better: messing about with the central coefficient is so 90's :-)"

Wow, Hrv. You know how to hurt a guy! 8)

In FLUENT, I did a similar thing as you do with .setReference...properly sizing the Ap and RHS, and zeroing out the neighbor coeffs to get away from scaling issues. You aren't supposed to be able to do that with FLUENT UDFs, but you can.

Best,

Mike

[hjasak]

Sorry sorry, really sorry Mike, I did not mean to come over all rude. Anyway, both you and I were doing CFD "in the 90's" and this cannot be helped. I'll buy you a beer when we next meet to make it up to you. It is nice to know you also did it the "proper way" as well :-)

Have a good weekend,

Hrv

[olwi]

My goodness, we're not here to relive the 90's, are we...

I've checked the code for fvMatrix.C/.H. It seems to me that setValues does what Hrvoje describes, but I have my doubts about setReference... setReference checks the needReference flag. How is that set?

Looking at icoFoam and accompanying example, pressure in cell 0 is set to 0.0 with setReference. But what if I have fixedValue conditions for pressure on an outlet?! Then that would be wrong. How is it disabled/enabled? We only want to fix a value in case we have only Neumann bc:s, making sure it's not floating all over the place.

/Ola


PS: Confession: When I did it in CFX years ago, I gladly did it the ugly way, messing with the central coefficient. It won't happen again. Ever.

[hjasak]

Consider 2 cases of incompressible flow, segregated solver: a duct and a lid-driven cavity. In a duct we specify the pressure at the outlet but a cavity has got walls all the way around so the pressure condition is zero gradient. In the first case, setting the reference for the pressure would be wrong (already fixed by the boundary) but in the second it is necessary.

Now generalise this: each fvPatchField provides a fixesValue() member, where fixedValue says yes and zeroGradient says no. Every new b.c. will provide this value. Thus, psi_.needReference() asks all boundary conditions whether they fix the value or not: if none do, the field needs reference.

Looking at the code in fvMatrix.C, I owe additional explanations. In the case of "weakly" referenced field value, I can be gentle with the matrix: at the expense of "fixed value" (of pressure) flapping about a bit, the diagonal is doubled and correction added to the source. This is very "soft-touch" and does the job well. The reason for this is that a strong way of forcing the fixed pressure value would cause a mass continuity error in the cell where the pressure is fixed.

For a strong way of fixing the value, have a look at the wall functions in the epsilon equation:

/home/hjasak/OpenFOAM/OpenFOAM-1.3/src/turbulenceModels/incompressible/wallFunct ions

in wallDissipationI.H, line 42:

epsEqn().setValues
(
p.faceCells(),
epsilon_.boundaryField()[patchi].patchInternalField()
);

The first lot gives you a bunch of indices (equations) and the second the values to be fixed. The business end lives in:
/home/hjasak/OpenFOAM/OpenFOAM-1.3/src/finiteVolume/fvMatrices/fvMatrix

fvMatrix.C, line 416 onwards.

Enjoy,

Hrv

[olwi]

Hmmm. Sounds like a smart solution. Remains getting my code work...

I tried printing p.needReference() and epot.needReference(), to check their contents. epot is my new variable with only zeroGradient bc:s, while pressure has a fixedValue at the outlet.

Problem is they are both 0 or both 1 depending on the order in which I call the functions!!! If p.needReference() is called first I get 0's, but if I call epot.needReference() first, then I get 1 both for epot and p... Can this be explained? Smells a static variable somewhere, but that couldn't be, could it?! Do I have to update bc:s or something to make it behave? Two field objects shouldn't affect each other, should they?

If this is actually what goes on, then it explains why my code doesn't work; epot gets no reference level, because the needReference-status is checked for p first in the code.

I must add that I run the cygwin 1.3 port, but I doubt that matters here, it has worked well in everything else.

/Ola

PS. I appreciate your thorough explanation of things, it really helps and motivates in exploring the internal workings of OF.

[olwi]

Hmmm. Sounds like a smart solution. Remains getting my code work...

I tried printing p.needReference() and epot.needReference(), to check their contents. epot is my new variable with only zeroGradient bc:s, while pressure has a fixedValue at the outlet.

Problem is they are both 0 or both 1 depending on the order in which I call the functions!!! If p.needReference() is called first I get 0's, but if I call epot.needReference() first, then I get 1 both for epot and p... Can this be explained? Smells a static variable somewhere, but that couldn't be, could it?! Do I have to update bc:s or something to make it behave? Two field objects shouldn't affect each other, should they?

If this is actually what goes on, then it explains why my code doesn't work; epot gets no reference level, because the needReference-status is checked for p first in the code.

I must add that I run the cygwin 1.3 port, but I doubt that matters here, it has worked well in everything else.

/Ola

PS. I appreciate your thorough explanation of things, it really helps and motivates in exploring the internal workings of OF.

[hjasak]

Actually, (I am sorry, but) you are right: this is the most idiotic and embarrasing long-standing bug in OpenFOAM, something I was sure does not exist. I have fixed it now - recompilation and testing is required.

A current implementation produces a really really bad effect: the same answer for all fields of each type, which is precisely what you see.

Sincere apologies + if you want a bug fix, please let me know. Full recompilation required.

Hrv

[olwi]

No need to be embarrased, even the sun has spots... If it's any comfort, I've filed three different new bugs in the stable Fluent release only the past four weeks. I'm probably one of their least favorite customers.

I've been impressed by the impeccable logic and design of OpenFOAM from the outset, and even more so now that I begin to poke around under the hood. It's a work of art, and you're doing a great job!

A bug fix would be great, but there's no immediate hurry. Do you prefer posting it on the forum? Otherwise here's my private email: ola.widlund@yahoo.com.

/Ola

[olwi]

Hi,

Just to close the thread: Hrv's bugfix fixed the problem (thanks again!). I now use the setReference functionality, and it all comes out ok. Thanks also to Michael and Mattijs for your input.

/Ola

[cliffoi]

Hi,
I realize that the thread is closed but have a query of my own. I have a volScalarField with only fixedGradient boundary conditions and would like to fix an arbitrary internal cell to a known value. I have tried the PEqn.setReference method but this seems to have no influence. Calling PEqn.needReference() before and after setReference is called returns true both times.
Could anybody suggest why this is not working?

Regards
Ivor

[mattijs]

Have you checked http://www.cfd-online.com/cgi-bin/Op...=9720#POST9720

[olwi]

Yes, the answer is in the posting of the 29 August, just scroll up a bit... ;)

You need the bug fix from Hrv (I can send it as well), and then compile the whole OpenFOAM. Hrv has a complete release from the 29 August on his public Mac archive on the net, but I'm not sure how official that is.

/Ola

[cliffoi]

Thanks Ola. Could you please email the bugfix to me. I'd appreciate it.

Regards
Ivor

[hjasak]

Hi Guys,

My OpenFOAM development version shows up regularly on

http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/release/

and you are free to use it. I do not intend to provide pre-compiled versions because it's a pain :-) The library is kept in sync with the other release but contains numerous additional solvers and tutorials, additional discretisation modules etc. as well as all the bug fixes I have done and collected.

As for being official, you may have noticed that OpenCFD Ltd. does not accept public code contributions (the sharp-eyed will also noticed that the list of contributors has also disappeared from the README file and having in mind that I started with OpenFOAM in 1993 and written well over 250k lines of code I am not too pleased). I have no intention of giving copyright on my work to someone else, and this situation is likely to continue for a while.

I've noticed there's no library there at the moment - will upload it as soon as it runs the test loop.

Hrv