diff options
Diffstat (limited to 'Source/WarpX.cpp')
| -rw-r--r-- | Source/WarpX.cpp | 201 |
1 files changed, 123 insertions, 78 deletions
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 0f57cbee1..1e48d659d 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -193,7 +193,7 @@ IntVect WarpX::filter_npass_each_dir(1); int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; -bool WarpX::do_nodal = false; +short WarpX::grid_type; amrex::IntVect m_rho_nodal_flag; int WarpX::do_similar_dm_pml = 1; @@ -908,9 +908,12 @@ WarpX::ReadParameters () pp_warpx.query("do_dynamic_scheduling", do_dynamic_scheduling); - pp_warpx.query("do_nodal", do_nodal); + // Integer that corresponds to the type of grid used in the simulation + // (collocated, staggered, hybrid) + grid_type = GetAlgorithmInteger(pp_warpx, "grid_type"); + // Use same shape factors in all directions, for gathering - if (do_nodal) galerkin_interpolation = false; + if (grid_type == GridType::Collocated) galerkin_interpolation = false; #ifdef WARPX_DIM_RZ // Only needs to be set with WARPX_DIM_RZ, otherwise defaults to 1 @@ -919,13 +922,47 @@ WarpX::ReadParameters () "The number of azimuthal modes (n_rz_azimuthal_modes) must be at least 1"); #endif - // If true, the current is deposited on a nodal grid and centered onto a staggered grid. - // Setting warpx.do_current_centering = 1 makes sense only if warpx.do_nodal = 0. Instead, - // if warpx.do_nodal = 1, Maxwell's equations are solved on a nodal grid and the current - // should not be centered onto a staggered grid. - if (WarpX::do_nodal == 0) + // Set default parameters with hybrid grid (parsed later below) + if (grid_type == GridType::Hybrid) + { + // Finite-order centering of fields (staggered to nodal) + // Default field gathering algorithm will be set below + field_centering_nox = 8; + field_centering_noy = 8; + field_centering_noz = 8; + // Finite-order centering of currents (nodal to staggered) + do_current_centering = true; + current_centering_nox = 8; + current_centering_noy = 8; + current_centering_noz = 8; + } + + // If true, the current is deposited on a nodal grid and centered onto + // a staggered grid. Setting warpx.do_current_centering=1 makes sense + // only if warpx.grid_type=hybrid. Instead, if warpx.grid_type=nodal or + // warpx.grid_type=staggered, Maxwell's equations are solved either on a + // collocated grid or on a staggered grid without current centering. + pp_warpx.query("do_current_centering", do_current_centering); + if (do_current_centering) { - pp_warpx.query("do_current_centering", do_current_centering); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + grid_type == GridType::Hybrid, + "warpx.do_current_centering=1 can be used only with warpx.grid_type=hybrid"); + + utils::parser::queryWithParser( + pp_warpx, "current_centering_nox", current_centering_nox); + utils::parser::queryWithParser( + pp_warpx, "current_centering_noy", current_centering_noy); + utils::parser::queryWithParser( + pp_warpx, "current_centering_noz", current_centering_noz); + + AllocateCenteringCoefficients(device_current_centering_stencil_coeffs_x, + device_current_centering_stencil_coeffs_y, + device_current_centering_stencil_coeffs_z, + current_centering_nox, + current_centering_noy, + current_centering_noz, + grid_type); } WARPX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -957,9 +994,9 @@ WarpX::ReadParameters () } if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { - // Force do_nodal=true (that is, not staggered) and - // use same shape factors in all directions, for gathering - do_nodal = true; + // Force grid_type=collocated (neither staggered nor hybrid) + // and use same shape factors in all directions for gathering + grid_type = GridType::Collocated; galerkin_interpolation = false; } #endif @@ -998,12 +1035,40 @@ WarpX::ReadParameters () "Vay deposition not implemented with multi-J algorithm"); } + // Query algo.field_gathering from input, set field_gathering_algo to + // "default" if not found (default defined in Utils/WarpXAlgorithmSelection.cpp) field_gathering_algo = GetAlgorithmInteger(pp_algo, "field_gathering"); - if (field_gathering_algo == GatheringAlgo::MomentumConserving) { - // Use same shape factors in all directions, for gathering - galerkin_interpolation = false; + + // Set default field gathering algorithm for hybrid grids (momentum-conserving) + std::string tmp_algo; + // - algo.field_gathering not found in the input + // - field_gathering_algo set to "default" above + // (default defined in Utils/WarpXAlgorithmSelection.cpp) + // - reset default value here for hybrid grids + if (pp_algo.query("field_gathering", tmp_algo) == false) + { + if (grid_type == GridType::Hybrid) + { + field_gathering_algo = GatheringAlgo::MomentumConserving; + } + } + // - algo.field_gathering found in the input + // - field_gathering_algo read above and set to user-defined value + else + { + if (grid_type == GridType::Hybrid) + { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + field_gathering_algo == GatheringAlgo::MomentumConserving, + "Hybrid grid (warpx.grid_type=hybrid) should be used only with " + "momentum-conserving field gathering algorithm " + "(algo.field_gathering=momentum-conserving)"); + } } + // Use same shape factors in all directions, for gathering + if (field_gathering_algo == GatheringAlgo::MomentumConserving) galerkin_interpolation = false; + em_solver_medium = GetAlgorithmInteger(pp_algo, "em_solver_medium"); if (em_solver_medium == MediumForEM::Macroscopic ) { macroscopic_solver_algo = GetAlgorithmInteger(pp_algo,"macroscopic_sigma_method"); @@ -1106,62 +1171,43 @@ WarpX::ReadParameters () ParmParse pp_interpolation("interpolation"); pp_interpolation.query("galerkin_scheme",galerkin_interpolation); + } - // Read order of finite-order centering of fields (staggered to nodal). - // Read this only if warpx.do_nodal = 0. Instead, if warpx.do_nodal = 1, - // Maxwell's equations are solved on a nodal grid and the electromagnetic - // forces are gathered from a nodal grid, hence the fields do not need to - // be centered onto a nodal grid. - if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving && - WarpX::do_nodal == 0) - { - utils::parser::queryWithParser( - pp_interpolation, "field_centering_nox", field_centering_nox); - utils::parser::queryWithParser( - pp_interpolation, "field_centering_noy", field_centering_noy); - utils::parser::queryWithParser( - pp_interpolation, "field_centering_noz", field_centering_noz); - } + { + ParmParse pp_warpx("warpx"); - // Read order of finite-order centering of currents (nodal to staggered) - if (WarpX::do_current_centering) + // If warpx.grid_type=staggered or warpx.grid_type=hybrid, + // and algo.field_gathering=momentum-conserving, the fields are solved + // on a staggered grid and centered onto a nodal grid for gathering. + // Instead, if warpx.grid_type=collocated, the momentum-conserving and + // energy conserving field gathering algorithms are equivalent (forces + // gathered from the collocated grid) and no fields centering occurs. + if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving && + WarpX::grid_type != GridType::Collocated) { utils::parser::queryWithParser( - pp_interpolation, "current_centering_nox", current_centering_nox); + pp_warpx, "field_centering_nox", field_centering_nox); utils::parser::queryWithParser( - pp_interpolation, "current_centering_noy", current_centering_noy); + pp_warpx, "field_centering_noy", field_centering_noy); utils::parser::queryWithParser( - pp_interpolation, "current_centering_noz", current_centering_noz); - } - - // Finite-order centering is not implemented with mesh refinement - // (note that when WarpX::do_nodal = 1 finite-order centering is not used anyways) - if (maxLevel() > 0 && WarpX::do_nodal == 0) - { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( - field_centering_nox == 2 && field_centering_noy == 2 && field_centering_noz == 2, - "High-order centering of fields (order > 2) is not implemented with mesh refinement"); - } + pp_warpx, "field_centering_noz", field_centering_noz); - if (WarpX::field_gathering_algo == GatheringAlgo::MomentumConserving && - WarpX::do_nodal == 0) - { AllocateCenteringCoefficients(device_field_centering_stencil_coeffs_x, device_field_centering_stencil_coeffs_y, device_field_centering_stencil_coeffs_z, field_centering_nox, field_centering_noy, - field_centering_noz); + field_centering_noz, + grid_type); } - if (WarpX::do_current_centering) + // Finite-order centering is not implemented with mesh refinement + // (note that when warpx.grid_type=collocated, finite-order centering is not used anyways) + if (maxLevel() > 0 && WarpX::grid_type != GridType::Collocated) { - AllocateCenteringCoefficients(device_current_centering_stencil_coeffs_x, - device_current_centering_stencil_coeffs_y, - device_current_centering_stencil_coeffs_z, - current_centering_nox, - current_centering_noy, - current_centering_noz); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + field_centering_nox == 2 && field_centering_noy == 2 && field_centering_noz == 2, + "High-order centering of fields (order > 2) is not implemented with mesh refinement"); } } @@ -1805,7 +1851,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d dx, do_subcycling, WarpX::use_fdtd_nci_corr, - do_nodal, + grid_type, do_moving_window, moving_window_dir, WarpX::nox, @@ -1920,7 +1966,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm G_nodal_flag = amrex::IntVect::TheCellVector(); // Overwrite nodal flags if necessary - if (do_nodal) { + if (grid_type == GridType::Collocated) { Ex_nodal_flag = IntVect::TheNodeVector(); Ey_nodal_flag = IntVect::TheNodeVector(); Ez_nodal_flag = IntVect::TheNodeVector(); @@ -2171,13 +2217,13 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm #endif } // ElectromagneticSolverAlgo::PSATD else { - m_fdtd_solver_fp[lev] = std::make_unique<FiniteDifferenceSolver>(electromagnetic_solver_id, dx, do_nodal); + m_fdtd_solver_fp[lev] = std::make_unique<FiniteDifferenceSolver>(electromagnetic_solver_id, dx, grid_type); } // // The Aux patch (i.e., the full solution) // - if (aux_is_nodal and !do_nodal) + if (aux_is_nodal and grid_type != GridType::Collocated) { // Create aux multifabs on Nodal Box Array BoxArray const nba = amrex::convert(ba,IntVect::TheNodeVector()); @@ -2266,11 +2312,11 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (do_divb_cleaning) { - if (do_nodal) + if (grid_type == GridType::Collocated) { AllocInitMultiFab(G_cp[lev], amrex::convert(cba, IntVect::TheUnitVector()), dm, ncomps, ngG, tag("G_cp"), 0.0_rt); } - else // do_nodal = 0 + else // grid_type=staggered or grid_type=hybrid { AllocInitMultiFab(G_cp[lev], amrex::convert(cba, IntVect::TheZeroVector()), dm, ncomps, ngG, tag("G_cp"), 0.0_rt); } @@ -2314,7 +2360,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm } // ElectromagneticSolverAlgo::PSATD else { m_fdtd_solver_cp[lev] = std::make_unique<FiniteDifferenceSolver>(electromagnetic_solver_id, cdx, - do_nodal); + grid_type); } } @@ -2401,7 +2447,7 @@ void WarpX::AllocLevelSpectralSolverRZ (amrex::Vector<std::unique_ptr<SpectralSo dm, n_rz_azimuthal_modes, noz_fft, - do_nodal, + grid_type, m_v_galilean, dx_vect, solver_dt, @@ -2455,7 +2501,7 @@ void WarpX::AllocLevelSpectralSolver (amrex::Vector<std::unique_ptr<SpectralSolv nox_fft, noy_fft, noz_fft, - do_nodal, + grid_type, m_v_galilean, m_v_comoving, dx_vect, @@ -2557,9 +2603,8 @@ WarpX::ComputeDivB (amrex::MultiFab& divB, int const dcomp, const std::array<const amrex::MultiFab* const, 3>& B, const std::array<amrex::Real,3>& dx) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!do_nodal, - "ComputeDivB not implemented with do_nodal." - "Shouldn't be too hard to make it general with class FiniteDifferenceSolver"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(grid_type != GridType::Collocated, + "ComputeDivB not implemented with warpx.grid_type=Collocated."); Real dxinv = 1._rt/dx[0], dyinv = 1._rt/dx[1], dzinv = 1._rt/dx[2]; @@ -2595,9 +2640,8 @@ WarpX::ComputeDivB (amrex::MultiFab& divB, int const dcomp, const std::array<const amrex::MultiFab* const, 3>& B, const std::array<amrex::Real,3>& dx, IntVect const ngrow) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!do_nodal, - "ComputeDivB not implemented with do_nodal." - "Shouldn't be too hard to make it general with class FiniteDifferenceSolver"); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(grid_type != GridType::Collocated, + "ComputeDivB not implemented with warpx.grid_type=collocated."); Real dxinv = 1._rt/dx[0], dyinv = 1._rt/dx[1], dzinv = 1._rt/dx[2]; @@ -2795,7 +2839,7 @@ WarpX::BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_ma #endif } -amrex::Vector<amrex::Real> WarpX::getFornbergStencilCoefficients(const int n_order, const bool nodal) +amrex::Vector<amrex::Real> WarpX::getFornbergStencilCoefficients(const int n_order, const short a_grid_type) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(n_order % 2 == 0, "n_order must be even"); @@ -2807,8 +2851,8 @@ amrex::Vector<amrex::Real> WarpX::getFornbergStencilCoefficients(const int n_ord // an overflow when evaluated numerically. One way to avoid the overflow is // to calculate the coefficients by recurrence. - // Coefficients for nodal (centered) finite-difference approximation - if (nodal == true) + // Coefficients for collocated (nodal) finite-difference approximation + if (a_grid_type == GridType::Collocated) { // First coefficient coeffs.at(0) = m * 2. / (m+1); @@ -2856,7 +2900,8 @@ void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector<amrex::Real> amrex::Gpu::DeviceVector<amrex::Real>& device_centering_stencil_coeffs_z, const int centering_nox, const int centering_noy, - const int centering_noz) + const int centering_noz, + const short a_grid_type) { // Vectors of Fornberg stencil coefficients amrex::Vector<amrex::Real> Fornberg_stencil_coeffs_x; @@ -2868,9 +2913,9 @@ void WarpX::AllocateCenteringCoefficients (amrex::Gpu::DeviceVector<amrex::Real> amrex::Vector<amrex::Real> host_centering_stencil_coeffs_y; amrex::Vector<amrex::Real> host_centering_stencil_coeffs_z; - Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(centering_nox, false); - Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(centering_noy, false); - Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(centering_noz, false); + Fornberg_stencil_coeffs_x = getFornbergStencilCoefficients(centering_nox, a_grid_type); + Fornberg_stencil_coeffs_y = getFornbergStencilCoefficients(centering_noy, a_grid_type); + Fornberg_stencil_coeffs_z = getFornbergStencilCoefficients(centering_noz, a_grid_type); host_centering_stencil_coeffs_x.resize(centering_nox); host_centering_stencil_coeffs_y.resize(centering_noy); |