# Writing own boundary condition

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

 March 11, 2008, 04:10 Hi, I'm quite new to openFoam #1 New Member   Christian Wesemeyer Join Date: Mar 2009 Posts: 19 Rep Power: 10 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

 March 11, 2008, 04:25 Look how parabolic velocity do #2 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,806 Rep Power: 25 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

 March 11, 2008, 05:01 Hi Hrvoje, thanks a lot for y #3 New Member   Christian Wesemeyer Join Date: Mar 2009 Posts: 19 Rep Power: 10 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

 March 11, 2008, 05:22 vectorField will give you a li #4 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,806 Rep Power: 25 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

 March 13, 2008, 06:00 Yes: your executable does not #6 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,806 Rep Power: 25 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

 March 13, 2008, 12:04 Hi, you're right, stupid as I #7 New Member   Christian Wesemeyer Join Date: Mar 2009 Posts: 19 Rep Power: 10 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::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/cwe-1.4.1/lib/linux64GccDPOpt/squareVelocityProfile.so" #5 Foam::fvPatchField >::adddictionaryConstructorToTable::New(Foam::fvPatch const&, Foam::DimensionedField, Foam::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/cwe-1.4.1/lib/linux64GccDPOpt/squareVelocityProfile.so" #6 Foam::fvPatchField >::New(Foam::fvPatch const&, Foam::DimensionedField, Foam::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam" #7 Foam::GeometricField, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::GeometricB oundaryField(Foam::fvBoundaryMesh const&, Foam::DimensionedField, Foam::volMesh> const&, Foam::dictionary const&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam" #8 Foam::GeometricField, Foam::fvPatchField, Foam::volMesh>::readField(Foam::Istream&) in "/home/wagner/cwe/OpenFOAM/OpenFOAM-1.4.1/applications/bin/linux64GccDPOpt/interFoam" #9 Foam::GeometricField, 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

 March 14, 2008, 06:21 I've looked through the source #8 New Member   Christian Wesemeyer Join Date: Mar 2009 Posts: 19 Rep Power: 10 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

 March 14, 2008, 09:01 Finally nearly everything work #9 New Member   Christian Wesemeyer Join Date: Mar 2009 Posts: 19 Rep Power: 10 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

 April 4, 2008, 05:22 To complete this post, I would #10 New Member   Christian Wesemeyer Join Date: Mar 2009 Posts: 19 Rep Power: 10 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

 April 4, 2008, 06:13 Hi Christian I have just be #11 Senior Member   Niels Gjoel Jacobsen Join Date: Mar 2009 Location: Deltares, Delft, The Netherlands Posts: 1,758 Rep Power: 29 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.

 April 4, 2008, 07:36 Hi Niels, thanks a lot for yo #12 New Member   Christian Wesemeyer Join Date: Mar 2009 Posts: 19 Rep Power: 10 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

 April 4, 2008, 08:22 Hi Christian, Sorry if I ca #13 Senior Member   Mark Olesen Join Date: Mar 2009 Location: http://olesenm.github.io/ Posts: 802 Rep Power: 21 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

 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 chen FLUENT 3 January 21, 2008 16:17 Adarsh FLUENT 0 March 30, 2007 11:11 plage OpenFOAM Running, Solving & CFD 4 October 3, 2006 12:21 Shukla Main CFD Forum 3 November 11, 2005 16:02 KKLAU FLUENT 1 April 6, 2004 23:58

All times are GMT -4. The time now is 04:32.