[superkelle]

Hi I wrote a fvOptions scalarCodedSource and I noticed that following does not work. I mean both Pout outputs does not give the same results:

const scalar xt =1;
const scalar yt =1;


scalar testRadius;
List< label> acessedCells; 

                   

forAll (C,i)
                {   
                    const scalar & Cx =C[i].component(0);
                    const scalar & Cy =C[i].component(1);
                    const scalar & Cz =C[i].component(2);


                    if ( (pow((Cx-xt),2) + pow((Cy-yt),2) ) <= pow(rt,2)) 
                    {   
                           acessedCells.append(i); 

                           testRadius = sqrt(sqr(Cx)+sqr(Cy));
                           Pout  << "testRadius: " << testRadius << endl;

                     }
                }
                
forAll (acessedCells,i)                    

{                    

    testRadius = sqrt(sqr(C[i].component(0))+sqr(C[i].component(1)));
    Pout  << "testRadius: " << testRadius << endl;
}

I think it has something to do with the numbering of the cells, and how I acess and store it. Does it maybe also matter that I used "renumberMesh" after decomposing?
Has anybody an idea what I did wrong?

[olesen]

In the first loop, you are looping over all of the cells, so 'i' is also your cell number.
In the second loop, you are looping over the cells that you previously found, so i=0 is the first one in that list, not the cell number, which would be assessed[i].

You can avoid some of these types of things with more modern structures. Using a range-for would allow you to iterate directly over assessed Cell.
Depending of the sparsity of cells that will test positive or negative (and depending on the expected number of resizes you might also consider using a 'bitSet'. This packs on/off positions into 32bit ints, which will be better for storage than storing the labels. You can also iterate over them directly.

BTW using a List and appending to it within the loop is an excellent way to completely kill performance. Each of those append() calls implies a malloc for resizing, plus copying all of the values (each time).
If you want a list for storing things, either a DynamicList where the resizing won't kill you, or a fully dimensioned List and track the indices by hand.

[superkelle]

In the first loop, you are looping over all of the cells, so 'i' is also your cell number.
In the second loop, you are looping over the cells that you previously found, so i=0 is the first one in that list, not the cell number, which would be assessed[i].

You can avoid some of these types of things with more modern structures. Using a range-for would allow you to iterate directly over assessed Cell.
Depending of the sparsity of cells that will test positive or negative (and depending on the expected number of resizes you might also consider using a 'bitSet'. This packs on/off positions into 32bit ints, which will be better for storage than storing the labels. You can also iterate over them directly.

BTW using a List and appending to it within the loop is an excellent way to completely kill performance. Each of those append() calls implies a malloc for resizing, plus copying all of the values (each time).

Originally Posted by olesen

If you want a list for storing things, either a DynamicList where the resizing won't kill you, or a fully dimensioned List and track the indices by hand.

You are right that was a stupid mistake. Ok besides it to be slow, it should work with :

forAll (acessedCells,j)                    
{                    

    testRadius = sqrt(sqr(C[acessedCells[j]]).component(0))+sqr(C[acessedCells[j]].component(1)));
    Pout  << "testRadius: " << testRadius << endl;
}

You do not want to know how I worked around it , you would scream because of the pain.


Never the less I also like to work with the standard library std, do you think std::vector<int> would be also a good option? And work with "push_back"? I suggest not because there you have also to resize it, don't I? For me, those simple containers are easier to handle.


I did not fully understand the bitSet suggestion, I try to explain how I understood it:


So when I work with a std::bitSet I would give it at the construction the size that is equal to the total cell count, right? So it is like a very simple field of bools to check if I should do some stuff in a second loop that again iterates over all cells ? Is that not also slow when I again iterate over all cells of the mesh?

[olesen]

Test reply - forum blocked?

[olesen]

Ok with quick reply (regular reply failed)



Sorry about the delay, but holidays came in between.
For the bitSet option, we need to be sure that we are speaking of the same things. There is no std::bitSet, it is called std::bitset but that won't help you much since its size is a template parameter and not dynamic. The bitSet that I refer to (Foam::bitSet) behaves largely like a boost::dynamic_bitset. Here is a reference for you (https://www.openfoam.com/documentati...1_1bitSet.html)

I understand your preference for regular C++ std structures vs OpenFOAM. This is a valid enough criticism for people starting to work with OpenFOAM (why not simply use standard C++?). Sometimes that answer is historic, sometimes it is due to different requirements. The std::vector is quite flexible and its push_back() method will be much, much better than using a Foam::List with append() since the std:.vector usually has 1.5x resizing factor for its internal storage area, which means that it will only do a realloc and copy every now and then for the push_back(), whereas the Foam::List has no additional storage overhead for its internal storage area and would thus do a realloc for every append(). However, for most of the OpenFOAM storage we prefer to have Foam::List since we do not want this additional flexibility and memory usage when handling large meshes. Additionally, with the Foam::List, Foam::UList etc, we have the ability to steal back storage in certain places in ways that can't be managed with std::vector. If we want have a bit of both (efficient resizable and 'stealable'), we would use a Foam:ynamicList for that. It currently uses a doubling strategy for its internals (with various implementation details to be discussed elsewhere) which makes it behave the most like std::vector. Since the OpenFOAM containers have the usual begin/end methods, you can use them in most C++ algorithms without a problem.

For your current problem, this is what one solution could look like:

const scalar xt = ...;  // location
const scalar yt = ...;  // location
const scalar acceptRad2 = ...;  // Accept radius^2

const vectorField& cc = mesh.cellCentres();
bitSet accessedCells(cc.size());


forAll(cc, celli)
{
    if ((sqr(cc[celli].x() - xt) + sqr(cc[celli].y() - yt)) < acceptRad2)
    {
        accessedCells.set(celli);
    }
}

// Later

for (const label celli : accessedCells)
{
    Pout<< celli << " radius " << hypot(cc[celli].x(), cc[celli].y()) << endl;
}

If you expect a reasonably number of cells to test positively (say 10-15%) this will be a fairly reasonable approach. The second loop using the bitSet should be fairly efficient since it only iterates across items that have been set 'on'. This allows it to skip over large sections of 'off' by testing the 32 bits at once.

If you only have a few items that you expect to test positively, then using a DynamicList (or std::vector) can give you better performance. In that case, your code would look like this:

const vectorField& cc = mesh.cellCentres();

// std::vector<label> accessedCells;
DynamicList<label> accessedCells;

accessedCells.reserve(0.05*cc.size());   // pre-reserve for 5% of mesh size (if we care)


forAll(cc, celli)
{
    if ((sqr(cc[celli].x() - xt) + sqr(cc[celli].y() - yt)) < acceptRad2)
    {
        // accessedCells.push_back(celli);
        accessedCells.append(celli);
    }
}

// Later

for (const label celli : accessedCells)
{
    Pout<< celli << " radius " << hypot(cc[celli].x(), cc[celli].y()) << endl;
}

Yet another possiblity (especially if you don't care much about order):

const vectorField& cc = mesh.cellCentres();

// std::unordered_set<label> accessedCells;
labelHashSet accessedCells;

forAll(cc, celli)
{
    if ((sqr(cc[celli].x() - xt) + sqr(cc[celli].y() - yt)) < acceptRad2)
    {
        // accessedCells.insert(celli);
        accessedCells.insert(celli);
    }
}

// Later

for (const label celli : accessedCells)
{
    Pout<< celli << " radius " << hypot(cc[celli].x(), cc[celli].y()) << endl;
}

For large numbers of insertions, the Foam::HashSet will generally outperform std::unordered_set.


Sorry that there is no single answer/solution, but these depend a bit on your problem.

/mark

[olesen]

Forgot to mention some other differences. The DynamicList append() method has many more variants than std::vector push_back(). For example,

labelList list1 = ...;
labelList list2 = ...;

// Later

DynamicList<label> dynList;
dynList.append(std::move(list1));
dynList.append(std::move(list2));

Best to scan the OpenFOAM documentation.

[superkelle]

Thank you very much for your time and effort. That was very helpful, I am always stunned how much passion many of the OF users in this forum put in their answers.

[superkelle]

I noticed that there is not something like a "bitSet" in OF7 of "The OpenFOAM Foundation", I'm sure that there is some equivalent, but I think before I spend another hour for searching I'll go with the std::vector. Have a look into the implemented OF containers of both forks:


https://cpp.openfoam.org/v7/dir_19d0...e3f56cbf7.html


https://www.openfoam.com/documentati...e3f56cbf7.html

[olesen]

Originally Posted by superkelle

I noticed that there is not something like a "bitSet" in OF7 of "The OpenFOAM Foundation", I'm sure that there is some equivalent, but I think before I spend another hour for searching I'll go with the std::vector.

No you are correct, the .org version does not have bitSet, the .com version does.
Go for std::vector if that works for you, but DynamicList would be the more OpenFOAM-like equivalent for, well, a dynamically resizable list.