Skip to content

Commit

Permalink
Merge pull request #35 from AMReX-Microelectronics/accurate_site_id_w…
Browse files Browse the repository at this point in the history 
…ith_trap_charges

Accurate site id with trap charges
  • Loading branch information
@saurabh-s-sawant
saurabh-s-sawant authored Oct 29, 2024
2 parents 1b9d849 + 6edb660 commit 7f3dec1
Show file tree
Hide file tree
Showing 15 changed files with 190 additions and 164 deletions.
9 changes: 9 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
Exec/*
!Exec/GNUmakefile
*.ipynb_checkpoints*
*.ex
*tmp_build_dir/
*.png
*.swp
perlscripts/
input_local/
1 change: 1 addition & 0 deletions Source/Input/MacroscopicProperties/ChargeDensitySource.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ class c_ChargeDensitySource
{
return m_container.Is_Occupation_Variable();
}
amrex::Real Get_mixing_factor() { return m_container.Get_mixing_factor(); }

void Deposit(amrex::MultiFab *const p_rho_mf);

Expand Down
21 changes: 20 additions & 1 deletion Source/Input/MacroscopicProperties/ChargeDensitySource.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,12 +77,31 @@ void c_ChargeDensitySource<ParticleContainerType>::Gather(
auto &par_potential = pti.get_potential();
auto *p_par_potential = par_potential.data();

auto &par_rel_diff = pti.get_relative_difference();
auto *p_par_rel_diff = par_rel_diff.data();

amrex::Real MF = Get_mixing_factor();
amrex::Real THRESHOLD_REL_DIFF = 1.;

amrex::ParallelFor(np,
[=] AMREX_GPU_DEVICE(int p)
{
p_par_potential[p] =
amrex::Real old_phi = p_par_potential[p];

amrex::Real new_phi =
CloudInCell::Gather_Trilinear(p_par[p].pos(),
plo, dx, phi);

p_par_rel_diff[p] = fabs((new_phi - old_phi) /
(new_phi + old_phi));

// if(p_par_rel_diff[p] < THRESHOLD_REL_DIFF) {
p_par_potential[p] =
new_phi * MF + old_phi * (1. - MF);
//}
// else {
// p_par_potential[p] = new_phi;
// }
});
}
}
31 changes: 6 additions & 25 deletions Source/Input/MacroscopicProperties/PointChargeContainer.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ struct PointChargeAttributes
charge_unit,
occupation,
potential,
rel_diff,
NUM
};
};
Expand Down Expand Up @@ -65,31 +66,6 @@ class c_PointChargeContainer
const amrex::Vector<amrex::Real> &charge_units = {},
const amrex::Vector<amrex::Real> &occupations = {});

// static void Set_StaticVar(amrex::Real& member_variable, amrex::Real
// value)
//{
// std::lock_guard<std::mutex> lock(mtx);
// member_variable = value;
// }

// struct GetV0 {
// AMREX_GPU_HOST_DEVICE int operator() () const { return V0; }
// amrex::Real V0;
// };
// static GetV0 Get_V0() { return GetV0{V0}; }

// struct GetEt {
// AMREX_GPU_HOST_DEVICE int operator() () const { return Et; }
// amrex::Real Et;
// };
// static GetEt Get_Et() { return GetEt{Et}; }
// struct GetChargeDensity {
// AMREX_GPU_HOST_DEVICE int operator() () const { return charge_density;
// } amrex::Real charge_density;
// };
// static GetChargeDensity Get_Charge_Density() { return
// GetChargeDensity{charge_density}; }

amrex::Real Get_V0() { return V0; }
amrex::Real Get_Et() { return Et; }
amrex::Real Get_mixing_factor() { return mixing_factor; }
Expand Down Expand Up @@ -164,6 +140,11 @@ class c_PointChargeContainer
return this->GetStructOfArrays().GetRealData(
PCA::realPA::occupation);
}

RealVector &get_relative_difference()
{
return this->GetStructOfArrays().GetRealData(PCA::realPA::rel_diff);
}
};

private:
Expand Down
109 changes: 75 additions & 34 deletions Source/Input/MacroscopicProperties/PointChargeContainer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@
#include <AMReX_ParmParse.H>
#include <AMReX_Parser.H>

#include <filesystem>
#include <iomanip>
#include <limits>

#include "../../Output/Output.H"
Expand Down Expand Up @@ -39,7 +41,21 @@ void c_PointChargeContainer::Define_OutputFile()
int step = rCode.get_step();
std::string main_output_foldername = rOutput.get_folder_name();
std::string pc_foldername = main_output_foldername + "/point_charge";
CreateDirectory(pc_foldername);

// if (std::filesystem::exists(pc_foldername)) {
// std::cout << "Directory exists: " << pc_foldername << "\n";
// } else {
// try {
// if (CreateDirectory(pc_foldername)) {
// } else {
// std::cout << "Failed to create directory, unknown
// reason.\n";
// }
// } catch (const std::filesystem::filesystem_error& e) {
// std::cerr << "Filesystem error: " << e.what() << "\n";
// }
// }a
UtilCreateCleanDirectory(pc_foldername, false);
pc_stepwise_filename = pc_foldername + "/total_charge_vs_step.dat";
outfile_pc_step.open(pc_stepwise_filename.c_str());
outfile_pc_step << "'step', ";
Expand All @@ -62,8 +78,8 @@ void c_PointChargeContainer::Write_OutputFile()
outfile_pc_step << std::setw(10) << step;
#ifdef USE_TRANSPORT
auto &rTransport = rCode.get_TransportSolver();
outfile_pc_step << std::setw(10) << rTransport.get_Vds()
<< std::setw(10) << rTransport.get_Vgs()
outfile_pc_step << std::setw(12) << rTransport.get_Vds()
<< std::setw(12) << rTransport.get_Vgs()
<< std::setw(10) << rTransport.get_Broyden_Step() - 1;
#endif
outfile_pc_step << std::setw(15) << Get_total_charge() << std::setw(15)
Expand Down Expand Up @@ -111,7 +127,9 @@ void c_PointChargeContainer::Check_PositionBounds(
}
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
GeomUtils::Is_ID_Within_Bounds(ID, minID, maxID) == true,
"Point charge " + std::to_string(s) + " is out of bounds.");
"Point charge " + std::to_string(s) + " is out of bounds. ID: [" +
std::to_string(ID[0]) + ", " + std::to_string(ID[1]) + ", " +
std::to_string(ID[2]) + "] ");
}
}

Expand All @@ -122,23 +140,23 @@ void c_PointChargeContainer::Read_PointCharges()
amrex::Vector<amrex::Real> vec_offset(AMREX_SPACEDIM, 0.0);
queryArrWithParser(pp_main, "offset", vec_offset, 0, AMREX_SPACEDIM);

amrex::Print() << "pc.offset: ";
for (const auto &v : vec_offset) amrex::Print() << v << " ";
amrex::Print() << "pc.offset / (nm): ";
for (const auto &v : vec_offset) amrex::Print() << v * 1e9 << " ";
amrex::Print() << "\n";

amrex::Vector<amrex::Real> vec_scaling(AMREX_SPACEDIM, 1.0);
queryArrWithParser(pp_main, "scaling", vec_scaling, 0, AMREX_SPACEDIM);

amrex::Print() << "pc.scaling: ";
for (const auto &v : vec_scaling) amrex::Print() << v << " ";
amrex::Print() << "pc.scaling / (nm): ";
for (const auto &v : vec_scaling) amrex::Print() << v * 1e9 << " ";
amrex::Print() << "\n";

amrex::Vector<amrex::Real> vec_max_bound(AMREX_SPACEDIM, 1.0);
amrex::Print() << "point charge max bound: ";
amrex::Print() << "max bound / (nm): ";
for (int d = 0; d < AMREX_SPACEDIM; ++d)
{
vec_max_bound[d] = vec_offset[d] + vec_scaling[d];
amrex::Print() << vec_max_bound[d] << " ";
amrex::Print() << vec_max_bound[d] * 1e9 << " ";
}
amrex::Print() << "\n";

Expand Down Expand Up @@ -255,6 +273,8 @@ void c_PointChargeContainer::Define(
real_attribs[PCA::realPA::charge_unit] = charge_units[p];
real_attribs[PCA::realPA::occupation] = occupations[p];
real_attribs[PCA::realPA::potential] = 0.0;
real_attribs[PCA::realPA::rel_diff] =
std::numeric_limits<double>::max();

std::pair<int, int> key{0, 0}; // {grid_index, tile_index}
int lev = 0;
Expand All @@ -265,19 +285,20 @@ void c_PointChargeContainer::Define(
particle_tile.push_back_real(real_attribs);

// Print the details of the newly added particle
amrex::Print() << "Defined particle ID: " << new_particle.id()
<< ", Position: (";
amrex::Print() << std::left << "ID: " << std::setw(8)
<< new_particle.id() << ", Pos/(nm): (";
for (int j = 0; j < AMREX_SPACEDIM; ++j)
{
amrex::Print() << new_particle.pos(j);
amrex::Print() << std::setw(12) << std::setprecision(8)
<< new_particle.pos(j) * 1e9;
if (j < AMREX_SPACEDIM - 1)
{
amrex::Print() << ", ";
}
}
amrex::Print() << "), charge_unit: " << charge_units[p]
<< ", occupation: " << occupations[p]
<< ", potential: 0."
amrex::Print() << ")"
<< ", q: " << std::setw(10) << charge_units[p]
<< ", occup.: " << std::setw(10) << occupations[p]
<< "\n";
}
}
Expand All @@ -292,8 +313,8 @@ void c_PointChargeContainer::Print_Container(bool print_positions)
{
int lev = 0;
Vector<int> particle_ids;
Vector<amrex::Real> charge_units, occupations, potentials, pos_x, pos_y,
pos_z;
Vector<amrex::Real> charge_units, occupations, potentials, rel_diff, pos_x,
pos_y, pos_z;

for (c_PointChargeContainer::ParIter pti(*this, lev); pti.isValid(); ++pti)
{
Expand All @@ -308,6 +329,7 @@ void c_PointChargeContainer::Print_Container(bool print_positions)
charge_units.push_back(soa_real[PCA::realPA::charge_unit][p]);
occupations.push_back(soa_real[PCA::realPA::occupation][p]);
potentials.push_back(soa_real[PCA::realPA::potential][p]);
rel_diff.push_back(soa_real[PCA::realPA::rel_diff][p]);
if (print_positions)
{
pos_x.push_back(p_par[p].pos(0));
Expand Down Expand Up @@ -341,14 +363,15 @@ void c_PointChargeContainer::Print_Container(bool print_positions)
// allocate global vectors
Vector<int> all_particle_ids;
Vector<amrex::Real> all_charge_units, all_occupations, all_potentials,
all_pos_x, all_pos_y, all_pos_z;
all_rel_diff, all_pos_x, all_pos_y, all_pos_z;

if (ParallelDescriptor::IOProcessor())
{
all_particle_ids.resize(total_particles);
all_charge_units.resize(total_particles);
all_occupations.resize(total_particles);
all_potentials.resize(total_particles);
all_rel_diff.resize(total_particles);

if (print_positions)
{
Expand All @@ -373,6 +396,9 @@ void c_PointChargeContainer::Print_Container(bool print_positions)
ParallelDescriptor::Gatherv(potentials.data(), local_num_particles,
all_potentials.data(), recvcounts, displs,
ParallelDescriptor::IOProcessorNumber());
ParallelDescriptor::Gatherv(rel_diff.data(), local_num_particles,
all_rel_diff.data(), recvcounts, displs,
ParallelDescriptor::IOProcessorNumber());

if (print_positions)
{
Expand All @@ -393,33 +419,44 @@ void c_PointChargeContainer::Print_Container(bool print_positions)
// Printing at root process and computing total charge
total_charge = 0.;
total_charge_units = 0;
amrex::Real THRESHOLD_REL_DIFF_TO_PRINT = 1.e-5;
if (ParallelDescriptor::IOProcessor())
{
amrex::Print() << "Point Charges: \n";
int np = all_charge_units.size();
amrex::Print() << "number of charges: " << np << "\n";
for (int p = 0; p < np; ++p)
{
amrex::Print() << "ID: " << all_particle_ids[p]
<< ", charge_unit: " << all_charge_units[p]
<< ", occupation: " << all_occupations[p]
<< ", potential: " << all_potentials[p];
if (print_positions)
if (all_rel_diff[p] > THRESHOLD_REL_DIFF_TO_PRINT)
{
amrex::Print()
<< ", Position: (" << all_pos_x[p] << ", " << all_pos_y[p];
<< std::left << "ID: " << std::setw(6)
<< all_particle_ids[p] << ", charge/(e): " << std::setw(6)
<< all_charge_units[p]
<< ", occupation: " << std::setprecision(8) << std::setw(12)
<< all_occupations[p]
<< ", phi/(V): " << std::setprecision(8) << std::setw(12)
<< all_potentials[p] << ", phi_rel_diff: " << std::setw(15)
<< std::scientific << all_rel_diff[p];

if (print_positions)
{
amrex::Print()
<< ", Position: (" << std::setw(10) << all_pos_x[p]
<< ", " << std::setw(10) << all_pos_y[p];
#if AMREX_SPACEDIM == 3
amrex::Print() << ", " << all_pos_z[p];
amrex::Print() << ", " << std::setw(10) << all_pos_z[p];
#endif
amrex::Print() << ")";
amrex::Print() << ")";
}
amrex::Print() << "\n";
}
amrex::Print() << "\n";
total_charge += all_charge_units[p] * all_occupations[p];
total_charge_units += all_charge_units[p];
}
amrex::Print() << "total charge: " << total_charge << "\n";
Write_OutputFile();
}
Write_OutputFile();
}

void c_PointChargeContainer::Compute_Occupation()
Expand Down Expand Up @@ -451,19 +488,23 @@ void c_PointChargeContainer::Compute_Occupation()
auto &par_charge_unit = pti.get_charge_unit();
auto *p_charge_unit = par_charge_unit.data();

auto &par_rel_diff = pti.get_relative_difference();
auto *p_par_rel_diff = par_rel_diff.data();

amrex::Real V0 = Get_V0();
amrex::Real Et = Get_Et();
amrex::Real MF = Get_mixing_factor();

// amrex::Real MF = 0.5;
// int THRESHOLD_REL_DIFF_OCCUPATION_COMPUT = 1.e-2;
amrex::ParallelFor(np,
[=] AMREX_GPU_DEVICE(int p)
{
// if(p_par_rel_diff[p] >
// THRESHOLD_REL_DIFF_OCCUPATION_COMPUT) {
amrex::Real argument =
-(p_potential[p] - V0) / Et;
amrex::Real new_occupation =
p_occupation[p] =
MathFunctions::Sigmoid(argument);
p_occupation[p] = p_occupation[p] * MF +
(1. - MF) * new_occupation;
//}
});
}
}
18 changes: 6 additions & 12 deletions Source/Solver/Transport/NEGF/CNT.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,23 +78,17 @@ class c_CNT : public c_NEGF_Common<ComplexType[NUM_MODES]>
MatrixBlock<BlkType> &gr, const ComplexType EmU) final;

public:
struct Get1DSiteID
static int get_1D_site_id(int par_id_local_to_NS)
{
AMREX_GPU_HOST_DEVICE int operator()(int global_par_id) const
{
return static_cast<int>(
amrex::Math::floor((global_par_id - 1) / type_id[0]));
}
amrex::Array<int, 2> type_id;
};

static Get1DSiteID get_1D_site_id() { return Get1DSiteID{type_id}; }
return static_cast<int>(
amrex::Math::floor((par_id_local_to_NS - 1) / type_id[0]));
}

struct GetAtomIDAtSite
{
AMREX_GPU_HOST_DEVICE int operator()(int global_par_id) const
AMREX_GPU_HOST_DEVICE int operator()(int par_id_local_to_NS) const
{
return static_cast<int>((global_par_id - 1) % type_id[0]);
return static_cast<int>((par_id_local_to_NS - 1) % type_id[0]);
}
amrex::Array<int, 2> type_id;
};
Expand Down
Loading

0 comments on commit 7f3dec1

Please sign in to comment.