diff options
| -rw-r--r-- | Docs/source/usage/examples.rst | 13 | ||||
| -rwxr-xr-x | Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py | 248 | ||||
| -rw-r--r-- | Examples/Physics_applications/capacitive_discharge/README.md | 5 | ||||
| -rwxr-xr-x | Examples/Physics_applications/capacitive_discharge/analysis_1d.py | 8 | ||||
| -rwxr-xr-x | Python/pywarpx/_libwarpx.py | 17 | ||||
| -rw-r--r-- | Regression/WarpX-tests.ini | 17 | ||||
| -rw-r--r-- | Source/Python/WarpXWrappers.H | 10 | ||||
| -rw-r--r-- | Source/Python/WarpXWrappers.cpp | 31 |
8 files changed, 346 insertions, 3 deletions
diff --git a/Docs/source/usage/examples.rst b/Docs/source/usage/examples.rst index 88be580fd..c21904161 100644 --- a/Docs/source/usage/examples.rst +++ b/Docs/source/usage/examples.rst @@ -69,13 +69,22 @@ Uniform plasma Capacitive discharge -------------------- -The Monte-Carlo collision (MCC) model can be used to simulate capacitive discharges between parallel plates. The implementation has been tested against the benchmark results from Turner et. al. in `Phys. Plasmas 20, 013507, 2013 <https://aip.scitation.org/doi/abs/10.1063/1.4775084>`_. The figure below shows a comparison of the ion density as calculated in WarpX (in June 2021) compared to the literature results (which can be found `here <https://aip.scitation.org/doi/suppl/10.1063/1.4775084>`__). An input file with parameters for the 1st case is given below except the total simulation steps have been reduced. +The Monte-Carlo collision (MCC) model can be used to simulate electron and ion collisions with a neutral background gas. In particular this can be used to study capacitive discharges between parallel plates. The implementation has been tested against the benchmark results from Turner et al. in `Phys. Plasmas 20, 013507, 2013 <https://aip.scitation.org/doi/abs/10.1063/1.4775084>`_. The figure below shows a comparison of the ion density as calculated in WarpX (in June 2021) compared to the literature results (which can be found `here <https://aip.scitation.org/doi/suppl/10.1063/1.4775084>`__). .. figure:: https://user-images.githubusercontent.com/40245517/123883650-4f614680-d8fe-11eb-9302-92a282191e8f.png :alt: MCC benchmark against Turner et. al. (2013). :width: 80% -* :download:`test case 1 <../../../Examples/Physics_applications/capacitive_discharge/inputs_2d>` +An input file to reproduce the benchmark calculations is linked below. +To run a given case ``-n``, from 1 to 4, execute: + + .. code-block:: bash + + python3 PICMI_inputs_1d.py -n 1 + +Once the simulation completes an output file ``avg_ion_density.npy`` will be created which can be compared to the literature results as in the plot above. Running case 1 on 4 processors takes roughly 20 minutes to complete. + +* :download:`input file <../../../Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py>` .. note:: diff --git a/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py new file mode 100755 index 000000000..c5e92a5b3 --- /dev/null +++ b/Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py @@ -0,0 +1,248 @@ +#!/usr/bin/env python3 +# +# --- Copyright 2021 Modern Electron +# --- Monte-Carlo Collision script to reproduce the benchmark tests from +# --- Turner et al. (2013) - https://doi.org/10.1063/1.4775084 + +import argparse +import sys + +import numpy as np +from pywarpx import callbacks, fields, picmi + +constants = picmi.constants + + +class CapacitiveDischargeExample(object): + '''The following runs a simulation of a parallel plate capacitor seeded + with a plasma in the spacing between the plates. A time varying voltage is + applied across the capacitor. The groups of 4 values below correspond to + the 4 cases simulated by Turner et al. (2013) in their benchmarks of + PIC-MCC codes. + ''' + + gap = 0.067 # m + + freq = 13.56e6 # Hz + voltage = [450.0, 200.0, 150.0, 120.0] # V + + gas_density = [9.64e20, 32.1e20, 96.4e20, 321e20] # m^-3 + gas_temp = 300.0 # K + m_ion = 6.67e-27 # kg + + plasma_density = [2.56e14, 5.12e14, 5.12e14, 3.84e14] # m^-3 + elec_temp = 30000.0 # K + + seed_nppc = 16 * np.array([32, 16, 8, 4]) + + nz = [128, 256, 512, 512] + + dt = 1.0 / (np.array([400, 800, 1600, 3200]) * freq) + + # Total simulation time in seconds + total_time = np.array([1280, 5120, 5120, 15360]) / freq + # Time (in seconds) between diagnostic evaluations + diag_interval = 32 / freq + + def __init__(self, n=0, test=False): + """Get input parameters for the specific case (n) desired.""" + self.test = test + + # Case specific input parameters + self.voltage = f"{self.voltage[n]}*sin(2*pi*{self.freq:.5e}*t)" + + self.gas_density = self.gas_density[n] + self.plasma_density = self.plasma_density[n] + self.seed_nppc = self.seed_nppc[n] + + self.nz = self.nz[n] + + self.dt = self.dt[n] + self.max_steps = int(self.total_time[n] / self.dt) + self.diag_steps = int(self.diag_interval / self.dt) + + if self.test: + self.max_steps = 20 + self.diag_steps = 5 + + self.rho_wrapper = None + self.ion_density_array = np.zeros(self.nz + 1) + + self.setup_run() + + def setup_run(self): + """Setup simulation components.""" + + ####################################################################### + # Set geometry and boundary conditions # + ####################################################################### + + self.grid = picmi.Cartesian1DGrid( + number_of_cells=[self.nz], + warpx_max_grid_size=128, + lower_bound=[0], + upper_bound=[self.gap], + lower_boundary_conditions=['dirichlet'], + upper_boundary_conditions=['dirichlet'], + lower_boundary_conditions_particles=['absorbing'], + upper_boundary_conditions_particles=['absorbing'], + warpx_potential_hi_z=self.voltage, + ) + + ####################################################################### + # Field solver # + ####################################################################### + + self.solver = picmi.ElectrostaticSolver( + grid=self.grid, method='Multigrid', required_precision=1e-12, + warpx_self_fields_verbosity=0 + ) + + ####################################################################### + # Particle types setup # + ####################################################################### + + self.electrons = picmi.Species( + particle_type='electron', name='electrons', + initial_distribution=picmi.UniformDistribution( + density=self.plasma_density, + rms_velocity=[np.sqrt(constants.kb * self.elec_temp / constants.m_e)]*3, + ) + ) + self.ions = picmi.Species( + particle_type='He', name='he_ions', + charge='q_e', mass=self.m_ion, + initial_distribution=picmi.UniformDistribution( + density=self.plasma_density, + rms_velocity=[np.sqrt(constants.kb * self.gas_temp / self.m_ion)]*3, + ) + ) + + ####################################################################### + # Collision initialization # + ####################################################################### + + cross_sec_direc = '../../../../warpx-data/MCC_cross_sections/He/' + mcc_electrons = picmi.MCCCollisions( + name='coll_elec', + species=self.electrons, + background_density=self.gas_density, + background_temperature=self.gas_temp, + background_mass=self.ions.mass, + scattering_processes={ + 'elastic' : { + 'cross_section' : cross_sec_direc+'electron_scattering.dat' + }, + 'excitation1' : { + 'cross_section': cross_sec_direc+'excitation_1.dat', + 'energy' : 19.82 + }, + 'excitation2' : { + 'cross_section': cross_sec_direc+'excitation_2.dat', + 'energy' : 20.61 + }, + 'ionization' : { + 'cross_section' : cross_sec_direc+'ionization.dat', + 'energy' : 24.55, + 'species' : self.ions + }, + } + ) + + mcc_ions = picmi.MCCCollisions( + name='coll_ion', + species=self.ions, + background_density=self.gas_density, + background_temperature=self.gas_temp, + scattering_processes={ + 'elastic' : { + 'cross_section' : cross_sec_direc+'ion_scattering.dat' + }, + 'back' : { + 'cross_section' : cross_sec_direc+'ion_back_scatter.dat' + }, + # 'charge_exchange' : { + # 'cross_section' : cross_sec_direc+'charge_exchange.dat' + # } + } + ) + + ####################################################################### + # Initialize simulation # + ####################################################################### + + self.sim = picmi.Simulation( + solver=self.solver, + time_step_size=self.dt, + max_steps=self.max_steps, + warpx_collisions=[mcc_electrons, mcc_ions], + warpx_load_balance_intervals=self.max_steps//5000, + verbose=self.test + ) + + self.sim.add_species( + self.electrons, + layout = picmi.GriddedLayout( + n_macroparticle_per_cell=[self.seed_nppc], grid=self.grid + ) + ) + self.sim.add_species( + self.ions, + layout = picmi.GriddedLayout( + n_macroparticle_per_cell=[self.seed_nppc], grid=self.grid + ) + ) + + ####################################################################### + # Add diagnostics for the CI test to be happy # + ####################################################################### + + field_diag = picmi.FieldDiagnostic( + name='diag1', + grid=self.grid, + period=0, + data_list=['rho_electrons', 'rho_he_ions'], + write_dir='.', + warpx_file_prefix='Python_background_mcc_1d_plt' + ) + self.sim.add_diagnostic(field_diag) + + def _get_rho_ions(self): + # deposit the ion density in rho_fp + self.sim.extension.depositChargeDensity('he_ions', 0) + + if self.rho_wrapper is None: + self.rho_wrapper = fields.RhoFPWrapper(0, False) + rho_data = self.rho_wrapper[Ellipsis][:,0] + self.ion_density_array += rho_data / constants.q_e / self.diag_steps + + def run_sim(self): + + self.sim.step(self.max_steps - self.diag_steps) + callbacks.installafterstep(self._get_rho_ions) + self.sim.step(self.diag_steps) + + if self.sim.extension.getMyProc() == 0: + np.save('avg_ion_density.npy', self.ion_density_array) + +########################## +# parse input parameters +########################## + +parser = argparse.ArgumentParser() +parser.add_argument( + '-t', '--test', help='toggle whether this script is run as a short CI test', + action='store_true', +) +parser.add_argument( + '-n', help='Test number to run (1 to 4)', required=False, type=int, + default=1 +) +args, left = parser.parse_known_args() +sys.argv = sys.argv[:1]+left + +if args.n < 1 or args.n > 4: + raise AttributeError('Test number must be an integer from 1 to 4.') + +run = CapacitiveDischargeExample(n=args.n-1, test=args.test) +run.run_sim() diff --git a/Examples/Physics_applications/capacitive_discharge/README.md b/Examples/Physics_applications/capacitive_discharge/README.md index 1b0deadd9..f5653b283 100644 --- a/Examples/Physics_applications/capacitive_discharge/README.md +++ b/Examples/Physics_applications/capacitive_discharge/README.md @@ -1,4 +1,7 @@ -The example in this directory is based on the benchmark cases from Turner et +The examples in this directory are based on the benchmark cases from Turner et al. (Phys. Plasmas 20, 013507, 2013). See the 'Examples' section in the documentation for previously computed comparisons between WarpX and the literature results. +The 1D PICMI input file can be used to reproduce the results from Turner et al. +for a given case, N, by executing: + `python3 PICMI_inputs_1d.py -n N` diff --git a/Examples/Physics_applications/capacitive_discharge/analysis_1d.py b/Examples/Physics_applications/capacitive_discharge/analysis_1d.py new file mode 100755 index 000000000..9a014e494 --- /dev/null +++ b/Examples/Physics_applications/capacitive_discharge/analysis_1d.py @@ -0,0 +1,8 @@ +#!/usr/bin/env python3 + +# Copyright 2021 Modern Electron + +import numpy as np + +density_data = np.load( 'avg_ion_density.npy' ) +assert np.isclose(np.mean(density_data), 2.53023e14) diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py index 105bfc37f..2475d8938 100755 --- a/Python/pywarpx/_libwarpx.py +++ b/Python/pywarpx/_libwarpx.py @@ -1048,6 +1048,23 @@ class LibWarpX(): ''' self.libwarpx_so.warpx_clearParticleBoundaryBuffer() + def depositChargeDensity(self, species_name, level): + ''' + + Deposit the specified species' charge density in rho_fp in order to + access that data via pywarpx.fields.RhoFPWrapper() + + Parameters + ---------- + + species_name : the species name that will be deposited. + level : Which AMR level to retrieve scraped particle data from. + + ''' + self.libwarpx_so.warpx_depositChargeDensity( + ctypes.c_char_p(species_name.encode('utf-8')), level + ) + def _get_mesh_field_list(self, warpx_func, level, direction, include_ghosts): """ Generic routine to fetch the list of field data arrays. diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 81a74825d..ea020eaff 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2884,6 +2884,23 @@ doVis = 0 compareParticles = 0 analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis.py +[Python_background_mcc_1d] +buildDir = . +inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_1d.py +runtime_params = +customRunCmd = python3 PICMI_inputs_1d.py --test +dim = 1 +addToCompileString = USE_PYTHON_MAIN=TRUE +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Physics_applications/capacitive_discharge/analysis_1d.py + [particle_absorption] buildDir = . inputFile = Examples/Modules/ParticleBoundaryProcess/inputs_absorption diff --git a/Source/Python/WarpXWrappers.H b/Source/Python/WarpXWrappers.H index 8d2ee4364..bd47a3ebb 100644 --- a/Source/Python/WarpXWrappers.H +++ b/Source/Python/WarpXWrappers.H @@ -110,6 +110,16 @@ extern "C" { void warpx_clearParticleBoundaryBuffer (); + /** + * \brief This function is used to deposit a given species' charge density + * in the rho_fp multifab which can then be accessed from python via + * pywarpx.fields.RhoFPWrapper() + * + * @param[in] species_name specifying the name of the species to deposit + * @param[in] lev mesh refinement level + */ + void warpx_depositChargeDensity (const char* species_name, int lev); + void warpx_ComputeDt (); void warpx_MoveWindow (int step, bool move_j); diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 7350bbbb6..061e80c6a 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -606,6 +606,37 @@ namespace particle_buffers.clearParticles(); } + void warpx_depositChargeDensity (const char* char_species_name, int lev) { + // this function is used to deposit a given species' charge density + // in the rho_fp multifab which can then be accessed from python via + // pywarpx.fields.RhoFPWrapper() + WarpX& warpx = WarpX::GetInstance(); + const auto & mypc = warpx.GetPartContainer(); + const std::string species_name(char_species_name); + auto & myspc = mypc.GetParticleContainerFromName(species_name); + auto * rho_fp = warpx.get_pointer_rho_fp(lev); + + if (rho_fp == nullptr) { + warpx.RecordWarning( + "WarpXWrappers", "rho_fp is not allocated", WarnPriority::low + ); + return; + } + + // reset rho before depositing + rho_fp->setVal(0.); + + for (WarpXParIter pti(myspc, lev); pti.isValid(); ++pti) + { + const long np = pti.numParticles(); + auto& wp = pti.GetAttribs(PIdx::w); + myspc.DepositCharge(pti, wp, nullptr, rho_fp, 0, 0, np, 0, lev, lev); + } +#ifdef WARPX_DIM_RZ + warpx.ApplyInverseVolumeScalingToChargeDensity(rho_fp, lev); +#endif + } + void warpx_ComputeDt () { WarpX& warpx = WarpX::GetInstance(); warpx.ComputeDt (); |