Long-term deviation between numerical and theoretical results in soliton propagation

[jinze-liu]

I am performing a nonlinear Schrödinger equation (NLSE) soliton simulation using WarpX. Initially, the numerical solution closely follows the theoretical prediction, but as time progresses, the discrepancy between the two grows. Eventually, the soliton shape begins to distort significantly. To address this, I have tried grid resolution，CFL and boundary_conditions. However, in all cases, the deviation still accumulates over time, and the soliton structure eventually deforms. I suspect that numerical dispersion or boundary effects may be influencing the simulation.

What could be the primary cause of this increasing discrepancy? Any insights or suggestions would be greatly appreciated!

Below are my input file settings and comparison plots of the numerical vs. theoretical results:

nx, ny, nz = 128, 32, 16384
Lx, Ly, Lz = 20, 1, 50 
max_grid_size = 512
blocking_factor = 32 
grid = picmi.Cartesian3DGrid(
    number_of_cells=[nx, ny, nz],
    lower_bound=[0, -Ly/2, -Lz/2],
    upper_bound=[Lx, Ly/2, Lz/2],
    lower_boundary_conditions=["periodic", "periodic", "periodic"],
    upper_boundary_conditions=["periodic", "periodic", "periodic"],
    lower_boundary_conditions_particles=["periodic", "periodic", "periodic"],
    upper_boundary_conditions_particles=["periodic", "periodic", "periodic"],
    warpx_max_grid_size=max_grid_size,
    warpx_blocking_factor=blocking_factor,
)

solver = picmi.ElectromagneticSolver(
    grid=grid,
    method = 'Yee',
    cfl=0.5
)

density = 7.9e8  
temperature = 2320  # in Kelvin

electrons = picmi.Species(
    particle_type='electron',
    name='electrons',
    initial_distribution=picmi.AnalyticDistribution(
        density_expression="7.9e8 ",
        upper_bound=[None, None, None],
        momentum_expressions = ["0.0", "0.0", "0.0"]
    )
)

sim = picmi.Simulation(
    solver=solver,
    max_steps=20000,
    verbose=1,
    particle_shape=3,
    warpx_use_filter = 0,
    warpx_current_deposition_algo='direct',
)

sim.add_species(electrons, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[1, 1, 10]))

external_Efield = picmi.AnalyticInitialField(
    Ex_expression="11712.8*cos(2*z)*(1 - 44998/(-1 + 150*cosh(0.299993*(z))))",
    Ey_expression="0.0",
    Ez_expression="0.0" 
)

external_Bfield = picmi.AnalyticInitialField(
    Bx_expression="0.0",  
    By_expression="0.0000139604*cos(2*z)*(1 - 44998/(-1 + 150*cosh(0.299993*(z))))",  
    Bz_expression="0.0"   
)

sim.add_applied_field(external_Efield)  # Correct method to add applied fields
sim.add_applied_field(external_Bfield)  # Add the magnetic field
# Define and add diagnostics
field_diag = picmi.FieldDiagnostic(
    name='field_diag',
    grid=grid,
    period=1000,
    data_list=['Ex', 'Ey', 'Bx', 'By', 'rho' ] ,  # 'Ex', 'Ey', 'Ez', 'Bx', 'By', 'Bz'
)

Z1field_probe_scat_line_rdiag = picmi.ReducedDiagnostic(
    diag_type='FieldProbe',
    name='0Ex-z',
    period=100,
    integrate=0,
    probe_geometry='Line',
    x_probe=Lx/2,
    y_probe=0,
    z_probe=-Lz/2,
    x1_probe=Lx/2,
    y1_probe=0,
    z1_probe=Lz/2,
    resolution=1000
)

sim.add_diagnostic(Z1field_probe_scat_line_rdiag)
sim.add_diagnostic(field_diag)
sim.step()

[jinze-liu]

After further testing, I found that the issue persists even in a simple wave propagation case. Initially, I set up a plane wave:

Ex_expression = "10000.0 * cos(2 * z)"

and tested without any particles, which showed stable propagation over long time steps. However, after introducing even a small number of electrons, the wave begins to decay significantly after a few hundred time steps.

To investigate further, I also added ions to ensure charge neutrality, but the issue remained.

And below is a comparison of Ex over time, showing how the wave decays unexpectedly:

nx, ny, nz = 128, 32, 16384
Lx, Ly, Lz = 20, 1, 50  
max_grid_size = 512
blocking_factor = 32 

# Define the simulation grid
grid = picmi.Cartesian3DGrid(
    number_of_cells=[nx, ny, nz],
    lower_bound=[0, -Ly/2, -Lz/2],
    upper_bound=[Lx, Ly/2, Lz/2],
    lower_boundary_conditions=["open", "open", "open"],
    upper_boundary_conditions=["open", "open", "open"],
    lower_boundary_conditions_particles=["open", "open", "open"],
    upper_boundary_conditions_particles=["open", "open", "open"],
    warpx_max_grid_size=max_grid_size,
    warpx_blocking_factor=blocking_factor,
)

# Electromagnetic solver
solver = picmi.ElectromagneticSolver(
    grid=grid,
    method='Yee',
    cfl=0.1
)

# External electric field
external_Efield = picmi.AnalyticInitialField(
    Ex_expression="10000.0 * cos(2 * z)",
    Ey_expression="0.0",
    Ez_expression="0.0" 
)

external_Bfield = picmi.AnalyticInitialField(
    Bx_expression="0.0",
    By_expression="0.0",
    Bz_expression="0.0"
)

# Initialize simulation
sim = picmi.Simulation(
    solver=solver,
    max_steps=20000,
    verbose=1,
    warpx_use_filter=0,
    warpx_current_deposition_algo="esirkepov",
)
# Add a small number of electrons
electrons = picmi.Species(
    particle_type='electron',
    name='electrons',
    initial_distribution=picmi.AnalyticDistribution(
        density_expression="1e6",
        upper_bound=[None, None, None],
        momentum_expressions=["0.0", "0.0", "0.0"]
    )
)
ions = picmi.Species(
    particle_type='He',
    name='ions',
    charge_state=+1,
    initial_distribution=picmi.AnalyticDistribution(
        density_expression="1e6",  
        momentum_expressions=["0.0", "0.0", "0.0"]
    )
)
sim.add_species(ions, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[1, 1, 10]))
sim.add_species(electrons, layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[1, 1, 10]))
# Apply external fields
sim.add_applied_field(external_Efield)
sim.add_applied_field(external_Bfield)

# Run simulation
sim.step()

[RemiLehe]

Sorry for the late reply!

For the case that you showed (soliton): my guess is that what you are observing might be related to the exact initialization of the simulation:

• physically, I think that the electrons that overlap with the soliton at t=0 should have a non-zero velocity at t=0 (to be consistent with the fact that they are feeling the fields of the soliton)
• however, in your simulation script above, the electrons have no velocity at t=0 (from the line momentum_expressions = ["0.0", "0.0", "0.0"])

Do you have an analytical expression for the velocity of the electrons as a function of z (similar to the expression that you have for E and B)? If yes, I would recommend incorporating this expression in your simulation script, for the initialization.

Regarding the second case ("plane wave"): thanks for your efforts in providing a simpler example that reproduces the issue. However, I am not sure that this example is representative, as it does not quite correspond to a plane wave, since the B field is 0 in this case. I would recommend initializing the B field consistently with the E field, and maybe using periodic boundary conditions.

[jinze-liu]

Thanks for your reply! I will carefully check the initial conditions in theory.