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

Writing own boundary condition

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

Reply
 
LinkBack Thread Tools Display Modes
Old   March 11, 2008, 04:10
Default Hi, I'm quite new to openFoam
  #1
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
Hi,
I'm quite new to openFoam and after some tests with the program I'm trying to write a boundary condition on my own.
I would like to implement the analytical solution for the velocity profile of a square duct as a boundary condition.
The program therefore should read in the coordinates of the points belonging to the patch where the boundary condition is applied to, calculate the velocity and then store the result.

So far I've searched the forum and figured out the example given in the wiki for the parabolic velocity profile.

My problem is, that I don't know how to read in the points and store the result for the velocity later on.

Thanks in advance

Christian
wese is offline   Reply With Quote

Old   March 11, 2008, 04:25
Default Look how parabolic velocity do
  #2
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,758
Rep Power: 21
hjasak will become famous soon enough
Look how parabolic velocity does it:

const vectorField& c = patch().Cf();


Now in variable c you have a vector of each face centre for all faces in your patch.

You can now do a calculation of the field; mine was something like

normal vector * peak value * (1 - sqr(local y)

You will change that based on the analytical solution that you need and maybe add different control parameters to allow you to scale/move the solution.

You then call operator= (equals) to assign the results to the patch values:

vectorField::operator=(n_*maxValue_*(1.0 - sqr(coord)));

Enjoy,

Hrv
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   March 11, 2008, 05:01
Default Hi Hrvoje, thanks a lot for y
  #3
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
Hi Hrvoje,
thanks a lot for your fast reply.
Can I than access the coordinates of the points by
{
c.x()
c.y()
c.z()
}
or how do I?
Is the operator applied to all 3 components of the vector?
If so, how do I just apply it on one component?

Thanks in advance
Christian
wese is offline   Reply With Quote

Old   March 11, 2008, 05:22
Default vectorField will give you a li
  #4
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,758
Rep Power: 21
hjasak will become famous soon enough
vectorField will give you a list of vectors for each face centre, each element of which contains the (x y z) coordinates. If you want to get eg. a y-component of fifth face centre, you can do:

c[4].y()

If you want a scalar field of the y-component for all face centres as a field, you do:

scalarField cy = c.component(vector::Y);

Clear?

Hrv
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   March 13, 2008, 03:56
Default Hi Hrvoje, sorry for my late
  #5
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
Hi Hrvoje,
sorry for my late reply, but I had to do some administrative stuff .

With your help I think, that I've cut the Gordian knot. To be honest, it has made me a bit weired, when I first had to think in vectors and tensors when writing source code, but now I was able to compile the project without any errors.

I've applied my new boundary condition to one of my inlets, but when I've tried to solve the problem with interFoam I receive this:

--> FOAM FATAL ERROR :
gradientInternalCoeffs cannot be called for a defaultFvPatchField (actual type squareVelocityProfile)
on patch Einlass_links of field U in file "/home/wagner/cwe/OpenFOAM/cwe-1.4.1/run/tutorials/tmischer_dec_neu7/0/U"
You are probably trying to solve for a field with a default boundary condition.

From function defaultFvPatchField<type>::gradientInternalCoeffs( ) const
in file fields/fvPatchFields/basic/default/defaultFvPatchField.C at line 694.
I've already searched my and your source code, but as far as I can see it looks nearly identical (despite of different names for the classes).

Do you have any idea, why I receive this error message?

Thanks in advance
Christian
wese is offline   Reply With Quote

Old   March 13, 2008, 06:00
Default Yes: your executable does not
  #6
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,758
Rep Power: 21
hjasak will become famous soon enough
Yes: your executable does not see the new boundary condition you implemented. Options:

- you did not instantiate the template, if it is a template
- you did not link the object/library you produced with the executable you are using
- you mis-spelled the name

Hope this helps,

Hrv
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   March 13, 2008, 12:04
Default Hi, you're right, stupid as I
  #7
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
Hi,
you're right, stupid as I am i've fotgotten to alter the ControlDict.
Shame on me ;).

Now it runs, but I receive new error messages:

#0 Foam::error::printStack(Foam:stream&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/lib/linux64GccDPOpt/libOpenFOAM.so"
#2 ?? in "/lib64/libc.so.6"
#3 Foam::squareVelocityProfileFvPatchVectorField::upd ateCoeffs() in "/home/wagner/cwe/OpenFOAM/cwe-1.4.1/lib/linux64GccDPOpt/squareVelocityProfile.so"
#4 Foam::squareVelocityProfileFvPatchVectorField::squ areVelocityProfileFvPatchVectorField(Foam::fvPatch const&, Foam::DimensionedField<foam::vector<double>, Foam::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/cwe-1.4.1/lib/linux64GccDPOpt/squareVelocityProfile.so"
#5 Foam::fvPatchField<foam::vector<double> >::adddictionaryConstructorToTable<foam::squarevel ocityprofilefvpatchvectorfield>::New(Foam::fvPatch const&, Foam::DimensionedField<foam::vector<double>, Foam::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/cwe-1.4.1/lib/linux64GccDPOpt/squareVelocityProfile.so"
#6 Foam::fvPatchField<foam::vector<double> >::New(Foam::fvPatch const&, Foam::DimensionedField<foam::vector<double>, Foam::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam"
#7 Foam::GeometricField<foam::vector<double>, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::GeometricB oundaryField(Foam::fvBoundaryMesh const&, Foam::DimensionedField<foam::vector<double>, Foam::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam"
#8 Foam::GeometricField<foam::vector<double>, Foam::fvPatchField, Foam::volMesh>::readField(Foam::Istream&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam"
#9 Foam::GeometricField<foam::vector<double>, Foam::fvPatchField, Foam::volMesh>::GeometricField(Foam::IOobject const&, Foam::fvMesh const&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam"
#10 main in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam"
#11 __libc_start_main in "/lib64/libc.so.6"
#12 Foam::regIOobject::readIfModified() in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam"
Floating point exception

I think, that #12 is a problem with my code (I've already tried to make my calculation easier and #12 disappeared).
I guess, that I've mad some wrong calculations, but what does the rest tell me?

To make it a bit easier, I've added my source code. Maybe I've made some other strange mistake.

The only input parameter the boundary condition expects is a scalar for "maxVelocity".

squareVelocityProfile.tar.gz

Thanks a lot for all of your help.

Christian
wese is offline   Reply With Quote

Old   March 14, 2008, 06:21
Default I've looked through the source
  #8
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
I've looked through the source code again and found several mistakes, which were responsible for the errors.

So far everything works, but I still have the problem, that the operator doesn't work, when I first create a scalar field, make some calculations with it (including the magnitude of some VectorFields) and than want to apply this in the operator (multiplied with a Vector).

When I just use the Vector with a scalar constant everything works fine, but with the scalarField not.

Must the scalarField be multiplied with a VectorField of same rang to get the program to work?

Thanks a lot in advance
Christian
wese is offline   Reply With Quote

Old   March 14, 2008, 09:01
Default Finally nearly everything work
  #9
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
Finally nearly everything works, but not I'm facing a new problem while writing this boundary condition.
See:
http://www.cfd-online.com/OpenFOAM_D...tml?1205499005
wese is offline   Reply With Quote

Old   April 4, 2008, 05:22
Default To complete this post, I would
  #10
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
To complete this post, I would like to tell, that I've finally solved the problem.
The new bc now also runs in parallel mode (the only requirement is, that the patch isn't split between different processors):
http://www.cfd-online.com/cgi-bin/OpenFOAM_Discus/show.cgi?1/7146
wese is offline   Reply With Quote

Old   April 4, 2008, 06:13
Default Hi Christian I have just be
  #11
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Deltares, Delft, The Netherlands
Posts: 1,603
Rep Power: 25
ngj will become famous soon enoughngj will become famous soon enough
Hi Christian

I have just been looking at the source code ... I see that you are using bounding box to find maximum and minimum of your patch.
I have had problems with bounding box in parallel as well, so you could try to specify your a and b vectors manually and out-comment the bounding box part, then from my experience you should be able to run in parallel without any limitations on where the splitting is made.

Have fun,

Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   April 4, 2008, 07:36
Default Hi Niels, thanks a lot for yo
  #12
New Member
 
Christian Wesemeyer
Join Date: Mar 2009
Posts: 19
Rep Power: 8
wese is on a distinguished road
Hi Niels,
thanks a lot for your reply.
My main problem with the boundingBox ind parallel was, that I received a lot of errors, if the box was of zero size.
The problem with the splitting through the patch is, that he takes the with and height of the patch as basis for the calculation.
I've tried this with simpleFoam and a square duct, where I've put the decomposition line just through my patch.
When I watched the output result I've seen the profile on the inlet but instead of a single one for the complete patch I've seen two for the "two halves" of the patch (one for every processor region).

I'll keep your idea in mind, maybe one day I'll need to implement it, but then I'll have to define the parameters for the boundary much more concrete to the specific mesh (because then I'll need to know and define the sizes of my inlet patch).

Have a nice weekend
Christian
wese is offline   Reply With Quote

Old   April 4, 2008, 08:22
Default Hi Christian, Sorry if I ca
  #13
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18
olesen will become famous soon enough
Hi Christian,

Sorry if I cannot add anything useful, but I'm interested in the reason why you cannot get the correct bounding box.

If I understood correctly, you are using the Foam::boundBox class and constructing it with the 'doReduce' option set to true but the reduce operations: maxOp<> and/or minOp<> are somehow failing when called from a slave node?

Are you constructing your own pointField as demand-driven data or indirectly via some other class method? If so, could it be that you are demanding similar information twice within the same calculation?

This might especially be true if you are doing something like this:

something() + boundingBox(myPointField(), true);

and if you've defined 'something()' to be extra smart about creating the pointField as well.
Eg,

myClass::myPointField()
{
if (!*ptsPtr_)
{
calcPointField()
}
return *ptsPtr_;
}
and

myClass::something()
{
if (!*ptsPtr_)
{
calcPointField()
cached_ = some reduce op
}
return cached_
}

You should be careful about this sort of code giving problems and resulting in strange values for 'cached_' that could depend on the calculation order.

Since your problem is not specific to boundary conditions, can you confirm that the same problem arises elsewhere?

/mark
olesen 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
Reading and Writing Boundary Conditions chen FLUENT 3 January 21, 2008 16:17
Boundary condition - FAN Adarsh FLUENT 0 March 30, 2007 11:11
Boundary condition of the third kind or Danckwertz boundary condition plage OpenFOAM Running, Solving & CFD 4 October 3, 2006 12:21
Slip Boundary Condition for Moving Boundary Shukla Main CFD Forum 3 November 11, 2005 16:02
Writing UDF for Robbins Bounday Condition KKLAU FLUENT 1 April 6, 2004 23:58


All times are GMT -4. The time now is 05:36.