Why avoid reseting particle positions in DFNFlow?

Hi Robert,

I can tell that Timos set updatePositions=false in the scripts he sent me lately. I actually tried to run the same simulation (hydraulic injection as described in Timos's paper in IJRRMS) with and without the updatePositions and I could see that the updatePosition=True produced no fluid driven induced cracks... I don't remember the exact reason for defining this function (see the comment related to the triangulation just above the function) but I think that you need updatePosition to be false, otherwise the triangulation will not record the fractured cells.

Timos or Bruno could probably be more specific.

Luc

Hi,
What justified the existence of a specific function was the need to make it do nothing... There was probably a better way than inserting one line in a duplicated code but at least this one did not touch the code of another engine.

Why this option was needed is a bit less clear. The idea of not triangulating was initially to gain cpu time for stiff materials, as you suggest Robert, nothing else IIRC. We subsequently forgot to write code compatible with updatePosition=true, I guess (showing that premature optimizations are not a good idea in general).
Possible problems I suspect are related to the algorithm scanning the edges of the triangulation to check which ones have a crack, and changing permeability accordingly. Not sure it works very well if the triangulation itslef is changing.
I don't see a fundamental reason why it should not work, it seems more like a blindspot of implementation.

Bruno

Hi,

The idea was to match the initial interactions with the related edges in order to prevent crack-flagged facets from appearing and disappearing.

Which can also be a matter of stability especially in the case of impermeable matrix.

Now, for stiff materials it can also be a speed-up trick (similar to removing the new contact detection from the collider...).

For larger deformations we should probably find another way to relate crack-flagged facets to cracks. Then there will be no need for this function. That's all I think.

Cheers,
Timos

Thank you Luc, Bruno, and Timos for responding to my question.

Hmmmm, I see now. OK so we are using a Delaunay triangulation to build the facet orientations, right? Extrapolating Timo's comment, small deformations should not change this triangulation, and we should be safe to update positions. As soon as a large deformation occurs, the Delaunay triangulation might change and one of our tricked facets might disappear, leading to a change back to poiseuille controlled permeability and possible instability in the solution of the linear system.

If this is truly the problem, and our solution is to lock the positions of the particles (w.r.t PFV), I don't think the DFNFlow scheme is truly poroelastic since we aren't allowing the change of mechanical volume of the pore to contribute to the pressure in the pore.

Yes Timos, the solution seems to be that we need a new way to relate crack-flagged facets to cracks. Actually, this raises another concern I've had, which is that the current method cannot trick facets near a broken bond that does not line up with an edge on the triangulation. I wonder if we could trick facets by looking for their intersection with a broken bond, instead of looking for edges that are broken bonds. This would be more or less an entire rewrite of DFNFlow and the intersection algorithm might be expensive...and it may still not solve the instability problems discussed in this thread :-(.

Bruno, do you see a solution to this that can be implemented solely on Yade's side of things, without diving into CGAL's libraries? It seems to me we could "force" CGAL to use the Delaunay triangulation once, followed by a continuous use of the same vertex IDs for the same cells during subsequent triangulations. I doubt this sort of functionality is included in CGAL. If we need to dive into CGAL, do you see an easy solution within CGAL or is it a likely headache? I only ask because I know you have your head much deeper into the CGAL side of things than me. I'm willing to implement the solution, but I want to make sure I am headed down the right path before starting.

Best,

Robert

>I don't think the DFNFlow scheme is truly poroelastic

No, no, that's fine. We are only speaking of re-meshing here. The volume changes are still calculated as usual because it uses the actual positions/velocities from O.bodies at each iteration, not the coordinates of the triangulated points.
In fact what happens now is what you describe in your last paragraph (in #4, using vertices ids).
Which stability problem are we speaking about??

>the current method cannot trick facets near a broken bond that does not line up with an edge on the triangulation.

True. At the same time, it is very unlikely that two particles close to each other will not share an edge - almost impossible unless you define very large radius for remote interactions (>~2x diameter).

> I wonder if we could trick facets by looking for their intersection with a broken bond,

A facet is a set of three edges, if there is no edge then there is no facet... or maybe I missed your point.
Cheers
Bruno

>The volume changes are still calculated as usual because it uses the actual positions/velocities from O.bodies at each iteration, not the coordinates of the triangulated points.

Thank you, Bruno. Do you mind checking my understanding really quick? I am almost certain that stopping the setPositionsBuffer is ruining the compressibility scheme:

The reason I said the volume changes are not calculated was that we are computing cell->info().volume() here [1], which uses the positionsBuffer. Next, we prepare cell->info().invVoidVolume() using cell->info().volume() here [2]. Finally, we set up the A matrix, and we consider compressibility due to mechanical volume changes [3]. If we only set the positionBuffer once (as we do in DFNFlow), the positions [4] used to compute cell->info().volume() will not change, the void volume will remain the same and the compressibility correction will not change in the solve stage.

Actually, just writing this I think I figured out an easy solution. We just need to keep updating the positionsBuffer, but we just add an easy flag for avoiding retriangulation if DFNFlow is activated. Currently, DFNFlow setPositionsBuffer is not allowing us to continue setting positions. I think we just need to get rid of it, let the FlowEngine setPositionsBuffer do its work and add some flag that prohibits the geometrical retriangulations.

>Which stability problem are we speaking about??

Is it possible that instability may occur when the A matrix changes drastically between time steps due to appearing/disappearing trickedpermeability facets?

>A facet is a set of three edges, if there is no edge then there is no facet... or maybe I missed your

Of course, but surely we could consider an algorithm that looks for an intersection between an interaction and a facet. For example: the interaction pokes through the center of the facet (instead of lining up with one of the facet's edges).

Best,

Robert

[1]https://github.com/yade/trunk/blob/master/pkg/pfv/FlowEngine.ipp.in#L602
[2]https://github.com/yade/trunk/blob/master/pkg/pfv/FlowEngine.ipp.in#L439
[3]https://github.com/yade/trunk/blob/master/lib/triangulation/FlowBoundingSphereLinSolv.ipp#L218
[4]https://github.com/yade/trunk/blob/master/pkg/pfv/DFNFlow.cpp#L182

OK, my testing leads me to believe that the deformations are not large enough in a hydraulic fracture simulation to change the Delaunay triangulation. Therefore, it should be safe to update positions the normal way. I've found that the remesh interval needs to be much smaller than if we do not update the particle positions (1/5th-ish). So I guess simulation speed was the only real issue, but the slow down is not due to setting positions, but instead due to the more frequent remesh requirement (and frequent solver building associated with it).

In comparison, the hydraulic fracture behavior difference is pronounced for fixed particle positions vs updated positions. For updated positions, the rate of injection pressure drop is about half of the rate observed for fixed particle positions. Further, the crack area is greater for the fixed positions by 5-10% compared to the updated positions. So the fracture is propagating faster at lower injection pressures for the fixed positions, and the fracture is propagating slower at higher injection pressures with updated positions. This makes sense, considering we aren't considering the impact of the changing void volume on pressures for the fractured cells when we fix positions. When we update positions, these void volumes increase in fractured cells, the pressures decrease slightly due to fluid compressibility, and the pore pressure at the tip will not be as high, leading to lower fluid induced forces and a slower propagation of the fracture. For fixed positions, a fractured cells void volume remains the same for the duration of the simulation. So as soon as the fluid pressure increases in a fractured cell, the pressure raises more quickly and the fluid induced forces follow, leading to quicker fracturing, and the fluid moves more quickly to the next cell. For updated positions, the fracture is propagating slower, the pressure is not relieved as quickly by the fast permeability changes enabled by fractured cells. In other words, the slow process of opening the fracture is accommodated when we update positions, whereas fixed particle positions yield stiffer and less accommodating fracture propagation behavior.

*The rate of injection pressure drop following breakdown.

Hi Robert,
I think you are right on the problem, sorry for not realizing earlier. updatePositions=false really leads to a one way coupling.
It is very unexpected, since I think the DFN engine had been validated in primary consolidation at some point. Maybe at that time the volume computation was still using bodies[k].pos instead of buffer[k], that would explain. Or maybe it was tested with updatePositions=true.

>the remesh interval needs to be much smaller than if we do not update the particle positions

Why does it need to be smaller?

I think the best fix would be to update positions always and re-triangulate never (given that updating "always" while re-triangulating "sometimes" is already possible with updatePositions=true). That was the original intention to make it possible and it should not be difficult to achieve. As you suggested the "updatePosition" condition should apply for retriangulation, not for updating volumes. This way we would use fresh positions to calculate updated volumes of tetrahedra - regardless of whether those fresh positions would lead to the same set of tetrahedra in a new triangulation.

There is in fact no need to change the code to achieve "update positions always and re-triangulate never". It should be enough to set defTolerance<0, meshUpdateInterval=1e20 (*), updatePositions=True. Or does it give problems?

The other question on maintaining a consistent fracture pattern after re-triangulating is less clear to me. I'm not sure how it would look to have a facet "intersecting" a cracked edge without having it as one of its edges. If there is a solution to this it is probably a bit more involved, and probably not necessary in most cases.

Last but important note, see https://bugs.launchpad.net/yade/+bug/1734653 for another (and only) known bug besides what we are discussing here.

Cheers

Bruno

(*) For the sake of homogeneity I'll make meshUpdateInterval<0 discard the condition as it does for defTolerance.

>Why does it [the remesh interval] need to be smaller?

Trial and error show that there is indeed a stability threshold associated with the remesh interval. I suspect that the compressibility term in question (the one that is essentially ignored when we updatePositions=False), is increasing the maximum eigenvalues of our A matrix (since adding a term like this can only increase the distortion of the system). Higher eigenvalues means a lower required time step. Since we only update the A matrix every N mechanical timesteps (dt), we are essentially asking the A matrix to accept a total-perturbation-timestep of N*dt. Remember, this scheme allows us to cheat a little and get sub total-perturbation-timestep forces, which compound to alter volumes and change compressibility term for next A matrix update. As the remesh interval decreases, we decrease the total-perturbation-timestep that we ask of each A matrix (since it is being updated more frequently), and satisfy the higher eigenvalues of our compressibility-augmented system. Disclaimer, this is just an educated guess.

We could calculate the maximum total-perturbation-timestep based on the eigenvalues, but it would be expensive. It could be worth just checking the eigenvalues/vectors once or twice during a simulation depending on how dynamic it is.

I do see that Catalano looked at using a parameter based on pore volume change to control the retriangulation frequency (Eq. 5.13 of [1]). Unfortunately, 1% does not work well for hydraulic fracturing simulations since we want to update the permeability as soon as cells are tricked. We don't want to wait around for a volume change to reach a threshold. I suppose we could use a smaller threshold, but still, I'd rather just update as soon as cells perms are tricked.

I should also mention here that Catalano's max FlowEngine timestep based on particle density, diameter, and fluid viscosity (Eq. 4.20 of [1]) seems to overestimate the time step by a couple orders of magnitude for my hydraulic fracture stimulations. Or perhaps I didn't implement it right...

>There is in fact no need to change the code to achieve "update positions always and re-triangulate never". It should be enough to set defTolerance<0, meshUpdateInterval=1e20 (*), updatePositions=True. Or does it give problems?

This would not work in the current code, since the remeshing also goes through and checks for new cell perms to trick as bonds break.

> I'm not sure how it would look to have a facet "intersecting" a cracked edge without having it as one of its edges. If there is a solution to this it is probably a bit more involved, and probably not necessary in most cases.

This is possible more often when we use extended interaction ranges for particles, and less often when an interaction is in a location that does not agree with the Delaunay triangulation.

>Last but important note, see https://bugs.launchpad.net/yade/+bug/1734653 for another (and only) known bug besides what we are discussing here.

I spent a good amount of time testing different methods for keeping a special positionsBuffer for triangulation only, and allowing the other positionsBuffers to keep doing their duty for updating volume changes. I ran into some multithreading headaches for some methods, and just annoying unexplainable bugs for others. Anyways, the current scheme works if we just updatePositions and remesh at a shorter interval...albeit at the expense of the Delaunay retriangulation only (how much does Delauany triangulation even cost...doesn't seem like much).

[1]Catalano, E. (2012). A pore-scale coupled hydromechanical model for biphasic granular media. Application to granular sediment hydrodynamics, 200.

On 11/27/2017 07:32 PM, Robert Caulk wrote:
>
> Trial and error show that there is indeed a stability threshold
> associated with the remesh interval. I suspect that the compressibility
> term in question (the one that is essentially ignored when we
> updatePositions=False), is increasing the maximum eigenvalues of our A
> matrix (since adding a term like this can only increase the distortion
> of the system). Higher eigenvalues means a lower required time step.
updatePositions=False means a one-way coupling, which is indeed more
stable than updatePositions=True (two-way).
However the stability condition is in terms of timestep, not meshing
frequency. The origin of the dependency vs. remeshing is unclear to me.

> I do see that Catalano looked at using a parameter based on pore volume
> change to control the retriangulation frequency (Eq. 5.13 of [1]).
Yes, but that is just an accuracy requirement, not related to stability.
We never found instabilites by not remeshing.
For instance the oedometer.py consolidation will work like a charm if
you make it run from begin to end with the same mesh.

> I suppose we could use a smaller threshold, but still, I'd
> rather just update as soon as cells perms are tricked.
You are right, "as soon as a crack occurs" is the best option. A problem
is that in large problems it implies remeshing  virtually at every
timestep, which is prohibitively expensive.

> I should also mention here that Catalano's max FlowEngine timestep based
> on particle density, diameter, and fluid viscosity (Eq. 4.20 of [1])
> seems to overestimate the time step by a couple orders of magnitude for
> my hydraulic fracture stimulations. Or perhaps I didn't implement it
> right...
In the settings of Catalano I believe the equation is correct, and
Oedometer.py should confirm that.
I would have to understand your settings in more details.
Permeable/impermeable matrix? compressible/incompressible fluid? update
positions or not?
Are you really in a viscosity dominated regime (else the timestep is as
determined by GlobalStiffnessTimeStepper)?

>
>> There is in fact no need to change the code to achieve "update positions always and re-triangulate never". It should be enough to set defTolerance<0, meshUpdateInterval=1e20 (*), updatePositions=True. Or does it give problems?
> This would not work in the current code, since the remeshing also goes
> through and checks for new cell perms to trick as bonds break
I was speaking of _not_ remeshing. :)
Just update cells volumes. It should happen if "defTolerance<0,
meshUpdateInterval=1e20, updatePos=True".
Could you try that?

> I spent a good amount of time testing different methods for keeping a
> special positionsBuffer for triangulation only, and allowing the other
> positionsBuffers to keep doing their duty for updating volume changes.
It seems we have diverging visions of this question. From now on let's
assume that we always have  updatePos=True.
In this situation positionsBuffer is updated at every iteration and it
is used to determine the volume changes at each iteration.
Whether or not remeshing occurs is a completely different question, even
if it uses the same "bufferized" data. The idea of having two distinct
buffers is still another independent aspect related to task parallelism,
and actually it provides only a limited speedup.

> I ran into some multithreading headaches for some methods, and just
> annoying unexplainable bugs for others. Anyways, the current scheme
> works if we just updatePositions and remesh at a shorter
> interval...albeit at the expense of the Delaunay retriangulation only
> (how much does Delauany triangulation even cost...doesn't seem like
> much).
Triangulation is fast compared to Cholesky factorization (CGAL also made
the triangulation parallel recently, which I did not try to use yet). It
means that tricking just one permeability without changing the mesh
(thus implying a full new Cholesky decomposition - unless advanced
methods enabling factors updates are applied) costs nearly the same as
remeshing (which also implies one Cholesky pass).

Could you show a basic example script with your typical settings?

Cheers
Bruno

Dear Bruno,

Thank you so much for taking the time to respond to this thread. Luckily, my original question has been more than answered and we can all now agree that updatePositions needs to be set to True if we want a true two way coupling. So that is good. But the remesh frequency instability is odd, and may require some more thought on my side. Anyways, after going through and trying to reply to your replies to my replies (LOL), I started to feel that these discussions are becoming too unwieldy for launchpad :-(. I guess we could try to set up a skype call or something if we felt the need to discuss this in depth. But here are my answers anyways.

>updatePositions=False means a one-way coupling, which is indeed more
>stable than updatePositions=True (two-way).
>However the stability condition is in terms of timestep, not meshing
>frequency. The origin of the dependency vs. remeshing is unclear to me.

Hmmm, so the idea that we are asking a single A matrix to accept a total perturbation that amounts to N*dt is not a problem? Well, perhaps it is an issue with the permeability tricking? If we increase pressure in our fracture, and many bonds break but the permeability is not tricked (because our remesh step is high), pore pressures will continue to increase as if the fracture has not propagated. Then, as soon as we do trick the permeability, we get a pressure gradient that is so high it requires a smaller time step to maintain stability. Maybe I am grabbing at straws now. I am not sure how else to explain the instability.

>You are right, "as soon as a crack occurs" is the best option. A problem
>is that in large problems it implies remeshing virtually at every
>timestep, which is prohibitively expensive.

I guess it depends entirely on the simulation. I've had some success with the method, it is currently an option in yade with flow.breakControlledRemesh=True. But if I have a few thousand broken bonds for a simulation that lasts 30000 steps, I would be refactorizing every 10 steps. That would be prohibitively expensive. I've thought about adding some logic to the current breakControlledRemesh so that it waits for 3-10 broken bonds before remeshing instead of 1.

>I was speaking of _not_ remeshing. :)
>Just update cells volumes. It should happen if "defTolerance<0,
>meshUpdateInterval=1e20, updatePos=True".
>Could you try that?

Oh I see, you are suggesting I try this as a way of narrowing down the source of the instability? I suppose I could, but I would need to turn off the injection flow rate, or else it would blow up. Without injection on, the model will be unconditionally stable. Maybe I do not understand.

>In this situation positionsBuffer is updated at every iteration and it
>is used to determine the volume changes at each iteration.
>Whether or not remeshing occurs is a completely different question, even
>if it uses the same "bufferized" data. The idea of having two distinct
>buffers is still another independent aspect related to task parallelism,
>and actually it provides only a limited speedup.

I was just talking about using a separately stored buffer for triangulation only, and then everything else (volume changes and total volumes) is computed using the updatedbuffer.

If I am reading the code correctly, in the situation where updatePostions=True, the RHS reflects volume changes at every timestep [1] but the total mechanical volume updates are only reflected in the A matrix every N flow.meshUpdateInterval steps [2]. This is where I thought the remesh frequency instability was originating when I was talking about the total perturbation timestep for the A matrix, but it seems that may not be the case.

Cheers,

Robert

[1]https://github.com/yade/trunk/blob/master/lib/triangulation/FlowBoundingSphereLinSolv.ipp#L328
[2]https://github.com/yade/trunk/blob/master/lib/triangulation/FlowBoundingSphereLinSolv.ipp#L218

Hello,

I am going to conduct yade simulations of fracture-induced volumetric expansion coupled with fluid flow.
The simulation will involve significant particle volume increase up to 50%, and one important question would be the propagation of fluid within the fractured medium.
In such simulation, the retriangulation with new particle positions together with DFNflow would be necessary.

Before going details into this issue on my side, I would like to ask kindly do you have any updates of ideas/codes in DFNflow concerning the retriangulation, e.g. is it necessary to relate crack-flagged facets to cracks?
Or maybe you can give me some suggestions about where to start, considering that I begin to learn to develop the source code.

Thank you very much in advance!
Best regards.
Lingran

Hello Lingran,

Happy to hear we have someone else to work on DFNFlow. Welcome :-)

In my opinion, you can do this with the current DFNFlow just fine. You simply need to set updatePositions=True and use a high remesh frequency (low meshUpdateInterval [1]). Warning: this will be computationally expensive.

As far as changing the particle sizes, you will need to do this with your constitutive law (JCFpm), not with DFNFlow.

>is it necessary to relate crack-flagged facets to cracks?

You mean: "do we need a NEW way to relate crack-flagged facets to cracks". The answer is, yes, only if your deformation is so large that the Delaunay triangulation is changing. In my experience with hydraulic fracture rock simulation, the Delaunay triangulation never changed. Your application sounds a bit theoretical. Maybe you can provide more information? What kind of stiff brittle (crack-able) material also experiences 50% solid volume changes?

Best,

Robert

[1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.FlowEngine.meshUpdateInterval

Hello Robert,

Thank you for your rapid feedback.

The simulation I talked about is related to rock hydration reaction, such as serpentinization, a geochemical metamorphic reaction that transfers Olivine to Serpentinite with the presense of water. The reaction is associated with significant changes in rock densities, volumes, mechanical properties, etc. Since the rock has low permeability and the volume expansion during the reaction tends to clog the initial pores, the progress of the reaction relies on the reaction-induced fractures which allow water to penetrate and contact with the unreacted materials. We have simulated this kind of reactions without the coupling with the fluid.

For the next step, we would like to couple the volumetric expansion with the fluid flow. With such large volumetric expansion of soild particles, the Delaunay triangulation will change. The retriangulation might give rise to some questions about cracks and permeability updates, as well as instability problems in the current DFNflow.

So my question is: how to propose a new way to relate crack-flagged facets to crack?

Lingran

Hi,
Just to summarize a few ideas after a discussion with Robert:
- the duplicated "setPositionsBuffer()" in DFN was definitely an experimental hack and it is unsafe to keep relying on this.
- removing this duplicate would let two-way-coupling work as it should (since the inherited setPositionsBuffer is not conditional).

What is really needed is 1/ a way to update permeability data based on cracks extension and aperture and, 2/ to rebuild the permeability matrix after that. To be sure no remeshing will occur it is enough to set FlowEngine's meshUpdateInterval and defTolerance negative.

Further: if permeability is only updated in the cracks (not in intact matrix) some performance could be gained by using cholmod update/downdate features (cholmod_updown), which updates factorization without going through the complete Cholesky algorithm (even changing the size of the problem for advancing cracks seems to be possible).

Bruno