Mesh-based Level Set Method Implementation

[NengLu]

See working cases (RTI model with free slip top) in https://github.com/NengLu/UW3_FreeSurf/tree/main/LevelSetMethod

TO DO:

• implementation in mesh deformed cases (free surface)
• test in parallel
• check function points_in_polygon
• reinitialization of level set function (Hillebrand et al., 2014)
   
  Reference:
• Hillebrand, B., Thieulot, C., Geenen, T., Van Den Berg, A. P., & Spakman, W. (2014). Using the level set method in geodynamical modeling of multi-material flows and Earth's free surface. Solid Earth, 5(2), 1087-1098.

[bknight1]

The results looks great @NengLu! Good work

[NengLu]

Thanks!@bknight1
I am trying sth with the free surface cases first.
It will be very helpful if you try the slab and detachment model (case d in Hillebrand2014), may need to be careful with the point in polygon function, I'm not sure if it works for the situation: the discontinue interface after the detachment happened.

[lmoresi]

If I were you @NengLu , I would try to use the kd-tree to set the initial value of the signed-distance function (phi) for the level set. In parallel, you need to be aware of that so-called fix that I checked in because 0.0 is a bad choice for phi in the absence of particles (i.e. it should be +/- 1).

[NengLu]

Thanks for the suggestion @lmoresi . Gona further to do some tests using kdtree instead scipy.distance.
About the phi, I considering it as [(phi<0, c1),(phi>=0, c2)] to include the condition phi=0. Here I calculated the phi by scipy.distance that have values (not just +- 1) as needed in the diffusion zone equation (Eq. (8) in (Hillebrand et al., 2014)).

[lmoresi]

Yes, that's the correct interpretation of the phi function so plus / minus infinity really but I know how much joy that brings ...

[lmoresi]

I was just surprised at how complicated the script became. The nice thing about level set methods is how simple they are when implemented. We can chat about it at the meeting tomorrow.

[NengLu]

Some updates (see details in https://github.com/NengLu/UW3_FreeSurf/tree/main/LevelSetMethod):

1. Using the kd-tree method instead scipy.distance in calculating the $\phi$ (sign-distance function)
2. I did to advect the $\phi$ is recacluating it based on the advected interface. (see MakerChainMethod.py)
3. Traditional way to advect the $\phi$ is using the equation (Eq. (4) in (Hillebrand et al., 2014)):
   $$\frac{\partial \phi}{\partial t}+\mathbf{v}\cdot \nabla \phi = 0$$
   When I tried to implement it with the uw3.advection solver, it had errors. I need help for the debugging @lmoresi .Thanks in advance! (see LevelSetMethod_fordebug.py)

[lmoresi]

The problem you are seeing is in this block. It is not related to the advection routine but it is triggered when we try to evaluate any expression. You have replaced the symbolic form of the mesh variable with something else (a piecewise function). What you need to do is evaluate the material function and assign this to the data array for the variable.

    material_fn = sympy.Piecewise((lightIndex,phi.sym[0] < 0.0,),(denseIndex, True))
    material_mesh.sym[0] = material_fn

@julesghub - we should make the .sym read only. It is not meant to be changed in this way.

[NengLu]

Thanks for pointing this out and it works after editing this material term (actually, I define this for plotting, it can be removed as it hasn't been used to define properties like density or viscosity).

The result looks pretty good. Here is the $\phi$ at the final step of the RTI model. (left is marker-chain, right is advection solver, and the colour back denotes the $\phi$ whose value is close to 0.).

But the advection solver is much slower than recalculating the $\phi$. Any tips to speed it up (already set the adv_phi.tolerance=1.0e-6)?

[lmoresi]

Don't worry about speed for the time being, worry about numerical stability etc. We can make things go faster once we are sure they work OK. You can always run in parallel !

For plotting, I suggest you change the scale to have zero in the centre (i.e. just let it be clim=[-A, A] where A is the max absolute distance - or some approximation to the largest distance in the simulation)