[mbfranklin@gmail.com]

Hello everybody,

I am working on implementing the conservative level set method presented by Olsson & Kreiss in their 2005 paper "A conservative level set method for two phase flow". I am currently using Scipy/Numpy for all of my numerical calculations.

I have been banging my head against the wall for months trying to debug my codes, and have yet to find any issues that are obviously incorrect, which brings me to this forum.

Has anybody else successfully implemented the level set method? I would love to have some conversation with anybody who is familiar with calculating normals, gradients, curvatures, etc. to bounce some ideas around, and perhaps share some numerical advice.

The first question that I have is this:

When calculating normal vectors from the gradients, one must divide by the magnitude of the gradient. In regions where the phase function is smooth, the gradients/magnitudes should go to zero, so dividing by this magnitude presents some numerical difficulty. Would you suggest

(1) adding a small number to the denominator (e.g. normal = gradient/(magnitude+eps) to avoid division by zero errors, or

(2) perform a check (e.g. if magnitude == 0: normal = 0)

More to come...any help would be MUCH appreciated!

[ata]

Hi
I think the first is better. Because with this approach you do not need an "if" and this can enhance speed of your code.
Enjoy
Good luck

Ata

[mbfranklin@gmail.com]

thanks for the reply! i'm new here, and really appreciate any help with my codes...looks like i came to the right place!

what value of eps would you suggest using in the calculation:

n = gradient/(magnitude + eps)

? i can set it to be any arbitrary value, and it has an affect on the normal vector field, including how far away it reaches from my circular interface. My initial phase-function varies smoothly from 1-to-0 over approximately 6 grid cells (very steep interface, but it's the value that Olsson & Kreiss report), and everywhere else it exponentially approaches zero, so at some point, the rate of exponential decrease towards the boundaries could possibly be on the order of eps. this can be seen in comparing two values of eps:

(1) eps = 1e-100 (extremely small number)

the normal field extends out completely to boundaries with magnitude =1.

(2) eps = 2.2e-16 (Python machine precision)

the normal field extends out about 2R from the circle of radius R, and the magnitude beyond that distance is =0.

so i don't like this approach, although it's what my advisor advised, too.

BUT using the "if statement" approach consistently extends the normal to the edge of the domain, and doesn't rely on some arbitrary constant to govern the normal vector field.

i'm so focused on my normal calculation, because it's my curvature field

kappa = divergence(normals)

which is giving me a real headache, although it should be simple!

more to come...

[ata]

Hi
Examine this two values:
Similar to OpenFOAM use
eps=1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
See OpenFOAM source code or similar to Henrik Rusche PhD Thesis "Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions" examine eps=10^-5
Enjoy
Good luck

Ata