[the_ichthyologist]

Hi all,
I'm new to OpenFOAM and have very limited knowledge of C++, I'm trying to modify the volFieldValue postProcess function in order to calculate 'weighted' percentiles (weighted with alpha.water - I'm using a modified interFoam - and cells volume). I'm using openfoam6.



I think I managed to compile a modified version of volFieldValue called "myvolFieldValue" without changing anything but adding a new operation (which doesn't do anything).


Now I'm trying to implement my function: basically I want to sort "values" of concentration field in a certain region of my mesh from lowest to highest, at the same time sort corresponding cell volumes and alphas with the same order; then sum over indices (dividing by the total ) until I tick the various percentiles (0.05 , 0.15 etc.).


I modified a quicksort algorithm in c++ but I can't define the right variable type for the arrays. It would be useful to reshape the arrays to simple vectors, I don't know if this is already done by volFieldValue when 'values', 'weightField' and 'V' are passed.


I attach all relevant files, this is the main error (for now):


fieldValues/myvolFieldValue/myvolFieldValueTemplates.C:230:14: error: cannot convert ‘const Foam::Field<double>’ to ‘const Foam::Field<double>*’ for argument ‘1’ to ‘void quickSort(const Foam::Field<double>*, const Foam::Field<double>*, const Foam::Field<double>*, int, int)’
quickSort(values, weightField, V, 0, n-1);


I also don't know if quicksort function definition is at the right place in the code and if I loop correctly over the fields once they are sorted.



I apologise for my poor knowledge, any help would be much appreciated


Giovanni

[olesen]

Originally Posted by the_ichthyologist

Hi all,
...
Now I'm trying to implement my function: basically I want to sort "values" of concentration field in a certain region of my mesh from lowest to highest, at the same time sort corresponding cell volumes and alphas with the same order; then sum over indices (dividing by the total ) until I tick the various percentiles (0.05 , 0.15 etc.).

I modified a quicksort algorithm in c++ but I can't define the right variable type for the arrays. It would be useful to reshape the arrays to simple vectors, I don't know if this is already done by volFieldValue when 'values', 'weightField' and 'V' are passed.

Instead of implementing/reimplementing quick sort, it makes more sense to leverage existing code. You could use the sortedOrder method directly, or with a hand-rolled comparator
https://develop.openfoam.com/Develop...ListOps.H#L190
When doing your sorted order, you use the 'primitiveField()' method on your volume field to obtain the underlying plain 'Field' (no dimensions, no boundaries). With this sorted order (the cell Ids in the sorted order), you can then walk through your other fields. I guess you will have to do all of the binning on each local process and then sum/merge the results across all processors.

[the_ichthyologist]

Dear Mark,
Thank you very much! That's been very helpful!

I managed to make it work with this code

scalarField vals(values.component(0));
scalarField myalpha(weightField.component(0)); // not used actually
scalarField myvol(V.component(0)); // not used actually

const double allconc=gSum(vals*weightField*V);

labelList order;
sortedOrder(vals,order);


const double target=0.50; // 50th percentile (median value)

int i=0;
double temp=0;

while (temp<target)
{
temp=temp+(vals[order[i]]*weightField[order[i]]*V[order[i]])/allconc;
i++;
}
scalar& res = setComponent(result, 0);
res=vals[order[i-1]];

which gives results close to what could be true but not really: I attach median and mean graphs, which are similar in trend but not in values. There is probably something wrong in the code which I can not detect, maybe it is something clear to you. I think the operations are actually made on the region, median values are below 1 as they should, but they are also wrong.



Moreover, the run crashes (after 2438 s of computational time though) for floating point exception, but this maybe is related to solver and mesh refinement issues.


I have also two other problems:
- I don't know how to give different percentiles without defining more operations, I tried by defining more components of 'result' but it does not seem to work
- More importantly, I really don't know how to run it on parallel: the distribution sorting needs to be done over all the region, it would be much more convenient to make calculations after recontructPar and reconstructParMesh (I'm using dynamic mesh, maybe this also gives some kind of issue), but I can't find how to do functionObject post-processing after the run, I guess it's something quite feasible and you probably know


Even partial solutions, guesses, or links to useful documentation/resources are welcome.


(PS I will update this reply : the solver gives some unreasonably high values of max concentration and unreasonable low (negative) values of minimum concentration, this could vanishes with a more refined mesh; however, running on parallel is still a big issue)



Giovanni

[olesen]

Take a look at the forces function object for some ideas about "binning" into a histogram. There is a far bit of code to sift through, but sure you'll find something useful.