From daba3143b4f340670637406e51acee793a2a0b87 Mon Sep 17 00:00:00 2001 From: Marco Acciarri Date: Fri, 6 Dec 2024 13:02:56 -0800 Subject: [PATCH] fixed some bugs and added checks --- .../FieldSolver/WarpXPushFieldsHybridPIC.cpp | 5 +- Source/Fluids/QdsmcParticleContainer.H | 5 ++ Source/Fluids/QdsmcParticleContainer.cpp | 55 +++++++++++++++++-- Source/Fluids/Qdsmc_K.H | 47 ++++++++++------ 4 files changed, 89 insertions(+), 23 deletions(-) diff --git a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp index f78fdf7933e..191d3ff606e 100644 --- a/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp +++ b/Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp @@ -181,6 +181,9 @@ void WarpX::HybridPICEvolveFields () m_hybrid_pic_model.get(), finest_level); + // Reset qdsmc particles positions to x0,y0,z0 and rest of attributes to 0 and redistribute + qdsmc_hybrid_electron_pc->ResetParticles(finest_level); + // Set fictitious electron particles velocities qdsmc_hybrid_electron_pc->SetV(finest_level, *m_fields.get(hybrid_electron_fl->name_mf_NU, Direction{0}, finest_level), @@ -193,7 +196,7 @@ void WarpX::HybridPICEvolveFields () *m_fields.get(FieldType::rho_fp, finest_level)); // Push fictitious electron particles - //qdsmc_hybrid_electron_pc->PushX(finest_level, dt[0]); + qdsmc_hybrid_electron_pc->PushX(finest_level, dt[0]); // Deposit entropy from qdsmc qdsmc_hybrid_electron_pc->DepositK(finest_level, *m_fields.get(hybrid_electron_fl->name_mf_K, finest_level)); diff --git a/Source/Fluids/QdsmcParticleContainer.H b/Source/Fluids/QdsmcParticleContainer.H index 0a0bcd066ae..3ccc6ce8ab8 100644 --- a/Source/Fluids/QdsmcParticleContainer.H +++ b/Source/Fluids/QdsmcParticleContainer.H @@ -105,6 +105,11 @@ public: // using linear interpolation void DepositK(int lev, amrex::MultiFab &Kfield); + // Reset particles positions to x0,y0,z0 + // and v, np, K to 0 + // then calls Redistribute() + void ResetParticles(int lev); + }; #endif // WARPX_QdsmcParticleContainer_H_ diff --git a/Source/Fluids/QdsmcParticleContainer.cpp b/Source/Fluids/QdsmcParticleContainer.cpp index 5b62f6dc5ac..5b06bc68979 100644 --- a/Source/Fluids/QdsmcParticleContainer.cpp +++ b/Source/Fluids/QdsmcParticleContainer.cpp @@ -272,7 +272,7 @@ QdsmcParticleContainer::SetV (int lev, amrex::Real vyp; amrex::Real vzp; - gather_vector_field_qdsmc(part_x0[ip], part_y0[ip], part_z0[ip], vxp, vyp, vzp, arrUxfield, arrUyfield, arrUzfield, plo, dinv); + gather_vector_field_qdsmc(part_x0[ip], part_y0[ip], part_z0[ip], vxp, vyp, vzp, arrUxfield, arrUyfield, arrUzfield, plo, dinv, box); part_vx[ip] = vxp; part_vy[ip] = vyp; @@ -307,8 +307,6 @@ QdsmcParticleContainer::SetK (int lev, amrex::Box box = amrex::convert( tilebox, ix_type_Kfield ); box.grow(Kfield.nGrowVect()); - const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - amrex::ParticleReal* const AMREX_RESTRICT part_x0 = attribs[QdsmcPIdx::x_node].dataPtr(); amrex::ParticleReal* const AMREX_RESTRICT part_y0 = attribs[QdsmcPIdx::y_node].dataPtr(); amrex::ParticleReal* const AMREX_RESTRICT part_z0 = attribs[QdsmcPIdx::z_node].dataPtr(); @@ -326,7 +324,7 @@ QdsmcParticleContainer::SetK (int lev, amrex::Real n_p; amrex::Real kn_p; - gather_density_entropy(part_x0[ip], part_y0[ip], part_z0[ip], n_p, kn_p, arrrhofield, arrKfield, plo, dinv, cell_volume); + gather_density_entropy(part_x0[ip], part_y0[ip], part_z0[ip], n_p, kn_p, arrrhofield, arrKfield, plo, dinv, cell_volume, box); part_np_real[ip] = n_p; part_entropy[ip] = kn_p; @@ -387,6 +385,48 @@ QdsmcParticleContainer::PushX (int lev, amrex::Real dt) } +void +QdsmcParticleContainer::ResetParticles(int lev) +{ + for (iterator pti(*this, lev); pti.isValid(); ++pti) + { + auto const np = pti.numParticles(); + auto& attribs = pti.GetStructOfArrays().GetRealData(); + + amrex::ParticleReal* const AMREX_RESTRICT part_x0 = attribs[QdsmcPIdx::x_node].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_y0 = attribs[QdsmcPIdx::y_node].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_z0 = attribs[QdsmcPIdx::z_node].dataPtr(); + + amrex::ParticleReal* const AMREX_RESTRICT part_x = attribs[QdsmcPIdx::x].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_y = attribs[QdsmcPIdx::y].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_z = attribs[QdsmcPIdx::z].dataPtr(); + + amrex::ParticleReal* const AMREX_RESTRICT part_vx = attribs[QdsmcPIdx::vx].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_vy = attribs[QdsmcPIdx::vy].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_vz = attribs[QdsmcPIdx::vz].dataPtr(); + + amrex::ParticleReal* const AMREX_RESTRICT part_entropy = attribs[QdsmcPIdx::entropy].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT part_np_real = attribs[QdsmcPIdx::np_real].dataPtr(); + + amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip) + { + part_x[ip] = part_x0[ip]; + part_y[ip] = part_y0[ip]; + part_z[ip] = part_z0[ip]; + + part_vx[ip] = 0; + part_vy[ip] = 0; + part_vz[ip] = 0; + + part_entropy[ip] = 0; + part_np_real[ip] = 0; + }); + } + + Redistribute(); +} + + // Generalize this function to --> DepositScalar void QdsmcParticleContainer::DepositK(int lev, amrex::MultiFab &Kfield) @@ -430,7 +470,7 @@ QdsmcParticleContainer::DepositK(int lev, amrex::MultiFab &Kfield) // avoid launching kernel for "empty" particles if(part_entropy[ip]>0) { - do_deposit_scalar(arrKField, part_x[ip], part_y[ip], part_z[ip], plo, dinv, part_entropy[ip]); + do_deposit_scalar(arrKField, part_x[ip], part_y[ip], part_z[ip], plo, dinv, part_entropy[ip], box); } }); } @@ -439,5 +479,8 @@ QdsmcParticleContainer::DepositK(int lev, amrex::MultiFab &Kfield) Kfield, 0, Kfield.nComp(), Kfield.nGrowVect(), Kfield.nGrowVect(), WarpX::do_single_precision_comms, period); - amrex::Gpu::streamSynchronize(); + // Add fill boundary ? + //Kfield.FillBoundary(Kfield.nGrowVect(), period); + + //amrex::Gpu::streamSynchronize(); } diff --git a/Source/Fluids/Qdsmc_K.H b/Source/Fluids/Qdsmc_K.H index dcebc1003f0..cb4e9bb5bef 100644 --- a/Source/Fluids/Qdsmc_K.H +++ b/Source/Fluids/Qdsmc_K.H @@ -32,18 +32,22 @@ void gather_density_entropy(const amrex::ParticleReal x0p, amrex::Array4 const& K_arr, amrex::GpuArray const& plo, const amrex::XDim3 & dinv, - const amrex::Real cell_volume) + const amrex::Real cell_volume, + const amrex::Box& bx) { - amrex::Real x = (x0p - plo[0])*dinv.x; - amrex::Real y = (y0p - plo[1])*dinv.y; - amrex::Real z = (z0p - plo[2])*dinv.z; + amrex::Real x = (x0p - plo[0])*dinv.x - 0.5; + amrex::Real y = (y0p - plo[1])*dinv.y - 0.5; + amrex::Real z = (z0p - plo[2])*dinv.z - 0.5; int i0p = std::floor(x); int j0p = std::floor(y); int k0p = std::floor(z); - n_p = rho_arr(i0p,j0p,k0p)*cell_volume/PhysConst::q_e; - kn_p = K_arr(i0p,j0p,k0p)*n_p; + amrex::IntVect c(i0p, j0p, k0p); + if (bx.contains(c)){ // Check if the cell is within the tile's box + n_p = rho_arr(i0p,j0p,k0p)*cell_volume/PhysConst::q_e; + kn_p = K_arr(i0p,j0p,k0p)*n_p; + } } @@ -64,19 +68,24 @@ void gather_vector_field_qdsmc(const amrex::ParticleReal x0p, amrex::Array4 const& field_y_arr, amrex::Array4 const& field_z_arr, amrex::GpuArray const& plo, - const amrex::XDim3 & dinv) + const amrex::XDim3 & dinv, + const amrex::Box& bx) { - amrex::Real x = (x0p - plo[0])*dinv.x; - amrex::Real y = (y0p - plo[1])*dinv.y; - amrex::Real z = (z0p - plo[2])*dinv.z; + amrex::Real x = (x0p - plo[0])*dinv.x - 0.5; + amrex::Real y = (y0p - plo[1])*dinv.y - 0.5; + amrex::Real z = (z0p - plo[2])*dinv.z - 0.5; int i0p = std::floor(x); int j0p = std::floor(y); int k0p = std::floor(z); - field_x_p = field_x_arr(i0p,j0p,k0p); - field_y_p = field_y_arr(i0p,j0p,k0p); - field_z_p = field_z_arr(i0p,j0p,k0p); + amrex::IntVect c(i0p, j0p, k0p); + if (bx.contains(c)){ // Check if the cell is within the tile's box + field_x_p = field_x_arr(i0p,j0p,k0p); + field_y_p = field_y_arr(i0p,j0p,k0p); + field_z_p = field_z_arr(i0p,j0p,k0p); + } + } @@ -112,7 +121,8 @@ void do_deposit_scalar(amrex::Array4 const& scalar_field, const amrex::ParticleReal zp, amrex::GpuArray const& plo, const amrex::XDim3 & dinv, - const amrex::ParticleReal valuep) + const amrex::ParticleReal valuep, + const amrex::Box& bx) { // calculate di, dj, dk for particle in cell amrex::Real x = (xp - plo[0])*dinv.x - 0.5; @@ -135,10 +145,15 @@ void do_deposit_scalar(amrex::Array4 const& scalar_field, for (int ioff = 0; ioff < 2; ++ioff){ for (int joff = 0; joff < 2; ++joff){ for (int koff = 0; koff < 2; ++koff){ - amrex::Gpu::Atomic::AddNoRet(&scalar_field(i0p+ioff,j0p+joff,k0p+koff), valuep*sx[ioff]*sy[joff]*sz[koff]); + + amrex::IntVect c(i0p+ioff, j0p+joff, k0p+koff); + if (bx.contains(c)){ // Check if the cell is within the tile's box + amrex::Gpu::Atomic::AddNoRet(&scalar_field(i0p+ioff,j0p+joff,k0p+koff), valuep*sx[ioff]*sy[joff]*sz[koff]); + } + } } } } -#endif //WARPX_Qdsmc_K_h +#endif //WARPX_Qdsmc_K_h \ No newline at end of file