CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Bugs

mapFields major bug

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

Like Tree3Likes
  • 2 Post By alchem
  • 1 Post By wyldckat

Reply
 
LinkBack Thread Tools Display Modes
Old   October 14, 2017, 05:53
Exclamation mapFields major bug
  #1
New Member
 
Join Date: May 2012
Posts: 4
Rep Power: 8
alchem is on a distinguished road
I tackled a really annoying bug of mapFields utility. mapFields is used when you need to interpolate fields from one mesh to another, particularly from crude mesh to refined one. In my case this utility sometimes throws this error:


Code:
Mapping fields for time 0.0035
 
 
     interpolating p
 #0  Foam::error::printStack(Foam::Ostream&) at ??:?
 #1  Foam::sigSegv::sigHandler(int) at ??:?
 #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
 #3  ? at ??:?
 #4  ? at ??:?
 #5  ? at ??:?
 #6  ? at ??:?
 #7  ? at ??:?
 #8  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
 #9  ? at ??:?
 Segmentation fault (core dumped)
Most of the times I thought that mesh I interpolate to is at fault. I mostly work with hexahedral or polyhedral meshes which is hard to create. After rebuilding mesh several times mapFields usually started to work. However one day I got similar error on tetrahedral mesh, simplest as mesh can be and rebuilding of that mesh with different sizes did not help. It became obvious that there is bug somwhere in mapFields utility.
After some time of searching through OpenFOAM source I found reason why segmentation fault throws:


src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C
Code:
void Foam::meshToMesh0::calculateInverseDistanceWeights() const
 ```
 label directCelli = -1;
 if (m < directHitTol || neighbours.empty())
 {
     directCelli = celli;
 }
 else
 {
     forAll(neighbours, ni)
     {
         scalar nm = mag(target - centreFrom[neighbours[ni]]);
         if (nm < directHitTol)
         {
             directCelli = neighbours[ni];
             break;
         }
     }
 }
 
 
 if (directCelli != -1)
 {
     // Direct hit
     invDistCoeffs[directCelli].setSize(1);
     invDistCoeffs[directCelli][0] = 1.0;
     V_ += fromMesh_.V()[cellAddressing_[directCelli]];
 }
 else
 {
 ...
 }
 ```
OpenFOAM uses reverse distance interpolation to interpolate field values. This interpolation will crash if point from which we interpolate and point interpolate to overlaps. So in OpenFOAM code this interpolation used only if distance between points is larger than directHitTol variable. If points is closer it is considered a direct hit and field value is copied from source point. At least as I imagined how it supposed to work.

In OpenFOAM implementation there is some non trivial logic. When we get into calculateInverseDistanceWeights function we know source cell we interpolate from. It is in cellAddressing_[celli]. So we calculate distance between centers of source and destination cells – variable m. Then we check if m is less than directHitTol. If it is, then it is a direct hit and we don’t need to do reverse distance interpolation, so we set invDistCoeffs[directCelli] with 1.0 where directCelli=celli and everything works fine. If distance m is bigger than directHitTol neighbors get checked. Thats where bug happens. If neighbor cell center for some reason is closer than directHitTol, invDistCoeffs[directCelli] will be initialized for neighbor cell while it should be initialized for celli we work with.

When invDistCoeffs get returned to interpolateField function it has some elements uninitialized and generate segmentation fault here:

src/sampling/meshToMesh0/meshToMesh0Templates.C

Code:
template<class Type, class CombineOp>
 void Foam::meshToMesh0::interpolateField
```
 const labelList& neighbours = cc[adr[celli]];
 const scalarList& w = weights[celli];
 
 
 Type f = fromVf[adr[celli]]*w[0];
 
 
 for (label ni = 1; ni < w.size(); ni++)
 {
     f += fromVf[neighbours[ni - 1]]*w[ni];
 }
 ```
w[0] is trying to access element of empty array.

There is another connected less critical bug.
directHitTol is hardcoded in meshToMesh0.C.


const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-05;


Yeah, just like that. If your mesh edge size is less than 1e-05 you will always get direct hit. Thus simplest fix to this bug is to set directHitTol lower than your edge length. Something like this


const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-10;


Oh, and there is another hard coded directHitTol in PatchToPatchInterpolation.C. I hope it is not affecting my calculations somehow...

I attached simple test case.

To check it execute following commands:
cd destination
mapFields ../source

P.S. Bug occurs on foundation OpenFOAM 3.1. OpenFOAM 5.0 has similar implementation of calculateInverseDistanceWeights so I think it will crash all the same
Attached Files
File Type: zip testCase.zip (173.5 KB, 8 views)
kindle and melabbassi like this.

Last edited by alchem; October 15, 2017 at 12:00.
alchem is offline   Reply With Quote

Old   December 21, 2017, 13:34
Default
  #2
New Member
 
Patrick H.
Join Date: Aug 2017
Location: Germany
Posts: 3
Rep Power: 3
patrex is on a distinguished road
The same is going to happen in OpenFOAM 4.1.

After I changed:

const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-05;

to:

const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-10;

like you said, it worked fine.

The bug is really annoying, because it's almost happend all the time at my meshes. So thank you for the examination and the provisional solution.
patrex is offline   Reply With Quote

Old   July 18, 2018, 06:22
Default
  #3
Member
 
Sebastian Trunk
Join Date: Mar 2015
Location: Erlangen, Germany
Posts: 57
Rep Power: 5
sisetrun is on a distinguished road
Thanks, worked for me as well...
sisetrun is offline   Reply With Quote

Old   July 23, 2018, 10:01
Default
  #4
New Member
 
Mohamed el Abbassi
Join Date: Oct 2016
Location: Delft, the Netherlands
Posts: 6
Rep Power: 4
melabbassi is on a distinguished road
Same issue with OF-v1712+. Your solution to reduce directHitTol worked. Thanks!
melabbassi is offline   Reply With Quote

Old   August 2, 2018, 09:34
Default
  #5
New Member
 
Join Date: Jun 2018
Posts: 10
Rep Power: 2
xshmuel is on a distinguished road
he guys I am getting more and more acquainted with openfoam. I have a question if you change the meshTOMesh0 file what parts do you then have to recompile? do you have to recompile the entire solver?

thanks for the help guys
xshmuel is offline   Reply With Quote

Old   August 2, 2018, 09:38
Default
  #6
New Member
 
Mohamed el Abbassi
Join Date: Oct 2016
Location: Delft, the Netherlands
Posts: 6
Rep Power: 4
melabbassi is on a distinguished road
Quote:
Originally Posted by xshmuel View Post
he guys I am getting more and more acquainted with openfoam. I have a question if you change the meshTOMesh0 file what parts do you then have to recompile? do you have to recompile the entire solver?

thanks for the help guys

You can run Allwmake again on the root directory. Only the modified files will be recompiled.
melabbassi is offline   Reply With Quote

Old   October 3, 2018, 23:54
Default
  #7
New Member
 
Sam
Join Date: Feb 2018
Posts: 4
Rep Power: 2
xevious is on a distinguished road
Quote:
Originally Posted by alchem View Post
I tackled a really annoying bug of mapFields utility. mapFields is used when you need to interpolate fields from one mesh to another, particularly from crude mesh to refined one. In my case this utility sometimes throws this error:


Code:
Mapping fields for time 0.0035
 
 
     interpolating p
 #0  Foam::error::printStack(Foam::Ostream&) at ??:?
 #1  Foam::sigSegv::sigHandler(int) at ??:?
 #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
 #3  ? at ??:?
 #4  ? at ??:?
 #5  ? at ??:?
 #6  ? at ??:?
 #7  ? at ??:?
 #8  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
 #9  ? at ??:?
 Segmentation fault (core dumped)
Most of the times I thought that mesh I interpolate to is at fault. I mostly work with hexahedral or polyhedral meshes which is hard to create. After rebuilding mesh several times mapFields usually started to work. However one day I got similar error on tetrahedral mesh, simplest as mesh can be and rebuilding of that mesh with different sizes did not help. It became obvious that there is bug somwhere in mapFields utility.
After some time of searching through OpenFOAM source I found reason why segmentation fault throws:


src/sampling/meshToMesh0/calculateMeshToMesh0Weights.C
Code:
void Foam::meshToMesh0::calculateInverseDistanceWeights() const
 ```
 label directCelli = -1;
 if (m < directHitTol || neighbours.empty())
 {
     directCelli = celli;
 }
 else
 {
     forAll(neighbours, ni)
     {
         scalar nm = mag(target - centreFrom[neighbours[ni]]);
         if (nm < directHitTol)
         {
             directCelli = neighbours[ni];
             break;
         }
     }
 }
 
 
 if (directCelli != -1)
 {
     // Direct hit
     invDistCoeffs[directCelli].setSize(1);
     invDistCoeffs[directCelli][0] = 1.0;
     V_ += fromMesh_.V()[cellAddressing_[directCelli]];
 }
 else
 {
 ...
 }
 ```
OpenFOAM uses reverse distance interpolation to interpolate field values. This interpolation will crash if point from which we interpolate and point interpolate to overlaps. So in OpenFOAM code this interpolation used only if distance between points is larger than directHitTol variable. If points is closer it is considered a direct hit and field value is copied from source point. At least as I imagined how it supposed to work.

In OpenFOAM implementation there is some non trivial logic. When we get into calculateInverseDistanceWeights function we know source cell we interpolate from. It is in cellAddressing_[celli]. So we calculate distance between centers of source and destination cells Ė variable m. Then we check if m is less than directHitTol. If it is, then it is a direct hit and we donít need to do reverse distance interpolation, so we set invDistCoeffs[directCelli] with 1.0 where directCelli=celli and everything works fine. If distance m is bigger than directHitTol neighbors get checked. Thats where bug happens. If neighbor cell center for some reason is closer than directHitTol, invDistCoeffs[directCelli] will be initialized for neighbor cell while it should be initialized for celli we work with.

When invDistCoeffs get returned to interpolateField function it has some elements uninitialized and generate segmentation fault here:

src/sampling/meshToMesh0/meshToMesh0Templates.C

Code:
template<class Type, class CombineOp>
 void Foam::meshToMesh0::interpolateField
```
 const labelList& neighbours = cc[adr[celli]];
 const scalarList& w = weights[celli];
 
 
 Type f = fromVf[adr[celli]]*w[0];
 
 
 for (label ni = 1; ni < w.size(); ni++)
 {
     f += fromVf[neighbours[ni - 1]]*w[ni];
 }
 ```
w[0] is trying to access element of empty array.

There is another connected less critical bug.
directHitTol is hardcoded in meshToMesh0.C.


const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-05;


Yeah, just like that. If your mesh edge size is less than 1e-05 you will always get direct hit. Thus simplest fix to this bug is to set directHitTol lower than your edge length. Something like this


const Foam::scalar Foam::meshToMesh0::directHitTol = 1e-10;


Oh, and there is another hard coded directHitTol in PatchToPatchInterpolation.C. I hope it is not affecting my calculations somehow...

I attached simple test case.

To check it execute following commands:
cd destination
mapFields ../source

P.S. Bug occurs on foundation OpenFOAM 3.1. OpenFOAM 5.0 has similar implementation of calculateInverseDistanceWeights so I think it will crash all the same



Hey, I tried your recommendation for reducing the directHitTol in the meshToMesh0.C file... i got the same error again. Which directory is the patchToPatchInterpolation.C file?
xevious is offline   Reply With Quote

Old   October 14, 2018, 09:29
Default
  #8
Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,453
Blog Entries: 41
Rep Power: 115
wyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of lightwyldckat is a glorious beacon of light
Greetings to all,

@xevious: You can find the file by running this command:
Code:
find $FOAM_SRC -iname patchToPatchInterpolation.C
The 'i' in '-iname' is so that it ignores the letter case.


@Everyone else: At first I was going to relay this bug report to the actual issue tracker here: https://bugs.openfoam.org - but I then first tested the provided cases with the utility mapFieldsPar, which worked without crashing... which as lead me to not yet report this bug.

So the story here is that the utility mapFieldsPar was designed to run in parallel and should not have the limitations that the current mapFields utility has.
In other words, mapFields was in fact replaced with mapFieldsPar sometime ago, but then they had to undo that because the utility mapFieldsPar also has a few limitations:

So my request is that, if possible, please try using mapFieldsPar as an alternative to the fix proposed by alchem. And if that still doesn't solve the problem, then please let us know, so that this issue can be reported at https://bugs.openfoam.org and subsequently fixed.

What I mean is that given that this is a bug in a part of the code that is meant to be discarded sometime in the near future, it's a bit annoying having to fix a bug in something that will not be of much use in the long term
However, if mapFieldsPar does not help solve this issue, then please let us know, so that it can be fixed, even if the code may be removed sometime in future releases.

Best regards,
Bruno
sisetrun likes this.
__________________
wyldckat is offline   Reply With Quote

Old   November 9, 2018, 11:03
Default weird behavior indeed
  #9
New Member
 
Join Date: Aug 2016
Posts: 11
Blog Entries: 68
Rep Power: 4
kindle is on a distinguished road
I found OpenFOAM/4.0-foss-2016b (easybuild) works fine while it is true that 2.3.x mapFields just cannot avoid even if I tried the above proposal

#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigSegv::sigHandler(int) at ??:?
#2 in "/lib64/libc.so.6"
#3
at ??:?
#4
at ??:?
#5
at ??:?
#6
at ??:?
#7
at ??:?
#8 __libc_start_main in "/lib64/libc.so.6"
#9
at ??:?
./mapFields.sh: line 8: 17987 Segmentation fault (core dumped)


Here are the details that I tried to modify these hard-coded const.

One in
Code:
OpenFOAM/sampling/.../patchToPatchInterpolation.C
another mentioned above
Code:
sampling/meshToMeshInterpolation/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
I re-compiled the corresponding lOpenFOAM and lsampling (without wclean)
I recompiled then the mapFields utility by wclean ; wmake

But it didn't work.

Then I change to OF-4 that I installed, it works just fine.

Quote:
Originally Posted by wyldckat View Post
Greetings to all,

@xevious: You can find the file by running this command:
Code:
find $FOAM_SRC -iname patchToPatchInterpolation.C
The 'i' in '-iname' is so that it ignores the letter case.


@Everyone else: At first I was going to relay this bug report to the actual issue tracker here: https://bugs.openfoam.org - but I then first tested the provided cases with the utility mapFieldsPar, which worked without crashing... which as lead me to not yet report this bug.

So the story here is that the utility mapFieldsPar was designed to run in parallel and should not have the limitations that the current mapFields utility has.
In other words, mapFields was in fact replaced with mapFieldsPar sometime ago, but then they had to undo that because the utility mapFieldsPar also has a few limitations:

So my request is that, if possible, please try using mapFieldsPar as an alternative to the fix proposed by alchem. And if that still doesn't solve the problem, then please let us know, so that this issue can be reported at https://bugs.openfoam.org and subsequently fixed.

What I mean is that given that this is a bug in a part of the code that is meant to be discarded sometime in the near future, it's a bit annoying having to fix a bug in something that will not be of much use in the long term
However, if mapFieldsPar does not help solve this issue, then please let us know, so that it can be fixed, even if the code may be removed sometime in future releases.

Best regards,
Bruno
kindle is offline   Reply With Quote

Old   November 10, 2018, 10:59
Default
  #10
New Member
 
Sam
Join Date: Feb 2018
Posts: 4
Rep Power: 2
xevious is on a distinguished road
Quote:
Originally Posted by wyldckat View Post
Greetings to all,

@xevious: You can find the file by running this command:
Code:
find $FOAM_SRC -iname patchToPatchInterpolation.C
The 'i' in '-iname' is so that it ignores the letter case.


@Everyone else: At first I was going to relay this bug report to the actual issue tracker here: https://bugs.openfoam.org - but I then first tested the provided cases with the utility mapFieldsPar, which worked without crashing... which as lead me to not yet report this bug.

So the story here is that the utility mapFieldsPar was designed to run in parallel and should not have the limitations that the current mapFields utility has.
In other words, mapFields was in fact replaced with mapFieldsPar sometime ago, but then they had to undo that because the utility mapFieldsPar also has a few limitations:

So my request is that, if possible, please try using mapFieldsPar as an alternative to the fix proposed by alchem. And if that still doesn't solve the problem, then please let us know, so that this issue can be reported at https://bugs.openfoam.org and subsequently fixed.

What I mean is that given that this is a bug in a part of the code that is meant to be discarded sometime in the near future, it's a bit annoying having to fix a bug in something that will not be of much use in the long term
However, if mapFieldsPar does not help solve this issue, then please let us know, so that it can be fixed, even if the code may be removed sometime in future releases.

Best regards,
Bruno



Hello, Thanks for the reply!


I tried all these methods, and still got this error with mapFieldsPar.


Create meshes

Source mesh size: 11419078 Target mesh size: 20750858


Consistently creating and mapping fields for time 0.08

Creating mesh-to-mesh addressing for region0 and region0 regions using cellVolumeWeight
[2] #0 Foam::error:rintStack(Foam::Ostream&) addr2line failed
[2] #1 Foam::sigSegv::sigHandler(int) addr2line failed
[2] #2 ? addr2line failed
[2] #3 Foam::meshToMesh::distributeAndMergeCells(Foam::ma pDistribute const&, Foam:olyMesh const&, Foam::globalIndex const&, Foam::Field<Foam::Vector<double> >&, Foam::List<Foam::face>&, Foam::List<int>&, Foam::List<int>&, Foam::List<int>&) const addr2line failed
[2] #4 Foam::meshToMesh::calculate(Foam::word const&, bool) addr2line failed
[2] #5 Foam::meshToMesh::constructNoCuttingPatches(Foam:: word const&, Foam::word const&, bool) addr2line failed
[2] #6 Foam::meshToMesh::meshToMesh(Foam:olyMesh const&, Foam:olyMesh const&, Foam::word const&, Foam::word const&, bool) addr2line failed
[2] #7 ?
[2] #8 ?
[2] #9 __libc_start_main addr2line failed
[2] #10 ?
[LAPTOP-DPO0FDTG:29404] *** Process received signal ***
[LAPTOP-DPO0FDTG:29404] Signal: Segmentation fault (11)
[LAPTOP-DPO0FDTG:29404] Signal code: (-6)
[LAPTOP-DPO0FDTG:29404] Failing at address: 0x3e8000072dc
[LAPTOP-DPO0FDTG:29404] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x354b0)[0x7fbd5f8754b0]
[LAPTOP-DPO0FDTG:29404] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x38)[0x7fbd5f875428]
[LAPTOP-DPO0FDTG:29404] [ 2] /lib/x86_64-linux-gnu/libc.so.6(+0x354b0)[0x7fbd5f8754b0]
[LAPTOP-DPO0FDTG:29404] [ 3] /opt/OpenFOAM/OpenFOAM-v1712/platforms/linux64Gcc63DPInt32Opt/lib/libsampling.so(_ZNK4Foam10meshToMesh23distributeAn dMergeCellsERKNS_13mapDistributeERKNS_8polyMeshERK NS_11globalIndexERNS_5FieldINS_6VectorIdEEEERNS_4L istINS_4faceEEERNSF_IiEESK_SK_+0xab8)[0x7fbd6421c148]
[LAPTOP-DPO0FDTG:29404] [ 4] /opt/OpenFOAM/OpenFOAM-v1712/platforms/linux64Gcc63DPInt32Opt/lib/libsampling.so(_ZN4Foam10meshToMesh9calculateERKNS _4wordEb+0x32c)[0x7fbd64201ddc]
[LAPTOP-DPO0FDTG:29404] [ 5] /opt/OpenFOAM/OpenFOAM-v1712/platforms/linux64Gcc63DPInt32Opt/lib/libsampling.so(_ZN4Foam10meshToMesh25constructNoCu ttingPatchesERKNS_4wordES3_b+0x2b2)[0x7fbd642040c2]
[LAPTOP-DPO0FDTG:29404] [ 6] /opt/OpenFOAM/OpenFOAM-v1712/platforms/linux64Gcc63DPInt32Opt/lib/libsampling.so(_ZN4Foam10meshToMeshC2ERKNS_8polyMe shES3_RKNS_4wordES6_b+0x107)[0x7fbd64204607]
[LAPTOP-DPO0FDTG:29404] [ 7] mapFieldsPar[0x44c351]
[LAPTOP-DPO0FDTG:29404] [ 8] mapFieldsPar[0x42b6be]
[LAPTOP-DPO0FDTG:29404] [ 9] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf0)[0x7fbd5f860830]
[LAPTOP-DPO0FDTG:29404] [10] mapFieldsPar[0x42c5e9]
[LAPTOP-DPO0FDTG:29404] *** End of error message ***
[0] #0 Foam::error:rintStack(Foam::Ostream&)[1] #0 Foam::error:rintStack(Foam::Ostream&)--------------------------------------------------------------------------
mpirun noticed that process rank 2 with PID 29404 on node LAPTOP-DPO0FDTG exited on signal 11 (Segmentation fault).
--------------------------------------------------------------------------
xevious 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
implementation of mapFields into parallel transient case simpomann OpenFOAM Pre-Processing 4 August 2, 2016 05:41
Personalization of mapFields and libsampling - Compilation issues saimat OpenFOAM Programming & Development 3 June 29, 2016 09:56
mapFields for 3D praveensrikanth91 OpenFOAM Pre-Processing 3 February 17, 2015 06:23
Strange random behaviour of mapFields blaise OpenFOAM Pre-Processing 0 November 3, 2014 10:37
The -parallel parameter of mapFields utility in OpenFOAM v2.3.0 shuoxue OpenFOAM Pre-Processing 1 April 28, 2014 06:59


All times are GMT -4. The time now is 22:14.