Implement CompressibleEulerPlasma equations

[warisa-r]

The system consists of two conservation laws, one for ions and the other for electrons. An example of these equations is also added.

The examples from the literature: An asymptotic preserving scheme for the two-fluid Euler–Poisson model in the quasineutral limit will be implemented.

• Modify the semidiscretization file to couple the Poisson equation
• Create 1st test case

[github-actions]

Review checklist

This checklist is meant to assist creators of PRs (to let them know what reviewers will typically look for) and reviewers (to guide them in a structured review process). Items do not need to be checked explicitly for a PR to be eligible for merging.

Purpose and scope

• The PR has a single goal that is clear from the PR title and/or description.
• All code changes represent a single set of modifications that logically belong together.
• No more than 500 lines of code are changed or there is no obvious way to split the PR into multiple PRs.

Code quality

• The code can be understood easily.
• Newly introduced names for variables etc. are self-descriptive and consistent with existing naming conventions.
• There are no redundancies that can be removed by simple modularization/refactoring.
• There are no leftover debug statements or commented code sections.
• The code adheres to our conventions and style guide, and to the Julia guidelines.

Documentation

• New functions and types are documented with a docstring or top-level comment.
• Relevant publications are referenced in docstrings (see example for formatting).
• Inline comments are used to document longer or unusual code sections.
• Comments describe intent ("why?") and not just functionality ("what?").
• If the PR introduces a significant change or new feature, it is documented in NEWS.md with its PR number.

Testing

• The PR passes all tests.
• New or modified lines of code are covered by tests.
• New or modified tests run in less then 10 seconds.

Performance

• There are no type instabilities or memory allocations in performance-critical parts.
• If the PR intent is to improve performance, before/after time measurements are posted in the PR.

Verification

• The correctness of the code was verified using appropriate tests.
• If new equations/methods are added, a convergence test has been run and the results
  are posted in the PR.

Created with ❤️ by the Trixi.jl community.

[DanielDoehring]

Good first effort!
Seeing this, I would in general opt for implementing the Multi-Ion Euler equations, similar to the Multi-Ion MHD equations. Then, we can derive the Electron-Ion as a special case for this.

[warisa-r]

Good first effort! Seeing this, I would in general opt for implementing the Multi-Ion Euler equations, similar to the Multi-Ion MHD equations. Then, we can derive the Electron-Ion as a special case for this.

Thank you for your detailed feedback! I have implemented the multi-ion system for compressible Euler equations. Is this now more in line with what you had in mind?

Running elixir_euler_multiion.jl yields

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'CompressibleEulerMultiIonEquations1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                197                run time:       8.55753000e-02 s
 Δt:             7.92328612e-03                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      2.00000000e+00 (100.000%)     time/DOF/rhs!:  1.84606481e-07 s
                                               PID:            3.91723582e-07 s
 #DOFs per field:            64                alloc'd memory:        158.932 MiB
 #elements:                  16

 Variable:       rho_1            rho_v1_1         rho_e_1          rho_2            rho_v1_2         rho_e_2          rho_3            rho_v1_3         rho_e_3
 L2 error:       2.31589449e-02   2.56244270e-02   2.69083460e-02   9.56221097e-03   1.90514680e-02   3.00279574e-02   2.31589449e-02   2.56244270e-02   2.69083460e-02
 Linf error:     3.31077549e-02   3.76854873e-02   4.56909872e-02   1.53340143e-02   3.08183559e-02   4.70978104e-02   3.31077549e-02   3.76854873e-02   4.56909872e-02
 Variable:       rho_1            v1_1             p_1              rho_2            v1_2             p_2              rho_3            v1_3             p_3
 L2 error prim.: 2.31589449e-02   3.35902165e-03   5.62572214e-03   9.56221097e-03   6.53763949e-03   1.18176493e-02   2.31589449e-02   3.35902165e-03   5.62572214e-03   
 Linf error pri.:3.31077549e-02   6.01291036e-03   1.08211138e-02   1.53340143e-02   1.15000182e-02   2.10914894e-02   3.31077549e-02   6.01291036e-03   1.08211138e-02
 ∑∂S/∂U ⋅ Uₜ :  -2.33315322e-03
────────────────────────────────────────────────────────────────────────────────────────────────────

P.S. The example also works with lax-friedrichs flux.

[DanielDoehring]

Yeah, looking good!

Now we need to find a way to make this system really coupled. Staying for the moment with the ion-only case, you could take a look at the work by Andres

trixi-framework#2213

and implement the ion-ion collision source terms. I am not sure if the testcase itself makes sense, or if you need for this the electron temperature/pressure business.

[warisa-r]

Yeah, looking good!

Now we need to find a way to make this system really coupled. Staying for the moment with the ion-only case, you could take a look at the work by Andres

trixi-framework#2213

and implement the ion-ion collision source terms. I am not sure if the testcase itself makes sense, or if you need for this the electron temperature/pressure business.

I can also implement the electrons' contribution. Based on this reference, it seems to not be much of a hassle since we only need additional input functions for electron pressure, temperature, and some constants to enable the ion-electron source term. However, I’d like to ask if we can simply adapt the existing source terms (with minor adjustments) from the multi-ion GLM-MHD equations for this purpose. Is what we are implementing the same with 13a), 13b), 13c) from this paper?

[DanielDoehring]

Taking a look at the PR by Andres, it seems the collision terms are purely hydrodynamic, so you should be able to directly re-use them for the Ion-Electron case.

What I am not sure about (maybe you need to read up on this) is the treatment of the electrons. As far as I understand, the electron pressure and temperature are simplifications that are based on neglecting the momentum equation for the electrons due to the significantly lower mass (and thus momentum) of the electrons compared to ions.
Additionally, also the continuity equation by assuming quasi-neutrality can be removed.

I really do not know if the testcase in the linked PR is still applicable for the full order model, for this I would need to consider those references. So I (have to) leave this assessment for the moment to you.

[warisa-r]

Taking a look at the PR by Andres, it seems the collision terms are purely hydrodynamic, so you should be able to directly re-use them for the Ion-Electron case.

What I am not sure about (maybe you need to read up on this) is the treatment of the electrons. As far as I understand, the electron pressure and temperature are simplifications that are based on neglecting the momentum equation for the electrons due to the significantly lower mass (and thus momentum) of the electrons compared to ions. Additionally, also the continuity equation by assuming quasi-neutrality can be removed.

I really do not know if the testcase in the linked PR is still applicable for the full order model, for this I would need to consider those references. So I (have to) leave this assessment for the moment to you.

I see. I will look it up and see what I can do. My last question at the moment is: what do you expect to see out of the result of running a successful simulation? Do you expect to see a temperature plot that is close to Figure 5d Ghosh et al. (the same one the PR you linked reproduced)?

[DanielDoehring]

Yes, so looking at that test you see that the magnetic fields are initialized with zero, and that they also do not appear in the source terms. So yes, I would hope that the temperature looks somewhat similar (depending on how much the explicit treatments of electrons change this).

[warisa-r]

The test case works for this implementation as well. This is the plot by our equations:

Looks the same as what the reference pull request (and the literature) produced. Like the reference PR, ignoring electron-ion collision terms also yields the same result of the temperature plot. Let me know if you want me to test out sth else! If everything looks OK then I'll proceed to clean up the code and extend the equations to higher dimensions (if it is necessary)

[DanielDoehring]

Awesome! Really good job!
Sounds good to me! For the moment 2D suffices for our purposes I think. In the meantime, if you need a break from programming you can maybe search for publications where people do actual simulations of the (originally intended) electron-ion system with Poisson coupling. It would be nice to have a meaningful testcase.

[warisa-r]

Got it!
As probably anticipated, the convergence test in 2D exhibited significantly higher error compared to the 1D case. Nevertheless, the 2D simulation successfully captured the temperature rise of both ion species, resulting in the same temperature plot.

The test case works for this implementation as well. This is the plot by our equations:

Looks the same as what the reference pull request (and the literature) produced. Like the reference PR, ignoring electron-ion collision terms also yields the same result of the temperature plot. Let me know if you want me to test out sth else! If everything looks OK then I'll proceed to clean up the code and extend the equations to higher dimensions (if it is necessary)

[warisa-r]

As of now, I will be working on the code under the assumption of a single ion species, since (as far as I understand) the current electron flux relies on a constant based on the ratio of m_e/ m_i. This term will become ambiguous for electron flux if we introduce multiple ion species. However, I will structure the code in a way that transitioning to a multi-species system will not be difficult in the future.

If I can find literature that reformulated the conservation laws so, that the ion flux is dependent on m_i / m_e (or if this is done by ourselves), I will change the flux function accordingly.

Let me know if I misunderstood anything!

[DanielDoehring]

As of now, I will be working on the code under the assumption of a single ion species, since (as far as I understand) the current electron flux relies on a constant based on the ratio of m_e/ m_i. This term will become ambiguous for electron flux if we introduce multiple ion species. However, I will structure the code in a way that transitioning to a multi-species system will not be difficult in the future.

If I can find literature that reformulated the conservation laws so, that the ion flux is dependent on m_i / m_e (or if this is done by ourselves), I will change the flux function accordingly.

Let me know if I misunderstood anything!

Yes exactly, only a single ion species, absolutely correct 👍

[warisa-r]

Here is the result of test case 1 that was run in the example

Compared to Fig 1. and Fig 2. in An asymptotic preserving scheme for the two-fluid Euler–Poisson model in the quasineutral limit, the sin waves that I obtained from my code have different amplitudes (in v1_ion) and there is a bit phase discrepancy (rho_electron, v1_electron). Maybe I did implement the initial and/or the boundary conditions incorrectly.

In addition: I, however, got an error message:

┌ Error: For periodic mesh non-periodic boundary conditions in x-direction are supplied.
└ @ Trixi C:\Users\user\Documents\Trixi.jl\src\semidiscretization\semidiscretization_hyperbolic.jl:234. 

But the example runs just fine. Do you know why this is the case? I tried to apply homogeneous Dirichlet boundary conditions to the Poisson equation and the periodic boundary conditions.

[DanielDoehring]

Just to have this posted somewhere: We should add at some point docstrings to these structs where we explicitly mention the paper we take the test cases from to familiarize people with the equations we use.

[DanielDoehring]

These parameters essentially govern how accurately the poisson equation is solved. You could experiment e.g. with stricter tolerance if this improves the result.

[warisa-r]

I see! I tried to adjust these a bit, including trying both lax friedrichs and hll flux. The result is still the same. Maybe I did sth wrong with the initial condition but I somehow cannot identify it.

[warisa-r]

I fixed some rather important bugs and read some more about the simulation. We do need a really small time step (delta_t < 10^-6) and really fine grid delta_x < 1e-4 (per instruction from the paper and also if these are two bigs the residual blows up to 2000-3000).

After fixing and trying this small grid size and time step, I found out that the simulation takes a very long time to run. Also, a more pressing matter is that the electric potential solver does not converge properly, leaving around 1e-3 to 1e-1 residual even with this set of parameters and increased maximum iterations.

My sincere apologies for the lousy report from before!

Let me know if I did something wrong or if you want me to try something! In the meantime, I'll look around to find if i can fix this issue or persist, run this on the cluster and send you the result despite the problem with the convergence of the electric potential solver.

@DanielDoehring

[DanielDoehring]

I see, thanks for the detailed report! The observations you report sound reasonable to me, so I think it is worth running these. If it is not feasible at the moment don't worry - I might get to it later. Also sorry for not reporting quickly to you - I am at the conference right now, i.e., long busy days + 7 hour timeshift. Anyways, great work as usual!

[warisa-r]

I see, thanks for the detailed report! The observations you report sound reasonable to me, so I think it is worth running these. If it is not feasible at the moment don't worry - I might get to it later. Also sorry for not reporting quickly to you - I am at the conference right now, i.e., long busy days + 7 hour timeshift. Anyways, great work as usual!

I wish you the very best of luck!

Another observation I want to write out and report is that I found that for electric potential solvers, we might not need the CFL number to be as small as the one we need for the time step of the Euler's equations (lower than 1e-6). This makes the simulation go faster. However, for certain time steps, I run into large residues, or it just takes much, much longer for the electric potential to converge. It might be worth it to experiment with these numbers before investing time in running the full simulation.