Skip to content

Implement rho in metal #44

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

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

Implement rho in metal #44

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
3f99d74
Implement rho in metal
prkkumar Oct 23, 2023
6454080
metal rho calculation -- incomplete
prkkumar Oct 26, 2023
bdc8bef
rho in metal -- TO DO fix compilation errors
prkkumar Oct 27, 2023
bc4b8d8
Use ParallelDescriptor::ReduceRealSum
prkkumar Oct 27, 2023
4f366b1
move avergae calculation outside MFIter
prkkumar Oct 27, 2023
f54bff8
Update Source/Solver/ChargeDensity.cpp
prkkumar Oct 30, 2023
2308296
fix initialization and main for metal charge
prkkumar Nov 2, 2023
9ad4f9c
treat metal as SC
prkkumar Nov 9, 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
28 changes: 28 additions & 0 deletions Source/FerroX.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,9 @@ AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> FerroX::FE_hi;
AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> FerroX::SC_hi;
AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> FerroX::Channel_hi;
AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> FerroX::Channel_lo;
AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> FerroX::Metal_hi;
AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> FerroX::Metal_lo;


// material parameters
AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_0;
Expand All @@ -184,6 +187,9 @@ AMREX_GPU_MANAGED amrex::Real FerroX::epsilonX_fe_tphase;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilonZ_fe;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_de;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_si;
AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_metal;
AMREX_GPU_MANAGED amrex::Real FerroX::metal_screening_length;
AMREX_GPU_MANAGED amrex::Real FerroX::metal_thickness;
AMREX_GPU_MANAGED amrex::Real FerroX::alpha; // alpha = 2*alpha_1
AMREX_GPU_MANAGED amrex::Real FerroX::beta; // beta = 4*alpha_11
AMREX_GPU_MANAGED amrex::Real FerroX::gamma; // gamma = 6*alpha_111
Expand Down Expand Up @@ -284,6 +290,16 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
pp.get("epsilonZ_fe",epsilonZ_fe);// epsilon_r for FE
pp.get("epsilon_de",epsilon_de);// epsilon_r for DE
pp.get("epsilon_si",epsilon_si);// epsilon_r for SC

epsilon_metal = 1.;
pp.query("epsilon_metal",epsilon_metal);// epsilon_r for metal

metal_screening_length = 1.e5;
pp.query("metal_screening_length",metal_screening_length);// metal_screening_length

metal_thickness = 0.;
pp.query("metal_thickness",metal_thickness);// metal_thickness

pp.get("alpha",alpha);
pp.get("beta",beta);
pp.get("gamma",FerroX::gamma);
Expand Down Expand Up @@ -415,6 +431,18 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
}
}

if (pp.queryarr("Metal_lo",temp)) {
for (int i=0; i<AMREX_SPACEDIM; ++i) {
Metal_lo[i] = temp[i];
}
}

if (pp.queryarr("Metal_hi",temp)) {
for (int i=0; i<AMREX_SPACEDIM; ++i) {
Metal_hi[i] = temp[i];
}
}

if (pp.queryarr("t_phase_lo",temp)) {
for (int i=0; i<AMREX_SPACEDIM; ++i) {
t_phase_lo[i] = temp[i];
Expand Down
6 changes: 6 additions & 0 deletions Source/FerroX_namespace.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,9 @@ namespace FerroX {
extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> SC_hi;
extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> Channel_hi;
extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> Channel_lo;
extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> Metal_hi;
extern AMREX_GPU_MANAGED amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> Metal_lo;


// material parameters
extern AMREX_GPU_MANAGED amrex::Real epsilon_0;
Expand All @@ -23,6 +26,9 @@ namespace FerroX {
extern AMREX_GPU_MANAGED amrex::Real epsilonZ_fe;
extern AMREX_GPU_MANAGED amrex::Real epsilon_de;
extern AMREX_GPU_MANAGED amrex::Real epsilon_si;
extern AMREX_GPU_MANAGED amrex::Real epsilon_metal;
extern AMREX_GPU_MANAGED amrex::Real metal_screening_length;
extern AMREX_GPU_MANAGED amrex::Real metal_thickness;
extern AMREX_GPU_MANAGED amrex::Real alpha; // alpha = 2*alpha_1
extern AMREX_GPU_MANAGED amrex::Real beta; // beta = 4*alpha_11
extern AMREX_GPU_MANAGED amrex::Real gamma; // gamma = 6*alpha_111
Expand Down
12 changes: 11 additions & 1 deletion Source/Solver/ChargeDensity.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,19 @@
using namespace amrex;
using namespace FerroX;


void ComputeRho(MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM> &P_old,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
const MultiFab& MaterialMask);
const MultiFab& MaterialMask,
const Geometry& geom,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi);

void Compute_P_Sum(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum);

void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count);

void Compute_P_av(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real &P_av_z);
180 changes: 146 additions & 34 deletions Source/Solver/ChargeDensity.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,20 +2,46 @@

// Compute rho in SC region for given phi
void ComputeRho(MultiFab& PoissonPhi,
Array<MultiFab, AMREX_SPACEDIM> &P_old,
MultiFab& rho,
MultiFab& e_den,
MultiFab& p_den,
const MultiFab& MaterialMask)
const MultiFab& MaterialMask,
const Geometry& geom,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_lo,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi)
{
// Calculate average Pr. We take the average only over the FE region
Real average_P_r = 0.;
Real total_P_r = 0.;
Real FE_index_counter = 0.;

Compute_P_av(P_old, total_P_r, MaterialMask, FE_index_counter, average_P_r);

//Calculate integrated electrode charge (Qe) based on eq 13 of https://pubs.aip.org/aip/jap/article/44/8/3379/6486/Depolarization-fields-in-thin-ferroelectric-films
Real FE_thickness = FE_hi[2] - FE_lo[2];
Real coth = (exp(2.*metal_thickness/metal_screening_length) + 1.0) / (exp(2.*metal_thickness/metal_screening_length) - 1.0);
Real csch = (2.*exp(metal_thickness/metal_screening_length)) / (exp(2.*metal_thickness/metal_screening_length) - 1.0);
Real numerator = 0.5 * FE_thickness * average_P_r / epsilonX_fe;
Real denominator = metal_screening_length/epsilon_de*(coth - csch + FE_thickness/(2.*epsilonX_fe));
Copy link
Collaborator

@jackieyao0114 jackieyao0114 Oct 30, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Real denominator = metal_screening_length/epsilon_de*(coth - csch + FE_thickness/(2.*epsilonX_fe));
Real denominator = metal_screening_length/epsilon_metal*(coth - csch + FE_thickness/(2.*epsilonX_fe));

Real Qe = -numerator/denominator;

// 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);

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>& hole_den_arr = p_den.array(mfi);
const Array4<Real>& e_den_arr = e_den.array(mfi);
const Array4<Real>& charge_den_arr = rho.array(mfi);
Expand All @@ -24,45 +50,60 @@ void ComputeRho(MultiFab& PoissonPhi,
const Array4<Real>& donor_den_arr = donor_den.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) {

if(use_Fermi_Dirac == 1){
//Fermi-Dirac
Real eta_n = q*(phi(i,j,k) - Ec)/(kb*T);
Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2)));
Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8));
Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0);
Real z = prob_lo[2] + (k+0.5) * dx[2];
Real z_metal = 0.;

e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n;

Real eta_p = q*(Ev - phi(i,j,k))/(kb*T);
Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2)));
Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8));
Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0);
if (mask(i,j,k) >= 2.0) {

hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p;
if (mask(i,j,k) == 4.0) { //Metal

if(z <= FE_lo[2]){
z_metal = std::abs(FE_lo[2] - (k + 0.5) * dx[2]);
} else if (z >= FE_hi[2]){
z_metal = std::abs((k + 0.5) * dx[2] - FE_hi[2]);
}
charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length);

} else {
//Maxwell-Boltzmann
Real n_0 = intrinsic_carrier_concentration;
Real p_0 = intrinsic_carrier_concentration;
hole_den_arr(i,j,k) = n_0*exp(-(q*phi(i,j,k))/(kb*T));
e_den_arr(i,j,k) = p_0*exp(q*phi(i,j,k)/(kb*T));
}

//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));

if(use_Fermi_Dirac == 1){
//Fermi-Dirac
Real eta_n = q*(phi(i,j,k) - Ec)/(kb*T);
Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2)));
Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8));
Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0);

e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n;

Real eta_p = q*(Ev - phi(i,j,k))/(kb*T);
Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2)));
Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8));
Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0);

hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p;
} else {
//Maxwell-Boltzmann
Real n_0 = intrinsic_carrier_concentration;
Real p_0 = intrinsic_carrier_concentration;
hole_den_arr(i,j,k) = n_0*exp(-(q*phi(i,j,k))/(kb*T));
e_den_arr(i,j,k) = p_0*exp(q*phi(i,j,k)/(kb*T));
}

//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 if (mask(i,j,k) == 2.0){ // 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 {

Expand All @@ -73,3 +114,74 @@ void ComputeRho(MultiFab& PoissonPhi,
}
}

void Compute_P_Sum(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum)
{

// Initialize to zero
sum = 0.;

ReduceOps<ReduceOpSum> reduce_op;

ReduceData<Real> reduce_data(reduce_op);
using ReduceTuple = typename decltype(reduce_data)::Type;

for (MFIter mfi(P[2],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Box& bx_grid = mfi.validbox();

auto const& fab = P[2].array(mfi);

reduce_op.eval(bx, reduce_data,
[=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
{
return {fab(i,j,k)};
});
}

sum = amrex::get<0>(reduce_data.value());
ParallelDescriptor::ReduceRealSum(sum);
}


void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count)
{

// Initialize to zero
count = 0.;

ReduceOps<ReduceOpSum> reduce_op;

ReduceData<Real> reduce_data(reduce_op);
using ReduceTuple = typename decltype(reduce_data)::Type;

for (MFIter mfi(MaterialMask, TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.tilebox();
const Box& bx_grid = mfi.validbox();

auto const& fab = MaterialMask.array(mfi);

reduce_op.eval(bx, reduce_data,
[=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
{
if(fab(i,j,k) == 0.) {
return {1.};
} else {
return {0.};
}

});
}

count = amrex::get<0>(reduce_data.value());
ParallelDescriptor::ReduceRealSum(count);
}

void Compute_P_av(const std::array<MultiFab, AMREX_SPACEDIM>& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real& P_av_z)
{
Compute_P_Sum(P, sum);
Compute_P_index_Sum(MaterialMask, count);
P_av_z = sum/count;
}

8 changes: 5 additions & 3 deletions Source/Solver/ElectrostaticSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,8 +102,8 @@ 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, P_old, rho, e_den, p_den, MaterialMask, geom, prob_lo, prob_hi);
//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 @@ -231,7 +231,7 @@ void InitializePermittivity(std::array<std::array<amrex::LinOpBCType,AMREX_SPACE

Real x = prob_lo[0] + (i+0.5) * dx[0];
Real y = prob_lo[1] + (j+0.5) * dx[1];
Real z = prob_lo[1] + (k+0.5) * dx[2];
Real z = prob_lo[2] + (k+0.5) * dx[2];

if(mask(i,j,k) == 0.0) {
beta(i,j,k) = epsilonX_fe * epsilon_0; //FE layer
Expand All @@ -244,6 +244,8 @@ void InitializePermittivity(std::array<std::array<amrex::LinOpBCType,AMREX_SPACE
beta(i,j,k) = epsilon_de * epsilon_0; //DE layer
} else if (mask(i,j,k) >= 2.0){
beta(i,j,k) = epsilon_si * epsilon_0; //SC layer
} else if (mask(i,j,k) == 4.0){
beta(i,j,k) = epsilon_metal * epsilon_0; //Metal layer
} else {
beta(i,j,k) = epsilon_de * epsilon_0; //Spacer is same as DE
}
Expand Down
2 changes: 2 additions & 0 deletions Source/Solver/Initialization.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -211,6 +211,8 @@ void InitializeMaterialMask(MultiFab& MaterialMask,
if (x <= Channel_hi[0] && x >= Channel_lo[0] && y <= Channel_hi[1] && y >= Channel_lo[1] && z <= Channel_hi[2] && z >= Channel_lo[2]){
mask(i,j,k) = 3.;
}
} else if (x <= Metal_hi[0] && x >= Metal_lo[0] && y <= Metal_hi[1] && y >= Metal_lo[1] && z <= Metal_hi[2] && z >= Metal_lo[2]) {
mask(i,j,k) = 4.;
} else {
mask(i,j,k) = 1.; //spacer is DE
}
Expand Down
12 changes: 6 additions & 6 deletions Source/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,8 +144,8 @@ void main_main (c_FerroX& rFerroX)
angle_theta.setVal(0.);

//Initialize material mask
InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi);
//InitializeMaterialMask(rFerroX, geom, MaterialMask);
//InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi);
InitializeMaterialMask(rFerroX, geom, MaterialMask);
if(Coordinate_Transformation == 1){
Initialize_tphase_Mask(rFerroX, geom, tphaseMask);
Initialize_Euler_angles(rFerroX, geom, angle_alpha, angle_beta, angle_theta);
Expand Down Expand Up @@ -302,7 +302,7 @@ void main_main (c_FerroX& rFerroX)
PoissonPhi.FillBoundary(geom.periodicity());

// Calculate rho from Phi in SC region
ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask);
ComputeRho(PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi);

if (contains_SC == 0) {
// no semiconductor region; set error to zero so the while loop terminates
Expand Down Expand Up @@ -423,7 +423,7 @@ void main_main (c_FerroX& rFerroX)
PoissonPhi.FillBoundary(geom.periodicity());

// Calculate rho from Phi in SC region
ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask);
ComputeRho(PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi);

if (contains_SC == 0) {
// no semiconductor region; set error to zero so the while loop terminates
Expand Down Expand Up @@ -493,7 +493,7 @@ void main_main (c_FerroX& rFerroX)
PoissonPhi.FillBoundary(geom.periodicity());

// Calculate rho from Phi in SC region
ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask);
ComputeRho(PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi);

if (contains_SC == 0) {
// no semiconductor region; set error to zero so the while loop terminates
Expand Down Expand Up @@ -652,7 +652,7 @@ void main_main (c_FerroX& rFerroX)
PoissonPhi.FillBoundary(geom.periodicity());

// Calculate rho from Phi in SC region
ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask);
ComputeRho(PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi);

if (contains_SC == 0) {
// no semiconductor region; set error to zero so the while loop terminates
Expand Down