|
| 1 | +#!/usr/bin/env python3 |
| 2 | + |
| 3 | +from pywarpx import picmi |
| 4 | + |
| 5 | +# Physical constants |
| 6 | +c = picmi.constants.c |
| 7 | +q_e = picmi.constants.q_e |
| 8 | + |
| 9 | +# We only run 100 steps for tests |
| 10 | +# Disable `max_step` below to run until the physical `stop_time`. |
| 11 | +max_step = 100 |
| 12 | +# time-scale with highly kinetic dynamics |
| 13 | +stop_time = 0.2e-12 |
| 14 | + |
| 15 | +# proper resolution for 30 n_c (dx<=3.33nm) incl. acc. length |
| 16 | +# (>=6x V100) |
| 17 | +# --> choose larger `max_grid_size` and `blocking_factor` for 1 to 8 grids per GPU accordingly |
| 18 | +#nx = 7488 |
| 19 | +#nz = 14720 |
| 20 | + |
| 21 | +# Number of cells |
| 22 | +nx = 384 |
| 23 | +nz = 512 |
| 24 | + |
| 25 | +# Domain decomposition (deactivate `warpx_numprocs` in `picmi.Simulation` for this to take effect) |
| 26 | +max_grid_size = 64 |
| 27 | +blocking_factor = 32 |
| 28 | + |
| 29 | +# Physical domain |
| 30 | +xmin = -7.5e-06 |
| 31 | +xmax = 7.5e-06 |
| 32 | +zmin = -5.0e-06 |
| 33 | +zmax = 25.0e-06 |
| 34 | + |
| 35 | +# Create grid |
| 36 | +grid = picmi.Cartesian2DGrid( |
| 37 | + number_of_cells=[nx, nz], |
| 38 | + lower_bound=[xmin, zmin], |
| 39 | + upper_bound=[xmax, zmax], |
| 40 | + lower_boundary_conditions=['open', 'open'], |
| 41 | + upper_boundary_conditions=['open', 'open'], |
| 42 | + lower_boundary_conditions_particles=['absorbing', 'absorbing'], |
| 43 | + upper_boundary_conditions_particles=['absorbing', 'absorbing'], |
| 44 | + warpx_max_grid_size=max_grid_size, |
| 45 | + warpx_blocking_factor=blocking_factor) |
| 46 | + |
| 47 | +# Particles: plasma parameters |
| 48 | +# critical plasma density |
| 49 | +nc = 1.742e27 # [m^-3] 1.11485e21 * 1.e6 / 0.8**2 |
| 50 | +# number density: "fully ionized" electron density as reference |
| 51 | +# [material 1] cryogenic H2 |
| 52 | +n0 = 30.0 # [n_c] |
| 53 | +# [material 2] liquid crystal |
| 54 | +# n0 = 192 |
| 55 | +# [material 3] PMMA |
| 56 | +# n0 = 230 |
| 57 | +# [material 4] Copper (ion density: 8.49e28/m^3; times ionization level) |
| 58 | +# n0 = 1400 |
| 59 | +plasma_density = n0 * nc |
| 60 | +preplasma_L = 0.05e-6 # [m] scale length (>0) |
| 61 | +preplasma_Lcut = 2.0e-6 # [m] hard cutoff from surface |
| 62 | +plasma_r0 = 2.5e-6 # [m] radius or half-thickness |
| 63 | +plasma_eps_z = 0.05e-6 # [m] small offset in z to make zmin, zmax interval larger than 2*(r0 + Lcut) |
| 64 | +plasma_creation_limit_z = plasma_r0 + preplasma_Lcut + plasma_eps_z # [m] upper limit in z for particle creation |
| 65 | + |
| 66 | +plasma_xmin = None |
| 67 | +plasma_ymin = None |
| 68 | +plasma_zmin = -plasma_creation_limit_z |
| 69 | +plasma_xmax = None |
| 70 | +plasma_ymax = None |
| 71 | +plasma_zmax = plasma_creation_limit_z |
| 72 | + |
| 73 | +density_expression_str = f'{plasma_density}*((abs(z)<={plasma_r0}) + (abs(z)<{plasma_r0}+{preplasma_Lcut}) * (abs(z)>{plasma_r0}) * exp(-(abs(z)-{plasma_r0})/{preplasma_L}))' |
| 74 | + |
| 75 | +slab_with_ramp_dist_hydrogen = picmi.AnalyticDistribution( |
| 76 | + density_expression=density_expression_str, |
| 77 | + lower_bound=[plasma_xmin, plasma_ymin, plasma_zmin], |
| 78 | + upper_bound=[plasma_xmax, plasma_ymax, plasma_zmax] |
| 79 | +) |
| 80 | + |
| 81 | +# thermal velocity spread for electrons in gamma*beta |
| 82 | +ux_th = .01 |
| 83 | +uz_th = .01 |
| 84 | + |
| 85 | +slab_with_ramp_dist_electrons = picmi.AnalyticDistribution( |
| 86 | + density_expression=density_expression_str, |
| 87 | + lower_bound=[plasma_xmin, plasma_ymin, plasma_zmin], |
| 88 | + upper_bound=[plasma_xmax, plasma_ymax, plasma_zmax], |
| 89 | + # if `momentum_expressions` and `momentum_spread_expressions` are unset, |
| 90 | + # a Gaussian momentum distribution is assumed given that `rms_velocity` has any non-zero elements |
| 91 | + rms_velocity=[c*ux_th, 0., c*uz_th] # thermal velocity spread in m/s |
| 92 | +) |
| 93 | + |
| 94 | +# TODO: add additional attributes orig_x and orig_z |
| 95 | +electrons = picmi.Species( |
| 96 | + particle_type='electron', |
| 97 | + name='electrons', |
| 98 | + initial_distribution=slab_with_ramp_dist_electrons, |
| 99 | +) |
| 100 | + |
| 101 | +# TODO: add additional attributes orig_x and orig_z |
| 102 | +hydrogen = picmi.Species( |
| 103 | + particle_type='proton', |
| 104 | + name='hydrogen', |
| 105 | + initial_distribution=slab_with_ramp_dist_hydrogen |
| 106 | +) |
| 107 | + |
| 108 | +# Laser |
| 109 | +# e_max = a0 * 3.211e12 / lambda_0[mu] |
| 110 | +# a0 = 16, lambda_0 = 0.8mu -> e_max = 64.22 TV/m |
| 111 | +e_max = 64.22e12 |
| 112 | +position_z = -4.0e-06 |
| 113 | +profile_t_peak = 50.e-15 |
| 114 | +profile_focal_distance = 4.0e-06 |
| 115 | +laser = picmi.GaussianLaser( |
| 116 | + wavelength=0.8e-06, |
| 117 | + waist=4.e-06, |
| 118 | + duration=30.e-15, |
| 119 | + focal_position=[0, 0, profile_focal_distance + position_z], |
| 120 | + centroid_position=[0, 0, position_z - c * profile_t_peak], |
| 121 | + propagation_direction=[0, 0, 1], |
| 122 | + polarization_direction=[1, 0, 0], |
| 123 | + E0=e_max, |
| 124 | + fill_in=False) |
| 125 | +laser_antenna = picmi.LaserAntenna( |
| 126 | + position=[0., 0., position_z], |
| 127 | + normal_vector=[0, 0, 1]) |
| 128 | + |
| 129 | +# Electromagnetic solver |
| 130 | +solver = picmi.ElectromagneticSolver( |
| 131 | + grid=grid, |
| 132 | + method='Yee', |
| 133 | + cfl=0.999, |
| 134 | + divE_cleaning=0, |
| 135 | + #warpx_pml_ncell=10 |
| 136 | +) |
| 137 | + |
| 138 | +# Diagnostics |
| 139 | +particle_diag = picmi.ParticleDiagnostic( |
| 140 | + name='Python_LaserIonAcc2d_plt', |
| 141 | + period=100, |
| 142 | + write_dir='./diags', |
| 143 | + warpx_format='openpmd', |
| 144 | + warpx_openpmd_backend='h5', |
| 145 | + # demonstration of a spatial and momentum filter |
| 146 | + warpx_plot_filter_function='(uz>=0) * (x<1.0e-6) * (x>-1.0e-6)' |
| 147 | +) |
| 148 | +# reduce resolution of output fields |
| 149 | +coarsening_ratio = [4, 4] |
| 150 | +ncell_field = [] |
| 151 | +for (ncell_comp, cr) in zip([nx,nz], coarsening_ratio): |
| 152 | + ncell_field.append(int(ncell_comp/cr)) |
| 153 | +field_diag = picmi.FieldDiagnostic( |
| 154 | + name='Python_LaserIonAcc2d_plt', |
| 155 | + grid=grid, |
| 156 | + period=100, |
| 157 | + number_of_cells=ncell_field, |
| 158 | + data_list=['B', 'E', 'J', 'rho', 'rho_electrons', 'rho_hydrogen'], |
| 159 | + write_dir='./diags', |
| 160 | + warpx_format='openpmd', |
| 161 | + warpx_openpmd_backend='h5' |
| 162 | +) |
| 163 | + |
| 164 | +particle_fw_diag = picmi.ParticleDiagnostic( |
| 165 | + name='openPMDfw', |
| 166 | + period=100, |
| 167 | + write_dir='./diags', |
| 168 | + warpx_format='openpmd', |
| 169 | + warpx_openpmd_backend='h5', |
| 170 | + warpx_plot_filter_function='(uz>=0) * (x<1.0e-6) * (x>-1.0e-6)' |
| 171 | +) |
| 172 | + |
| 173 | +particle_bw_diag = picmi.ParticleDiagnostic( |
| 174 | + name='openPMDbw', |
| 175 | + period=100, |
| 176 | + write_dir='./diags', |
| 177 | + warpx_format='openpmd', |
| 178 | + warpx_openpmd_backend='h5', |
| 179 | + warpx_plot_filter_function='(uz<0)' |
| 180 | +) |
| 181 | + |
| 182 | +# histograms with 2.0 degree acceptance angle in fw direction |
| 183 | +# 2 deg * pi / 180 : 0.03490658503 rad |
| 184 | +# half-angle +/- : 0.017453292515 rad |
| 185 | +histuH_rdiag = picmi.ReducedDiagnostic( |
| 186 | + diag_type='ParticleHistogram', |
| 187 | + name='histuH', |
| 188 | + period=100, |
| 189 | + species=hydrogen, |
| 190 | + bin_number=1000, |
| 191 | + bin_min=0.0, |
| 192 | + bin_max=0.474, # 100 MeV protons |
| 193 | + histogram_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)', |
| 194 | + filter_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)') |
| 195 | + |
| 196 | +histue_rdiag = picmi.ReducedDiagnostic( |
| 197 | + diag_type='ParticleHistogram', |
| 198 | + name='histue', |
| 199 | + period=100, |
| 200 | + species=electrons, |
| 201 | + bin_number=1000, |
| 202 | + bin_min=0.0, |
| 203 | + bin_max=197.0, # 100 MeV electrons |
| 204 | + histogram_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)', |
| 205 | + filter_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)') |
| 206 | + |
| 207 | +# just a test entry to make sure that the histogram filter is purely optional: |
| 208 | +# this one just records uz of all hydrogen ions, independent of their pointing |
| 209 | +histuzAll_rdiag = picmi.ReducedDiagnostic( |
| 210 | + diag_type='ParticleHistogram', |
| 211 | + name='histuzAll', |
| 212 | + period=100, |
| 213 | + species=hydrogen, |
| 214 | + bin_number=1000, |
| 215 | + bin_min=-0.474, |
| 216 | + bin_max=0.474, |
| 217 | + histogram_function='uz') |
| 218 | + |
| 219 | +field_probe_z_rdiag = picmi.ReducedDiagnostic( |
| 220 | + diag_type='FieldProbe', |
| 221 | + name='FieldProbe_Z', |
| 222 | + period=100, |
| 223 | + integrate=0, |
| 224 | + probe_geometry='Line', |
| 225 | + x_probe=0.0, |
| 226 | + z_probe=-5.0e-6, |
| 227 | + x1_probe=0.0, |
| 228 | + z1_probe=25.0e-6, |
| 229 | + resolution=3712) |
| 230 | + |
| 231 | +field_probe_scat_point_rdiag = picmi.ReducedDiagnostic( |
| 232 | + diag_type='FieldProbe', |
| 233 | + name='FieldProbe_ScatPoint', |
| 234 | + period=1, |
| 235 | + integrate=0, |
| 236 | + probe_geometry='Point', |
| 237 | + x_probe=0.0, |
| 238 | + z_probe=15.0e-6) |
| 239 | + |
| 240 | +field_probe_scat_line_rdiag = picmi.ReducedDiagnostic( |
| 241 | + diag_type='FieldProbe', |
| 242 | + name='FieldProbe_ScatLine', |
| 243 | + period=100, |
| 244 | + integrate=1, |
| 245 | + probe_geometry='Line', |
| 246 | + x_probe=-2.5e-6, |
| 247 | + z_probe=15.0e-6, |
| 248 | + x1_probe=2.5e-6, |
| 249 | + z1_probe=15e-6, |
| 250 | + resolution=201) |
| 251 | + |
| 252 | +load_balance_costs_rdiag = picmi.ReducedDiagnostic( |
| 253 | + diag_type='LoadBalanceCosts', |
| 254 | + name='LBC', |
| 255 | + period=100) |
| 256 | + |
| 257 | +# Set up simulation |
| 258 | +sim = picmi.Simulation( |
| 259 | + solver=solver, |
| 260 | + max_time=stop_time, # need to remove `max_step` to run this far |
| 261 | + verbose=1, |
| 262 | + particle_shape='cubic', |
| 263 | + warpx_numprocs=[1, 2], # deactivate `numprocs` for dynamic load balancing |
| 264 | + warpx_use_filter=1, |
| 265 | + warpx_load_balance_intervals=100, |
| 266 | + warpx_load_balance_costs_update='heuristic' |
| 267 | +) |
| 268 | + |
| 269 | +# Add plasma electrons |
| 270 | +sim.add_species( |
| 271 | + electrons, |
| 272 | + layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[2,2]) |
| 273 | + # for more realistic simulations, try to avoid that macro-particles represent more than 1 n_c |
| 274 | + #layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[4,8]) |
| 275 | +) |
| 276 | + |
| 277 | +# Add hydrogen ions |
| 278 | +sim.add_species( |
| 279 | + hydrogen, |
| 280 | + layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[2,2]) |
| 281 | + # for more realistic simulations, try to avoid that macro-particles represent more than 1 n_c |
| 282 | + #layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[4,8]) |
| 283 | +) |
| 284 | + |
| 285 | +# Add laser |
| 286 | +sim.add_laser( |
| 287 | + laser, |
| 288 | + injection_method=laser_antenna) |
| 289 | + |
| 290 | +# Add full diagnostics |
| 291 | +sim.add_diagnostic(particle_diag) |
| 292 | +sim.add_diagnostic(field_diag) |
| 293 | +sim.add_diagnostic(particle_fw_diag) |
| 294 | +sim.add_diagnostic(particle_bw_diag) |
| 295 | +# Add reduced diagnostics |
| 296 | +sim.add_diagnostic(histuH_rdiag) |
| 297 | +sim.add_diagnostic(histue_rdiag) |
| 298 | +sim.add_diagnostic(histuzAll_rdiag) |
| 299 | +sim.add_diagnostic(field_probe_z_rdiag) |
| 300 | +sim.add_diagnostic(field_probe_scat_point_rdiag) |
| 301 | +sim.add_diagnostic(field_probe_scat_line_rdiag) |
| 302 | +sim.add_diagnostic(load_balance_costs_rdiag) |
| 303 | +# TODO: make ParticleHistogram2D available |
| 304 | + |
| 305 | +# Write input file that can be used to run with the compiled version |
| 306 | +sim.write_input_file(file_name='inputs_2d_picmi') |
| 307 | + |
| 308 | +# Initialize inputs and WarpX instance |
| 309 | +sim.initialize_inputs() |
| 310 | +sim.initialize_warpx() |
| 311 | + |
| 312 | +# Advance simulation until last time step |
| 313 | +sim.step(max_step) |
0 commit comments