Moving Window & Embedded Boundaries

[wwshunan]

Hi,
I want to simulate the particle beam between the electrodes. The electrode is imported as an embedded boundary from a STL file. The electrodes are very long, so the moving window was used in the simulation. However I found the embedded boundary also moved with the moving window. How can I solve problem? My code is given below:

from pywarpx import picmi, warpx
from scipy.constants import c, e
from wrdis import readDis 
from pywarpx.callbacks import installcallback
from pywarpx.particle_containers import ParticleContainerWrapper
import pywarpx
import matplotlib.pyplot as plt
import numpy as np

data = readDis('part_rfq.dst')
part_num = int(data.shape[0])

# 设置位置
mass = 938.272
freq = 202.56e6
I = 70e-3
q = I / freq
gamma = data[:, 5] / mass + 1
beta = (1 - gamma**-2)**0.5
lm = c / freq
x = data[:, 0]*0.01
y = data[:, 2]*0.01
z = -data[:, 4]*beta*lm/(2*np.pi)
z -= z.max() + 0.0248536
ux = gamma * data[:, 1] * beta * c 
uy = gamma * data[:, 3] * beta * c 
uz = gamma * beta * c  

constants = picmi.constants
number_of_steps = 50    # Number of simulation steps

# 设置网格参数
nx = 64
ny = 64
nz = 64
xmin = -0.02
xmax = 0.02 
ymin = -0.02
ymax = 0.02
zmin = z.min()
zmax = z.max()
z_range = zmax - zmin
zmin -= 0.2 * z_range
zmax += 0.2 * z_range

v0 = np.mean(beta * c)


# 创建网格
grid = picmi.Cartesian3DGrid(
    number_of_cells=[nx, ny, nz],
    lower_bound=[xmin, ymin, zmin],
    upper_bound=[xmax, ymax, zmax],
    lower_boundary_conditions=['neumann', 'neumann', 'neumann'],
    upper_boundary_conditions=['neumann', 'neumann', 'neumann'],
    lower_boundary_conditions_particles=["absorbing", "absorbing", "absorbing"],
    upper_boundary_conditions_particles=["absorbing", "absorbing", "absorbing"],
    moving_window_velocity=[0.0, 0.0, v0],
)

# 导入 STL 文件
eb = picmi.EmbeddedBoundary(
    stl_file='electrode.stl',
    #potential="(67000.0-(((y*y-x*x)>0)*1.0)*134000.0))*t"
    potential=f"67000.0*sin(2*{np.pi*freq}*t)-(((y*y-x*x)>0)*1.0)*134000.0*sin(2*{np.pi*freq}*t)"
    #potential = "2.0 * (y*y > x*x)"
)

# Define the proton species with in-memory particle data
protons = picmi.Species(
    particle_type='proton',  # Uses m_p and e from scipy.constants
    name='protons',
    initial_distribution=picmi.ParticleListDistribution(
        x=x, y=y, z=z,
        ux=ux, uy=uy, uz=uz,  # Convert momentum to velocity (m/s)
        weight=np.ones_like(x) * q / (e * len(x))    # Uniform weight, total charge = e
    )
)


field_diag = picmi.FieldDiagnostic(
    name='diag',
    grid = grid,
    period = 10,
    data_list = ['E', 'part_per_cell', 'phi'],
    warpx_format = 'plotfile',
    write_dir = './results',
    warpx_file_prefix = 'plt')

particle_diag = picmi.ParticleDiagnostic(
    name='diag',
    period=10,
    species=[protons],
    data_list = ['position', 'weighting', 'momentum'],
    warpx_format='plotfile'
)

solver = picmi.ElectrostaticSolver(
    grid=grid, method="Multigrid", required_precision=1e-7
)

# 创建模拟对象
sim = picmi.Simulation(
    warpx_embedded_boundary=eb,
    time_step_size=1/freq/30,    # Adjust time step (in seconds)
    solver=solver,
    max_steps=10
)
sim.add_species(protons, layout=None)
#sim.add_species(protons, layout=picmi.PseudoRandomLayout(n_macroparticles=len(x)))
#sim.add_species(protons, layout=picmi.PseudoRandomLayout(n_macroparticles=number_sim_particles))
sim.add_diagnostic(field_diag)
sim.add_diagnostic(particle_diag)

# 创建粒子容器包装器
proton_container = ParticleContainerWrapper('protons')

# 回调函数
last_time = [0.0]
last_velocity = [v0]
current_zmin = [zmin]
current_zmax = [zmax]

def adjust_moving_window():
    from pywarpx._libwarpx import libwarpx
    current_step = libwarpx.libwarpx_so.get_instance().getistep(0)
    current_time = libwarpx.libwarpx_so.get_instance().gett_new(0)
    
    # 计算时间步长
    dt = current_time - last_time[0] if current_step > 1 else 0.0
    last_time[0] = current_time
    
    # 获取速度和位置
    uz_arrays = proton_container.get_particle_uz(level=0, copy_to_host=True)
    z_arrays = proton_container.get_particle_z(level=0, copy_to_host=True)
    
    if uz_arrays and z_arrays:
        uz_combined = np.concatenate(uz_arrays)
        z_combined = np.concatenate(z_arrays)
        uz_avg = np.mean(uz_combined) if len(uz_combined) > 0 else v0
        z_avg = np.mean(z_combined) if len(z_combined) > 0 else 0
        z_min_bunch = np.min(z_combined)
        z_max_bunch = np.max(z_combined)
        
        # 调整窗口速度
        window_velocity = uz_avg
        z_window_center = (current_zmin[0] + current_zmax[0]) / 2
        z_deviation = z_avg - z_window_center
        if abs(z_deviation) > 0.02 and dt > 0:  # 位置校正
            correction = z_deviation / dt
            window_velocity += correction
        
        # 动态调整网格大小
        z_range = z_max_bunch - z_min_bunch
        buffer = 0.2 * z_range if z_range > 0 else 0.02  # 20% 缓冲或最小 2 cm
        new_zmin = z_min_bunch - buffer
        new_zmax = z_max_bunch + buffer
        
        # 更新网格边界
        if abs(new_zmin - current_zmin[0]) > 0.01 or abs(new_zmax - current_zmax[0]) > 0.01:  # 避免频繁小调整
            current_zmin[0] = new_zmin
            current_zmax[0] = new_zmax
            warpx.__setattr__('prob_lo', [xmin, ymin, new_zmin])
            warpx.__setattr__('prob_hi', [xmax, ymax, new_zmax])
        
        warpx.__setattr__('moving_window_v', window_velocity / c)
        warpx.__setattr__('moving_window_dir', 'z')
        
        # 计算能量
        gamma = np.sqrt(1 + (uz_combined/c)**2)
        energy = (gamma - 1) * mass * c**2 / e  # eV
        avg_energy = np.mean(energy) if len(energy) > 1 else 0
        
        if current_step % 10 == 0:
            print(f"Step {current_step}, Time {current_time:.2e} s, Window velocity: {window_velocity:.2e} m/s, "
                  f"z_avg: {z_avg:.3f} m, z_range: {z_range:.3f} m, Grid: [{new_zmin:.3f}, {new_zmax:.3f}], "
                  f"Avg Energy: {avg_energy:.2e} eV")
            if z_min_bunch < new_zmin + 0.01 or z_max_bunch > new_zmax - 0.01:
                print(f"Warning: Particles near grid boundaries (z_min_bunch: {z_min_bunch:.3f}, z_max_bunch: {z_max_bunch:.3f})")
    else:
        warpx.__setattr__('moving_window_v', last_velocity[0] / c)
        last_velocity[0] = last_velocity[0]

# 回调函数：动态调整移动窗口速度
#def adjust_moving_window():
#    from pywarpx._libwarpx import libwarpx
#    current_step = libwarpx.libwarpx_so.get_instance().getistep(0)
#    current_time = libwarpx.libwarpx_so.get_instance().gett_new(0)
#    
#    # 获取粒子速度（列表形式，需合并）
#    uz_arrays = proton_container.get_particle_uz(level=0, copy_to_host=True)
#    if uz_arrays:  # 检查是否为空
#        uz_combined = np.concatenate(uz_arrays)
#        uz_avg = np.mean(uz_combined) if len(uz_combined) > 0 else v0
#    else:
#        uz_avg = v0  # 若无粒子，保持初始速度
#    
#    warpx.__setattr__('moving_window_v', uz_avg / c)
#    warpx.__setattr__('moving_window_dir', 'z')
#    if current_step != "unknown" and current_step % 10 == 0:  # 每 10 步打印
#        print(f"Step {current_step}, Time {current_time:.2e} s, Window velocity: {uz_avg:.2e} m/s")
#    elif current_step == "unknown":
#        print(f"Step unknown, Window velocity: {uz_avg:.2e} m/s")

# 注册回调
installcallback('afterstep', adjust_moving_window)  # 在每步结束后调整速度
# 运行模拟
sim.step(number_of_steps)

[ax3l]

Thank you for the report!

However I found the embedded boundary also moved with the moving window.

I think we do not have moving window support for embedded boundaries yet. We still need to add support for that.

[wwshunan]

Thank you for your answer. Could you tell me the approach to achieve this function? I'll try to use AI to assist me to modify the code. The things needed to be finished in the code for the embedded boundary in the moving window should be similar to those in the static window. Meanwhile, I would also like to refer to the code files for embedding boundaries. If you could tell me the relevant file names, that would be even better.