Hello, I have an equation f
 

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

One way to do it, though it is
 

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.

Check e.g. icoFoam.C (around l
 

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

pEqn.setReference(pRefCell, pRefValue);

Thanks for the quicke response
 

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! http://www.cfd-online.com/OpenFOAM_D...part/happy.gif

/Ola

Mattijs's comment is actually
 

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

Forgot to say: this will work
 

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

"Mattijs's comment is actually
 

"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

Sorry sorry, really sorry Mike
 

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

My goodness, we're not here to
 

My goodness, we're not here to relive the 90's, are we... http://www.cfd-online.com/OpenFOAM_D...part/happy.gif

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.

Consider 2 cases of incompress
 

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

Hmmm. Sounds like a smart solu
 

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.

Hmmm. Sounds like a smart solu
 

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.

Actually, (I am sorry, but) yo
 

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

No need to be embarrased, even
 

No need to be embarrased, even the sun has spots... http://www.cfd-online.com/OpenFOAM_D...part/happy.gif 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

Hi, Just to close the threa
 

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

Hi, I realize that the thread
 

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

Have you checked http://www.c
 

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

Yes, the answer is in the post
 

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

Thanks Ola. Could you please e
 

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

Regards
Ivor

Hi Guys, My OpenFOAM develo
 

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

Dear Hvorje! It seems for m
 

Dear Hvorje!

It seems for me that your "soft-touch" value setting algorithm has problems with finite, non-zero gradients on boundaries. At least i got a problem with PISO pressure correction algorithm applied to incompressible NS equation (the same as icoFoam) but with volumetric forces. I started with an incompressible buoyancy force, so my "corrected" boundary condition for pressure is snGrad(p)=nf*(0,0,-alpha*g*(T-T0)) on the whole border of domain. Things went well until walls was heated in way which introduces non-zero pressure normal gradient on boundary. Perhaps, I can fix single point (p) value in a "hard" way, however than I can expect (as I got from your post above) the problem with pEqn.flux() algorithm. From the other hand the construction

phi-=fvc::interpolate(rUA*fvc::grad(p))&Mesh.Sf()

instead of

phi-=pEqn.flux()


at the end of PISO loop seems to do some work on well-posed structured mesh (results looks fine for "cubical" domain), however I am afraid this may be not so good on arbitrary unstructured grid. I don't know, it's just my guess about la Place operator which has different numeric discretization as interpolate(grad()), so that both constructions may not be the same.

My questions are

1. Do really setReference() function fails in the cases with non-zero gradients (but still only gradients) on boundary?

2. If yes, how difficult is to correct mass flux (what I mean is to restore/repair a matrix after solve(), so that .flux() method work again) in a single cell where field value is set in a "hard" way? May you (or someone else) can post (or point to location) where I can get information about how exactly flux method works, so I can fix this myself?

thanks a lot for any help!

regards,

/Normunds

Dear Hrvoje and all! I am r
 

Dear Hrvoje and all!

I am really sorry about wrong spelling of your name as well about the content of my previous message - it seems I was wrong, the function setReference() do it's job well, it was not a cause of my problem.

I just recognized (to be more precise it's still my guess) that by adding volume force to momentum conservation equation (as it is suggested somewhere on this board)



fvVectorMatrix UEqn (
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
== volForce
);

solve(UEqn == - fvc::grad(p))


one may end up with mathematically ill defined problem for p due to differences in numeric implementation of operators. In other words integral of normal p derivative (p flux) over all boundary surfaces of solution domain must be exactly consumed by the volumetric source term in the pressure equation or pressure field equation has no solution. The latest seems to be a case when volForce is added to momentum equation in way shown above.

It looks like I have to do

fvVectorMatrix UEqn (
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == - fvc::grad(p)+volForce)

and below in PISO loop

U = rUA*UEqn.H();

phi = (fvc::interpolate(U) & mesh.Sf())
+ (fvc::ddtPhiCorr(rUA, U, phi));

adjustPhi(phi, U, p); // I am not sure about this !!!!

U += rUA*volForce;
phi += (fvc::interpolate(rUA*volForce) & mesh.Sf());


or something like that so that non-zero p integral flux from boundary is exactly compensated by volForce yielded terms in div(phi). However, this is not anymore related to topic of this thread, my apologies to Ola.

Once again my apologies to all!

/Normunds

would you please post the file
 

would you please post the file that should be modified (copied from the above URL) for the above bug to fixed. I use OF V1.3. Thanks.

Best regards,
Maka.

Hi! im new in OpenFoam and tr
 

Hi!
im new in OpenFoam and trying to create a solver for a simple problem like: "div(-k grad (p)=f" with k=id (k might be heterogeneous later). First, how to discretize the source term f when its constant(f=1)? Second, I dont know how to simply impose Neumann boundary conditions. Can somebody help me please!!!
Thx!
stephane.

Hi all. I want to fix the val
 

Hi all.
I want to fix the value of several cells using setValues while equation being solved.
So I think that I have to select cells with cellSet,
and then use setValues with the equation matrix in solver.

I read sets in solver file after cellSet as follows and it seems working fine.

const cellSet selectedCells
(
IOobject
(
"selectedCells",
"constant/polyMesh/sets",
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);

But, setValues function requires not cellSet but labelList....

Thus, my quastions are
1. How can I make labelList from cellSet?
(I think both are similar very much. Am I wrong?)
2. Does anybody think of better way? or Did anybody do the same thing?

regards,
Naoki

Hi all I want to fix the valu
 

Hi all
I want to fix the value of several cells using setValues while equation being solved.
So I think that I have to select cells with cellSet,
and then use setValues with the equation matrix in solver.

I read sets in solver file after cellSet as follows and it seems working fine.

const cellSet selectedCells
(
IOobject
(
"selectedCells",
"constant/polyMesh/sets",
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);

But, setValues function requires not cellSet but labelList....

Thus, my quastions are
1. How can I make labelList from cellSet?
(I think both are similar very much. Am I wrong?)
2. Does anybody think of better way? or Did anybody do the same thing?

regards,
Naoki

I 've solved the problem. I'm
 

I 've solved the problem.
I'm sorry for the silly question.
Here is a part of the code to change cellSet to labelList and do setValue with labelList.

---
cellSet newSet(mesh, setName);
labelList refCells;
refCells = newSet.toc();
scalarField refField(refCells.size());
...
...
someEqn().setValue(refCells,refField);
---
Actually setValue function helps me greatly.

regards,
Naoki

Dear friends,

I am facing the problem you solved in this thread. I ask you a question: maybe you can help.
I want to impose a fixed value of a temperature in a certain volume of my domain.
How can I do this? Should I use the fvMatrix::setValues() method?

Could you briefly explain and what I should practically do?

Thanks a lot,

Samuele