From 24e3acabb41d915b7d5d61f633e676590a0ca15e Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 24 May 2019 14:58:43 -0700 Subject: Update Python wrapper to handle RZ version --- Source/Python/WarpXWrappers.cpp | 46 ++++++++++++++++++++++------------------- 1 file changed, 25 insertions(+), 21 deletions(-) (limited to 'Source/Python/WarpXWrappers.cpp') diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 16d7cd841..d39873c67 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -10,11 +10,14 @@ namespace { - double** getMultiFabPointers(const amrex::MultiFab& mf, int *num_boxes, int *ngrow, int **shapes) + double** getMultiFabPointers(const amrex::MultiFab& mf, int *num_boxes, int *ncomps, int *ngrow, int **shapes) { + *ncomps = mf.nComp(); *ngrow = mf.nGrow(); *num_boxes = mf.local_size(); - *shapes = (int*) malloc(AMREX_SPACEDIM * (*num_boxes) * sizeof(int)); + int shapesize = AMREX_SPACEDIM; + if (mf.nComp() > 1) shapesize += 1; + *shapes = (int*) malloc(shapesize * (*num_boxes) * sizeof(int)); double** data = (double**) malloc((*num_boxes) * sizeof(double*)); int i = 0; @@ -24,8 +27,9 @@ namespace for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi, ++i ) { data[i] = (double*) mf[mfi].dataPtr(); for (int j = 0; j < AMREX_SPACEDIM; ++j) { - (*shapes)[AMREX_SPACEDIM*i+j] = mf[mfi].box().length(j); + (*shapes)[shapesize*i+j] = mf[mfi].box().length(j); } + if (mf.nComp() > 1) (*shapes)[shapesize*i+2] = mf.nComp(); } return data; } @@ -197,9 +201,9 @@ extern "C" } double** warpx_getEfield(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getEfield(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getEfieldLoVects(int lev, int direction, @@ -209,9 +213,9 @@ extern "C" } double** warpx_getEfieldCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getEfieldCPLoVects(int lev, int direction, @@ -221,9 +225,9 @@ extern "C" } double** warpx_getEfieldFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getEfieldFPLoVects(int lev, int direction, @@ -233,9 +237,9 @@ extern "C" } double** warpx_getBfield(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getBfield(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getBfieldLoVects(int lev, int direction, @@ -245,9 +249,9 @@ extern "C" } double** warpx_getBfieldCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getBfieldCPLoVects(int lev, int direction, @@ -257,9 +261,9 @@ extern "C" } double** warpx_getBfieldFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getBfieldFPLoVects(int lev, int direction, @@ -269,9 +273,9 @@ extern "C" } double** warpx_getCurrentDensity(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getcurrent(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getCurrentDensityLoVects(int lev, int direction, @@ -281,9 +285,9 @@ extern "C" } double** warpx_getCurrentDensityCP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getCurrentDensityCPLoVects(int lev, int direction, @@ -293,9 +297,9 @@ extern "C" } double** warpx_getCurrentDensityFP(int lev, int direction, - int *return_size, int *ngrow, int **shapes) { + int *return_size, int *ncomps, int *ngrow, int **shapes) { auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ngrow, shapes); + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); } int* warpx_getCurrentDensityFPLoVects(int lev, int direction, -- cgit v1.2.3 From 6fcf8b97a6bb432b08c456a275d551d53b998be6 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Mon, 26 Aug 2019 09:59:33 -0700 Subject: Clean up of RZ mutlimode, mostly removing imaginary part of mode 0 --- .../langmuir_PICMI_rz_multimode_analyze.py | 4 +- Source/Diagnostics/FieldIO.cpp | 4 +- Source/Evolve/WarpXEvolveEM.cpp | 8 ++-- Source/FieldSolver/WarpXPushFieldsEM.cpp | 16 +++---- Source/FieldSolver/WarpX_FDTD.H | 52 +++++++++++----------- Source/Particles/Deposition/CurrentDeposition.H | 12 ++--- Source/Particles/Gather/FieldGather.H | 24 +++++----- Source/Python/WarpXWrappers.cpp | 2 +- Source/WarpX.cpp | 13 +++--- 9 files changed, 67 insertions(+), 68 deletions(-) (limited to 'Source/Python/WarpXWrappers.cpp') diff --git a/Examples/Tests/Langmuir/langmuir_PICMI_rz_multimode_analyze.py b/Examples/Tests/Langmuir/langmuir_PICMI_rz_multimode_analyze.py index 2a7b317a6..6e63040fe 100644 --- a/Examples/Tests/Langmuir/langmuir_PICMI_rz_multimode_analyze.py +++ b/Examples/Tests/Langmuir/langmuir_PICMI_rz_multimode_analyze.py @@ -145,8 +145,8 @@ Ex_sim_modes = ExWrapper()[...] Ez_sim_modes = EzWrapper()[...] # Sum the real components to get the field along x-axis (theta = 0) -Er_sim = np.sum(Ex_sim_modes[:,:,::2], axis=2) -Ez_sim = np.sum(Ez_sim_modes[:,:,::2], axis=2) +Er_sim = Ex_sim_modes[:,:,0] + np.sum(Ex_sim_modes[:,:,1::2], axis=2) +Ez_sim = Ez_sim_modes[:,:,0] + np.sum(Ez_sim_modes[:,:,1::2], axis=2) # The analytical solutions Er_th = calcEr(zz[:-1,:], rr[:-1,:] + dr/2., k0, w0, wp, t0, [epsilon0, epsilon1, epsilon2]) diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index c5b7a9185..82378edeb 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -219,7 +219,7 @@ ConstructTotalRZField(std::array< std::unique_ptr, 3 >& mf_total, MultiFab::Copy(*mf_total[0], *vector_field[0], 0, 0, 1, vector_field[0]->nGrowVect()); MultiFab::Copy(*mf_total[1], *vector_field[1], 0, 0, 1, vector_field[1]->nGrowVect()); MultiFab::Copy(*mf_total[2], *vector_field[2], 0, 0, 1, vector_field[2]->nGrowVect()); - for (int ic=2 ; ic < vector_field[0]->nComp() ; ic += 2) { + for (int ic=1 ; ic < vector_field[0]->nComp() ; ic += 2) { MultiFab::Add(*mf_total[0], *vector_field[0], ic, 0, 1, vector_field[0]->nGrowVect()); MultiFab::Add(*mf_total[1], *vector_field[1], ic, 0, 1, vector_field[1]->nGrowVect()); MultiFab::Add(*mf_total[2], *vector_field[2], ic, 0, 1, vector_field[2]->nGrowVect()); @@ -448,7 +448,7 @@ WarpX::AverageAndPackFields ( Vector& varnames, // imaginary part of mode 0 for code symmetry, even though // it is always zero. int modes_factor = 1; - if (n_rz_azimuthal_modes > 1) modes_factor = 2*n_rz_azimuthal_modes + 1; + /* if (n_rz_azimuthal_modes > 1) modes_factor = 2*n_rz_azimuthal_modes; */ // Count how many different fields should be written (ncomp) const int ncomp = fields_to_plot.size() diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 7573b572e..5313fb418 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -356,7 +356,7 @@ WarpX::OneStep_sub1 (Real curtime) RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); NodalSyncJ(fine_lev, PatchType::fine); - ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*n_rz_azimuthal_modes); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, rho_fp[fine_lev]->nComp()); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -383,7 +383,7 @@ WarpX::OneStep_sub1 (Real curtime) PushParticlesandDepose(coarse_lev, curtime); StoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, n_rz_azimuthal_modes); + AddRhoFromFineLevelandSumBoundary(coarse_lev, 0, rho_cp[coarse_lev]->nComp()/2); EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]); EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf); @@ -410,7 +410,7 @@ WarpX::OneStep_sub1 (Real curtime) RestrictRhoFromFineToCoarsePatch(fine_lev); ApplyFilterandSumBoundaryJ(fine_lev, PatchType::fine); NodalSyncJ(fine_lev, PatchType::fine); - ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, 2*n_rz_azimuthal_modes); + ApplyFilterandSumBoundaryRho(fine_lev, PatchType::fine, 0, rho_fp[fine_lev]->nComp()/2); NodalSyncRho(fine_lev, PatchType::fine, 0, 2); EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]); @@ -436,7 +436,7 @@ WarpX::OneStep_sub1 (Real curtime) // by only half a coarse step (second half) RestoreCurrent(coarse_lev); AddCurrentFromFineLevelandSumBoundary(coarse_lev); - AddRhoFromFineLevelandSumBoundary(coarse_lev, n_rz_azimuthal_modes, n_rz_azimuthal_modes); + AddRhoFromFineLevelandSumBoundary(coarse_lev, rho_cp[coarse_lev]->nComp()/2, rho_cp[coarse_lev]->nComp()/2); EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]); FillBoundaryE(fine_lev, PatchType::coarse); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 2af6e9e24..b876f2751 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -627,14 +627,14 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Note that Jr(i==0) is at 1/2 dr. if (rmin == 0. && 0 <= i && i < ngJ) { + Jr_arr(i,j,0,2*imode-1) -= ifact*Jr_arr(-1-i,j,0,2*imode-1); Jr_arr(i,j,0,2*imode) -= ifact*Jr_arr(-1-i,j,0,2*imode); - Jr_arr(i,j,0,2*imode+1) -= ifact*Jr_arr(-1-i,j,0,2*imode+1); } // Apply the inverse volume scaling // Since Jr is not cell centered in r, no need for distinction // between on axis and off-axis factors + Jr_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); Jr_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); - Jr_arr(i,j,0,2*imode+1) /= (2.*MathConst::pi*r); } }, [=] AMREX_GPU_DEVICE (int i, int j, int k) @@ -661,18 +661,18 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Jt is located on the boundary if (rmin == 0. && 0 < i && i <= ngJ) { + Jt_arr(i,j,0,2*imode-1) += ifact*Jt_arr(-i,j,0,2*imode-1); Jt_arr(i,j,0,2*imode) += ifact*Jt_arr(-i,j,0,2*imode); - Jt_arr(i,j,0,2*imode+1) += ifact*Jt_arr(-i,j,0,2*imode+1); } // Apply the inverse volume scaling // Jt is forced to zero on axis. if (r == 0.) { + Jt_arr(i,j,0,2*imode-1) = 0.; Jt_arr(i,j,0,2*imode) = 0.; - Jt_arr(i,j,0,2*imode+1) = 0.; } else { + Jt_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); Jt_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); - Jt_arr(i,j,0,2*imode+1) /= (2.*MathConst::pi*r); } } }, @@ -700,18 +700,18 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // to the cells above the axis. // Jz is located on the boundary if (rmin == 0. && 0 < i && i <= ngJ) { + Jz_arr(i,j,0,2*imode-1) += ifact*Jz_arr(-i,j,0,2*imode-1); Jz_arr(i,j,0,2*imode) += ifact*Jz_arr(-i,j,0,2*imode); - Jz_arr(i,j,0,2*imode+1) += ifact*Jz_arr(-i,j,0,2*imode+1); } // Apply the inverse volume scaling if (r == 0.) { // Verboncoeur JCP 164, 421-427 (2001) : corrected volume on axis + Jz_arr(i,j,0,2*imode-1) /= (MathConst::pi*dr/3.); Jz_arr(i,j,0,2*imode) /= (MathConst::pi*dr/3.); - Jz_arr(i,j,0,2*imode+1) /= (MathConst::pi*dr/3.); } else { + Jz_arr(i,j,0,2*imode-1) /= (2.*MathConst::pi*r); Jz_arr(i,j,0,2*imode) /= (2.*MathConst::pi*r); - Jz_arr(i,j,0,2*imode+1) /= (2.*MathConst::pi*r); } } diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index 650e222e1..51e8c46a5 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -28,21 +28,21 @@ void warpx_push_bx_yee(int i, int j, int k, Array4 const& Bx, // (due to the 1/r terms). The following expressions regularize // these divergences by assuming, on axis : // Ez/r = 0/r + dEz/dr - Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) - imode*dtsdx*Ez(i+1,j,0,2*imode+1) & + Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i+1,j,0,2*imode) & + + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); + Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i+1,j,0,2*imode-1) & + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); - Bx(i,j,0,2*imode+1) = Bx(i,j,0,2*imode+1) + imode*dtsdx*Ez(i+1,j,0,2*imode) & - + dtsdz*(Ey(i,j+1,0,2*imode+1) - Ey(i,j,0,2*imode+1)); } else { + Bx(i,j,0,2*imode-1) = 0.; Bx(i,j,0,2*imode) = 0.; - Bx(i,j,0,2*imode+1) = 0.; } } else { // Br(i,j,m) = Br(i,j,m) + I*m*dt*Ez(i,j,m)/r + dtsdz*(Et(i,j+1,m) - Et(i,j,m)) const Real r = rmin*dxinv + i; - Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) - imode*dtsdx*Ez(i,j,0,2*imode+1)/r & + Bx(i,j,0,2*imode-1) = Bx(i,j,0,2*imode-1) - imode*dtsdx*Ez(i,j,0,2*imode)/r & + + dtsdz*(Ey(i,j+1,0,2*imode-1) - Ey(i,j,0,2*imode-1)); + Bx(i,j,0,2*imode) = Bx(i,j,0,2*imode) + imode*dtsdx*Ez(i,j,0,2*imode-1)/r & + dtsdz*(Ey(i,j+1,0,2*imode) - Ey(i,j,0,2*imode)); - Bx(i,j,0,2*imode+1) = Bx(i,j,0,2*imode+1) + imode*dtsdx*Ez(i,j,0,2*imode)/r & - + dtsdz*(Ey(i,j+1,0,2*imode+1) - Ey(i,j,0,2*imode+1)); } } #endif @@ -64,10 +64,10 @@ void warpx_push_by_yee(int i, int j, int k, Array4 const& By, #if (defined WARPX_DIM_RZ) for (int imode=1 ; imode < nmodes ; imode++) { // Bt(i,j,m) = Bt(i,j,m) + dtsdr*(Ez(i+1,j,m) - Ez(i,j,m)) - dtsdz*(Er(i,j+1,m) - Er(i,j,m)) + By(i,j,0,2*imode-1) += + dtsdx*(Ez(i+1,j ,0,2*imode-1) - Ez(i,j,0,2*imode-1)) + - dtsdz*(Ex(i ,j+1,0,2*imode-1) - Ex(i,j,0,2*imode-1)); By(i,j,0,2*imode) += + dtsdx*(Ez(i+1,j ,0,2*imode) - Ez(i,j,0,2*imode)) - - dtsdz*(Ex(i ,j+1,0,2*imode) - Ex(i,j,0,2*imode)); - By(i,j,0,2*imode+1) += + dtsdx*(Ez(i+1,j ,0,2*imode+1) - Ez(i,j,0,2*imode+1)) - - dtsdz*(Ex(i ,j+1,0,2*imode+1) - Ex(i,j,0,2*imode+1)); + - dtsdz*(Ex(i ,j+1,0,2*imode) - Ex(i,j,0,2*imode)); } #endif #endif @@ -90,8 +90,8 @@ void warpx_push_bz_yee(int i, int j, int k, Array4 const& Bz, Bz(i,j,0,0) += - dtsdx*(ru*Ey(i+1,j,0,0) - rd*Ey(i,j,0,0)); for (int imode=1 ; imode < nmodes ; imode++) { // Bz(i,j,m) = Bz(i,j,m) - dtsdr*(ru*Et(i+1,j,m) - rd*Et(i,j,m)) - I*m*dt*Er(i,j,m)/r - Bz(i,j,0,2*imode) += - dtsdx*(ru*Ey(i+1,j,0,2*imode) - rd*Ey(i,j,0,2*imode)) + imode*dtsdx*Ex(i,j,0,2*imode+1)/r; - Bz(i,j,0,2*imode+1) += - dtsdx*(ru*Ey(i+1,j,0,2*imode+1) - rd*Ey(i,j,0,2*imode+1)) - imode*dtsdx*Ex(i,j,0,2*imode)/r; + Bz(i,j,0,2*imode-1) += - dtsdx*(ru*Ey(i+1,j,0,2*imode-1) - rd*Ey(i,j,0,2*imode-1)) + imode*dtsdx*Ex(i,j,0,2*imode)/r; + Bz(i,j,0,2*imode) += - dtsdx*(ru*Ey(i+1,j,0,2*imode) - rd*Ey(i,j,0,2*imode)) - imode*dtsdx*Ex(i,j,0,2*imode-1)/r; } #endif } @@ -113,10 +113,10 @@ void warpx_push_ex_yee(int i, int j, int k, Array4 const& Ex, const Real r = rmin*dxinv+ i + 0.5; for (int imode=1 ; imode < nmodes ; imode++) { // Er(i,j,m) = Er(i,j,m) - I*m*dt*Bz(i,j,m)/r - dtsdz*(Bt(i,j,m) - Bt(i,j-1,m)) - mudt*Jr(i,j,m) - Ex(i,j,0,2*imode) += - dtsdz_c2*(By(i,j,0,2*imode) - By(i,j-1,0,2*imode)) + imode*dtsdx_c2*Bz(i,j,0,2*imode+1)/r + Ex(i,j,0,2*imode-1) += - dtsdz_c2*(By(i,j,0,2*imode-1) - By(i,j-1,0,2*imode-1)) + imode*dtsdx_c2*Bz(i,j,0,2*imode)/r + - mu_c2_dt*Jx(i,j,0,2*imode-1); + Ex(i,j,0,2*imode) += - dtsdz_c2*(By(i,j,0,2*imode) - By(i,j-1,0,2*imode)) - imode*dtsdx_c2*Bz(i,j,0,2*imode-1)/r - mu_c2_dt*Jx(i,j,0,2*imode); - Ex(i,j,0,2*imode+1) += - dtsdz_c2*(By(i,j,0,2*imode+1) - By(i,j-1,0,2*imode+1)) - imode*dtsdx_c2*Bz(i,j,0,2*imode)/r - - mu_c2_dt*Jx(i,j,0,2*imode+1); } #endif #endif @@ -156,21 +156,21 @@ void warpx_push_ey_yee(int i, int j, int k, Array4 const& Ey, // Er(r=0,m=1) = 0.5*[Er(r=dr/2,m=1) + Er(r=-dr/2,m=1)] // And using the rule applying for the guards cells // Er(r=-dr/2,m=1) = Er(r=dr/2,m=1). Thus: Et(i,j,m) = -i*Er(i,j,m) - Ey(i,j,0,2*imode) = Ex(i,j,0,2*imode+1); - Ey(i,j,0,2*imode+1) = -Ex(i,j,0,2*imode); + Ey(i,j,0,2*imode-1) = Ex(i,j,0,2*imode); + Ey(i,j,0,2*imode) = -Ex(i,j,0,2*imode-1); } else { // Etheta should remain 0 on axis, for modes different than m=1 + Ey(i,j,0,2*imode-1) = 0.; Ey(i,j,0,2*imode) = 0.; - Ey(i,j,0,2*imode+1) = 0.; } } else { // Et(i,j,m) = Et(i,j,m) - dtsdr*(Bz(i,j,m) - Bz(i-1,j,m)) + dtsdz*(Br(i,j,m) - Br(i,j-1,m)) - mudt*Jt(i,j,m) + Ey(i,j,0,2*imode-1) += - dtsdx_c2 * (Bz(i,j,0,2*imode-1) - Bz(i-1,j,0,2*imode-1)) + + dtsdz_c2 * (Bx(i,j,0,2*imode-1) - Bx(i,j-1,0,2*imode-1)) + - mu_c2_dt * Jy(i,j,0,2*imode-1); Ey(i,j,0,2*imode) += - dtsdx_c2 * (Bz(i,j,0,2*imode) - Bz(i-1,j,0,2*imode)) + dtsdz_c2 * (Bx(i,j,0,2*imode) - Bx(i,j-1,0,2*imode)) - mu_c2_dt * Jy(i,j,0,2*imode); - Ey(i,j,0,2*imode+1) += - dtsdx_c2 * (Bz(i,j,0,2*imode+1) - Bz(i-1,j,0,2*imode+1)) - + dtsdz_c2 * (Bx(i,j,0,2*imode+1) - Bx(i,j-1,0,2*imode+1)) - - mu_c2_dt * Jy(i,j,0,2*imode+1); } } #endif @@ -200,19 +200,19 @@ void warpx_push_ez_yee(int i, int j, int k, Array4 const& Ez, } for (int imode=1 ; imode < nmodes ; imode++) { if (i == 0 && rmin == 0) { + Ez(i,j,0,2*imode-1) = 0.; Ez(i,j,0,2*imode) = 0.; - Ez(i,j,0,2*imode+1) = 0.; } else { const Real r = rmin*dxinv + i + 0.5; const Real ru = 1. + 0.5/(rmin*dxinv + i); const Real rd = 1. - 0.5/(rmin*dxinv + i); // Ez(i,j,m) = Ez(i,j,m) + dtsdr*(ru*Bt(i,j,m) - rd*Bt(i-1,j,m)) + I*m*dt*Br(i,j,m)/r - mudt*Jz(i,j,m) + Ez(i,j,0,2*imode-1) += + dtsdx_c2 * (ru*By(i,j,0,2*imode-1) - rd*By(i-1,j,0,2*imode-1)) + - imode*dtsdx_c2*Bx(i,j,0,2*imode)/r + - mu_c2_dt * Jz(i,j,0,2*imode-1); Ez(i,j,0,2*imode) += + dtsdx_c2 * (ru*By(i,j,0,2*imode) - rd*By(i-1,j,0,2*imode)) - - imode*dtsdx_c2*Bx(i,j,0,2*imode+1)/r + + imode*dtsdx_c2*Bx(i,j,0,2*imode-1)/r - mu_c2_dt * Jz(i,j,0,2*imode); - Ez(i,j,0,2*imode+1) += + dtsdx_c2 * (ru*By(i,j,0,2*imode+1) - rd*By(i-1,j,0,2*imode+1)) - + imode*dtsdx_c2*Bx(i,j,0,2*imode)/r - - mu_c2_dt * Jz(i,j,0,2*imode+1); } } #endif diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index f736d6af0..068a6e190 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -398,8 +398,8 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes const Complex djr_cmplx = 2.*sdxi*xy_mid; - amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.real()); - amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode+1), djr_cmplx.imag()); + amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real()); + amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag()); xy_mid = xy_mid*xy_mid0; } #endif @@ -420,8 +420,8 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, // The minus sign comes from the different convention with respect to Davidson et al. const Complex djt_cmplx = -2.*I*(lo.x+i_new-1 + i + xmin*dxi)*wq*invdtdx/(double)imode* (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old)); - amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.real()); - amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode+1), djt_cmplx.imag()); + amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); + amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag()); xy_new = xy_new*xy_new0; xy_mid = xy_mid*xy_mid0; xy_old = xy_old*xy_old0; @@ -439,8 +439,8 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes const Complex djz_cmplx = 2.*sdzk*xy_mid; - amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.real()); - amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode+1), djz_cmplx.imag()); + amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real()); + amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag()); xy_mid = xy_mid*xy_mid0; } #endif diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index eea637aa9..d8d1d78ef 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -156,8 +156,8 @@ void doGatherShapeN(const amrex::Real * const xp, // Gather field on particle Eyp[i] from field on grid ey_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ - const amrex::Real dEy = (+ ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode)*xy.real() - - ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode+1)*xy.imag()); + const amrex::Real dEy = (+ ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode-1)*xy.real() + - ey_arr(lo.x+j+ix, lo.y+l+iz, 0, 2*imode)*xy.imag()); Eyp[ip] += sx[ix]*sz[iz]*dEy; } } @@ -165,11 +165,11 @@ void doGatherShapeN(const amrex::Real * const xp, // Gather field on particle Bzp[i] from field on grid bz_arr for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - const amrex::Real dEx = (+ ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.real() - - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode+1)*xy.imag()); + const amrex::Real dEx = (+ ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode-1)*xy.real() + - ex_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.imag()); Exp[ip] += sx0[ix]*sz[iz]*dEx; - const amrex::Real dBz = (+ bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.real() - - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode+1)*xy.imag()); + const amrex::Real dBz = (+ bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode-1)*xy.real() + - bz_arr(lo.x+j0+ix, lo.y+l +iz, 0, 2*imode)*xy.imag()); Bzp[ip] += sx0[ix]*sz[iz]*dBz; } } @@ -177,19 +177,19 @@ void doGatherShapeN(const amrex::Real * const xp, // Gather field on particle Bxp[i] from field on grid bx_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order; ix++){ - const amrex::Real dEz = (+ ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.real() - - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode+1)*xy.imag()); + const amrex::Real dEz = (+ ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode-1)*xy.real() + - ez_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.imag()); Ezp[ip] += sx[ix]*sz0[iz]*dEz; - const amrex::Real dBx = (+ bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.real() - - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode+1)*xy.imag()); + const amrex::Real dBx = (+ bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode-1)*xy.real() + - bx_arr(lo.x+j+ix, lo.y+l0 +iz, 0, 2*imode)*xy.imag()); Bxp[ip] += sx[ix]*sz0[iz]*dBx; } } // Gather field on particle Byp[i] from field on grid by_arr for (int iz=0; iz<=depos_order-lower_in_v; iz++){ for (int ix=0; ix<=depos_order-lower_in_v; ix++){ - const amrex::Real dBy = (+ by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode)*xy.real() - - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode+1)*xy.imag()); + const amrex::Real dBy = (+ by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode-1)*xy.real() + - by_arr(lo.x+j0+ix, lo.y+l0+iz, 0, 2*imode)*xy.imag()); Byp[ip] += sx0[ix]*sz0[iz]*dBy; } } diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 3ed4830f5..a60efd498 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -29,7 +29,7 @@ namespace for (int j = 0; j < AMREX_SPACEDIM; ++j) { (*shapes)[shapesize*i+j] = mf[mfi].box().length(j); } - if (mf.nComp() > 1) (*shapes)[shapesize*i+2] = mf.nComp(); + if (mf.nComp() > 1) (*shapes)[shapesize*i+AMREX_SPACEDIM] = mf.nComp(); } return data; } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 0b6e93a95..8e8c25f4d 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -787,13 +787,12 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm { #if defined WARPX_DIM_RZ - if (n_rz_azimuthal_modes > 1) { - // There is a real and imaginary component for each mode - ncomps = n_rz_azimuthal_modes*2; - } else { - // With only mode 0, only reals are used - ncomps = 1; - } + // With RZ multimode, there is a real and imaginary component + // for each mode, except mode 0 which is purely real + // Component 0 is mode 0. + // Odd components are the real parts. + // Even components are the imaginary parts. + ncomps = n_rz_azimuthal_modes*2 - 1; #endif // -- cgit v1.2.3