Skip to content

Commit

Permalink
fixed bugs in fluid container qdsmc functions, added safety n_floor c…
Browse files Browse the repository at this point in the history 
…heck to Ke and Ue calculation
  • Loading branch information
@marcoacc95
marcoacc95 committed Nov 29, 2024
1 parent 4191b7d commit e6fae85
Show file tree
Hide file tree
Showing 8 changed files with 206 additions and 132 deletions.
1 change: 1 addition & 0 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -201,6 +201,7 @@ public:
amrex::GpuArray<int, 3> Ez_IndexType;

bool m_solve_electron_energy_equation;

};

/**
Expand Down
20 changes: 4 additions & 16 deletions Source/FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -235,26 +235,14 @@ void HybridPICModel::InitData ()
}

const auto elec_temp = m_elec_temp;
// Initialize electron temperature multifab
// This assumes an uniform electron temperature in the domain
// ...
// 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(*warpx.m_fields.get("fluid_temperature_electrons_hybrid", warpx.finestLevel()), TilingIfNotGPU()); mfi.isValid(); ++mfi )
{
Array4<Real> const& Te = warpx.m_fields.get("fluid_temperature_electrons_hybrid", warpx.finestLevel())->array(mfi);
const Box& tilebox = mfi.tilebox();

ParallelFor(tilebox, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
Te(i, j, k) = elec_temp;
});
}
// Initialize electron temperature multifab
// This assumes an uniform electron temperature in the domain
warpx.m_fields.get("fluid_temperature_electrons_hybrid", warpx.finestLevel())->setVal(elec_temp);

// Fill Boundaries in electron temperature multifab
warpx.m_fields.get("fluid_temperature_electrons_hybrid", warpx.finestLevel())->FillBoundary(warpx.Geom(warpx.finestLevel()).periodicity());

}

void HybridPICModel::GetCurrentExternal ()
Expand Down
29 changes: 14 additions & 15 deletions Source/FieldSolver/WarpXPushFieldsHybridPIC.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,36 +171,35 @@ void WarpX::HybridPICEvolveFields ()

// Initialize Ke = ne^(1-gamma)*Te
// Move up, shouldn't Ke be defined with both ne and Te at n instead of ne at n+1 ?
hybrid_electron_fl->HybridInitializeKe(m_fields, m_hybrid_pic_model->m_gamma , finest_level);
hybrid_electron_fl->HybridInitializeKe(m_fields, m_hybrid_pic_model->m_gamma, m_hybrid_pic_model->m_n_floor, finest_level);

// Update Ue in electron fluid container
// check at what step this should be calculated, here Ue is at n+1
// maybe move up to be consistent with Ke initialization
hybrid_electron_fl->HybridInitializeUe(m_fields,
*current_fp_temp[finest_level][0],
*current_fp_temp[finest_level][1],
*current_fp_temp[finest_level][2],
m_hybrid_pic_model.get(), finest_level);
current_fp_temp[finest_level],
m_hybrid_pic_model.get(),
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),
*m_fields.get(hybrid_electron_fl->name_mf_NU, Direction{1}, finest_level),
*m_fields.get(hybrid_electron_fl->name_mf_NU, Direction{2}, finest_level));
//qdsmc_hybrid_electron_pc->SetV(finest_level,
// *m_fields.get(hybrid_electron_fl->name_mf_NU, Direction{0}, finest_level),
// *m_fields.get(hybrid_electron_fl->name_mf_NU, Direction{1}, finest_level),
// *m_fields.get(hybrid_electron_fl->name_mf_NU, Direction{2}, finest_level));

// Set fictitious electron particles entropy
qdsmc_hybrid_electron_pc->SetK(finest_level,
*m_fields.get(hybrid_electron_fl->name_mf_K, finest_level),
*m_fields.get(FieldType::rho_fp, finest_level));
//qdsmc_hybrid_electron_pc->SetK(finest_level,
// *m_fields.get(hybrid_electron_fl->name_mf_K, finest_level),
// *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));
//qdsmc_hybrid_electron_pc->DepositK(finest_level, *m_fields.get(hybrid_electron_fl->name_mf_K, finest_level));

// Update Te after QDSMC solver:
hybrid_electron_fl->HybridUpdateTe(m_fields, m_hybrid_pic_model->m_gamma, finest_level);
//hybrid_electron_fl->HybridUpdateTe(m_fields, m_hybrid_pic_model->m_gamma, m_hybrid_pic_model->m_n_floor, finest_level);
}

// Calculate the electron pressure at t=n+1
Expand Down
12 changes: 4 additions & 8 deletions Source/Fluids/QdsmcParticleContainer.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,10 @@
#include <AMReX_AmrCoreFwd.H>
#include <AMReX_Vector.H>

#include "Particles/ParticleBoundaries.H"
#include <ablastr/fields/MultiFabRegister.H>
//#include "Particles/ParticleBoundaries.H"
//#include <ablastr/fields/MultiFabRegister.H>

#include "QdsmcParticleContainer_fwd.H"
//#include "QdsmcParticleContainer_fwd.H"

/**
* This enumerated struct is used to index the qdsmc fictitious
Expand All @@ -38,10 +38,6 @@ struct QdsmcPIdx
z, z_node,
vx, vy, vz,
entropy, np_real,

#ifdef WARPX_DIM_RZ
theta, ///< RZ needs all three position components
#endif
nattribs
};
};
Expand Down Expand Up @@ -74,7 +70,7 @@ public:

//! similar to WarpXParticleContainer::AddNParticles

void AddNParticles (int lev, amrex::Long n,
void AddNParticles (int lev, amrex::Long np,
amrex::Vector<amrex::ParticleReal> const & x,
amrex::Vector<amrex::ParticleReal> const & y,
amrex::Vector<amrex::ParticleReal> const & z);
Expand Down
163 changes: 107 additions & 56 deletions Source/Fluids/QdsmcParticleContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@
#include <AMReX_ParticleUtil.H>
#include <AMReX_Utility.H>

#include <AMReX_Print.H>

#include <algorithm>
#include <cmath>
Expand All @@ -63,7 +64,6 @@ QdsmcParticleContainer::QdsmcParticleContainer (AmrCore* amr_core)
: ParticleContainerPureSoA<QdsmcPIdx::nattribs, 0>(amr_core->GetParGDB())
{
SetParticleSize();
//InitParticles(0); // only level 0 is used in HybridSolver, commented since this is already called in WarpX.cpp
}


Expand All @@ -78,68 +78,121 @@ QdsmcParticleContainer::AddNParticles (int lev, amrex::Long n,
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(y.size() == n,"y.size() != # of qdsmc particles to add");
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(z.size() == n,"z.size() != # of qdsmc particles to add");

if (n <= 0){
Redistribute();
return;
}
// have to resize here, not in the constructor because grids have not
// been built when constructor was called.
reserveData();
resizeData();

long ibegin = 0;
long iend = n;

const int myproc = amrex::ParallelDescriptor::MyProc();
const int nprocs = amrex::ParallelDescriptor::NProcs();
const auto navg = n/nprocs;
const auto nleft = n - navg * nprocs;
if (myproc < nleft) {
ibegin = myproc*(navg+1);
iend = ibegin + navg+1;
} else {
ibegin = myproc*navg + nleft;
iend = ibegin + navg;
}


//if (np <= 0){
// Redistribute();
// return;
//}

//std::ofstream ofs("xyz_InitParticles.txt", std::ofstream::out); // REMOVE
//for(int i=0; i<np; i++){ //REMOVE
// Print(ofs) << " x " << x[i] << " y " << y[i] << " z " << z[i] << std::endl; // REMOVE
//} //REMOVE
//ofs.close(); // REMOVE
// up to here, x,y,z vectors have the correct values

// Add to grid 0 and tile 0
// Redistribute() will move them to proper places.
auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0);

// Creates a temporary tile to obtain data from simulation. This data
// is then coppied to the permament tile which is stored on the particle
// (particle_tile).
using PinnedTile = typename ContainerLike<amrex::PinnedArenaAllocator>::ParticleTileType;

PinnedTile pinned_tile;
pinned_tile.define(NumRuntimeRealComps(), NumRuntimeIntComps());

for (int i = 0; i < n; i++)
const std::size_t np = iend-ibegin;

std::ofstream ofsID("IDs.txt", std::ofstream::out); // REMOVE
for (auto i = ibegin; i < iend; ++i)
{
auto & idcpu_data = pinned_tile.GetStructOfArrays().GetIdCPUData();
idcpu_data.push_back(amrex::SetParticleIDandCPU(ParticleType::NextID(), ParallelDescriptor::MyProc()));
amrex::Long current_id = ParticleType::NextID();
idcpu_data.push_back(amrex::SetParticleIDandCPU(current_id, ParallelDescriptor::MyProc()));
Print(ofsID) << " i " << i <<" current_id " << current_id << std::endl; // REMOVE
}
ofsID.close(); // REMOVE

// write Real attributes (SoA) to particle initialized zero
DefineAndReturnParticleTile(0, 0, 0);
//
// Check here if the idcpu_data is correct and compare with physical particle container
// ...
// ...

std::ofstream ofs0("some_params.txt", std::ofstream::out); // REMOVE
Print(ofs0) << " ibegin " << ibegin << " iend " << iend << " np " << np << " NReal " << NReal << std::endl; // REMOVE
Print(ofs0) << " NumRuntimeRealComps " << NumRuntimeRealComps() << " NumRuntimeIntComps " << NumRuntimeIntComps() << std::endl; // REMOVE
ofs0.close(); // REMOVE

//std::ofstream ofs("xyz_InitParticles.txt", std::ofstream::out); // REMOVE
//for(int i=0; i<np; i++){ //REMOVE
// Print(ofs) << " x " << x[i] << " y " << y[i] << " z " << z[i] << std::endl; // REMOVE
//} //REMOVE
//ofs.close(); // REMOVE

// for RZ write theta value
#ifdef WARPX_DIM_RZ
pinned_tile.push_back_real(QdsmcPIdx::theta, n, 0.0);
#endif
#if !defined (WARPX_DIM_1D_Z)
pinned_tile.push_back_real(QdsmcPIdx::x, x);
pinned_tile.push_back_real(QdsmcPIdx::x_node, x);
pinned_tile.push_back_real(QdsmcPIdx::x, x.data() + ibegin, x.data() + iend);
pinned_tile.push_back_real(QdsmcPIdx::x_node, x.data() + ibegin, x.data() + iend);
#endif
#if defined (WARPX_DIM_3D)
pinned_tile.push_back_real(QdsmcPIdx::y, y);
pinned_tile.push_back_real(QdsmcPIdx::y_node, y);
pinned_tile.push_back_real(QdsmcPIdx::y, y.data() + ibegin, y.data() + iend);
pinned_tile.push_back_real(QdsmcPIdx::y_node, y.data() + ibegin, y.data() + iend);
#endif
pinned_tile.push_back_real(QdsmcPIdx::z, z);
pinned_tile.push_back_real(QdsmcPIdx::z_node, z);
pinned_tile.push_back_real(QdsmcPIdx::vx, n, 0.0);
pinned_tile.push_back_real(QdsmcPIdx::vy, n, 0.0);
pinned_tile.push_back_real(QdsmcPIdx::vz, n, 0.0);
pinned_tile.push_back_real(QdsmcPIdx::entropy, n, 0.0);
pinned_tile.push_back_real(QdsmcPIdx::np_real, n, 0.0);

const auto old_np = particle_tile.numParticles();
const auto new_np = old_np + pinned_tile.numParticles();
pinned_tile.push_back_real(QdsmcPIdx::z, z.data() + ibegin, z.data() + iend);
pinned_tile.push_back_real(QdsmcPIdx::z_node, z.data() + ibegin, z.data() + iend);
pinned_tile.push_back_real(QdsmcPIdx::vx, np, 0.0_prt);
pinned_tile.push_back_real(QdsmcPIdx::vy, np, 0.0_prt);
pinned_tile.push_back_real(QdsmcPIdx::vz, np, 0.0_prt);
pinned_tile.push_back_real(QdsmcPIdx::entropy, np, 0.0_prt);
pinned_tile.push_back_real(QdsmcPIdx::np_real, np, 0.0_prt);

if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){
DefineAndReturnParticleTile(0, 0, 0);
}

pinned_tile.resize(np);

auto old_np = particle_tile.numParticles();
auto new_np = old_np + pinned_tile.numParticles();
particle_tile.resize(new_np);
amrex::copyParticles(
particle_tile, pinned_tile, 0, old_np, pinned_tile.numParticles());
particle_tile, pinned_tile, 0, old_np, pinned_tile.numParticles()
);

// Move particles to their appropriate tiles
Redistribute();

// Remove particles that are inside the embedded boundaries here
/**
*
*
*
*/
// Remove particles that are inside the embedded boundaries
//#ifdef AMREX_USE_EB
// if (EB::enabled()) {
// auto & warpx = WarpX::GetInstance();
// scrapeParticlesAtEB(
// *this,
// warpx.m_fields.get_mr_levels(FieldType::distance_to_eb, warpx.finestLevel()),
// ParticleBoundaryProcess::Absorb());
// deleteInvalidParticles();
// }
//#endif
}


Expand All @@ -152,7 +205,7 @@ QdsmcParticleContainer::InitParticles (int lev)

const amrex::Real* dx = warpx.Geom(lev).CellSize();

// Read this from domain ?
// Number of cells in each direction
int nx = (probhi[0] - problo[0])/dx[0];
int ny = (probhi[1] - problo[1])/dx[1];
int nz = (probhi[2] - problo[2])/dx[2];
Expand All @@ -164,6 +217,13 @@ QdsmcParticleContainer::InitParticles (int lev)
amrex::Vector<amrex::ParticleReal> ypos;
amrex::Vector<amrex::ParticleReal> zpos;

//std::ofstream ofs1("variables_InitParticles.txt", std::ofstream::out); // REMOVE
//Print(ofs1) << " problo[0] " << problo[0] << " problo[1] " << problo[1] << " problo[2] " << problo[2] << std::endl; // REMOVE
//Print(ofs1) << " probhi[0] " << probhi[0] << " probhi[1] " << probhi[1] << " probhi[2] " << probhi[2] << std::endl; // REMOVE
//Print(ofs1) << " dx[0] " << dx[0] << " dx[1] " << dx[1] << " dx[2] " << dx[2] << std::endl; // REMOVE
//Print(ofs1) << " nx " << nx << " ny " << ny << " nz " << nz << std::endl; // REMOVE
//ofs1.close(); // REMOVE

// for now, only one MPI rank adds fictitious particles
// this is done only once in an entire simulation
if (ParallelDescriptor::IOProcessor())
Expand All @@ -174,30 +234,21 @@ QdsmcParticleContainer::InitParticles (int lev)
{
for ( int k = 0; k <= nz; k++)
{
// cell centered
//amrex::Real x = problo[0] + (i+0.5)*dx[0];
//amrex::Real y = problo[1] + (j+0.5)*dx[1];
//amrex::Real z = problo[2] + (k+0.5)*dx[2];

// nodal
amrex::Real x = problo[0] + i*dx[0];
amrex::Real y = problo[1] + j*dx[1];
amrex::Real z = problo[2] + k*dx[2];

if( x <= probhi[0] && y <= probhi[1] && z <= probhi[2])
{
xpos.push_back(x);
ypos.push_back(y);
zpos.push_back(z);

n_to_add++;
}
amrex::ParticleReal x = problo[0] + i*dx[0];
amrex::ParticleReal y = problo[1] + j*dx[1];
amrex::ParticleReal z = problo[2] + k*dx[2];

xpos.push_back(x);
ypos.push_back(y);
zpos.push_back(z);

n_to_add++;
}
}
}
}

AddNParticles (0, n_to_add, xpos, ypos, zpos);
AddNParticles(0, n_to_add, xpos, ypos, zpos);
}


Expand Down Expand Up @@ -338,7 +389,7 @@ QdsmcParticleContainer::PushX (int lev, amrex::Real dt)
});
}

Redistribute();
//Redistribute();
// search for maximum part_dx/part_dy/part_dz and assert if larger than dx/dy/dz
// ...
// ...
Expand Down
8 changes: 4 additions & 4 deletions Source/Fluids/WarpXFluidContainer.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,16 +134,16 @@ public:
void DepositCharge (ablastr::fields::MultiFabRegister& m_fields, amrex::MultiFab &rho, int lev, int icomp = 0);

// HybridUpdateUe calculates Ue from Je = J_ampere - Ji and saves in NU for hybrid solver
void HybridInitializeUe (ablastr::fields::MultiFabRegister& fields,
amrex::MultiFab &ji_x, amrex::MultiFab &ji_y, amrex::MultiFab &ji_z,
void HybridInitializeUe (ablastr::fields::MultiFabRegister& m_fields,
ablastr::fields::VectorField const& Ji,
HybridPICModel const* hybrid_model, int lev);

// Sets Ke in fluid container MF to be used by qdsmc particles
void HybridInitializeKe (ablastr::fields::MultiFabRegister& m_fields, amrex::Real gamma, int lev);
void HybridInitializeKe (ablastr::fields::MultiFabRegister& m_fields, amrex::Real gamma, amrex::Real n_floor, int lev);

// UpdateTe takes the already interpolated Ke values from the grid after pushing the fictitious particles
// and calculates Te at n+1
void HybridUpdateTe (ablastr::fields::MultiFabRegister& m_fields, amrex::Real gamma, int lev);
void HybridUpdateTe (ablastr::fields::MultiFabRegister& m_fields, amrex::Real gamma, amrex::Real n_floor, int lev);

[[nodiscard]] amrex::Real getCharge () const {return charge;}
[[nodiscard]] amrex::Real getMass () const {return mass;}
Expand Down
Loading

0 comments on commit e6fae85

Please sign in to comment.