some doubts in simulating coulomb explosion

[RayleighZZ]

Hello devs:

Recently, I am using epoch-4.17.13 to simulate coulomb explosion, gcc version 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04), 15 cores. In short, the physical process I want to reproduce is about many deuterons are presented in a very small spherical area and the simulation area is a cube with a side length 30 times the radius of the spherical deuteron system (Attached my input.deck). The spherical deuteron system is of solid density and theoretically, deuterons will diffuse within tens of femtoseconds and be accelerated to 10^-16J, due to the coulomb repulsion between them. But according to my simulation results, when the simulation time is one picosecond, the average kinetic energy of the deuterons is still almost equal to the initial moment 10^-21J, and the deuteron system has almost no diffusion (Attached my results, at the initial moment and one picosecond, the
deuteron number density distribution in the x-y plane where the center of the deuteron system is located). The average weight of peusoparticles is almost 1.

I analyzed the reasons. Theoretically, the maximum electric field strength in the simulation area can reach 10^12 V/m, but when I tested results in 1 femtosecond, the maximum electric field strength of the system can only reach 10^7 V/m.
So, what is the reason for this and how can I solve it?
And, is there something wrong with my input.deck?
Thank you very much!
My post analysis script for matlab is also attached. You could copy that for analysis.
inputdeck.txt


matlab_script.txt

[TomGoffrey]

Your initial conditions are not charge neutral - but you start with E=0. This is inconsistent.

EPOCH solves Maxwell's equations in a manner consistent with Gauss' law, but only if your initial conditions are also consistent with Gauss' law. The discussion in #158 might be useful.

The particles in EPOCH only couple to the EM fields via current deposition, and that's the electric field you are observing.

[TomGoffrey]

No problem.

[RayleighZZ]

Hi Tom,
@TomGoffrey After our last discussion, I tried some methods to modify the initial electric field, but I didn't succeed. The initial electric field I want to add is spatially varying in the whole simulation area. I found a subroutine set_initial_values in src/housekeeping/setup.F90, the initial electric field in each direction is set to 0 there. I tried to modify this subroutine, but it seems that the modifications are only for the area controlled by a certain processor due to the MPI in Epoch (I am using epoch-4.17.13, gcc version 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04), 15 cores). I'm not very familiar with MPI. Are there some ways to set a spatially varying initial electric field? Should I make some modifications in src/fields.F90?
Looking forward to your reply. Thank you very much!

[TomGoffrey]

Perhaps it would be easier for you to use an external tool to calculate the initial field, and supply the data to EPOCH in a file.

You can do this by specifying ex = 'ex.dat' in the fields block. The file should be raw binary, stored in column-major order.

[RayleighZZ]

Ok. I'll try this. And, whether the electric field set by fields.block is only loaded into the simulation area at the initial moment, or is it loaded every moment? Do I need to consider the timing control of the applied electric field?

[RayleighZZ]

The applied electric field set by fields.block is loaded every timestep. I should consider the time control of the applied electric field.