[kmooney]

Howdy Foamers,

I'm currently working on adding a load balancer to run with dynamicRefineFvMesh. I feel that adding such a capability will vastly improve the capabilities of openFoam in terms of scaling and efficiency. I'm actually rather surprised that is has yet to be implemented considering dynamic adaptive mesh refinement is such a useful tool. For anyone not familiar with what I'm talking about check out the interDyMFoam damBreak tutorial in 2.1.x (I think the 1.6-ext has an issue).

Currently, if you perform a re-balance mid-run, the refinementHistory data is lost and all cells are considered un-refined. By the next time step the refiner starts refining cells it shouldn't be and everything goes to heck.

My question that I'm putting out there is:

1: Is anyone out there familiar with the oct-tree storage style of the 'refinementHistory' structures? At first I thought I understood it but after performing some simple 2 cell refinement studies I realize I'm missing something.

2: Has anyone worked with the agglomeration maps that are built into the 'decomposer' functions? From what I can tell, if we move a parent cell and all of its oct-tree refined children together during a balance, we can maintain the validity of the refinementHistory data.
see: Foam::refinementHistory::distribute(mapDistributeP olyMesh) in OpenFOAM/src/dynamicMesh/lnInclude/refinementHistory.H in OF2.1.x. If we can work with the refinementHistory correctly, we can build the agglomeration map, and force Scotch to keep the cells together.

I know there are a lot of people out there that need to be able to run large scale AMR simulations. I feel like if just a few of us stick our heads together for a few days we'll be able to set this up and really expand and accelerate our work.

Cheers!!

Kyle

P.S. Sorry I spelled agglomeration wrong in the title of the thread...

[wyldckat]

Greetings Kyle,

I've fixed the typo in the title

As for the issue at hand, I know that mturcios777 was working a lot lately with dynamic refinement and triggered on a few limitations that have already been fixed. Check his posts here on the forum, as well as his reports on the bug tracker:

• http://www.openfoam.org/mantisbt/view.php?id=653
• http://www.openfoam.org/mantisbt/view.php?id=638

From one of the comments on the tracker, you'll see he should be back in about a week from now.

Good luck!
Bruno

[kmooney]

Hi Bruno,

Thanks for the info! I actually ran into the same issues that he reported just recently when I tried to combine a moving mesh case with adaptive hex-refinement. Until now I though it was a bug in my code but it appears it's been addressed in the foam source.

[tomislav_maric]

Originally Posted by kmooney

1: Is anyone out there familiar with the oct-tree storage style of the 'refinementHistory' structures? At first I thought I understood it but after performing some simple 2 cell refinement studies I realize I'm missing something.

Hi Kyle,

I'm trying to figure this out as well. Here's how far I've got:

It seems that both the visibleCells and splitCells relate to the current mesh. If the cell is visible by the refinement engine, it is divided by the hexRef8 cutter. This is recognized by the fact that the visibleCell label is different than -1. This suggests that the visibleCells index list is as long as the mesh.cells() for the current mesh state. I have checked this and it is right.

Now, the index of the visibleCells corresponds to the splitCells if it is different than -1. Inside split cells on this index, there lies a parent index and a label list of the sub-cells.

Parent index:
The index of the split cells will be -1 if it has no parent. If it has a parent, it means that the refinement layer is greater than 1, and the parent_ index (public attribute) will correspond to another element of splitCells (not the unrefined mesh).

Label list:
These are the indices (again to splitCells) of the sub-cells of the split cell.

So, if e.g. we would have a refined cell in our mesh, with a second level of refinement, the visibleCells index for it would be different than -1.

In case of two levels of refinement for a single cell,

the splitCells may look like this:

splitCells
(
3 () // 0 Unrefined subcell: parent 3
3 () // 1 Unrefined subcell: parent 3
3 () // 2 -//-
-1 (0 1 2 5 6 7 8 9) // 3 Refined cell, top level of refinement, no parent
3 (10 11 12 13 14 15 16 17) // 4 Refined subcell, level 2, parent 3
3 () // 5 Unrefined subcell: parent 3
3 () // 6 -//-
3 () // 7 -//-
3 () // 8 -//-
9 () // 9 -//-
4 () // 10 Unrefined subcell: parent 4
4 () // 11 -//-
4 () // 12  -//-
4 () // 13  -//-
4 () // 14  -//-
4 () // 15  -//-
4 () // 16  -//-
4 () // 17  -//-
)

Coolio. Okay, what I would need is the label of the parent cell coming from the mesh before refinement is executed, and this is not kept by the refinementHistory. It is generated by the mesh cutter, and obtained as an autoPtr<mapPolyMesh> within the dynamicRefineFvMesh::update() method. There is no access to the map available from the mesh. The map is computed during refinement and dies directly after... So, nothing else to do but to copy paste code from the ::update() and shove my mapper in there in a new class. Copy pasting is the root of all debugging evil, so hence this post.. any suggestions?

Tomislav

[kmooney]

Do you feel that given a current hex8 refined mesh and some refinement history, we would be able to create a routine like this:

//given some current mesh's cell index cellI

label topLevelParentIndex = returnTopLevelParent(cellI);

What I think we need to do to perform the load balance is to bin up all the cells that originate from the same top level parent from the start of the run. Lets call these t=0 top level cells superParents. If we can agglomerate cells by superParent bins, we could instruct ptScotch to keep each bin of cells together through the load balance. From what I gather from the refinementHistory source, we can then balance the history trees using the same fvDistribute map as long as the trees are not broken through the transfer.

I've tried to quickly write a recursive returnTopLevelParent(cellI) function but I'm not getting the correct child->superParent maps yet.

What do you think?

[tomislav_maric]

Originally Posted by kmooney

Do you feel that given a current hex8 refined mesh and some refinement history, we would be able to create a routine like this:

//given some current mesh's cell index cellI

label topLevelParentIndex = returnTopLevelParent(cellI);

What I think we need to do to perform the load balance is to bin up all the cells that originate from the same top level parent from the start of the run. Lets call these t=0 top level cells superParents. If we can agglomerate cells by superParent bins, we could instruct ptScotch to keep each bin of cells together through the load balance. From what I gather from the refinementHistory source, we can then balance the history trees using the same fvDistribute map as long as the trees are not broken through the transfer.

I've tried to quickly write a recursive returnTopLevelParent(cellI) function but I'm not getting the correct child->superParent maps yet.

What do you think?

No chance getting this from refinementHistory, hexRef8, or anything other, than the mapPolyMesh. I managed to inherit from dynamicRefineFvMesh and make it store an autoPtr<mapPolyMesh> and provide it (via a ref of course) by a public member function. I can send you the code without a problem, or better yet, let's set up a git repo on bitbucket for this!

This info allows me to actually address the true parent of the refined cell, not the index from splitCells. I need this for the geometrical mapping procedure for the volume fraction field, that I got to work today.

Since this mesh is storing an instance of the mapPolyMesh, we can log the changes from the upper layer of the code using another class then, and construct the bins you are talking about and store them there as the mesh evolves, if this is what you mean...

In the same way my geomInterface class now talks to the mesh, asks about the cell labels of the refined cells in the mesh before refinement is done, and does the mapping..

[kmooney]

You don't think that a recursive style call of

refinementHistory::parentIndex(cellI)

would return a superParent index?

[tomislav_maric]

Originally Posted by kmooney

You don't think that a recursive style call of

refinementHistory::parentIndex(cellI)

would return a superParent index?

It will for the current refinement:

        label parentIndex(const label cellI) const
        {
            label index = visibleCells_[cellI];

            if (index < 0)
            {
                FatalErrorIn("refinementHistory::parentIndex(const label)")
                    << "Cell " << cellI << " is not visible"
                    << abort(FatalError);
            }
            return splitCells_[index].parent_;
        }

At the end, you are given a label from splitCells, and from refinementHistory.H:

All refinement history. Used in unrefinement.
- visibleCells: valid for the current mesh and contains per cell -1 (cell unrefined) or an index into splitCells_.

• splitCells: for every split contains the parent (also index into splitCells) and optionally a subsplit as 8 indices into splitCells. Note that the numbers in splitCells are not cell labels, they are purely indices into splitCells

[tomislav_maric]

I have some data that maps to the old mesh and doesn't get hit by the dynamicFvMesh::update(), so I tried to address the elements of the storage using the parentIndex(label), and setting values initially to like -100, to see what I am setting. As expected from the splitCells comments, I was either experiencing segfaults by faulty indices (greater than the size of the List) or setting some unrelated cell values all over the old data, definitely not the parent cells of the old, unrefined mesh.

But now I'm using the mapPolyMesh for this, since this information is actually used by the polyMesh::updateMesh(const mapPolyMesh&) to map the fields that are registered to the mesh. And this works great... The only trick I used is to inherit from dynamicRefineFvMesh, store an autoPtr<mapPolyMesh> and access it.

[kmooney]

Whelp, set up a bit bucket! I have a redistribution case set up for testing. I might be able to drop your code into mine and try to get ptScotch to distribute the way we want!

[kmooney]

...this is all based on the assumption that the refinementHistory::distribute(const mapDistributePolyMesh&) function is indeed working correctly. If its not working right we have a whole other can of worms to deal with.

[tomislav_maric]

It's a deal! I'll drop you an invite to the repository, and we can start from there.

Oh, my bet is that there'll be a truck full of worms to deal with.. I wold gladly deal with only just a can...

[tomislav_maric]

Hi,

it seems that the dynamicRefineFvMesh rewrites the connectivity info of a mesh when the unrefine method is called. This is a problem, since there is another layer of operations present between the old state, and the new state. Let's imagine we have multiple refinement steps that look like this:

U -> R0 -> R1 -> RN

where U marks the un-refined mesh and RN marks the refined mesh with a refinement level of N. If we store a mapPolyMesh for each refinement level, then getting back from RN to U is possible. However, sometimes, unrefinement is called:

U -> R0 -> R1 -> R2 U2 -> R3 -> ... -> RN (optionally UN)

where U is the unrefinement of cells. Unrefinement changes the mapPolyMesh of the refined mesh and messes up the connectivity. This is what I see when I try to address data from the step "R3" to "R2". There is no direct connection, since U2 basically involves a creation and a destruction of a mesh object (all the lists get expanded/shrinked, and the addressing of involved cells, faces points is taken over). We need to consider this as well....

T.