The implementation of the combined boundary condition for the Poisson solver in WarpX #6019
Unanswered
ChenCuiPlasma
asked this question in
Q&A
Replies: 3 comments 20 replies
-
|
Thanks for the questions! The multigrid Poisson solver and the embedded boundaries are enabled primarily by AMReX, while the RZ geometry and the electrostatic solver are more like WarpX-centric features. Let's ping some developers from both sides who may be familiar with these features. @WeiqunZhang @atmyers @jlvay @dpgrote @JustinRayAngus @RevathiJambunathan Pinging also @RemiLehe for future reference. |
Beta Was this translation helpful? Give feedback.
3 replies
-
|
The solution that is proposed with a conductor right at the boundary might be the simplest way to go. @WeiqunZhang, @atmyers, do you have any clue of why the multigrid solver is failing to converge in that case? |
Beta Was this translation helpful? Give feedback.
-
|
Probably because the solver treats it as a singular problem because all domain boundaries are set to Neumann. |
Beta Was this translation helpful? Give feedback.
-
|
Dear Dr. Zoni, Dr. Vay, and Dr. Zhang, Thank you so much for your help and efforts in resolving my question. I really appreciate your efforts and attention on this problem. I have attached my newest output after adopting the modification from Dr. Zhang to the thread below. Best, |
Beta Was this translation helpful? Give feedback.
-
|
Could you give this a try? If it works, we can think about how to add this to WarpX. |
Beta Was this translation helpful? Give feedback.
6 replies
-
|
You don't need to use EB, if you apply my changes. |
Beta Was this translation helpful? Give feedback.
-
|
I am also assuming you are compiling with |
Beta Was this translation helpful? Give feedback.
-
|
Thanks, Dr. Zhang. Now I removed the EB, and the new output looks: |
Beta Was this translation helpful? Give feedback.
-
Yes, I am compiling with 'WarpX_DIMS=RZ'. |
Beta Was this translation helpful? Give feedback.
-
|
That's a different issue now. What are your dx and dz? If they are very different, the solver may have trouble converging. If that's the case, try to make the two numbers closer. |
Beta Was this translation helpful? Give feedback.
-
|
Hi Dr. Zhang, Thanks for your suggestions. In my current code, the dx are set to be equal to dz and equals to 1.66*10^-5 m. I was thinking it might be my plasma injection issue, so I modified my plasma injection condition and make sure to inject two equal species with the same conditions (but just charge sign different). However, the MLMG still fails to converge but it now fails to converge in a few steps. |
Beta Was this translation helpful? Give feedback.
11 replies
-
|
Do you have a test case I can try? |
Beta Was this translation helpful? Give feedback.
-
|
Hi Dr. Zhang, please see the attached code as a test case. I also attached the scenarios for the errors and the corresopnding output after the code body. In this test case, using the pywarpx 25.06, if we assign the boundary of the embedded object to be (zmin, zmax);(rmin, rmax)= (-2,0);(0,10), then the MLMG solver can not converge. If we slightly changed the zmax to be 0.00001, then the MLMG solver converge successfully. Please see the line 176-179 of the code ( # PICMI Python Script for WarpX
# Translated from the in-house axisymmetric espic code.
import math
from pywarpx import callbacks, picmi
from scipy.constants import c, e, epsilon_0, m_e, m_p, k
import scipy.constants as constants
import numpy as np
from pywarpx import picmi, particle_containers, libwarpx
def inertial_length(n: float):
return constants.c / plasma_freq(n)
def J_to_eV(energy: float):
return energy * 6.242e18
def plasma_freq(n: float):
return np.sqrt(n * constants.e**2 / (constants.m_e * constants.epsilon_0))
def cyclotron_freq(m: float, B: float):
return constants.e * B / m
def thermal_velocity(T: float, m: float):
T_j = T / 6.242e18 # temperature in joules
return np.sqrt(2 * T_j / m)
def ion_sound_velocity(T_e: float, T_i: float, m_i: float, Z: int = 1):
T_ij = T_i / 6.242e18 # ion temperature in joules
T_ej = T_e / 6.242e18 # electron temperature in joules
return np.sqrt((Z * T_ej + T_ij) / m_i)
def debye_length(T: float, n: float):
return 7430 * np.sqrt(T / n)
def larmor_radius(m: float, v_perp: float, q: float, B: float):
return m * v_perp / q / B
class FluxMaxwellian_ZInjector(object):
def __init__(
self,
species: picmi.Species,
T: float,
weight: float,
nparts: int,
zmin: float,
zmax: float,
rmin: float,
rmax: float,
):
self.species = species
self.T = T
self.weight = weight
self.nparts = nparts
self.sigma = thermal_velocity(self.T, self.species.mass)
self.zmin = zmin
self.zmax = zmax
self.rmin = rmin
self.rmax = rmax
def inject_parts(self):
nprocs = libwarpx.amr.ParallelDescriptor.NProcs()
nparts_per_proc = int(self.nparts / nprocs)
r = self.rmax * np.sqrt(np.random.rand(nparts_per_proc))
theta = 2 * np.pi * np.random.rand(nparts_per_proc)
x_pos = r * np.cos(theta)
y_pos = r * np.sin(theta)
z_pos = np.random.uniform(self.zmin, self.zmax, nparts_per_proc)
vx_vals = np.random.normal(0, self.sigma, nparts_per_proc) / np.sqrt(2)
vy_vals = np.random.normal(0, self.sigma, nparts_per_proc) / np.sqrt(2)
vz_vals = self.sigma * np.sqrt(
-2.0 * np.log(1.0 - np.random.rand(nparts_per_proc))
)
part_wrapper = particle_containers.ParticleContainerWrapper(self.species.name)
part_wrapper.add_particles(
x=x_pos, y=y_pos, z=z_pos, ux=vx_vals, uy=vy_vals, uz=vz_vals, w=self.weight
)
# =============================================================================
# 1. Conversion rule from the normalized units to SI units
# =============================================================================
# Nomraliing Reference Constants
n0 = 1.0e18 # m^-3
Te0 = 5 # eV
Ti0 = 0.5 # eV
e0 = abs(e) # Elementary charge in C
m0 = m_e # Electron mass in kg
eps0 = epsilon_0 # vacuum permittivity in F/m
k_B = k # Boltzmann constant in J/K
T_e0_K = Te0 * e / k_B
lambda_D0 = math.sqrt((eps0 * k_B * T_e0_K / (n0 * e0**2))) # debye length
omega_pe0 = math.sqrt((n0 * e0**2 / (eps0 * m0)))
v_te0 = math.sqrt(k_B * T_e0_K / m0)
C_s0_proton = math.sqrt((k_B * T_e0_K) / m_p) # Sound speed in m/s
print(f"--- Conversion Factors ---")
print(f"Reference Density n_e0 = {n0:.8e} m^-3")
print(f"Reference Temperature T_e0 = {Te0:.8f} eV")
print(f"Debye Length lambda_D0 = {lambda_D0:.8e} m")
print(f"Plasma Frequency omega_pe0 = {omega_pe0:.8e} rad/s")
print(f"Thermal Velocity v_te0 = {v_te0:.8e} m/s")
print(f"Proton Acoustic C_s0 = {C_s0_proton:.8e} m/s")
print(f"-----------------------------")
# =============================================================================
# 2. Normalized Simulation Parameters (from espix .in files)
# =============================================================================
# Here are the unitless values from the input files.
# From specProp.in
SPEC_NUM = 2
SPEC_1_MASS_NORM = 1.0 # electron mass normalized to itself
SPEC_2_MASS_NORM = m_p / m_e # proton mass normalized to electron mass
SPEC_1_CHARGE_NORM = -1.0 # electron charge normalized to itself
SPEC_2_CHARGE_NORM = 1.0 # proton charge normalized to electron charge
CS0_NORM = (
C_s0_proton / math.sqrt(SPEC_2_MASS_NORM / (m_p / m_e)) / v_te0
) # Sound speed normalized to thermal velocity
# From simulationCtrl.in
DT_NORM = 0.1
END_STEP = 100000
DIAG_EVERY_STEP = 1000
# # domain decomposition (WarpX wants one entry per coordinate)
# MPI_Z_NUM = 4
# MPI_R_NUM = 1
# From meshGeom.in
DOMAIN_Z_MIN_NORM = 0.0
DOMAIN_Z_MAX_NORM = 200.0
DOMAIN_R_MIN_NORM = 0.0
DOMAIN_R_MAX_NORM = 200.0
NZ = 200
NR = 200
FIELD_BOUNDARY_CONDITION = ["neumann", "neumann", "none", "neumann"]
PART_SPEC_1_BOUNDARY_CONDITION = [
"Reflecting",
"Reflecting",
"none",
"Reflecting",
]
PART_SPEC_2_BOUNDARY_CONDITION = ["Absorbing", "Absorbing", "none", "Absorbing"]
# From beams_in_domain.in
BEAM_RADIUS_NORM = 5.0
BEAM_CENTER_Z_NORM = 0.0
BEAM_CENTER_R_NORM = 5.0
BEAM_DENSITY_NORM = 1.0
BEAM_ELEC_TEMP_NORM = 1.0
BEAM_ION_TEMP_NORM = 0.025
MACH_NUM = 4
# BEAM_DRIFT_VEL_NORM = 0.0233370481
BEAM_DRIFT_VEL_NORM = MACH_NUM * CS0_NORM
# From objectsInDomain.in
OBJECTS_NUM = 1
OBJECTS_POTENTIAL_FIX = True
OBJECTS_POTENTIAL_NORM = 0.0 # Normalized potential for objects
OBJECTS_Z_MIN_NORM = -2.0
OBJECTS_Z_MAX_NORM = 0.00000
OBJECTS_R_MIN_NORM = 0.0
OBJECTS_R_MAX_NORM = 10.0
# From field_solver.in
FIELD_SOLVER_TYPE = "Multigrid"
FIELD_SOLVER_CONVERGENCE = 1e-6 # For iterative solvers, the convergence criteria
FIELD_SOLVER_MAX_ITER = 200 # Maximum number of iterations for iterative solvers
FIELD_SOLVER_NEWTON_CONVERGENCE = 1e-8 # Convergence criteria for Newton's method
FIELD_SOLVER_MAX_NEWTON_ITER = 40 # Maximum number of iterations for Newton's method
# From weight.part, not sure, may need adjustment
NPPC_E = 10
NPPC_I = 10
# =============================================================================
# 3. Conversion into SI Simulation Parameters for WarpX
# =============================================================================
# Here we apply the conversion factors to the normalized values.
# Time step: t = t_tilde / omega_pe0
DT_SI = DT_NORM / omega_pe0
# Domain size: x = x_tilde * lambda_D0
DOMAIN_Z_MIN_SI = DOMAIN_Z_MIN_NORM * lambda_D0 # Effectively Z
DOMAIN_R_MIN_SI = DOMAIN_R_MIN_NORM * lambda_D0
DOMAIN_Z_MAX_SI = DOMAIN_Z_MAX_NORM * lambda_D0 # Effectively R
DOMAIN_R_MAX_SI = DOMAIN_R_MAX_NORM * lambda_D0
# Beam properties
# Radius: x = x_tilde * lambda_D0
BEAM_RADIUS_SI = BEAM_RADIUS_NORM * lambda_D0
BEAM_CENTER_Z_SI = BEAM_CENTER_Z_NORM * lambda_D0
BEAM_CENTER_R_SI = BEAM_CENTER_R_NORM * lambda_D0
# Density: n = n_tilde * n_e0
BEAM_DENSITY_SI = BEAM_DENSITY_NORM * n0
# Temperature is T_tilde * T_e0. Since T_e0 = 1 eV, T_norm = T_eV.
BEAM_ELEC_TEMP_SI = BEAM_ELEC_TEMP_NORM * Te0 * e0 / k_B
BEAM_ION_TEMP_SI = BEAM_ION_TEMP_NORM * Te0 * e0 / k_B
# Velocity: v = v_tilde * v_te0
BEAM_DRIFT_VEL_SI = BEAM_DRIFT_VEL_NORM * v_te0
# Species mass and charge are already in SI
SPEC_1_MASS_SI = SPEC_1_MASS_NORM * m0
SPEC_2_MASS_SI = SPEC_2_MASS_NORM * m0
SPEC_1_CHARGE_SI = SPEC_1_CHARGE_NORM * e0
SPEC_2_CHARGE_SI = SPEC_2_CHARGE_NORM * e0
# Objects information
OBJECTS_Z_MIN_SI = OBJECTS_Z_MIN_NORM * lambda_D0
OBJECTS_Z_MAX_SI = OBJECTS_Z_MAX_NORM * lambda_D0
OBJECTS_R_MIN_SI = OBJECTS_R_MIN_NORM * lambda_D0
OBJECTS_R_MAX_SI = OBJECTS_R_MAX_NORM * lambda_D0
OBJECTS_POTENTIAL_SI = OBJECTS_POTENTIAL_NORM / (Te0 / e0)
# =============================================================================
# 4. Grid and Solver (using SI values)
# =============================================================================
grid = picmi.CylindricalGrid(
nr=NR, # Number of radial cells
nz=NZ, # Number of axial cells
n_azimuthal_modes=1, # truly axisymmetric
rmin=DOMAIN_R_MIN_SI,
rmax=DOMAIN_R_MAX_SI,
zmin=DOMAIN_Z_MIN_SI,
zmax=DOMAIN_Z_MAX_SI,
bc_rmin=FIELD_BOUNDARY_CONDITION[2],
bc_rmax=FIELD_BOUNDARY_CONDITION[3],
bc_zmin=FIELD_BOUNDARY_CONDITION[0],
bc_zmax=FIELD_BOUNDARY_CONDITION[1],
# Here use the ion particle boundary conditions as the global boundary conditions to simplify the setup.
# We will override these for each species later.
bc_rmin_particles=PART_SPEC_2_BOUNDARY_CONDITION[2],
bc_rmax_particles=PART_SPEC_2_BOUNDARY_CONDITION[3],
bc_zmin_particles=PART_SPEC_2_BOUNDARY_CONDITION[0],
bc_zmax_particles=PART_SPEC_2_BOUNDARY_CONDITION[1],
)
electrostatic_solver = picmi.ElectrostaticSolver(
grid=grid,
method=FIELD_SOLVER_TYPE,
required_precision=FIELD_SOLVER_CONVERGENCE,
maximum_iterations=FIELD_SOLVER_MAX_ITER,
)
# =============================================================================
# 5. Embedded boundary: cylindrical conductor z belongs to [−2,0] m, r<=r_cond
# =============================================================================
# Conditions defining the solid (conductor) region:
# r >= OBJECTS_R_MIN_SI
# r <= OBJECTS_R_MAX_SI
# z >= OBJECTS_Z_MIN_SI
# z <= OBJECTS_Z_MAX_SI
# Inside conductor = all four ≥ 0
impl_fn = (
f"min( min(x - {OBJECTS_R_MIN_SI}, {OBJECTS_R_MAX_SI} - x), "
f"min(z - {OBJECTS_Z_MIN_SI}, {OBJECTS_Z_MAX_SI} - z) )"
)
conductor_pot = f"{OBJECTS_POTENTIAL_SI}" # Potential of the conductor in SI units
emitter_eb = picmi.EmbeddedBoundary(
implicit_function=impl_fn,
potential=conductor_pot, # fixed potential relative to reference
)
# =============================================================================
# 6. Species Definition and Injection Setup (using SI Uints)
# =============================================================================
# elec_part_bc = picmi.ParticleBoundaryConditions(
# x_lower_boundary=PART_SPEC_1_BOUNDARY_CONDITION[0],
# x_upper_boundary=PART_SPEC_1_BOUNDARY_CONDITION[1],
# y_lower_boundary=PART_SPEC_1_BOUNDARY_CONDITION[2],
# y_upper_boundary=PART_SPEC_1_BOUNDARY_CONDITION[3],
# )
# ion_part_bc = picmi.ParticleBoundaryConditions(
# x_lower_boundary=PART_SPEC_2_BOUNDARY_CONDITION[0],
# x_upper_boundary=PART_SPEC_2_BOUNDARY_CONDITION[1],
# y_lower_boundary=PART_SPEC_2_BOUNDARY_CONDITION[2],
# y_upper_boundary=PART_SPEC_2_BOUNDARY_CONDITION[3],
# )
_face_keys = ["r_lo", "r_hi", "z_lo", "z_hi"]
_reflect_kw = {
"r_lo": "warpx_reflection_model_xlo",
"r_hi": "warpx_reflection_model_xhi",
"z_lo": "warpx_reflection_model_zlo",
"z_hi": "warpx_reflection_model_zhi",
}
_grid_default = {
"r_lo": "none",
"r_hi": "Absorb",
"z_lo": "Absorb",
"z_hi": "Absorb",
}
def build_species_kwargs(spec_bc):
"""Return a dict of reflection-model overrides for this species.
spec_bc is a 4-element list [r_lo, r_hi, z_lo, z_hi]."""
kw = {}
for idx, bc in enumerate(spec_bc):
face = _face_keys[idx]
grid_bc = _grid_default[face]
if bc == grid_bc:
continue # same as grid → no override
# Translate desired bc into reflection probability
if bc in {"Reflect", "Reflecting"}:
prob = "1.0" # always reflect
elif bc in {"Absorb", "Absorbing"}:
prob = "0.0" # never reflect -> absorb
elif bc in {"None", "none"}:
prob = "1.0"
else:
raise ValueError(f"Unsupported BC '{bc}' for face {face}")
kw[_reflect_kw[face]] = prob
return kw
electrons = picmi.Species(
name="electrons",
charge=SPEC_1_CHARGE_SI,
mass=SPEC_1_MASS_SI,
# warpx_reflection_model_zlo="1.0",
# particle_boundary_conditions=elec_part_bc,
# **_spec_1_kwargs,
warpx_reflection_model_zlo="0.0", # Reflecting at z_lo
warpx_reflection_model_zhi="0.0", # Absorbing at z_hi
warpx_reflection_model_xhi="0.0", # Reflecting at r_lo
)
ions = picmi.Species(
name="protons",
charge=SPEC_2_CHARGE_SI,
mass=SPEC_2_MASS_SI,
# particle_boundary_conditions=ion_part_bc,
# **_spec_2_kwargs,
warpx_reflection_model_zlo="0.0", # Reflecting at z_lo
warpx_reflection_model_zhi="0.0", # Absorbing at z_hi
warpx_reflection_model_xhi="0.0", # Reflecting at r_lo
)
electron_injector = FluxMaxwellian_ZInjector(
species=electrons,
T=Te0, # eV
weight=1.0, # or match to your needs
nparts=10000, # total number to inject per timestep (will be split across procs)
zmin=OBJECTS_Z_MAX_SI,
zmax=OBJECTS_Z_MAX_SI, # small layer above the conductor
rmin=OBJECTS_R_MIN_SI,
rmax=OBJECTS_R_MAX_SI,
)
ion_injector = FluxMaxwellian_ZInjector(
species=ions,
T=Ti0,
weight=1.0,
nparts=10000,
zmin=OBJECTS_Z_MAX_SI,
zmax=OBJECTS_Z_MAX_SI, # small layer above the conductor
rmin=OBJECTS_R_MIN_SI,
rmax=OBJECTS_R_MAX_SI,
)
# Layout:
layout_e = picmi.PseudoRandomLayout(n_macroparticles_per_cell=NPPC_E)
layout_i = picmi.PseudoRandomLayout(n_macroparticles_per_cell=NPPC_I)
# =============================================================================
# 7. Diagnostics
# =============================================================================
field_diag = picmi.FieldDiagnostic(
name="diag1",
grid=grid,
period=DIAG_EVERY_STEP,
data_list=["Er", "Ez", "phi", "rho", "rho_electrons", "rho_protons"],
warpx_format="plotfile",
)
part_diag = picmi.ParticleDiagnostic(
name="diag2",
period=DIAG_EVERY_STEP,
species=[electrons, ions],
warpx_format="plotfile",
)
reduced_diag = picmi.ReducedDiagnostic(
name="ParticleNumber", period=DIAG_EVERY_STEP, diag_type="ParticleNumber"
)
# =============================================================================
# 8. Simulation Setup
# =============================================================================
simulation = picmi.Simulation(
solver=electrostatic_solver,
time_step_size=DT_SI,
max_steps=END_STEP,
# warpx_embedded_boundary=emitter_eb,
verbose=1,
)
simulation.add_species(electrons, layout=layout_e)
simulation.add_species(ions, layout=layout_i)
simulation.add_diagnostic(field_diag)
simulation.add_diagnostic(part_diag)
simulation.add_diagnostic(reduced_diag)
simulation.initialize_inputs()
simulation.initialize_warpx()
# =============================================================================
# 9. Run the simulation
# =============================================================================
def inject_both():
electron_injector.inject_parts()
ion_injector.inject_parts()
callbacks.installparticleinjection(inject_both)
# then just run
simulation.step(END_STEP)If we assign the zmax to be just 0, , with the line If we changed the OBJECTS_Z_MAX_NORM to be 0.00001, with the line Now we changed the code based on what Dr. Zhang suggested, and recompile the WarpX code (Version 25.07) using MacOS clang 19.1, then we comment out the line Thank you very much! Best, |
Beta Was this translation helpful? Give feedback.
-
|
I have tried to look into this, but I am not familiar with the python stuff. I can reproduce your failure. But I have no ideas what you are talking about. For example, what the following means.
Unless you can provide C++ approach, I don't really know how to help. |
Beta Was this translation helpful? Give feedback.
-
|
@WeiqunZhang If you add Write input file that can be used to run with the compiled versionsimulation.write_input_file(file_name="inputs_picmi") before simulation.step(END_STEP) then you can look at the WarpX input script in inputs_picmi. |
Beta Was this translation helpful? Give feedback.
-
|
But it has some python callback functions. |
Beta Was this translation helpful? Give feedback.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Uh oh!
There was an error while loading. Please reload this page.
Uh oh!
There was an error while loading. Please reload this page.
-
Dear WarpX team,
I hope this message finds you well.
I’m modeling an electric-propulsion thruster plume with the axisymmetric R-Z electrostatic solver (WarpX 25.06, PICMI / pywarpx).
The physical setup requires a piece-wise boundary (i.e., a combined bonudary condition) at z = 0:
(e.g. Rmax = 200 mm, rc = 10 mm).
In our previous simulation using our in-house PIC code, we imposed this by:
When I mimic this in WarpX with an embedded conductor at the same location, the multigrid Poisson solver fails to converge; it appears the solver still treats the entire z = 0 plane as Neumann. Shifting the conductor slightly inside the domain restores convergence, but that distorts the physical geometry.
If it is possible, may I ask a few questions? Thanks very much in advance.
Questions
Thank you very much for your time and efforts in helping me with this issue and for developing WarpX.
Best regards,
Chen Cui
Beta Was this translation helpful? Give feedback.
All reactions