From 7d26194e42dde770dda3d870d6e9e7f500b05d57 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 29 Feb 2024 15:51:05 -0800 Subject: [PATCH] Add test for particle boundary interaction --- .../PICMI_inputs_rz.py | 165 ++++++++++++++++++ .../particle_boundary_interaction/analysis.py | 51 ++++++ Regression/WarpX-tests.ini | 21 +++ 3 files changed, 237 insertions(+) create mode 100644 Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py create mode 100644 Examples/Tests/particle_boundary_interaction/analysis.py diff --git a/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py b/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py new file mode 100644 index 00000000000..aced3fe3a87 --- /dev/null +++ b/Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python3 +# +# --- Input file for particle-boundary interaction testing in RZ. +# --- This input is a simple case of reflection +# --- of one electron on the surface of a sphere. +import numpy as np + +from pywarpx import callbacks, particle_containers, picmi + +########################## +# numerics parameters +########################## + +dt = 1.0e-11 + +# --- Nb time steps + +max_steps = 23 +diagnostic_interval = 1 + +# --- grid + +nr = 64 +nz= 64 + +rmin = 0.0 +rmax = 2 +zmin = -2 +zmax = 2 + +########################## +# numerics components +########################## + +grid = picmi.CylindricalGrid( + number_of_cells = [nr, nz], + n_azimuthal_modes = 1, + lower_bound = [rmin, zmin], + upper_bound = [rmax, zmax], + lower_boundary_conditions = ['none', 'dirichlet'], + upper_boundary_conditions = ['dirichlet', 'dirichlet'], + lower_boundary_conditions_particles = ['absorbing', 'reflecting'], + upper_boundary_conditions_particles = ['absorbing', 'reflecting'] +) + + +solver = picmi.ElectrostaticSolver( + grid=grid, method='Multigrid', + warpx_absolute_tolerance=1e-7 +) + +embedded_boundary = picmi.EmbeddedBoundary( + implicit_function="-(x**2+y**2+z**2-radius**2)", + radius = 0.2 +) + +########################## +# physics components +########################## + +#one particle +e_dist=picmi.ParticleListDistribution(x=0.0, y=0.0, z=-0.25, ux=0.5e10, uy=0.0, uz=1.0e10, weight=1) + +electrons = picmi.Species( + particle_type='electron', name='electrons', initial_distribution=e_dist, warpx_save_particles_at_eb=1 +) + +########################## +# diagnostics +########################## + +field_diag = picmi.FieldDiagnostic( + name = 'diag1', + grid = grid, + period = diagnostic_interval, + data_list = ['Er', 'Ez', 'phi', 'rho','rho_electrons'], + warpx_format = 'openpmd', + write_dir = '.', + warpx_file_prefix = 'particle_boundary_interaction_plt' +) + +part_diag = picmi.ParticleDiagnostic(name = 'diag1', + period = diagnostic_interval, + species = [electrons], + warpx_format = 'openpmd', + write_dir = '.', + warpx_file_prefix = 'particle_boundary_interaction_plt' +) + +########################## +# simulation setup +########################## + +sim = picmi.Simulation( + solver=solver, + time_step_size = dt, + max_steps = max_steps, + warpx_embedded_boundary=embedded_boundary, + warpx_amrex_the_arena_is_managed=1, +) + +sim.add_species( + electrons, + layout = picmi.GriddedLayout( + n_macroparticle_per_cell=[10, 1, 1], grid=grid + ) +) +sim.add_diagnostic(part_diag) +sim.add_diagnostic(field_diag) + +sim.initialize_inputs() +sim.initialize_warpx() + +########################## +# python particle data access +########################## + +def mirror_reflection(): + buffer = particle_containers.ParticleBoundaryBufferWrapper() #boundary buffer + scraping=(len(buffer.get_particle_boundary_buffer("electrons", 'eb', 'deltaTimeScraped', 0))>0) #bouleen saying if particles where scraped + + if scraping: #otherwise np.concatenate doesnt work + delta_t = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'deltaTimeScraped', 0)) + + + #STEP 1: extract the different parameters of the boundary buffer (normal, time, position) + r = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'x', 0)) + theta = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'theta', 0)) + z = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'z', 0)) + x= r*np.cos(theta) #from RZ coordinates to 3D coordinates + y= r*np.sin(theta) + ux = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'ux', 0)) + uy = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'uy', 0)) + uz = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'uz', 0)) + w = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'w', 0)) + nx = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'nx', 0)) + ny = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'ny', 0)) + nz = np.concatenate(buffer.get_particle_boundary_buffer("electrons", 'eb', 'nz', 0)) + + + #STEP 2: use these parameters to inject particle from the same position in the plasma + elect_pc = particle_containers.ParticleContainerWrapper('electrons') #general particle container + + ####this part is specific to the case of simple reflection. + un=ux*nx+uy*ny+uz*nz + ux_reflect=-2*un*nx+ux #for a "mirror reflection" u(sym)=-2(u.n)n+u + uy_reflect=-2*un*ny+uy + uz_reflect=-2*un*nz+uz + elect_pc.add_particles( + x=x + (dt-delta_t)*ux_reflect, y=y + (dt-delta_t)*uy_reflect, z=z + (dt-delta_t)*uz_reflect, + ux=ux_reflect, uy=uy_reflect, uz=uz_reflect, + w=w + ) #adds the particle in the general particle container at the next step + #### Can be modified depending to the model of interaction. + + buffer.clear_buffer() #reinitialise the boundary buffer + +callbacks.installafterstep(mirror_reflection) #mirror_reflection is called at the next step + # using the new particle container modified at the last step + +########################## +# simulation run +########################## + +sim.step(max_steps) #the whole process is done "max_steps" times diff --git a/Examples/Tests/particle_boundary_interaction/analysis.py b/Examples/Tests/particle_boundary_interaction/analysis.py new file mode 100644 index 00000000000..f4cebbedfb9 --- /dev/null +++ b/Examples/Tests/particle_boundary_interaction/analysis.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python + +""" +This script tests the last coordinate after adding an electron. +The sphere is centered on O and has a radius of 0.2 (EB) +The electron is initially at: (0,0,-0.25) and moves with a velocity: +(0.5e10,0,1.0e10) with a time step of 1e-11. +An input file PICMI_inputs_rz.py is used. +""" +import os +import sys + +import numpy as np +from openpmd_viewer import OpenPMDTimeSeries +import yt + +yt.funcs.mylog.setLevel(0) +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + +# Open plotfile specified in command line +filename = sys.argv[1] +test_name = os.path.split(os.getcwd())[1] +checksumAPI.evaluate_checksum(test_name, filename, output_format='openpmd') + +ts = OpenPMDTimeSeries('./diags/diag1/') + +it=ts.iterations +x,y,z=ts.get_particle(['x','y','z'], species='electrons', iteration=it[-1]) + +# Analytical results calculated +x_analytic=0.03532 +y_analytic=0.00000 +z_analytic=-0.20531 + +print('NUMERICAL coordinates of the point of contact:') +print('x=%5.5f, y=%5.5f, z=%5.5f' % (x[0], y[0], z[0])) +print('\n') +print('ANALYTICAL coordinates of the point of contact:') +print('x=%5.5f, y=%5.5f, z=%5.5f' % (x_analytic, y_analytic, z_analytic)) + +tolerance=1e-5 + +diff_x=np.abs((x[0]-x_analytic)/x_analytic) +diff_z=np.abs((z[0]-z_analytic)/z_analytic) + +print('\n') +print("percentage error for x = %5.4f %%" %(diff_x *100)) +print("percentage error for z = %5.4f %%" %(diff_z *100)) + +assert (diff_x < tolerance) and (y[0] < 1e-8) and (diff_z < tolerance), 'Test particle_boundary_interaction did not pass' diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 128495195da..9c0437de5d3 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -4664,6 +4664,27 @@ doVis = 0 compareParticles = 1 analysisRoutine = Examples/Tests/energy_conserving_thermal_plasma/analysis.py +[particle_boundary_interaction] +buildDir = . +inputFile = Examples/Tests/particle_boundary_interaction/PICMI_inputs_rz.py +runtime_params = +customRunCmd = python3 PICMI_inputs_rz.py +dim = 2 +addToCompileString = USE_PYTHON_MAIN=TRUE USE_RZ=TRUE +cmakeSetupOpts = -DWarpX_DIMS="RZ" -DWarpX_EB=ON -DWarpX_PYTHON=ON +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons +outputFile = particle_boundary_interaction_plt +analysisRoutine = Examples/Tests/particle_boundary_interaction/analysis.py + [focusing_gaussian_beam] buildDir = . inputFile = Examples/Tests/gaussian_beam/inputs_focusing_beam