diff options
Diffstat (limited to 'Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp')
| -rw-r--r-- | Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp | 338 |
1 files changed, 338 insertions, 0 deletions
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp new file mode 100644 index 000000000..9e04045d0 --- /dev/null +++ b/Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp @@ -0,0 +1,338 @@ +/* Copyright 2023 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Roelof Groenewald (TAE Technologies) + * + * License: BSD-3-Clause-LBNL + */ + +#include "HybridPICModel.H" + +using namespace amrex; + +HybridPICModel::HybridPICModel ( int nlevs_max ) +{ + ReadParameters(); + AllocateMFs(nlevs_max); +} + +void HybridPICModel::ReadParameters () +{ + ParmParse pp_hybrid("hybrid_pic_model"); + + // The B-field update is subcycled to improve stability - the number + // of sub steps can be specified by the user (defaults to 50). + utils::parser::queryWithParser(pp_hybrid, "substeps", m_substeps); + + // The hybrid model requires an electron temperature, reference density + // and exponent to be given. These values will be used to calculate the + // electron pressure according to p = n0 * Te * (n/n0)^gamma + utils::parser::queryWithParser(pp_hybrid, "gamma", m_gamma); + if (!utils::parser::queryWithParser(pp_hybrid, "elec_temp", m_elec_temp)) { + Abort("hybrid_pic_model.elec_temp must be specified when using the hybrid solver"); + } + bool n0_ref_given = utils::parser::queryWithParser(pp_hybrid, "n0_ref", m_n0_ref); + if (m_gamma != 1.0 && !n0_ref_given) { + Abort("hybrid_pic_model.n0_ref should be specified if hybrid_pic_model.gamma != 1"); + } + + pp_hybrid.query("plasma_resistivity(rho)", m_eta_expression); + utils::parser::queryWithParser(pp_hybrid, "n_floor", m_n_floor); + + // convert electron temperature from eV to J + m_elec_temp *= PhysConst::q_e; +} + +void HybridPICModel::AllocateMFs (int nlevs_max) +{ + electron_pressure_fp.resize(nlevs_max); + rho_fp_temp.resize(nlevs_max); + current_fp_temp.resize(nlevs_max); + current_fp_ampere.resize(nlevs_max); +} + +void HybridPICModel::AllocateLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm, + const int ncomps, const IntVect& ngJ, const IntVect& ngRho, + const IntVect& jx_nodal_flag, + const IntVect& jy_nodal_flag, + const IntVect& jz_nodal_flag, + const IntVect& rho_nodal_flag) +{ + // set human-readable tag for each MultiFab + auto const tag = [lev]( std::string tagname ) { + tagname.append("[l=").append(std::to_string(lev)).append("]"); + return tagname; + }; + + auto & warpx = WarpX::GetInstance(); + + // The "electron_pressure_fp" multifab stores the electron pressure calculated + // from the specified equation of state. + // The "rho_fp_temp" multifab is used to store the ion charge density + // interpolated or extrapolated to appropriate timesteps. + // The "current_fp_temp" multifab is used to store the ion current density + // interpolated or extrapolated to appropriate timesteps. + // The "current_fp_ampere" multifab stores the total current calculated as + // the curl of B. + warpx.AllocInitMultiFab(electron_pressure_fp[lev], amrex::convert(ba, rho_nodal_flag), + dm, ncomps, ngRho, tag("electron_pressure_fp"), 0.0_rt); + + warpx.AllocInitMultiFab(rho_fp_temp[lev], amrex::convert(ba, rho_nodal_flag), + dm, ncomps, ngRho, tag("rho_fp_temp"), 0.0_rt); + + warpx.AllocInitMultiFab(current_fp_temp[lev][0], amrex::convert(ba, jx_nodal_flag), + dm, ncomps, ngJ, tag("current_fp_temp[x]"), 0.0_rt); + warpx.AllocInitMultiFab(current_fp_temp[lev][1], amrex::convert(ba, jy_nodal_flag), + dm, ncomps, ngJ, tag("current_fp_temp[y]"), 0.0_rt); + warpx.AllocInitMultiFab(current_fp_temp[lev][2], amrex::convert(ba, jz_nodal_flag), + dm, ncomps, ngJ, tag("current_fp_temp[z]"), 0.0_rt); + + warpx.AllocInitMultiFab(current_fp_ampere[lev][0], amrex::convert(ba, jx_nodal_flag), + dm, ncomps, ngJ, tag("current_fp_ampere[x]"), 0.0_rt); + warpx.AllocInitMultiFab(current_fp_ampere[lev][1], amrex::convert(ba, jy_nodal_flag), + dm, ncomps, ngJ, tag("current_fp_ampere[y]"), 0.0_rt); + warpx.AllocInitMultiFab(current_fp_ampere[lev][2], amrex::convert(ba, jz_nodal_flag), + dm, ncomps, ngJ, tag("current_fp_ampere[z]"), 0.0_rt); +} + +void HybridPICModel::ClearLevel (int lev) +{ + electron_pressure_fp[lev].reset(); + rho_fp_temp[lev].reset(); + for (int i = 0; i < 3; ++i) { + current_fp_temp[lev][i].reset(); + current_fp_ampere[lev][i].reset(); + } +} + +void HybridPICModel::InitData () +{ + m_resistivity_parser = std::make_unique<amrex::Parser>( + utils::parser::makeParser(m_eta_expression, {"rho"})); + m_eta = m_resistivity_parser->compile<1>(); + + auto & warpx = WarpX::GetInstance(); + + // Get the grid staggering of the fields involved in calculating E + amrex::IntVect Jx_stag = warpx.getcurrent_fp(0,0).ixType().toIntVect(); + amrex::IntVect Jy_stag = warpx.getcurrent_fp(0,1).ixType().toIntVect(); + amrex::IntVect Jz_stag = warpx.getcurrent_fp(0,2).ixType().toIntVect(); + amrex::IntVect Bx_stag = warpx.getBfield_fp(0,0).ixType().toIntVect(); + amrex::IntVect By_stag = warpx.getBfield_fp(0,1).ixType().toIntVect(); + amrex::IntVect Bz_stag = warpx.getBfield_fp(0,2).ixType().toIntVect(); + amrex::IntVect Ex_stag = warpx.getEfield_fp(0,0).ixType().toIntVect(); + amrex::IntVect Ey_stag = warpx.getEfield_fp(0,1).ixType().toIntVect(); + amrex::IntVect Ez_stag = warpx.getEfield_fp(0,2).ixType().toIntVect(); + + // copy data to device + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + Jx_IndexType[idim] = Jx_stag[idim]; + Jy_IndexType[idim] = Jy_stag[idim]; + Jz_IndexType[idim] = Jz_stag[idim]; + Bx_IndexType[idim] = Bx_stag[idim]; + By_IndexType[idim] = By_stag[idim]; + Bz_IndexType[idim] = Bz_stag[idim]; + Ex_IndexType[idim] = Ex_stag[idim]; + Ey_IndexType[idim] = Ey_stag[idim]; + Ez_IndexType[idim] = Ez_stag[idim]; + } + + // Below we set all the unused dimensions to have nodal values for J, B & E + // since these values will be interpolated onto a nodal grid - if this is + // not done the Interp function returns nonsense values. +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_1D_Z) + Jx_IndexType[2] = 1; + Jy_IndexType[2] = 1; + Jz_IndexType[2] = 1; + Bx_IndexType[2] = 1; + By_IndexType[2] = 1; + Bz_IndexType[2] = 1; + Ex_IndexType[2] = 1; + Ey_IndexType[2] = 1; + Ez_IndexType[2] = 1; +#endif +#if defined(WARPX_DIM_1D_Z) + Jx_IndexType[1] = 1; + Jy_IndexType[1] = 1; + Jz_IndexType[1] = 1; + Bx_IndexType[1] = 1; + By_IndexType[1] = 1; + Bz_IndexType[1] = 1; + Ex_IndexType[1] = 1; + Ey_IndexType[1] = 1; + Ez_IndexType[1] = 1; +#endif +} + +void HybridPICModel::CalculateCurrentAmpere ( + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield, + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths) +{ + auto& warpx = WarpX::GetInstance(); + for (int lev = 0; lev <= warpx.finestLevel(); ++lev) + { + CalculateCurrentAmpere(Bfield[lev], edge_lengths[lev], lev); + } +} + +void HybridPICModel::CalculateCurrentAmpere ( + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths, + const int lev) +{ + WARPX_PROFILE("WarpX::CalculateCurrentAmpere()"); + + auto& warpx = WarpX::GetInstance(); + warpx.get_pointer_fdtd_solver_fp(lev)->CalculateCurrentAmpere( + current_fp_ampere[lev], Bfield, edge_lengths, lev + ); + + // we shouldn't apply the boundary condition to J since J = J_i - J_e but + // the boundary correction was already applied to J_i and the B-field + // boundary ensures that J itself complies with the boundary conditions, right? + // ApplyJfieldBoundary(lev, Jfield[0].get(), Jfield[1].get(), Jfield[2].get()); + for (int i=0; i<3; i++) current_fp_ampere[lev][i]->FillBoundary(warpx.Geom(lev).periodicity()); +} + +void HybridPICModel::HybridPICSolveE ( + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> & Efield, + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Jfield, + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& Bfield, + amrex::Vector<std::unique_ptr<amrex::MultiFab>> const& rhofield, + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3>> const& edge_lengths, + DtType a_dt_type) +{ + auto& warpx = WarpX::GetInstance(); + for (int lev = 0; lev <= warpx.finestLevel(); ++lev) + { + HybridPICSolveE( + Efield[lev], Jfield[lev], Bfield[lev], rhofield[lev], + edge_lengths[lev], lev, a_dt_type + ); + } +} + +void HybridPICModel::HybridPICSolveE ( + std::array< std::unique_ptr<amrex::MultiFab>, 3> & Efield, + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield, + std::unique_ptr<amrex::MultiFab> const& rhofield, + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths, + const int lev, DtType a_dt_type) +{ + WARPX_PROFILE("WarpX::HybridPICSolveE()"); + + HybridPICSolveE( + Efield, Jfield, Bfield, rhofield, edge_lengths, lev, + PatchType::fine, a_dt_type + ); + if (lev > 0) + { + amrex::Abort(Utils::TextMsg::Err( + "HybridPICSolveE: Only one level implemented for hybrid-PIC solver.")); + } +} + +void HybridPICModel::HybridPICSolveE ( + std::array< std::unique_ptr<amrex::MultiFab>, 3> & Efield, + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Jfield, + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& Bfield, + std::unique_ptr<amrex::MultiFab> const& rhofield, + std::array< std::unique_ptr<amrex::MultiFab>, 3> const& edge_lengths, + const int lev, PatchType patch_type, DtType a_dt_type) +{ + auto& warpx = WarpX::GetInstance(); + + // Solve E field in regular cells + // The first half step uses t=n quantities, the second half t=n+1/2 + // quantities and the full step uses t=n+1 quantities + if (a_dt_type == DtType::FirstHalf) { + warpx.get_pointer_fdtd_solver_fp(lev)->HybridPICSolveE( + Efield, current_fp_ampere[lev], + current_fp_temp[lev], Bfield, + rho_fp_temp[lev], + electron_pressure_fp[lev], + edge_lengths, lev, this, a_dt_type + ); + } + else if (a_dt_type == DtType::SecondHalf) { + warpx.get_pointer_fdtd_solver_fp(lev)->HybridPICSolveE( + Efield, current_fp_ampere[lev], + Jfield, Bfield, + rho_fp_temp[lev], + electron_pressure_fp[lev], + edge_lengths, lev, this, a_dt_type + ); + } + else { + warpx.get_pointer_fdtd_solver_fp(lev)->HybridPICSolveE( + Efield, current_fp_ampere[lev], + current_fp_temp[lev], Bfield, + rhofield, + electron_pressure_fp[lev], + edge_lengths, lev, this, a_dt_type + ); + } + + warpx.ApplyEfieldBoundary(lev, patch_type); +} + +void HybridPICModel::CalculateElectronPressure(DtType a_dt_type) +{ + auto& warpx = WarpX::GetInstance(); + for (int lev = 0; lev <= warpx.finestLevel(); ++lev) + { + CalculateElectronPressure(lev, a_dt_type); + } +} + +void HybridPICModel::CalculateElectronPressure(const int lev, DtType a_dt_type) +{ + WARPX_PROFILE("WarpX::CalculateElectronPressure()"); + + auto& warpx = WarpX::GetInstance(); + // The full step uses rho^{n+1}, otherwise use the old or averaged + // charge density. + if (a_dt_type == DtType::Full) { + FillElectronPressureMF( + electron_pressure_fp[lev], warpx.get_pointer_rho_fp(lev) + ); + } else { + FillElectronPressureMF( + electron_pressure_fp[lev], rho_fp_temp[lev].get() + ); + } + warpx.ApplyElectronPressureBoundary(lev, PatchType::fine); + electron_pressure_fp[lev]->FillBoundary(warpx.Geom(lev).periodicity()); +} + + +void HybridPICModel::FillElectronPressureMF ( + std::unique_ptr<amrex::MultiFab> const& Pe_field, + amrex::MultiFab* const& rho_field ) +{ + const auto n0_ref = m_n0_ref; + const auto elec_temp = m_elec_temp; + const auto gamma = m_gamma; + + // Loop through the grids, and over the tiles within each grid +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(*Pe_field, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + // Extract field data for this grid/tile + Array4<Real const> const& rho = rho_field->const_array(mfi); + Array4<Real> const& Pe = Pe_field->array(mfi); + + // Extract tileboxes for which to loop + const Box& tilebox = mfi.tilebox(); + + ParallelFor(tilebox, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Pe(i, j, k) = ElectronPressure::get_pressure( + n0_ref, elec_temp, gamma, rho(i, j, k) + ); + }); + } +} |