Spontaneous Radiation Rises as Resolution increases

[Xuan0533]

I got trouble with the resolution setup. It is found that as the transverse resolution increases by increasing number of nodes keep, the spontaneous radiation continuous to increase with no end.

A specific image file is attached:

I would really appreciate if I can get some help!

[Xuan0533]

The input file:

&MDATA
  qOneD                  = .false.
  qScaled                = .false.  ! Using unscaled (S.I.) inputs
  sElectronThreshold     = 0.05
  qFocussing             = .true.   ! if focussing is included in the transverse plane
  qUndEnds               = .false.
  ioutInfo=3
  meshType  = 0          ! Full mesh
  sFModelLengthZ2        = 1.5674743331935322e-15
  sFModelLengthX         = 0.00025
  sFModelLengthY         = 0.00025
  nodesPerLambdar        = 20  ! <- Nodes per resonant wavelength in the temporal direction
  iNumNodesX             = 201
  iNumNodesY             = 201
  qFMesh_G = .false.            ! <- Auto-setup of the mesh: 
  srho                   = 0.00473540795492875
  saw                    = 3.0
  sgamma_r               = 8727.533
  lambda_w               = 0.03
  zundType               = 'helical'
  sZ0                    = 0.0
  stepsPerPeriod         = 30  ! Number of steps per wiggler period
  iWriteNthSteps         = 1500
  iWriteIntNthSteps      = 30
  qInitWrLat             = .false.
  qDumpEnd               = .false.
  Lattfile= 'FACET.latt'
  beam_file = 'FACET.beam'
/

The Beam file:

&NBLIST
nbeams = 1
dtype = 'simple'
/

&BLIST
sSigmaE = 1.1181e-5, 1.1181e-5, 1.1129e-16, 6.1339e-6, 6.1339e-6, 0.00065187
sLenE = 1.1181e-4, 1.1181e-4, 1.1129e-15, 6.1339e-5, 6.1339e-5, 0.0065187
bcenter = 0.0000
gammaf = 1.0
Ipk = 11964.77   ! Peak Current
alphax = 0.0   ! Twiss params...
alphay = 0.0
emitx = 7.3104e-11  ! geometric emittances...
emity = 7.3104e-11
TrLdMeth = 2
iMPsZ2PerWave = 20
nseqparts = 1000

The Lattice file:

QU -1.240 1.240
UN  'helical'  25   1.0   0.0   30   1.0   1.0   0.0    0.0
QU 0.620 -0.620
UN  'helical'  25   1.0   0.0   30   1.0   1.0   0.0    0.0
QU -1.240 1.240

[mightylorenzo]

Hmm..there's one or two possibilities here...

First of all, just a minor point, are you filtering the radiation field before plotting the power? This may all be low-frequency content, but I can't see a power spectrum so I don't know...usually it's OK since there is a high-pass filtering that Puffin itself does when running in 3D (defaulted at a 1/3rd of the resonant frequency, but it can be adjusted), but it can also depend on the temporal profile of the beam if it's short enough. Usually I'd actually just band-pass the field around the resonant frequency and also the harmonics separately if necessary, but it depends on what you're looking for and what scenario you're investigating.

Another thing is, the radiation field mesh uses an absorbing boundary of 16 nodes surrounding the outer edge on all sides. They are included in the mesh setup in the input file, so if you tell puffin to setup 100 nodes in the x-direction for the mesh, the inner 68 will be non-absorbing nodes surrounded by 16 absorbing nodes on each side (32 absorbing nodes in total). This number of absorbing nodes is currently fixed in puffin. If you have a lot of diffraction, then you may see power being sucked out of the system through the absorbing boundaries. Changing the grid spacing while keeping the length constant would change the length of the absorbing boundary, resulting in more or less power being absorbed by the boundaries.

I've not had a chance to look properly at the numbers, but if you have a lot of diffraction (Rayleigh range comparable to the undulator period) you may also find that calculating the diffraction step every period is not good enough. Having said that, I seem to recall that in practice this didn't seem to affect the results significantly even when the Rayleigh range was pretty much equal to the undulator period, but I may be remembering wrong :) . You can adjust the diffraction step size with the 'sDiffFrac' variable in the main input file. It is the length of the diffraction step expressed in units of the undulator period, defaulted to 1. You could try reducing it. Worst case is calculating a diffraction step at the end of every RK4 temporal-longitudinal step (sDiffFrac = 0.03333333333333333)...this will massively increase your simulation time...so I hope it's not that...

The last thing I can think of is simply an under-sampling of the beam compared to the radiation grid. My quick (maybe wrong) calc from looking at the input file looks like there's roughly 10,000 transverse grid elements (100 in x and 100 in y) and only 1000 randomly distributed (OK, well it's a Halton sequence but almost :) ) macroparticles for each slice in the transverse direction....for context, the way I usually set it up is to have the beam contained within around 20 or so nodes in the centre of the grid in each of x and y when using 1000 macroparticles in each transverse slice. In the example file you gave, the beam may not be being properly sampled by the field, if that makes any sense. This interplay between the beam and grid sampling is crucial, rather than the absolute individual sampling of the beam and field grid. Maybe my numbers for the file here are incorrect...but as I say I'd aim to have about 20 or so nodes in the centre of the grid in each of x and y sampling the beam here (usually 5 or 6 times the standard deviation if it's a Gaussian beam), with the outer nodes there to allow for the field diffraction (usually 64 or 128 total, but depends on the amount of diffraction you need to model). With a Gaussian in x and y that's fine, adding more doesn't give you much, if any, benefit (certainly not worth the extra compute tbh). OTOH maybe your focussing system is causing the beam to stretch/contract significantly in x or y across the length of the undulator but hopefully this can point you in the right direction.

I'd suggest testing this in periodic mode if you're not already - you'll be able test all this relatively quickly to pin down/eliminate the transverse variations you're investigating before moving to a full-fat simulation.

[Xuan0533]

Thank you very much for the suggestions.

I have experimented with changing various setup parameters such as sDiffFrac, sFiltFrac, and nseqparts, but in all cases, I observe a similar phenomenon: the spontaneous radiation power increases as the transverse resolution increases.

For reference, I am using the 3D example input file for LCLS, with two modifications:

• qElectronFieldCoupling = .false.

• sperwaves = 5

I increase the resolution by scaling both iNumNodes and iRedNodes by factors of 2, 3, 4, 5, etc. For each case, maximum power and the 3D spectrum of radiation at z = 3m is shown below:

[Xuan0533]

The 3D spectrum code I use is:

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

h5_fld_filename = out_dir +  'lcls_aperp_90.h5'
with h5py.File(h5_fld_filename, 'r') as f:
    for item in f.keys():
        print(f"Item: {item}")
        print(f"Data: {f[item]}")
        # Print attributes
        for attr_name in f[item].attrs:
            print(f"  Attribute: {attr_name} = {f[item].attrs[attr_name]}")
            
with h5py.File(h5_fld_filename, 'r') as f:
    #aperp = f['aperp'][:]
    gamma_r = f['/runInfo'].attrs['gamma_r']
    lc = f['/runInfo'].attrs['Lc']
    lg = f['/runInfo'].attrs['Lg']
    kappa = f['/runInfo'].attrs['kappa']
    Ex_SI = (f['aperp'][0]*gamma_r*electron_mass*c**2)/(e*lg*kappa)
    Ey_SI = -(f['aperp'][1]*gamma_r*electron_mass*c**2)/(e*lg*kappa)
    mesh_lb = f['/globalLimits'].attrs['vsLowerBounds']
    mesh_ub = f['/globalLimits'].attrs['vsUpperBounds']
    xlamds = f['runInfo'].attrs['lambda_r']
    xlamd = f['runInfo'].attrs['lambda_w']
    z2Length = (f['runInfo'].attrs['nZ2']) * f['runInfo'].attrs['sLengthOfElmZ2'] * lc
    dz = f['runInfo'].attrs['sLengthOfElmZ2'] * lc`

nz, nx, ny = Ex_SI.shape
mesh_z = np.linspace(mesh_lb[0], mesh_ub[0], nz+1)*lc
mesh_x = np.linspace(mesh_lb[1], mesh_ub[1], nx+1)*np.sqrt(lc*lg)
mesh_y = np.linspace(mesh_lb[2], mesh_ub[2], ny+1)*np.sqrt(lc*lg)
dz = np.diff(mesh_z)[0]
dx = np.diff(mesh_x)[0]
dy = np.diff(mesh_y)[0]
kx = np.fft.fftshift(np.fft.fftfreq(nx, dx))
ky = np.fft.fftshift(np.fft.fftfreq(ny, dy))
kz = np.fft.fftshift(np.fft.fftfreq(nz, dz))
freq_r = c/xlamds

# reciprocal spacings
dkx, dky, dkz = 1 / (nx * dx), 1 / (ny * dy), 1 / (nz * dz)
#dkx, dky, dkz = np.diff(kx)[0], np.diff(ky)[0], np.diff(kz)[0]

Kz, Kx, Ky = np.meshgrid(kz, kx, ky, indexing="ij")
K = np.sqrt(Kx**2 + Ky**2 + Kz**2)
photon_energy_3d = K.flatten() * c * h / e  # photon energy, eV
dEk = (np.abs(FFT(Ex_SI).flatten())**2 + np.abs(FFT(Ey_SI).flatten())**2) * dkx * dky * dkz * epsilon_0   # Spectral Energy Density

bins = 50
hist, bin_edges = np.histogram(photon_energy_3d, bins=bins, weights=dEk) # Computes the histogram of photon energies (photon_energy_3d) with weighting each entry by the spectral energy density (dEk).
hist_x = bin_edges[:-1] + np.diff(bin_edges) / 2  # Centers of the histogram bin
hist_width = np.diff(bin_edges)  # Width of the histogram bin
hist_y = hist / hist_width  # Normalizes the histogram counts (hist) by the bin widths to get the spectral energy at each photon energies (photon_energy_3d)

fig, ax = plt.subplots()
ax.bar(hist_x, hist_y * 1e6, hist_width)
ax.axvline(x = (freq_r*h/e), color="k", linestyle = '--')
ax.set_xlabel("Photon Energy (eV)")
ax.set_ylabel("Spectral Energy (µJ/eV)")
ax2 = ax.twinx()
ax2.plot(hist_x, np.cumsum(hist) * 1e6,  color="r")
ax2.set_ylabel("Cumulative Energy (µJ)", color="r")
ax2.tick_params(axis='y', colors='red')
ax2.set_ylim(0, None)

[mightylorenzo]

Thanks for the work and the plots!

I think the problem here is that you are increasing the transverse resolution of the mesh without increasing the sampling resolution of the beam. As I said above, it's not just the node spacing you need to worry about, it's how the mesh is 'sampling' the beam. Doing the same simulation, I see the same behaviour, the spontaneous emission power increases. But if, when doubling the field resolution, I also double nseqparts (doubling the number of electrons in the transverse plane), then I get identical results to the non-doubled example.

What you are observing (I think) is the field oversampling the beam in the transverse plane, so that the charge density is seen by the field to have an artificial filamentation, exciting many coherent transverse modes. To ensure you avoid this, you need to ensure the node spacing is sufficiently large so that the mesh sees the beam as continuous in the transverse plane - unless the filamentation is genuine!! 😀

[Xuan0533]

Thanks for your help! I just tried simulation with double the resolution and double nesqparts, and now, the synchrotron radiation issue is gone! The oversampling makes sense to me.

May I ask about what should be a correct large node spacing, like how many macroparticles should be put per cell? Based on your suggestion, for ebeam contained in 20 or so nodes should have at least 1000 slice macroparticles for the field to sample the beam correctly, so it is like 50 macripaticles per cell. Is this correct?

Actually, I am using PUFFIN with h5 beam file with equal weight particles. For this case, how should I do convergence test to find the proper mesh setup?

A couple of follow‑ups, please:

1. Based on pervious discussion, the ebeam contained in 20 nodes should have at least 50 macroparticles per cell (i.e. 1 000 per slice in total), is this reasonable? If not, may I ask about what correct 'large' node spacing should be?

2. Since I’m loading an equal‑weight HDF5 beam, what’s the best strategy to do convergence test to find proper mesh setup?

Thanks again for your guidance!