View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0003458OpenFOAMContributionpublic2020-02-20 10:312020-03-24 16:54
Reporterkevnor Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityN/A
Status resolvedResolutionunable to reproduce 
Product Versiondev 
Fixed in Versiondev 
Summary0003458: meshToMesh0::cellAddresses speed up
DescriptionIn the moderate-sized (~20M cell) mesh mappings I've been testing on non-convex geometries, the computation of mesh-to-mesh cell addressing in meshToMesh0::cellAddresses can be sped up significantly (in some cases by > 100x) by updating the current position 'curCell' in the source mesh to correspond to the cell returned by the octree search whenever this is performed. If this isn't done, the curCell in the source mesh can get 'stuck' at an edge, with the result that octree search is used for all the following points that are not visible from curCell, at great expense.

I've attached a patch that implements this.

Also, I'm not sure that I understand the point of the "rescue" logic in the same routine and this is also causing a slow-down in our cases! It seems that it merely avoids checking neighbour cells if the closest cell-centre is on the boundary and just goes straight to octree search. However, when meshes are graded near the boundary (as they often are), it is actually fairly common for the point to lie in a neighbour cell. It's relatively cheap to check these neighbours and expensive to do a full octree search...

A second patch is included that implements both this and the curCell update.

As I said, disabling the rescue logic results in speed-ups for the cases we're looking at, but I guess that there must be situations where it helps for it to have been introduced in the first place?! It seems to have been there since at least version 2.
TagsNo tags attached.

Activities

kevnor

2020-02-20 10:31

reporter  

meshToMesh0_noRescue.patch (7,244 bytes)    
diff --git src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
index e635732..1bb18bd 100644
--- src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
+++ src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
@@ -49,38 +49,6 @@ void Foam::meshToMesh0::calcAddressing()
     const cellList& fromCells = fromMesh_.cells();
     const pointField& fromPoints = fromMesh_.points();
 
-    // In an attempt to preserve the efficiency of linear search and prevent
-    // failure, a RESCUE mechanism will be set up. Here, we shall mark all
-    // cells next to the solid boundaries. If such a cell is found as the
-    // closest, the relationship between the origin and cell will be examined.
-    // If the origin is outside the cell, a global n-squared search is
-    // triggered.
-
-    // SETTING UP RESCUE
-
-    // visit all boundaries and mark the cell next to the boundary.
-
-    if (debug)
-    {
-        InfoInFunction << "Setting up rescue" << endl;
-    }
-
-    List<bool> boundaryCell(fromCells.size(), false);
-
-    // set reference to boundary
-    const polyPatchList& patchesFrom = fromMesh_.boundaryMesh();
-
-    forAll(patchesFrom, patchi)
-    {
-        // get reference to cells next to the boundary
-        const labelUList& bCells = patchesFrom[patchi].faceCells();
-
-        forAll(bCells, facei)
-        {
-            boundaryCell[bCells[facei]] = true;
-        }
-    }
-
     treeBoundBox meshBb(fromPoints);
 
     scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size())));
@@ -118,7 +86,6 @@ void Foam::meshToMesh0::calcAddressing()
         cellAddressing_,
         toMesh_.cellCentres(),
         fromMesh_,
-        boundaryCell,
         oc
     );
 
@@ -135,7 +102,6 @@ void Foam::meshToMesh0::calcAddressing()
                 boundaryAddressing_[patchi],
                 toPatch.faceCentres(),
                 fromMesh_,
-                boundaryCell,
                 oc
             );
         }
@@ -189,7 +155,7 @@ void Foam::meshToMesh0::calcAddressing()
 
                 scalar distSqr = sqr(wallBb.mag());
 
-                forAll(toPatch, toi)
+                for(label toi=0; toi<toPatch.size(); ++toi)
                 {
                     boundaryAddressing_[patchi][toi] = oc.findNearest
                     (
@@ -214,11 +180,10 @@ void Foam::meshToMesh0::cellAddresses
     labelList& cellAddressing_,
     const pointField& points,
     const fvMesh& fromMesh,
-    const List<bool>& boundaryCell,
     const indexedOctree<treeDataCell>& oc
 ) const
 {
-    // the implemented search method is a simple neighbour array search.
+    // The implemented search method is a simple neighbour array search.
     // It starts from a cell zero, searches its neighbours and finds one
     // which is nearer to the target point than the current position.
     // The location of the "current position" is reset to that cell and
@@ -275,64 +240,55 @@ void Foam::meshToMesh0::cellAddresses
         }
         else
         {
-            // If curCell is a boundary cell then the point maybe either outside
-            // the domain or in an other region of the doamin, either way use
-            // the octree search to find it.
-            if (boundaryCell[curCell])
+            // Not in nearest cell: search the neighbours
+            bool found = false;
+
+            // set the current list of neighbouring cells
+            const labelList& neighbours = cc[curCell];
+
+            forAll(neighbours, nI)
             {
-                cellAddressing_[toI] = oc.findInside(p);
+                // search through all the neighbours.
+                if (fromMesh.pointInCell(p, neighbours[nI]))
+                {
+                    cellAddressing_[toI] = neighbours[nI];
+                    found = true;
+                    break;
+                }
             }
-            else
+
+            if (!found)
             {
-                // If not on the boundary search the neighbours
-                bool found = false;
+                // If still not found search the neighbour-neighbours
 
                 // set the current list of neighbouring cells
                 const labelList& neighbours = cc[curCell];
 
                 forAll(neighbours, nI)
                 {
-                    // search through all the neighbours.
-                    // If point is in neighbour reset current cell
-                    if (fromMesh.pointInCell(p, neighbours[nI]))
-                    {
-                        cellAddressing_[toI] = neighbours[nI];
-                        found = true;
-                        break;
-                    }
-                }
-
-                if (!found)
-                {
-                    // If still not found search the neighbour-neighbours
+                    // set the current list of neighbour-neighbouring cells
+                    const labelList& nn = cc[neighbours[nI]];
 
-                    // set the current list of neighbouring cells
-                    const labelList& neighbours = cc[curCell];
-
-                    forAll(neighbours, nI)
+                    forAll(nn, nI)
                     {
-                        // set the current list of neighbour-neighbouring cells
-                        const labelList& nn = cc[neighbours[nI]];
-
-                        forAll(nn, nI)
+                        // search through all the neighbours.
+                        if (fromMesh.pointInCell(p, nn[nI]))
                         {
-                            // search through all the neighbours.
-                            // If point is in neighbour reset current cell
-                            if (fromMesh.pointInCell(p, nn[nI]))
-                            {
-                                cellAddressing_[toI] = nn[nI];
-                                found = true;
-                                break;
-                            }
+                            cellAddressing_[toI] = nn[nI];
+                            found = true;
+                            break;
                         }
-                        if (found) break;
                     }
+                    if (found) break;
                 }
+            }
 
-                if (!found)
-                {
-                    // Still not found so us the octree
-                    cellAddressing_[toI] = oc.findInside(p);
+            if (!found)
+            {
+                // Still not found so use the octree
+                cellAddressing_[toI] = oc.findInside(p);
+                if(cellAddressing_[toI] != -1) {
+                    curCell = cellAddressing_[toI];
                 }
             }
         }
diff --git src/sampling/meshToMesh0/meshToMesh0.H src/sampling/meshToMesh0/meshToMesh0.H
index b20e7f0..5a134f8 100644
--- src/sampling/meshToMesh0/meshToMesh0.H
+++ src/sampling/meshToMesh0/meshToMesh0.H
@@ -110,7 +110,6 @@ class meshToMesh0
             labelList& cells,
             const pointField& points,
             const fvMesh& fromMesh,
-            const List<bool>& boundaryCell,
             const indexedOctree<treeDataCell>& oc
         ) const;
meshToMesh0_noRescue.patch (7,244 bytes)   
meshToMesh0_updateCurCell.patch (1,083 bytes)    
diff --git src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
index e635732..a9b31fa 100644
--- src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
+++ src/sampling/meshToMesh0/calculateMeshToMesh0Addressing.C
@@ -281,6 +281,9 @@ void Foam::meshToMesh0::cellAddresses
             if (boundaryCell[curCell])
             {
                 cellAddressing_[toI] = oc.findInside(p);
+                if(cellAddressing_[toI] != -1) {
+                    curCell = cellAddressing_[toI];
+		}
             }
             else
             {
@@ -331,8 +334,11 @@ void Foam::meshToMesh0::cellAddresses
 
                 if (!found)
                 {
-                    // Still not found so us the octree
+                    // Still not found so use the octree
                     cellAddressing_[toI] = oc.findInside(p);
+                    if(cellAddressing_[toI] != -1) {                        
+                        curCell = cellAddressing_[toI];
+		    }
                 }
             }
         }
meshToMesh0_updateCurCell.patch (1,083 bytes)   

henry

2020-02-21 15:35

manager   ~0011205

Last edited: 2020-02-21 15:35

 For what range of cases have you tested these patches?

kevnor

2020-02-21 17:32

reporter   ~0011208

So far I've only tested it on 4 different cases here (plus those tutorial cases that use mapFields, though those are all rather small). Our cases include tet and hex meshes between 2 million and 30 million cells and a variety of different geometries. I get widely varying speedups for the different cases, but the patched versions are faster than the unpatched for all the larger cases. I can post more exact figures if its useful?

Unfortunately I'm not able to upload the meshes I've tested so far here, but I will try and put together some simple blockMesh / snappyHexMesh meshes that demonstrate the issue shortly!

I understand this is going to need further testing to make sure that it doesn't break anything -- in particular, as I mentioned above, I guess that the 'rescue' logic must have been put in for a good reason.

henry

2020-03-20 21:20

manager   ~0011264

Last edited: 2020-03-20 21:32

 On the complex geometry cases I have tested meshToMesh0_updateCurCell.patch gives a slight improvement in performance but meshToMesh0_noRescue.patch gives a very significant drop in performance. The results of the mapping are the same for the current version and both patches.

henry

2020-03-20 21:24

manager   ~0011265

Can you provide any case which shows the benefit of the either or both of the two patches?

henry

2020-03-24 16:54

manager   ~0011275

I have committed the "curCell" optimisation as it shows a minor improvement in performance for some cases.

commit 4be01b4e70e4bab877346cce52900a59fd6bfeac

However removing the "rescue" logic has a significant detrimental effect for all the cases I have tested and I have not committed it.
We have not found any cases for which either change provides significant benefit and we need you to provide a means to reproduce your findings in order to work further on this.

Issue History

Date Modified Username Field Change
2020-02-20 10:31 kevnor New Issue
2020-02-20 10:31 kevnor File Added: meshToMesh0_noRescue.patch
2020-02-20 10:31 kevnor File Added: meshToMesh0_updateCurCell.patch
2020-02-21 15:35 henry Note Added: 0011205
2020-02-21 15:35 henry Note Edited: 0011205
2020-02-21 17:32 kevnor Note Added: 0011208
2020-03-20 21:20 henry Note Added: 0011264
2020-03-20 21:24 henry Note Added: 0011265
2020-03-20 21:32 henry Note Edited: 0011264
2020-03-24 16:54 henry Assigned To => henry
2020-03-24 16:54 henry Status new => resolved
2020-03-24 16:54 henry Resolution open => unable to reproduce
2020-03-24 16:54 henry Fixed in Version => dev
2020-03-24 16:54 henry Note Added: 0011275