Skip to content
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

drift-diffusion model for carrier transport #38

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

drift-diffusion model for carrier transport #38

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
feade4f
drift-diffusion model for carrier transport
prkkumar Aug 2, 2023
be91814
fix time integration of continuity equation
prkkumar Aug 2, 2023
00aeeba
update newton's method dF_dPhi for drift-diffusion
prkkumar Aug 9, 2023
221ef30
update newtons method everywhere
prkkumar Aug 9, 2023
612beeb
add dFdPhi function with drift-diffusion
prkkumar Aug 11, 2023
32130bb
DD fuxes
prkkumar Aug 11, 2023
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
16 changes: 16 additions & 0 deletions Source/FerroX.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,12 @@ AMREX_GPU_MANAGED amrex::Real FerroX::T;
AMREX_GPU_MANAGED amrex::Real FerroX::acceptor_doping;
AMREX_GPU_MANAGED amrex::Real FerroX::donor_doping;

//Drift-Diffusion
AMREX_GPU_MANAGED amrex::Real FerroX::electron_mobility;
AMREX_GPU_MANAGED amrex::Real FerroX::hole_mobility;
AMREX_GPU_MANAGED amrex::Real FerroX::electron_affinity;
AMREX_GPU_MANAGED amrex::Real FerroX::bandgap;

// P and Phi Bc
AMREX_GPU_MANAGED amrex::Real FerroX::lambda;
AMREX_GPU_MANAGED amrex::GpuArray<int, AMREX_SPACEDIM> FerroX::P_BC_flag_lo;
Expand Down Expand Up @@ -440,6 +446,16 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
donor_doping = 0.0;
pp.query("donor_doping",donor_doping);

//Drift-Diffusion
electron_mobility = 0.0;
pp.query("electron_mobility",electron_mobility);
hole_mobility = 0.0;
pp.query("hole_mobility",hole_mobility);
electron_affinity = 0.0;
pp.query("electron_affinity",electron_affinity);
bandgap = 0.0;
pp.query("bandgap",bandgap);

Coordinate_Transformation = 0;
pp.query("Coordinate_Transformation",Coordinate_Transformation);

Expand Down
4 changes: 4 additions & 0 deletions Source/FerroX_namespace.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,10 @@ namespace FerroX {
extern AMREX_GPU_MANAGED amrex::Real T;
extern AMREX_GPU_MANAGED amrex::Real acceptor_doping;
extern AMREX_GPU_MANAGED amrex::Real donor_doping;
extern AMREX_GPU_MANAGED amrex::Real electron_mobility;
extern AMREX_GPU_MANAGED amrex::Real hole_mobility;
extern AMREX_GPU_MANAGED amrex::Real electron_affinity;
extern AMREX_GPU_MANAGED amrex::Real bandgap;

// P and Phi Bc
extern AMREX_GPU_MANAGED amrex::Real lambda;
Expand Down
10 changes: 10 additions & 0 deletions Source/Solver/ChargeDensity.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,3 +12,13 @@ void ComputeRho(MultiFab& PoissonPhi,
MultiFab& p_den,
const MultiFab& MaterialMask);

void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM> &Jn,
Array<MultiFab, AMREX_SPACEDIM> &Jp,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& e_den_old,
MultiFab& p_den_old,
const Geometry& geom,
const MultiFab& MaterialMask);
175 changes: 175 additions & 0 deletions Source/Solver/ChargeDensity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include "ChargeDensity.H"
#include "DerivativeAlgorithm.H"

// Compute rho in SC region for given phi
void ComputeRho(MultiFab& PoissonPhi,
Expand Down Expand Up @@ -70,3 +71,177 @@ void ComputeRho(MultiFab& PoissonPhi,
}
}

// Drift-Diffusion : Compute rho in SC region for given phi, Jn, and Jp
void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM> &Jn,
Array<MultiFab, AMREX_SPACEDIM> &Jp,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& e_den_old,
MultiFab& p_den_old,
const Geometry& geom,
const MultiFab& MaterialMask)
{
MultiFab Ec_mf(rho.boxArray(), rho.DistributionMap(), 1, 0);
MultiFab Ev_mf(rho.boxArray(), rho.DistributionMap(), 1, 0);

//Calculate Ec and Ev
for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.validbox();

const Array4<Real>& phi = PoissonPhi.array(mfi);
const Array4<Real>& Ec_arr = Ec_mf.array(mfi);
const Array4<Real>& Ev_arr = Ev_mf.array(mfi);
const Array4<Real const>& mask = MaterialMask.array(mfi);

amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{

if (mask(i,j,k) >= 2.0) {
//electrons
Ec_arr(i,j,k) = -q*phi(i,j,k) - electron_affinity;
Ev_arr(i,j,k) = Ec_arr(i,j,k) - bandgap;
} else {
Ec_arr(i,j,k) = 0.0;
Ev_arr(i,j,k) = 0.0;
}
});
}

// Calculate currents Jn and Jp
for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.validbox();

// extract dx from the geometry object
GpuArray<Real,AMREX_SPACEDIM> dx = geom.CellSizeArray();

// Calculate charge density from Phi, Nc, Nv, Ec, and Ev

const Array4<Real>& hole_den_arr = p_den.array(mfi);
const Array4<Real>& e_den_arr = e_den.array(mfi);
const Array4<Real>& phi = PoissonPhi.array(mfi);
const Array4<Real const>& mask = MaterialMask.array(mfi);

const Array4<Real>& Ec_arr = Ec_mf.array(mfi);
const Array4<Real>& Ev_arr = Ev_mf.array(mfi);
const Array4<Real>& Jnx_arr = Jn[0].array(mfi);
const Array4<Real>& Jny_arr = Jn[1].array(mfi);
const Array4<Real>& Jnz_arr = Jn[2].array(mfi);

const Array4<Real>& Jpx_arr = Jp[0].array(mfi);
const Array4<Real>& Jpy_arr = Jp[1].array(mfi);
const Array4<Real>& Jpz_arr = Jp[2].array(mfi);

amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{

if (mask(i,j,k) >= 2.0) {

//electrons
Real gradEc_x = DFDx(Ec_arr, i, j, k, dx);
Real gradEc_y = DFDy(Ec_arr, i, j, k, dx);
Real gradEc_z = DFDz(Ec_arr, i, j, k, dx);

Real gradn_x = DFDx(e_den_arr, i, j, k, dx);
Real gradn_y = DFDy(e_den_arr, i, j, k, dx);
Real gradn_z = DFDz(e_den_arr, i, j, k, dx);

Jnx_arr(i,j,k) = electron_mobility*e_den_arr(i,j,k)*gradEc_x + kb*T*electron_mobility*gradn_x;
Jny_arr(i,j,k) = electron_mobility*e_den_arr(i,j,k)*gradEc_y + kb*T*electron_mobility*gradn_y;
Jnz_arr(i,j,k) = electron_mobility*e_den_arr(i,j,k)*gradEc_z + kb*T*electron_mobility*gradn_z;

//holes
Real gradEv_x = DFDx(Ev_arr, i, j, k, dx);
Real gradEv_y = DFDy(Ev_arr, i, j, k, dx);
Real gradEv_z = DFDz(Ev_arr, i, j, k, dx);

Real gradp_x = DFDx(hole_den_arr, i, j, k, dx);
Real gradp_y = DFDy(hole_den_arr, i, j, k, dx);
Real gradp_z = DFDz(hole_den_arr, i, j, k, dx);

Jpx_arr(i,j,k) = hole_mobility*hole_den_arr(i,j,k)*gradEv_x - kb*T*hole_mobility*gradp_x;
Jpy_arr(i,j,k) = hole_mobility*hole_den_arr(i,j,k)*gradEv_y - kb*T*hole_mobility*gradp_y;
Jpz_arr(i,j,k) = hole_mobility*hole_den_arr(i,j,k)*gradEv_z - kb*T*hole_mobility*gradp_z;

} else {

Jnx_arr(i,j,k) = 0.0;
Jny_arr(i,j,k) = 0.0;
Jnz_arr(i,j,k) = 0.0;

Jpx_arr(i,j,k) = 0.0;
Jpy_arr(i,j,k) = 0.0;
Jpz_arr(i,j,k) = 0.0;

}
});
}

// loop over boxes
for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.validbox();

// extract dx from the geometry object
GpuArray<Real,AMREX_SPACEDIM> dx = geom.CellSizeArray();

// Calculate charge density from Phi, Nc, Nv, Ec, and Ev
MultiFab acceptor_den(rho.boxArray(), rho.DistributionMap(), 1, 0);
MultiFab donor_den(rho.boxArray(), rho.DistributionMap(), 1, 0);
MultiFab div_Jn(rho.boxArray(), rho.DistributionMap(), 1, 0);
MultiFab div_Jp(rho.boxArray(), rho.DistributionMap(), 1, 0);

const Array4<Real>& hole_den_arr = p_den.array(mfi);
const Array4<Real>& e_den_arr = e_den.array(mfi);
const Array4<Real>& hole_den_old_arr = p_den.array(mfi);
const Array4<Real>& e_den_old_arr = e_den.array(mfi);
const Array4<Real>& charge_den_arr = rho.array(mfi);
const Array4<Real>& phi = PoissonPhi.array(mfi);
const Array4<Real>& acceptor_den_arr = acceptor_den.array(mfi);
const Array4<Real>& donor_den_arr = donor_den.array(mfi);
const Array4<Real const>& mask = MaterialMask.array(mfi);

const Array4<Real>& Jnx_arr = Jn[0].array(mfi);
const Array4<Real>& Jny_arr = Jn[1].array(mfi);
const Array4<Real>& Jnz_arr = Jn[2].array(mfi);
const Array4<Real>& div_Jn_arr = div_Jn.array(mfi);

const Array4<Real>& Jpx_arr = Jp[0].array(mfi);
const Array4<Real>& Jpy_arr = Jp[1].array(mfi);
const Array4<Real>& Jpz_arr = Jp[2].array(mfi);
const Array4<Real>& div_Jp_arr = div_Jp.array(mfi);

amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{

if (mask(i,j,k) >= 2.0) {

div_Jn_arr(i,j,k) = DFDx(Jnx_arr, i,j,k,dx) + DFDy(Jny_arr, i,j,k,dx) + DFDz(Jnz_arr, i,j,k,dx);
div_Jp_arr(i,j,k) = DFDx(Jpx_arr, i,j,k,dx) + DFDy(Jpy_arr, i,j,k,dx) + DFDz(Jpz_arr, i,j,k,dx);

e_den_arr(i,j,k) = e_den_old_arr(i,j,k) + dt*div_Jn_arr(i,j,k);
hole_den_arr(i,j,k) = hole_den_old_arr(i,j,k) - dt*div_Jp_arr(i,j,k);

//If in channel, set acceptor doping, else (Source/Drain) set donor doping
if (mask(i,j,k) == 3.0) {
acceptor_den_arr(i,j,k) = acceptor_doping;
donor_den_arr(i,j,k) = 0.0;
} else { // Source / Drain
acceptor_den_arr(i,j,k) = 0.0;
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));

} else {

charge_den_arr(i,j,k) = 0.0;

}
});
}
}

9 changes: 9 additions & 0 deletions Source/Solver/DerivativeAlgorithm.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,15 @@ using namespace FerroX;
return (F(i,j+1,k) - F(i,j-1,k))/(2.*dx[1]);
}

/**
* Perform first derivative dF/dz */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
static amrex::Real DFDz (
amrex::Array4<amrex::Real> const& F,
int const i, int const j, int const k, amrex::GpuArray<amrex::Real, 3> dx) {
return (F(i,j,k+1) - F(i,j,k-1))/(2.*dx[2]);
}

/**
* Perform first derivative dP/dx */
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand Down
28 changes: 23 additions & 5 deletions Source/Solver/ElectrostaticSolver.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,18 +39,36 @@ void InitializePermittivity(std::array<std::array<amrex::LinOpBCType,AMREX_SPACE
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi);

void dF_dPhi(MultiFab& alpha_cc,
MultiFab& PoissonRHS,
MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM>& P_old,
MultiFab& PoissonRHS,
MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM>& P_old,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& MaterialMask,
MultiFab& MaterialMask,
MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta,
const Geometry& geom,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi);

void dF_dPhi_DD(MultiFab& alpha_cc,
MultiFab& PoissonRHS,
MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM>& P_old,
Array<MultiFab, AMREX_SPACEDIM>& Jn,
Array<MultiFab, AMREX_SPACEDIM>& Jp,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& e_den_old,
MultiFab& p_den_old,
MultiFab& MaterialMask,
MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta,
const Geometry& geom,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi);


void ComputePoissonRHS_Newton(MultiFab& PoissonRHS,
MultiFab& PoissonPhi,
MultiFab& alpha_cc);
Expand Down
38 changes: 37 additions & 1 deletion Source/Solver/ElectrostaticSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,13 +102,49 @@ void dF_dPhi(MultiFab& alpha_cc,
PoissonPhi_plus_delta.plus(delta, 0, 1, 0);

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

//Compute RHS of Poisson equation
ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom);

MultiFab::LinComb(alpha_cc, 1./delta, PoissonRHS_phi_plus_delta, 0, -1./delta, PoissonRHS, 0, 0, 1, 0);
}

void dF_dPhi_DD(MultiFab& alpha_cc,
MultiFab& PoissonRHS,
MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM>& P_old,
Array<MultiFab, AMREX_SPACEDIM>& Jn,
Array<MultiFab, AMREX_SPACEDIM>& Jp,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
MultiFab& e_den_old,
MultiFab& p_den_old,
MultiFab& MaterialMask,
MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta,
const Geometry& geom,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi)

{

MultiFab PoissonPhi_plus_delta(PoissonPhi.boxArray(), PoissonPhi.DistributionMap(), 1, 0);
MultiFab PoissonRHS_phi_plus_delta(PoissonRHS.boxArray(), PoissonRHS.DistributionMap(), 1, 0);

MultiFab::Copy(PoissonPhi_plus_delta, PoissonPhi, 0, 0, 1, 0);
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_DriftDiffusion(PoissonPhi_plus_delta, Jn, Jp, rho, e_den, p_den, e_den_old, p_den_old, geom, MaterialMask);

//Compute RHS of Poisson equation
ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom);

MultiFab::LinComb(alpha_cc, 1./delta, PoissonRHS_phi_plus_delta, 0, -1./delta, PoissonRHS, 0, 0, 1, 0);
}

void ComputePoissonRHS_Newton(MultiFab& PoissonRHS,
MultiFab& PoissonPhi,
MultiFab& alpha_cc)
Expand Down
Loading