Level set convection - Zalesak sphere - problems with new algebraic stabilization element

[ddiezrod]

Description I've created a 3D zalesak sphere test to check how well the different level set options work. For this test, I am doing half a turn with 200 time steps and checking how well the original shape is maintained.

From what I see. the original supg levelset is behaving much bettter than the new agebraic stabilization element.

Left: SUPG, Right: New algebraic stabilization Original shape: Quarter turn Half turn

In both cases, the test is run with bfecc. We can clearly see that shape is much better preserved in the SUPG element and that we even have some oscillations/undershoots in the "air" region in the new element.

I have two questions regarding this:

• Is there any time step limitation for the new element? This could explain why its behaving much worse in this case.
• In the original paper where the formulation of the new element is described, there is a "w" parameter that seems to have an important effect on the results. I did not see this parameter in the implemented formulation, was this parameter maybe somehow simplified?

@rubenzorrilla @mrhashemi

[mrhashemi]

Please try using a much smaller CFL number (start with 0.1 for example) for the algebraic stabilization. It uses an explicit formulation for higher-order (generally nonlinear) terms and therefore, the time-step should be restricted.

[ddiezrod]

Please try using a much smaller CFL number (start with 0.1 for example) for the algebraic stabilization. It uses an explicit formulation for higher-order (generally nonlinear) terms and therefore, the time-step should be restricted.

OK, I'll try that. Does this mean that in general CFL should be lower than 1? In our problems CFL is usually 10, so if thats the case I guess we just cannot use this formulation.

[mrhashemi]

Please try using a much smaller CFL number (start with 0.1 for example) for the algebraic stabilization. It uses an explicit formulation for higher-order (generally nonlinear) terms and therefore, the time-step should be restricted.

OK, I'll try that. Does this mean that in general CFL should be lower than 1? In our problems CFL is usually 10, so if thats the case I guess we just cannot use this formulation.

You can, by using internal substeps. dividing for example each step to 20 or 30 substeps (only for the convection). Otherwise, you should implement some nonlinear formulation and pay the price (of the computational cost) by doing Newton-Raphson or fixed-point iterations.

[ddiezrod]

Ok @mrhashemi , thanks for your answer. I will upload the results using CFL =1 in both cases.

[mrhashemi]

Ok @mrhashemi , thanks for your answer. I will upload the results using CFL =1 in both cases.

• Please first try CFL=0.1 for the case with HO term ("include_anti_diffusivity_terms" : true) and then, increase it (anyway it should be smaller than 1 for this explicit formulation).

• Use CFL~1 for algebraic stabilization without HO terms ("include_anti_diffusivity_terms" : false).

• For the std. cross-wind stabilized element, using Crank-Nicolson, there is no limitation on CFL.

[ddiezrod]

@mrhashemi I ran the tests with cfl = 0.95 (using substepping) and the results do not look much better, although instabilities have disappeared.

left (with antidiffusive) right (without)

I could not run with cfl = 0.1 (with cfl 1 this is taking already something like 10 h)

[mrhashemi]

With the current explicit anti-diffusion terms, it is impossible to use CFL~1.

Could you please let me know how you calculate h for determining CFL? I assume you let the code do that and if I am not wrong, the code uses an average element size (1/|grad N|). Since the elements are not of equal length at all the edges, setting CFl=0.95 is a little bit tricky.

Let me try it myself and let you know.

[ddiezrod]

With the current explicit anti-diffusion terms, it is impossible to use CFL~1.

Could you please let me know how you calculate h for determining CFL? I assume you let the code do that and if I am not wrong, the code uses an average element size (1/|grad N|). Since the elements are not of equal length at all the edges, setting CFl=0.95 is a little bit tricky.

Let me try it myself and let you know.

Yes, Im just letting the code compute its own substeps so CFL = 0.95 is used.

[mrhashemi]

@ddiezrod thank you for raising this issue. I did some further tests: despite the promising result of the currently implemented algebraic stabilization with anti-diffusive terms in 2D (using a rather fine mesh), it is still too diffusive in 3D using a very coarse mesh. The solution is to substitute the average distance gradient with the exact value and make the limiter sharper (increase the power). However, I have not yet completely tested the side effects and cannot implement these modifications. I should also investigate the time marching approach. Anyway, the current implementation would work fine in 3D if the mesh is fine enough and CFL is small.

Just for comparison, these results are obtained without BFECC (CFL is very small and the mesh is coarse): 1- Algebraic stabilization with anti-diffusion terms:

2- SUPG + CW:

Initial:

[ddiezrod]

HI @mrhashemi , thanks for your comment. I wouldn't say the mesh is very coarse in my case (about 2.7 M elements) but it still looks very diffusive even using cfl = 0.5. I'm not very familiar with the formulation but I can try to test something if you want. Also, if you are interested. I can send you the test by mail, it seems its too big to upload it here.

[mrhashemi]

Hi @ddiezrod The mesh you use is quite fine :+1: I think you can reduce the number of elements to ~500K. The CFL I used is in the order of 0.1 (counting on h = the shortest edge of the elements).

If you agree, I can create a draft PR and add a couple of additional setting parameters. Then, it would be easier to run more tests.

[ddiezrod]

Hi @ddiezrod The mesh you use is quite fine 👍 I think you can reduce the number of elements to ~500K. The CFL I used is in the order of 0.1 (counting on h = the shortest edge of the elements).

If you agree, I can create a draft PR and add a couple of additional setting parameters. Then, it would be easier to run more tests.

That'd be great, thanks!

[RiccardoRossi]

news on this? i think it is quite important!

[mrhashemi]

news on this? i think it is quite important!

I'll first make some main settings optional in the convection process (the most important one is the acuteness of the limiter separately for BFECC and for the more accurate algebraically stabilized scheme). This would not significantly improve the CFL limitation. The next step would be to update the gradient (using a first-order scheme, maybe the semi-Lagrangian method).

[ddiezrod]

I ran again this test with the latest changes, it looks much better (although I think I also fixed a test in my code that could be important). Both cases were run with tetha = 0.5, no BFECC and CFL = 0.8. Results with SUPG still look superior to those with algebraic stabilization. (Left algebraic, right SUPG after a half turn)

[mrhashemi]

I ran again this test with the latest changes, it looks much better (although I think I also fixed a test in my code that could be important). Both cases were run with tetha = 0.5, no BFECC and CFL = 0.8. Results with SUPG still look superior to those with algebraic stabilization. (Left algebraic, right SUPG after a half turn)

Is the SUPG result obtained with crosswind stabilization (0.7)? Compared to the first test that you performed with BFECC (posted above), it seems that BFECC worsens the SUPG+CW, which is quite surprising since, in 2D tests, it gives highly diffusive results. Anyway, I do not think there is any further room to improve what you've obtained here for this 3D case. What I can suggest is to run and compare the results for the standard 2D Zalesak test (with three shapes).

[ddiezrod]

I used CW=0.0. I did a couple of modifications in the test so to compare with bfecc I should run the test again, I can run them again using bfecc if that info is interesting to you. I could try also the 2D, but we are only interested in 3D resulst at the moment...

[mrhashemi]

I used CW=0.0. I did a couple of modifications in the test so to compare with bfecc I should run the test again, I can run them again using bfecc if that info is interesting to you. I could try also the 2D, but we are only interested in 3D resulst at the moment...

Ok, now, it makes sense. Without CW, the method is not stable. There is a high chance of instabilities with pure SUPG. So, it is better to set it to the default value even though the result for this test seems fine. In case you are focused on 3D cases, maybe I can find a benchmark.

[ddiezrod]

I ran again the test using the default value (CW = 0.7) and result is now much better for algebraic stab. (left algebraic, right supg )

I have to say that we normally use cw = 0.0 and I cannot say results are always great, but we do not see instabilities very often., even with high CFLs

[mrhashemi]

I ran again the test using the default value (CW = 0.7) and result is now much better for algebraic stab. (left algebraic, right supg )

I have to say that we normally use cw = 0.0 and I cannot say results are always great, but we do not see instabilities very often., even with high CFLs

Maybe I used an incorrect term. With theta=0.5 and pure SUPG you can choose very large CFLs. The problem is with not complying with DMP (it is not monotonic) and the result would be noisy with over/undershoots once a steep gradient occurs.

[ddiezrod]

Maybe I used an incorrect term. With theta=0.5 and pure SUPG you can choose very large CFLs. The problem is with not complying with DMP (it is not monotonic) and the result would be noisy with over/undershoots once a steep gradient occurs.

Got it! Thanks for the all the explanations! Is there any plan for making algebraic stab. work with higher CFLs or you are just interested in low CFLs for the moment?

[mrhashemi]

Maybe I used an incorrect term. With theta=0.5 and pure SUPG you can choose very large CFLs. The problem is with not complying with DMP (it is not monotonic) and the result would be noisy with over/undershoots once a steep gradient occurs.

Got it! Thanks for the all the explanations! Is there any plan for making algebraic stab. work with higher CFLs or you are just interested in low CFLs for the moment?

I have not tested the maximum possible CFL. Maybe it works up to CFL~1.5 (or even larger) with theta=0.5. I'd appreciate it if you can test it for this 3D case. There are other approaches (mostly somehow implicit) that might work for large CFLs. At the moment, I have no plan for doing further implementations. The only modification that I will implement (in soon future hopefully) is a modified BFECC algorithm. It is slightly more costly than the current algorithm, but further improves the results.

[ddiezrod]

I have tested and at least until CFL = 5 it is working fine. With CFL = 10 you start seeing instabilities. I will close this issue for the moment, I think we can consider that algebraic stabilization works fine for low values of CFL, but its not recommended to use in general free surfaces cases due to this limitation. Also, although it might be more unstable in some cases, the supg element with cwd = 0.0 is less diffusive.