Skip to content

add nucleation parser and sro params #52

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
prkkumar wants to merge 6 commits into
base: development
Choose a base branch
from
Open

add nucleation parser and sro params #52

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
7b51e24
add nucleation parser and sro params
Mar 29, 2024
491809e
option for both complete and partial ionization
Mar 29, 2024
c4417f5
uniform initialization of P
Apr 30, 2024
667ccb8
add warning message for the new initialization
Apr 30, 2024
47caecc
noise
May 10, 2024
d50c934
noise
Jun 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 26 additions & 9 deletions Source/FerroX.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,12 +212,15 @@ AMREX_GPU_MANAGED amrex::Real FerroX::intrinsic_carrier_concentration;
AMREX_GPU_MANAGED int FerroX::use_Fermi_Dirac;
AMREX_GPU_MANAGED int FerroX::use_work_function;
AMREX_GPU_MANAGED amrex::Real FerroX::metal_work_function;
AMREX_GPU_MANAGED int FerroX::use_sro;
AMREX_GPU_MANAGED int FerroX::complete_ionization;

// P and Phi Bc
AMREX_GPU_MANAGED amrex::Real FerroX::lambda;
AMREX_GPU_MANAGED amrex::GpuArray<int, AMREX_SPACEDIM> FerroX::P_BC_flag_lo;
AMREX_GPU_MANAGED amrex::GpuArray<int, AMREX_SPACEDIM> FerroX::P_BC_flag_hi;
AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> FerroX::Remnant_P;
AMREX_GPU_MANAGED amrex::Real FerroX::noise_amplitude;

//problem type : initialization of P for 2D/3D/convergence problems
AMREX_GPU_MANAGED int FerroX::prob_type;
Expand Down Expand Up @@ -358,6 +361,9 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
num_Vapp_max = 1;
pp.query("num_Vapp_max",num_Vapp_max);

noise_amplitude = 0.0;
pp.query("noise_amplitude",noise_amplitude);

//stack dimensions in 3D. This is an alternate way of initializing the device geometry, which works in simpler scenarios.
//A more general way of initializing device geometry is accomplished through masks which use function parsers

Expand Down Expand Up @@ -447,16 +453,23 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
}
}

// For Silicon:
// Nc = 2.8e25 m^-3
// Nv = 1.83e25 m^-3
// Band gap Eg = 1.12eV
use_sro = 0;
pp.query("use_sro",use_sro);

if (use_sro == 0){
//For Silicon:
Nc = 2.8e25; // m^-3
Nv = 1.83e25; //m^-3
bandgap = 1.12; //eV
affinity = 4.05; //eV
// 1eV = 1.602e-19 J

Nc = 2.8e25;
Nv = 1.83e25;
bandgap = 1.12; //eV
affinity = 4.05; //eV
} else {
//for SRO:
Nc = 4.640317153991026e+26;
Nv = 1.9697595861578862e+30;
bandgap = 1.25; //eV
affinity = 7.2; //eV
}
q = 1.602e-19;
kb = 1.38e-23; // Boltzmann constant
T = 300; // Room Temp
Expand All @@ -466,6 +479,10 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
donor_doping = 0.0;
pp.query("donor_doping",donor_doping);

complete_ionization = 1;
pp.query("complete_ionization",complete_ionization);

//For partial ionization:
//The most common acceptor dopant in bulk Si is boron (B), which has Ea = 44 meV
//The most common donors in bulk Si are phosphorus (P) and arsenic (As),
//which have ionization energies of Ed = 46 meV and 54 meV, respectively.
Expand Down
3 changes: 3 additions & 0 deletions Source/FerroX_namespace.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ namespace FerroX {
extern AMREX_GPU_MANAGED int use_Fermi_Dirac;
extern AMREX_GPU_MANAGED int use_work_function;
extern AMREX_GPU_MANAGED amrex::Real metal_work_function;
extern AMREX_GPU_MANAGED int use_sro;
extern AMREX_GPU_MANAGED int complete_ionization;


// P and Phi Bc
Expand All @@ -59,6 +61,7 @@ namespace FerroX {
extern AMREX_GPU_MANAGED amrex::GpuArray<int, AMREX_SPACEDIM> P_BC_flag_lo;
extern AMREX_GPU_MANAGED amrex::GpuArray<int, AMREX_SPACEDIM> P_BC_flag_hi;
extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> Remnant_P;
extern AMREX_GPU_MANAGED amrex::Real noise_amplitude;


//problem type : initialization of P for 2D/3D/convergence problems
Expand Down
5 changes: 4 additions & 1 deletion Source/Solver/ChargeDensity.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,8 @@ void ComputeRho(MultiFab& PoissonPhi,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
const MultiFab& MaterialMask);
MultiFab& acceptor_den,
MultiFab& donor_den,
const MultiFab& MaterialMask,
const Geometry& geom);

48 changes: 32 additions & 16 deletions Source/Solver/ChargeDensity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,15 +16,12 @@ void ComputeRho(MultiFab& PoissonPhi,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
const MultiFab& MaterialMask)
MultiFab& acceptor_den,
MultiFab& donor_den,
const MultiFab& MaterialMask,
const Geometry& geom)
{

//Define acceptor and donor multifabs for doping and fill them with zero.
MultiFab acceptor_den(rho.boxArray(), rho.DistributionMap(), 1, 0);
MultiFab donor_den(rho.boxArray(), rho.DistributionMap(), 1, 0);
acceptor_den.setVal(0.);
donor_den.setVal(0.);

// loop over boxes
for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi)
{
Expand All @@ -47,7 +44,12 @@ void ComputeRho(MultiFab& PoissonPhi,

//Following: http://dx.doi.org/10.1063/1.4825209

amrex::Real Ef = 0.0;
amrex::Real Ef;
if (mask(i,j,k) == 2.0){
Ef = 0.0;
} else {
Ef = -q*Phi_Bc_hi;
}
amrex::Real Eg = bandgap;
amrex::Real Chi = affinity;
amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc/Nv)/q;
Expand All @@ -62,8 +64,8 @@ void ComputeRho(MultiFab& PoissonPhi,
//g_D is the donor ground state degeneracy factor and is equal to 2
//because a donor level can accept one electron with either spin or can have no electron when filled.

amrex::Real g_A = 4.0;
amrex::Real g_D = 2.0;
amrex::Real g_A = 4.0;
amrex::Real g_D = 2.0;

amrex::Real Ea = acceptor_ionization_energy;
amrex::Real Ed = donor_ionization_energy;
Expand All @@ -76,18 +78,27 @@ void ComputeRho(MultiFab& PoissonPhi,
e_den_arr(i,j,k) = Nc*FD_half(eta_n);
hole_den_arr(i,j,k) = Nv*FD_half(eta_p);

acceptor_den_arr(i,j,k) = acceptor_doping/(1.0 + g_A*exp((-q*Ef + q*Ea + q*phi_ref - q*Chi - q*Eg - q*phi(i,j,k))/(kb*T)));
donor_den_arr(i,j,k) = donor_doping/(1.0 + g_D*exp( (q*Ef + q*Ed - q*phi_ref + q*Chi + q*phi(i,j,k)) / (kb*T) ));
if (complete_ionization == 0){
acceptor_den_arr(i,j,k) = acceptor_doping/(1.0 + g_A*exp((-q*Ef + q*Ea + q*phi_ref - q*Chi - q*Eg - q*phi(i,j,k))/(kb*T)));
donor_den_arr(i,j,k) = donor_doping/(1.0 + g_D*exp( (q*Ef + q*Ed - q*phi_ref + q*Chi + q*phi(i,j,k)) / (kb*T) ));
} else {
acceptor_den_arr(i,j,k) = acceptor_doping;
donor_den_arr(i,j,k) = donor_doping;
}

} else {

//Maxwell-Boltzmann
e_den_arr(i,j,k) = Nc*exp( -(Ec_corr - q*Ef) / (kb*T) );
hole_den_arr(i,j,k) = Nv*exp( -(q*Ef - Ev_corr) / (kb*T) );

acceptor_den_arr(i,j,k) = acceptor_doping/(1.0 + g_A*exp((-q*Ef + q*Ea + q*phi_ref - q*Chi - q*Eg - q*phi(i,j,k))/(kb*T)));
donor_den_arr(i,j,k) = donor_doping/(1.0 + g_D*exp( (q*Ef + q*Ed - q*phi_ref + q*Chi + q*phi(i,j,k)) / (kb*T) ));


if (complete_ionization == 0){
acceptor_den_arr(i,j,k) = acceptor_doping/(1.0 + g_A*exp((-q*Ef + q*Ea + q*phi_ref - q*Chi - q*Eg - q*phi(i,j,k))/(kb*T)));
donor_den_arr(i,j,k) = donor_doping/(1.0 + g_D*exp( (q*Ef + q*Ed - q*phi_ref + q*Chi + q*phi(i,j,k)) / (kb*T) ));
} else {
acceptor_den_arr(i,j,k) = acceptor_doping;
donor_den_arr(i,j,k) = donor_doping;
}
}

charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k));
Expand All @@ -99,5 +110,10 @@ void ComputeRho(MultiFab& PoissonPhi,
}
});
}
rho.FillBoundary(geom.periodicity());
e_den.FillBoundary(geom.periodicity());
p_den.FillBoundary(geom.periodicity());
acceptor_den.FillBoundary(geom.periodicity());
donor_den.FillBoundary(geom.periodicity());
}

8 changes: 4 additions & 4 deletions Source/Solver/DerivativeAlgorithm.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,8 +193,8 @@ using namespace FerroX;

} else if (P_BC_flag_lo[2] == 1){

Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda);
return (dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]); // dP/dz = P_lo/lambda;
Real F_lo = F(i,j,k)/(1 + dx[2]/2/(-lambda));
return (dx[2]*F_lo/(-lambda) - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]); // dP/dz = P_lo/lambda;

// Real F_lo = (9. * F(i,j,k) - F(i,j,k+1)) / (3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative
// return -(dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]);// dP/dz = P_lo/lambda;
Expand Down Expand Up @@ -407,8 +407,8 @@ using namespace FerroX;

} else if (P_BC_flag_lo[2] == 1){

Real F_lo = F(i,j,k)/(1 + dx[2]/2/lambda);
return (-dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];//dPdz = P_lo/lambda;
Real F_lo = F(i,j,k)/(1 + dx[2]/2/(-lambda));
return (-dx[2]*F_lo/(-lambda) - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];//dPdz = P_lo/lambda;

// Real F_lo = (9. * F(i,j,k) - F(i,j,k+1)) / (3. * dx[2] / lambda + 8.); // derived with 2nd order one-sided 1st derivative
// return (-dx[2]*F_lo/lambda - F(i,j,k) + F(i,j,k+1))/dx[2]/dx[2];// dPdz = P_lo/lambda;
Expand Down
4 changes: 3 additions & 1 deletion Source/Solver/ElectrostaticSolver.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,8 @@ void dF_dPhi(MultiFab& alpha_cc,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& acceptor_den,
MultiFab& donor_den,
MultiFab& MaterialMask,
MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta,
const Geometry& geom,
Expand All @@ -61,4 +63,4 @@ void SetPoissonBC(c_FerroX& rFerroX, std::array<std::array<amrex::LinOpBCType,AM

void Fill_Constant_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& PoissonPhi);
void Fill_FunctionBased_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time);

void SetNucleation(Array<MultiFab, AMREX_SPACEDIM> &P_old, MultiFab& MaterialMask, const amrex::GpuArray<int, AMREX_SPACEDIM>& n_cell);
58 changes: 57 additions & 1 deletion Source/Solver/ElectrostaticSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
#include <AMReX.H>
#include "ElectrostaticSolver.H"
#include "DerivativeAlgorithm.H"
#include "ChargeDensity.H"
Expand Down Expand Up @@ -87,6 +88,8 @@ void dF_dPhi(MultiFab& alpha_cc,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& acceptor_den,
MultiFab& donor_den,
MultiFab& MaterialMask,
MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta,
const Geometry& geom,
Expand All @@ -102,7 +105,7 @@ void dF_dPhi(MultiFab& alpha_cc,
PoissonPhi_plus_delta.plus(delta, 0, 1, 0);

// Calculate rho from Phi in SC region
ComputeRho(PoissonPhi_plus_delta, rho, e_den, p_den, MaterialMask);
ComputeRho(PoissonPhi_plus_delta, rho, e_den, p_den, acceptor_den, donor_den, MaterialMask, geom);

//Compute RHS of Poisson equation
ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom);
Expand Down Expand Up @@ -521,3 +524,56 @@ void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray<int, AMREX_SPACEDIM>
PoissonPhi.FillBoundary(geom.periodicity());
}


void SetNucleation(Array<MultiFab, AMREX_SPACEDIM> &P_old, MultiFab& NucleationMask, const amrex::GpuArray<int, AMREX_SPACEDIM>& n_cell)
{
int seed = random_seed;

int nprocs = ParallelDescriptor::NProcs();

if (prob_type == 1) {
amrex::InitRandom(seed , nprocs, seed ); // give all MPI ranks the same seed
} else {
amrex::InitRandom(seed+ParallelDescriptor::MyProc(), nprocs, seed+ParallelDescriptor::MyProc()); // give all MPI ranks a different seed
}

int nrand = n_cell[0]*n_cell[2];
amrex::Gpu::ManagedVector<Real> rngs(nrand, 0.0);

// generate random numbers on the host
for (int i=0; i<nrand; ++i) {
//rngs[i] = amrex::RandomNormal(0.,1.); // zero mean, unit variance
rngs[i] = amrex::Random(); // uniform [0,1] option
}

for (MFIter mfi(P_old[0]); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();

const Array4<Real> &pOld_p = P_old[0].array(mfi);
const Array4<Real> &pOld_q = P_old[1].array(mfi);
const Array4<Real> &pOld_r = P_old[2].array(mfi);
const Array4<Real>& mask = NucleationMask.array(mfi);

Real* rng = rngs.data();

amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::RandomEngine const& engine) noexcept
{
if (mask(i,j,k) == 1.) {
if (prob_type == 1) { //2D

pOld_p(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[0]*noise_amplitude;
pOld_q(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[1]*noise_amplitude;
pOld_r(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[2]*noise_amplitude;

} else if (prob_type == 2) { //3D

pOld_p(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[0]*noise_amplitude;
pOld_q(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[1]*noise_amplitude;
pOld_r(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[2]*noise_amplitude;
}
}
});
}

}
3 changes: 3 additions & 0 deletions Source/Solver/Initialization.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ void InitializePandRho(Array<MultiFab, AMREX_SPACEDIM> &P_old,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& acceptor_den,
MultiFab& donor_den,
const MultiFab& MaterialMask,
const MultiFab& tphaseMask,
const amrex::GpuArray<int, AMREX_SPACEDIM>& n_cell,
Expand All @@ -26,5 +28,6 @@ void InitializeMaterialMask(MultiFab& MaterialMask,

void InitializeMaterialMask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& MaterialMask);
void Initialize_tphase_Mask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& tphaseMask);
void Initialize_nucleation_Mask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& NucleationMask);
void Initialize_Euler_angles(c_FerroX& rFerroX, const Geometry& geom, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta);

Loading