Electrostatic solver failing to converge with very high tolerance

[archermarx]

Hi all,

I'm running a 2D simulation using this PICMI input

import numpy as np
from math import sqrt, ceil, floor
from pywarpx import callbacks, fields, libwarpx, particle_containers, picmi
from periodictable import elements

# constants
m_p = picmi.constants.m_p
m_e = picmi.constants.m_e
q_e = picmi.constants.q_e
ep0 = picmi.constants.ep0

####################################################################
#                     CONFIGURABLE OPTIONS                         #
####################################################################

verbose = True             # Whether to use verbose output
num_grids = 1              # Number of subgrids to decompose domain into
particles_per_cell = 100    # Number of particles per cell
dt = 1e-12
seed = 11235813             # Random seed

L_azimuthal = 1e-3        # Simulation domain length (m)
L_axial = 1e-3     # Virtual axial length
max_time = 1e-6    # Max time (s)
num_diags = 1000     # Number of diagnostic outputs
n0 = 1e17           # Plasma density
B0 = 2e-2           # Magnetic field strength (T)
species = "Xe"      # Ion species
T_e = 2.0          # Electron temperature (eV)
T_i = 0.2           # Ion temperature (eV)
V_d = 200           # discharge voltage

# number of cells
Nx = 32 # axial cells
Nz = 32 # azimuthal cells

####################################################################
#                        DERIVED VALUES                            #
####################################################################

m_i = elements.symbol(species).mass * m_p       # Ion mass
ve_rms = sqrt(q_e * T_e / m_e)                  # Electron rms speed
vi_rms = sqrt(q_e * T_i / m_i)                  # Ion rms speed
u_i = sqrt(2 * q_e * V_d / m_i) / 2             # ion directed velocity

max_steps = ceil(max_time / dt)                 # Maximum simulation steps
diag_inter_time = max_time / num_diags          # Interval between diagnostic outputs (s)
diag_inter_iter = round(diag_inter_time / dt)   # Interval between diagnostic outputs (iters)

####################################################################
#                       SIMULATION SETUP                           #
####################################################################

# Grid
grid = picmi.Cartesian2DGrid(
    number_of_cells = [Nx, Nz],
    warpx_max_grid_size = Nx * Nz,
    warpx_blocking_factor = 1,
    lower_bound = [0, 0],
    upper_bound = [L_axial, L_azimuthal],
    lower_boundary_conditions = ['dirichlet', 'periodic'],
    upper_boundary_conditions = ['dirichlet', 'periodic'],
    lower_boundary_conditions_particles = ['absorbing', 'periodic'],
    upper_boundary_conditions_particles = ['absorbing', 'periodic'],
    warpx_potential_lo_x = 0.0,
    warpx_potential_hi_x = V_d,
)

# Field solver
solver = picmi.ElectrostaticSolver(
    grid=grid,
)

# Initialize simulation
sim = picmi.Simulation(
    solver = solver,
    time_step_size = dt,
    max_time = max_time,
    verbose = verbose,
    warpx_use_filter = True,
    warpx_field_gathering_algo = 'energy-conserving',
    warpx_serialize_initial_conditions = True,
    warpx_random_seed = seed,
    warpx_sort_intervals = 500
)
solver.sim = sim

# Applied fields
external_field = picmi.ConstantAppliedField(
    By = B0
)
sim.add_applied_field(external_field)

# Particles
gridded_layout = picmi.GriddedLayout(n_macroparticle_per_cell = [particles_per_cell, particles_per_cell], grid = grid)

# electron distribution
electron_distribution = picmi.UniformDistribution(
    density = n0, rms_velocity = [ve_rms, ve_rms, ve_rms], directed_velocity = [-u_i, 0.0, 0.0]
)

# electron definition
electrons = picmi.Species(
    particle_type = 'electron', name = 'electrons', initial_distribution = [electron_distribution]
)

sim.add_species(electrons, layout = [gridded_layout])

# ion distribution
ion_distribution = picmi.UniformDistribution(
    density = n0, rms_velocity = [vi_rms, vi_rms, vi_rms], directed_velocity = [u_i, 0.0, 0.0]
)

ions = picmi.Species(
    particle_type = species, name = 'ions', mass = m_i, charge = 'q_e',
    initial_distribution = [ion_distribution]
)

sim.add_species(ions, layout = [gridded_layout])

diag_name = "ecdi_2d"

# Diagnostics
field_diag = picmi.FieldDiagnostic(
    name = diag_name,
    grid = grid,
    period = diag_inter_iter,
    data_list = ['rho_ions', 'rho_electrons', 'Ez', 'j'],
    write_dir = 'diags',
    warpx_format = "openpmd"
)
sim.add_diagnostic(field_diag)

particle_diag = picmi.ParticleDiagnostic(
    name = diag_name,
    period = diag_inter_iter,
    write_dir = 'diags',
    warpx_format = "openpmd"
)
sim.add_diagnostic(particle_diag)

# Initialize simulation
sim.initialize_inputs()
sim.initialize_warpx()

####################################################################
#                        RUN SIMULATION                            #
####################################################################
sim.step(max_steps)

On the first iteration, the electrostatic solver fails to converge, and the residual is quite large:

MLMG: Iteration 193 Fine resid/bnorm = 14.75470384
MLMG: Iteration 194 Fine resid/bnorm = 17.49102104
MLMG: Iteration 195 Fine resid/bnorm = 17.19491238
MLMG: Iteration 196 Fine resid/bnorm = 17.72385067
MLMG: Iteration 197 Fine resid/bnorm = 17.77645028
MLMG: Iteration 198 Fine resid/bnorm = 17.1501666
MLMG: Iteration 199 Fine resid/bnorm = 14.85552749
MLMG: Iteration 200 Fine resid/bnorm = 22.35247512
MLMG: Failed to converge after 200 iterations. resid, resid/bnorm = 0.000286684874, 22.35247512
amrex::Abort::0::MLMG failed. !!!
SIGABRT

The full output file is here: error.txt

Any idea what's going on? I don't think I'm trying to do anything too complex.

[ax3l]

Hi @archermarx,

It looks like your problem has a hard time converging in the Poisson solver... You can try to increase the maximum number of allowed steps or the tolerance when the result is accepted (maximum_iterations and required_precision): https://warpx.readthedocs.io/en/latest/usage/python.html#pywarpx.picmi.ElectrostaticSolver

[archermarx]

Hi,

I realize that, but I'm more wondering why this problem is so difficult that it fails on the first iteration? I have dirichlet boundary conditions and an initially-uniform plasma. That doesn't seem to me like a problem that should be challenging for a poisson solver. I'm curious whether I'm doing something wrong here.

On Mon, Jan 29, 2024, 7:34 PM Axel Huebl @.***> wrote:

Hi @archermarx https://github.com/archermarx,

It looks like your problem has a hard time converging in the Poisson solver... You can try to increase the maximum number of allowed steps or the tolerance when the result is accepted (maximum_iterations and required_precision):

https://warpx.readthedocs.io/en/latest/usage/python.html#pywarpx.picmi.ElectrostaticSolver

— Reply to this email directly, view it on GitHub https://github.com/ECP-WarpX/WarpX/issues/4641#issuecomment-1915833088, or unsubscribe https://github.com/notifications/unsubscribe-auth/AOP3QGUB2LEQ6O2T3GPP4O3YRA5YRAVCNFSM6AAAAABCMLRN2SVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSMJVHAZTGMBYHA . You are receiving this because you were mentioned.Message ID: @.***>

[WeiqunZhang]

Could you provide an inputs file that can be used for C++ code?

https://warpx.readthedocs.io/en/latest/usage/python.html#pywarpx.picmi.Simulation.write_input_file

[archermarx]

Certainly

# algo
algo.field_gathering = energy-conserving
algo.particle_shape = 1

# amr
amr.blocking_factor = 1
amr.max_grid_size = 1024
amr.max_level = 0
amr.n_cell = 32 32

# boundary
boundary.field_hi = pec periodic
boundary.field_lo = pec periodic
boundary.particle_hi = absorbing periodic
boundary.particle_lo = absorbing periodic
boundary.potential_hi_x = 200
boundary.potential_lo_x = 0.0

# diagnostics
diagnostics.diags_names = ecdi_2d

# ecdi_2d
ecdi_2d.diag_type = Full
ecdi_2d.fields_to_plot = Ez rho_electrons rho_ions
ecdi_2d.file_prefix = diags/ecdi_2d
ecdi_2d.format = openpmd
ecdi_2d.intervals = 1000
ecdi_2d.species = electrons ions
ecdi_2d.write_species = 1

# electrons
electrons.charge = -q_e
electrons.dist0.density = 1e+17
electrons.dist0.injection_style = nuniformpercell
electrons.dist0.momentum_distribution_type = gaussian
electrons.dist0.num_particles_per_cell_each_dim = 100 100
electrons.dist0.profile = constant
electrons.dist0.ux_m = -2.8491467279331137e-05
electrons.dist0.ux_th = 0.0019783585031909016
electrons.dist0.uy_m = 0.0
electrons.dist0.uy_th = 0.0019783585031909016
electrons.dist0.uz_m = 0.0
electrons.dist0.uz_th = 0.0019783585031909016
electrons.initialize_self_fields = 0
electrons.injection_sources = dist0
electrons.mass = m_e

# geometry
geometry.dims = 2
geometry.prob_hi = 0.001 0.001
geometry.prob_lo = 0 0

# ions
ions.charge = q_e
ions.dist0.density = 1e+17
ions.dist0.injection_style = nuniformpercell
ions.dist0.momentum_distribution_type = gaussian
ions.dist0.num_particles_per_cell_each_dim = 100 100
ions.dist0.profile = constant
ions.dist0.ux_m = 2.8491467279331137e-05
ions.dist0.ux_th = 1.2741771523059082e-06
ions.dist0.uy_m = 0.0
ions.dist0.uy_th = 1.2741771523059082e-06
ions.dist0.uz_m = 0.0
ions.dist0.uz_th = 1.2741771523059082e-06
ions.initialize_self_fields = 0
ions.injection_sources = dist0
ions.mass = 2.196035502270312e-25

# particles
particles.B_ext_particle_init_style = constant
particles.B_external_particle = 0.0 0.02 0.0
particles.species_names = electrons ions

# stop_time
stop_time = 1e-06

# warpx
warpx.const_dt = 1e-12
warpx.do_electrostatic = labframe
warpx.random_seed = 11235813
warpx.serialize_initial_conditions = 1
warpx.sort_intervals = 500
warpx.use_filter = 1
warpx.verbose = 1

inputs_2d.txt

[WeiqunZhang]

The issue is mlmg.setAlwaysUseBNorm(always_use_bnorm), where always_use_bnorm is true. Here $b$ means the right hand side. For a linear system $L(phi) = b$, the residual is $r = b - L(phi)$. The iterative solver stops when the max norm of $r$ is below a certain target. If always_use_bnorm is true, the target is set to relative tolerance times the max norm of $b$. The issue is $b$ is very small in this case. The initial residual is about 2e11, but the inital rhs is at 1.e-7. The solver is able to reduce the residual to 0.02 in 6 iterations, but that residual is still much higher than the rhs.

If I make the following change,

diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H
index eaeadfbb2..9476bbf0c 100644
--- a/Source/ablastr/fields/PoissonSolver.H
+++ b/Source/ablastr/fields/PoissonSolver.H
@@ -253,7 +253,7 @@ computePhi (amrex::Vector<amrex::MultiFab*> const & rho,
         amrex::MLMG mlmg(linop); // actual solver defined here
         mlmg.setVerbose(verbosity);
         mlmg.setMaxIter(max_iters);
-        mlmg.setAlwaysUseBNorm(always_use_bnorm);
+        mlmg.setAlwaysUseBNorm(false);

         // Solve Poisson equation at lev
         mlmg.solve( {phi[lev]}, {rho[lev]},

and run with warpx.self_fields_required_precision=1.e-10, the test seems to run fine.

@RemiLehe What do you suggestion?

[archermarx]

Hey, thanks for looking into this! If always_use_bnorm is set to false, then what convergence criterion is used?

Also, in the mean-time before any changes to the code are made, do you have any suggestions for how to address this? Perhaps initializing the particles pseudo-randomly instead of using a griddedlayout might cause the RHS to be large enough to avoid this problem?

[WeiqunZhang]

It will then use the max of the initial residual and the rhs.

How to address this issue is more a question for @RemiLehe. The solver will never be able to reduce the residual from 2.e11 to something below 1.e-7.

[ax3l]

Note if we change this: same logic in Source/ablastr/fields/PoissonSolver.H and Source/ablastr/fields/VectorPoissonSolver.H

[archermarx]

A side note, if I initialize the particles pseudo-randomly, the RHS is big enough that the solver can converge if i set a tolerance of like 1e-10. For variance reduction purposes, it would be nice to be able to use the uniform distribution.

[ax3l]

X-ref: from @roelof-groenewald - introduced in #2410

[ax3l]

Discussed in developer meeting today:

• Probably setting default back to AMReX default (alwaysUseBNorm=false) and adding an ABLASTR runtime option for cases such as #2410

[roelof-groenewald]

@ax3l Wouldn't it be better to leave the default as is now and allow an ABLASTR runtime option to overwrite it to have always_use_bnorm = false, since this seems to be an unusual case?