CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

How force fixed value of variable in one cell

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   August 25, 2006, 12:11
Default Hello, I have an equation f
  #1
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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
olwi is offline   Reply With Quote

Old   August 25, 2006, 12:17
Default One way to do it, though it is
  #2
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 124
Rep Power: 12
mprinkey will become famous soon enough
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.
mprinkey is offline   Reply With Quote

Old   August 25, 2006, 12:30
Default Check e.g. icoFoam.C (around l
  #3
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 16
mattijs is on a distinguished road
Check e.g. icoFoam.C (around line 86):

pEqn.setReference(pRefCell, pRefValue);
mattijs is offline   Reply With Quote

Old   August 25, 2006, 12:46
Default Thanks for the quicke response
  #4
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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
olwi is offline   Reply With Quote

Old   August 25, 2006, 13:07
Default Mattijs's comment is actually
  #5
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
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
pyt likes this.
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   August 25, 2006, 13:10
Default Forgot to say: this will work
  #6
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   August 25, 2006, 13:16
Default "Mattijs's comment is actually
  #7
Senior Member
 
Michael Prinkey
Join Date: Mar 2009
Location: Pittsburgh PA
Posts: 124
Rep Power: 12
mprinkey will become famous soon enough
"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
mprinkey is offline   Reply With Quote

Old   August 25, 2006, 14:46
Default Sorry sorry, really sorry Mike
  #8
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   August 28, 2006, 08:03
Default My goodness, we're not here to
  #9
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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.
olwi is offline   Reply With Quote

Old   August 29, 2006, 02:44
Default Consider 2 cases of incompress
  #10
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   August 29, 2006, 06:59
Default Hmmm. Sounds like a smart solu
  #11
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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 is offline   Reply With Quote

Old   August 29, 2006, 07:15
Default Hmmm. Sounds like a smart solu
  #12
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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 is offline   Reply With Quote

Old   August 29, 2006, 08:24
Default Actually, (I am sorry, but) yo
  #13
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   August 29, 2006, 09:05
Default No need to be embarrased, even
  #14
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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 is offline   Reply With Quote

Old   August 31, 2006, 03:57
Default Hi, Just to close the threa
  #15
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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
olwi is offline   Reply With Quote

Old   November 8, 2006, 03:38
Default Hi, I realize that the thread
  #16
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
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
cliffoi is offline   Reply With Quote

Old   November 8, 2006, 06:27
Default Have you checked http://www.c
  #17
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 16
mattijs is on a distinguished road
Have you checked http://www.cfd-online.com/cgi-bin/Op...=9720#POST9720
mattijs is offline   Reply With Quote

Old   November 8, 2006, 06:39
Default Yes, the answer is in the post
  #18
Member
 
Ola Widlund
Join Date: Mar 2009
Location: Sweden
Posts: 87
Rep Power: 8
olwi is on a distinguished road
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
olwi is offline   Reply With Quote

Old   November 8, 2006, 06:46
Default Thanks Ola. Could you please e
  #19
Member
 
Ivor Clifford
Join Date: Mar 2009
Location: Switzerland
Posts: 91
Rep Power: 8
cliffoi is on a distinguished road
Thanks Ola. Could you please email the bugfix to me. I'd appreciate it.

Regards
Ivor
cliffoi is offline   Reply With Quote

Old   November 8, 2006, 07:55
Default Hi Guys, My OpenFOAM develo
  #20
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,761
Rep Power: 21
hjasak will become famous soon enough
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
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
How force solver to update variable with nearly diagonal matrix olwi OpenFOAM Running, Solving & CFD 0 February 29, 2008 04:52
NS with variable body force-Momentum conserved? CFDtoy Main CFD Forum 0 January 25, 2008 00:47
*VARIABLE CELL SIZE IN LES p.cooper Main CFD Forum 1 July 28, 2002 02:27
How to add a fixed force on momentum equation sun Phoenics 3 January 19, 2002 07:31
Defining a new cell variable Ulrich FLUENT 1 May 25, 2000 16:04


All times are GMT -4. The time now is 18:17.