diff options
| -rw-r--r-- | Python/pywarpx/PGroup.py | 8 | ||||
| -rwxr-xr-x | Python/pywarpx/_libwarpx.py | 280 | ||||
| -rw-r--r-- | Python/pywarpx/picmi.py | 2 | ||||
| -rw-r--r-- | Source/Diagnostics/FieldIO.cpp | 16 | ||||
| -rw-r--r-- | Source/Diagnostics/WarpXIO.cpp | 3 | ||||
| -rw-r--r-- | Source/Evolve/WarpXEvolveEM.cpp | 19 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 2 | ||||
| -rw-r--r-- | Source/FortranInterface/WarpX_f.H | 5 | ||||
| -rw-r--r-- | Source/FortranInterface/WarpX_picsar.F90 | 437 | ||||
| -rw-r--r-- | Source/Parallelization/WarpXComm.cpp | 102 | ||||
| -rw-r--r-- | Source/Parallelization/WarpXRegrid.cpp | 66 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 2 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 4 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 1 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 44 | ||||
| -rw-r--r-- | Source/Python/WarpXWrappers.cpp | 46 | ||||
| -rw-r--r-- | Source/WarpX.H | 4 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 112 |
18 files changed, 705 insertions, 448 deletions
diff --git a/Python/pywarpx/PGroup.py b/Python/pywarpx/PGroup.py index 68b37740d..48e68ceb5 100644 --- a/Python/pywarpx/PGroup.py +++ b/Python/pywarpx/PGroup.py @@ -83,6 +83,10 @@ class PGroup(object): return _libwarpx.get_particle_y(self.ispecie)[self.igroup] yp = property(getyp) + def getrp(self): + return _libwarpx.get_particle_r(self.ispecie)[self.igroup] + rp = property(getrp) + def getzp(self): return _libwarpx.get_particle_z(self.ispecie)[self.igroup] zp = property(getzp) @@ -136,6 +140,10 @@ class PGroup(object): return _libwarpx.get_particle_Bz(self.ispecie)[self.igroup] bz = property(getbz) + def gettheta(self): + return _libwarpx.get_particle_theta(self.ispecie)[self.igroup] + theta = property(gettheta) + class PGroups(object): def __init__(self, ispecie=0): self.ispecie = ispecie diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py index 4c3283b97..5b7e14dbb 100755 --- a/Python/pywarpx/_libwarpx.py +++ b/Python/pywarpx/_libwarpx.py @@ -399,7 +399,10 @@ def get_particle_x(species_number): ''' structs = get_particle_structs(species_number) - return [struct['x'] for struct in structs] + if geometry_dim == '3d' or geometry_dim == '2d': + return [struct['x'] for struct in structs] + elif geometry_dim == 'rz': + return [struct['x']*np.cos(theta) for struct, theta in zip(structs, get_particle_theta(species_number))] def get_particle_y(species_number): @@ -410,7 +413,26 @@ def get_particle_y(species_number): ''' structs = get_particle_structs(species_number) - return [struct['y'] for struct in structs] + if geometry_dim == '3d' or geometry_dim == '2d': + return [struct['y'] for struct in structs] + elif geometry_dim == 'rz': + return [struct['x']*np.sin(theta) for struct, theta in zip(structs, get_particle_theta(species_number))] + + +def get_particle_r(species_number): + ''' + + Return a list of numpy arrays containing the particle 'r' + positions on each tile. + + ''' + structs = get_particle_structs(species_number) + if geometry_dim == 'rz': + return [struct['x'] for struct in structs] + elif geometry_dim == '3d': + return [np.sqrt(struct['x']**2 + struct['y']**2) for struct in structs] + elif geometry_dim == '2d': + raise Exception('get_particle_r: There is no r coordinate with 2D Cartesian') def get_particle_z(species_number): @@ -556,6 +578,53 @@ def get_particle_Bz(species_number): return get_particle_arrays(species_number, 9) +def get_particle_theta(species_number): + ''' + + Return a list of numpy arrays containing the particle + theta on each tile. + + ''' + + if geometry_dim == 'rz': + return get_particle_arrays(species_number, 10) + elif geometry_dim == '3d': + return [np.arctan2(struct['y'], struct['x']) for struct in structs] + elif geometry_dim == '2d': + raise Exception('get_particle_r: There is no theta coordinate with 2D Cartesian') + + +def _get_mesh_field_list(warpx_func, level, direction, include_ghosts): + """ + Generic routine to fetch the list of field data arrays. + """ + shapes = _LP_c_int() + size = ctypes.c_int(0) + ncomps = ctypes.c_int(0) + ngrow = ctypes.c_int(0) + data = warpx_func(level, direction, + ctypes.byref(size), ctypes.byref(ncomps), + ctypes.byref(ngrow), ctypes.byref(shapes)) + ng = ngrow.value + grid_data = [] + shapesize = dim + if ncomps.value > 1: + shapesize += 1 + for i in range(size.value): + shape = tuple([shapes[shapesize*i + d] for d in range(shapesize)]) + # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. + arr = np.ctypeslib.as_array(data[i], shape[::-1]).T + arr.setflags(write=1) + if include_ghosts: + grid_data.append(arr) + else: + grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) + + _libc.free(shapes) + _libc.free(data) + return grid_data + + def get_mesh_electric_field(level, direction, include_ghosts=True): ''' @@ -582,28 +651,7 @@ def get_mesh_electric_field(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getEfield(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getEfield, level, direction, include_ghosts) def get_mesh_electric_field_cp(level, direction, include_ghosts=True): @@ -631,28 +679,7 @@ def get_mesh_electric_field_cp(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getEfieldCP(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getEfieldCP, level, direction, include_ghosts) def get_mesh_electric_field_fp(level, direction, include_ghosts=True): @@ -680,28 +707,7 @@ def get_mesh_electric_field_fp(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getEfieldFP(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getEfieldFP, level, direction, include_ghosts) def get_mesh_magnetic_field(level, direction, include_ghosts=True): @@ -730,28 +736,7 @@ def get_mesh_magnetic_field(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getBfield(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getBfield, level, direction, include_ghosts) def get_mesh_magnetic_field_cp(level, direction, include_ghosts=True): @@ -779,28 +764,7 @@ def get_mesh_magnetic_field_cp(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getBfieldCP(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getBfieldCP, level, direction, include_ghosts) def get_mesh_magnetic_field_fp(level, direction, include_ghosts=True): @@ -828,28 +792,7 @@ def get_mesh_magnetic_field_fp(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getBfieldFP(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getBfieldFP, level, direction, include_ghosts) def get_mesh_current_density(level, direction, include_ghosts=True): @@ -876,28 +819,7 @@ def get_mesh_current_density(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getCurrentDensity(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getCurrentDensity, level, direction, include_ghosts) def get_mesh_current_density_cp(level, direction, include_ghosts=True): @@ -925,28 +847,7 @@ def get_mesh_current_density_cp(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getCurrentDensityCP(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityCP, level, direction, include_ghosts) def get_mesh_current_density_fp(level, direction, include_ghosts=True): @@ -974,28 +875,7 @@ def get_mesh_current_density_fp(level, direction, include_ghosts=True): ''' assert(level == 0) - - shapes = _LP_c_int() - size = ctypes.c_int(0) - ngrow = ctypes.c_int(0) - data = libwarpx.warpx_getCurrentDensityFP(level, direction, - ctypes.byref(size), ctypes.byref(ngrow), - ctypes.byref(shapes)) - ng = ngrow.value - grid_data = [] - for i in range(size.value): - shape = tuple([shapes[dim*i + d] for d in range(dim)]) - # --- The data is stored in Fortran order, hence shape is reversed and a transpose is taken. - arr = np.ctypeslib.as_array(data[i], shape[::-1]).T - arr.setflags(write=1) - if include_ghosts: - grid_data.append(arr) - else: - grid_data.append(arr[[slice(ng, -ng) for _ in range(dim)]]) - - _libc.free(shapes) - _libc.free(data) - return grid_data + return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityFP, level, direction, include_ghosts) def _get_mesh_array_lovects(level, direction, include_ghosts=True, getarrayfunc=None): diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index f242f8589..151de388d 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -256,13 +256,13 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): assert self.lower_bound[0] >= 0., Exception('Lower radial boundary must be >= 0.') assert self.bc_rmin != 'periodic' and self.bc_rmax != 'periodic', Exception('Radial boundaries can not be periodic') - assert self.n_azimuthal_modes is None or self.n_azimuthal_modes == 1, Exception('Only one azimuthal mode supported') # Geometry pywarpx.geometry.coord_sys = 1 # RZ pywarpx.geometry.is_periodic = '0 %d'%(self.bc_zmin=='periodic') # Is periodic? pywarpx.geometry.prob_lo = self.lower_bound # physical domain pywarpx.geometry.prob_hi = self.upper_bound + pywarpx.warpx.nmodes = self.n_azimuthal_modes if self.moving_window_velocity is not None and np.any(np.not_equal(self.moving_window_velocity, 0.)): pywarpx.warpx.do_moving_window = 1 diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index e3d44d1fc..15eac3449 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -200,8 +200,14 @@ PackPlotDataPtrs (Vector<const MultiFab*>& pmf, pmf[1] = data[1].get(); pmf[2] = data[2].get(); #elif (AMREX_SPACEDIM == 2) - pmf[0] = data[0].get(); - pmf[1] = data[2].get(); + if (data[0]->nComp() > 1) { + // Only grabs the first component + pmf[0] = new MultiFab(*data[0], amrex::make_alias, 0, 1); + pmf[1] = new MultiFab(*data[2], amrex::make_alias, 0, 1); + } else { + pmf[0] = data[0].get(); + pmf[1] = data[2].get(); + } #endif } @@ -543,8 +549,8 @@ WriteRawField( const MultiFab& F, const DistributionMapping& dm, VisMF::Write(F, prefix); } else { // Copy original MultiFab into one that does not have guard cells - MultiFab tmpF( F.boxArray(), dm, 1, 0); - MultiFab::Copy(tmpF, F, 0, 0, 1, 0); + MultiFab tmpF( F.boxArray(), dm, F.nComp(), 0); + MultiFab::Copy(tmpF, F, 0, 0, F.nComp(), 0); VisMF::Write(tmpF, prefix); } @@ -566,7 +572,7 @@ WriteZeroRawField( const MultiFab& F, const DistributionMapping& dm, std::string prefix = amrex::MultiFabFileFullPrefix(lev, filename, level_prefix, field_name); - MultiFab tmpF(F.boxArray(), dm, 1, ng); + MultiFab tmpF(F.boxArray(), dm, F.nComp(), ng); tmpF.setVal(0.); VisMF::Write(tmpF, prefix); } diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 6d646dd42..59ba04c5f 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -580,7 +580,8 @@ WarpX::WritePlotFile () const if (F_fp[lev]) WriteRawField( *F_fp[lev], dm, raw_pltname, level_prefix, "F_fp", lev, plot_raw_fields_guards); if (plot_rho) { // Use the component 1 of `rho_fp`, i.e. rho_new for time synchronization - MultiFab rho_new(*rho_fp[lev], amrex::make_alias, 1, 1); + // If nComp > 1, this is the upper half of the list of components. + MultiFab rho_new(*rho_fp[lev], amrex::make_alias, rho_fp[lev]->nComp()/2, rho_fp[lev]->nComp()/2); WriteRawField( rho_new, dm, raw_pltname, level_prefix, "rho_fp", lev, plot_raw_fields_guards); } } diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index ad7c7d840..dc6cfddeb 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -477,8 +477,23 @@ WarpX::ComputeDt () if (maxwell_fdtd_solver_id == 0) { // CFL time step Yee solver #ifdef WARPX_RZ - // Derived semi-analytically by R. Lehe - deltat = cfl * 1./( std::sqrt((1+0.2105)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); + // In the rz case, the Courant limit has been evaluated + // semi-analytically by R. Lehe, and resulted in the following + // coefficients. For an explanation, see (not officially published) + // www.normalesup.org/~lehe/Disp_relation_Circ.pdf + // NB : Here the coefficient for m=1 as compared to this document, + // as it was observed in practice that this coefficient was not + // high enough (The simulation became unstable). + Real multimode_coeffs[6] = { 0.2105, 1.0, 3.5234, 8.5104, 15.5059, 24.5037 }; + Real multimode_alpha; + if (nmodes < 7) { + // Use the table of the coefficients + multimode_alpha = multimode_coeffs[nmodes-1]; + } else { + // Use a realistic extrapolation + multimode_alpha = (nmodes - 1)*(nmodes - 1) - 0.4; + } + deltat = cfl * 1./( std::sqrt((1+multimode_alpha)/(dx[0]*dx[0]) + 1./(dx[1]*dx[1])) * PhysConst::c ); #else deltat = cfl * 1./( std::sqrt(AMREX_D_TERM( 1./(dx[0]*dx[0]), + 1./(dx[1]*dx[1]), diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index c53e13f8f..fc4fb902b 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -109,6 +109,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) tbx.loVect(), tbx.hiVect(), tby.loVect(), tby.hiVect(), tbz.loVect(), tbz.hiVect(), + &nmodes, BL_TO_FORTRAN_3D((*Ex)[mfi]), BL_TO_FORTRAN_3D((*Ey)[mfi]), BL_TO_FORTRAN_3D((*Ez)[mfi]), @@ -271,6 +272,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), tez.loVect(), tez.hiVect(), + &nmodes, BL_TO_FORTRAN_3D((*Ex)[mfi]), BL_TO_FORTRAN_3D((*Ey)[mfi]), BL_TO_FORTRAN_3D((*Ez)[mfi]), diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 52996a60a..ba36c44ab 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -109,6 +109,7 @@ extern "C" amrex::Real* jx, const long* jx_ng, const int* jx_ntot, amrex::Real* jy, const long* jy_ng, const int* jy_ntot, amrex::Real* jz, const long* jz_ng, const int* jz_ntot, + const long* nmodes, const long* np, const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp, const amrex::Real* uxp, const amrex::Real* uyp,const amrex::Real* uzp, @@ -124,6 +125,7 @@ extern "C" amrex::Real* jx, const long* jx_ng, const int* jx_ntot, amrex::Real* jy, const long* jy_ng, const int* jy_ntot, amrex::Real* jz, const long* jz_ng, const int* jz_ntot, + const long* nmodes, const amrex::Real* rmin, const amrex::Real* dr); @@ -143,6 +145,7 @@ extern "C" const amrex::Real* bxg, const int* bxg_lo, const int* bxg_hi, const amrex::Real* byg, const int* byg_lo, const int* byg_hi, const amrex::Real* bzg, const int* bzg_lo, const int* bzg_hi, + const long* nmodes, const int* ll4symtry, const int* l_lower_order_in_v, const int* l_nodal, const long* lvect, const long* field_gathe_algo); @@ -201,6 +204,7 @@ extern "C" const int* xlo, const int* xhi, const int* ylo, const int* yhi, const int* zlo, const int* zhi, + const long* nmodes, BL_FORT_FAB_ARG_3D(ex), BL_FORT_FAB_ARG_3D(ey), BL_FORT_FAB_ARG_3D(ez), @@ -221,6 +225,7 @@ extern "C" const int* xlo, const int* xhi, const int* ylo, const int* yhi, const int* zlo, const int* zhi, + const long* nmodes, const BL_FORT_FAB_ARG_3D(ex), const BL_FORT_FAB_ARG_3D(ey), const BL_FORT_FAB_ARG_3D(ez), diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index c17e8861b..808135c61 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -89,6 +89,7 @@ contains ex,ey,ez,bx,by,bz,ixyzmin,xmin,ymin,zmin,dx,dy,dz,nox,noy,noz, & exg,exg_lo,exg_hi,eyg,eyg_lo,eyg_hi,ezg,ezg_lo,ezg_hi, & bxg,bxg_lo,bxg_hi,byg,byg_lo,byg_hi,bzg,bzg_lo,bzg_hi, & + nmodes, & ll4symtry,l_lower_order_in_v, l_nodal,& lvect,field_gathe_algo) & bind(C, name="warpx_geteb_energy_conserving") @@ -100,12 +101,22 @@ contains integer, intent(in) :: ixyzmin(AMREX_SPACEDIM) real(amrex_real), intent(in) :: xmin,ymin,zmin,dx,dy,dz integer(c_long), intent(in) :: field_gathe_algo - integer(c_long), intent(in) :: np,nox,noy,noz + integer(c_long), intent(in) :: np,nmodes,nox,noy,noz integer(c_int), intent(in) :: ll4symtry,l_lower_order_in_v, l_nodal integer(c_long),intent(in) :: lvect real(amrex_real), intent(in), dimension(np) :: xp,yp,zp real(amrex_real), intent(out), dimension(np) :: ex,ey,ez,bx,by,bz +#ifdef WARPX_RZ + real(amrex_real),intent(in):: exg(exg_lo(1):exg_hi(1),exg_lo(2):exg_hi(2),2*nmodes) + real(amrex_real),intent(in):: eyg(eyg_lo(1):eyg_hi(1),eyg_lo(2):eyg_hi(2),2*nmodes) + real(amrex_real),intent(in):: ezg(ezg_lo(1):ezg_hi(1),ezg_lo(2):ezg_hi(2),2*nmodes) + real(amrex_real),intent(in):: bxg(bxg_lo(1):bxg_hi(1),bxg_lo(2):bxg_hi(2),2*nmodes) + real(amrex_real),intent(in):: byg(byg_lo(1):byg_hi(1),byg_lo(2):byg_hi(2),2*nmodes) + real(amrex_real),intent(in):: bzg(bzg_lo(1):bzg_hi(1),bzg_lo(2):bzg_hi(2),2*nmodes) +#else real(amrex_real),intent(in):: exg(*), eyg(*), ezg(*), bxg(*), byg(*), bzg(*) +#endif + logical(pxr_logical) :: pxr_ll4symtry, pxr_l_lower_order_in_v, pxr_l_nodal ! Compute the number of valid cells and guard cells @@ -114,6 +125,11 @@ contains exg_nguards(AMREX_SPACEDIM), eyg_nguards(AMREX_SPACEDIM), ezg_nguards(AMREX_SPACEDIM), & bxg_nguards(AMREX_SPACEDIM), byg_nguards(AMREX_SPACEDIM), bzg_nguards(AMREX_SPACEDIM) +#ifdef WARPX_RZ + complex(amrex_real), allocatable, dimension(:,:,:) :: erg_c, etg_c, ezg_c, brg_c, btg_c, bzg_c + integer :: alloc_status +#endif + pxr_ll4symtry = ll4symtry .eq. 1 pxr_l_lower_order_in_v = l_lower_order_in_v .eq. 1 pxr_l_nodal = l_nodal .eq. 1 @@ -131,6 +147,41 @@ contains byg_nvalid = byg_lo + byg_hi - 2_c_long*ixyzmin + 1_c_long bzg_nvalid = bzg_lo + bzg_hi - 2_c_long*ixyzmin + 1_c_long +#ifdef WARPX_RZ + if (nmodes > 1) then + + allocate(erg_c(exg_lo(1):exg_hi(1),exg_lo(2):exg_hi(2),nmodes), & + etg_c(eyg_lo(1):eyg_hi(1),eyg_lo(2):eyg_hi(2),nmodes), & + ezg_c(ezg_lo(1):ezg_hi(1),ezg_lo(2):ezg_hi(2),nmodes), & + brg_c(bxg_lo(1):bxg_hi(1),bxg_lo(2):bxg_hi(2),nmodes), & + btg_c(byg_lo(1):byg_hi(1),byg_lo(2):byg_hi(2),nmodes), & + bzg_c(bzg_lo(1):bzg_hi(1),bzg_lo(2):bzg_hi(2),nmodes), stat=alloc_status) + if (alloc_status /= 0) then + print*,"Error: warpx_geteb_energy_conserving: complex arrays could not be allocated" + stop + endif + + erg_c(:,:,:) = cmplx(exg(:,:,1::2), exg(:,:,2::2), amrex_real) + etg_c(:,:,:) = cmplx(eyg(:,:,1::2), eyg(:,:,2::2), amrex_real) + ezg_c(:,:,:) = cmplx(ezg(:,:,1::2), ezg(:,:,2::2), amrex_real) + brg_c(:,:,:) = cmplx(bxg(:,:,1::2), bxg(:,:,2::2), amrex_real) + btg_c(:,:,:) = cmplx(byg(:,:,1::2), byg(:,:,2::2), amrex_real) + bzg_c(:,:,:) = cmplx(bzg(:,:,1::2), bzg(:,:,2::2), amrex_real) + + call geteb2dcirc_energy_conserving_generic(np, xp, yp, zp, ex, ey, ez, bx, by, bz, & + xmin, zmin, dx, dz, nmodes, nox, noz, & + pxr_l_lower_order_in_v, pxr_l_nodal, & + erg_c, exg_nguards, exg_nvalid, & + etg_c, eyg_nguards, eyg_nvalid, & + ezg_c, ezg_nguards, ezg_nvalid, & + brg_c, bxg_nguards, bxg_nvalid, & + btg_c, byg_nguards, byg_nvalid, & + bzg_c, bzg_nguards, bzg_nvalid) + + deallocate(erg_c, etg_c, ezg_c, brg_c, btg_c, bzg_c) + + else +#endif CALL WRPX_PXR_GETEB_ENERGY_CONSERVING(np,xp,yp,zp, & ex,ey,ez,bx,by,bz,xmin,ymin,zmin,dx,dy,dz,nox,noy,noz, & exg,exg_nguards,exg_nvalid,& @@ -141,6 +192,9 @@ contains bzg,bzg_nguards,bzg_nvalid,& pxr_ll4symtry, pxr_l_lower_order_in_v, pxr_l_nodal, & lvect, field_gathe_algo ) +#ifdef WARPX_RZ + endif +#endif end subroutine warpx_geteb_energy_conserving @@ -308,17 +362,24 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> @param[in] charge_depo_algo algorithm choice for the charge deposition !> subroutine warpx_current_deposition( & - jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot, & + jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot,nmodes, & np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q,xmin,ymin,zmin,dt,dx,dy,dz,nox,noy,noz,& lvect,current_depo_algo) & bind(C, name="warpx_current_deposition") integer, intent(in) :: jx_ntot(AMREX_SPACEDIM), jy_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM) integer(c_long), intent(in) :: jx_ng, jy_ng, jz_ng + integer(c_long), intent(IN) :: nmodes integer(c_long), intent(IN) :: np integer(c_long), intent(IN) :: nox,noy,noz +#ifdef WARPX_RZ + real(amrex_real), intent(IN OUT):: jx(jx_ntot(1),jx_ntot(2),2*nmodes) + real(amrex_real), intent(IN OUT):: jy(jy_ntot(1),jy_ntot(2),2*nmodes) + real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2*nmodes) +#else real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) +#endif real(amrex_real), intent(IN) :: q real(amrex_real), intent(IN) :: dx,dy,dz real(amrex_real), intent(IN) :: dt @@ -329,6 +390,13 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n integer(c_long), intent(IN) :: lvect integer(c_long), intent(IN) :: current_depo_algo +#ifdef WARPX_RZ + logical(pxr_logical) :: l_particles_weight = .true. + integer(c_long) :: type_rz_depose = 1 + complex(amrex_real), allocatable, dimension(:,:,:) :: jr_c, jt_c, jz_c + integer :: alloc_status +#endif + ! Compute the number of valid cells and guard cells integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) @@ -350,13 +418,53 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n nox,noy,noz,current_depo_algo) ! Dimension 2 #elif (AMREX_SPACEDIM==2) - CALL WRPX_PXR_CURRENT_DEPOSITION( & +#ifdef WARPX_RZ + if (nmodes > 1) then + + allocate(jr_c(jx_ntot(1),jx_ntot(2),nmodes), & + jt_c(jy_ntot(1),jy_ntot(2),nmodes), & + jz_c(jz_ntot(1),jz_ntot(2),nmodes), stat=alloc_status) + if (alloc_status /= 0) then + print*,"Error: warpx_current_deposition: complex arrays could not be allocated" + stop + endif + + jr_c = 0._amrex_real + jt_c = 0._amrex_real + jz_c = 0._amrex_real + + CALL pxr_depose_jrjtjz_esirkepov_n_2d_circ( & + jr_c,jx_nguards,jx_nvalid, & + jt_c,jy_nguards,jy_nvalid, & + jz_c,jz_nguards,jz_nvalid, & + nmodes, & + np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, & + xmin,zmin,dt,dx,dz, & + nox,noz,l_particles_weight,type_rz_depose) + + jx(:,:,1::2) = jx(:,:,1::2) + real(jr_c) + jx(:,:,2::2) = jx(:,:,2::2) + aimag(jr_c) + jy(:,:,1::2) = jy(:,:,1::2) + real(jt_c) + jy(:,:,2::2) = jy(:,:,2::2) + aimag(jt_c) + jz(:,:,1::2) = jz(:,:,1::2) + real(jz_c) + jz(:,:,2::2) = jz(:,:,2::2) + aimag(jz_c) + + deallocate(jr_c) + deallocate(jt_c) + deallocate(jz_c) + + else +#endif + CALL WRPX_PXR_CURRENT_DEPOSITION( & jx,jx_nguards,jx_nvalid, & jy,jy_nguards,jy_nvalid, & jz,jz_nguards,jz_nvalid, & np,xp,yp,zp,uxp,uyp,uzp,gaminv,w,q, & xmin,zmin,dt,dx,dz,nox,noz,lvect, & current_depo_algo) +#ifdef WARPX_RZ + endif +#endif #endif end subroutine warpx_current_deposition @@ -367,45 +475,87 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> Applies the inverse volume scaling for RZ current deposition !> !> @details - !> The scaling is done for single mode only + !> The scaling is done for all modes ! - !> @param[inout] jx,jy,jz current arrays - !> @param[in] jx_ntot,jy_ntot,jz_ntot vectors with total number of + !> @param[inout] jr,jt,jz current arrays + !> @param[in] jr_ntot,jt_ntot,jz_ntot vectors with total number of !> cells (including guard cells) along each axis for each current - !> @param[in] jx_ng,jy_ng,jz_ng vectors with number of guard cells along each + !> @param[in] jr_ng,jt_ng,jz_ng vectors with number of guard cells along each !> axis for each current !> @param[in] rmin tile grid minimum radius !> @param[in] dr radial space discretization steps !> subroutine warpx_current_deposition_rz_volume_scaling( & - jx,jx_ng,jx_ntot,jy,jy_ng,jy_ntot,jz,jz_ng,jz_ntot, & - rmin,dr) & + jr,jr_ng,jr_ntot,jt,jt_ng,jt_ntot,jz,jz_ng,jz_ntot, & + nmodes,rmin,dr) & bind(C, name="warpx_current_deposition_rz_volume_scaling") - integer, intent(in) :: jx_ntot(AMREX_SPACEDIM), jy_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM) - integer(c_long), intent(in) :: jx_ng, jy_ng, jz_ng - real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) + integer, intent(in) :: jr_ntot(AMREX_SPACEDIM), jt_ntot(AMREX_SPACEDIM), jz_ntot(AMREX_SPACEDIM) + integer(c_long), intent(in) :: jr_ng, jt_ng, jz_ng + integer(c_long), intent(in) :: nmodes + real(amrex_real), intent(IN OUT):: jr(jr_ntot(1),jr_ntot(2),2*nmodes) + real(amrex_real), intent(IN OUT):: jt(jt_ntot(1),jt_ntot(2),2*nmodes) + real(amrex_real), intent(IN OUT):: jz(jz_ntot(1),jz_ntot(2),2*nmodes) real(amrex_real), intent(IN) :: rmin, dr #ifdef WARPX_RZ + + complex(amrex_real), allocatable, dimension(:,:,:) :: jr_c, jt_c, jz_c + integer :: alloc_status + integer(c_long) :: type_rz_depose = 1 -#endif + ! Compute the number of valid cells and guard cells - integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & - jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) - jx_nvalid = jx_ntot - 2*jx_ng - jy_nvalid = jy_ntot - 2*jy_ng + integer(c_long) :: jr_nvalid(AMREX_SPACEDIM), jt_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & + jr_nguards(AMREX_SPACEDIM), jt_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) + jr_nvalid = jr_ntot - 2*jr_ng + jt_nvalid = jt_ntot - 2*jt_ng jz_nvalid = jz_ntot - 2*jz_ng - jx_nguards = jx_ng - jy_nguards = jy_ng + jr_nguards = jr_ng + jt_nguards = jt_ng jz_nguards = jz_ng -#ifdef WARPX_RZ - CALL WRPX_PXR_RZ_VOLUME_SCALING_J( & - jx,jx_nguards,jx_nvalid, & - jy,jy_nguards,jy_nvalid, & + if (nmodes > 1) then + + allocate(jr_c(jr_ntot(1),jr_ntot(2),nmodes), & + jt_c(jt_ntot(1),jt_ntot(2),nmodes), & + jz_c(jz_ntot(1),jz_ntot(2),nmodes), stat=alloc_status) + if (alloc_status /= 0) then + print*,"Error: warpx_current_deposition_rz_volume_scaling: complex arrays could not be allocated" + stop + endif + + jr_c = cmplx(jr(:,:,1::2), jr(:,:,2::2), amrex_real) + jt_c = cmplx(jt(:,:,1::2), jt(:,:,2::2), amrex_real) + jz_c = cmplx(jz(:,:,1::2), jz(:,:,2::2), amrex_real) + + CALL apply_2dcirc_volume_scaling_j( & + jr_c, jr_nguards, jr_nvalid, & + jt_c, jt_nguards, jt_nvalid, & + jz_c, jz_nguards, jz_nvalid, & + nmodes, & + rmin,dr, & + type_rz_depose) + + jr(:,:,1::2) = real(jr_c) + jr(:,:,2::2) = aimag(jr_c) + jt(:,:,1::2) = real(jt_c) + jt(:,:,2::2) = aimag(jt_c) + jz(:,:,1::2) = real(jz_c) + jz(:,:,2::2) = aimag(jz_c) + + deallocate(jr_c) + deallocate(jt_c) + deallocate(jz_c) + + else + CALL apply_rz_volume_scaling_j( & + jr,jr_nguards,jr_nvalid, & + jt,jt_nguards,jt_nvalid, & jz,jz_nguards,jz_nvalid, & rmin,dr,type_rz_depose) + endif + #endif end subroutine warpx_current_deposition_rz_volume_scaling @@ -541,6 +691,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> @param[in] dtsdx, dtsdy, dtsdz factors c**2 * dt/(dx, dy, dz) subroutine warpx_push_evec( & xlo, xhi, ylo, yhi, zlo, zhi, & + nmodes, & ex, exlo, exhi, & ey, eylo, eyhi, & ez, ezlo, ezhi, & @@ -560,31 +711,115 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jxlo(BL_SPACEDIM), jxhi(BL_SPACEDIM), jylo(BL_SPACEDIM), jyhi(BL_SPACEDIM), & jzlo(BL_SPACEDIM), jzhi(BL_SPACEDIM) - real(amrex_real), intent(IN OUT):: ex(*), ey(*), ez(*) + integer(c_long), intent(in) :: nmodes - real(amrex_real), intent(IN):: bx(*), by(*), bz(*), jx(*), jy(*), jz(*) +#ifdef WARPX_RZ + real(amrex_real), intent(IN OUT):: ex(exlo(1):exhi(1),exlo(2):exhi(2),2*nmodes) + real(amrex_real), intent(IN OUT):: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),2*nmodes) + real(amrex_real), intent(IN OUT):: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),2*nmodes) + real(amrex_real), intent(IN):: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),2*nmodes) + real(amrex_real), intent(IN):: by(bylo(1):byhi(1),bylo(2):byhi(2),2*nmodes) + real(amrex_real), intent(IN):: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2*nmodes) + real(amrex_real), intent(IN):: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2),2*nmodes) + real(amrex_real), intent(IN):: jy(jylo(1):jyhi(1),jylo(2):jyhi(2),2*nmodes) + real(amrex_real), intent(IN):: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2),2*nmodes) +#else + real(amrex_real), intent(IN OUT):: ex(*) + real(amrex_real), intent(IN OUT):: ey(*) + real(amrex_real), intent(IN OUT):: ez(*) + real(amrex_real), intent(IN):: bx(*) + real(amrex_real), intent(IN):: by(*) + real(amrex_real), intent(IN):: bz(*) + real(amrex_real), intent(IN):: jx(*) + real(amrex_real), intent(IN):: jy(*) + real(amrex_real), intent(IN):: jz(*) +#endif real(amrex_real), intent(IN) :: mudt, dtsdx, dtsdy, dtsdz real(amrex_real), intent(IN) :: xmin, dx - CALL WRPX_PXR_PUSH_EVEC(& - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi,& - ey, eylo, eyhi,& - ez, ezlo, ezhi,& - bx, bxlo, bxhi,& - by, bylo, byhi,& - bz, bzlo, bzhi,& - jx, jxlo, jxhi,& - jy, jylo, jyhi,& - jz, jzlo, jzhi,& - mudt, dtsdx, dtsdy, dtsdz & #ifdef WARPX_RZ - ,xmin,dx & + complex(amrex_real), allocatable, dimension(:,:,:) :: er_c, et_c, ez_c + complex(amrex_real), allocatable, dimension(:,:,:) :: br_c, bt_c, bz_c + complex(amrex_real), allocatable, dimension(:,:,:) :: jr_c, jt_c, jz_c + integer :: alloc_status + + if (nmodes == 1) then +#endif + CALL WRPX_PXR_PUSH_EVEC(& + xlo, xhi, ylo, yhi, zlo, zhi, & + ex, exlo, exhi,& + ey, eylo, eyhi,& + ez, ezlo, ezhi,& + bx, bxlo, bxhi,& + by, bylo, byhi,& + bz, bzlo, bzhi,& + jx, jxlo, jxhi,& + jy, jylo, jyhi,& + jz, jzlo, jzhi,& + mudt, dtsdx, dtsdy, dtsdz & +#ifdef WARPX_RZ + ,xmin, dx & #endif - ) + ) +#ifdef WARPX_RZ + else + + allocate(er_c(exlo(1):exhi(1),exlo(2):exhi(2),nmodes), & + et_c(eylo(1):eyhi(1),eylo(2):eyhi(2),nmodes), & + ez_c(ezlo(1):ezhi(1),ezlo(2):ezhi(2),nmodes), & + br_c(bxlo(1):bxhi(1),bxlo(2):bxhi(2),nmodes), & + bt_c(bylo(1):byhi(1),bylo(2):byhi(2),nmodes), & + bz_c(bzlo(1):bzhi(1),bzlo(2):bzhi(2),nmodes), & + jr_c(jxlo(1):jxhi(1),jxlo(2):jxhi(2),nmodes), & + jt_c(jylo(1):jyhi(1),jylo(2):jyhi(2),nmodes), & + jz_c(jzlo(1):jzhi(1),jzlo(2):jzhi(2),nmodes), stat=alloc_status) + if (alloc_status /= 0) then + print*,"Error: warpx_push_evec: complex arrays could not be allocated" + stop + endif + + er_c = cmplx(ex(:,:,1::2), ex(:,:,2::2), amrex_real) + et_c = cmplx(ey(:,:,1::2), ey(:,:,2::2), amrex_real) + ez_c = cmplx(ez(:,:,1::2), ez(:,:,2::2), amrex_real) + br_c = cmplx(bx(:,:,1::2), bx(:,:,2::2), amrex_real) + bt_c = cmplx(by(:,:,1::2), by(:,:,2::2), amrex_real) + bz_c = cmplx(bz(:,:,1::2), bz(:,:,2::2), amrex_real) + jr_c = cmplx(jx(:,:,1::2), jx(:,:,2::2), amrex_real) + jt_c = cmplx(jy(:,:,1::2), jy(:,:,2::2), amrex_real) + jz_c = cmplx(jz(:,:,1::2), jz(:,:,2::2), amrex_real) + + CALL pxrpush_emrz_evec_multimode(& + xlo, xhi, ylo, yhi, zlo, zhi, & + nmodes, & + er_c, exlo, exhi,& + et_c, eylo, eyhi,& + ez_c, ezlo, ezhi,& + br_c, bxlo, bxhi,& + bt_c, bylo, byhi,& + bz_c, bzlo, bzhi,& + jr_c, jxlo, jxhi,& + jt_c, jylo, jyhi,& + jz_c, jzlo, jzhi,& + mudt, dtsdx, dtsdy, dtsdz, xmin, dx & + ) + + ! Only E needs to be copied back + ex(:,:,1::2) = real(er_c) + ex(:,:,2::2) = aimag(er_c) + ey(:,:,1::2) = real(et_c) + ey(:,:,2::2) = aimag(et_c) + ez(:,:,1::2) = real(ez_c) + ez(:,:,2::2) = aimag(ez_c) + + deallocate(er_c, et_c, ez_c) + deallocate(br_c, bt_c, bz_c) + deallocate(jr_c, jt_c, jz_c) + + endif +#endif end subroutine warpx_push_evec ! _________________________________________________________________ @@ -603,6 +838,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n !> @param[in] dtsdx, dtsdy, dtsdz factors 0.5 * dt/(dx, dy, dz) subroutine warpx_push_bvec( & xlo, xhi, ylo, yhi, zlo, zhi, & + nmodes, & ex, exlo, exhi, & ey, eylo, eyhi, & ez, ezlo, ezhi, & @@ -619,41 +855,110 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n bylo(BL_SPACEDIM), byhi(BL_SPACEDIM), bzlo(BL_SPACEDIM), bzhi(BL_SPACEDIM), & maxwell_fdtd_solver_id - real(amrex_real), intent(IN OUT):: ex(*), ey(*), ez(*) + integer(c_long), intent(in) :: nmodes - real(amrex_real), intent(IN):: bx(*), by(*), bz(*) +#ifdef WARPX_RZ + real(amrex_real), intent(IN):: ex(exlo(1):exhi(1),exlo(2):exhi(2),2*nmodes) + real(amrex_real), intent(IN):: ey(eylo(1):eyhi(1),eylo(2):eyhi(2),2*nmodes) + real(amrex_real), intent(IN):: ez(ezlo(1):ezhi(1),ezlo(2):ezhi(2),2*nmodes) + real(amrex_real), intent(IN OUT):: bx(bxlo(1):bxhi(1),bxlo(2):bxhi(2),2*nmodes) + real(amrex_real), intent(IN OUT):: by(bylo(1):byhi(1),bylo(2):byhi(2),2*nmodes) + real(amrex_real), intent(IN OUT):: bz(bzlo(1):bzhi(1),bzlo(2):bzhi(2),2*nmodes) +#else + real(amrex_real), intent(IN):: ex(*) + real(amrex_real), intent(IN):: ey(*) + real(amrex_real), intent(IN):: ez(*) + real(amrex_real), intent(IN OUT):: bx(*) + real(amrex_real), intent(IN OUT):: by(*) + real(amrex_real), intent(IN OUT):: bz(*) +#endif real(amrex_real), intent(IN) :: dtsdx, dtsdy, dtsdz real(amrex_real), intent(IN) :: xmin, dx - IF (maxwell_fdtd_solver_id .eq. 0) THEN - ! Yee FDTD solver - CALL WRPX_PXR_PUSH_BVEC( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - bx, bxlo, bxhi, & - by, bylo, byhi, & - bz, bzlo, bzhi, & - dtsdx,dtsdy,dtsdz & #ifdef WARPX_RZ - ,xmin,dx & + complex(amrex_real), allocatable, dimension(:,:,:) :: er_c, et_c, ez_c + complex(amrex_real), allocatable, dimension(:,:,:) :: br_c, bt_c, bz_c + integer :: alloc_status #endif - ) - ELSE IF (maxwell_fdtd_solver_id .eq. 1) THEN - ! Cole-Karkkainen FDTD solver - CALL WRPX_PXR_PUSH_BVEC_CKC( & - xlo, xhi, ylo, yhi, zlo, zhi, & - ex, exlo, exhi, & - ey, eylo, eyhi, & - ez, ezlo, ezhi, & - bx, bxlo, bxhi, & - by, bylo, byhi, & - bz, bzlo, bzhi, & - dtsdx,dtsdy,dtsdz) - ENDIF + + if (nmodes == 1) then + + IF (maxwell_fdtd_solver_id .eq. 0) THEN + ! Yee FDTD solver + CALL WRPX_PXR_PUSH_BVEC( & + xlo, xhi, ylo, yhi, zlo, zhi, & + ex, exlo, exhi, & + ey, eylo, eyhi, & + ez, ezlo, ezhi, & + bx, bxlo, bxhi, & + by, bylo, byhi, & + bz, bzlo, bzhi, & + dtsdx,dtsdy,dtsdz & +#ifdef WARPX_RZ + ,xmin,dx & +#endif + ) + ELSE IF (maxwell_fdtd_solver_id .eq. 1) THEN + ! Cole-Karkkainen FDTD solver + CALL WRPX_PXR_PUSH_BVEC_CKC( & + xlo, xhi, ylo, yhi, zlo, zhi, & + ex, exlo, exhi, & + ey, eylo, eyhi, & + ez, ezlo, ezhi, & + bx, bxlo, bxhi, & + by, bylo, byhi, & + bz, bzlo, bzhi, & + dtsdx,dtsdy,dtsdz) + ENDIF + +#ifdef WARPX_RZ + else + + allocate(er_c(exlo(1):exhi(1),exlo(2):exhi(2),nmodes), & + et_c(eylo(1):eyhi(1),eylo(2):eyhi(2),nmodes), & + ez_c(ezlo(1):ezhi(1),ezlo(2):ezhi(2),nmodes), & + br_c(bxlo(1):bxhi(1),bxlo(2):bxhi(2),nmodes), & + bt_c(bylo(1):byhi(1),bylo(2):byhi(2),nmodes), & + bz_c(bzlo(1):bzhi(1),bzlo(2):bzhi(2),nmodes), stat=alloc_status) + if (alloc_status /= 0) then + print*,"Error: warpx_push_bvec: complex arrays could not be allocated" + stop + endif + + er_c = cmplx(ex(:,:,1::2), ex(:,:,2::2), amrex_real) + et_c = cmplx(ey(:,:,1::2), ey(:,:,2::2), amrex_real) + ez_c = cmplx(ez(:,:,1::2), ez(:,:,2::2), amrex_real) + br_c = cmplx(bx(:,:,1::2), bx(:,:,2::2), amrex_real) + bt_c = cmplx(by(:,:,1::2), by(:,:,2::2), amrex_real) + bz_c = cmplx(bz(:,:,1::2), bz(:,:,2::2), amrex_real) + + CALL pxrpush_emrz_bvec_multimode(& + xlo, xhi, ylo, yhi, zlo, zhi, & + nmodes, & + er_c, exlo, exhi,& + et_c, eylo, eyhi,& + ez_c, ezlo, ezhi,& + br_c, bxlo, bxhi,& + bt_c, bylo, byhi,& + bz_c, bzlo, bzhi,& + dtsdx, dtsdy, dtsdz, xmin, dx & + ) + + ! Only B needs to be copied back + bx(:,:,1::2) = real(br_c) + bx(:,:,2::2) = aimag(br_c) + by(:,:,1::2) = real(bt_c) + by(:,:,2::2) = aimag(bt_c) + bz(:,:,1::2) = real(bz_c) + bz(:,:,2::2) = aimag(bz_c) + + deallocate(er_c, et_c, ez_c) + deallocate(br_c, bt_c, bz_c) + +#endif + endif end subroutine warpx_push_bvec ! _________________________________________________________________ diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 5c9fa144f..67a6557b7 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -59,24 +59,24 @@ WarpX::UpdateAuxilaryData () // B field { - MultiFab dBx(Bfield_cp[lev][0]->boxArray(), dm, 1, ng); - MultiFab dBy(Bfield_cp[lev][1]->boxArray(), dm, 1, ng); - MultiFab dBz(Bfield_cp[lev][2]->boxArray(), dm, 1, ng); + MultiFab dBx(Bfield_cp[lev][0]->boxArray(), dm, Bfield_cp[lev][0]->nComp(), ng); + MultiFab dBy(Bfield_cp[lev][1]->boxArray(), dm, Bfield_cp[lev][1]->nComp(), ng); + MultiFab dBz(Bfield_cp[lev][2]->boxArray(), dm, Bfield_cp[lev][2]->nComp(), ng); dBx.setVal(0.0); dBy.setVal(0.0); dBz.setVal(0.0); - dBx.ParallelCopy(*Bfield_aux[lev-1][0], 0, 0, 1, ng, ng, crse_period); - dBy.ParallelCopy(*Bfield_aux[lev-1][1], 0, 0, 1, ng, ng, crse_period); - dBz.ParallelCopy(*Bfield_aux[lev-1][2], 0, 0, 1, ng, ng, crse_period); + dBx.ParallelCopy(*Bfield_aux[lev-1][0], 0, 0, Bfield_aux[lev-1][0]->nComp(), ng, ng, crse_period); + dBy.ParallelCopy(*Bfield_aux[lev-1][1], 0, 0, Bfield_aux[lev-1][1]->nComp(), ng, ng, crse_period); + dBz.ParallelCopy(*Bfield_aux[lev-1][2], 0, 0, Bfield_aux[lev-1][2]->nComp(), ng, ng, crse_period); if (Bfield_cax[lev][0]) { - MultiFab::Copy(*Bfield_cax[lev][0], dBx, 0, 0, 1, ng); - MultiFab::Copy(*Bfield_cax[lev][1], dBy, 0, 0, 1, ng); - MultiFab::Copy(*Bfield_cax[lev][2], dBz, 0, 0, 1, ng); + MultiFab::Copy(*Bfield_cax[lev][0], dBx, 0, 0, Bfield_cax[lev][0]->nComp(), ng); + MultiFab::Copy(*Bfield_cax[lev][1], dBy, 0, 0, Bfield_cax[lev][1]->nComp(), ng); + MultiFab::Copy(*Bfield_cax[lev][2], dBz, 0, 0, Bfield_cax[lev][2]->nComp(), ng); } - MultiFab::Subtract(dBx, *Bfield_cp[lev][0], 0, 0, 1, ng); - MultiFab::Subtract(dBy, *Bfield_cp[lev][1], 0, 0, 1, ng); - MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, 1, ng); + MultiFab::Subtract(dBx, *Bfield_cp[lev][0], 0, 0, Bfield_cp[lev][0]->nComp(), ng); + MultiFab::Subtract(dBy, *Bfield_cp[lev][1], 0, 0, Bfield_cp[lev][1]->nComp(), ng); + MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, Bfield_cp[lev][2]->nComp(), ng); const Real* dx = Geom(lev-1).CellSize(); const int refinement_ratio = refRatio(lev-1)[0]; @@ -134,24 +134,24 @@ WarpX::UpdateAuxilaryData () // E field { - MultiFab dEx(Efield_cp[lev][0]->boxArray(), dm, 1, ng); - MultiFab dEy(Efield_cp[lev][1]->boxArray(), dm, 1, ng); - MultiFab dEz(Efield_cp[lev][2]->boxArray(), dm, 1, ng); + MultiFab dEx(Efield_cp[lev][0]->boxArray(), dm, Efield_cp[lev][0]->nComp(), ng); + MultiFab dEy(Efield_cp[lev][1]->boxArray(), dm, Efield_cp[lev][1]->nComp(), ng); + MultiFab dEz(Efield_cp[lev][2]->boxArray(), dm, Efield_cp[lev][2]->nComp(), ng); dEx.setVal(0.0); dEy.setVal(0.0); dEz.setVal(0.0); - dEx.ParallelCopy(*Efield_aux[lev-1][0], 0, 0, 1, ng, ng, crse_period); - dEy.ParallelCopy(*Efield_aux[lev-1][1], 0, 0, 1, ng, ng, crse_period); - dEz.ParallelCopy(*Efield_aux[lev-1][2], 0, 0, 1, ng, ng, crse_period); + dEx.ParallelCopy(*Efield_aux[lev-1][0], 0, 0, Efield_aux[lev-1][0]->nComp(), ng, ng, crse_period); + dEy.ParallelCopy(*Efield_aux[lev-1][1], 0, 0, Efield_aux[lev-1][1]->nComp(), ng, ng, crse_period); + dEz.ParallelCopy(*Efield_aux[lev-1][2], 0, 0, Efield_aux[lev-1][2]->nComp(), ng, ng, crse_period); if (Efield_cax[lev][0]) { - MultiFab::Copy(*Efield_cax[lev][0], dEx, 0, 0, 1, ng); - MultiFab::Copy(*Efield_cax[lev][1], dEy, 0, 0, 1, ng); - MultiFab::Copy(*Efield_cax[lev][2], dEz, 0, 0, 1, ng); + MultiFab::Copy(*Efield_cax[lev][0], dEx, 0, 0, Efield_cax[lev][0]->nComp(), ng); + MultiFab::Copy(*Efield_cax[lev][1], dEy, 0, 0, Efield_cax[lev][1]->nComp(), ng); + MultiFab::Copy(*Efield_cax[lev][2], dEz, 0, 0, Efield_cax[lev][2]->nComp(), ng); } - MultiFab::Subtract(dEx, *Efield_cp[lev][0], 0, 0, 1, ng); - MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, 1, ng); - MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, 1, ng); + MultiFab::Subtract(dEx, *Efield_cp[lev][0], 0, 0, Efield_cp[lev][0]->nComp(), ng); + MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, Efield_cp[lev][1]->nComp(), ng); + MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, Efield_cp[lev][2]->nComp(), ng); const int refinement_ratio = refRatio(lev-1)[0]; #ifdef _OPEMP @@ -199,8 +199,8 @@ WarpX::UpdateAuxilaryData () FArrayBox& aux = (*Efield_aux[lev][idim])[mfi]; FArrayBox& fp = (*Efield_fp[lev][idim])[mfi]; const Box& bx = aux.box(); - aux.copy(fp, bx, 0, bx, 0, 1); - aux.plus(efab[idim], bx, bx, 0, 0, 1); + aux.copy(fp, bx, 0, bx, 0, Efield_fp[lev][idim]->nComp()); + aux.plus(efab[idim], bx, bx, 0, 0, Efield_fp[lev][idim]->nComp()); } } } @@ -388,7 +388,7 @@ WarpX::SyncCurrent () // (potentially large) stencil of the multi-pass bilinear filter. j_fp[lev][idim].reset(new MultiFab(current_fp[lev][idim]->boxArray(), current_fp[lev][idim]->DistributionMap(), - 1, ng)); + current_fp[lev][idim]->nComp(), ng)); // Apply the filter to current_fp, store the result in j_fp. bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); // Then swap j_fp and current_fp @@ -405,7 +405,7 @@ WarpX::SyncCurrent () for (int idim = 0; idim < 3; ++idim) { j_cp[lev][idim].reset(new MultiFab(current_cp[lev][idim]->boxArray(), current_cp[lev][idim]->DistributionMap(), - 1, ng)); + current_cp[lev][idim]->nComp(), ng)); bilinear_filter.ApplyStencil(*j_cp[lev][idim], *current_cp[lev][idim]); std::swap(j_cp[lev][idim], current_cp[lev][idim]); } @@ -417,7 +417,7 @@ WarpX::SyncCurrent () for (int idim = 0; idim < 3; ++idim) { j_buf[lev][idim].reset(new MultiFab(current_buf[lev][idim]->boxArray(), current_buf[lev][idim]->DistributionMap(), - 1, ng)); + current_buf[lev][idim]->nComp(), ng)); bilinear_filter.ApplyStencil(*j_buf[lev][idim], *current_buf[lev][idim]); std::swap(*j_buf[lev][idim], *current_buf[lev][idim]); } @@ -448,16 +448,16 @@ WarpX::SyncCurrent () const MultiFab* ccz = current_cp[lev+1][2].get(); if (current_buf[lev+1][0]) { - MultiFab::Add(*current_buf[lev+1][0], *current_cp[lev+1][0], 0, 0, 1, ngsrc); - MultiFab::Add(*current_buf[lev+1][1], *current_cp[lev+1][1], 0, 0, 1, ngsrc); - MultiFab::Add(*current_buf[lev+1][2], *current_cp[lev+1][2], 0, 0, 1, ngsrc); + MultiFab::Add(*current_buf[lev+1][0], *current_cp[lev+1][0], 0, 0, current_cp[lev+1][0]->nComp(), ngsrc); + MultiFab::Add(*current_buf[lev+1][1], *current_cp[lev+1][1], 0, 0, current_cp[lev+1][1]->nComp(), ngsrc); + MultiFab::Add(*current_buf[lev+1][2], *current_cp[lev+1][2], 0, 0, current_cp[lev+1][2]->nComp(), ngsrc); ccx = current_buf[lev+1][0].get(); ccy = current_buf[lev+1][1].get(); ccz = current_buf[lev+1][2].get(); } - current_fp[lev][0]->copy(*ccx,0,0,1,ngsrc,ngdst,period,FabArrayBase::ADD); - current_fp[lev][1]->copy(*ccy,0,0,1,ngsrc,ngdst,period,FabArrayBase::ADD); - current_fp[lev][2]->copy(*ccz,0,0,1,ngsrc,ngdst,period,FabArrayBase::ADD); + current_fp[lev][0]->copy(*ccx,0,0,current_fp[lev][0]->nComp(),ngsrc,ngdst,period,FabArrayBase::ADD); + current_fp[lev][1]->copy(*ccy,0,0,current_fp[lev][1]->nComp(),ngsrc,ngdst,period,FabArrayBase::ADD); + current_fp[lev][2]->copy(*ccz,0,0,current_fp[lev][2]->nComp(),ngsrc,ngdst,period,FabArrayBase::ADD); } // Sum up coarse patch @@ -478,7 +478,7 @@ WarpX::SyncCurrent () // current_fp has right number of ghost cells. std::swap(j_fp[lev][idim], current_fp[lev][idim]); // Then copy the interior of j_fp to current_fp. - MultiFab::Copy(*current_fp[lev][idim], *j_fp[lev][idim], 0, 0, 1, 0); + MultiFab::Copy(*current_fp[lev][idim], *j_fp[lev][idim], 0, 0, j_fp[lev][idim]->nComp(), 0); // current_fp has right number of ghost cells and // correct filtered values here. } @@ -487,7 +487,7 @@ WarpX::SyncCurrent () { for (int idim = 0; idim < 3; ++idim) { std::swap(j_cp[lev][idim], current_cp[lev][idim]); - MultiFab::Copy(*current_cp[lev][idim], *j_cp[lev][idim], 0, 0, 1, 0); + MultiFab::Copy(*current_cp[lev][idim], *j_cp[lev][idim], 0, 0, j_cp[lev][idim]->nComp(), 0); } } for (int lev = 1; lev <= finest_level; ++lev) @@ -495,7 +495,7 @@ WarpX::SyncCurrent () for (int idim = 0; idim < 3; ++idim) { if (j_buf[lev][idim]) { std::swap(j_buf[lev][idim], current_buf[lev][idim]); - MultiFab::Copy(*current_buf[lev][idim], *j_buf[lev][idim], 0, 0, 1, 0); + MultiFab::Copy(*current_buf[lev][idim], *j_buf[lev][idim], 0, 0, j_buf[lev][idim]->nComp(), 0); } } } @@ -543,7 +543,7 @@ WarpX::SyncCurrent (const std::array<const amrex::MultiFab*,3>& fine, ffab.resize(fbx); fbx &= (*fine[idim])[mfi].box(); ffab.setVal(0.0); - ffab.copy((*fine[idim])[mfi], fbx, 0, fbx, 0, 1); + ffab.copy((*fine[idim])[mfi], fbx, 0, fbx, 0, fine[idim]->nComp()); WRPX_SYNC_CURRENT(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_ANYD((*crse[idim])[mfi]), BL_TO_FORTRAN_ANYD(ffab), @@ -731,7 +731,7 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) if (use_filter) { IntVect ng = j[idim]->nGrowVect(); ng += bilinear_filter.stencil_length_each_dir-1; - MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), 1, ng); + MultiFab jf(j[idim]->boxArray(), j[idim]->DistributionMap(), j[idim]->nComp(), ng); bilinear_filter.ApplyStencil(jf, *j[idim]); WarpXSumGuardCells(*(j[idim]), jf, period); } else { @@ -764,7 +764,7 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) const auto& period = Geom(lev).periodicity(); for (int idim = 0; idim < 3; ++idim) { MultiFab mf(current_fp[lev][idim]->boxArray(), - current_fp[lev][idim]->DistributionMap(), 1, 0); + current_fp[lev][idim]->DistributionMap(), current_fp[lev][idim]->nComp(), 0); mf.setVal(0.0); if (use_filter && current_buf[lev+1][idim]) { @@ -772,16 +772,16 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) IntVect ng = current_cp[lev+1][idim]->nGrowVect(); ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jfc(current_cp[lev+1][idim]->boxArray(), - current_cp[lev+1][idim]->DistributionMap(), 1, ng); + current_cp[lev+1][idim]->DistributionMap(), current_cp[lev+1][idim]->nComp(), ng); bilinear_filter.ApplyStencil(jfc, *current_cp[lev+1][idim]); // buffer patch of fine level MultiFab jfb(current_buf[lev+1][idim]->boxArray(), - current_buf[lev+1][idim]->DistributionMap(), 1, ng); + current_buf[lev+1][idim]->DistributionMap(), current_buf[lev+1][idim]->nComp(), ng); bilinear_filter.ApplyStencil(jfb, *current_buf[lev+1][idim]); - MultiFab::Add(jfb, jfc, 0, 0, 1, ng); - mf.ParallelAdd(jfb, 0, 0, 1, ng, IntVect::TheZeroVector(), period); + MultiFab::Add(jfb, jfc, 0, 0, jfc.nComp(), ng); + mf.ParallelAdd(jfb, 0, 0, jfb.nComp(), ng, IntVect::TheZeroVector(), period); WarpXSumGuardCells(*current_cp[lev+1][idim], jfc, period); } @@ -791,29 +791,29 @@ WarpX::AddCurrentFromFineLevelandSumBoundary (int lev) IntVect ng = current_cp[lev+1][idim]->nGrowVect(); ng += bilinear_filter.stencil_length_each_dir-1; MultiFab jf(current_cp[lev+1][idim]->boxArray(), - current_cp[lev+1][idim]->DistributionMap(), 1, ng); + current_cp[lev+1][idim]->DistributionMap(), current_cp[lev+1][idim]->nComp(), ng); bilinear_filter.ApplyStencil(jf, *current_cp[lev+1][idim]); - mf.ParallelAdd(jf, 0, 0, 1, ng, IntVect::TheZeroVector(), period); + mf.ParallelAdd(jf, 0, 0, jf.nComp(), ng, IntVect::TheZeroVector(), period); WarpXSumGuardCells(*current_cp[lev+1][idim], jf, period); } else if (current_buf[lev+1][idim]) // but no filter { MultiFab::Copy(*current_buf[lev+1][idim], - *current_cp [lev+1][idim], 0, 0, 1, + *current_cp [lev+1][idim], 0, 0, current_cp [lev+1][idim]->nComp(), current_cp[lev+1][idim]->nGrow()); - mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, 1, + mf.ParallelAdd(*current_buf[lev+1][idim], 0, 0, current_buf[lev+1][idim]->nComp(), current_buf[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); WarpXSumGuardCells(*(current_cp[lev+1][idim]), period); } else // no filter, no buffer { - mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, 1, + mf.ParallelAdd(*current_cp[lev+1][idim], 0, 0, current_cp[lev+1][idim]->nComp(), current_cp[lev+1][idim]->nGrowVect(), IntVect::TheZeroVector(), period); WarpXSumGuardCells(*(current_cp[lev+1][idim]), period); } - MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, 1, 0); + MultiFab::Add(*current_fp[lev][idim], mf, 0, 0, current_fp[lev][idim]->nComp(), 0); } NodalSyncJ(lev, PatchType::fine); NodalSyncJ(lev+1, PatchType::coarse); diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index 8d7873041..eb119d4a2 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -46,21 +46,21 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa { const IntVect& ng = Bfield_fp[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_fp[lev][idim]->boxArray(), - dm, 1, ng)); - pmf->Redistribute(*Bfield_fp[lev][idim], 0, 0, 1, ng); + dm, Bfield_fp[lev][idim]->nComp(), ng)); + pmf->Redistribute(*Bfield_fp[lev][idim], 0, 0, Bfield_fp[lev][idim]->nComp(), ng); Bfield_fp[lev][idim] = std::move(pmf); } { const IntVect& ng = Efield_fp[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_fp[lev][idim]->boxArray(), - dm, 1, ng)); - pmf->Redistribute(*Efield_fp[lev][idim], 0, 0, 1, ng); + dm, Efield_fp[lev][idim]->nComp(), ng)); + pmf->Redistribute(*Efield_fp[lev][idim], 0, 0, Efield_fp[lev][idim]->nComp(), ng); Efield_fp[lev][idim] = std::move(pmf); } { const IntVect& ng = current_fp[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(current_fp[lev][idim]->boxArray(), - dm, 1, ng)); + dm, current_fp[lev][idim]->nComp(), ng)); current_fp[lev][idim] = std::move(pmf); current_fp_owner_masks[lev][idim] = std::move(current_fp[lev][idim]->OwnerMask(period)); } @@ -68,7 +68,7 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa { const IntVect& ng = current_store[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(current_store[lev][idim]->boxArray(), - dm, 1, ng)); + dm, current_store[lev][idim]->nComp(), ng)); // no need to redistribute current_store[lev][idim] = std::move(pmf); } @@ -77,8 +77,8 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa if (F_fp[lev] != nullptr) { const IntVect& ng = F_fp[lev]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(F_fp[lev]->boxArray(), - dm, 1, ng)); - pmf->Redistribute(*F_fp[lev], 0, 0, 1, ng); + dm, F_fp[lev]->nComp(), ng)); + pmf->Redistribute(*F_fp[lev], 0, 0, F_fp[lev]->nComp(), ng); F_fp[lev] = std::move(pmf); } @@ -96,8 +96,8 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa if (lev == 0) { for (int idim = 0; idim < 3; ++idim) { - Bfield_aux[lev][idim].reset(new MultiFab(*Bfield_fp[lev][idim], amrex::make_alias, 0, 1)); - Efield_aux[lev][idim].reset(new MultiFab(*Efield_fp[lev][idim], amrex::make_alias, 0, 1)); + Bfield_aux[lev][idim].reset(new MultiFab(*Bfield_fp[lev][idim], amrex::make_alias, 0, Bfield_aux[lev][idim]->nComp())); + Efield_aux[lev][idim].reset(new MultiFab(*Efield_fp[lev][idim], amrex::make_alias, 0, Efield_aux[lev][idim]->nComp())); } } else { for (int idim=0; idim < 3; ++idim) @@ -105,15 +105,15 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa { const IntVect& ng = Bfield_aux[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_aux[lev][idim]->boxArray(), - dm, 1, ng)); - // pmf->Redistribute(*Bfield_aux[lev][idim], 0, 0, 1, ng); + dm, Bfield_aux[lev][idim]->nComp(), ng)); + // pmf->Redistribute(*Bfield_aux[lev][idim], 0, 0, Bfield_aux[lev][idim]->nComp(), ng); Bfield_aux[lev][idim] = std::move(pmf); } { const IntVect& ng = Efield_aux[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_aux[lev][idim]->boxArray(), - dm, 1, ng)); - // pmf->Redistribute(*Efield_aux[lev][idim], 0, 0, 1, ng); + dm, Efield_aux[lev][idim]->nComp(), ng)); + // pmf->Redistribute(*Efield_aux[lev][idim], 0, 0, Efield_aux[lev][idim]->nComp(), ng); Efield_aux[lev][idim] = std::move(pmf); } } @@ -127,21 +127,21 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa { const IntVect& ng = Bfield_cp[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_cp[lev][idim]->boxArray(), - dm, 1, ng)); - pmf->Redistribute(*Bfield_cp[lev][idim], 0, 0, 1, ng); + dm, Bfield_cp[lev][idim]->nComp(), ng)); + pmf->Redistribute(*Bfield_cp[lev][idim], 0, 0, Bfield_cp[lev][idim]->nComp(), ng); Bfield_cp[lev][idim] = std::move(pmf); } { const IntVect& ng = Efield_cp[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_cp[lev][idim]->boxArray(), - dm, 1, ng)); - pmf->Redistribute(*Efield_cp[lev][idim], 0, 0, 1, ng); + dm, Efield_cp[lev][idim]->nComp(), ng)); + pmf->Redistribute(*Efield_cp[lev][idim], 0, 0, Efield_cp[lev][idim]->nComp(), ng); Efield_cp[lev][idim] = std::move(pmf); } { const IntVect& ng = current_cp[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>( new MultiFab(current_cp[lev][idim]->boxArray(), - dm, 1, ng)); + dm, current_cp[lev][idim]->nComp(), ng)); current_cp[lev][idim] = std::move(pmf); current_cp_owner_masks[lev][idim] = std::move( current_cp[lev][idim]->OwnerMask(cperiod)); @@ -151,8 +151,8 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa if (F_cp[lev] != nullptr) { const IntVect& ng = F_cp[lev]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(F_cp[lev]->boxArray(), - dm, 1, ng)); - pmf->Redistribute(*F_cp[lev], 0, 0, 1, ng); + dm, F_cp[lev]->nComp(), ng)); + pmf->Redistribute(*F_cp[lev], 0, 0, F_cp[lev]->nComp(), ng); F_cp[lev] = std::move(pmf); } @@ -173,24 +173,24 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa { const IntVect& ng = Bfield_cax[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Bfield_cax[lev][idim]->boxArray(), - dm, 1, ng)); - // pmf->ParallelCopy(*Bfield_cax[lev][idim], 0, 0, 1, ng, ng); + dm, Bfield_cax[lev][idim]->nComp(), ng)); + // pmf->ParallelCopy(*Bfield_cax[lev][idim], 0, 0, Bfield_cax[lev][idim]->nComp(), ng, ng); Bfield_cax[lev][idim] = std::move(pmf); } if (Efield_cax[lev][idim]) { const IntVect& ng = Efield_cax[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(Efield_cax[lev][idim]->boxArray(), - dm, 1, ng)); - // pmf->ParallelCopy(*Efield_cax[lev][idim], 0, 0, 1, ng, ng); + dm, Efield_cax[lev][idim]->nComp(), ng)); + // pmf->ParallelCopy(*Efield_cax[lev][idim], 0, 0, Efield_cax[lev][idim]->nComp(), ng, ng); Efield_cax[lev][idim] = std::move(pmf); } if (current_buf[lev][idim]) { const IntVect& ng = current_buf[lev][idim]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(current_buf[lev][idim]->boxArray(), - dm, 1, ng)); - // pmf->ParallelCopy(*current_buf[lev][idim], 0, 0, 1, ng, ng); + dm, current_buf[lev][idim]->nComp(), ng)); + // pmf->ParallelCopy(*current_buf[lev][idim], 0, 0, current_buf[lev][idim]->nComp(), ng, ng); current_buf[lev][idim] = std::move(pmf); } } @@ -198,24 +198,24 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa { const IntVect& ng = charge_buf[lev]->nGrowVect(); auto pmf = std::unique_ptr<MultiFab>(new MultiFab(charge_buf[lev]->boxArray(), - dm, 1, ng)); - // pmf->ParallelCopy(*charge_buf[lev][idim], 0, 0, 1, ng, ng); + dm, charge_buf[lev]->nComp(), ng)); + // pmf->ParallelCopy(*charge_buf[lev][idim], 0, 0, charge_buf[lev]->nComp(), ng, ng); charge_buf[lev] = std::move(pmf); } if (current_buffer_masks[lev]) { const IntVect& ng = current_buffer_masks[lev]->nGrowVect(); auto pmf = std::unique_ptr<iMultiFab>(new iMultiFab(current_buffer_masks[lev]->boxArray(), - dm, 1, ng)); - // pmf->ParallelCopy(*current_buffer_masks[lev], 0, 0, 1, ng, ng); + dm, current_buffer_masks[lev]->nComp(), ng)); + // pmf->ParallelCopy(*current_buffer_masks[lev], 0, 0, current_buffer_masks[lev]->nComp(), ng, ng); current_buffer_masks[lev] = std::move(pmf); } if (gather_buffer_masks[lev]) { const IntVect& ng = gather_buffer_masks[lev]->nGrowVect(); auto pmf = std::unique_ptr<iMultiFab>(new iMultiFab(gather_buffer_masks[lev]->boxArray(), - dm, 1, ng)); - // pmf->ParallelCopy(*gather_buffer_masks[lev], 0, 0, 1, ng, ng); + dm, gather_buffer_masks[lev]->nComp(), ng)); + // pmf->ParallelCopy(*gather_buffer_masks[lev], 0, 0, gather_buffer_masks[lev]->nComp(), ng, ng); gather_buffer_masks[lev] = std::move(pmf); } } diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 6d618c096..f39a2a36d 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -299,7 +299,7 @@ MultiParticleContainer::GetChargeDensity (int lev, bool local) std::unique_ptr<MultiFab> rho = allcontainers[0]->GetChargeDensity(lev, true); for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) { std::unique_ptr<MultiFab> rhoi = allcontainers[i]->GetChargeDensity(lev, true); - MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow()); + MultiFab::Add(*rho, *rhoi, 0, 0, rho->nComp(), rho->nGrow()); } if (!local) { const Geometry& gm = allcontainers[0]->Geom(lev); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 212084e64..e1b012464 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1096,6 +1096,7 @@ PhysicalParticleContainer::FieldGather (int lev, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), + &WarpX::nmodes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); @@ -1396,6 +1397,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*bxfab), BL_TO_FORTRAN_ANYD(*byfab), BL_TO_FORTRAN_ANYD(*bzfab), + &WarpX::nmodes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); @@ -1493,6 +1495,7 @@ PhysicalParticleContainer::Evolve (int lev, BL_TO_FORTRAN_ANYD(*cbxfab), BL_TO_FORTRAN_ANYD(*cbyfab), BL_TO_FORTRAN_ANYD(*cbzfab), + &WarpX::nmodes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); } @@ -1825,6 +1828,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), + &WarpX::nmodes, &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index fd1b2dfb5..919cd9f0f 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -427,6 +427,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, BL_TO_FORTRAN_ANYD(bxfab), BL_TO_FORTRAN_ANYD(byfab), BL_TO_FORTRAN_ANYD(bzfab), + &WarpX::nmodes, &ll4symtry, &l_lower_order_in_v, &WarpX::do_nodal, &lvect_fieldgathe, &WarpX::field_gathering_algo); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 9791eee80..681758c45 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -324,9 +324,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, tby.grow(ngJ); tbz.grow(ngJ); - local_jx[thread_num].resize(tbx); - local_jy[thread_num].resize(tby); - local_jz[thread_num].resize(tbz); + local_jx[thread_num].resize(tbx, jx.nComp()); + local_jy[thread_num].resize(tby, jy.nComp()); + local_jz[thread_num].resize(tbz, jz.nComp()); jx_ptr = local_jx[thread_num].dataPtr(); jy_ptr = local_jy[thread_num].dataPtr(); @@ -407,6 +407,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &np_current, m_xp[thread_num].dataPtr(), m_yp[thread_num].dataPtr(), @@ -424,6 +425,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &xyzmin[0], &dx[0]); #endif } @@ -433,9 +435,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - jx[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); - jy[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, 1); - jz[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, 1); + jx[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, local_jx[thread_num].nComp()); + jy[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, local_jy[thread_num].nComp()); + jz[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, local_jz[thread_num].nComp()); BL_PROFILE_VAR_STOP(blp_accumulate); #endif @@ -465,9 +467,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, tby.grow(ngJ); tbz.grow(ngJ); - local_jx[thread_num].resize(tbx); - local_jy[thread_num].resize(tby); - local_jz[thread_num].resize(tbz); + local_jx[thread_num].resize(tbx, jx.nComp()); + local_jy[thread_num].resize(tby, jy.nComp()); + local_jz[thread_num].resize(tbz, jz.nComp()); jx_ptr = local_jx[thread_num].dataPtr(); jy_ptr = local_jy[thread_num].dataPtr(); @@ -549,6 +551,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &ncrse, m_xp[thread_num].dataPtr() +np_current, m_yp[thread_num].dataPtr() +np_current, @@ -567,6 +570,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr, &ngJ, jxntot.getVect(), jy_ptr, &ngJ, jyntot.getVect(), jz_ptr, &ngJ, jzntot.getVect(), + &WarpX::nmodes, &xyzmin[0], &dx[0]); #endif } @@ -576,9 +580,9 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - (*cjx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); - (*cjy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, 1); - (*cjz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, 1); + (*cjx)[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, local_jx[thread_num].nComp()); + (*cjy)[pti].atomicAdd(local_jy[thread_num], tby, tby, 0, 0, local_jy[thread_num].nComp()); + (*cjz)[pti].atomicAdd(local_jz[thread_num], tbz, tbz, 0, 0, local_jz[thread_num].nComp()); BL_PROFILE_VAR_STOP(blp_accumulate); #endif @@ -612,11 +616,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array<Real, 3>& xyzmin = xyzmin_tile; #ifdef AMREX_USE_GPU - data_ptr = (*rhomf)[pti].dataPtr(icomp); + data_ptr = (*rhomf)[pti].dataPtr(icomp*(rhomf->nComp()/2)); auto rholen = (*rhomf)[pti].length(); #else tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box); + local_rho[thread_num].resize(tile_box, rhomf->nComp()); data_ptr = local_rho[thread_num].dataPtr(); auto rholen = local_rho[thread_num].length(); @@ -655,7 +659,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); + (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp*(rhomf->nComp()/2), (rhomf->nComp()/2)); BL_PROFILE_VAR_STOP(blp_accumulate); #endif @@ -674,7 +678,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #else tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); tile_box.grow(ngRho); - local_rho[thread_num].resize(tile_box); + local_rho[thread_num].resize(tile_box, crhomf->nComp()); data_ptr = local_rho[thread_num].dataPtr(); auto rholen = local_rho[thread_num].length(); @@ -715,7 +719,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); - (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); + (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp*(crhomf->nComp()/2), (crhomf->nComp()/2)); BL_PROFILE_VAR_STOP(blp_accumulate); #endif @@ -770,7 +774,7 @@ WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, BoxArray coarsened_fine_BA = fine_BA; coarsened_fine_BA.coarsen(m_gdb->refRatio(lev)); - MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, 1, 0); + MultiFab coarsened_fine_data(coarsened_fine_BA, fine_dm, rho[lev+1]->nComp(), 0); coarsened_fine_data.setVal(0.0); IntVect ratio(AMREX_D_DECL(2, 2, 2)); // FIXME @@ -801,7 +805,7 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) const int ng = WarpX::nox; - auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,1,ng)); + auto rho = std::unique_ptr<MultiFab>(new MultiFab(nba,dm,WarpX::ncomps,ng)); rho->setVal(0.0); #ifdef _OPENMP @@ -829,7 +833,7 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); const std::array<Real, 3>& xyzmin = xyzmin_tile; tile_box.grow(ng); - rho_loc.resize(tile_box); + rho_loc.resize(tile_box, rho->nComp()); rho_loc = 0.0; data_ptr = rho_loc.dataPtr(); auto rholen = rho_loc.length(); diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 3c1a930b3..10e5ed8dd 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -10,11 +10,14 @@ namespace { - double** getMultiFabPointers(const amrex::MultiFab& mf, int *num_boxes, int *ngrow, int **shapes) + double** getMultiFabPointers(const amrex::MultiFab& mf, int *num_boxes, int *ncomps, int *ngrow, int **shapes) { + *ncomps = mf.nComp(); *ngrow = mf.nGrow(); *num_boxes = mf.local_size(); - *shapes = (int*) malloc(AMREX_SPACEDIM * (*num_boxes) * sizeof(int)); + int shapesize = AMREX_SPACEDIM; + if (mf.nComp() > 1) shapesize += 1; + *shapes = (int*) malloc(shapesize * (*num_boxes) * sizeof(int)); double** data = (double**) malloc((*num_boxes) * sizeof(double*)); int i = 0; @@ -24,8 +27,9 @@ namespace for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi, ++i ) { data[i] = (double*) mf[mfi].dataPtr(); for (int j = 0; j < AMREX_SPACEDIM; ++j) { - (*shapes)[AMREX_SPACEDIM*i+j] = mf[mfi].box().length(j); + (*shapes)[shapesize*i+j] = mf[mfi].box().length(j); } + if (mf.nComp() > 1) (*shapes)[shapesize*i+2] = mf.nComp(); } return data; } @@ -197,9 +201,9 @@ extern "C" } double** warpx_getEfield(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getEfield(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getEfieldLoVects(int lev, int direction, @@ -209,9 +213,9 @@ extern "C" } double** warpx_getEfieldCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getEfieldCPLoVects(int lev, int direction, @@ -221,9 +225,9 @@ extern "C" } double** warpx_getEfieldFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getEfieldFPLoVects(int lev, int direction, @@ -233,9 +237,9 @@ extern "C" } double** warpx_getBfield(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getBfield(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getBfieldLoVects(int lev, int direction, @@ -245,9 +249,9 @@ extern "C" } double** warpx_getBfieldCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getBfieldCPLoVects(int lev, int direction, @@ -257,9 +261,9 @@ extern "C" } double** warpx_getBfieldFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getBfieldFPLoVects(int lev, int direction, @@ -269,9 +273,9 @@ extern "C" } double** warpx_getCurrentDensity(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getcurrent(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getCurrentDensityLoVects(int lev, int direction, @@ -281,9 +285,9 @@ extern "C" } double** warpx_getCurrentDensityCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getCurrentDensityCPLoVects(int lev, int direction, @@ -293,9 +297,9 @@ extern "C" } double** warpx_getCurrentDensityFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getCurrentDensityFPLoVects(int lev, int direction, diff --git a/Source/WarpX.H b/Source/WarpX.H index 35b072142..903e01770 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -91,6 +91,10 @@ public: static long noy; static long noz; + // Number of modes for the RZ multimode version + static long nmodes; + static long ncomps; + static bool use_fdtd_nci_corr; static int l_lower_order_in_v; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 3d7f7dcc5..93653040d 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -41,6 +41,9 @@ long WarpX::field_gathering_algo = 1; long WarpX::particle_pusher_algo = 0; int WarpX::maxwell_fdtd_solver_id = 0; +long WarpX::nmodes = 1; +long WarpX::ncomps = 1; + long WarpX::nox = 1; long WarpX::noy = 1; long WarpX::noz = 1; @@ -463,6 +466,10 @@ WarpX::ReadParameters () // Use same shape factors in all directions, for gathering l_lower_order_in_v = false; } + + // Only needs to be set with WARPX_RZ, otherwise defaults to 1. + pp.query("nmodes", nmodes); + } { @@ -682,20 +689,31 @@ void WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, int ngF) { + +#if defined WARPX_RZ + if (nmodes > 1) { + // There is a real and imaginary component for each mode + ncomps = nmodes*2; + } else { + // With only mode 0, only reals are used + ncomps = 1; + } +#endif + // // The fine patch // - Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,1,ngE)); - Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,1,ngE)); - Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,1,ngE)); + Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE)); - Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,1,ngE)); - Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,1,ngE)); - Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,1,ngE)); + Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); - current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ)); - current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); - current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); + current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ)); + current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ)); + current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ)); const auto& period = Geom(lev).periodicity(); current_fp_owner_masks[lev][0] = std::move(current_fp[lev][0]->OwnerMask(period)); @@ -704,25 +722,25 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (do_dive_cleaning || plot_rho) { - rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho)); + rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); rho_fp_owner_masks[lev] = std::move(rho_fp[lev]->OwnerMask(period)); } if (do_subcycling == 1 && lev == 0) { - current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,1,ngJ)); - current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,1,ngJ)); - current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,1,ngJ)); + current_store[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ)); + current_store[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ)); + current_store[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ)); } if (do_dive_cleaning) { - F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,1, ngF)); + F_fp[lev].reset (new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,ncomps, ngF)); } #ifdef WARPX_USE_PSATD else { - rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2,ngRho)); + rho_fp[lev].reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); rho_fp_owner_masks[lev] = std::move(rho_fp[lev]->OwnerMask(period)); } #endif @@ -733,19 +751,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (lev == 0) { for (int idir = 0; idir < 3; ++idir) { - Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, 1)); - Bfield_aux[lev][idir].reset(new MultiFab(*Bfield_fp[lev][idir], amrex::make_alias, 0, 1)); + Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, ncomps)); + Bfield_aux[lev][idir].reset(new MultiFab(*Bfield_fp[lev][idir], amrex::make_alias, 0, ncomps)); } } else { - Bfield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,1,ngE)); - Bfield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,1,ngE)); - Bfield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,1,ngE)); + Bfield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE)); - Efield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,1,ngE)); - Efield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,1,ngE)); - Efield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,1,ngE)); + Efield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); } // @@ -757,19 +775,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm cba.coarsen(refRatio(lev-1)); // Create the MultiFabs for B - Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,1,ngE)); - Bfield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,1,ngE)); - Bfield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,1,ngE)); + Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE)); // Create the MultiFabs for E - Efield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,1,ngE)); - Efield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,1,ngE)); - Efield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE)); + Efield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); // Create the MultiFabs for the current - current_cp[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ)); - current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); - current_cp[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); + current_cp[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ)); + current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ)); + current_cp[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ)); const auto& cperiod = Geom(lev-1).periodicity(); current_cp_owner_masks[lev][0] = std::move(current_cp[lev][0]->OwnerMask(cperiod)); @@ -777,17 +795,17 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm current_cp_owner_masks[lev][2] = std::move(current_cp[lev][2]->OwnerMask(cperiod)); if (do_dive_cleaning || plot_rho){ - rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); + rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); rho_cp_owner_masks[lev] = std::move(rho_cp[lev]->OwnerMask(cperiod)); } if (do_dive_cleaning) { - F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,1, ngF)); + F_cp[lev].reset (new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,ncomps, ngF)); } #ifdef WARPX_USE_PSATD else { - rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); + rho_cp[lev].reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); rho_cp_owner_masks[lev] = std::move(rho_cp[lev]->OwnerMask(cperiod)); } #endif @@ -803,28 +821,28 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (n_field_gather_buffer > 0) { // Create the MultiFabs for B - Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,1,ngE)); - Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,1,ngE)); - Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,1,ngE)); + Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE)); // Create the MultiFabs for E - Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,1,ngE)); - Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,1,ngE)); - Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,1,ngE)); + Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); - gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); + gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) ); // Gather buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. } if (n_current_deposition_buffer > 0) { - current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,1,ngJ)); - current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,1,ngJ)); - current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,1,ngJ)); + current_buf[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ)); + current_buf[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ)); + current_buf[lev][2].reset( new MultiFab(amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ)); if (do_dive_cleaning || plot_rho) { - charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2,ngRho)); + charge_buf[lev].reset( new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()),dm,2*ncomps,ngRho)); } - current_buffer_masks[lev].reset( new iMultiFab(ba, dm, 1, 1) ); + current_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) ); // Current buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. } |