Skip to content

Commit

Permalink
fixed some bugs and added checks
Browse files Browse the repository at this point in the history
  • Loading branch information
@marcoacc95
marcoacc95 committed Dec 6, 2024
1 parent cdd2752 commit daba314
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 23 deletions.
5 changes: 4 additions & 1 deletion Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -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));
Expand Down
5 changes: 5 additions & 0 deletions Source/Fluids/QdsmcParticleContainer.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_
55 changes: 49 additions & 6 deletions Source/Fluids/QdsmcParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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();
Expand All @@ -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;
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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);
}
});
}
Expand All @@ -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();
}
47 changes: 31 additions & 16 deletions Source/Fluids/Qdsmc_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,18 +32,22 @@ void gather_density_entropy(const amrex::ParticleReal x0p,
amrex::Array4<amrex::Real const> const& K_arr,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> 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;
}
}


Expand All @@ -64,19 +68,24 @@ void gather_vector_field_qdsmc(const amrex::ParticleReal x0p,
amrex::Array4<amrex::Real const> const& field_y_arr,
amrex::Array4<amrex::Real const> const& field_z_arr,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> 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);
}

}


Expand Down Expand Up @@ -112,7 +121,8 @@ void do_deposit_scalar(amrex::Array4<amrex::Real> const& scalar_field,
const amrex::ParticleReal zp,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> 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;
Expand All @@ -135,10 +145,15 @@ void do_deposit_scalar(amrex::Array4<amrex::Real> 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

0 comments on commit daba314

Please sign in to comment.