diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 6a86e1fe10d..6b21758f6e8 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1747,8 +1747,8 @@ are applied to the grid directly. In particular, these fields can be seen in the One can refer to input files in ``Examples/Tests/LoadExternalField`` for more information. Regarding how to prepare the openPMD data file, one can refer to the `openPMD-example-datasets `__. - Note that if both `B_ext_grid_init_style` and `E_ext_grid_init_style` are set to - `read_from_file`, the openPMD file specified by `warpx.read_fields_from_path` + Note that if both ``B_ext_grid_init_style`` and ``E_ext_grid_init_style`` are set to + ``read_from_file``, the openPMD file specified by ``warpx.read_fields_from_path`` should contain both B and E external fields data. * ``warpx.E_external_grid`` & ``warpx.B_external_grid`` (list of `3 floats`) @@ -1804,6 +1804,21 @@ are applied to the particles directly, at each timestep. As a results, these fie Note that the position is defined in Cartesian coordinates, as a function of (x,y,z), even for RZ. + * ``read_from_file``: load the external field from an openPMD file. + An additional parameter, indicating the path of an openPMD data file, ``particles.read_fields_from_path`` + must be specified, from which the external E field data can be loaded into WarpX. + One can refer to input files in ``Examples/Tests/LoadExternalField`` for more information. + Regarding how to prepare the openPMD data file, one can refer to + the `openPMD-example-datasets `__. + Note that if both ``B_ext_particle_init_style`` and ``E_ext_particle_init_style`` are set to + ``read_from_file``, the openPMD file specified by ``particles.read_fields_from_path`` + should contain both B and E external fields data. + + .. note:: + + When using ``read_from_file``, the fields loaded from the file will be interpolated + to the resolution of the grid used for the simulation. + * ``repeated_plasma_lens``: apply a series of plasma lenses. The properties of the lenses are defined in the lab frame by the input parameters: diff --git a/Examples/Tests/LoadExternalField/PICMI_inputs_3d.py b/Examples/Tests/LoadExternalField/PICMI_inputs_3d_grid_fields.py similarity index 96% rename from Examples/Tests/LoadExternalField/PICMI_inputs_3d.py rename to Examples/Tests/LoadExternalField/PICMI_inputs_3d_grid_fields.py index 849d4d2be08..d128d9c10e0 100644 --- a/Examples/Tests/LoadExternalField/PICMI_inputs_3d.py +++ b/Examples/Tests/LoadExternalField/PICMI_inputs_3d_grid_fields.py @@ -93,7 +93,7 @@ species=[ions], data_list = ['ux', 'uy', 'uz', 'x', 'y', 'z', 'weighting'], write_dir='.', - warpx_file_prefix='Python_LoadExternalField3D_plt' + warpx_file_prefix='Python_LoadExternalGridField3D_plt' ) field_diag = picmi.FieldDiagnostic( name='diag1', @@ -101,7 +101,7 @@ period=300, data_list = ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Jx', 'Jy', 'Jz'], write_dir='.', - warpx_file_prefix='Python_LoadExternalField3D_plt' + warpx_file_prefix='Python_LoadExternalGridField3D_plt' ) ################################# diff --git a/Examples/Tests/LoadExternalField/PICMI_inputs_3d_particle_fields.py b/Examples/Tests/LoadExternalField/PICMI_inputs_3d_particle_fields.py new file mode 100644 index 00000000000..7bf7c5a084c --- /dev/null +++ b/Examples/Tests/LoadExternalField/PICMI_inputs_3d_particle_fields.py @@ -0,0 +1,140 @@ +#!/usr/bin/env python3 +# +# --- Input file for loading initial field from openPMD file. + +from pywarpx import picmi + +constants = picmi.constants + +################################# +####### GENERAL PARAMETERS ###### +################################# + +max_steps = 300 + +max_grid_size = 40 +nx = max_grid_size +ny = max_grid_size +nz = max_grid_size + +xmin = -1 +xmax = 1 +ymin = xmin +ymax = xmax +zmin = 0 +zmax = 5 + +number_per_cell = 200 + +################################# +############ NUMERICS ########### +################################# + +verbose = 1 +dt = 4.4e-7 +use_filter = 0 + +# Order of particle shape factors +particle_shape = 1 + +################################# +############ PLASMA ############# +################################# + +ion_dist = picmi.ParticleListDistribution( + x=0.0, + y=0.2, + z=2.5, + ux=9.5e-05*constants.c, + uy=0.0*constants.c, + uz=1.34e-4*constants.c, + weight=1.0 + ) + +ions = picmi.Species( + particle_type='H', + name='proton', charge='q_e',mass="m_p", + warpx_do_not_deposit=1, + initial_distribution=ion_dist + ) + +################################# +######## INITIAL FIELD ########## +################################# + +applied_field = picmi.LoadAppliedField( + read_fields_from_path="../../../../openPMD-example-datasets/example-femm-3d.h5", + load_E=False + ) + +################################# +###### GRID AND SOLVER ########## +################################# + +grid = picmi.Cartesian3DGrid( + number_of_cells=[nx, ny, nz], + warpx_max_grid_size=max_grid_size, + lower_bound=[xmin, ymin, zmin], + upper_bound=[xmax, ymax, zmax], + lower_boundary_conditions=['dirichlet', 'dirichlet', 'dirichlet'], + upper_boundary_conditions=['dirichlet', 'dirichlet', 'dirichlet'], + lower_boundary_conditions_particles=['absorbing', 'absorbing', 'absorbing'], + upper_boundary_conditions_particles=['absorbing', 'absorbing', 'absorbing'] + ) +solver = picmi.ElectrostaticSolver(grid=grid) + +################################# +######### DIAGNOSTICS ########### +################################# + +particle_diag = picmi.ParticleDiagnostic( + name='diag1', + period=300, + species=[ions], + data_list = ['ux', 'uy', 'uz', 'x', 'y', 'z', 'weighting'], + write_dir='.', + warpx_file_prefix='Python_LoadExternalParticleField3D_plt' + ) +field_diag = picmi.FieldDiagnostic( + name='diag1', + grid=grid, + period=300, + data_list = ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Jx', 'Jy', 'Jz'], + write_dir='.', + warpx_file_prefix='Python_LoadExternalParticleField3D_plt' + ) + +################################# +####### SIMULATION SETUP ######## +################################# + +sim = picmi.Simulation( + solver=solver, + max_steps=max_steps, + verbose=verbose, + warpx_serialize_initial_conditions=False, + warpx_grid_type='collocated', + warpx_do_dynamic_scheduling=False, + warpx_use_filter=use_filter, + time_step_size=dt, + particle_shape=particle_shape + ) + +sim.add_applied_field(applied_field) + +sim.add_species( + ions, + layout=picmi.PseudoRandomLayout( + n_macroparticles_per_cell=number_per_cell, grid=grid + ) + ) + +sim.add_diagnostic(field_diag) +sim.add_diagnostic(particle_diag) + +################################# +##### SIMULATION EXECUTION ###### +################################# + +#sim.write_input_file('PICMI_inputs_3d') +sim.step(max_steps) diff --git a/Examples/Tests/LoadExternalField/inputs_rz b/Examples/Tests/LoadExternalField/inputs_rz_grid_fields similarity index 100% rename from Examples/Tests/LoadExternalField/inputs_rz rename to Examples/Tests/LoadExternalField/inputs_rz_grid_fields diff --git a/Examples/Tests/LoadExternalField/inputs_rz_particle_fields b/Examples/Tests/LoadExternalField/inputs_rz_particle_fields new file mode 100644 index 00000000000..b76d4cb7efc --- /dev/null +++ b/Examples/Tests/LoadExternalField/inputs_rz_particle_fields @@ -0,0 +1,65 @@ +warpx.serialize_initial_conditions = 0 +warpx.do_dynamic_scheduling = 0 +particles.do_tiling = 0 + +particles.B_ext_particle_init_style = "read_from_file" +particles.read_fields_from_path = "../../../../openPMD-example-datasets/example-femm-thetaMode.h5" + +warpx.grid_type = collocated +warpx.do_electrostatic = labframe + +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 300 +amr.n_cell = 40 40 +warpx.numprocs = 1 1 +amr.max_level = 0 +geometry.dims = RZ + +geometry.prob_lo = 0.0 0.0 +geometry.prob_hi = 1.0 5.0 + +################################# +###### Boundary Condition ####### +################################# +boundary.field_lo = none pec +boundary.field_hi = pec pec +boundary.potential_lo_x = 0 +boundary.potential_hi_x = 0 +boundary.potential_lo_y = 0 +boundary.potential_hi_y = 0 +boundary.potential_lo_z = 0 +boundary.potential_hi_z = 0 + +################################# +############ NUMERICS ########### +################################# +warpx.serialize_initial_conditions = 1 +warpx.verbose = 1 +warpx.const_dt = 4.40917904849092e-7 +warpx.use_filter = 0 + +# Order of particle shape factors +algo.particle_shape = 1 + +################################# +############ PLASMA ############# +################################# +particles.species_names = proton +proton.injection_style = "SingleParticle" +proton.single_particle_pos = 0.0 0.2 2.5 +proton.single_particle_u = 9.506735958279367e-05 0.0 0.00013435537232359165 +proton.single_particle_weight = 1.0 +proton.do_not_deposit = 1 +proton.mass = m_p +proton.charge = q_e + +# Diagnostics +diagnostics.diags_names = diag1 chk +diag1.intervals = 300 +diag1.diag_type = Full + +chk.intervals = 150 +chk.diag_type = Full +chk.format = checkpoint diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 06cd58be6ef..d4203e95509 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -1527,6 +1527,13 @@ def applied_field_initialize_inputs(self): expression = pywarpx.my_constants.mangle_expression(expression, self.mangle_dict) pywarpx.warpx.__setattr__(f'B{sdir}_external_grid_function(x,y,z)', expression) +class LoadAppliedField(picmistandard.PICMI_LoadAppliedField): + def applied_field_initialize_inputs(self): + pywarpx.particles.read_fields_from_path = self.read_fields_from_path + if self.load_E: + pywarpx.particles.E_ext_particle_init_style = 'read_from_file' + if self.load_B: + pywarpx.particles.B_ext_particle_init_style = 'read_from_file' class ConstantAppliedField(picmistandard.PICMI_ConstantAppliedField): def applied_field_initialize_inputs(self): diff --git a/Regression/Checksum/benchmarks_json/LoadExternalFieldRZ.json b/Regression/Checksum/benchmarks_json/LoadExternalFieldRZGrid.json similarity index 100% rename from Regression/Checksum/benchmarks_json/LoadExternalFieldRZ.json rename to Regression/Checksum/benchmarks_json/LoadExternalFieldRZGrid.json diff --git a/Regression/Checksum/benchmarks_json/LoadExternalFieldRZParticles.json b/Regression/Checksum/benchmarks_json/LoadExternalFieldRZParticles.json new file mode 100644 index 00000000000..dc4ffa3beb2 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/LoadExternalFieldRZParticles.json @@ -0,0 +1,22 @@ +{ + "lev=0": { + "Br": 0.6380326770459374, + "Bt": 0.0, + "Bz": 4.850698699951099, + "Er": 0.0, + "Et": 0.0, + "Ez": 0.0, + "jr": 0.0, + "jt": 0.0, + "jz": 0.0 + }, + "proton": { + "particle_momentum_x": 1.2945379007125463e-23, + "particle_momentum_y": 7.901629565307941e-23, + "particle_momentum_z": 2.0004592432172922e-23, + "particle_position_x": 0.12402004585603355, + "particle_position_y": 4.363249204658946, + "particle_theta": 0.22200801011337173, + "particle_weight": 1.0 + } +} \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Python_LoadExternalField3D.json b/Regression/Checksum/benchmarks_json/Python_LoadExternalGridField3D.json similarity index 100% rename from Regression/Checksum/benchmarks_json/Python_LoadExternalField3D.json rename to Regression/Checksum/benchmarks_json/Python_LoadExternalGridField3D.json diff --git a/Regression/Checksum/benchmarks_json/Python_LoadExternalParticleField3D.json b/Regression/Checksum/benchmarks_json/Python_LoadExternalParticleField3D.json new file mode 100644 index 00000000000..a34b2f136a6 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/Python_LoadExternalParticleField3D.json @@ -0,0 +1,22 @@ +{ + "lev=0": { + "Bx": 29.8448958859583, + "By": 29.844895885958305, + "Bz": 197.55124310523018, + "Ex": 0.0, + "Ey": 0.0, + "Ez": 0.0, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0 + }, + "proton": { + "particle_momentum_x": 1.2118953253556959e-23, + "particle_momentum_y": 7.822512598009088e-23, + "particle_momentum_z": 2.2761898274347412e-23, + "particle_position_x": 0.12238072371983327, + "particle_position_y": 0.009653941154598592, + "particle_position_z": 4.317544769595657, + "particle_weight": 1.0 + } +} \ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 4972d453f5f..ad498d63945 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2266,10 +2266,29 @@ doVis = 0 compareParticles = 0 analysisRoutine = Examples/Tests/resampling/analysis_leveling_thinning.py -[LoadExternalFieldRZ] +[LoadExternalFieldRZGrid] buildDir = . -inputFile = Examples/Tests/LoadExternalField/inputs_rz -runtime_params = warpx.abort_on_warning_threshold=medium chk.file_prefix=LoadExternalFieldRZ_chk chk.file_min_digits=5 +inputFile = Examples/Tests/LoadExternalField/inputs_rz_grid_fields +runtime_params = warpx.abort_on_warning_threshold=medium chk.file_prefix=LoadExternalFieldRZGrid_chk chk.file_min_digits=5 +dim = 2 +addToCompileString = USE_RZ=TRUE +cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_OPENPMD=ON +restartTest = 1 +restartFileNum = 150 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = proton +analysisRoutine = Examples/Tests/LoadExternalField/analysis_rz.py + +[LoadExternalFieldRZParticles] +buildDir = . +inputFile = Examples/Tests/LoadExternalField/inputs_rz_particle_fields +runtime_params = warpx.abort_on_warning_threshold=medium chk.file_prefix=LoadExternalFieldRZParticles_chk chk.file_min_digits=5 dim = 2 addToCompileString = USE_RZ=TRUE cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_OPENPMD=ON @@ -3465,11 +3484,31 @@ compileTest = 0 doVis = 0 analysisRoutine = Examples/analysis_default_openpmd_regression.py -[Python_LoadExternalField3D] +[Python_LoadExternalGridField3D] buildDir = . -inputFile = Examples/Tests/LoadExternalField/PICMI_inputs_3d.py +inputFile = Examples/Tests/LoadExternalField/PICMI_inputs_3d_grid_fields.py runtime_params = -customRunCmd = python3 PICMI_inputs_3d.py +customRunCmd = python3 PICMI_inputs_3d_grid_fields.py +dim = 3 +addToCompileString = USE_PYTHON_MAIN=TRUE +cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_APP=OFF -DWarpX_PYTHON=ON -DWarpX_OPENPMD=ON +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = proton +analysisRoutine = Examples/Tests/LoadExternalField/analysis_3d.py + +[Python_LoadExternalParticleField3D] +buildDir = . +inputFile = Examples/Tests/LoadExternalField/PICMI_inputs_3d_particle_fields.py +runtime_params = +customRunCmd = python3 PICMI_inputs_3d_particle_fields.py dim = 3 addToCompileString = USE_PYTHON_MAIN=TRUE cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_APP=OFF -DWarpX_PYTHON=ON -DWarpX_OPENPMD=ON diff --git a/Source/Initialization/ExternalField.cpp b/Source/Initialization/ExternalField.cpp index 51a5d821148..d86c0a484bf 100644 --- a/Source/Initialization/ExternalField.cpp +++ b/Source/Initialization/ExternalField.cpp @@ -175,7 +175,6 @@ ExternalFieldParams::ExternalFieldParams(const amrex::ParmParse& pp_warpx) // // External fields from file // - if (E_ext_grid_type == ExternalFieldType::read_from_file || B_ext_grid_type == ExternalFieldType::read_from_file){ const std::string read_fields_from_path="./"; diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index cba5729d142..59e2b78774c 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1368,6 +1368,21 @@ void WarpX::CheckKnownIssues() void WarpX::LoadExternalFieldsFromFile (int const lev) { + // External fields from file are currently not compatible with the moving window + // In order to support the moving window, the MultiFab containing the external + // fields should be updated every time the window moves. + if ( (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) || + (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) || + (mypc->m_B_ext_particle_s == "read_from_file") || + (mypc->m_E_ext_particle_s == "read_from_file") ) { + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + WarpX::do_moving_window == 0, + "External fields from file are not compatible with the moving window." ); + } + + // External grid fields + if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) { #if defined(WARPX_DIM_RZ) WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, @@ -1392,6 +1407,41 @@ WarpX::LoadExternalFieldsFromFile (int const lev) ReadExternalFieldFromFile(m_p_ext_field_params->external_fields_path, Efield_fp_external[lev][0].get(), "E", "x"); ReadExternalFieldFromFile(m_p_ext_field_params->external_fields_path, Efield_fp_external[lev][1].get(), "E", "y"); ReadExternalFieldFromFile(m_p_ext_field_params->external_fields_path, Efield_fp_external[lev][2].get(), "E", "z"); +#endif + } + + // External particle fields + + if (mypc->m_B_ext_particle_s == "read_from_file") { + std::string external_fields_path; + const amrex::ParmParse pp_particles("particles"); + pp_particles.get("read_fields_from_path", external_fields_path ); +#if defined(WARPX_DIM_RZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, + "External field reading is not implemented for more than one RZ mode (see #3829)"); + ReadExternalFieldFromFile(external_fields_path, B_external_particle_field[lev][0].get(), "B", "r"); + ReadExternalFieldFromFile(external_fields_path, B_external_particle_field[lev][1].get(), "B", "t"); + ReadExternalFieldFromFile(external_fields_path, B_external_particle_field[lev][2].get(), "B", "z"); +#else + ReadExternalFieldFromFile(external_fields_path, B_external_particle_field[lev][0].get(), "B", "x"); + ReadExternalFieldFromFile(external_fields_path, B_external_particle_field[lev][1].get(), "B", "y"); + ReadExternalFieldFromFile(external_fields_path, B_external_particle_field[lev][2].get(), "B", "z"); +#endif + } + if (mypc->m_E_ext_particle_s == "read_from_file") { + std::string external_fields_path; + const amrex::ParmParse pp_particles("particles"); + pp_particles.get("read_fields_from_path", external_fields_path ); +#if defined(WARPX_DIM_RZ) + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(n_rz_azimuthal_modes == 1, + "External field reading is not implemented for more than one RZ mode (see #3829)"); + ReadExternalFieldFromFile(external_fields_path, E_external_particle_field[lev][0].get(), "E", "r"); + ReadExternalFieldFromFile(external_fields_path, E_external_particle_field[lev][1].get(), "E", "t"); + ReadExternalFieldFromFile(external_fields_path, E_external_particle_field[lev][2].get(), "E", "z"); +#else + ReadExternalFieldFromFile(external_fields_path, E_external_particle_field[lev][0].get(), "E", "x"); + ReadExternalFieldFromFile(external_fields_path, E_external_particle_field[lev][1].get(), "E", "y"); + ReadExternalFieldFromFile(external_fields_path, E_external_particle_field[lev][2].get(), "E", "z"); #endif } } diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index dcf5fbe5621..e7df489236e 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -18,6 +18,7 @@ #include "Utils/WarpXProfilerWrapper.H" #include "WarpXComm_K.H" #include "WarpXSumGuardCells.H" +#include "Particles/MultiParticleContainer.H" #include #include @@ -59,6 +60,21 @@ WarpX::UpdateAuxilaryData () } else { UpdateAuxilaryDataStagToNodal(); } + + // When loading particle fields from file: add the external fields: + for (int lev = 0; lev <= finest_level; ++lev) { + if (mypc->m_E_ext_particle_s == "read_from_file") { + amrex::MultiFab::Add(*Efield_aux[lev][0], *E_external_particle_field[lev][0], 0, 0, E_external_particle_field[lev][0]->nComp(), guard_cells.ng_FieldGather); + amrex::MultiFab::Add(*Efield_aux[lev][1], *E_external_particle_field[lev][1], 0, 0, E_external_particle_field[lev][1]->nComp(), guard_cells.ng_FieldGather); + amrex::MultiFab::Add(*Efield_aux[lev][2], *E_external_particle_field[lev][2], 0, 0, E_external_particle_field[lev][2]->nComp(), guard_cells.ng_FieldGather); + } + if (mypc->m_B_ext_particle_s == "read_from_file") { + amrex::MultiFab::Add(*Bfield_aux[lev][0], *B_external_particle_field[lev][0], 0, 0, B_external_particle_field[lev][0]->nComp(), guard_cells.ng_FieldGather); + amrex::MultiFab::Add(*Bfield_aux[lev][1], *B_external_particle_field[lev][1], 0, 0, B_external_particle_field[lev][0]->nComp(), guard_cells.ng_FieldGather); + amrex::MultiFab::Add(*Bfield_aux[lev][2], *B_external_particle_field[lev][2], 0, 0, B_external_particle_field[lev][0]->nComp(), guard_cells.ng_FieldGather); + } + } + } void @@ -293,6 +309,30 @@ WarpX::UpdateAuxilaryDataStagToNodal () void WarpX::UpdateAuxilaryDataSameType () { + // Update aux field, including guard cells, up to ng_FieldGather + const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; + + // Level 0: Copy from fine to aux + // Note: in some configurations, Efield_aux/Bfield_aux and Efield_fp/Bfield_fp are simply aliases to the + // same MultiFab object. MultiFab::Copy operation automatically detects this and does nothing in this case. + if (WarpX::fft_do_time_averaging) + { + MultiFab::Copy(*Efield_aux[0][0], *Efield_avg_fp[0][0], 0, 0, Efield_aux[0][0]->nComp(), ng_src); + MultiFab::Copy(*Efield_aux[0][1], *Efield_avg_fp[0][1], 0, 0, Efield_aux[0][1]->nComp(), ng_src); + MultiFab::Copy(*Efield_aux[0][2], *Efield_avg_fp[0][2], 0, 0, Efield_aux[0][2]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][0], *Bfield_avg_fp[0][0], 0, 0, Bfield_aux[0][0]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][1], *Bfield_avg_fp[0][1], 0, 0, Bfield_aux[0][1]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][2], *Bfield_avg_fp[0][2], 0, 0, Bfield_aux[0][2]->nComp(), ng_src); + } + else + { + MultiFab::Copy(*Efield_aux[0][0], *Efield_fp[0][0], 0, 0, Efield_aux[0][0]->nComp(), ng_src); + MultiFab::Copy(*Efield_aux[0][1], *Efield_fp[0][1], 0, 0, Efield_aux[0][1]->nComp(), ng_src); + MultiFab::Copy(*Efield_aux[0][2], *Efield_fp[0][2], 0, 0, Efield_aux[0][2]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][0], *Bfield_fp[0][0], 0, 0, Bfield_aux[0][0]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][1], *Bfield_fp[0][1], 0, 0, Bfield_aux[0][1]->nComp(), ng_src); + MultiFab::Copy(*Bfield_aux[0][2], *Bfield_fp[0][2], 0, 0, Bfield_aux[0][2]->nComp(), ng_src); + } for (int lev = 1; lev <= finest_level; ++lev) { const amrex::Periodicity& crse_period = Geom(lev-1).periodicity(); @@ -308,8 +348,6 @@ WarpX::UpdateAuxilaryDataSameType () dBy.setVal(0.0); dBz.setVal(0.0); - // Guard cells may not be up to date beyond ng_FieldGather - const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; // Copy Bfield_aux to the dB MultiFabs, using up to ng_src (=ng_FieldGather) guard // cells from Bfield_aux and filling up to ng (=nGrow) guard cells in the dB MultiFabs @@ -379,8 +417,6 @@ WarpX::UpdateAuxilaryDataSameType () dEy.setVal(0.0); dEz.setVal(0.0); - // Guard cells may not be up to date beyond ng_FieldGather - const amrex::IntVect& ng_src = guard_cells.ng_FieldGather; // Copy Efield_aux to the dE MultiFabs, using up to ng_src (=ng_FieldGather) guard // cells from Efield_aux and filling up to ng (=nGrow) guard cells in the dE MultiFabs ablastr::utils::communication::ParallelCopy(dEx, *Efield_aux[lev - 1][0], 0, 0, diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index 52586a5c88c..d33bd17ccdd 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -192,6 +192,12 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) { RemakeMultiFab(Efield_fp_external[lev][idim], true); } + if (mypc->m_B_ext_particle_s == "read_from_file") { + RemakeMultiFab(B_external_particle_field[lev][idim], true); + } + if (mypc->m_E_ext_particle_s == "read_from_file") { + RemakeMultiFab(E_external_particle_field[lev][idim], true); + } RemakeMultiFab(current_fp[lev][idim], false); RemakeMultiFab(current_store[lev][idim], false); if (current_deposition_algo == CurrentDepositionAlgo::Vay) { diff --git a/Source/Particles/Gather/GetExternalFields.cpp b/Source/Particles/Gather/GetExternalFields.cpp index 8e3c679aa3a..bb55f79f394 100644 --- a/Source/Particles/Gather/GetExternalFields.cpp +++ b/Source/Particles/Gather/GetExternalFields.cpp @@ -83,6 +83,18 @@ GetExternalEBField::GetExternalEBField (const WarpXParIter& a_pti, long a_offset m_repeated_plasma_lens_strengths_B = mypc.d_repeated_plasma_lens_strengths_B.data(); } + // When the external particle fields are read from file, + // the external fields are not added directly inside the gather kernel. + // (Hence of `None`, which ensures that the gather kernel is compiled without support + // for external fields.) Instead, the external fields are added to the MultiFab + // Efield_aux and Bfield_aux before the particles gather from these MultiFab. + if (mypc.m_E_ext_particle_s == "read_from_file") { + m_Etype = ExternalFieldInitType::None; + } + if (mypc.m_B_ext_particle_s == "read_from_file") { + m_Btype = ExternalFieldInitType::None; + } + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_Etype != Unknown, "Unknown E_ext_particle_init_style"); WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_Btype != Unknown, "Unknown B_ext_particle_init_style"); diff --git a/Source/WarpX.H b/Source/WarpX.H index 55d0f9ec3e5..0a64c550a27 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -1518,6 +1518,8 @@ private: // Same as Bfield_fp/Efield_fp for reading external field data amrex::Vector, 3 > > Efield_fp_external; amrex::Vector, 3 > > Bfield_fp_external; + amrex::Vector, 3 > > E_external_particle_field; + amrex::Vector, 3 > > B_external_particle_field; //! EB: Lengths of the mesh edges amrex::Vector, 3 > > m_edge_lengths; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 2e6532b0ea0..82ce6dda438 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -344,6 +344,8 @@ WarpX::WarpX () // Same as Bfield_fp/Efield_fp for reading external field data Bfield_fp_external.resize(1); Efield_fp_external.resize(1); + B_external_particle_field.resize(1); + E_external_particle_field.resize(1); m_edge_lengths.resize(nlevs_max); m_face_areas.resize(nlevs_max); @@ -2562,23 +2564,35 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(Efield_aux[lev][1], nba, dm, ncomps, ngEB, lev, "Efield_aux[y]", 0.0_rt); AllocInitMultiFab(Efield_aux[lev][2], nba, dm, ncomps, ngEB, lev, "Efield_aux[z]", 0.0_rt); } else if (lev == 0) { - if (!WarpX::fft_do_time_averaging) { - // In this case, the aux grid is simply an alias of the fp grid - AliasInitMultiFab(Efield_aux[lev][0], *Efield_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][1], *Efield_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][2], *Efield_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); - - AliasInitMultiFab(Bfield_aux[lev][0], *Bfield_fp[lev][0], 0, ncomps, lev, "Bfield_aux[x]", 0.0_rt); - AliasInitMultiFab(Bfield_aux[lev][1], *Bfield_fp[lev][1], 0, ncomps, lev, "Bfield_aux[y]", 0.0_rt); - AliasInitMultiFab(Bfield_aux[lev][2], *Bfield_fp[lev][2], 0, ncomps, lev, "Bfield_aux[z]", 0.0_rt); - } else { - AliasInitMultiFab(Efield_aux[lev][0], *Efield_avg_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][1], *Efield_avg_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); - AliasInitMultiFab(Efield_aux[lev][2], *Efield_avg_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); - + if (WarpX::fft_do_time_averaging) { AliasInitMultiFab(Bfield_aux[lev][0], *Bfield_avg_fp[lev][0], 0, ncomps, lev, "Bfield_aux[x]", 0.0_rt); AliasInitMultiFab(Bfield_aux[lev][1], *Bfield_avg_fp[lev][1], 0, ncomps, lev, "Bfield_aux[y]", 0.0_rt); AliasInitMultiFab(Bfield_aux[lev][2], *Bfield_avg_fp[lev][2], 0, ncomps, lev, "Bfield_aux[z]", 0.0_rt); + + AliasInitMultiFab(Efield_aux[lev][0], *Efield_avg_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][1], *Efield_avg_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][2], *Efield_avg_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); + } else { + if (mypc->m_B_ext_particle_s == "read_from_file") { + AllocInitMultiFab(Bfield_aux[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[x]"); + AllocInitMultiFab(Bfield_aux[lev][1], amrex::convert(ba, By_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[y]"); + AllocInitMultiFab(Bfield_aux[lev][2], amrex::convert(ba, Bz_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[z]"); + } else { + // In this case, the aux grid is simply an alias of the fp grid (most common case in WarpX) + AliasInitMultiFab(Bfield_aux[lev][0], *Bfield_fp[lev][0], 0, ncomps, lev, "Bfield_aux[x]", 0.0_rt); + AliasInitMultiFab(Bfield_aux[lev][1], *Bfield_fp[lev][1], 0, ncomps, lev, "Bfield_aux[y]", 0.0_rt); + AliasInitMultiFab(Bfield_aux[lev][2], *Bfield_fp[lev][2], 0, ncomps, lev, "Bfield_aux[z]", 0.0_rt); + } + if (mypc->m_E_ext_particle_s == "read_from_file") { + AllocInitMultiFab(Efield_aux[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, lev, "Efield_aux[x]"); + AllocInitMultiFab(Efield_aux[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, lev, "Efield_aux[y]"); + AllocInitMultiFab(Efield_aux[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, lev, "Efield_aux[z]"); + } else { + // In this case, the aux grid is simply an alias of the fp grid (most common case in WarpX) + AliasInitMultiFab(Efield_aux[lev][0], *Efield_fp[lev][0], 0, ncomps, lev, "Efield_aux[x]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][1], *Efield_fp[lev][1], 0, ncomps, lev, "Efield_aux[y]", 0.0_rt); + AliasInitMultiFab(Efield_aux[lev][2], *Efield_fp[lev][2], 0, ncomps, lev, "Efield_aux[z]", 0.0_rt); + } } } else { AllocInitMultiFab(Bfield_aux[lev][0], amrex::convert(ba, Bx_nodal_flag), dm, ncomps, ngEB, lev, "Bfield_aux[x]"); @@ -2590,6 +2604,44 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(Efield_aux[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, lev, "Efield_aux[z]"); } + // The external fields that are read from file + if (m_p_ext_field_params->B_ext_grid_type == ExternalFieldType::read_from_file) { + // These fields will be added directly to the grid, i.e. to fp, and need to match the index type + AllocInitMultiFab(Bfield_fp_external[lev][0], amrex::convert(ba, Bfield_fp[lev][0]->ixType()), + dm, ncomps, ngEB, lev, "Bfield_fp_external[x]", 0.0_rt); + AllocInitMultiFab(Bfield_fp_external[lev][1], amrex::convert(ba, Bfield_fp[lev][1]->ixType()), + dm, ncomps, ngEB, lev, "Bfield_fp_external[y]", 0.0_rt); + AllocInitMultiFab(Bfield_fp_external[lev][2], amrex::convert(ba, Bfield_fp[lev][2]->ixType()), + dm, ncomps, ngEB, lev, "Bfield_fp_external[z]", 0.0_rt); + } + if (mypc->m_B_ext_particle_s == "read_from_file") { + // These fields will be added to the fields that the particles see, and need to match the index type + AllocInitMultiFab(B_external_particle_field[lev][0], amrex::convert(ba, Bfield_aux[lev][0]->ixType()), + dm, ncomps, ngEB, lev, "B_external_particle_field[x]", 0.0_rt); + AllocInitMultiFab(B_external_particle_field[lev][1], amrex::convert(ba, Bfield_aux[lev][1]->ixType()), + dm, ncomps, ngEB, lev, "B_external_particle_field[y]", 0.0_rt); + AllocInitMultiFab(B_external_particle_field[lev][2], amrex::convert(ba, Bfield_aux[lev][2]->ixType()), + dm, ncomps, ngEB, lev, "B_external_particle_field[z]", 0.0_rt); + } + if (m_p_ext_field_params->E_ext_grid_type == ExternalFieldType::read_from_file) { + // These fields will be added directly to the grid, i.e. to fp, and need to match the index type + AllocInitMultiFab(Efield_fp_external[lev][0], amrex::convert(ba, Efield_fp[lev][0]->ixType()), + dm, ncomps, ngEB, lev, "Efield_fp_external[x]", 0.0_rt); + AllocInitMultiFab(Efield_fp_external[lev][1], amrex::convert(ba, Efield_fp[lev][1]->ixType()), + dm, ncomps, ngEB, lev, "Efield_fp_external[y]", 0.0_rt); + AllocInitMultiFab(Efield_fp_external[lev][2], amrex::convert(ba, Efield_fp[lev][2]->ixType()), + dm, ncomps, ngEB, lev, "Efield_fp_external[z]", 0.0_rt); + } + if (mypc->m_E_ext_particle_s == "read_from_file") { + // These fields will be added to the fields that the particles see, and need to match the index type + AllocInitMultiFab(E_external_particle_field[lev][0], amrex::convert(ba, Efield_aux[lev][0]->ixType()), + dm, ncomps, ngEB, lev, "E_external_particle_field[x]", 0.0_rt); + AllocInitMultiFab(E_external_particle_field[lev][1], amrex::convert(ba, Efield_aux[lev][1]->ixType()), + dm, ncomps, ngEB, lev, "E_external_particle_field[y]", 0.0_rt); + AllocInitMultiFab(E_external_particle_field[lev][2], amrex::convert(ba, Efield_aux[lev][2]->ixType()), + dm, ncomps, ngEB, lev, "E_external_particle_field[z]", 0.0_rt); + } + // // The coarse patch // diff --git a/requirements.txt b/requirements.txt index 60f14e1282d..bb0d8ebb704 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,7 +4,7 @@ periodictable~=1.5 # PICMI # note: don't forget to update the version in Docs/requirements.txt, too -picmistandard==0.28.0 +picmistandard==0.29.0 # for development against an unreleased PICMI version, use: #picmistandard @ git+https://github.com/picmi-standard/picmi.git#subdirectory=PICMI_Python