diff options
12 files changed, 402 insertions, 112 deletions
diff --git a/Examples/Tests/galilean/analysis.py b/Examples/Tests/galilean/analysis.py index 898ac1435..9fe6bab72 100755 --- a/Examples/Tests/galilean/analysis.py +++ b/Examples/Tests/galilean/analysis.py @@ -28,6 +28,7 @@ filename = sys.argv[1] # Parse some input arguments from output file 'warpx_used_inputs' current_correction = False time_averaging = False +periodic_single_box = False warpx_used_inputs = open('./warpx_used_inputs', 'r').read() if re.search('geometry.dims\s*=\s*2', warpx_used_inputs): dims = '2D' @@ -39,6 +40,8 @@ if re.search('psatd.current_correction\s*=\s*1', warpx_used_inputs): current_correction = True if re.search('psatd.do_time_averaging\s*=\s*1', warpx_used_inputs): time_averaging = True +if re.search('psatd.periodic_single_box_fft\s*=\s*1', warpx_used_inputs): + periodic_single_box = True ds = yt.load(filename) @@ -58,21 +61,31 @@ tol_charge = 1e-9 if dims == '2D': if not current_correction: energy_ref = 35657.41657683263 - if current_correction: + if current_correction and periodic_single_box: energy_ref = 35024.0275199999 + if current_correction and not periodic_single_box: + energy_ref = 35675.25563324745 + tol_energy = 2e-8 + tol_charge = 2e-4 if time_averaging: energy_ref = 26208.04843478073 tol_energy = 1e-6 elif dims == 'RZ': if not current_correction: energy_ref = 191002.6526271543 - if current_correction: + if current_correction and periodic_single_box: energy_ref = 472779.70801323955 + if current_correction and not periodic_single_box: + energy_ref = 511671.4108624746 + tol_charge = 2e-4 elif dims == '3D': if not current_correction: energy_ref = 661285.098907683 - if current_correction: + if current_correction and periodic_single_box: energy_ref = 856783.3007547935 + if current_correction and not periodic_single_box: + energy_ref = 875307.5138913819 + tol_charge = 1e-2 if time_averaging: energy_ref = 14.564631643496 tol_energy = 1e-4 diff --git a/Regression/Checksum/benchmarks_json/galilean_2d_psatd_current_correction.json b/Regression/Checksum/benchmarks_json/galilean_2d_psatd_current_correction.json index ce7e37dd1..8d43e0b51 100644 --- a/Regression/Checksum/benchmarks_json/galilean_2d_psatd_current_correction.json +++ b/Regression/Checksum/benchmarks_json/galilean_2d_psatd_current_correction.json @@ -1,31 +1,31 @@ { "electrons": { - "particle_momentum_x": 1.5446768224262931e-21, + "particle_momentum_x": 1.5536145943675821e-21, "particle_momentum_y": 0.0, - "particle_momentum_z": 1.7807675479103523e-16, - "particle_position_x": 1267465.6959779875, - "particle_position_y": 15724303.998550693, + "particle_momentum_z": 1.7807674983977113e-16, + "particle_position_x": 1267465.7449785147, + "particle_position_y": 15724304.003312223, "particle_weight": 1.6888332018290936e+18 }, "ions": { - "particle_momentum_x": 2.609389259207186e-18, + "particle_momentum_x": 2.6369229634007033e-18, "particle_momentum_y": 0.0, - "particle_momentum_z": 3.269760999127174e-13, - "particle_position_x": 1267465.6476913819, - "particle_position_y": 15724303.99882355, + "particle_momentum_z": 3.2697610879105043e-13, + "particle_position_x": 1267465.8546859026, + "particle_position_y": 15724304.003084136, "particle_weight": 1.6888332018290936e+18 }, "lev=0": { "Bx": 0.0, - "By": 0.002955847949540036, + "By": 0.0030245915907163154, "Bz": 0.0, - "Ex": 893534.7419589404, + "Ex": 914579.058236283, "Ey": 0.0, - "Ez": 65657.88089087335, - "divE": 1689738.9457807848, - "jx": 219.89616041000212, + "Ez": 69297.0731582175, + "divE": 1695716.5213981706, + "jx": 211.85815442883973, "jy": 0.0, - "jz": 4459.160944658457, - "rho": 1.4961265980553603e-05 + "jz": 4474.2194399497175, + "rho": 1.5014190692786408e-05 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_2d_psatd_current_correction_psb.json b/Regression/Checksum/benchmarks_json/galilean_2d_psatd_current_correction_psb.json new file mode 100644 index 000000000..ce7e37dd1 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/galilean_2d_psatd_current_correction_psb.json @@ -0,0 +1,31 @@ +{ + "electrons": { + "particle_momentum_x": 1.5446768224262931e-21, + "particle_momentum_y": 0.0, + "particle_momentum_z": 1.7807675479103523e-16, + "particle_position_x": 1267465.6959779875, + "particle_position_y": 15724303.998550693, + "particle_weight": 1.6888332018290936e+18 + }, + "ions": { + "particle_momentum_x": 2.609389259207186e-18, + "particle_momentum_y": 0.0, + "particle_momentum_z": 3.269760999127174e-13, + "particle_position_x": 1267465.6476913819, + "particle_position_y": 15724303.99882355, + "particle_weight": 1.6888332018290936e+18 + }, + "lev=0": { + "Bx": 0.0, + "By": 0.002955847949540036, + "Bz": 0.0, + "Ex": 893534.7419589404, + "Ey": 0.0, + "Ez": 65657.88089087335, + "divE": 1689738.9457807848, + "jx": 219.89616041000212, + "jy": 0.0, + "jz": 4459.160944658457, + "rho": 1.4961265980553603e-05 + } +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_3d_psatd_current_correction.json b/Regression/Checksum/benchmarks_json/galilean_3d_psatd_current_correction.json index 47f4e9c2f..16814f47b 100644 --- a/Regression/Checksum/benchmarks_json/galilean_3d_psatd_current_correction.json +++ b/Regression/Checksum/benchmarks_json/galilean_3d_psatd_current_correction.json @@ -1,33 +1,33 @@ { "electrons": { - "particle_momentum_x": 7.768538631989314e-22, - "particle_momentum_y": 7.844994340141792e-22, - "particle_momentum_z": 8.903837510927089e-17, - "particle_position_x": 158433.32170594894, - "particle_position_y": 158432.8515695265, - "particle_position_z": 5891662.9548929185, + "particle_momentum_x": 7.799828234013091e-22, + "particle_momentum_y": 7.853620523635122e-22, + "particle_momentum_z": 8.903838009194013e-17, + "particle_position_x": 158433.52653925028, + "particle_position_y": 158432.74997779692, + "particle_position_z": 5891662.962223673, "particle_weight": 2.041377132710917e+18 }, "ions": { - "particle_momentum_x": 1.3137653484757431e-18, - "particle_momentum_y": 1.3110225003256574e-18, - "particle_momentum_z": 1.6348803844352492e-13, - "particle_position_x": 158433.312978349, - "particle_position_y": 158432.84896819026, - "particle_position_z": 5891662.955099393, + "particle_momentum_x": 1.3150842145882957e-18, + "particle_momentum_y": 1.3043330137743231e-18, + "particle_momentum_z": 1.634880568866417e-13, + "particle_position_x": 158433.58054870443, + "particle_position_y": 158432.80156002127, + "particle_position_z": 5891662.961766562, "particle_weight": 2.041377132710917e+18 }, "lev=0": { - "Bx": 0.006449275564525033, - "By": 0.0064783778061885235, - "Bz": 0.0006158841538190647, - "Ex": 1950801.24799352, - "Ey": 1945623.9590479229, - "Ez": 150385.8857169562, - "divE": 6191274.556649665, - "jx": 505.8903778073368, - "jy": 510.11178001328403, - "jz": 16346.473505698516, - "rho": 5.481870772518483e-05 + "Bx": 0.0068073634628550315, + "By": 0.006697549542915319, + "Bz": 0.0005743262065396818, + "Ex": 2018478.8211576312, + "Ey": 2052642.246798101, + "Ez": 149912.8262920749, + "divE": 6615408.753505244, + "jx": 490.81504811208214, + "jy": 497.1861446484984, + "jz": 17469.798807376228, + "rho": 5.862037378889882e-05 } }
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_3d_psatd_current_correction_psb.json b/Regression/Checksum/benchmarks_json/galilean_3d_psatd_current_correction_psb.json new file mode 100644 index 000000000..47f4e9c2f --- /dev/null +++ b/Regression/Checksum/benchmarks_json/galilean_3d_psatd_current_correction_psb.json @@ -0,0 +1,33 @@ +{ + "electrons": { + "particle_momentum_x": 7.768538631989314e-22, + "particle_momentum_y": 7.844994340141792e-22, + "particle_momentum_z": 8.903837510927089e-17, + "particle_position_x": 158433.32170594894, + "particle_position_y": 158432.8515695265, + "particle_position_z": 5891662.9548929185, + "particle_weight": 2.041377132710917e+18 + }, + "ions": { + "particle_momentum_x": 1.3137653484757431e-18, + "particle_momentum_y": 1.3110225003256574e-18, + "particle_momentum_z": 1.6348803844352492e-13, + "particle_position_x": 158433.312978349, + "particle_position_y": 158432.84896819026, + "particle_position_z": 5891662.955099393, + "particle_weight": 2.041377132710917e+18 + }, + "lev=0": { + "Bx": 0.006449275564525033, + "By": 0.0064783778061885235, + "Bz": 0.0006158841538190647, + "Ex": 1950801.24799352, + "Ey": 1945623.9590479229, + "Ez": 150385.8857169562, + "divE": 6191274.556649665, + "jx": 505.8903778073368, + "jy": 510.11178001328403, + "jz": 16346.473505698516, + "rho": 5.481870772518483e-05 + } +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json b/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json index 73cd09635..3dd8810fa 100644 --- a/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json +++ b/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction.json @@ -1,30 +1,30 @@ { "electrons": { - "particle_momentum_x": 7.005566900701696e-22, - "particle_momentum_y": 2.741432867448382e-22, - "particle_momentum_z": 8.903838637330695e-17, - "particle_position_x": 633733.8772916717, - "particle_position_y": 7862151.998406768, - "particle_theta": 51362.087114626614, + "particle_momentum_x": 7.040321706649852e-22, + "particle_momentum_y": 2.6350019627856963e-22, + "particle_momentum_z": 8.903841205118235e-17, + "particle_position_x": 633733.2382222114, + "particle_position_y": 7862152.022952429, + "particle_theta": 51150.20809917081, "particle_weight": 1.0261080645329302e+20 }, "ions": { - "particle_momentum_x": 1.3125873166811053e-18, - "particle_momentum_y": 2.7371989851752584e-22, - "particle_momentum_z": 1.6348804524259806e-13, - "particle_position_x": 633733.7424153683, - "particle_position_y": 7862151.997142274, - "particle_theta": 51470.93289837983, + "particle_momentum_x": 1.309068869366174e-18, + "particle_momentum_y": 2.6235118222356472e-22, + "particle_momentum_z": 1.6348804184310653e-13, + "particle_position_x": 633733.1816905431, + "particle_position_y": 7862151.995500945, + "particle_theta": 51448.117157711284, "particle_weight": 1.0261080645329302e+20 }, "lev=0": { - "By": 0.0017275299193154866, - "Ex": 520734.8193541992, - "Ey": 130008.65587120061, - "Ez": 36540.80980315481, - "divE": 1228895.3465882093, - "jx": 103.9586429924576, - "jz": 3245.258754763504, - "rho": 1.088087020096811e-05 + "By": 0.0017182252289440264, + "Ex": 518347.23478549306, + "Ey": 137106.07262027927, + "Ez": 35918.446653926556, + "divE": 1227626.4690642208, + "jx": 159.1245447865701, + "jz": 49252.95943703725, + "rho": 1.0869618748586227e-05 } -} +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction_psb.json b/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction_psb.json new file mode 100644 index 000000000..73cd09635 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/galilean_rz_psatd_current_correction_psb.json @@ -0,0 +1,30 @@ +{ + "electrons": { + "particle_momentum_x": 7.005566900701696e-22, + "particle_momentum_y": 2.741432867448382e-22, + "particle_momentum_z": 8.903838637330695e-17, + "particle_position_x": 633733.8772916717, + "particle_position_y": 7862151.998406768, + "particle_theta": 51362.087114626614, + "particle_weight": 1.0261080645329302e+20 + }, + "ions": { + "particle_momentum_x": 1.3125873166811053e-18, + "particle_momentum_y": 2.7371989851752584e-22, + "particle_momentum_z": 1.6348804524259806e-13, + "particle_position_x": 633733.7424153683, + "particle_position_y": 7862151.997142274, + "particle_theta": 51470.93289837983, + "particle_weight": 1.0261080645329302e+20 + }, + "lev=0": { + "By": 0.0017275299193154866, + "Ex": 520734.8193541992, + "Ey": 130008.65587120061, + "Ez": 36540.80980315481, + "divE": 1228895.3465882093, + "jx": 103.9586429924576, + "jz": 3245.258754763504, + "rho": 1.088087020096811e-05 + } +} diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 6263db152..8d611da80 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2372,7 +2372,7 @@ compareParticles = 1 particleTypes = electrons ions analysisRoutine = Examples/Tests/galilean/analysis.py -[galilean_2d_psatd_current_correction] +[galilean_2d_psatd_current_correction_psb] buildDir = . inputFile = Examples/Tests/galilean/inputs_2d runtime_params = psatd.periodic_single_box_fft=1 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE @@ -2390,6 +2390,24 @@ compareParticles = 1 particleTypes = electrons ions analysisRoutine = Examples/Tests/galilean/analysis.py +[galilean_2d_psatd_current_correction] +buildDir = . +inputFile = Examples/Tests/galilean/inputs_2d +runtime_params = psatd.periodic_single_box_fft=0 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE amr.max_grid_size=64 amr.blocking_factor=64 +dim = 2 +addToCompileString = USE_PSATD=TRUE +cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons ions +analysisRoutine = Examples/Tests/galilean/analysis.py + [galilean_2d_psatd_hybrid] buildDir = . inputFile = Examples/Tests/galilean/inputs_2d_hybrid @@ -2444,7 +2462,7 @@ compareParticles = 1 particleTypes = electrons ions analysisRoutine = Examples/Tests/galilean/analysis.py -[galilean_rz_psatd_current_correction] +[galilean_rz_psatd_current_correction_psb] buildDir = . inputFile = Examples/Tests/galilean/inputs_rz runtime_params = psatd.periodic_single_box_fft=1 psatd.current_correction=1 electrons.random_theta=0 ions.random_theta=0 @@ -2462,6 +2480,24 @@ compareParticles = 1 particleTypes = electrons ions analysisRoutine = Examples/Tests/galilean/analysis.py +[galilean_rz_psatd_current_correction] +buildDir = . +inputFile = Examples/Tests/galilean/inputs_rz +runtime_params = psatd.periodic_single_box_fft=0 psatd.current_correction=1 electrons.random_theta=0 ions.random_theta=0 amr.max_grid_size=32 amr.blocking_factor=32 +dim = 2 +addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack +cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_PSATD=ON +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons ions +analysisRoutine = Examples/Tests/galilean/analysis.py + [galilean_3d_psatd] buildDir = . inputFile = Examples/Tests/galilean/inputs_3d @@ -2480,7 +2516,7 @@ compareParticles = 1 particleTypes = electrons ions analysisRoutine = Examples/Tests/galilean/analysis.py -[galilean_3d_psatd_current_correction] +[galilean_3d_psatd_current_correction_psb] buildDir = . inputFile = Examples/Tests/galilean/inputs_3d runtime_params = warpx.numprocs=1 1 1 psatd.periodic_single_box_fft=1 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE @@ -2498,6 +2534,24 @@ compareParticles = 1 particleTypes = electrons ions analysisRoutine = Examples/Tests/galilean/analysis.py +[galilean_3d_psatd_current_correction] +buildDir = . +inputFile = Examples/Tests/galilean/inputs_3d +runtime_params = warpx.numprocs=1 1 2 psatd.periodic_single_box_fft=0 psatd.update_with_rho=0 psatd.current_correction=1 diag1.fields_to_plot=Ex Ey Ez Bx By Bz jx jy jz rho divE +dim = 3 +addToCompileString = USE_PSATD=TRUE +cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PSATD=ON +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons ions +analysisRoutine = Examples/Tests/galilean/analysis.py + [averaged_galilean_2d_psatd] buildDir = . inputFile = Examples/Tests/averaged_galilean/inputs_avg_2d diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index b7a8c64f1..cc71e78d5 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -398,19 +398,9 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("afterdeposition"); - // Synchronize J and rho: filter, exchange boundary, interpolate across levels. - // With Vay current deposition, the current deposited at this point is not yet - // the actual current J. This is computed later in WarpX::PushPSATD, by calling - // WarpX::PSATDVayDeposition. The function SyncCurrent is called after that, - // instead of here, so that we synchronize the correct current. - // With current centering, the nodal current is deposited in 'current_fp_nodal': - // SyncCurrent stores the result of its centering into 'current_fp' and then - // performs both filtering, if used, and exchange of guard cells. - if (WarpX::current_deposition_algo != CurrentDepositionAlgo::Vay) - { - SyncCurrent(current_fp, current_cp); - } - SyncRho(); + // Synchronize J and rho: + // filter (if used), exchange guard cells, interpolate across MR levels + SyncCurrentAndRho(); // At this point, J is up-to-date inside the domain, and E and B are // up-to-date including enough guard cells for first step of the field @@ -495,6 +485,46 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("afterEsolve"); } +void WarpX::SyncCurrentAndRho () +{ + if (maxwell_solver_id == MaxwellSolverAlgo::PSATD) + { + if (fft_periodic_single_box) + { + // With periodic single box, synchronize J and rho here, + // even with current correction or Vay deposition + if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + // TODO Replace current_cp with current_cp_vay once Vay deposition is implemented with MR + SyncCurrent(current_fp_vay, current_cp); + SyncRho(); + } + else + { + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + } + else // no periodic single box + { + // Without periodic single box, synchronize J and rho here, + // except with current correction or Vay deposition: + // in these cases, synchronize later (in WarpX::PushPSATD) + if (current_correction == false && + current_deposition_algo != CurrentDepositionAlgo::Vay) + { + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + } + } + else // FDTD + { + SyncCurrent(current_fp, current_cp); + SyncRho(); + } +} + void WarpX::OneStep_multiJ (const amrex::Real cur_time) { diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6d604bc15..7499ee140 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -270,7 +270,8 @@ WarpX::PSATDBackwardTransformG () void WarpX::PSATDForwardTransformJ ( const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp, - const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp) + const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp, + const bool apply_kspace_filter) { SpectralFieldIndex Idx; int idx_jx, idx_jy, idx_jz; @@ -299,7 +300,7 @@ void WarpX::PSATDForwardTransformJ ( #ifdef WARPX_DIM_RZ // Apply filter in k space if needed - if (WarpX::use_kspace_filter) + if (use_kspace_filter && apply_kspace_filter) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -311,6 +312,8 @@ void WarpX::PSATDForwardTransformJ ( } } } +#else + amrex::ignore_unused(apply_kspace_filter); #endif } @@ -349,8 +352,10 @@ void WarpX::PSATDBackwardTransformJ ( void WarpX::PSATDForwardTransformRho ( const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_fp, const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_cp, - const int icomp, const int dcomp) + const int icomp, const int dcomp, const bool apply_kspace_filter) { + if (charge_fp[0] == nullptr) return; + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; // Select index in k space @@ -368,7 +373,7 @@ void WarpX::PSATDForwardTransformRho ( #ifdef WARPX_DIM_RZ // Apply filter in k space if needed - if (WarpX::use_kspace_filter) + if (use_kspace_filter && apply_kspace_filter) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -380,6 +385,8 @@ void WarpX::PSATDForwardTransformRho ( } } } +#else + amrex::ignore_unused(apply_kspace_filter); #endif } @@ -635,46 +642,114 @@ WarpX::PushPSATD () "PushFieldsEM: PSATD solver selected but not built")); #else - PSATDForwardTransformEB(Efield_fp, Bfield_fp, Efield_cp, Bfield_cp); - - amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp = - (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) ? current_fp_vay : current_fp; + if (fft_periodic_single_box) + { + if (current_correction) + { + // FFT of J and rho + PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new - PSATDForwardTransformJ(J_fp, current_cp); + // Correct J in k-space + PSATDCurrentCorrection(); - // Do rho FFTs only if needed - if (WarpX::update_with_rho || WarpX::current_correction || WarpX::do_dive_cleaning) - { - PSATDForwardTransformRho(rho_fp, rho_cp, 0,0); // rho old - PSATDForwardTransformRho(rho_fp, rho_cp, 1,1); // rho new + // Inverse FFT of J + PSATDBackwardTransformJ(current_fp, current_cp); + } + else if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + // FFT of D and rho (if used) + // TODO Replace current_cp with current_cp_vay once Vay deposition is implemented with MR + PSATDForwardTransformJ(current_fp_vay, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new + + // Compute J from D in k-space + PSATDVayDeposition(); + + // Inverse FFT of J, subtract cumulative sums of D + PSATDBackwardTransformJ(current_fp, current_cp); + // TODO Cumulative sums need to be fixed with periodic single box + PSATDSubtractCurrentPartialSumsAvg(); + + // FFT of J after subtraction of cumulative sums + PSATDForwardTransformJ(current_fp, current_cp); + } + else // no current correction, no Vay deposition + { + // FFT of J and rho (if used) + PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new + } } - - // Correct the current in Fourier space so that the continuity equation is satisfied, and - // transform back to real space so that the current correction is reflected in the diagnostics - if (WarpX::current_correction) + else // no periodic single box { - PSATDCurrentCorrection(); - PSATDBackwardTransformJ(current_fp, current_cp); - } + if (current_correction) + { + // FFT of J and rho +#ifdef WARPX_DIM_RZ + // In RZ geometry, do not apply filtering here, since it is + // applied in the subsequent calls to these functions (below) + const bool apply_kspace_filter = false; + PSATDForwardTransformJ(current_fp, current_cp, apply_kspace_filter); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0, apply_kspace_filter); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1, apply_kspace_filter); // rho new +#else + PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new +#endif - // Compute the current in Fourier space according to the Vay deposition scheme, and - // transform back to real space so that the Vay deposition is reflected in the diagnostics - if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) - { - PSATDVayDeposition(); - PSATDBackwardTransformJ(current_fp, current_cp); - PSATDSubtractCurrentPartialSumsAvg(); - SyncCurrent(current_fp, current_cp); + // Correct J in k-space + PSATDCurrentCorrection(); + + // Inverse FFT of J + PSATDBackwardTransformJ(current_fp, current_cp); + + // Synchronize J and rho + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + else if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + // FFT of D + PSATDForwardTransformJ(current_fp_vay, current_cp); + + // Compute J from D in k-space + PSATDVayDeposition(); + + // Inverse FFT of J, subtract cumulative sums of D + PSATDBackwardTransformJ(current_fp, current_cp); + PSATDSubtractCurrentPartialSumsAvg(); + + // Synchronize J and rho (if used) + SyncCurrent(current_fp, current_cp); + SyncRho(); + } + + // FFT of J and rho (if used) PSATDForwardTransformJ(current_fp, current_cp); + PSATDForwardTransformRho(rho_fp, rho_cp, 0, 0); // rho old + PSATDForwardTransformRho(rho_fp, rho_cp, 1, 1); // rho new } + // FFT of E and B + PSATDForwardTransformEB(Efield_fp, Bfield_fp, Efield_cp, Bfield_cp); + #ifdef WARPX_DIM_RZ if (pml_rz[0]) pml_rz[0]->PushPSATD(0); #endif + // FFT of F and G if (WarpX::do_dive_cleaning) PSATDForwardTransformF(); if (WarpX::do_divb_cleaning) PSATDForwardTransformG(); + + // Update E, B, F, and G in k-space PSATDPushSpectralFields(); + + // Inverse FFT of E, B, F, and G PSATDBackwardTransformEB(Efield_fp, Bfield_fp, Efield_cp, Bfield_cp); if (WarpX::fft_do_time_averaging) PSATDBackwardTransformEBavg(Efield_avg_fp, Bfield_avg_fp, Efield_avg_cp, Bfield_avg_cp); diff --git a/Source/WarpX.H b/Source/WarpX.H index c16f30408..7a051e30d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -647,6 +647,13 @@ public: void FillBoundaryAux (int lev, amrex::IntVect ng); /** + * \brief Synchronize J and rho: + * filter (if used), exchange guard cells, interpolate across MR levels. + * Contains separate calls to WarpX::SyncCurrent and WarpX::SyncRho. + */ + void SyncCurrentAndRho (); + + /** * \brief Apply filter and sum guard cells across MR levels. * If current centering is used, center the current from a nodal grid * to a staggered grid. For each MR level beyond level 0, interpolate @@ -1506,10 +1513,13 @@ private: * storing the fine patch current to be transformed * \param J_cp Vector of three-dimensional arrays (for each level) * storing the coarse patch current to be transformed + * \param[in] apply_kspace_filter Control whether to apply filtering + * (only used in RZ geometry to avoid double filtering) */ void PSATDForwardTransformJ ( const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_fp, - const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp); + const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>,3>>& J_cp, + const bool apply_kspace_filter=true); /** * \brief Backward FFT of J on all mesh refinement levels @@ -1531,11 +1541,13 @@ private: * \param charge_cp Vector (for each level) storing the coarse patch charge to be transformed * \param[in] icomp index of fourth component (0 for rho_old, 1 for rho_new) * \param[in] dcomp index of spectral component (0 for rho_old, 1 for rho_new) + * \param[in] apply_kspace_filter Control whether to apply filtering + * (only used in RZ geometry to avoid double filtering) */ void PSATDForwardTransformRho ( const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_fp, const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& charge_cp, - const int icomp, const int dcomp); + const int icomp, const int dcomp, const bool apply_kspace_filter=true); /** * \brief Copy rho_new to rho_old in spectral space diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 12747b7b9..86961dcaf 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1194,6 +1194,13 @@ WarpX::ReadParameters () "Option algo.current_deposition=vay must be used with psatd.periodic_single_box_fft=0."); } + if (current_deposition_algo == CurrentDepositionAlgo::Vay) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + current_correction == false, + "Options algo.current_deposition=vay and psatd.current_correction=1 cannot be combined together."); + } + // Auxiliary: boosted_frame = true if warpx.gamma_boost is set in the inputs amrex::ParmParse pp_warpx("warpx"); const bool boosted_frame = pp_warpx.query("gamma_boost", gamma_boost); @@ -1328,10 +1335,15 @@ WarpX::ReadParameters () } } - // Fill guard cells with backward FFTs if Vay current deposition is used - if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Vay) + // Without periodic single box, fill guard cells with backward FFTs, + // with current correction or Vay deposition + if (fft_periodic_single_box == false) { - WarpX::m_fill_guards_current = amrex::IntVect(1); + if (current_correction || + current_deposition_algo == CurrentDepositionAlgo::Vay) + { + WarpX::m_fill_guards_current = amrex::IntVect(1); + } } } |