[sbaffini]

I'm all but an expert in this field but, as a practictioner in need of a positivity preserving scheme, I gathered the following information:

1) A 1D, first order positively preserving scheme can be made nD and high order (EDIT: I now see this is similar to a work you already cited, so you probably know this as well)

https://deepblue.lib.umich.edu/bitst...pdf?sequence=1

2) The scheme that Einfeldt defines as HLLEM in "On Godunov-type methods near low densities" is proven by himself to be positively preserving and has been later used as a basis for a positivity preserving entropy fix for a Roe scheme, e.g.,

https://perso.ensta-paris.fr/~pelant...nti_oxford.pdf

Altough, I see there might be misalignment in nomenclature across references

3) My personal experience is with a preconditioned Roe scheme using an entropy fix based on the preconditioned HLLE+ scheme

https://arc.aiaa.org/doi/10.2514/1.12176

and it is positivity preserving at 2nd order for a shock tube on an unstructured grid (that is, an actual 3D tube with a finite cross section). Of course, a working limiter is needed here.

May I suggest you to use the very test in the Einfeldt paper to check your scheme?

[dappoli]

Originally Posted by sbaffini

Altough, I see there might be misalignment in nomenclature across references

I agree, for instance looking at the paper you linked for the HLLE+, it seems that I have been using the HLLE+ all the time, where while the HLLEM scheme does not use the Roe average for the contact discontinuity velocity.

However in the paper I followed originally:

https://www.sciencedirect.com/scienc...092?via%3Dihub (page 3) it uses the Roe averages.



Anyway I will try Einfeldt test to see what happens

[Eifoehn4]

Originally Posted by dappoli

First of all thank you for the reply.

I think we are talking about the same thing, for me a positive preserving scheme is a scheme that generates "physically admissible solutions from physical data" [1], or in other words it keeps positive pressure and density if the initial condition are positive. For the Euler equations the physically admissible state in this case is defined as:


being U the conserved variable vector and m the momentum along the 3 directions. For example Roe scheme is not positive preserving and neither is Osher's. Einfeldt [2] demonstrated that the Rusanov and HLLE are positive preserving, but for example HLLEM is not. All this holds true for 1st order schemes, once 2nd order schemes are taken into account, from what I found there is no proof of satisfying the positive preserving property. That is why the Muscl reconstruction near low density fails also for Rusanov. Clearly if the timestep is strongly reduced this is in general postponed, but inevitable. Having said this, I am new to the study and application of Riemann solvers so I might have said something wrong.



Going back to Noh:

You are right, but being focused on other problems, I forgot that I am using a variation of the MINMOD slope limiter described in Kurganov-Tadmor [3] (the scheme I am using which is equal to a MUSCL Rusanov) defined as:


which allows to reduce the diffusivity of the minmod by tuning the parameter . Setting the parameter to 1 recovers the original minmod. I had it at 1.4. Attached you can find the solution with the parameter = 1. Now the solution is sharp, although the HLLEM keeps behaving strangely at the centre of the domain (note that the simulation is the same, I just cut the plot). Do you have any idea of why it behaves like that at the origin?



Finally, I coded the "hurricane" test described in [4] where a uniform domain is initialized with a vortex like velocity. This creates a vacuum point in the centre of the domain and if the scheme is not positive preserving, eventually it will fail. Of course 2nd order MUSCL Rusanov, HLLE and HLLEM fail in this, but first order HLLE and Rusanov do not. What I did in the end is to use a posteriori approach:

I advance 1 timestep,

I check for negative pressure, temperature or NANs ,

If I detect errors, I reject the step and re do it with 1st order Rusanov in the troubled cell and its neighbours.



I already tried doing this with little success, but I was not considering the neighbouring cells and just using the Rusanov flux in the troubled one. This eventually does not work. Now I will have to deal with some special cases such as when the troubled cell is on the boundary between processors (I use domain decomposition with MPI).

Attached you can see a video of the hurricane like simulation with the HLLEM a polsteriori corrected scheme. Sorry for the low quality, the solution is smooth in reality but the forum has a very low size limit



Now I will have to try this on my problem, probably with a careful optimization of the minmod limiter and neighbour radius to avoid failure.


[1] https://link.springer.com/chapter/10...-011-5169-6_16
[2] https://www.sciencedirect.com/scienc...21999191902113
[3] https://www.math.umd.edu/~tadmor/pub...or.JCP-00I.pdf
[4] https://www.math.hkust.edu.hk/~makxu...D-inviscid.pdf

Yes I think so. Now it makes sense, the simulation looks similar now. The HLLEM (HLLE+) for the Noh test might be a normal behavior.

I have also testest the Hurrican example, the high speed example, with reconstruction (TVD, MINMOD) + Rusanov and even HLLC and had no problems. Your video seems unsymmetric, the solution should be symmetric.

[arjun]

May i ask, how much is lowest density you are running into?

[Eifoehn4]

Originally Posted by arjun

May i ask, how much is lowest density you are running into?

Depends on the mesh resolution and the settings. With linear TVD reconstruction (MINMOD), CFL=0.5, explicit Euler O(2) and HLLC I get following values for the hurricane example (high speed):

[arjun]

Originally Posted by Eifoehn4

Depends on the mesh resolution and the settings. With linear TVD reconstruction (MINMOD), CFL=0.5, explicit Euler O(2) and HLLC I get following values for the hurricane example (high speed):

Thank you. I was just curious. I actually have worked on cases with very low density and I decided against Riemann type solvers. I decided to go pressure based solver route. With that in my case i could go as low as 1E-5 without any problems, and i can go till 1E-7 if i push it.

[Eifoehn4]

Originally Posted by arjun

Thank you. I was just curious. I actually have worked on cases with very low density and I decided against Riemann type solvers. I decided to go pressure based solver route. With that in my case i could go as low as 1E-5 without any problems, and i can go till 1E-7 if i push it.

Sounds great.

We should keep in mind, that the example is designed in such a way, that it contains a physical singularity in the center (vaccum point). Every consistent, convergent and positivy preserving solver should reach a density of without much effort (under mesh refinement) until round-off errors come into play.

The validity of the Euler Equations under very high Knudsen numbers is a different topic.

[arjun]

Originally Posted by Eifoehn4

Sounds great.

We should keep in mind, that the example is designed in such a way, that it contains a physical singularity in the center (vaccum point). Every consistent, convergent and positivy preserving solver should reach a density of without much effort (under mesh refinement) until round-off errors come into play.

The validity of the Euler Equations under very high Knudsen numbers is a different topic.

Can you point out which test you are talking about. Sorry a lot going on in this thread so for me hard to keep track.



I can try with pressure based solver and see what happens. I expect it to break down too but things can't improve if they don't break.

[Eifoehn4]

Originally Posted by arjun

Can you point out which test you are talking about. Sorry a lot going on in this thread so for me hard to keep track.



I can try with pressure based solver and see what happens. I expect it to break down too but things can't improve if they don't break.

The numbers are given for the hurricane example, reference [4], see the previous answers.

[FMDenaro]

Originally Posted by dappoli

These are the equations:




which are integrated for electrons and ions where q is the charge ( e, the fundamental charge), B the magnetic field, in my case perpendicular to the domain,
the total energy ( Note that T is expressed in eV here), the energy injection source term and the elctric potential obtained from:






At the moment I am trying also with a Strang splitting, so that I solve for the homogeneous system which is identical to Euler. The eigenvalues are the usual eigenvalues. In the algorithm I advance one dt the single species and then compute the new potential for the next timestep.

If you are not familiar with plasma, just as a reference, the problem is complicated by the broad range of timescales. For instance, for Xenon m_i = 2.18e-25 kg and m_e = 9.1e-31 kg, so that there is a huge difference in the timescales of ions and electrons. Thermal velocities for electrons are in the order of 1e6 m/s and the temperature is 1-2 order of magnitude larger than the one of ions. Moreover, to solve the non-neutrality of the problem, the debye length and plasma frequency must be resolved, limiting the problem to very fine meshes and timesteps in the 0.5-10 ps order.



I know that in the pressure-density-velocity the TVD is not granted, but more in general the positivity of schemes such as the Rusanov is not granted either for multidimensional system with reconstruction > order 1. The negative temperature sometimes occurs in the face reconstruction, but that is easily solvable (or at least more easily with a dedicated limiter), but sometimes it just pops up in the cell centre. After I advance in time the conservative variables, recovering the temperature in the cell centre my give rise to a negative value. Note that in this problem it is not necessary to know the temperature at cell centre, but if one has to compute collision frequencies ecc. it is mandatory.



I know that there is a work from Linde and Roe where multidimensional schemes of order >1 are made positive preserving using a dedicated limiter, which however is somehow not straightforward for the pressure. Indeed they state that the pressure must be limited iteratively to guarantee the positivity. I would try to avoid this given the fact that the timestep is already incredibly small and time can be an issue.



https://link.springer.com/chapter/10...-011-5169-6_16


By the way, does it make sense for you that i can see a discontinuity in the solution (again in the temperature, but I guess only because I am dividing by the discontinuous density) where I use 1st order rusanov instead of 2nd order HLLEM? Should the interface be smooth? Because I might have an error in the code if that is the case

I am not an expert with plasma physics and I never worked on your problem but I would ask you how do you get the Poisson equations for phi from the second equation, what are the BC.s you prescribed and how do you solve.

If I understand correctly, the solved function is a "source" term in the energy (third equation).


What is more, did you evaluate if you have a stiff hyperbolic problem due to the coupling?

[Eifoehn4]

Originally Posted by FMDenaro

I am not an expert with plasma physics and I never worked on your problem but I would ask you how do you get the Poisson equations for phi from the second equation, what are the BC.s you prescribed and how do you solve.

If I understand correctly, the solved function is a "source" term in the energy (third equation).


What is more, did you evaluate if you have a stiff hyperbolic problem due to the coupling?

The model originally stems from the coupled Euler-Maxwell equations, in the limit of electro-statics. Therefore, the potential is the electro-static potenial. The gradients (electric field) exerts forces (Lorentz force) on the electron or ion species, which carry charge, which is some how similar to pressure for the incompressible Navier-Stokes equations. The model is well-known in the semi-conductor industrie.

[FMDenaro]

Originally Posted by Eifoehn4

The model originally stems from the coupled Euler-Maxwell equations, in the limit of electro-statics. Therefore, the potential is the electro-static potenial. The gradients (electric field) exerts forces (Lorentz force) on the electron or ion species, which carry charge, which is some how similar to pressure for the incompressible Navier-Stokes equations. The model is well-known in the semi-conductor industrie.

Ok, thanks, that confirms the reason of my question.
From application of the divergence operator on the second equation ?
What about the BCs and the discretization of the Poisson problem?

[Eifoehn4]

Originally Posted by FMDenaro

Ok, thanks, that confirms the reason of my question.
From application of the divergence operator on the second equation ?
What about the BCs and the discretization of the Poisson problem?

The potential rests on the assumption that the electric field is irrotational:



The boundary conditions and discretization may be choosen independantly.

[FMDenaro]

Originally Posted by Eifoehn4

The potential rests on the assumption that the electric field is irrotational:



The boundary conditions and discretization may be choosen independantly.

Could you please better detail? I mean, apply the divergence operator to the second equation and write down the Neumann BCs.
Is that sufficient to ensure the compatibility condition for the existance of a solution?

[Eifoehn4]

Originally Posted by FMDenaro

Could you please better detail? I mean, apply the divergence operator to the second equation and write down the Neumann BCs.
Is that sufficient to ensure the compatibility condition for the existance of a solution?

To provide a little clarity: The model here is an extention of the standard Euler equation with additional external forces (source terms). Here, the additional external force is the electric/(magnetic) field, either supplied by the full Maxwell equations or by the Coloumb law via the electro-static potential.

You may solve the compressible Euler-Maxwell or Euler-Poisson equation in a pressure based fashion rather than a density based fashion. In this case you would end up with an additional Poisson like equation, however for pressure, including additional source terms. Except external forces, there should be no difference to handle BCs in order to fullfill the compatibility condition for pressure.

The former electro-static Poisson equation should of course also satisfy the compatibilty condition in order to be unique.

[dappoli]

Hi everyone, sorry for not replying for a while.

Eifoehn4 is right, the system comes from the Euler-Maxwell in the electrostatic limit. In my case boundary conditions for Poisson are simple dirichlet.

Originally Posted by arjun

May i ask, how much is lowest density you are running into?

In my specific problem I got number densities as low as 0.7 m^-3 while other regions of the domain are at 10^15, so very low. It turns out that even with a more diffusive scheme sometimes the simulation fails, in particular when I have a sudden diffusion of density that expands in the low density region.

This creates a non physical shock with extreme temperatures (to reproduce the problem try to do a free expansion in quasi-vacuum). In order to avoid this I decided to try another approach and coded a "vacuum tracking" algorithm which follows the interface between absolute vacuum (just 0 ) and solves Riemann problems with vacuum. This works fine especially in the free expansion case since the approximate/exact riemann solution is available (see Munz 1994), however once I have source terms, such as magnetic confinement ecc, it is not anymore free expansion and I think I have to modify the riemann solution to account for this. In particular I guess I have to remove the rarefaction wave since nothing can expand beyond the magnetic line (free expansion always assumes a rarefaction wave at the interface, being the contact discontinuity impossible to obtain due to the pressure gradient). I am not really sure on how to modify the flux function to account for this but I will look into it.



I never thought of using pressure based solver, I honestly know only the one like Piso or Simple, but I always skip them for this kind of problems. Might be due to the fact that I studied on LeVeque and he just goes with Riemann problems. I will look into them as well, thanks

[arjun]

Please check the Toro's book as that is the only book that kind of talk about Riemann solvers with Vacuum.

I thought about the appraoch you talk about detecting Vacuum and this is why i was reading the Toro's book for it.

The major reason why I gave up on Riemann solver was that i had a calculation with Vacuum pump (2 phase flow) and Riemann solver did not give me good flow profile in rotating geometry while pressure based solver did give me good profile. In this pump low pressures are generated so was exploring the idea of detecting vacuum and then going forward with modifying Riemann solver.

The unreliability of the Riemann solver made me give up on it.



I point out that the solver was validated for standard test problem so it was not inaccurate just that when rotation was involved things did not add up.

I actually wrote the Riemann solver in long double precision because of low densities involved.


Anyway keep us updated and also check the Toro's book. I do work on same problem actually and will be spending next 2 to 3 weeks on it (i resume where i left off). I have electron transport and now i will add ion transport with it.

Originally Posted by dappoli

Hi everyone, sorry for not replying for a while.

Eifoehn4 is right, the system comes from the Euler-Maxwell in the electrostatic limit. In my case boundary conditions for Poisson are simple dirichlet.





In my specific problem I got number densities as low as 0.7 m^-3 while other regions of the domain are at 10^15, so very low. It turns out that even with a more diffusive scheme sometimes the simulation fails, in particular when I have a sudden diffusion of density that expands in the low density region.

This creates a non physical shock with extreme temperatures (to reproduce the problem try to do a free expansion in quasi-vacuum). In order to avoid this I decided to try another approach and coded a "vacuum tracking" algorithm which follows the interface between absolute vacuum (just 0 ) and solves Riemann problems with vacuum. This works fine especially in the free expansion case since the approximate/exact riemann solution is available (see Munz 1994), however once I have source terms, such as magnetic confinement ecc, it is not anymore free expansion and I think I have to modify the riemann solution to account for this. In particular I guess I have to remove the rarefaction wave since nothing can expand beyond the magnetic line (free expansion always assumes a rarefaction wave at the interface, being the contact discontinuity impossible to obtain due to the pressure gradient). I am not really sure on how to modify the flux function to account for this but I will look into it.



I never thought of using pressure based solver, I honestly know only the one like Piso or Simple, but I always skip them for this kind of problems. Might be due to the fact that I studied on LeVeque and he just goes with Riemann problems. I will look into them as well, thanks

[Eifoehn4]

Originally Posted by arjun

Please check the Toro's book as that is the only book that kind of talk about Riemann solvers with Vacuum.

I thought about the appraoch you talk about detecting Vacuum and this is why i was reading the Toro's book for it.

The major reason why I gave up on Riemann solver was that i had a calculation with Vacuum pump (2 phase flow) and Riemann solver did not give me good flow profile in rotating geometry while pressure based solver did give me good profile. In this pump low pressures are generated so was exploring the idea of detecting vacuum and then going forward with modifying Riemann solver.

The unreliability of the Riemann solver made me give up on it.



I point out that the solver was validated for standard test problem so it was not inaccurate just that when rotation was involved things did not add up.

I actually wrote the Riemann solver in long double precision because of low densities involved.


Anyway keep us updated and also check the Toro's book. I do work on same problem actually and will be spending next 2 to 3 weeks on it (i resume where i left off). I have electron transport and now i will add ion transport with it.

Did you consider All-Mach Riemann solvers? Without any fix, Riemann solvers become quite diffusive in the low Mach regime. The reason for this is that the jump terms (responsible for numerical diffusion) depend on all characteristics, also the highest, e.g. the speed of sound.

[arjun]

Originally Posted by Eifoehn4

Did you consider All-Mach Riemann solvers? Without any fix, Riemann solvers become quite diffusive in the low Mach regime. The reason for this is that the jump terms (responsible for numerical diffusion) depend on all characteristics, also the highest, e.g. the speed of sound.

Thank you. I will check it once i have some time to spare. I used hllc in that sim.