diff options
Diffstat (limited to 'Source')
50 files changed, 1078 insertions, 136 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 441acb743..47e4cce62 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -184,7 +184,9 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { +#if (AMREX_SPACEDIM >= 2) int jdim = (idim+1) % AMREX_SPACEDIM; +#endif #if (AMREX_SPACEDIM == 3) int kdim = (idim+2) % AMREX_SPACEDIM; #endif @@ -198,6 +200,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n { direct_faces.push_back(kv.first); } +#if (AMREX_SPACEDIM >= 2) else if (amrex::grow(grid_box, jdim, ncell).intersects(box)) { side_faces.push_back(kv.first); @@ -227,8 +230,10 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n { corners.push_back(kv.first); } +#endif } +#if (AMREX_SPACEDIM >= 2) for (auto gid : corners) { const Box& grid_box = grids[gid]; @@ -261,6 +266,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n amrex::Abort("SigmaBox::SigmaBox(): corners, how did this happen?\n"); } } +#endif #if (AMREX_SPACEDIM == 3) for (auto gid : side_side_edges) @@ -302,6 +308,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n } #endif +#if (AMREX_SPACEDIM >= 2) for (auto gid : side_faces) { const Box& grid_box = grids[gid]; @@ -317,6 +324,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n"); } } +#endif for (auto gid : direct_faces) { @@ -369,7 +377,12 @@ SigmaBox::ComputePMLFactorsB (const Real* a_dx, Real dt) N[idim] = sigma_star[idim].size(); dx[idim] = a_dx[idim]; } - amrex::ParallelFor(amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), + amrex::ParallelFor( +#if (AMREX_SPACEDIM >= 2) + amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), +#else + N[0], +#endif [=] AMREX_GPU_DEVICE (int i) noexcept { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { @@ -398,7 +411,12 @@ SigmaBox::ComputePMLFactorsE (const Real* a_dx, Real dt) N[idim] = sigma[idim].size(); dx[idim] = a_dx[idim]; } - amrex::ParallelFor(amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), + amrex::ParallelFor( +#if (AMREX_SPACEDIM >= 2) + amrex::max(AMREX_D_DECL(N[0],N[1],N[2])), +#else + N[0], +#endif [=] AMREX_GPU_DEVICE (int i) noexcept { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { @@ -524,6 +542,8 @@ PML::PML (const int lev, const BoxArray& grid_ba, const DistributionMapping& /*g IntVect ngFFT = IntVect(ngFFt_x, ngFFt_y, ngFFt_z); #elif (AMREX_SPACEDIM == 2) IntVect ngFFT = IntVect(ngFFt_x, ngFFt_z); +#elif (AMREX_SPACEDIM == 1) + IntVect ngFFT = IntVect(ngFFt_z); #endif // Set the number of guard cells to the maximum of each field diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index 9c7435b6c..4458f1e2a 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -80,6 +80,9 @@ PEC::ApplyPECtoEfield (std::array<amrex::MultiFab*, 3> Efield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 0; PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -90,6 +93,9 @@ PEC::ApplyPECtoEfield (std::array<amrex::MultiFab*, 3> Efield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 1; PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -100,6 +106,9 @@ PEC::ApplyPECtoEfield (std::array<amrex::MultiFab*, 3> Efield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 2; PEC::SetEfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -163,6 +172,9 @@ PEC::ApplyPECtoBfield (std::array<amrex::MultiFab*, 3> Bfield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 0; PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -173,6 +185,9 @@ PEC::ApplyPECtoBfield (std::array<amrex::MultiFab*, 3> Bfield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 1; PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, @@ -183,6 +198,9 @@ PEC::ApplyPECtoBfield (std::array<amrex::MultiFab*, 3> Bfield, const int lev, #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) amrex::ignore_unused(k); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(j,k); +#endif amrex::IntVect iv(AMREX_D_DECL(i,j,k)); const int icomp = 2; PEC::SetBfieldOnPEC(icomp, domain_lo, domain_hi, iv, n, diff --git a/Source/BoundaryConditions/WarpX_PML_kernels.H b/Source/BoundaryConditions/WarpX_PML_kernels.H index ca6b594ac..c434e09f3 100644 --- a/Source/BoundaryConditions/WarpX_PML_kernels.H +++ b/Source/BoundaryConditions/WarpX_PML_kernels.H @@ -27,6 +27,12 @@ void warpx_damp_pml_ex (int i, int j, int k, Array4<Real> const& Ex, int xlo, int ylo, int zlo, const bool dive_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Ex, Ex_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -98,6 +104,12 @@ void warpx_damp_pml_ey (int i, int j, int k, Array4<Real> const& Ey, int xlo, int ylo, int zlo, const bool dive_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Ey, Ey_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo, dive_cleaning); @@ -166,6 +178,12 @@ void warpx_damp_pml_ez (int i, int j, int k, Array4<Real> const& Ez, int xlo, int ylo, int zlo, const bool dive_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Ez, Ez_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, dive_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -237,6 +255,12 @@ void warpx_damp_pml_bx (int i, int j, int k, Array4<Real> const& Bx, int xlo, int ylo, int zlo, const bool divb_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Bx, Bx_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -308,6 +332,12 @@ void warpx_damp_pml_by (int i, int j, int k, Array4<Real> const& By, int xlo, int ylo, int zlo, const bool divb_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, By, By_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo, divb_cleaning); @@ -376,6 +406,12 @@ void warpx_damp_pml_bz (int i, int j, int k, Array4<Real> const& Bz, int xlo, int ylo, int zlo, const bool divb_cleaning) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, Bz, Bz_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo, divb_cleaning); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); @@ -446,6 +482,12 @@ void warpx_damp_pml_scalar (int i, int j, int k, Array4<Real> const& arr, const Real* const sigma_star_fac_z, int xlo, int ylo, int zlo) { +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(i, j, k, arr, arr_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z, + sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, xlo, ylo, zlo); + amrex::Abort("PML not implemented in Cartesian 1D geometry"); +#endif + #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo); diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index 2f2231c0c..236516f07 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -301,6 +301,7 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) / dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]) ) ); // Take the max of 0 and num_zcells_lab int Nz_lab = std::max( 0, num_zcells_lab ); +#if (AMREX_SPACEDIM >= 2) // Number of lab-frame cells in x-direction at level, lev const int num_xcells_lab = static_cast<int>( floor ( ( diag_dom.hi(0) - diag_dom.lo(0) ) @@ -308,6 +309,7 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) ) ); // Take the max of 0 and num_ycells_lab int Nx_lab = std::max( 0, num_xcells_lab); +#endif #if (AMREX_SPACEDIM == 3) // Number of lab-frame cells in the y-direction at level, lev const int num_ycells_lab = static_cast<int>( floor ( @@ -317,8 +319,10 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) // Take the max of 0 and num_xcells_lab int Ny_lab = std::max( 0, num_ycells_lab ); m_snapshot_ncells_lab[i_buffer] = {Nx_lab, Ny_lab, Nz_lab}; -#else +#elif (AMREX_SPACEDIM == 2) m_snapshot_ncells_lab[i_buffer] = {Nx_lab, Nz_lab}; +#else + m_snapshot_ncells_lab[i_buffer] = amrex::IntVect({Nz_lab}); #endif } diff --git a/Source/Diagnostics/BackTransformedDiagnostic.cpp b/Source/Diagnostics/BackTransformedDiagnostic.cpp index 17b59d3ae..c2d255d06 100644 --- a/Source/Diagnostics/BackTransformedDiagnostic.cpp +++ b/Source/Diagnostics/BackTransformedDiagnostic.cpp @@ -160,8 +160,10 @@ namespace // Create a 3D, nx x ny x nz dataspace. #if (AMREX_SPACEDIM == 3) hsize_t dims[3] = {nx, ny, nz}; -#else +#elif (AMREX_SPACEDIM == 2) hsize_t dims[3] = {nx, nz}; +#else + hsize_t dims[3] = {nz}; #endif hid_t grid_space = H5Screate_simple(AMREX_SPACEDIM, dims, NULL); @@ -422,11 +424,15 @@ namespace shift[0] = lo_x; shift[1] = lo_y; shift[2] = lo_z; -#else +#elif (AMREX_SPACEDIM == 2) hsize_t slab_offsets[2], slab_dims[2]; int shift[2]; shift[0] = lo_x; shift[1] = lo_z; +#else + hsize_t slab_offsets[1], slab_dims[1]; + int shift[1]; + shift[0] = lo_z; #endif hid_t slab_dataspace; @@ -584,13 +590,19 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, m_dz_lab_ = PhysConst::c * m_dt_boost_ * m_inv_beta_boost_ * m_inv_gamma_boost_; m_inv_dz_lab_ = 1.0_rt / m_dz_lab_; int Nz_lab = static_cast<unsigned>((zmax_lab - zmin_lab) * m_inv_dz_lab_); +#if (AMREX_SPACEDIM >= 2) int Nx_lab = geom.Domain().length(0); +#endif #if (AMREX_SPACEDIM == 3) int Ny_lab = geom.Domain().length(1); IntVect prob_ncells_lab = {Nx_lab, Ny_lab, Nz_lab}; -#else +#elif (AMREX_SPACEDIM == 2) // Ny_lab = 1; IntVect prob_ncells_lab = {Nx_lab, Nz_lab}; +#else + // Nx_lab = 1; + // Ny_lab = 1; + IntVect prob_ncells_lab({Nz_lab}); #endif writeMetaData(); @@ -640,7 +652,6 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, for (int i = 0; i < N_slice_snapshots; ++i) { - IntVect slice_ncells_lab ; // To construct LabFrameSlice(), the location of lo() and hi() of the // reduced diag is computed using the user-defined values of the @@ -657,6 +668,7 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, ( (1._rt+m_beta_boost_)*m_gamma_boost_); auto Nz_slice_lab = static_cast<int>( (zmax_slice_lab - zmin_slice_lab) * m_inv_dz_lab_); +#if (AMREX_SPACEDIM >= 2) auto Nx_slice_lab = static_cast<int>( (current_slice_hi[0] - current_slice_lo[0] ) / geom.CellSize(0)); @@ -664,6 +676,7 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, // if the x-dimension is reduced, increase total_cells by 1 // to be consistent with the number of cells created for the output. if (Nx_lab != Nx_slice_lab) Nx_slice_lab++; +#endif #if (AMREX_SPACEDIM == 3) auto Ny_slice_lab = static_cast<int>( (current_slice_hi[1] - current_slice_lo[1]) / @@ -672,9 +685,11 @@ BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab, // if the y-dimension is reduced, increase total_cells by 1 // to be consistent with the number of cells created for the output. if (Ny_lab != Ny_slice_lab) Ny_slice_lab++; - slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab}; + amrex::IntVect slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab}; +#elif (AMREX_SPACEDIM == 2) + amrex::IntVect slice_ncells_lab = {Nx_slice_lab, Nz_slice_lab}; #else - slice_ncells_lab = {Nx_slice_lab, Nz_slice_lab}; + amrex::IntVect slice_ncells_lab({Nz_slice_lab}); #endif IntVect slice_lo(AMREX_D_DECL(0,0,0)); @@ -1213,7 +1228,11 @@ createLabFrameDirectories() { const auto lo = lbound(m_buff_box_); for (int comp = 0; comp < m_ncomp_to_dump_; ++comp) { output_create_field(m_file_name, m_mesh_field_names_[comp], +#if ( AMREX_SPACEDIM >= 2 ) m_prob_ncells_lab_[0], +#else + 1, +#endif #if ( AMREX_SPACEDIM == 3 ) m_prob_ncells_lab_[1], #else @@ -1387,8 +1406,10 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, const int icomp = field_map_ptr[n]; #if (AMREX_SPACEDIM == 3) buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp); -#else +#elif (AMREX_SPACEDIM == 2) buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp); +#else + buf_arr(k_lab,j,k,n) = tmp_arr(i,j,k,icomp); #endif } ); @@ -1424,8 +1445,10 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, const int icomp = field_map_ptr[n]; #if (AMREX_SPACEDIM == 3) buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp); -#else +#elif (AMREX_SPACEDIM == 2) buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp); +#else + buf_arr(k_lab,j,k,n) = tmp_arr(i,j,k,icomp); #endif }); } @@ -1531,8 +1554,10 @@ AddPartDataToParticleBuffer( int* const AMREX_RESTRICT IndexLocation = IndexForPartCopy.dataPtr(); // Compute extent of the reduced domain +/- user-defined physical width +#if (AMREX_SPACEDIM >= 2) Real const xmin = m_diag_domain_lab_.lo(0)-m_particle_slice_dx_lab_; Real const xmax = m_diag_domain_lab_.hi(0)+m_particle_slice_dx_lab_; +#endif #if (AMREX_SPACEDIM == 3) Real const ymin = m_diag_domain_lab_.lo(1)-m_particle_slice_dx_lab_; Real const ymax = m_diag_domain_lab_.hi(1)+m_particle_slice_dx_lab_; @@ -1544,8 +1569,11 @@ AddPartDataToParticleBuffer( [=] AMREX_GPU_DEVICE(int i) { Flag[i] = 0; +#if (AMREX_SPACEDIM >= 2) if ( x_temp[i] >= (xmin) && - x_temp[i] <= (xmax) ) { + x_temp[i] <= (xmax) ) +#endif + { #if (AMREX_SPACEDIM == 3) if (y_temp[i] >= (ymin) && y_temp[i] <= (ymax) ) diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index 9df599fa7..6fef39ce8 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -109,8 +109,10 @@ Diagnostics::BaseReadParameters () if (warpx.do_moving_window) { #if (AMREX_SPACEDIM == 3) amrex::Vector<int> dim_map {0, 1, 2}; -#else +#elif (AMREX_SPACEDIM == 2) amrex::Vector<int> dim_map {0, 2}; +#else + amrex::Vector<int> dim_map {2}; #endif if (warpx.boost_direction[ dim_map[warpx.moving_window_dir] ] == 1) { // Convert user-defined lo and hi for diagnostics to account for boosted-frame diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index 3490c4755..4bcf233fb 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -496,6 +496,11 @@ FullDiagnostics::MovingWindowAndGalileanDomainShift (int step) new_lo[1] = current_lo[1] + warpx.m_galilean_shift[2]; new_hi[1] = current_hi[1] + warpx.m_galilean_shift[2]; } +#elif (AMREX_SPACEDIM == 1 ) + { + new_lo[0] = current_lo[0] + warpx.m_galilean_shift[2]; + new_hi[0] = current_hi[0] + warpx.m_galilean_shift[2]; + } #endif // Update RealBox of geometry with galilean-shifted boundary for (int lev = 0; lev < nmax_lev; ++lev) { diff --git a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp index 4896176c2..0c2cf7e95 100644 --- a/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp +++ b/Source/Diagnostics/ReducedDiags/BeamRelevant.cpp @@ -60,6 +60,16 @@ BeamRelevant::BeamRelevant (std::string rd_name) // 12,13: emittance x,z // 14: charge m_data.resize(15, 0.0_rt); +#elif (defined WARPX_DIM_1D_Z) + // 0 : mean z + // 1,2,3 : mean px,py,pz + // 4 : gamma + // 5 : rms z + // 6,7,8 : rms px,py,pz + // 9 : rms gamma + // 10 : emittance z + // 11 : charge + m_data.resize(12, 0.0_rt); #endif if (ParallelDescriptor::IOProcessor()) @@ -112,6 +122,23 @@ BeamRelevant::BeamRelevant (std::string rd_name) ofs << "[" << c++ << "]emittance_x(m)"; ofs << m_sep; ofs << "[" << c++ << "]emittance_z(m)"; ofs << m_sep; ofs << "[" << c++ << "]charge(C)"; ofs << std::endl; +#elif (defined WARPX_DIM_1D_Z) + int c = 0; + ofs << "#"; + ofs << "[" << c++ << "]step()"; ofs << m_sep; + ofs << "[" << c++ << "]time(s)"; ofs << m_sep; + ofs << "[" << c++ << "]z_mean(m)"; ofs << m_sep; + ofs << "[" << c++ << "]px_mean(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]py_mean(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]pz_mean(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]gamma_mean()"; ofs << m_sep; + ofs << "[" << c++ << "]z_rms(m)"; ofs << m_sep; + ofs << "[" << c++ << "]px_rms(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]py_rms(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]pz_rms(kg*m/s)"; ofs << m_sep; + ofs << "[" << c++ << "]gamma_rms()"; ofs << m_sep; + ofs << "[" << c++ << "]emittance_z(m)"; ofs << m_sep; + ofs << "[" << c++ << "]charge(C)"; ofs << std::endl; #endif // close file ofs.close(); @@ -143,6 +170,8 @@ void BeamRelevant::ComputeDiags (int step) int const index_z = 2; #elif (defined WARPX_DIM_XZ || defined WARPX_DIM_RZ) int const index_z = 1; +#elif (defined WARPX_DIM_1D_Z) + int const index_z = 0; #endif // loop over species @@ -250,10 +279,16 @@ void BeamRelevant::ComputeDiags (int step) const ParticleReal p_pos0 = p.pos(0); const ParticleReal p_w = p.rdata(PIdx::w); -#if (defined WARPX_DIM_RZ) +#if (defined WARPX_DIM_1D_Z) + const ParticleReal p_x = 0.0; + const ParticleReal p_y = 0.0; +#elif (defined WARPX_DIM_RZ) const ParticleReal p_theta = p.rdata(PIdx::theta); const ParticleReal p_x = p_pos0*std::cos(p_theta); const ParticleReal p_y = p_pos0*std::sin(p_theta); +#elif (defined WARPX_DIM_XZ) + const ParticleReal p_x = p_pos0; + const ParticleReal p_y = 0.0; #else const ParticleReal p_pos1 = p.pos(1); const ParticleReal p_x = p_pos0; @@ -351,6 +386,20 @@ void BeamRelevant::ComputeDiags (int step) m_data[13] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c; m_data[14] = charge; amrex::ignore_unused(y_mean, y_ms, yuy); +#elif (defined WARPX_DIM_1D_Z) + m_data[0] = z_mean; + m_data[1] = ux_mean * m; + m_data[2] = uy_mean * m; + m_data[3] = uz_mean * m; + m_data[4] = gm_mean; + m_data[5] = std::sqrt(z_ms); + m_data[6] = std::sqrt(ux_ms) * m; + m_data[7] = std::sqrt(uy_ms) * m; + m_data[8] = std::sqrt(uz_ms) * m; + m_data[9] = std::sqrt(gm_ms); + m_data[10] = std::sqrt(z_ms*uz_ms-zuz*zuz) / PhysConst::c; + m_data[11] = charge; + amrex::ignore_unused(x_mean, x_ms, xux, y_mean, y_ms, yuy); #endif } // end loop over species diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp index 378a11989..d83669ce5 100644 --- a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp @@ -99,7 +99,9 @@ void FieldEnergy::ComputeDiags (int step) // get cell size Geometry const & geom = warpx.Geom(lev); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + auto dV = geom.CellSize(0); +#elif (AMREX_SPACEDIM == 2) auto dV = geom.CellSize(0) * geom.CellSize(1); #elif (AMREX_SPACEDIM == 3) auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2); diff --git a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp index c835feafd..1a0c9029d 100644 --- a/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldMomentum.cpp @@ -186,7 +186,9 @@ void FieldMomentum::ComputeDiags (int step) // Get cell size amrex::Geometry const & geom = warpx.Geom(lev); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + auto dV = geom.CellSize(0); +#elif (AMREX_SPACEDIM == 2) auto dV = geom.CellSize(0) * geom.CellSize(1); #elif (AMREX_SPACEDIM == 3) auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2); diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index 648617719..adaa3d722 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -189,7 +189,9 @@ bool FieldProbe::ProbeInDomain () const * value in a position array will be the y value, but in the case of 2D, prob_lo[1] * and prob_hi[1] refer to z. This is a result of warpx.Geom(lev). */ -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + return z_probe >= prob_lo[1] && z_probe < prob_hi[1]; +#elif (AMREX_SPACEDIM == 2) return x_probe >= prob_lo[0] && x_probe < prob_hi[0] && z_probe >= prob_lo[1] && z_probe < prob_hi[1]; #else @@ -251,7 +253,10 @@ void FieldProbe::ComputeDiags (int step) { const auto cell_size = gm.CellSizeArray(); const int i_probe = static_cast<int>(amrex::Math::floor((x_probe - prob_lo[0]) / cell_size[0])); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const int j_probe = 0; + const int k_probe = 0; +#elif (AMREX_SPACEDIM == 2) const int j_probe = static_cast<int>(amrex::Math::floor((z_probe - prob_lo[1]) / cell_size[1])); const int k_probe = 0; #elif(AMREX_SPACEDIM == 3) diff --git a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp index c30d5820f..f092d53e3 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp @@ -185,6 +185,8 @@ void ParticleExtrema::ComputeDiags (int step) int const index_z = 2; #elif (defined WARPX_DIM_XZ || defined WARPX_DIM_RZ) int const index_z = 1; +#elif (defined WARPX_DIM_1D_Z) + int const index_z = 0; #endif // loop over species @@ -211,6 +213,8 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::cos(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMin(xmin); +#elif (defined WARPX_DIM_1D_Z) + Real xmin = 0.0_rt; #else Real xmin = ReduceMin( myspc, [=] AMREX_GPU_HOST_DEVICE (const PType& p) @@ -224,6 +228,8 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::cos(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMax(xmax); +#elif (defined WARPX_DIM_1D_Z) + Real xmax = 0.0_rt; #else Real xmax = ReduceMax( myspc, [=] AMREX_GPU_HOST_DEVICE (const PType& p) @@ -237,7 +243,7 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::sin(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMin(ymin); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) Real ymin = 0.0_rt; #else Real ymin = ReduceMin( myspc, @@ -252,7 +258,7 @@ void ParticleExtrema::ComputeDiags (int step) [=] AMREX_GPU_HOST_DEVICE (const PType& p) { return p.pos(0)*std::sin(p.rdata(PIdx::theta)); }); ParallelDescriptor::ReduceRealMax(ymax); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) Real ymax = 0.0_rt; #else Real ymax = ReduceMax( myspc, diff --git a/Source/Evolve/WarpXComputeDt.cpp b/Source/Evolve/WarpXComputeDt.cpp index 392ef81b5..08d12dde6 100644 --- a/Source/Evolve/WarpXComputeDt.cpp +++ b/Source/Evolve/WarpXComputeDt.cpp @@ -38,7 +38,9 @@ WarpX::ComputeDt () if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) { // Computation of dt for spectral algorithm // (determined by the minimum cell size in all directions) -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + deltat = cfl * dx[0] / PhysConst::c; +#elif (AMREX_SPACEDIM == 2) deltat = cfl * std::min(dx[0], dx[1]) / PhysConst::c; #else deltat = cfl * std::min(dx[0], std::min(dx[1], dx[2])) / PhysConst::c; @@ -84,10 +86,13 @@ WarpX::PrintDtDxDyDz () for (int lev=0; lev <= max_level; lev++) { const amrex::Real* dx_lev = geom[lev].CellSize(); amrex::Print() << "Level " << lev << ": dt = " << dt[lev] +#if (defined WARPX_DIM_1D_Z) + << " ; dz = " << dx_lev[0] << '\n'; +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) << " ; dx = " << dx_lev[0] -#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) << " ; dz = " << dx_lev[1] << '\n'; #elif (defined WARPX_DIM_3D) + << " ; dx = " << dx_lev[0] << " ; dy = " << dx_lev[1] << " ; dz = " << dx_lev[2] << '\n'; #endif diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index d04a54ae3..770887cfe 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -253,7 +253,6 @@ WarpX::Evolve (int numsteps) bool move_j = is_synchronized; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. - int num_moved = MoveWindow(step+1, move_j); mypc->ContinuousFluxInjection(dt[0]); diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index 5f33bbfb4..582890aa8 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -400,7 +400,9 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab> // Set the value of beta amrex::Array<amrex::Real,AMREX_SPACEDIM> beta_solver = -# if (AMREX_SPACEDIM==2) +# if (AMREX_SPACEDIM==1) + {{ beta[2] }}; // beta_x and beta_z +# elif (AMREX_SPACEDIM==2) {{ beta[0], beta[2] }}; // beta_x and beta_z # else {{ beta[0], beta[1], beta[2] }}; @@ -463,7 +465,13 @@ WarpX::computePhiCartesian (const amrex::Vector<std::unique_ptr<amrex::MultiFab> if (do_electrostatic == ElectrostaticSolverAlgo::LabFrame) { for (int lev = 0; lev <= max_level; ++lev) { -#if (AMREX_SPACEDIM==2) +#if (AMREX_SPACEDIM==1) + mlmg.getGradSolution( + {amrex::Array<amrex::MultiFab*,1>{ + get_pointer_Efield_fp(lev, 2) + }} + ); +#elif (AMREX_SPACEDIM==2) mlmg.getGradSolution( {amrex::Array<amrex::MultiFab*,2>{ get_pointer_Efield_fp(lev, 0),get_pointer_Efield_fp(lev, 2) @@ -578,21 +586,28 @@ WarpX::computeE (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > #endif for ( MFIter mfi(*phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - const Real inv_dx = 1._rt/dx[0]; #if (AMREX_SPACEDIM == 3) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dy = 1._rt/dx[1]; const Real inv_dz = 1._rt/dx[2]; -#else +#elif (AMREX_SPACEDIM == 2) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dz = 1._rt/dx[1]; +#else + const Real inv_dz = 1._rt/dx[0]; #endif +#if (AMREX_SPACEDIM >= 2) const Box& tbx = mfi.tilebox( E[lev][0]->ixType().toIntVect() ); +#endif #if (AMREX_SPACEDIM == 3) const Box& tby = mfi.tilebox( E[lev][1]->ixType().toIntVect() ); #endif const Box& tbz = mfi.tilebox( E[lev][2]->ixType().toIntVect() ); const auto& phi_arr = phi[lev]->array(mfi); +#if (AMREX_SPACEDIM >= 2) const auto& Ex_arr = (*E[lev][0])[mfi].array(); +#endif #if (AMREX_SPACEDIM == 3) const auto& Ey_arr = (*E[lev][1])[mfi].array(); #endif @@ -631,7 +646,7 @@ WarpX::computeE (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > +(beta_y*beta_z-1)*inv_dz*( phi_arr(i,j,k+1)-phi_arr(i,j,k) ); } ); -#else +#elif (AMREX_SPACEDIM == 2) amrex::ParallelFor( tbx, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Ex_arr(i,j,k) += @@ -646,6 +661,14 @@ WarpX::computeE (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > +(beta_y*beta_z-1)*inv_dz*( phi_arr(i,j+1,k)-phi_arr(i,j,k) ); } ); +#else + amrex::ParallelFor( tbz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Ez_arr(i,j,k) += + +(beta_y*beta_z-1)*inv_dz*( phi_arr(i+1,j,k)-phi_arr(i,j,k) ); + } + ); + ignore_unused(beta_x); #endif } } @@ -680,12 +703,15 @@ WarpX::computeB (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > #endif for ( MFIter mfi(*phi[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - const Real inv_dx = 1._rt/dx[0]; #if (AMREX_SPACEDIM == 3) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dy = 1._rt/dx[1]; const Real inv_dz = 1._rt/dx[2]; -#else +#elif (AMREX_SPACEDIM == 2) + const Real inv_dx = 1._rt/dx[0]; const Real inv_dz = 1._rt/dx[1]; +#else + const Real inv_dz = 1._rt/dx[0]; #endif const Box& tbx = mfi.tilebox( B[lev][0]->ixType().toIntVect() ); const Box& tby = mfi.tilebox( B[lev][1]->ixType().toIntVect() ); @@ -728,7 +754,7 @@ WarpX::computeB (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > + phi_arr(i+1,j+1,k)-phi_arr(i,j+1,k))); } ); -#else +#elif (AMREX_SPACEDIM == 2) amrex::ParallelFor( tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Bx_arr(i,j,k) += inv_c * ( @@ -746,6 +772,18 @@ WarpX::computeB (amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> > +beta_y*inv_dx*( phi_arr(i+1,j,k)-phi_arr(i,j,k) )); } ); +#else + amrex::ParallelFor( tbx, tby, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Bx_arr(i,j,k) += inv_c * ( + -beta_y*inv_dz*( phi_arr(i+1,j,k)-phi_arr(i,j,k) )); + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + By_arr(i,j,k) += inv_c * ( + +beta_x*inv_dz*(phi_arr(i+1,j,k)-phi_arr(i,j,k))); + } + ); + ignore_unused(beta_z,tbz,Bz_arr); #endif } } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp index 02b712679..46b32a8ff 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp @@ -165,9 +165,11 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( #else // Calculate relevant coefficients +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) amrex::Real const cdt_over_dx = PhysConst::c*dt*m_h_stencil_coefs_x[0]; amrex::Real const coef1_x = (1._rt - cdt_over_dx)/(1._rt + cdt_over_dx); amrex::Real const coef2_x = 2._rt*cdt_over_dx/(1._rt + cdt_over_dx) / PhysConst::c; +#endif #ifdef WARPX_DIM_3D amrex::Real const cdt_over_dy = PhysConst::c*dt*m_h_stencil_coefs_y[0]; amrex::Real const coef1_y = (1._rt - cdt_over_dy)/(1._rt + cdt_over_dy); @@ -177,8 +179,10 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( amrex::Real const coef1_z = (1._rt - cdt_over_dz)/(1._rt + cdt_over_dz); amrex::Real const coef2_z = 2._rt*cdt_over_dz/(1._rt + cdt_over_dz) / PhysConst::c; +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) bool const apply_lo_x = (field_boundary_lo[0] == FieldBoundaryType::Absorbing_SilverMueller); bool const apply_hi_x = (field_boundary_hi[0] == FieldBoundaryType::Absorbing_SilverMueller); +#endif #ifdef WARPX_DIM_3D bool const apply_lo_y = (field_boundary_lo[1] == FieldBoundaryType::Absorbing_SilverMueller); bool const apply_hi_y = (field_boundary_hi[1] == FieldBoundaryType::Absorbing_SilverMueller); @@ -201,10 +205,14 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // Extract field data for this grid/tile Array4<Real> const& Ex = Efield[0]->array(mfi); Array4<Real> const& Ey = Efield[1]->array(mfi); +#ifndef WARPX_DIM_1D_Z Array4<Real> const& Ez = Efield[2]->array(mfi); +#endif Array4<Real> const& Bx = Bfield[0]->array(mfi); Array4<Real> const& By = Bfield[1]->array(mfi); +#ifndef WARPX_DIM_1D_Z Array4<Real> const& Bz = Bfield[2]->array(mfi); +#endif // Extract the tileboxes for which to loop Box tbx = mfi.tilebox(Bfield[0]->ixType().toIntVect()); @@ -244,18 +252,27 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // At the -z boundary (innermost guard cell) if ( apply_lo_z && ( j==domain_box.smallEnd(1)-1 ) ) Bx(i,j,k) = coef1_z * Bx(i,j,k) + coef2_z * Ey(i,j+1,k); +#elif WARPX_DIM_1D_Z + // At the +z boundary (innermost guard cell) + if ( apply_hi_z && ( i==domain_box.bigEnd(0)+1 ) ) + Bx(i,j,k) = coef1_z * Bx(i,j,k) - coef2_z * Ey(i,j,k); + // At the -z boundary (innermost guard cell) + if ( apply_lo_z && ( i==domain_box.smallEnd(0)-1 ) ) + Bx(i,j,k) = coef1_z * Bx(i,j,k) + coef2_z * Ey(i+1,j,k); #endif }, // Apply Boundary condition to By [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) // At the +x boundary (innermost guard cell) if ( apply_hi_x && ( i==domain_box.bigEnd(0)+1 ) ) By(i,j,k) = coef1_x * By(i,j,k) - coef2_x * Ez(i,j,k); // At the -x boundary (innermost guard cell) if ( apply_lo_x && ( i==domain_box.smallEnd(0)-1 ) ) By(i,j,k) = coef1_x * By(i,j,k) + coef2_x * Ez(i+1,j,k); +#endif #ifdef WARPX_DIM_3D // At the +z boundary (innermost guard cell) if ( apply_hi_z && ( k==domain_box.bigEnd(2)+1 ) ) @@ -270,18 +287,27 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // At the -z boundary (innermost guard cell) if ( apply_lo_z && ( j==domain_box.smallEnd(1)-1 ) ) By(i,j,k) = coef1_z * By(i,j,k) - coef2_z * Ex(i,j+1,k); +#elif WARPX_DIM_1D_Z + // At the +z boundary (innermost guard cell) + if ( apply_hi_z && ( i==domain_box.bigEnd(0)+1 ) ) + By(i,j,k) = coef1_z * By(i,j,k) + coef2_z * Ex(i,j,k); + // At the -z boundary (innermost guard cell) + if ( apply_lo_z && ( i==domain_box.smallEnd(0)-1 ) ) + By(i,j,k) = coef1_z * By(i,j,k) - coef2_z * Ex(i+1,j,k); #endif }, // Apply Boundary condition to Bz [=] AMREX_GPU_DEVICE (int i, int j, int k){ +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) // At the +x boundary (innermost guard cell) if ( apply_hi_x && ( i==domain_box.bigEnd(0)+1 ) ) Bz(i,j,k) = coef1_x * Bz(i,j,k) + coef2_x * Ey(i,j,k); // At the -x boundary (innermost guard cell) if ( apply_lo_x && ( i==domain_box.smallEnd(0)-1 ) ) Bz(i,j,k) = coef1_x * Bz(i,j,k) - coef2_x * Ey(i+1,j,k); +#endif #ifdef WARPX_DIM_3D // At the +y boundary (innermost guard cell) if ( apply_hi_y && ( j==domain_box.bigEnd(1)+1 ) ) @@ -289,6 +315,8 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( // At the -y boundary (innermost guard cell) if ( apply_lo_y && ( j==domain_box.smallEnd(1)-1 ) ) Bz(i,j,k) = coef1_y * Bz(i,j,k) + coef2_y * Ex(i,j+1,k); +#elif WARPX_DIM_1D_Z + ignore_unused(i,j,k); #endif } ); diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp index ac19984cd..318081159 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp @@ -176,7 +176,6 @@ void FiniteDifferenceSolver::EvolveECartesian ( - T_Algo::DownwardDx(Bz, coefs_x, n_coefs_x, i, j, k) + T_Algo::DownwardDz(Bx, coefs_z, n_coefs_z, i, j, k) - PhysConst::mu0 * jy(i, j, k) ); - }, [=] AMREX_GPU_DEVICE (int i, int j, int k){ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H index 216723553..d7946c97e 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianCKCAlgorithm.H @@ -71,6 +71,13 @@ struct CartesianCKCAlgorithm { constexpr Real gammax=0._rt, gammay=0._rt, gammaz=0._rt; constexpr Real betaxy=0._rt, betazy=0._rt, betayx=0._rt, betayz=0._rt; constexpr Real alphay=0._rt; +#elif defined WARPX_DIM_1D_Z + Real const alphaz = inv_dz; + // Other coefficients are 0 in 1D Cartesian + // (and will actually not be used in the stencil) + constexpr Real gammax=0._rt, gammay=0._rt, gammaz=0._rt; + constexpr Real betaxy=0._rt, betazy=0._rt, betayx=0._rt, betayz=0._rt, betaxz=0._rt, betazx=0._rt; + constexpr Real alphay=0._rt, alphax=0._rt; #endif // Store the coefficients in array `stencil_coefs`, in prescribed order @@ -124,14 +131,16 @@ struct CartesianCKCAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) amrex::Real const alphax = coefs_x[1]; -#if defined WARPX_DIM_3D - amrex::Real const betaxy = coefs_x[2]; -#endif amrex::Real const betaxz = coefs_x[3]; +#endif #if defined WARPX_DIM_3D + amrex::Real const betaxy = coefs_x[2]; amrex::Real const gammax = coefs_x[4]; #endif + #if defined WARPX_DIM_3D return alphax * (F(i+1,j ,k ,ncomp) - F(i, j, k ,ncomp)) + betaxy * (F(i+1,j+1,k ,ncomp) - F(i ,j+1,k ,ncomp) @@ -146,6 +155,9 @@ struct CartesianCKCAlgorithm { return alphax * (F(i+1,j ,k ,ncomp) - F(i, j, k ,ncomp)) + betaxz * (F(i+1,j+1,k ,ncomp) - F(i ,j+1,k ,ncomp) + F(i+1,j-1,k ,ncomp) - F(i ,j-1,k ,ncomp)); +#elif (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 #endif } @@ -158,8 +170,14 @@ struct CartesianCKCAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else amrex::Real const inv_dx = coefs_x[0]; return inv_dx*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); +#endif } /** @@ -186,10 +204,10 @@ struct CartesianCKCAlgorithm { + F(i-1,j+1,k+1,ncomp) - F(i-1,j ,k+1,ncomp) + F(i+1,j+1,k-1,ncomp) - F(i+1,j ,k-1,ncomp) + F(i-1,j+1,k-1,ncomp) - F(i-1,j ,k-1,ncomp)); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #else amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); @@ -210,7 +228,7 @@ struct CartesianCKCAlgorithm { amrex::Real const inv_dy = coefs_y[0]; amrex::ignore_unused(n_coefs_y); return inv_dy*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); return 0._rt; // 2D Cartesian: derivative along y is 0 @@ -233,7 +251,9 @@ struct CartesianCKCAlgorithm { amrex::ignore_unused(n_coefs_z); Real const alphaz = coefs_z[1]; +#if (defined WARPX_DIM_3D || WARPX_DIM_XZ) Real const betazx = coefs_z[2]; +#endif #if defined WARPX_DIM_3D Real const betazy = coefs_z[3]; Real const gammaz = coefs_z[4]; @@ -252,6 +272,8 @@ struct CartesianCKCAlgorithm { return alphaz * (F(i ,j+1,k ,ncomp) - F(i ,j ,k ,ncomp)) + betazx * (F(i+1,j+1,k ,ncomp) - F(i+1,j ,k ,ncomp) + F(i-1,j+1,k ,ncomp) - F(i-1,j ,k ,ncomp)); +#elif (defined WARPX_DIM_1D_Z) + return alphaz * (F(i+1 ,j,k ,ncomp) - F(i ,j ,k ,ncomp)); #endif } @@ -269,6 +291,8 @@ struct CartesianCKCAlgorithm { return inv_dz*( F(i,j,k,ncomp) - F(i,j,k-1,ncomp) ); #elif (defined WARPX_DIM_XZ) return inv_dz*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return inv_dz*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); #endif } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H index 5c1d68687..01d83804b 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianNodalAlgorithm.H @@ -74,8 +74,13 @@ struct CartesianNodalAlgorithm { int const i, int const j, int const k, int const ncomp=0 ) { using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + ignore_unused(i, j, k, coefs_x, ncomp, F); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else Real const inv_dx = coefs_x[0]; return 0.5_rt*inv_dx*( F(i+1,j,k,ncomp) - F(i-1,j,k,ncomp) ); +#endif } /** @@ -88,8 +93,14 @@ struct CartesianNodalAlgorithm { amrex::Real const * const coefs_x, int const n_coefs_x, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + ignore_unused(i, j, k, coefs_x, n_coefs_x, ncomp, F); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else return UpwardDx( F, coefs_x, n_coefs_x, i, j, k ,ncomp); // For CartesianNodalAlgorithm, UpwardDx and DownwardDx are equivalent +#endif } /** @@ -106,9 +117,9 @@ struct CartesianNodalAlgorithm { #if defined WARPX_DIM_3D Real const inv_dy = coefs_y[0]; return 0.5_rt*inv_dy*( F(i,j+1,k,ncomp) - F(i,j-1,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_1D_Z) ignore_unused(i, j, k, coefs_y, ncomp, F); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #endif } @@ -142,6 +153,8 @@ struct CartesianNodalAlgorithm { return 0.5_rt*inv_dz*( F(i,j,k+1,ncomp) - F(i,j,k-1,ncomp) ); #elif (defined WARPX_DIM_XZ) return 0.5_rt*inv_dz*( F(i,j+1,k,ncomp) - F(i,j-1,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return 0.5_rt*inv_dz*( F(i+1,j,k,ncomp) - F(i-1,j,k,ncomp) ); #endif } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H index 4fe813754..f2c8553d5 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H @@ -71,8 +71,14 @@ struct CartesianYeeAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else amrex::Real const inv_dx = coefs_x[0]; return inv_dx*( F(i+1,j,k,ncomp) - F(i,j,k,ncomp) ); +#endif } /** @@ -84,8 +90,14 @@ struct CartesianYeeAlgorithm { amrex::Real const * const coefs_x, int const /*n_coefs_x*/, int const i, int const j, int const k, int const ncomp=0 ) { + using namespace amrex; +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(F, coefs_x, i, j, k, ncomp); + return 0._rt; // 1D Cartesian: derivative along x is 0 +#else amrex::Real const inv_dx = coefs_x[0]; return inv_dx*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); +#endif } /** @@ -101,10 +113,10 @@ struct CartesianYeeAlgorithm { Real const inv_dy = coefs_y[0]; amrex::ignore_unused(n_coefs_y); return inv_dy*( F(i,j+1,k,ncomp) - F(i,j,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #else amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); @@ -125,10 +137,10 @@ struct CartesianYeeAlgorithm { Real const inv_dy = coefs_y[0]; amrex::ignore_unused(n_coefs_y); return inv_dy*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); -#elif (defined WARPX_DIM_XZ) +#elif (defined WARPX_DIM_XZ || WARPX_DIM_1D_Z) amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); - return 0._rt; // 2D Cartesian: derivative along y is 0 + return 0._rt; // 1D and 2D Cartesian: derivative along y is 0 #else amrex::ignore_unused(F, coefs_y, n_coefs_y, i, j, k, ncomp); @@ -149,6 +161,8 @@ struct CartesianYeeAlgorithm { return inv_dz*( F(i,j,k+1,ncomp) - F(i,j,k,ncomp) ); #elif (defined WARPX_DIM_XZ) return inv_dz*( F(i,j+1,k,ncomp) - F(i,j,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return inv_dz*( F(i+1,j,k,ncomp) - F(i,j,k,ncomp) ); #endif } @@ -167,6 +181,8 @@ struct CartesianYeeAlgorithm { return inv_dz*( F(i,j,k,ncomp) - F(i,j,k-1,ncomp) ); #elif (defined WARPX_DIM_XZ) return inv_dz*( F(i,j,k,ncomp) - F(i,j-1,k,ncomp) ); +#elif (defined WARPX_DIM_1D_Z) + return inv_dz*( F(i,j,k,ncomp) - F(i-1,j,k,ncomp) ); #endif } diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 2ac27227d..6f1ebad78 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -32,6 +32,10 @@ void warpx_computedivb(int i, int j, int k, int dcomp, divB(i,j,0,dcomp) = (Bx(i+1,j ,0) - Bx(i,j,0))*dxinv + (Bz(i ,j+1,0) - Bz(i,j,0))*dzinv; amrex::ignore_unused(k, By, dyinv); +#elif defined WARPX_DIM_1D_Z + divB(i,0,0,dcomp) = (Bz(i+1,0 ,0) - Bz(i,0,0))*dzinv; + amrex::ignore_unused(j, Bx, dxinv); + amrex::ignore_unused(k, By, dyinv); #elif defined WARPX_DIM_RZ const amrex::Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); const amrex::Real rd = 1. - 0.5/(rmin*dxinv + i + 0.5); diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index 8a16e2a9b..601bd9ffc 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -80,9 +80,17 @@ void BilinearFilter::ComputeStencils(){ stencil_z.resize( 1u + npass_each_dir[1] ); compute_stencil(stencil_x, npass_each_dir[0]); compute_stencil(stencil_z, npass_each_dir[1]); +#elif (AMREX_SPACEDIM == 1) + // npass_each_dir = npass_z + stencil_z.resize( 1u + npass_each_dir[0] ); + compute_stencil(stencil_z, npass_each_dir[0]); #endif slen = stencil_length_each_dir.dim3(); #if (AMREX_SPACEDIM == 2) slen.z = 1; #endif +#if (AMREX_SPACEDIM == 1) + slen.y = 1; + slen.z = 1; +#endif } diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 654a1aa56..b72d6cec6 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -92,7 +92,9 @@ void Filter::DoFilter (const Box& tbx, Array4<Real > const& dst, int scomp, int dcomp, int ncomp) { +#if (AMREX_SPACEDIM >= 2) amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); +#endif #if (AMREX_SPACEDIM == 3) amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); #endif @@ -129,7 +131,7 @@ void Filter::DoFilter (const Box& tbx, dst(i,j,k,dcomp+n) = d; }); -#else +#elif (AMREX_SPACEDIM == 2) AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, { Real d = 0.0; @@ -155,6 +157,32 @@ void Filter::DoFilter (const Box& tbx, dst(i,j,k,dcomp+n) = d; }); +#elif (AMREX_SPACEDIM == 1) + AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, + { + Real d = 0.0; + + // Pad source array with zeros beyond ghost cells + // for out-of-bound accesses due to large-stencil operations + const auto src_zeropad = [src] (const int jj, const int kk, const int ll, const int nn) noexcept + { + return src.contains(jj,kk,ll) ? src(jj,kk,ll,nn) : 0.0_rt; + }; + + for (int iz=0; iz < slen_local.z; ++iz){ + for (int iy=0; iy < slen_local.y; ++iy){ + for (int ix=0; ix < slen_local.x; ++ix){ + Real sss = sz[iy]; + d += sss*( src_zeropad(i-ix,j,k,scomp+n) + +src_zeropad(i+ix,j,k,scomp+n)); + } + } + } + + dst(i,j,k,dcomp+n) = d; + }); +#else + amrex::Abort("Filter not implemented for the current geometry!"); #endif } @@ -247,7 +275,9 @@ void Filter::DoFilter (const Box& tbx, const auto lo = amrex::lbound(tbx); const auto hi = amrex::ubound(tbx); // tmp and dst are of type Array4 (Fortran ordering) +#if (AMREX_SPACEDIM >= 2) amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); +#endif #if (AMREX_SPACEDIM == 3) amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); #endif @@ -267,8 +297,10 @@ void Filter::DoFilter (const Box& tbx, for (int ix=0; ix < slen.x; ++ix){ #if (AMREX_SPACEDIM == 3) Real sss = sx[ix]*sy[iy]*sz[iz]; -#else +#elif (AMREX_SPACEDIM == 2) Real sss = sx[ix]*sz[iy]; +#else + Real sss = sz[ix]; #endif // 3 nested loop on 3D array for (int k = lo.z; k <= hi.z; ++k) { @@ -284,11 +316,16 @@ void Filter::DoFilter (const Box& tbx, +tmp(i+ix,j-iy,k+iz,scomp+n) +tmp(i-ix,j+iy,k+iz,scomp+n) +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else +#elif (AMREX_SPACEDIM == 2) dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) +tmp(i+ix,j-iy,k,scomp+n) +tmp(i-ix,j+iy,k,scomp+n) +tmp(i+ix,j+iy,k,scomp+n)); +#elif (AMREX_SPACEDIM == 1) + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j,k,scomp+n) + +tmp(i+ix,j,k,scomp+n)); +#else + amrex::Abort("Filter not implemented for the current geometry!"); #endif } } diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 224bfb946..fa88eb854 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -23,6 +23,7 @@ using namespace amrex; NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set, amrex::Real cdtodz, bool nodal_gather){ +#if (AMREX_SPACEDIM >= 2) // Store parameters into class data members m_coeff_set = coeff_set; m_cdtodz = cdtodz; @@ -35,10 +36,15 @@ NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set, amrex::Real cdto stencil_length_each_dir = {1,5}; slen = {1,5,1}; #endif +#else + amrex::ignore_unused(coeff_set, cdtodz, nodal_gather); + amrex::Abort("NCIGodfreyFilter not implemented in 1D!"); +#endif } void NCIGodfreyFilter::ComputeStencils(){ +#if (AMREX_SPACEDIM >= 2) using namespace warpx::nci_godfrey; // Sanity checks: filter length shoulz be 5 in z @@ -128,4 +134,7 @@ void NCIGodfreyFilter::ComputeStencils(){ Gpu::copyAsync(Gpu::hostToDevice,h_stencil_z.begin(),h_stencil_z.end(),stencil_z.begin()); Gpu::synchronize(); +#else + amrex::Abort("NCIGodfreyFilter not implemented in 1D!"); +#endif } diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index af6a7a208..5710cab9a 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -81,10 +81,16 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) // NOTE: When periodic boundaries are used, default injection range is set to mother grid dimensions. const amrex::Geometry& geom = WarpX::GetInstance().Geom(0); if( geom.isPeriodic(0) ) { +# ifndef WARPX_DIM_1D_Z xmin = geom.ProbLo(0); xmax = geom.ProbHi(0); +# else + zmin = geom.ProbLo(0); + zmax = geom.ProbHi(0); +# endif } +# ifndef WARPX_DIM_1D_Z if( geom.isPeriodic(1) ) { # ifndef WARPX_DIM_3D zmin = geom.ProbLo(1); @@ -94,6 +100,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) ymax = geom.ProbHi(1); # endif } +# endif # ifdef WARPX_DIM_3D if( geom.isPeriodic(2) ) { @@ -246,16 +253,18 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) flux_normal_axis = 0; } #else +# ifndef WARPX_DIM_1D_Z if (flux_normal_axis_string == "x" || flux_normal_axis_string == "X") { flux_normal_axis = 0; } +# endif #endif #ifdef WARPX_DIM_3D - else if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") { + if (flux_normal_axis_string == "y" || flux_normal_axis_string == "Y") { flux_normal_axis = 1; } #endif - else if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") { + if (flux_normal_axis_string == "z" || flux_normal_axis_string == "Z") { flux_normal_axis = AMREX_SPACEDIM-1; } #ifdef WARPX_DIM_3D @@ -263,8 +272,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) #else # ifdef WARPX_DIM_RZ std::string flux_normal_axis_help = "'r' or 'z'."; -# else +# elif WARPX_DIM_XZ std::string flux_normal_axis_help = "'x' or 'z'."; +# else + std::string flux_normal_axis_help = "'z'."; # endif #endif AMREX_ALWAYS_ASSERT_WITH_MESSAGE(flux_normal_axis >= 0, @@ -282,7 +293,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) } else if (injection_style == "nuniformpercell") { // Note that for RZ, three numbers are expected, r, theta, and z. // For 2D, only two are expected. The third is overwritten with 1. -#if WARPX_DIM_XZ + // For 1D, only one is expected. The second and third are overwritten with 1. +#if defined(WARPX_DIM_1D_Z) + constexpr int num_required_ppc_each_dim = 1; +#elif defined(WARPX_DIM_XZ) constexpr int num_required_ppc_each_dim = 2; #else constexpr int num_required_ppc_each_dim = 3; @@ -292,6 +306,10 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) #if WARPX_DIM_XZ num_particles_per_cell_each_dim.push_back(1); #endif +#if WARPX_DIM_1D_Z + num_particles_per_cell_each_dim.push_back(1); // overwrite 2nd number with 1 + num_particles_per_cell_each_dim.push_back(1); // overwrite 3rd number with 1 +#endif #if WARPX_DIM_RZ if (WarpX::n_rz_azimuthal_modes > 1) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 41fda5567..f92fe80aa 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -86,11 +86,13 @@ WarpX::PostProcessBaseGrids (BoxArray& ba0) const klo += domlo[2]; khi += domlo[2]; #endif +#if (AMREX_SPACEDIM >= 2) for (int j = 0; j < numprocs[1]; ++j) { int jlo = (j < extra[1]) ? j*(sz[1]+1) : (j*sz[1]+extra[1]); int jhi = (j < extra[1]) ? jlo+(sz[1]+1)-1 : jlo+sz[1]-1; jlo += domlo[1]; jhi += domlo[1]; +#endif for (int i = 0; i < numprocs[0]; ++i) { int ilo = (i < extra[0]) ? i*(sz[0]+1) : (i*sz[0]+extra[0]); int ihi = (i < extra[0]) ? ilo+(sz[0]+1)-1 : ilo+sz[0]-1; @@ -358,6 +360,7 @@ WarpX::computeMaxStepBoostAccelerator(const amrex::Geometry& a_geom){ void WarpX::InitNCICorrector () { +#if !(defined WARPX_DIM_1D_Z) if (WarpX::use_fdtd_nci_corr) { for (int lev = 0; lev <= max_level; ++lev) @@ -367,8 +370,10 @@ WarpX::InitNCICorrector () amrex::Real dz, cdtodz; if (AMREX_SPACEDIM == 3){ dz = dx[2]; - }else{ + }else if(AMREX_SPACEDIM == 2){ dz = dx[1]; + }else{ + dz = dx[0]; } cdtodz = PhysConst::c * dt[lev] / dz; @@ -385,6 +390,7 @@ WarpX::InitNCICorrector () nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); } } +#endif } void @@ -663,13 +669,20 @@ WarpX::InitializeExternalFieldsOnGridUsingParser ( #endif // Shift required in the x-, y-, or z- position // depending on the index type of the multifab +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt; amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0._rt; amrex::Real fac_z = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; #else + amrex::Real fac_x = (1._rt - x_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1._rt - x_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1._rt - x_nodal_flag[2]) * dx_lev[2] * 0.5_rt; @@ -682,13 +695,20 @@ WarpX::InitializeExternalFieldsOnGridUsingParser ( #ifdef AMREX_USE_EB if(geom_data_y(i, j, k)<=0) return; #endif +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt; amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0._rt; amrex::Real fac_z = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; #elif (AMREX_SPACEDIM==3) + amrex::Real fac_x = (1._rt - y_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1._rt - y_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1._rt - y_nodal_flag[2]) * dx_lev[2] * 0.5_rt; @@ -701,13 +721,20 @@ WarpX::InitializeExternalFieldsOnGridUsingParser ( #ifdef AMREX_USE_EB if(geom_data_z(i, j, k)<=0) return; #endif +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0._rt; + amrex::Real y = 0._rt; + amrex::Real fac_z = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; + amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt; amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0._rt; amrex::Real fac_z = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real z = j*dx_lev[1] + real_box.lo(1) + fac_z; #elif (AMREX_SPACEDIM==3) + amrex::Real fac_x = (1._rt - z_nodal_flag[0]) * dx_lev[0] * 0.5_rt; + amrex::Real x = i*dx_lev[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1._rt - z_nodal_flag[1]) * dx_lev[1] * 0.5_rt; amrex::Real y = j*dx_lev[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1._rt - z_nodal_flag[2]) * dx_lev[2] * 0.5_rt; diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 5df0635e7..e974e5764 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -58,7 +58,9 @@ include $(AMREX_HOME)/Tools/GNUMake/Make.defs ifeq ($(DIM),3) DEFINES += -DWARPX_DIM_3D -else +endif + +ifeq ($(DIM),2) ifeq ($(USE_RZ),TRUE) DEFINES += -DWARPX_DIM_RZ else @@ -66,6 +68,10 @@ else endif endif +ifeq ($(DIM),1) + DEFINES += -DWARPX_DIM_1D_Z +endif + -include Make.package include $(WARPX_HOME)/Source/Make.package include $(WARPX_HOME)/Source/BoundaryConditions/Make.package diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp index 44848c915..b72e13cdc 100644 --- a/Source/Parallelization/GuardCellManager.cpp +++ b/Source/Parallelization/GuardCellManager.cpp @@ -115,6 +115,9 @@ guardCellManager::Init ( #elif (AMREX_SPACEDIM == 2) ng_alloc_EB = IntVect(ngx,ngz); ng_alloc_J = IntVect(ngJx,ngJz); +#elif (AMREX_SPACEDIM == 1) + ng_alloc_EB = IntVect(ngz); + ng_alloc_J = IntVect(ngJz); #endif // TODO Adding one cell for rho should not be necessary, given that the number of guard cells @@ -195,6 +198,8 @@ guardCellManager::Init ( IntVect ngFFT = IntVect(ngFFt_x, ngFFt_y, ngFFt_z); #elif (AMREX_SPACEDIM == 2) IntVect ngFFT = IntVect(ngFFt_x, ngFFt_z); +#elif (AMREX_SPACEDIM == 1) + IntVect ngFFT = IntVect(ngFFt_z); #endif // All boxes should have the same number of guard cells, to avoid temporary parallel copies: diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 2d1e1f04a..ecaddb570 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -979,7 +979,9 @@ WarpX::ApplyFilterandSumBoundaryJ (int lev, PatchType patch_type) IntVect ng_depos_J = get_ng_depos_J(); if (WarpX::do_current_centering) { -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + ng_depos_J[0] += WarpX::current_centering_noz / 2; +#elif (AMREX_SPACEDIM == 2) ng_depos_J[0] += WarpX::current_centering_nox / 2; ng_depos_J[1] += WarpX::current_centering_noz / 2; #elif (AMREX_SPACEDIM == 3) diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index 04226e971..5995243f5 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -20,17 +20,17 @@ void warpx_interp (int j, int k, int l, { using namespace amrex; - // NOTE Indices (j,k,l) in the following refer to (x,z,-) in 2D and (x,y,z) in 3D + // NOTE Indices (j,k,l) in the following refer to (z,-,-) in 1D, (x,z,-) in 2D, and (x,y,z) in 3D // Refinement ratio const int rj = rr[0]; - const int rk = rr[1]; - const int rl = (AMREX_SPACEDIM == 2) ? 1 : rr[2]; + const int rk = (AMREX_SPACEDIM == 1) ? 1 : rr[1]; + const int rl = (AMREX_SPACEDIM <= 2) ? 1 : rr[2]; // Staggering (0: cell-centered; 1: nodal) const int sj = arr_stag[0]; - const int sk = arr_stag[1]; - const int sl = (AMREX_SPACEDIM == 2) ? 0 : arr_stag[2]; + const int sk = (AMREX_SPACEDIM == 1) ? 0 : arr_stag[1]; + const int sl = (AMREX_SPACEDIM <= 2) ? 0 : arr_stag[2]; // Number of points used for interpolation from coarse grid to fine grid const int nj = (sj == 0) ? 1 : 2; @@ -87,7 +87,21 @@ void warpx_interp_nd_bfield_x (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real bg = owx * Bxg(jg ,0,0) + + wx * Bxg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real bc = owx * Bxc(jg ,0,0) + + wx * Bxc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5_rt*(Bxf_zeropad(j-1,0,0) + Bxf_zeropad(j,0,0)); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real bg = owx * owy * Bxg(jg ,kg ,0) @@ -163,7 +177,21 @@ void warpx_interp_nd_bfield_y (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real bg = owx * Byg(jg ,0,0) + + wx * Byg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real bc = owx * Byc(jg ,0,0) + + wx * Byc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5_rt*(Byf_zeropad(j-1,0,0) + Byf_zeropad(j,0,0)); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real bg = owx * owy * Byg(jg ,kg ,0) @@ -241,7 +269,21 @@ void warpx_interp_nd_bfield_z (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real bg = owx * Bzg(jg ,0,0) + + wx * Bzg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real bc = owx * Bzc(jg ,0,0) + + wx * Bzc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real bf = Bzf_zeropad(j,0,0); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real bg = owx * owy * Bzg(jg ,kg ,0) @@ -318,7 +360,21 @@ void warpx_interp_nd_efield_x (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real eg = owx * Exg(jg ,0,0) + + wx * Exg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real ec = owx * Exc(jg ,0,0) + + wx * Exc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real ef = Exf_zeropad(j,0,0); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real eg = owx * owy * Exg(jg ,kg ,0) @@ -394,7 +450,23 @@ void warpx_interp_nd_efield_y (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) + + +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real eg = owx * Eyg(jg ,0,0) + + wx * Eyg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real ec = owx * Eyc(jg ,0,0) + + wx * Eyc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real ef = Eyf_zeropad(j,0,0); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real eg = owx * owy * Eyg(jg ,kg ,0) @@ -469,7 +541,21 @@ void warpx_interp_nd_efield_z (int j, int k, int l, Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; Real owy = 1.0_rt-wy; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + + // interp from coarse nodal to fine nodal + Real eg = owx * Ezg(jg ,0,0) + + wx * Ezg(jg+1,0,0); + + // interp from coarse staggered to fine nodal + Real ec = owx * Ezc(jg ,0,0) + + wx * Ezc(jg+1,0,0); + + // interp from fine staggered to fine nodal + Real ef = 0.5_rt*(Ezf_zeropad(j-1,0,0) + Ezf_zeropad(j,0,0)); + amrex::ignore_unused(owy); + +#elif (AMREX_SPACEDIM == 2) // interp from coarse nodal to fine nodal Real eg = owx * owy * Ezg(jg ,kg ,0) @@ -568,6 +654,9 @@ void warpx_interp (const int j, }; // Avoid compiler warnings +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(nox, noy, stencil_coeffs_x, stencil_coeffs_y); +#endif #if (AMREX_SPACEDIM == 2) amrex::ignore_unused(noy, stencil_coeffs_y); #endif @@ -586,32 +675,43 @@ void warpx_interp (const int j, // Staggering (s = 0 if cell-centered, s = 1 if nodal) const int sj = (dst_nodal) ? src_stag[0] : dst_stag[0]; +#if (AMREX_SPACEDIM >= 2) const int sk = (dst_nodal) ? src_stag[1] : dst_stag[1]; +#endif #if (AMREX_SPACEDIM == 3) const int sl = (dst_nodal) ? src_stag[2] : dst_stag[2]; #endif // Interpolate along j,k,l only if source MultiFab is staggered along j,k,l const bool interp_j = (sj == 0); +#if (AMREX_SPACEDIM >= 2) const bool interp_k = (sk == 0); +#endif #if (AMREX_SPACEDIM == 3) const bool interp_l = (sl == 0); #endif +#if (AMREX_SPACEDIM == 1) + const int noj = noz; +#elif (AMREX_SPACEDIM == 2) const int noj = nox; -#if (AMREX_SPACEDIM == 2) const int nok = noz; #elif (AMREX_SPACEDIM == 3) + const int noj = nox; const int nok = noy; const int nol = noz; #endif // Additional normalization factor const amrex::Real wj = (interp_j) ? 0.5_rt : 1.0_rt; +#if (AMREX_SPACEDIM == 1) + constexpr amrex::Real wk = 1.0_rt; + constexpr amrex::Real wl = 1.0_rt; +#elif (AMREX_SPACEDIM == 2) const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt; -#if (AMREX_SPACEDIM == 2) constexpr amrex::Real wl = 1.0_rt; #elif (AMREX_SPACEDIM == 3) + const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt; const amrex::Real wl = (interp_l) ? 0.5_rt : 1.0_rt; #endif @@ -620,11 +720,17 @@ void warpx_interp (const int j, const int jmax = (interp_j) ? j + noj/2 + shift - 1 : j; // Min and max for interpolation loop along k +#if (AMREX_SPACEDIM == 1) + // k = 0 always + const int kmin = k; + const int kmax = k; +#else const int kmin = (interp_k) ? k - nok/2 + shift : k; const int kmax = (interp_k) ? k + nok/2 + shift - 1 : k; +#endif // Min and max for interpolation loop along l -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM <= 2) // l = 0 always const int lmin = l; const int lmax = l; @@ -696,6 +802,11 @@ void warpx_interp (const int j, else // PSATD (finite-order interpolation) { +#if (AMREX_SPACEDIM == 1) + amrex::Abort("PSATD not implemented in 1D"); +#endif + +#if (AMREX_SPACEDIM >= 2) // 1D not implemented for PSATD amrex::Real cj = 1.0_rt; amrex::Real ck = 1.0_rt; amrex::Real cl = 1.0_rt; @@ -725,6 +836,7 @@ void warpx_interp (const int j, } } } +#endif //1D } dst_arr(j,k,l) = wj * wk * wl * res; diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index fa3a0effc..98da517d6 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -257,7 +257,9 @@ public: const amrex::Real dt = WarpX::GetInstance().getdt(lev); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); -#if defined WARPX_DIM_XZ +#if defined WARPX_DIM_1D_Z + auto dV = geom.CellSize(0); +#elif defined WARPX_DIM_XZ auto dV = geom.CellSize(0) * geom.CellSize(1); #elif defined WARPX_DIM_RZ amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box @@ -406,7 +408,9 @@ public: const amrex::Real dt = WarpX::GetInstance().getdt(lev); amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev); -#if defined WARPX_DIM_XZ +#if defined WARPX_DIM_1D_Z + auto dV = geom.CellSize(0); +#elif defined WARPX_DIM_XZ auto dV = geom.CellSize(0) * geom.CellSize(1); #elif defined WARPX_DIM_RZ amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index d651b1b6c..84db47ba5 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -58,16 +58,22 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) const bool do_ionization = ion_lev; - const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dzi = 1.0_rt/dx[2]; +#if (AMREX_SPACEDIM == 1) + const amrex::Real invvol = dzi; +#endif #if (AMREX_SPACEDIM == 2) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real invvol = dxi*dzi; #elif (defined WARPX_DIM_3D) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dyi = 1.0_rt/dx[1]; const amrex::Real invvol = dxi*dyi*dzi; #endif +#if (AMREX_SPACEDIM >= 2) const amrex::Real xmin = xyzmin[0]; +#endif #if (defined WARPX_DIM_3D) const amrex::Real ymin = xyzmin[1]; #endif @@ -105,6 +111,8 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, GetPosition(ip, xp, yp, zp); // --- Compute shape factors + Compute_shape_factor< depos_order > const compute_shape_factor; +#if (AMREX_SPACEDIM >= 2) // x direction // Get particle position in grid coordinates #if (defined WARPX_DIM_RZ) @@ -128,13 +136,12 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, // i: leftmost grid point that the particle touches amrex::Real sx[depos_order + 1] = {0._rt}; int i = 0; - Compute_shape_factor< depos_order > const compute_shape_factor; if (rho_type[0] == NODE) { i = compute_shape_factor(sx, x); } else if (rho_type[0] == CELL) { i = compute_shape_factor(sx, x - 0.5_rt); } - +#endif //AMREX_SPACEDIM >= 2 #if (defined WARPX_DIM_3D) // y direction const amrex::Real y = (yp - ymin)*dyi; @@ -157,6 +164,13 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, } // Deposit charge into rho_arr +#if (defined WARPX_DIM_1D_Z) + for (int iz=0; iz<=depos_order; iz++){ + amrex::Gpu::Atomic::AddNoRet( + &rho_arr(lo.x+k+iz, 0, 0, 0), + sz[iz]*wq); + } +#endif #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 6c330e388..640ba4514 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -79,16 +79,22 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) const bool do_ionization = ion_lev; - const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dzi = 1.0_rt/dx[2]; +#if (AMREX_SPACEDIM == 1) + const amrex::Real invvol = dzi; +#endif #if (AMREX_SPACEDIM == 2) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real invvol = dxi*dzi; #elif (defined WARPX_DIM_3D) + const amrex::Real dxi = 1.0_rt/dx[0]; const amrex::Real dyi = 1.0_rt/dx[1]; const amrex::Real invvol = dxi*dyi*dzi; #endif +#if (AMREX_SPACEDIM >= 2) const amrex::Real xmin = xyzmin[0]; +#endif #if (defined WARPX_DIM_3D) const amrex::Real ymin = xyzmin[1]; #endif @@ -164,6 +170,8 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, const amrex::Real wqz = wq*invvol*vz; // --- Compute shape factors + Compute_shape_factor< depos_order > const compute_shape_factor; +#if (AMREX_SPACEDIM >= 2) // x direction // Get particle position after 1/2 push back in position #if (defined WARPX_DIM_RZ) @@ -181,7 +189,6 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, double sx_cell[depos_order + 1] = {0.}; int j_node = 0; int j_cell = 0; - Compute_shape_factor< depos_order > const compute_shape_factor; if (jx_type[0] == NODE || jy_type[0] == NODE || jz_type[0] == NODE) { j_node = compute_shape_factor(sx_node, xmid); } @@ -202,6 +209,7 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, int const j_jx = ((jx_type[0] == NODE) ? j_node : j_cell); int const j_jy = ((jy_type[0] == NODE) ? j_node : j_cell); int const j_jz = ((jz_type[0] == NODE) ? j_node : j_cell); +#endif //AMREX_SPACEDIM >= 2 #if (defined WARPX_DIM_3D) // y direction @@ -258,6 +266,19 @@ void doDepositionShapeN(const GetParticlePosition& GetPosition, int const l_jz = ((jz_type[zdir] == NODE) ? l_node : l_cell); // Deposit current into jx_arr, jy_arr and jz_arr +#if (defined WARPX_DIM_1D_Z) + for (int iz=0; iz<=depos_order; iz++){ + amrex::Gpu::Atomic::AddNoRet( + &jx_arr(lo.x+l_jx+iz, 0, 0, 0), + sz_jx[iz]*wqx); + amrex::Gpu::Atomic::AddNoRet( + &jy_arr(lo.x+l_jy+iz, 0, 0, 0), + sz_jy[iz]*wqy); + amrex::Gpu::Atomic::AddNoRet( + &jz_arr(lo.x+l_jz+iz, 0, 0, 0), + sz_jz[iz]*wqz); + } +#endif #if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ @@ -367,11 +388,15 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) bool const do_ionization = ion_lev; +#if !(defined WARPX_DIM_1D_Z) Real const dxi = 1.0_rt / dx[0]; -#if !(defined WARPX_DIM_RZ) +#endif +#if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z) Real const dtsdx0 = dt*dxi; #endif +#if !(defined WARPX_DIM_1D_Z) Real const xmin = xyzmin[0]; +#endif #if (defined WARPX_DIM_3D) Real const dyi = 1.0_rt / dx[1]; Real const dtsdy0 = dt*dyi; @@ -389,6 +414,9 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Real const invdtdx = 1.0_rt / (dt*dx[2]); Real const invdtdz = 1.0_rt / (dt*dx[0]); Real const invvol = 1.0_rt / (dx[0]*dx[2]); +#elif (defined WARPX_DIM_1D_Z) + Real const invdtdz = 1.0_rt / (dt*dx[0]); + Real const invvol = 1.0_rt / (dx[2]); #endif #if (defined WARPX_DIM_RZ) @@ -396,8 +424,10 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, #endif Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c ); +#if !(defined WARPX_DIM_1D_Z) Real constexpr one_third = 1.0_rt / 3.0_rt; Real constexpr one_sixth = 1.0_rt / 6.0_rt; +#endif // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr #if defined(WARPX_USE_GPUCLOCK) @@ -429,7 +459,9 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, ParticleReal xp, yp, zp; GetPosition(ip, xp, yp, zp); +#if !(defined WARPX_DIM_1D_Z) Real const wqx = wq*invdtdx; +#endif #if (defined WARPX_DIM_3D) Real const wqy = wq*invdtdy; #endif @@ -476,9 +508,11 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, double const x_old = (rp_old - xmin)*dxi; #else // Keep these double to avoid bug in single precision +#if !(defined WARPX_DIM_1D_Z) double const x_new = (xp - xmin)*dxi; double const x_old = x_new - dtsdx0*uxp[ip]*gaminv; #endif +#endif #if (defined WARPX_DIM_3D) // Keep these double to avoid bug in single precision double const y_new = (yp - ymin)*dyi; @@ -492,6 +526,9 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; #elif (defined WARPX_DIM_XZ) Real const vy = uyp[ip]*gaminv; +#elif (defined WARPX_DIM_1D_Z) + Real const vx = uxp[ip]*gaminv; + Real const vy = uyp[ip]*gaminv; #endif // Shape factor arrays @@ -499,8 +536,10 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, // to possibly hold the factor for the old particle // which can be at a different grid location. // Keep these double to avoid bug in single precision +#if !(defined WARPX_DIM_1D_Z) double sx_new[depos_order + 3] = {0.}; double sx_old[depos_order + 3] = {0.}; +#endif #if (defined WARPX_DIM_3D) // Keep these double to avoid bug in single precision double sy_new[depos_order + 3] = {0.}; @@ -516,8 +555,10 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Compute_shape_factor< depos_order > compute_shape_factor; Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; +#if !(defined WARPX_DIM_1D_Z) const int i_new = compute_shape_factor(sx_new+1, x_new); const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new); +#endif #if (defined WARPX_DIM_3D) const int j_new = compute_shape_factor(sy_new+1, y_new); const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new); @@ -526,9 +567,11 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new); // computes min/max positions of current contributions +#if !(defined WARPX_DIM_1D_Z) int dil = 1, diu = 1; if (i_old < i_new) dil = 0; if (i_old > i_new) diu = 0; +#endif #if (defined WARPX_DIM_3D) int djl = 1, dju = 1; if (j_old < j_new) djl = 0; @@ -636,6 +679,23 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, #endif } } +#elif (defined WARPX_DIM_1D_Z) + + for (int k=dkl; k<=depos_order+2-dku; k++) { + amrex::Real sdxi = 0._rt; + sdxi += wq*vx*invvol*0.5_rt*(sz_old[k] + sz_new[k]); + amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+k_new-1+k, 0, 0, 0), sdxi); + } + for (int k=dkl; k<=depos_order+2-dku; k++) { + amrex::Real sdyj = 0._rt; + sdyj += wq*vy*invvol*0.5_rt*(sz_old[k] + sz_new[k]); + amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+k_new-1+k, 0, 0, 0), sdyj); + } + for (int k=dkl; k<=depos_order+1-dku; k++) { + amrex::Real sdzk = 0._rt; + sdzk += wqz*(sz_old[k] - sz_new[k]); + amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k_new-1+k, 0, 0, 0), sdzk); + } #endif } ); @@ -700,11 +760,18 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, amrex::Abort("Vay deposition not implemented in RZ geometry"); #endif +#if (defined WARPX_DIM_1D_Z) + amrex::ignore_unused(GetPosition, + wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab, + np_to_depose, dt, dx, xyzmin, lo, q, n_rz_azimuthal_modes); + amrex::Abort("Vay deposition not implemented in cartesian 1D geometry"); +#endif + #if !defined(AMREX_USE_GPU) amrex::ignore_unused(cost, load_balance_costs_update_algo); #endif -#if !(defined WARPX_DIM_RZ) +#if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z) amrex::ignore_unused(n_rz_azimuthal_modes); // If ion_lev is a null pointer, then do_ionization=0, else do_ionization=1 @@ -1001,6 +1068,6 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, amrex::The_Managed_Arena()->free(cost_real); } # endif -#endif // #if !(defined WARPX_DIM_RZ) +#endif // #if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z) } #endif // CURRENTDEPOSITION_H_ diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 562f5b242..23a27f330 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -67,17 +67,25 @@ void doGatherShapeN (const amrex::ParticleReal xp, amrex::ignore_unused(yp); #endif +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(xp,yp); +#endif + #ifndef WARPX_DIM_RZ amrex::ignore_unused(n_rz_azimuthal_modes); #endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 3) const amrex::Real dxi = 1.0_rt/dx[0]; +#endif const amrex::Real dzi = 1.0_rt/dx[2]; #if (AMREX_SPACEDIM == 3) const amrex::Real dyi = 1.0_rt/dx[1]; #endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 3) const amrex::Real xmin = xyzmin[0]; +#endif #if (AMREX_SPACEDIM == 3) const amrex::Real ymin = xyzmin[1]; #endif @@ -88,6 +96,11 @@ void doGatherShapeN (const amrex::ParticleReal xp, constexpr int CELL = amrex::IndexType::CELL; // --- Compute shape factors + + Compute_shape_factor< depos_order > const compute_shape_factor; + Compute_shape_factor<depos_order - galerkin_interpolation > const compute_shape_factor_galerkin; + +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 3) // x direction // Get particle position #ifdef WARPX_DIM_RZ @@ -110,8 +123,6 @@ void doGatherShapeN (const amrex::ParticleReal xp, int j_cell = 0; int j_node_v = 0; int j_cell_v = 0; - Compute_shape_factor< depos_order > const compute_shape_factor; - Compute_shape_factor<depos_order - galerkin_interpolation > const compute_shape_factor_galerkin; if ((ey_type[0] == NODE) || (ez_type[0] == NODE) || (bx_type[0] == NODE)) { j_node = compute_shape_factor(sx_node, x); } @@ -136,6 +147,7 @@ void doGatherShapeN (const amrex::ParticleReal xp, int const j_bx = ((bx_type[0] == NODE) ? j_node : j_cell ); int const j_by = ((by_type[0] == NODE) ? j_node_v : j_cell_v); int const j_bz = ((bz_type[0] == NODE) ? j_node_v : j_cell_v); +#endif #if (AMREX_SPACEDIM == 3) // y direction @@ -214,7 +226,33 @@ void doGatherShapeN (const amrex::ParticleReal xp, // AMREX_SPACEDIM nested loops because the deposition // order can differ for each component of each field // when galerkin_interpolation is set to 1 -#if (AMREX_SPACEDIM == 2) + +#if (AMREX_SPACEDIM == 1) + // Gather field on particle Eyp from field on grid ey_arr + // Gather field on particle Exp from field on grid ex_arr + // Gather field on particle Bzp from field on grid bz_arr + for (int iz=0; iz<=depos_order; iz++){ + Eyp += sz_ey[iz]* + ey_arr(lo.x+l_ey+iz, 0, 0, 0); + Exp += sz_ex[iz]* + ex_arr(lo.x+l_ex+iz, 0, 0, 0); + Bzp += sz_bz[iz]* + bz_arr(lo.x+l_bz+iz, 0, 0, 0); + } + + // Gather field on particle Byp from field on grid by_arr + // Gather field on particle Ezp from field on grid ez_arr + // Gather field on particle Bxp from field on grid bx_arr + for (int iz=0; iz<=depos_order-galerkin_interpolation; iz++){ + Ezp += sz_ez[iz]* + ez_arr(lo.x+l_ez+iz, 0, 0, 0); + Bxp += sz_bx[iz]* + bx_arr(lo.x+l_bx+iz, 0, 0, 0); + Byp += sz_by[iz]* + by_arr(lo.x+l_by+iz, 0, 0, 0); + } + +#elif (AMREX_SPACEDIM == 2) // Gather field on particle Eyp from field on grid ey_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index 559489f34..e7c47a6f4 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -134,6 +134,12 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_nvec[1] == amrex::Real(0), "Laser propagation direction must be 0 along y in 2D"); #endif +#ifdef WARPX_DIM_1D_Z + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_nvec[0] == amrex::Real(0), + "Laser propagation direction must be 0 along x in 1D"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_nvec[1] == amrex::Real(0), + "Laser propagation direction must be 0 along y in 2D"); +#endif // Plane normal Real s = 1.0_rt / std::sqrt(m_nvec[0]*m_nvec[0] + m_nvec[1]*m_nvec[1] + m_nvec[2]*m_nvec[2]); @@ -168,9 +174,12 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) m_u_X = m_p_X; m_u_Y = m_p_Y; -#else +#elif (defined WARPX_DIM_XZ) m_u_X = CrossProduct({0., 1., 0.}, m_nvec); m_u_Y = {0., 1., 0.}; +#elif (defined WARPX_DIM_1D_Z) + m_u_X = {1., 0., 0.}; + m_u_Y = {0., 1., 0.}; #endif m_laser_injection_box= Geom(0).ProbDomain(); @@ -192,7 +201,10 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, // Sanity checks int dir = WarpX::moving_window_dir; std::vector<Real> windir(3, 0.0); -#if (AMREX_SPACEDIM==2) +#if (AMREX_SPACEDIM==1) + windir[2] = 1.0; + amrex::ignore_unused(dir); +#elif (AMREX_SPACEDIM==2) windir[2*dir] = 1.0; #else windir[dir] = 1.0; @@ -243,8 +255,10 @@ LaserParticleContainer::ContinuousInjection (const RealBox& injection_box) // Convert updated_position to Real* to use RealBox::contains(). #if (AMREX_SPACEDIM == 3) const Real* p_pos = m_updated_position.dataPtr(); -#else +#elif (AMREX_SPACEDIM == 2) const Real p_pos[2] = {m_updated_position[0], m_updated_position[2]}; +#else + const Real p_pos[1] = {m_updated_position[2]}; #endif if ( injection_box.contains(p_pos) ){ // Update laser_injection_box with current value @@ -278,6 +292,13 @@ LaserParticleContainer::UpdateContinuousInjectionPosition (Real dt) // which has 3 components, for both 2D and 3D simulations. m_updated_position[2*dir] -= WarpX::beta_boost * WarpX::boost_direction[2*dir] * PhysConst::c * dt; +#elif ( AMREX_SPACEDIM == 1 ) + // In 1D, dir=0 corresponds to z + // This needs to be converted in order to index `boost_direction` + // which has 3 components, for 1D, 2D, and 3D simulations. + m_updated_position[2] -= WarpX::beta_boost * + WarpX::boost_direction[2] * PhysConst::c * dt; + amrex::ignore_unused(dir); #endif } } @@ -317,12 +338,13 @@ LaserParticleContainer::InitData (int lev) m_position = m_updated_position; } +#if (AMREX_SPACEDIM >= 2) auto Transform = [&](int const i, int const j) -> Vector<Real>{ #if (AMREX_SPACEDIM == 3) return { m_position[0] + (S_X*(Real(i)+0.5_rt))*m_u_X[0] + (S_Y*(Real(j)+0.5_rt))*m_u_Y[0], m_position[1] + (S_X*(Real(i)+0.5_rt))*m_u_X[1] + (S_Y*(Real(j)+0.5_rt))*m_u_Y[1], m_position[2] + (S_X*(Real(i)+0.5_rt))*m_u_X[2] + (S_Y*(Real(j)+0.5_rt))*m_u_Y[2] }; -#else +#elif (AMREX_SPACEDIM == 2) amrex::ignore_unused(j); # if (defined WARPX_DIM_RZ) return { m_position[0] + (S_X*(Real(i)+0.5_rt)), @@ -335,18 +357,21 @@ LaserParticleContainer::InitData (int lev) # endif #endif }; +#endif // Given the "lab" frame coordinates, return the real coordinates in the laser plane coordinates auto InverseTransform = [&](const Vector<Real>& pos) -> Vector<Real>{ #if (AMREX_SPACEDIM == 3) return {m_u_X[0]*(pos[0]-m_position[0])+m_u_X[1]*(pos[1]-m_position[1])+m_u_X[2]*(pos[2]-m_position[2]), m_u_Y[0]*(pos[0]-m_position[0])+m_u_Y[1]*(pos[1]-m_position[1])+m_u_Y[2]*(pos[2]-m_position[2])}; -#else +#elif (AMREX_SPACEDIM == 2) # if (defined WARPX_DIM_RZ) return {pos[0]-m_position[0], 0.0_rt}; # else return {m_u_X[0]*(pos[0]-m_position[0])+m_u_X[2]*(pos[2]-m_position[2]), 0.0_rt}; # endif +#else + return {m_u_X[2]*(pos[2]-m_position[2]), 0.0_rt}; #endif }; @@ -374,11 +399,14 @@ LaserParticleContainer::InitData (int lev) compute_min_max(prob_hi[0], prob_lo[1], prob_hi[2]); compute_min_max(prob_lo[0], prob_hi[1], prob_hi[2]); compute_min_max(prob_hi[0], prob_hi[1], prob_hi[2]); -#else +#elif (AMREX_SPACEDIM == 2) compute_min_max(prob_lo[0], 0.0, prob_lo[1]); compute_min_max(prob_hi[0], 0.0, prob_lo[1]); compute_min_max(prob_lo[0], 0.0, prob_hi[1]); compute_min_max(prob_hi[0], 0.0, prob_hi[1]); +#else + compute_min_max(0.0, 0.0, prob_lo[0]); + compute_min_max(0.0, 0.0, prob_hi[0]); #endif } @@ -404,8 +432,10 @@ LaserParticleContainer::InitData (int lev) } } } -#else +#elif (AMREX_SPACEDIM == 2) BoxArray plane_ba { Box {IntVect(plane_lo[0],0), IntVect(plane_hi[0],0)} }; +#else + BoxArray plane_ba { Box {IntVect(0), IntVect(0)} }; #endif amrex::Vector<amrex::Real> particle_x, particle_y, particle_z, particle_w; @@ -419,11 +449,17 @@ LaserParticleContainer::InitData (int lev) const Box& bx = plane_ba[i]; for (IntVect cell = bx.smallEnd(); cell <= bx.bigEnd(); bx.next(cell)) { +#if (AMREX_SPACEDIM >= 2) const Vector<Real>& pos = Transform(cell[0], cell[1]); +#else + const Vector<Real>& pos = { 0.0_rt, 0.0_rt, m_position[2] }; +#endif #if (AMREX_SPACEDIM == 3) const Real* x = pos.dataPtr(); -#else +#elif (AMREX_SPACEDIM == 2) const Real x[2] = {pos[0], pos[2]}; +#else + const Real x[1] = {pos[2]}; #endif if (m_laser_injection_box.contains(x)) { @@ -583,6 +619,7 @@ LaserParticleContainer::Evolve (int lev, } } + if (rho && ! skip_deposition) { int* AMREX_RESTRICT ion_lev = nullptr; DepositCharge(pti, wp, ion_lev, rho, 1, 0, @@ -636,7 +673,7 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const Sy = std::min(std::min(dx[0]/(std::abs(m_u_Y[0])+eps), dx[1]/(std::abs(m_u_Y[1])+eps)), dx[2]/(std::abs(m_u_Y[2])+eps)); -#else +#elif (AMREX_SPACEDIM == 2) # if (defined WARPX_DIM_RZ) Sx = dx[0]; # else @@ -644,6 +681,10 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const dx[2]/(std::abs(m_u_X[2])+eps)); # endif Sy = 1.0; +#else + Sx = 1.0; + Sy = 1.0; + amrex::ignore_unused(eps); #endif } @@ -694,6 +735,7 @@ LaserParticleContainer::calculate_laser_plane_coordinates (const WarpXParIter& p { const auto GetPosition = GetParticlePosition(pti); +#if (AMREX_SPACEDIM >= 2) Real tmp_u_X_0 = m_u_X[0]; Real tmp_u_X_2 = m_u_X[2]; Real tmp_position_0 = m_position[0]; @@ -705,6 +747,7 @@ LaserParticleContainer::calculate_laser_plane_coordinates (const WarpXParIter& p Real tmp_u_Y_2 = m_u_Y[2]; Real tmp_position_1 = m_position[1]; #endif +#endif amrex::ParallelFor( np, @@ -725,6 +768,9 @@ LaserParticleContainer::calculate_laser_plane_coordinates (const WarpXParIter& p tmp_u_X_0 * (x - tmp_position_0) + tmp_u_X_2 * (z - tmp_position_2); pplane_Yp[i] = 0.; +#elif (AMREX_SPACEDIM == 1) + pplane_Xp[i] = 0.; + pplane_Yp[i] = 0.; #endif } ); @@ -792,7 +838,9 @@ LaserParticleContainer::update_laser_particle (WarpXParIter& pti, // Push the the particle positions ParticleReal x, y, z; GetPosition(i, x, y, z); +#if !(defined WARPX_DIM_1D_Z) x += vx * dt; +#endif #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) y += vy * dt; #endif diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index b877708f9..70254aadd 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -357,6 +357,10 @@ MultiParticleContainer::ReadParameters () AMREX_ALWAYS_ASSERT_WITH_MESSAGE(WarpX::use_fdtd_nci_corr==0, "ERROR: use_fdtd_nci_corr is not supported in RZ"); #endif +#ifdef WARPX_DIM_1D_Z + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(WarpX::use_fdtd_nci_corr==0, + "ERROR: use_fdtd_nci_corr is not supported in 1D"); +#endif ParmParse pp_lasers("lasers"); pp_lasers.queryarr("names", lasers_names); @@ -1303,6 +1307,9 @@ MultiParticleContainer::doQEDSchwinger () #ifdef WARPX_DIM_RZ amrex::Abort("Schwinger process not implemented in rz geometry"); #endif +#ifdef WARPX_DIM_1D_Z + amrex::Abort("Schwinger process not implemented in 1D geometry"); +#endif // Get cell volume. In 2D the transverse size is // chosen by the user in the input file. diff --git a/Source/Particles/ParticleBoundaries_K.H b/Source/Particles/ParticleBoundaries_K.H index f33b29c17..48be4e65a 100644 --- a/Source/Particles/ParticleBoundaries_K.H +++ b/Source/Particles/ParticleBoundaries_K.H @@ -82,7 +82,10 @@ namespace ApplyParticleBoundaries { */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void - apply_boundaries (amrex::ParticleReal& x, amrex::Real xmin, amrex::Real xmax, + apply_boundaries ( +#ifndef WARPX_DIM_1D_Z + amrex::ParticleReal& x, amrex::Real xmin, amrex::Real xmax, +#endif #ifdef WARPX_DIM_3D amrex::ParticleReal& y, amrex::Real ymin, amrex::Real ymax, #endif @@ -96,10 +99,12 @@ namespace ApplyParticleBoundaries { bool change_sign_uy = false; bool change_sign_uz = false; +#ifndef WARPX_DIM_1D_Z apply_boundary(x, xmin, xmax, change_sign_ux, particle_lost, boundaries.xmin_bc, boundaries.xmax_bc, boundaries.reflection_model_xlo(-ux), boundaries.reflection_model_xhi(ux), engine); +#endif #ifdef WARPX_DIM_3D apply_boundary(y, ymin, ymax, change_sign_uy, particle_lost, boundaries.ymin_bc, boundaries.ymax_bc, diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp index 64efcefeb..d8856c297 100644 --- a/Source/Particles/ParticleBoundaryBuffer.cpp +++ b/Source/Particles/ParticleBoundaryBuffer.cpp @@ -69,12 +69,17 @@ ParticleBoundaryBuffer::ParticleBoundaryBuffer () for (int ispecies = 0; ispecies < numSpecies(); ++ispecies) { amrex::ParmParse pp_species(getSpeciesNames()[ispecies]); +#if AMREX_SPACEDIM == 1 + pp_species.query("save_particles_at_zlo", m_do_boundary_buffer[0][ispecies]); + pp_species.query("save_particles_at_zhi", m_do_boundary_buffer[1][ispecies]); +#elif AMREX_SPACEDIM == 2 pp_species.query("save_particles_at_xlo", m_do_boundary_buffer[0][ispecies]); pp_species.query("save_particles_at_xhi", m_do_boundary_buffer[1][ispecies]); -#if AMREX_SPACEDIM == 2 pp_species.query("save_particles_at_zlo", m_do_boundary_buffer[2][ispecies]); pp_species.query("save_particles_at_zhi", m_do_boundary_buffer[3][ispecies]); #else + pp_species.query("save_particles_at_xlo", m_do_boundary_buffer[0][ispecies]); + pp_species.query("save_particles_at_xhi", m_do_boundary_buffer[1][ispecies]); pp_species.query("save_particles_at_ylo", m_do_boundary_buffer[2][ispecies]); pp_species.query("save_particles_at_yhi", m_do_boundary_buffer[3][ispecies]); pp_species.query("save_particles_at_zlo", m_do_boundary_buffer[4][ispecies]); diff --git a/Source/Particles/ParticleCreation/SmartCreate.H b/Source/Particles/ParticleCreation/SmartCreate.H index 926168985..957efb2c4 100644 --- a/Source/Particles/ParticleCreation/SmartCreate.H +++ b/Source/Particles/ParticleCreation/SmartCreate.H @@ -46,13 +46,17 @@ struct SmartCreate const int cpu = 0, const int id = 0) const noexcept { - prt.m_aos[i_prt].pos(0) = x; #if (AMREX_SPACEDIM == 3) + prt.m_aos[i_prt].pos(0) = x; prt.m_aos[i_prt].pos(1) = y; prt.m_aos[i_prt].pos(2) = z; #elif (AMREX_SPACEDIM == 2) + prt.m_aos[i_prt].pos(0) = x; prt.m_aos[i_prt].pos(1) = z; amrex::ignore_unused(y); +#else + prt.m_aos[i_prt].pos(0) = z; + amrex::ignore_unused(x,y); #endif prt.m_aos[i_prt].cpu() = cpu; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c331b7a61..e34c4c28b 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -145,7 +145,7 @@ namespace pos.x = lo_corner[0] + (iv[0]+r.x)*dx[0]; pos.y = lo_corner[1] + (iv[1]+r.y)*dx[1]; pos.z = lo_corner[2] + (iv[2]+r.z)*dx[2]; -#else +#elif (AMREX_SPACEDIM == 2) pos.x = lo_corner[0] + (iv[0]+r.x)*dx[0]; pos.y = 0.0_rt; #if defined WARPX_DIM_XZ @@ -154,6 +154,10 @@ namespace // Note that for RZ, r.y will be theta pos.z = lo_corner[1] + (iv[1]+r.z)*dx[1]; #endif +#else + pos.x = 0.0_rt; + pos.y = 0.0_rt; + pos.z = lo_corner[0] + (iv[0]+r.z)*dx[0]; #endif return pos; } @@ -185,7 +189,9 @@ namespace ) noexcept { p.pos(0) = 0._rt; +#if (AMREX_SPACEDIM >= 2) p.pos(1) = 0._rt; +#endif #if (AMREX_SPACEDIM == 3) p.pos(2) = 0._rt; #endif @@ -283,7 +289,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If old particle positions should be saved add the needed components pp_species_name.query("save_previous_position", m_save_previous_position); if (m_save_previous_position) { +#if (AMREX_SPACEDIM >= 2) AddRealComp("prev_x"); +#endif #if (AMREX_SPACEDIM == 3) AddRealComp("prev_y"); #endif @@ -435,6 +443,11 @@ PhysicalParticleContainer::AddGaussianBeam ( const Real x = amrex::RandomNormal(x_m, x_rms); constexpr Real y = 0._prt; const Real z = amrex::RandomNormal(z_m, z_rms); +#elif (defined WARPX_DIM_1D_Z) + const Real weight = q_tot/(npart*charge*x_rms*y_rms); + constexpr Real x = 0._prt; + constexpr Real y = 0._prt; + const Real z = amrex::RandomNormal(z_m, z_rms); #endif if (plasma_injector->insideBounds(x, y, z) && std::abs( x - x_m ) < x_cut * x_rms && @@ -506,18 +519,20 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot, openPMD::ParticleSpecies ps = it.particles.begin()->second; auto const npart = ps["position"]["x"].getExtent()[0]; +#if !(defined WARPX_DIM_1D_Z) std::shared_ptr<ParticleReal> ptr_x = ps["position"]["x"].loadChunk<ParticleReal>(); double const position_unit_x = ps["position"]["x"].unitSI(); +#endif std::shared_ptr<ParticleReal> ptr_z = ps["position"]["z"].loadChunk<ParticleReal>(); double const position_unit_z = ps["position"]["z"].unitSI(); std::shared_ptr<ParticleReal> ptr_ux = ps["momentum"]["x"].loadChunk<ParticleReal>(); double const momentum_unit_x = ps["momentum"]["x"].unitSI(); std::shared_ptr<ParticleReal> ptr_uz = ps["momentum"]["z"].loadChunk<ParticleReal>(); double const momentum_unit_z = ps["momentum"]["z"].unitSI(); -# ifndef WARPX_DIM_XZ +# if !(defined(WARPX_DIM_XZ) || defined(WARPX_DIM_1D_Z)) std::shared_ptr<ParticleReal> ptr_y = ps["position"]["y"].loadChunk<ParticleReal>(); double const position_unit_y = ps["position"]["y"].unitSI(); -# endif +#endif std::shared_ptr<ParticleReal> ptr_uy = nullptr; double momentum_unit_y = 1.0; if (ps["momentum"].contains("y")) { @@ -547,13 +562,17 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot, } for (auto i = decltype(npart){0}; i<npart; ++i){ +#if !(defined WARPX_DIM_1D_Z) ParticleReal const x = ptr_x.get()[i]*position_unit_x; +#else + ParticleReal const x = 0.0_prt; +#endif ParticleReal const z = ptr_z.get()[i]*position_unit_z+z_shift; -# if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) +#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) ParticleReal const y = ptr_y.get()[i]*position_unit_y; -# else +#else ParticleReal const y = 0.0_prt; -# endif +#endif if (plasma_injector->insideBounds(x, y, z)) { ParticleReal const ux = ptr_ux.get()[i]*momentum_unit_x/PhysConst::m_e; ParticleReal const uz = ptr_uz.get()[i]*momentum_unit_z/PhysConst::m_e; @@ -711,6 +730,8 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) scale_fac = dx[0]*dx[1]*dx[2]/num_ppc; #elif AMREX_SPACEDIM==2 scale_fac = dx[0]*dx[1]/num_ppc; +#elif AMREX_SPACEDIM==1 + scale_fac = dx[0]/num_ppc; #endif } @@ -839,9 +860,12 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pcounts[index] = num_ppc; } } -#if (AMREX_SPACEDIM != 3) +#if (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); #endif +#if (AMREX_SPACEDIM == 1) + amrex::ignore_unused(j,k); +#endif }); // Max number of new particles. All of them are created, @@ -954,7 +978,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) ); continue; } -#else +#elif (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi @@ -965,6 +989,17 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) ); continue; } +#else + amrex::ignore_unused(j,k); + if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { + ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi +#ifdef WARPX_QED + ,loc_has_quantum_sync, p_optical_depth_QSR + ,loc_has_breit_wheeler, p_optical_depth_BW +#endif + ); + continue; + } #endif // Save the x and y values to use in the insideBounds checks. @@ -1110,6 +1145,8 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) #endif p.pos(0) = xb; p.pos(1) = pos.z; +#else + p.pos(0) = pos.z; #endif } }); @@ -1287,8 +1324,10 @@ PhysicalParticleContainer::AddPlasmaFlux (int lev, amrex::Real dt) pcounts[index] = num_ppc_int; } } -#if (AMREX_SPACEDIM != 3) +#if (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); +#elif (AMREX_SPACEDIM == 1) + amrex::ignore_unused(j,k); #endif }); @@ -1402,12 +1441,18 @@ PhysicalParticleContainer::AddPlasmaFlux (int lev, amrex::Real dt) p.id() = -1; continue; } -#else +#elif (AMREX_SPACEDIM == 2) amrex::ignore_unused(k); if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { p.id() = -1; continue; } +#else + amrex::ignore_unused(j,k); + if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { + p.id() = -1; + continue; + } #endif // Save the x and y values to use in the insideBounds checks. @@ -1491,6 +1536,8 @@ PhysicalParticleContainer::AddPlasmaFlux (int lev, amrex::Real dt) #endif p.pos(0) = xb; p.pos(1) = pos.z; +#else + p.pos(0) = pos.z; #endif } }); @@ -1768,7 +1815,9 @@ PhysicalParticleContainer::applyNCIFilter ( const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::noz)}); +#elif (AMREX_SPACEDIM == 2) const Box& tbox = amrex::grow(box,{static_cast<int>(WarpX::nox), static_cast<int>(WarpX::noz)}); #else @@ -1836,7 +1885,13 @@ PhysicalParticleContainer::SplitParticles (int lev) long np_split; if(split_type==0) { - np_split = (AMREX_SPACEDIM == 3) ? 8 : 4; + if (AMREX_SPACEDIM == 3){ + np_split = 8; + } else if (AMREX_SPACEDIM == 2){ + np_split = 4; + } else { + np_split = 2; + } } else { np_split = 2*AMREX_SPACEDIM; } @@ -1876,7 +1931,20 @@ PhysicalParticleContainer::SplitParticles (int lev) // If particle is tagged, split it and put the // split particles in local arrays psplit_x etc. np_split_to_add += np_split; -#if (AMREX_SPACEDIM==2) +#if (AMREX_SPACEDIM==1) + // Split particle in two along z axis + // 2 particles in 1d, split_type doesn't matter? Discuss with Remi + for (int ishift = -1; ishift < 2; ishift +=2 ){ + // Add one particle with offset in z + psplit_x.push_back( xp ); + psplit_y.push_back( yp ); + psplit_z.push_back( zp + ishift*split_offset[2] ); + psplit_ux.push_back( uxp[i] ); + psplit_uy.push_back( uyp[i] ); + psplit_uz.push_back( uzp[i] ); + psplit_w.push_back( wp[i]/np_split ); + } +#elif (AMREX_SPACEDIM==2) if (split_type==0){ // Split particle in two along each diagonals // 4 particles in 2d @@ -2133,7 +2201,9 @@ PhysicalParticleContainer::GetParticleSlice ( WARPX_PROFILE("PhysicalParticleContainer::GetParticleSlice()"); // Assume that the boost in the positive z direction. -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + AMREX_ALWAYS_ASSERT(direction == 0); +#elif (AMREX_SPACEDIM == 2) AMREX_ALWAYS_ASSERT(direction == 1); #else AMREX_ALWAYS_ASSERT(direction == 2); @@ -2428,7 +2498,11 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, ParticleReal* y_old = nullptr; ParticleReal* z_old = nullptr; if (save_previous_position) { +#if (AMREX_SPACEDIM >= 2) x_old = pti.GetAttribs(particle_comps["prev_x"]).dataPtr(); +#else + amrex::ignore_unused(x_old); +#endif #if (AMREX_SPACEDIM == 3) y_old = pti.GetAttribs(particle_comps["prev_y"]).dataPtr(); #else @@ -2465,7 +2539,9 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, getPosition(ip, xp, yp, zp); if (save_previous_position) { +#if (AMREX_SPACEDIM >= 2) x_old[ip] = xp; +#endif #if (AMREX_SPACEDIM == 3) y_old[ip] = yp; #endif diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index b0f1257f5..4df1ef3c1 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -35,10 +35,14 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p, x = p.pos(0); y = p.pos(1); z = p.pos(2); -#else +#elif WARPX_DIM_XZ x = p.pos(0); y = amrex::ParticleReal(0.0); z = p.pos(1); +#else + x = amrex::ParticleReal(0.0); + y = amrex::ParticleReal(0.0); + z = p.pos(0); #endif } @@ -55,6 +59,9 @@ struct GetParticlePosition const RType* m_theta = nullptr; #elif (AMREX_SPACEDIM == 2) static constexpr RType m_y_default = RType(0.0); +#elif (AMREX_SPACEDIM == 1) + static constexpr RType m_x_default = RType(0.0); + static constexpr RType m_y_default = RType(0.0); #endif GetParticlePosition () = default; @@ -93,10 +100,14 @@ struct GetParticlePosition x = p.pos(0); y = p.pos(1); z = p.pos(2); -#else +#elif WARPX_DIM_XZ x = p.pos(0); y = m_y_default; z = p.pos(1); +#else + x = m_x_default; + y = m_y_default; + z = p.pos(0); #endif } @@ -117,10 +128,14 @@ struct GetParticlePosition x = p.pos(0); y = p.pos(1); z = p.pos(2); -#else +#elif WARPX_DIM_XZ x = p.pos(0); y = m_y_default; z = p.pos(1); +#else + x = m_x_default; + y = m_y_default; + z = p.pos(0); #endif } }; @@ -158,6 +173,9 @@ struct SetParticlePosition #if defined(WARPX_DIM_XZ) amrex::ignore_unused(y); #endif +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(x,y); +#endif #ifdef WARPX_DIM_RZ m_theta[i] = std::atan2(y, x); m_structs[i].pos(0) = std::sqrt(x*x + y*y); @@ -166,9 +184,11 @@ struct SetParticlePosition m_structs[i].pos(0) = x; m_structs[i].pos(1) = y; m_structs[i].pos(2) = z; -#else +#elif WARPX_DIM_XZ m_structs[i].pos(0) = x; m_structs[i].pos(1) = z; +#else + m_structs[i].pos(0) = z; #endif } @@ -182,6 +202,9 @@ struct SetParticlePosition #if defined(WARPX_DIM_XZ) amrex::ignore_unused(y); #endif +#if defined(WARPX_DIM_1D_Z) + amrex::ignore_unused(x,y); +#endif #ifdef WARPX_DIM_RZ m_structs[i].pos(0) = x; m_theta[i] = y; @@ -190,9 +213,11 @@ struct SetParticlePosition m_structs[i].pos(0) = x; m_structs[i].pos(1) = y; m_structs[i].pos(2) = z; -#else +#elif WARPX_DIM_XZ m_structs[i].pos(0) = x; m_structs[i].pos(1) = z; +#else + m_structs[i].pos(0) = z; #endif } }; diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 8ebb0b49a..39b12cd67 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -29,7 +29,11 @@ void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::Parti // Compute inverse Lorentz factor const amrex::Real inv_gamma = 1._rt/std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)*inv_c2); // Update positions over one time step +#if (AMREX_SPACEDIM >= 2) x += ux * inv_gamma * dt; +#else + amrex::ignore_unused(x); +#endif #if (AMREX_SPACEDIM == 3) || (defined WARPX_DIM_RZ) // RZ pushes particles in 3D y += uy * inv_gamma * dt; #else diff --git a/Source/Particles/Pusher/UpdatePositionPhoton.H b/Source/Particles/Pusher/UpdatePositionPhoton.H index f52eb1c75..5e958c2c1 100644 --- a/Source/Particles/Pusher/UpdatePositionPhoton.H +++ b/Source/Particles/Pusher/UpdatePositionPhoton.H @@ -32,7 +32,11 @@ void UpdatePositionPhoton( const amrex::Real c_over_umod = (u_norm == 0._rt) ? 0._rt: PhysConst::c/u_norm; // Update positions over one time step +#if (AMREX_SPACEDIM >= 2) x += ux * c_over_umod * dt; +#else + amrex::ignore_unused(x); +#endif #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) // RZ pushes particles in 3D y += uy * c_over_umod * dt; #else diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 554d8fee4..f7cb27b28 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -109,12 +109,16 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) // The boundary conditions are read in in ReadBCParams but a child class // can allow these value to be overwritten if different boundary // conditions are desired for a specific species +#ifndef WARPX_DIM_1D_Z m_boundary_conditions.SetBoundsX(WarpX::particle_boundary_lo[0], WarpX::particle_boundary_hi[0]); +#endif #ifdef WARPX_DIM_3D m_boundary_conditions.SetBoundsY(WarpX::particle_boundary_lo[1], WarpX::particle_boundary_hi[1]); m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[2], WarpX::particle_boundary_hi[2]); -#else +#elif WARPX_DIM_XZ || WARPX_DIM_RZ m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[1], WarpX::particle_boundary_hi[1]); +#else + m_boundary_conditions.SetBoundsZ(WarpX::particle_boundary_lo[0], WarpX::particle_boundary_hi[0]); #endif m_boundary_conditions.BuildReflectionModelParsers(); } @@ -212,6 +216,9 @@ WarpXParticleContainer::AddNParticles (int /*lev*/, p.pos(0) = x[i]; #endif p.pos(1) = z[i]; +#else //AMREX_SPACEDIM == 1 + amrex::ignore_unused(x,y); + p.pos(0) = z[i]; #endif pinned_tile.push_back(p); @@ -325,7 +332,9 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, // deposit part of its current in a neighboring box. However, this should catch particles // traveling many cells away, for example with algorithms that allow for large time steps. -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::noz/2)); +#elif (AMREX_SPACEDIM == 2) const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::nox/2), static_cast<int>(WarpX::noz/2)); #elif (AMREX_SPACEDIM == 3) @@ -608,7 +617,9 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, // are not trivial, this check might be too strict and we might need to relax it, as currently // done for the current deposition. -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::noz/2+1)); +#elif (AMREX_SPACEDIM == 2) const amrex::IntVect shape_extent = amrex::IntVect(static_cast<int>(WarpX::nox/2+1), static_cast<int>(WarpX::noz/2+1)); #elif (AMREX_SPACEDIM == 3) @@ -1132,8 +1143,10 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ { auto GetPosition = GetParticlePosition(pti); auto SetPosition = SetParticlePosition(pti); +#ifndef WARPX_DIM_1D_Z const Real xmin = Geom(lev).ProbLo(0); const Real xmax = Geom(lev).ProbHi(0); +#endif #ifdef WARPX_DIM_3D const Real ymin = Geom(lev).ProbLo(1); const Real ymax = Geom(lev).ProbHi(1); @@ -1163,7 +1176,10 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ // Note that for RZ, (x, y, z) is actually (r, theta, z). bool particle_lost = false; - ApplyParticleBoundaries::apply_boundaries(x, xmin, xmax, + ApplyParticleBoundaries::apply_boundaries( +#ifndef WARPX_DIM_1D_Z + x, xmin, xmax, +#endif #ifdef WARPX_DIM_3D y, ymin, ymax, #endif diff --git a/Source/Utils/CoarsenIO.cpp b/Source/Utils/CoarsenIO.cpp index 2cb97b0e1..e2ea419ad 100644 --- a/Source/Utils/CoarsenIO.cpp +++ b/Source/Utils/CoarsenIO.cpp @@ -41,24 +41,36 @@ CoarsenIO::Loop ( MultiFab& mf_dst, GpuArray<int,3> cr; // coarsening ratio sf[0] = stag_src[0]; +#if (AMREX_SPACEDIM == 1) + sf[1] = 0; +#else sf[1] = stag_src[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sf[2] = 0; #elif (AMREX_SPACEDIM == 3) sf[2] = stag_src[2]; #endif sc[0] = stag_dst[0]; +#if (AMREX_SPACEDIM == 1) + sc[1] = 0; +#else sc[1] = stag_dst[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sc[2] = 0; #elif (AMREX_SPACEDIM == 3) sc[2] = stag_dst[2]; #endif cr[0] = crse_ratio[0]; +#if (AMREX_SPACEDIM == 1) + cr[1] = 1; +#else cr[1] = crse_ratio[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) cr[2] = 1; #elif (AMREX_SPACEDIM == 3) cr[2] = crse_ratio[2]; diff --git a/Source/Utils/CoarsenMR.cpp b/Source/Utils/CoarsenMR.cpp index 80562d542..b701af088 100644 --- a/Source/Utils/CoarsenMR.cpp +++ b/Source/Utils/CoarsenMR.cpp @@ -30,24 +30,36 @@ CoarsenMR::Loop ( MultiFab& mf_dst, GpuArray<int,3> cr; // coarsening ratio sf[0] = stag_src[0]; +#if (AMREX_SPACEDIM == 1) + sf[1] = 0; +#else sf[1] = stag_src[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sf[2] = 0; #elif (AMREX_SPACEDIM == 3) sf[2] = stag_src[2]; #endif sc[0] = stag_dst[0]; +#if (AMREX_SPACEDIM == 1) + sc[1] = 0; +#else sc[1] = stag_dst[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) sc[2] = 0; #elif (AMREX_SPACEDIM == 3) sc[2] = stag_dst[2]; #endif cr[0] = crse_ratio[0]; +#if (AMREX_SPACEDIM == 1) + cr[1] = 1; +#else cr[1] = crse_ratio[1]; -#if (AMREX_SPACEDIM == 2) +#endif +#if (AMREX_SPACEDIM == 2 || AMREX_SPACEDIM == 1) cr[2] = 1; #elif (AMREX_SPACEDIM == 3) cr[2] = crse_ratio[2]; diff --git a/Source/Utils/Interpolate_K.H b/Source/Utils/Interpolate_K.H index 9c8e0a2b4..f452f1acf 100644 --- a/Source/Utils/Interpolate_K.H +++ b/Source/Utils/Interpolate_K.H @@ -18,11 +18,16 @@ void interp (int j, int k, int l, Real const wx = static_cast<Real>(type[0]) * static_cast<amrex::Real>(j-jg*r_ratio) * rr; Real const owx = 1.0_rt-wx; +#if (AMREX_SPACEDIM >= 2) int const kg = amrex::coarsen(k,r_ratio); Real const wy = static_cast<Real>(type[1]) * static_cast<amrex::Real>(k-kg*r_ratio) * rr; Real const owy = 1.0_rt-wy; +#endif -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + fine(j,k,l) = owx * crse(jg ,0,0) + + wx * crse(jg+1,0,0); +#elif (AMREX_SPACEDIM == 2) fine(j,k,l) = owx * owy * crse(jg ,kg ,0) + owx * wy * crse(jg ,kg+1,0) + wx * owy * crse(jg+1,kg ,0) diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 600ad180e..3bb48d1a5 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -63,6 +63,12 @@ WarpX::UpdatePlasmaInjectionPosition (amrex::Real a_dt) // This needs to be converted in order to index `boost_direction` // which has 3 components, for both 2D and 3D simulations. WarpX::boost_direction[2*dir] * PhysConst::c * a_dt; +#elif ( AMREX_SPACEDIM == 1 ) + // In 1D, dir=0 corresponds to z + // This needs to be converted in order to index `boost_direction` + // which has 3 components, for 1D, 2D, and 3D simulations. + WarpX::boost_direction[2] * PhysConst::c * a_dt; + amrex::ignore_unused(dir); #endif } } @@ -365,13 +371,20 @@ WarpX::shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, int num_shift, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { // Compute x,y,z co-ordinates based on index type of mf +#if (AMREX_SPACEDIM==1) + amrex::Real x = 0.0; + amrex::Real y = 0.0; + amrex::Real fac_z = (1.0 - mf_type[0]) * dx[0]*0.5; + amrex::Real z = i*dx[0] + real_box.lo(0) + fac_z; +#elif (AMREX_SPACEDIM==2) amrex::Real fac_x = (1.0 - mf_type[0]) * dx[0]*0.5; amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; -#if (AMREX_SPACEDIM==2) amrex::Real y = 0.0; amrex::Real fac_z = (1.0 - mf_type[1]) * dx[1]*0.5; amrex::Real z = j*dx[1] + real_box.lo(1) + fac_z; #else + amrex::Real fac_x = (1.0 - mf_type[0]) * dx[0]*0.5; + amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; amrex::Real fac_y = (1.0 - mf_type[1]) * dx[1]*0.5; amrex::Real y = j*dx[1] + real_box.lo(1) + fac_y; amrex::Real fac_z = (1.0 - mf_type[2]) * dx[2]*0.5; @@ -417,6 +430,11 @@ WarpX::ShiftGalileanBoundary () m_v_galilean[0]*time_shift, std::numeric_limits<amrex::Real>::quiet_NaN(), m_v_galilean[2]*time_shift }; +#elif (AMREX_SPACEDIM == 1) + m_galilean_shift = { + std::numeric_limits<Real>::quiet_NaN(), + std::numeric_limits<Real>::quiet_NaN(), + m_v_galilean[2]*time_shift }; #endif #if (AMREX_SPACEDIM == 3) @@ -430,8 +448,13 @@ WarpX::ShiftGalileanBoundary () new_hi[0] = current_hi[0] + m_galilean_shift[0]; new_lo[1] = current_lo[1] + m_galilean_shift[2]; new_hi[1] = current_hi[1] + m_galilean_shift[2]; - } - #endif + } +#elif (AMREX_SPACEDIM == 1) + { + new_lo[0] = current_lo[0] + m_galilean_shift[2]; + new_hi[0] = current_hi[0] + m_galilean_shift[2]; + } +#endif time_of_last_gal_shift = cur_time; ResetProbDomain(amrex::RealBox(new_lo, new_hi)); diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 06f0a7da5..d77a35b3a 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -170,8 +170,10 @@ void ConvertLabParamsToBoost() #if (AMREX_SPACEDIM == 3) Vector<int> dim_map {0, 1, 2}; -#else +#elif (AMREX_SPACEDIM == 2) Vector<int> dim_map {0, 2}; +#else + Vector<int> dim_map {2}; #endif for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) @@ -219,6 +221,8 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) amrex::Array<amrex::Real,3> galilean_shift = { 0., 0., 0., }; #elif (AMREX_SPACEDIM == 2) amrex::Array<amrex::Real,3> galilean_shift = { 0., std::numeric_limits<Real>::quiet_NaN(), 0., } ; +#elif (AMREX_SPACEDIM == 1) + amrex::Array<amrex::Real,3> galilean_shift = {std::numeric_limits<Real>::quiet_NaN(), std::numeric_limits<Real>::quiet_NaN(), 0., } ; #endif const amrex::Real zmin_box = WarpX::LowerCorner(bx, galilean_shift, lev)[2]; const amrex::Real zmax_box = WarpX::UpperCorner(bx, lev)[2]; @@ -226,8 +230,10 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) // Get box lower index in the z direction #if (AMREX_SPACEDIM==3) const int lo_ind = bx.loVect()[2]; -#else +#elif (AMREX_SPACEDIM==2) const int lo_ind = bx.loVect()[1]; +#else + const int lo_ind = bx.loVect()[0]; #endif // Check if box intersect with [zmin, zmax] if ( (zmax>zmin_box && zmin<=zmax_box) ){ @@ -237,8 +243,10 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept{ #if (AMREX_SPACEDIM==3) const Real z_gridpoint = zmin_box+(k-lo_ind)*dz; -#else +#elif (AMREX_SPACEDIM==2) const Real z_gridpoint = zmin_box+(j-lo_ind)*dz; +#else + const Real z_gridpoint = zmin_box+(i-lo_ind)*dz; #endif if ( (z_gridpoint >= zmin) && (z_gridpoint < zmax) ) { arr(i,j,k) = 0.; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index f41de8b58..5f9106388 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -659,7 +659,9 @@ WarpX::ReadParameters () Vector<int> parse_filter_npass_each_dir(AMREX_SPACEDIM,1); queryArrWithParser(pp_warpx, "filter_npass_each_dir", parse_filter_npass_each_dir, 0, AMREX_SPACEDIM); filter_npass_each_dir[0] = parse_filter_npass_each_dir[0]; +#if (AMREX_SPACEDIM >= 2) filter_npass_each_dir[1] = parse_filter_npass_each_dir[1]; +#endif #if (AMREX_SPACEDIM == 3) filter_npass_each_dir[2] = parse_filter_npass_each_dir[2]; #endif @@ -1432,7 +1434,9 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + amrex::RealVect dx({WarpX::CellSize(lev)[2]}); +#elif (AMREX_SPACEDIM == 2) amrex::RealVect dx = {WarpX::CellSize(lev)[0], WarpX::CellSize(lev)[2]}; #elif (AMREX_SPACEDIM == 3) amrex::RealVect dx = {WarpX::CellSize(lev)[0], WarpX::CellSize(lev)[1], WarpX::CellSize(lev)[2]}; @@ -1493,7 +1497,18 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm amrex::IntVect F_nodal_flag, G_nodal_flag; // Set nodal flags -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + // AMReX convention: x = missing dimension, y = missing dimension, z = only dimension + Ex_nodal_flag = IntVect(1); + Ey_nodal_flag = IntVect(1); + Ez_nodal_flag = IntVect(0); + Bx_nodal_flag = IntVect(0); + By_nodal_flag = IntVect(0); + Bz_nodal_flag = IntVect(1); + jx_nodal_flag = IntVect(1); + jy_nodal_flag = IntVect(1); + jz_nodal_flag = IntVect(0); +#elif (AMREX_SPACEDIM == 2) // AMReX convention: x = first dimension, y = missing dimension, z = second dimension Ex_nodal_flag = IntVect(0,1); Ey_nodal_flag = IntVect(1,1); @@ -2041,7 +2056,7 @@ WarpX::CellSize (int lev) #elif (AMREX_SPACEDIM == 2) return { dx[0], 1.0, dx[1] }; #else - static_assert(AMREX_SPACEDIM != 1, "1D is not supported"); + return { 1.0, 1.0, dx[0] }; #endif } @@ -2065,6 +2080,9 @@ WarpX::LowerCorner(const Box& bx, std::array<amrex::Real,3> galilean_shift, int #elif (AMREX_SPACEDIM == 2) return { xyzmin[0] + galilean_shift[0], std::numeric_limits<Real>::lowest(), xyzmin[1] + galilean_shift[2] }; + +#elif (AMREX_SPACEDIM == 1) + return { std::numeric_limits<Real>::lowest(), std::numeric_limits<Real>::lowest(), xyzmin[0] + galilean_shift[2] }; #endif } @@ -2077,6 +2095,8 @@ WarpX::UpperCorner(const Box& bx, int lev) return { xyzmax[0], xyzmax[1], xyzmax[2] }; #elif (AMREX_SPACEDIM == 2) return { xyzmax[0], std::numeric_limits<Real>::max(), xyzmax[1] }; +#elif (AMREX_SPACEDIM == 1) + return { std::numeric_limits<Real>::max(), std::numeric_limits<Real>::max(), xyzmax[0] }; #endif } @@ -2276,7 +2296,19 @@ WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_ma const auto hi = amrex::ubound( tbx ); Array4<int> msk = buffer_mask.array(); Array4<int const> gmsk = guard_mask.array(); -#if (AMREX_SPACEDIM == 2) +#if (AMREX_SPACEDIM == 1) + int k = lo.z; + int j = lo.y; + for (int i = lo.x; i <= hi.x; ++i) { + setnull = false; + // If gmsk=0 for any neighbor within ng cells, current cell is in the buffer region + for (int ii = i-ng; ii <= i+ng; ++ii) { + if ( gmsk(ii,j,k) == 0 ) setnull = true; + } + if ( setnull ) msk(i,j,k) = 0; + else msk(i,j,k) = 1; + } +#elif (AMREX_SPACEDIM == 2) int k = lo.z; for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { |