Thomas-Fermi: Not Charge Conserving

[ax3l]

The Thomas-Fermi ionization creates either electrons with too high weighting or too many of them.

See the comparison from the new LCT foil example in #2008

#!/usr/bin/env python
#

import matplotlib.pyplot as plt
import h5py as h5
import numpy as np

starts = {
    "random": "004",
    "quiet": "005"
}
start = starts["random"]
sims = {
    "BSI": "lctFoil-" + start + "-BSI",
    "ADK": "lctFoil-" + start + "-ADK",
    "TF": "lctFoil-" + start + "-TF",
    "combined": "lctFoil-" + start + "-BSI-ADK-TF"
}
dt = 4.91356e-18 # s
dx = 800.e-9 / 384. # m
step = 750

fig = plt.figure()
cax = fig.add_axes([0.87, 0.1, 0.01, 0.75])
ax1 = fig.add_axes([0.1, 0.5, 0.35, 0.35])
ax2 = fig.add_axes([0.5, 0.5, 0.35, 0.35])
ax3 = fig.add_axes([0.1, 0.1, 0.35, 0.35])
ax4 = fig.add_axes([0.5, 0.1, 0.35, 0.35])
axes = [ax1, ax2, ax3, ax4]

def get_nZ(flds, species):
    r = flds[species + "_chargeDensity"]
    d = r[()] * r.attrs["unitSI"] / 1.602e-19 * 1.e-6 # elements / cm^3
    return d

def plot_sim(ax, sim):

    f = h5.File("../" + sims[sim] + "/simOutput/h5/simData_" + str(step) + ".h5" , "r")
    ne = get_nZ(f["/data/" + str(step) + "/fields/"], "e")
    nH = get_nZ(f["/data/" + str(step) + "/fields/"], "H")
    nC = get_nZ(f["/data/" + str(step) + "/fields/"], "C")
    nN = get_nZ(f["/data/" + str(step) + "/fields/"], "N")

    d = ne + nH + nC + nN
    ax.set_title(sim)
    return ax.imshow(
      #np.abs(d),
      d,
      #cmap='CMRmap_r',
      cmap='RdBu',
      origin="lower",
      aspect="auto",
      interpolation="nearest",
      #vmin=0.,
      #vmax=2.e24,
      vmin=-3.e22, vmax=3.e22,
      extent=[0., dx * 1.e6 * d.shape[0], 0., dx * 1.e6 * d.shape[1]]
    )
    f.close()

i = 0
for sim in sims:
    print(sim)
    im = plot_sim(axes[i], sim)
    i += 1

#fig.colorbar(im, cax=cax, label=r'$n_e$ [$q_e \cdot$ cm$^{-3}$]')
#fig.colorbar(im, cax=cax, label=r'$n_{Z,C}$ [$q_e \cdot$ cm$^{-3}$]')
fig.colorbar(im, cax=cax, label=r'$n_{Z,\Sigma{e,H,C,N}}$ [$q_e \cdot$ cm$^{-3}$]')
fig.suptitle("time = {:5.3f} fs".format(dt * step * 1.e15))
ax3.set_xlabel(r'$x$ [$\mu$m]')
ax4.set_xlabel(r'$x$ [$\mu$m]')
ax1.set_ylabel(r'$y$ [$\mu$m]')
ax3.set_ylabel(r'$y$ [$\mu$m]')
plt.show()

Comparisons

(laser from bottom plane)

Quiet Start

n_e

n_Z,C

sum(n_Z)

Random Start

n_e

n_Z,C

sum(n_Z)

[n01r]

I can reproduce this result.

I checked with only carbon now and found that the weighting seems to be correct but the number of macroelectrons is way too large!

Update:
In AlgorithmThomasFermi.hpp I then set newBoundElectrons to an arbitrary constant value (3.0). With that I wanted to exclude the possibility that with a fixed number of electrons to create the wrong number is actually spawned. I observed that the number of electrons created is the same as I'd expect.

[ax3l]

Current state of the FoilLCT example in dev as of 69011eff40f4a7e828f18cbae86343349f01d6ec (step 2000 with default settings, e.g. random start):