diff options
12 files changed, 338 insertions, 217 deletions
diff --git a/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd.json b/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd.json index 641457e04..0589f213f 100644 --- a/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd.json +++ b/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd.json @@ -2,32 +2,32 @@ "electrons": { "particle_cpu": 0.0, "particle_id": 2147516416.0, - "particle_momentum_x": 2.532371210270441e-21, - "particle_momentum_y": 2.699519499032127e-21, - "particle_momentum_z": 1.780767630883076e-16, - "particle_position_x": 4.055885890661150e+05, - "particle_position_y": 2.012710911909795e+07, + "particle_momentum_x": 2.532315559872013e-21, + "particle_momentum_y": 2.6999645401303122e-21, + "particle_momentum_z": 1.7807676308672323e-16, + "particle_position_x": 405588.5890661151, + "particle_position_y": 20127109.119097948, "particle_weight": 6.917460794691972e+17 }, "ions": { "particle_cpu": 0.0, "particle_id": 6442483712.0, - "particle_momentum_x": 2.609374760931462e-18, - "particle_momentum_y": 2.617947652774628e-18, - "particle_momentum_z": 3.269760998813916e-13, - "particle_position_x": 405588.43859727494, - "particle_position_y": 20127109.117509525, + "particle_momentum_x": 2.6093747608241028e-18, + "particle_momentum_y": 2.61794770341563e-18, + "particle_momentum_z": 3.269760999095121e-13, + "particle_position_x": 405588.4386175681, + "particle_position_y": 20127109.117484942, "particle_weight": 6.917460794691972e+17 }, "lev=0": { - "Bx": 2.836797963533862e-03, - "By": 1.587110142217534e-03, - "Bz": 7.993802082803127e-03, - "Ex": 4.830436110421224e+05, - "Ey": 2.565122701641558e+06, - "Ez": 4.462554901421999e+04, - "jx": 2.196458000289692e+02, - "jy": 9.854589204209451e+02, - "jz": 4.147941958892408e+03 + "Bx": 0.002836797963533812, + "By": 0.0015871101422500937, + "Bz": 0.007993802082802991, + "Ex": 483043.6110477062, + "Ey": 2565122.7016415303, + "Ez": 44625.54901861734, + "jx": 219.64580002836755, + "jy": 985.4589204209067, + "jz": 4147.941958882783 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json index 473d072bc..2df4d7266 100644 --- a/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json +++ b/Regression/Checksum/benchmarks_json/averaged_galilean_2d_psatd_hybrid.json @@ -2,9 +2,9 @@ "electrons": { "particle_cpu": 32768.0, "particle_id": 1073774592.0, - "particle_momentum_x": 2.6991465744589784e-21, - "particle_momentum_y": 2.3548757050745318e-21, - "particle_momentum_z": 1.7807675042822086e-16, + "particle_momentum_x": 2.6987960602750912e-21, + "particle_momentum_y": 2.3583158815448188e-21, + "particle_momentum_z": 1.7807675042654014e-16, "particle_position_x": 405588.6263336367, "particle_position_y": 20127109.123049334, "particle_weight": 6.917460794691972e+17 @@ -12,20 +12,20 @@ "ions": { "particle_cpu": 32768.0, "particle_id": 3221258240.0, - "particle_momentum_x": 2.6368707684152848e-18, - "particle_momentum_y": 2.6255385742132203e-18, - "particle_momentum_z": 3.269761087904483e-13, + "particle_momentum_x": 2.6368707637661792e-18, + "particle_momentum_y": 2.6255383353049537e-18, + "particle_momentum_z": 3.2697610879044984e-13, "particle_position_x": 405588.6872874654, "particle_position_y": 20127109.122936208, "particle_weight": 6.917460794691972e+17 }, "lev=0": { - "Bx": 0.002203710933841355, - "By": 0.004460064079572663, - "Bz": 0.00636094468142458, - "Ex": 1416289.0559212451, + "Bx": 0.002232320213633506, + "By": 0.004459938776538334, + "Bz": 0.006374599526320067, + "Ex": 1416337.8255897765, "Ey": 2097441.857526919, - "Ez": 322604.31292830594, + "Ez": 322659.2369253079, "jx": 434.8052109005167, "jy": 844.9945365202007, "jz": 5016.6719265439715 diff --git a/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd.json b/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd.json index 69e5e47ee..a93cd9e02 100644 --- a/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd.json +++ b/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd.json @@ -2,8 +2,8 @@ "electrons": { "particle_cpu": 0.0, "particle_id": 536887296.0, - "particle_momentum_x": 3.235085756863653e-20, - "particle_momentum_y": 3.119669915093045e-20, + "particle_momentum_x": 3.235370594842269e-20, + "particle_momentum_y": 3.116556644010673e-20, "particle_momentum_z": 8.904212408978889e-17, "particle_position_x": 158433.3678817281, "particle_position_y": 158432.165536343, @@ -13,8 +13,8 @@ "ions": { "particle_cpu": 0.0, "particle_id": 1610629120.0, - "particle_momentum_x": 1.3144416374780757e-18, - "particle_momentum_y": 1.3116365676113791e-18, + "particle_momentum_x": 1.314439943956008e-18, + "particle_momentum_y": 1.311628100016085e-18, "particle_momentum_z": 1.6348803463772493e-13, "particle_position_x": 158433.36794743972, "particle_position_y": 158432.1305546865, @@ -32,4 +32,4 @@ "jy": 69843.47286996675, "jz": 21538.61066907263 } -}
\ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json b/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json index 545d4cdf7..9881b0b72 100644 --- a/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json +++ b/Regression/Checksum/benchmarks_json/averaged_galilean_3d_psatd_hybrid.json @@ -2,9 +2,9 @@ "electrons": { "particle_cpu": 0.0, "particle_id": 536887296.0, - "particle_momentum_x": 3.1129132496015666e-20, - "particle_momentum_y": 3.0017627740521483e-20, - "particle_momentum_z": 8.904171965856788e-17, + "particle_momentum_x": 3.069060217914455e-20, + "particle_momentum_y": 2.955864813654792e-20, + "particle_momentum_z": 8.904171964713278e-17, "particle_position_x": 158433.3693892244, "particle_position_y": 158432.16743352715, "particle_position_z": 15724303.927902682, @@ -13,21 +13,21 @@ "ions": { "particle_cpu": 0.0, "particle_id": 1610629120.0, - "particle_momentum_x": 1.3144073041863281e-18, - "particle_momentum_y": 1.3115588299779562e-18, - "particle_momentum_z": 1.6348803504860865e-13, + "particle_momentum_x": 1.3143704893850579e-18, + "particle_momentum_y": 1.311559649208915e-18, + "particle_momentum_z": 1.6348803504847822e-13, "particle_position_x": 158433.36794743096, "particle_position_y": 158432.13055422663, "particle_position_z": 15724303.986784957, "particle_weight": 4.082754265421834e+18 }, "lev=0": { - "Bx": 0.056921656761459924, - "By": 0.058934240410496434, - "Bz": 0.5206831875786353, - "Ex": 171806975.3287552, - "Ey": 164770789.98409316, - "Ez": 1036142.4032854013, + "Bx": 0.05811677943915807, + "By": 0.06003789952009811, + "Bz": 0.53300620716022, + "Ex": 173236565.20945996, + "Ey": 165890551.4406476, + "Ez": 1037164.1567886411, "jx": 76766.28402952384, "jy": 75076.8695684449, "jz": 20606.225852789266 diff --git a/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json b/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json index 80d8f60ac..06c4b8367 100644 --- a/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json +++ b/Regression/Checksum/benchmarks_json/comoving_hybrid_2d.json @@ -2,43 +2,43 @@ "beam": { "particle_cpu": 0.0, "particle_id": 466643.0, - "particle_momentum_x": 0.0000000000000000005940643790503620617444, - "particle_momentum_y": 0.0000000000000000004207935220641095129183, - "particle_momentum_z": 0.0000000000000000090690824780397199942459, - "particle_position_x": 0.0013704760088849208231331910923245231970, - "particle_position_y": 0.1486369608860869107047619763761758804321, + "particle_momentum_x": 5.935248102846229e-19, + "particle_momentum_y": 4.2079455139390086e-19, + "particle_momentum_z": 9.069157096209629e-18, + "particle_position_x": 0.0013709782662344476, + "particle_position_y": 0.14863704108231895, "particle_weight": 2911663983235.947 }, "electrons": { "particle_cpu": 155246.0, "particle_id": 26701298988.0, - "particle_momentum_x": 0.0000000000000000009186987708931365023716, - "particle_momentum_y": 0.0000000000000000008099466808209477065805, - "particle_momentum_z": 0.0000000000000005716177443199551787605836, - "particle_position_x": 6.6343937669634316378619587339926511049271, - "particle_position_y": 58.08977454115495, + "particle_momentum_x": 9.169991045389105e-19, + "particle_momentum_y": 8.134032712971527e-19, + "particle_momentum_z": 5.716870930487574e-16, + "particle_position_x": 6.634385105820577, + "particle_position_y": 58.08979694901561, "particle_weight": 3.212795210423339e+18 }, "ions": { "particle_cpu": 155246.0, "particle_id": 26830093836.0, - "particle_momentum_x": 0.0000000000000000020044237826247600848882, - "particle_momentum_y": 0.0000000000000000008198157862412493956856, - "particle_momentum_z": 1.0089822624453072e-12, - "particle_position_x": 6.622297128431678, - "particle_position_y": 58.08877603287398, + "particle_momentum_x": 2.0028127909408492e-18, + "particle_momentum_y": 8.233828374969269e-19, + "particle_momentum_z": 1.0089821936116652e-12, + "particle_position_x": 6.6222971327574776, + "particle_position_y": 58.088776020644396, "particle_weight": 3.212794492131793e+18 }, "lev=0": { - "Bx": 1301466.3798134964890778064727783203125000000000, - "By": 4006270.2276948434300720691680908203125000000000, - "Bz": 203201.1633351982163731008768081665039062500000, - "Ex": 1203026764868191.75, - "Ey": 407303968803835.0625, - "Ez": 217578142282980.46875, - "jx": 35138099060953276.0, - "jy": 23083112114752576.0, - "jz": 413107254378017984.0, - "rho": 1350522452.8903205394744873046875 + "Bx": 1305305.431740205, + "By": 4004955.3295348613, + "Bz": 203986.21157566475, + "Ex": 1202527548976529.5, + "Ey": 408729463466964.56, + "Ez": 217606643896681.66, + "jx": 3.5156420344831892e+16, + "jy": 2.330006936826398e+16, + "jz": 4.1305946927278195e+17, + "rho": 1350581221.909034 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json b/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json index 5bfebde72..4dc35e93f 100644 --- a/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json +++ b/Regression/Checksum/benchmarks_json/galilean_hybrid_2d.json @@ -2,43 +2,43 @@ "beam": { "particle_cpu": 0.0, "particle_id": 447573.0, - "particle_momentum_x": 5.389906975731078e-19, - "particle_momentum_y": 4.0113509524430207e-19, - "particle_momentum_z": 8.088485048828958e-18, - "particle_position_x": 0.0014081774396371826, - "particle_position_y": 0.14121931180937125, + "particle_momentum_x": 5.384525633353412e-19, + "particle_momentum_y": 4.011357939321231e-19, + "particle_momentum_z": 8.08848605215067e-18, + "particle_position_x": 0.0014094341950741979, + "particle_position_y": 0.1412193837490303, "particle_weight": 2780592292672.2705 }, "electrons": { "particle_cpu": 155383.0, "particle_id": 26721507984.0, - "particle_momentum_x": 7.338476660439086e-19, - "particle_momentum_y": 5.65821279562348e-19, - "particle_momentum_z": 5.696860669282423e-16, - "particle_position_x": 6.639871516463294, - "particle_position_y": 58.21048773913985, - "particle_weight": 3.2159037448973373e+18 + "particle_momentum_x": 9.809438398217426e-19, + "particle_momentum_y": 5.687153095622055e-19, + "particle_momentum_z": 5.697506964559476e-16, + "particle_position_x": 6.639908741842812, + "particle_position_y": 58.210506112500646, + "particle_weight": 3.215903744897337e+18 }, "ions": { "particle_cpu": 155428.0, "particle_id": 26814790002.0, - "particle_momentum_x": 1.6705682162320084e-18, - "particle_momentum_y": 5.723976856025395e-19, - "particle_momentum_z": 1.0101678467302113e-12, - "particle_position_x": 6.630049969920412, - "particle_position_y": 58.2491476342914, + "particle_momentum_x": 1.906616117322184e-18, + "particle_momentum_y": 5.75419216528982e-19, + "particle_momentum_z": 1.0101677927464587e-12, + "particle_position_x": 6.630049949959092, + "particle_position_y": 58.24914762101138, "particle_weight": 3.21662600232034e+18 }, "lev=0": { - "Bx": 1220834.3602474995, - "By": 3462824.6001510955, - "Bz": 194280.5148985213, - "Ex": 1039574206585584.9, - "Ey": 383018547917553.56, - "Ez": 200106779561626.4, - "jx": 3.097332988821253e+16, - "jy": 2.215463603536144e+16, - "jz": 2.921753560471727e+17, - "rho": 994167865.3632829 + "Bx": 1225109.2628475595, + "By": 3748536.8460838096, + "Bz": 195216.2190878017, + "Ex": 1092766647326914.6, + "Ey": 384575347595082.3, + "Ez": 343130879909148.75, + "jx": 3.51564554699967e+16, + "jy": 2.2405130712647624e+16, + "jz": 3.0976785868961254e+17, + "rho": 1054402097.6298534 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json index ad0acf8f7..f3c0b0a48 100644 --- a/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json +++ b/Regression/Checksum/benchmarks_json/pml_psatd_dive_divb_cleaning.json @@ -1,11 +1,11 @@ { "lev=0": { - "Bx": 1.3610687381384125e-08, - "By": 1.361341814453787e-08, - "Bz": 1.361623444282223e-08, - "Ex": 18.733019975918303, - "Ey": 18.73132517981054, - "Ez": 18.731454226529113, + "Bx": 1.55895366473558e-08, + "By": 1.5015206116838625e-08, + "Bz": 1.516725065646249e-08, + "Ex": 19.175719716905533, + "Ey": 19.25088294791974, + "Ez": 19.244587842285785, "rho": 0.00017713857500552797 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/pml_x_psatd.json b/Regression/Checksum/benchmarks_json/pml_x_psatd.json index e3eb32102..f845c7199 100644 --- a/Regression/Checksum/benchmarks_json/pml_x_psatd.json +++ b/Regression/Checksum/benchmarks_json/pml_x_psatd.json @@ -1,12 +1,12 @@ { "lev=0": { - "Bx": 0.0000000115934922595108653400133868777849, - "By": 0.0000000151105874509448789476968307921370, - "Bz": 0.0000000089855317729980750970015230638492, - "Ex": 3.9386851408461347467948598932707682251930, - "Ey": 4.1349222519754311733208851364906877279282, - "Ez": 3.3085920402662600814380766678368672728539, - "divE": 188916.91920602522441186010837554931640625, - "rho": 1.6719346122590767e-06 + "Bx": 1.2110998243414779e-08, + "By": 1.5545485705380978e-08, + "Bz": 9.134211013787592e-09, + "Ex": 4.050591183040409, + "Ey": 4.248614715015242, + "Ez": 3.3401282934200776, + "divE": 188920.76108729007, + "rho": 1.6719346125882405e-06 } -} +}
\ No newline at end of file diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 9c4e786de..138476f88 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -219,6 +219,15 @@ SpectralFieldData::BackwardTransform( const int lev, const bool is_nodal_z = mf.is_nodal(1); #endif + const int si = (is_nodal_x) ? 1 : 0; +#if (AMREX_SPACEDIM == 2) + const int sj = (is_nodal_z) ? 1 : 0; + const int sk = 0; +#elif (AMREX_SPACEDIM == 3) + const int sj = (is_nodal_y) ? 1 : 0; + const int sk = (is_nodal_z) ? 1 : 0; +#endif + // Loop over boxes for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) @@ -256,51 +265,49 @@ SpectralFieldData::BackwardTransform( const int lev, // Copy field into temporary array tmp_arr(i,j,k) = spectral_field_value; }); - } // Perform Fourier transform from `tmpSpectralField` to `tmpRealField` AnyFFT::Execute(backward_plan[mfi]); - // Copy the temporary field `tmpRealField` to the real-space field `mf` - // (only in the valid cells ; not in the guard cells) - // Normalize (divide by 1/N) since the FFT+IFFT results in a factor N + // Copy the temporary field tmpRealField to the real-space field mf and + // normalize, dividing by N, since (FFT + inverse FFT) results in a factor N { - Array4<Real> mf_arr = mf[mfi].array(); - Array4<const Real> tmp_arr = tmpRealField[mfi].array(); - // Normalization: divide by the number of points in realspace - // (includes the guard cells) - const Box realspace_bx = tmpRealField[mfi].box(); - const Real inv_N = 1./realspace_bx.numPts(); - - if (m_periodic_single_box) { - // Enforce periodicity on the nodes, by using modulo in indices - // This is because `tmp_arr` is cell-centered while `mf_arr` can be nodal - int const nx = realspace_bx.length(0); - int const ny = realspace_bx.length(1); -#if (AMREX_SPACEDIM == 3) - int const nz = realspace_bx.length(2); -#else - int constexpr nz = 1; + amrex::Box const& mf_box = (m_periodic_single_box) ? mfi.validbox() : mfi.fabbox(); + amrex::Array4<amrex::Real> mf_arr = mf[mfi].array(); + amrex::Array4<const amrex::Real> tmp_arr = tmpRealField[mfi].array(); + + const amrex::Real inv_N = 1._rt / tmpRealField[mfi].box().numPts(); + + // Total number of cells, including ghost cells (nj represents ny in 3D and nz in 2D) + const int ni = mf_box.length(0); + const int nj = mf_box.length(1); +#if (AMREX_SPACEDIM == 2) + constexpr int nk = 1; +#elif (AMREX_SPACEDIM == 3) + const int nk = mf_box.length(2); #endif - ParallelFor( - mfi.validbox(), - /* GCC 8.1-8.2 work-around (ICE): - * named capture in nonexcept lambda needed for modulo operands - * https://godbolt.org/z/ppbAzd - */ - [mf_arr, i_comp, inv_N, tmp_arr, nx, ny, nz] - AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i%nx, j%ny, k%nz); - }); - } else { - ParallelFor( mfi.validbox(), - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Copy and normalize field - mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k); - }); - } - + // Lower bound of the box (lo_j represents lo_y in 3D and lo_z in 2D) + const int lo_i = amrex::lbound(mf_box).x; + const int lo_j = amrex::lbound(mf_box).y; +#if (AMREX_SPACEDIM == 2) + constexpr int lo_k = 0; +#elif (AMREX_SPACEDIM == 3) + const int lo_k = amrex::lbound(mf_box).z; +#endif + // Loop over cells within full box, including ghost cells + ParallelFor(mf_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Assume periodicity and set the last outer guard cell equal to the first one: + // this is necessary in order to get the correct value along a nodal direction, + // because the last point along a nodal direction is always discarded when FFTs + // are computed, as the real-space box is always cell-centered. + const int ii = (i == lo_i + ni - si) ? lo_i : i; + const int jj = (j == lo_j + nj - sj) ? lo_j : j; + const int kk = (k == lo_k + nk - sk) ? lo_k : k; + // Copy and normalize field + mf_arr(i,j,k,i_comp) = inv_N * tmp_arr(ii,jj,kk); + }); } if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 222a6ae6f..672eaec0c 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -11,6 +11,7 @@ #include "BoundaryConditions/WarpX_PML_kernels.H" #include "BoundaryConditions/PML_current.H" #include "WarpX_FDTD.H" +#include "WarpXPushFieldsEM_K.H" #ifdef BL_USE_SENSEI_INSITU # include <AMReX_AmrMeshInSituBridge.H> @@ -161,8 +162,16 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) { PushPSATDSinglePatch( lev, *spectral_solver_cp[lev], Efield_cp[lev], Bfield_cp[lev], Efield_avg_cp[lev], Bfield_avg_cp[lev], current_cp[lev], rho_cp[lev] ); } - if (use_damp_fields_in_z_guard) { - DampFieldsInGuards( Efield_fp[lev], Bfield_fp[lev] ); + + // Damp the fields in the guard cells along z + if (use_damp_fields_in_z_guard) + { + DampFieldsInGuards(Efield_fp[lev], Bfield_fp[lev]); + + if (WarpX::fft_do_time_averaging) + { + DampFieldsInGuards(Efield_avg_fp[lev], Bfield_avg_fp[lev]); + } } #endif } @@ -419,7 +428,6 @@ WarpX::DampFieldsInGuards(std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield for ( amrex::MFIter mfi(*Efield[0], amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi ) { - amrex::Array4<amrex::Real> const& Ex_arr = Efield[0]->array(mfi); amrex::Array4<amrex::Real> const& Ey_arr = Efield[1]->array(mfi); amrex::Array4<amrex::Real> const& Ez_arr = Efield[2]->array(mfi); @@ -427,74 +435,72 @@ WarpX::DampFieldsInGuards(std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield amrex::Array4<amrex::Real> const& By_arr = Bfield[1]->array(mfi); amrex::Array4<amrex::Real> const& Bz_arr = Bfield[2]->array(mfi); - // Get the tilebox from Efield so that it includes the guard cells. - amrex::Box tilebox = (*Efield[0])[mfi].box(); - int const nz_tile = tilebox.bigEnd(zdir); + // Get the tileboxes from Efield and Bfield so that they include the guard cells + // and take the staggering of each MultiFab into account + amrex::Box tex = amrex::convert((*Efield[0])[mfi].box(), Efield[0]->ixType().toIntVect()); + amrex::Box tey = amrex::convert((*Efield[1])[mfi].box(), Efield[1]->ixType().toIntVect()); + amrex::Box tez = amrex::convert((*Efield[2])[mfi].box(), Efield[2]->ixType().toIntVect()); + amrex::Box tbx = amrex::convert((*Bfield[0])[mfi].box(), Bfield[0]->ixType().toIntVect()); + amrex::Box tby = amrex::convert((*Bfield[1])[mfi].box(), Bfield[1]->ixType().toIntVect()); + amrex::Box tbz = amrex::convert((*Bfield[2])[mfi].box(), Bfield[2]->ixType().toIntVect()); + + // Get smallEnd of tileboxes + const int tex_smallEnd_z = tex.smallEnd(zdir); + const int tey_smallEnd_z = tey.smallEnd(zdir); + const int tez_smallEnd_z = tez.smallEnd(zdir); + const int tbx_smallEnd_z = tbx.smallEnd(zdir); + const int tby_smallEnd_z = tby.smallEnd(zdir); + const int tbz_smallEnd_z = tbz.smallEnd(zdir); + + // Get bigEnd of tileboxes + const int tex_bigEnd_z = tex.bigEnd(zdir); + const int tey_bigEnd_z = tey.bigEnd(zdir); + const int tez_bigEnd_z = tez.bigEnd(zdir); + const int tbx_bigEnd_z = tbx.bigEnd(zdir); + const int tby_bigEnd_z = tby.bigEnd(zdir); + const int tbz_bigEnd_z = tbz.bigEnd(zdir); // Box for the whole simulation domain amrex::Box const& domain = Geom(0).Domain(); int const nz_domain = domain.bigEnd(zdir); - if (tilebox.smallEnd(zdir) < 0) { - - // Apply damping factor in guards cells below the lower end of the domain - int const nz_guard = -tilebox.smallEnd(zdir); + // Set the tileboxes so that they only cover the lower/upper half of the guard cells + constrain_tilebox_to_guards(tex, zdir, nz_domain, tex_smallEnd_z, tex_bigEnd_z); + constrain_tilebox_to_guards(tey, zdir, nz_domain, tey_smallEnd_z, tey_bigEnd_z); + constrain_tilebox_to_guards(tez, zdir, nz_domain, tez_smallEnd_z, tez_bigEnd_z); + constrain_tilebox_to_guards(tbx, zdir, nz_domain, tbx_smallEnd_z, tbx_bigEnd_z); + constrain_tilebox_to_guards(tby, zdir, nz_domain, tby_smallEnd_z, tby_bigEnd_z); + constrain_tilebox_to_guards(tbz, zdir, nz_domain, tbz_smallEnd_z, tbz_bigEnd_z); - // Set so the box only covers the lower half of the guard cells - tilebox.setBig(zdir, -nz_guard/2-1); - - amrex::ParallelFor(tilebox, Efield[0]->nComp(), - [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + amrex::ParallelFor( + tex, Efield[0]->nComp(), [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) { -#if (AMREX_SPACEDIM == 3) - amrex::Real zcell = static_cast<amrex::Real>(k + nz_guard); -#else - amrex::Real zcell = static_cast<amrex::Real>(j + nz_guard); -#endif - const amrex::Real phase = MathConst::pi*zcell/nz_guard; - const amrex::Real sin_phase = std::sin(phase); - const amrex::Real damp_factor = sin_phase*sin_phase; - - Ex_arr(i,j,k,icomp) *= damp_factor; - Ey_arr(i,j,k,icomp) *= damp_factor; - Ez_arr(i,j,k,icomp) *= damp_factor; - Bx_arr(i,j,k,icomp) *= damp_factor; - By_arr(i,j,k,icomp) *= damp_factor; - Bz_arr(i,j,k,icomp) *= damp_factor; - - }); - - } - else if (nz_tile > nz_domain) { - - // Apply damping factor in guards cells above the upper end of the domain - int nz_guard = nz_tile - nz_domain; - - // Set so the box only covers the upper half of the guard cells - tilebox.setSmall(zdir, nz_domain + nz_guard/2 + 1); - - amrex::ParallelFor(tilebox, Efield[0]->nComp(), - [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + damp_field_in_guards(Ex_arr, i, j, k, icomp, zdir, nz_domain, tex_smallEnd_z, tex_bigEnd_z); + }, + tey, Efield[1]->nComp(), [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) { -#if (AMREX_SPACEDIM == 3) - amrex::Real zcell = static_cast<amrex::Real>(nz_tile - k); -#else - amrex::Real zcell = static_cast<amrex::Real>(nz_tile - j); -#endif - const amrex::Real phase = MathConst::pi*zcell/nz_guard; - const amrex::Real sin_phase = std::sin(phase); - const amrex::Real damp_factor = sin_phase*sin_phase; - - Ex_arr(i,j,k,icomp) *= damp_factor; - Ey_arr(i,j,k,icomp) *= damp_factor; - Ez_arr(i,j,k,icomp) *= damp_factor; - Bx_arr(i,j,k,icomp) *= damp_factor; - By_arr(i,j,k,icomp) *= damp_factor; - Bz_arr(i,j,k,icomp) *= damp_factor; - - }); + damp_field_in_guards(Ey_arr, i, j, k, icomp, zdir, nz_domain, tey_smallEnd_z, tey_bigEnd_z); + }, + tez, Efield[2]->nComp(), [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + { + damp_field_in_guards(Ez_arr, i, j, k, icomp, zdir, nz_domain, tez_smallEnd_z, tez_bigEnd_z); + } + ); - } + amrex::ParallelFor( + tbx, Bfield[0]->nComp(), [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + { + damp_field_in_guards(Bx_arr, i, j, k, icomp, zdir, nz_domain, tbx_smallEnd_z, tbx_bigEnd_z); + }, + tby, Bfield[1]->nComp(), [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + { + damp_field_in_guards(By_arr, i, j, k, icomp, zdir, nz_domain, tby_smallEnd_z, tby_bigEnd_z); + }, + tbz, Bfield[2]->nComp(), [=] AMREX_GPU_DEVICE (int i, int j, int k, int icomp) + { + damp_field_in_guards(Bz_arr, i, j, k, icomp, zdir, nz_domain, tbz_smallEnd_z, tbz_bigEnd_z); + } + ); } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM_K.H b/Source/FieldSolver/WarpXPushFieldsEM_K.H new file mode 100644 index 000000000..7f6ca5f69 --- /dev/null +++ b/Source/FieldSolver/WarpXPushFieldsEM_K.H @@ -0,0 +1,109 @@ +/* Copyright 2019-2020 + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WarpXPushFieldsEM_K_h +#define WarpXPushFieldsEM_K_h + +#include "Utils/WarpXConst.H" +#include <AMReX.H> + +/* + * \brief Set a tilebox so that it only covers the lower/upper half of the guard cells. + * + * \param[in,out] tb tilebox to be modified + * \param[in] dir direction where the tilebox smallEnd/bigEnd is modified + * \param[in] n_domain number of cells in the whole simulation domain + * \param[in] tb_smallEnd smallEnd of the tilebox + * \param[in] tb_bigEnd bigEnd of the tilebox + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void constrain_tilebox_to_guards( + amrex::Box& tb, + const int dir, + const int n_domain, + const int tb_smallEnd, + const int tb_bigEnd) +{ + using namespace amrex; + + const int n_tile = tb_bigEnd; + + if (tb_smallEnd < 0) + { + const int n_guard = -tb_smallEnd; + tb.setBig(dir, 0 - n_guard/2 - 1); + } + else if (n_tile > n_domain) + { + const int n_guard = n_tile - n_domain; + tb.setSmall(dir, n_domain + n_guard/2 + 1); + } +} + +/* + * \brief Damp a given field in the guard cells along a given direction + * + * \param[in,out] mf_arr array that contains the field values to be damped + * \oaram[in] i index along x + * \oaram[in] j index along y (in 3D) or z (in 2D/RZ) + * \oaram[in] k index along z (in 3D, \c k = 0 in 2D/RZ) + * \param[in] icomp index along the fourth component of the array + * \param]in] dir direction where the field will be damped + * \param[in] n_domain number of cells in the whole simulation domain + * \param[in] tb_smallEnd smallEnd of the current tilebox + * \param[in] tb_bigEnd bigEnd of the current tilebox + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void damp_field_in_guards( + amrex::Array4<amrex::Real> const& mf_arr, + const int i, + const int j, + const int k, + const int icomp, + const int dir, + const int n_domain, + const int tb_smallEnd, + const int tb_bigEnd) +{ + using namespace amrex; + + // dir = 0: idx = i (x) + // dir = 1: idx = j (y in 3D, z in 2D/RZ) + // dir = 2: idx = k (z in 3D) + const int idx = ((dir == 0) ? i : ((dir == 1) ? j : k)); + + const int n_tile = tb_bigEnd; + + if (tb_smallEnd < 0) + { + // Apply damping factor in guards cells below the lower end of the domain + const int n_guard = -tb_smallEnd; + + const amrex::Real cell = static_cast<amrex::Real>(idx + n_guard); + + const amrex::Real phase = MathConst::pi * cell / n_guard; + const amrex::Real sin_phase = std::sin(phase); + const amrex::Real damp_factor = sin_phase * sin_phase; + + mf_arr(i,j,k,icomp) *= damp_factor; + } + else if (n_tile > n_domain) + { + // Apply damping factor in guards cells above the upper end of the domain + const int n_guard = n_tile - n_domain; + + const amrex::Real cell = static_cast<amrex::Real>(n_tile - idx); + + const amrex::Real phase = MathConst::pi * cell / n_guard; + const amrex::Real sin_phase = std::sin(phase); + const amrex::Real damp_factor = sin_phase * sin_phase; + + mf_arr(i,j,k,icomp) *= damp_factor; + } +} + +#endif diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 2629da00a..895db9fd2 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -959,13 +959,12 @@ WarpX::ReadParameters () "psatd.update_with_rho must be equal to 1 for comoving PSATD"); } -# ifdef WARPX_DIM_RZ - if (!Geom(0).isPeriodic(1)) { + constexpr int zdir = AMREX_SPACEDIM - 1; + if (!Geom(0).isPeriodic(zdir)) + { use_damp_fields_in_z_guard = true; } pp_psatd.query("use_damp_fields_in_z_guard", use_damp_fields_in_z_guard); -# endif - } // for slice generation // |