[wese]

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

[hjasak]

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

[wese]

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

[hjasak]

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

[wese]

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

[hjasak]

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

[wese]

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]

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]

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]

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

[ngj]

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

[wese]

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

[olesen]

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