# How force fixed value of variable in one cell

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

 August 25, 2006, 12:11 Hello, I have an equation f #1 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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

 August 25, 2006, 12:17 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 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.

 August 25, 2006, 12:30 Check e.g. icoFoam.C (around l #3 Super Moderator   Mattijs Janssens Join Date: Mar 2009 Posts: 1,416 Rep Power: 16 Check e.g. icoFoam.C (around line 86): pEqn.setReference(pRefCell, pRefValue);

 August 25, 2006, 12:46 Thanks for the quicke response #4 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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

 August 25, 2006, 13:07 Mattijs's comment is actually #5 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,761 Rep Power: 21 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

 August 25, 2006, 13:10 Forgot to say: this will work #6 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,761 Rep Power: 21 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

 August 25, 2006, 13:16 "Mattijs's comment is actually #7 Senior Member   Michael Prinkey Join Date: Mar 2009 Location: Pittsburgh PA Posts: 124 Rep Power: 12 "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

 August 25, 2006, 14:46 Sorry sorry, really sorry Mike #8 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,761 Rep Power: 21 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

 August 28, 2006, 08:03 My goodness, we're not here to #9 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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.

 August 29, 2006, 02:44 Consider 2 cases of incompress #10 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,761 Rep Power: 21 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

 August 29, 2006, 06:59 Hmmm. Sounds like a smart solu #11 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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.

 August 29, 2006, 07:15 Hmmm. Sounds like a smart solu #12 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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.

 August 29, 2006, 08:24 Actually, (I am sorry, but) yo #13 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,761 Rep Power: 21 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

 August 29, 2006, 09:05 No need to be embarrased, even #14 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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

 August 31, 2006, 03:57 Hi, Just to close the threa #15 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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

 November 8, 2006, 03:38 Hi, I realize that the thread #16 Member   Ivor Clifford Join Date: Mar 2009 Location: Switzerland Posts: 91 Rep Power: 8 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

 November 8, 2006, 06:27 Have you checked http://www.c #17 Super Moderator   Mattijs Janssens Join Date: Mar 2009 Posts: 1,416 Rep Power: 16 Have you checked http://www.cfd-online.com/cgi-bin/Op...=9720#POST9720

 November 8, 2006, 06:39 Yes, the answer is in the post #18 Member   Ola Widlund Join Date: Mar 2009 Location: Sweden Posts: 87 Rep Power: 8 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

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

 November 8, 2006, 07:55 Hi Guys, My OpenFOAM develo #20 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,761 Rep Power: 21 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

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post olwi OpenFOAM Running, Solving & CFD 0 February 29, 2008 04:52 CFDtoy Main CFD Forum 0 January 25, 2008 00:47 p.cooper Main CFD Forum 1 July 28, 2002 02:27 sun Phoenics 3 January 19, 2002 07:31 Ulrich FLUENT 1 May 25, 2000 16:04

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