From 3f99d74098bf8ed88829bc8d033afaf35b5bf8af Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Mon, 23 Oct 2023 11:50:44 -0700 Subject: [PATCH 1/8] Implement rho in metal --- Source/FerroX.cpp | 24 +++++++++ Source/FerroX_namespace.H | 5 ++ Source/Solver/ChargeDensity.cpp | 76 ++++++++++++++++----------- Source/Solver/ElectrostaticSolver.cpp | 2 + Source/Solver/Initialization.cpp | 2 + 5 files changed, 77 insertions(+), 32 deletions(-) diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 680074f..ba4ba07 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -176,6 +176,9 @@ AMREX_GPU_MANAGED amrex::GpuArray FerroX::FE_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::SC_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Channel_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Channel_lo; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::Metal_hi; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::Metal_lo; + // material parameters AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_0; @@ -184,6 +187,8 @@ 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::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 @@ -284,6 +289,13 @@ void InitializeFerroXNamespace(const amrex::GpuArray SC_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Channel_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Channel_lo; + extern AMREX_GPU_MANAGED amrex::GpuArray Metal_hi; + extern AMREX_GPU_MANAGED amrex::GpuArray Metal_lo; + // material parameters extern AMREX_GPU_MANAGED amrex::Real epsilon_0; @@ -23,6 +26,8 @@ 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 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 diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index 7c93c24..66335ca 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -2,6 +2,7 @@ // Compute rho in SC region for given phi void ComputeRho(MultiFab& PoissonPhi, + Array &P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, @@ -16,6 +17,10 @@ void ComputeRho(MultiFab& PoissonPhi, MultiFab acceptor_den(rho.boxArray(), rho.DistributionMap(), 1, 0); MultiFab donor_den(rho.boxArray(), rho.DistributionMap(), 1, 0); + const Array4 &pOld_p = P_old[0].array(mfi); + const Array4 &pOld_q = P_old[1].array(mfi); + const Array4 &pOld_r = P_old[2].array(mfi); + const Array4& hole_den_arr = p_den.array(mfi); const Array4& e_den_arr = e_den.array(mfi); const Array4& charge_den_arr = rho.array(mfi); @@ -30,39 +35,46 @@ void ComputeRho(MultiFab& PoissonPhi, 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); - - 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; + if (mask(i,j,k) == 4.0) { //Metal + + charge_den_arr(i,j,k) = //TODO; + } 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 { diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index a1cac53..8c9e9f9 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -244,6 +244,8 @@ void InitializePermittivity(std::array= 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 } diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 82cf1a0..2f82452 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -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 } From 6454080959336ab05c8eb10a4091c15039ef16e4 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Thu, 26 Oct 2023 14:00:36 -0700 Subject: [PATCH 2/8] metal rho calculation -- incomplete --- Source/Solver/ChargeDensity.cpp | 34 ++++++++++++++++++++++++++++++++- Source/main.cpp | 4 ++-- 2 files changed, 35 insertions(+), 3 deletions(-) diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index 66335ca..f565660 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -29,6 +29,36 @@ void ComputeRho(MultiFab& PoissonPhi, const Array4& donor_den_arr = donor_den.array(mfi); const Array4& mask = MaterialMask.array(mfi); + // Calculate average P + //Calculate Qe based on eq 13 + //rhs = Qe/screening_length*exp(-z/screening_lenth) + + + // Calculate average Pr. We take the average only over the FE region + Real average_P_r = 0.; + Real total_P_r = 0.; + int FE_index_counter = 0; + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + if (mask(i,j,k) == 0.0){ + total_P_r += pOld_r(i,j,k); + FE_index_counter += 1; + } + }); + + average_P_r = total_P_r/FE_index_counter; + + //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 numerator = 0.5 * FE_thickness * average_P_r / epsilonX_fe; + Real denominator = metal_screening_length/ + Real Qe = + + //amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + //{ + // Real coth = + // Real csch + //}); amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { @@ -37,7 +67,9 @@ void ComputeRho(MultiFab& PoissonPhi, if (mask(i,j,k) == 4.0) { //Metal - charge_den_arr(i,j,k) = //TODO; + Real z_metal_fe_lo + Real z_metal_fe_lo + charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z/metal_screening_length); } else { diff --git a/Source/main.cpp b/Source/main.cpp index 6cbcf9d..d14b13f 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -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); From bdc8bef71ddc10d1e3e0de7b6978160cd4081977 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Thu, 26 Oct 2023 17:03:47 -0700 Subject: [PATCH 3/8] rho in metal -- TO DO fix compilation errors --- Source/FerroX.cpp | 4 +++ Source/FerroX_namespace.H | 1 + Source/Solver/ChargeDensity.H | 8 +++-- Source/Solver/ChargeDensity.cpp | 50 +++++++++++++++------------ Source/Solver/ElectrostaticSolver.cpp | 6 ++-- Source/main.cpp | 8 ++--- 6 files changed, 45 insertions(+), 32 deletions(-) diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index ba4ba07..d8c435d 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -189,6 +189,7 @@ 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 @@ -296,6 +297,9 @@ void InitializeFerroXNamespace(const amrex::GpuArray &P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - const MultiFab& MaterialMask); - + const MultiFab& MaterialMask, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi); diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index f565660..a099ad2 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -6,13 +6,19 @@ void ComputeRho(MultiFab& PoissonPhi, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - const MultiFab& MaterialMask) + const MultiFab& MaterialMask, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi) { // loop over boxes for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); + // extract dx from the geometry object + GpuArray 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); @@ -29,47 +35,45 @@ void ComputeRho(MultiFab& PoissonPhi, const Array4& donor_den_arr = donor_den.array(mfi); const Array4& mask = MaterialMask.array(mfi); - // Calculate average P - //Calculate Qe based on eq 13 - //rhs = Qe/screening_length*exp(-z/screening_lenth) - - // Calculate average Pr. We take the average only over the FE region Real average_P_r = 0.; Real total_P_r = 0.; int FE_index_counter = 0; - amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (mask(i,j,k) == 0.0){ - total_P_r += pOld_r(i,j,k); - FE_index_counter += 1; + if (mask(i, j, k) == 0.0) { + total_P_r += pOld_r(i, j, k); + FE_index_counter += 1; } }); - - average_P_r = total_P_r/FE_index_counter; + + average_P_r = total_P_r / static_cast(FE_index_counter); //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 z_metal = 0.; + 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/ - Real Qe = + Real denominator = metal_screening_length/epsilon_de*(coth - csch + FE_thickness/(2.*epsilonX_fe)); + Real Qe = -numerator/denominator; - //amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - //{ - // Real coth = - // Real csch - //}); - amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real z = prob_lo[2] + (k+0.5) * dx[2]; + if (mask(i,j,k) >= 2.0) { if (mask(i,j,k) == 4.0) { //Metal - Real z_metal_fe_lo - Real z_metal_fe_lo - charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z/metal_screening_length); + 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 { diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 8c9e9f9..d4f46b8 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -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); @@ -231,7 +231,7 @@ void InitializePermittivity(std::array Date: Fri, 27 Oct 2023 14:42:23 -0700 Subject: [PATCH 4/8] Use ParallelDescriptor::ReduceRealSum --- Source/Solver/ChargeDensity.H | 6 +++ Source/Solver/ChargeDensity.cpp | 85 ++++++++++++++++++++++++++++----- 2 files changed, 80 insertions(+), 11 deletions(-) diff --git a/Source/Solver/ChargeDensity.H b/Source/Solver/ChargeDensity.H index 8b03aa5..ab29aec 100644 --- a/Source/Solver/ChargeDensity.H +++ b/Source/Solver/ChargeDensity.H @@ -16,3 +16,9 @@ void ComputeRho(MultiFab& PoissonPhi, const Geometry& geom, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); + +void Compute_P_Sum(const std::array& P, Real& sum); + +void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count); + +void Compute_P_av(const std::array& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real &P_av_z); diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index a099ad2..52f6bce 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -38,21 +38,12 @@ void ComputeRho(MultiFab& PoissonPhi, // Calculate average Pr. We take the average only over the FE region Real average_P_r = 0.; Real total_P_r = 0.; - int FE_index_counter = 0; - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (mask(i, j, k) == 0.0) { - total_P_r += pOld_r(i, j, k); - FE_index_counter += 1; - } - }); + Real FE_index_counter = 0.; - average_P_r = total_P_r / static_cast(FE_index_counter); + 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 z_metal = 0.; 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; @@ -63,6 +54,7 @@ void ComputeRho(MultiFab& PoissonPhi, { Real z = prob_lo[2] + (k+0.5) * dx[2]; + Real z_metal = 0.; if (mask(i,j,k) >= 2.0) { @@ -121,3 +113,74 @@ void ComputeRho(MultiFab& PoissonPhi, } } +void Compute_P_Sum(const std::array& P, Real& sum) +{ + + // Initialize to zero + sum = 0.; + + ReduceOps reduce_op; + + ReduceData 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 reduce_op; + + ReduceData 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& 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; +} + From 4f366b1f2f8506bb31c9efdb4e81b343e2c64ec9 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Fri, 27 Oct 2023 15:22:33 -0700 Subject: [PATCH 5/8] move avergae calculation outside MFIter --- Source/Solver/ChargeDensity.cpp | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index 52f6bce..d332f45 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -11,6 +11,21 @@ void ComputeRho(MultiFab& PoissonPhi, const amrex::GpuArray& prob_lo, const amrex::GpuArray& 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)); + Real Qe = -numerator/denominator; + // loop over boxes for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) { @@ -35,20 +50,6 @@ void ComputeRho(MultiFab& PoissonPhi, const Array4& donor_den_arr = donor_den.array(mfi); const Array4& mask = MaterialMask.array(mfi); - // 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)); - Real Qe = -numerator/denominator; amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { From f54bff8dafb4dfdec4d8aba412a44c984d436295 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar <89051199+prkkumar@users.noreply.github.com> Date: Mon, 30 Oct 2023 15:42:00 -0700 Subject: [PATCH 6/8] Update Source/Solver/ChargeDensity.cpp Co-authored-by: Zhi (Jackie) Yao <58234082+jackieyao0114@users.noreply.github.com> --- Source/Solver/ChargeDensity.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index d332f45..7139f11 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -23,7 +23,7 @@ void ComputeRho(MultiFab& PoissonPhi, 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)); + Real denominator = metal_screening_length/epsilon_de*(coth - csch) + FE_thickness/(2.*epsilonX_fe); Real Qe = -numerator/denominator; // loop over boxes From 230829617f0596f9062796b8e0ef996312760f31 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Thu, 2 Nov 2023 11:22:03 -0700 Subject: [PATCH 7/8] fix initialization and main for metal charge --- Source/Solver/ChargeDensity.cpp | 6 ++ Source/Solver/ElectrostaticSolver.cpp | 2 +- Source/Solver/Initialization.cpp | 93 +++++++++++++++++-------- Source/Utils/FerroXUtils/FerroXUtil.cpp | 2 +- Source/main.cpp | 16 ++--- 5 files changed, 80 insertions(+), 39 deletions(-) diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index 7139f11..12de4d1 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -25,6 +25,11 @@ void ComputeRho(MultiFab& PoissonPhi, 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); Real Qe = -numerator/denominator; + + amrex::Print() << "ChargeDen average_P_r = " << average_P_r << "\n"; + amrex::Print() << "ChargeDen numerator = " << numerator << "\n"; + amrex::Print() << "ChargeDen denominator = " << denominator << "\n"; + amrex::Print() << "ChargeDen Qe = " << Qe << "\n"; // loop over boxes for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) @@ -66,6 +71,7 @@ void ComputeRho(MultiFab& PoissonPhi, } else if (z >= FE_hi[2]){ z_metal = std::abs((k + 0.5) * dx[2] - FE_hi[2]); } + // amrex::Print() << "Qe = " << Qe << "\n"; charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); } else { diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index d4f46b8..f903d9f 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -242,7 +242,7 @@ void InitializePermittivity(std::array= 2.0){ + } else if (mask(i,j,k) == 2.0 || mask(i,j,k) == 3.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 diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 2f82452..bdeb288 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -1,4 +1,5 @@ #include "Initialization.H" +#include "ChargeDensity.H" #include "Utils/eXstaticUtils/eXstaticUtil.H" #include "../../Utils/SelectWarpXUtils/WarpXUtil.H" @@ -59,6 +60,25 @@ void InitializePandRho(Array &P_old, rngs[i] = amrex::Random(); // uniform [0,1] option } + Real average_P_r = 1.e-6; + //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)); + Real Qe = -numerator/denominator; + + amrex::Print() << "average_P_r = " << average_P_r << "\n"; + amrex::Print() << "numerator = " << numerator << "\n"; + amrex::Print() << "denominator = " << denominator << "\n"; + amrex::Print() << "Qe = " << Qe << "\n"; + // loop over boxes for (MFIter mfi(rho); mfi.isValid(); ++mfi) { @@ -131,45 +151,60 @@ void InitializePandRho(Array &P_old, amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real z = prob_lo[2] + (k+0.5) * dx[2]; + Real z_metal = 0.; + //SC region if (mask(i,j,k) >= 2.0) { - if(use_Fermi_Dirac == 1){ - - //Approximate FD integral - Real Phi = 0.5*(Ec + Ev); //eV - Real eta_n = q*(Phi - 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.0))); - 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; + if (mask(i,j,k) == 4.0) { //Metal - Real eta_p = q*(Ev - Phi)/(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.0))); - 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(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]); + } + //if(Qe != 0.) amrex::Print() << "initialization : Qe = " << Qe << "\n"; + charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); - hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; - } else { + if(use_Fermi_Dirac == 1){ + + //Approximate FD integral + Real Phi = 0.5*(Ec + Ev); //eV + Real eta_n = q*(Phi - 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.0))); + 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; - hole_den_arr(i,j,k) = intrinsic_carrier_concentration; - e_den_arr(i,j,k) = intrinsic_carrier_concentration; + Real eta_p = q*(Ev - Phi)/(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.0))); + 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 { + + hole_den_arr(i,j,k) = intrinsic_carrier_concentration; + e_den_arr(i,j,k) = intrinsic_carrier_concentration; + + } + + //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 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)); - + }); } for (int i = 0; i < 3; i++){ diff --git a/Source/Utils/FerroXUtils/FerroXUtil.cpp b/Source/Utils/FerroXUtils/FerroXUtil.cpp index 0237195..32e614c 100644 --- a/Source/Utils/FerroXUtils/FerroXUtil.cpp +++ b/Source/Utils/FerroXUtils/FerroXUtil.cpp @@ -25,7 +25,7 @@ void FerroX_Util::Contains_sc(MultiFab& MaterialMask, bool& contains_SC) for (auto k = lo.z; k <= hi.z; ++k) { for (auto j = lo.y; j <= hi.y; ++j) { for (auto i = lo.x; i <= hi.x; ++i) { - if (mask(i,j,k) >= 2.0) { + if (mask(i,j,k) == 2.0 || mask(i,j,k) == 3.0) { has_SC = 1; } } diff --git a/Source/main.cpp b/Source/main.cpp index 30b1f9f..e84a4de 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -285,9 +285,9 @@ void main_main (c_FerroX& rFerroX) //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); @@ -405,9 +405,9 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_new_pre, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); @@ -475,9 +475,9 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_new, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); @@ -634,9 +634,9 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); From 9ad4f9c1cbd82915515c455664ad75312b0793e1 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Thu, 9 Nov 2023 14:19:14 -0800 Subject: [PATCH 8/8] treat metal as SC --- Source/Solver/ChargeDensity.cpp | 23 ++++++++++++++++------- Source/Solver/Initialization.cpp | 29 ++++++++++++++++++----------- 2 files changed, 34 insertions(+), 18 deletions(-) diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index 12de4d1..f7777a1 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -66,13 +66,22 @@ void ComputeRho(MultiFab& PoissonPhi, 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]); - } - // amrex::Print() << "Qe = " << Qe << "\n"; - charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); + //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]); + //} + // // amrex::Print() << "Qe = " << Qe << "\n"; + //charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); + + //Treat Metal as SC + 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)); + 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)); } else { diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index bdeb288..5e1b8de 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -74,10 +74,10 @@ void InitializePandRho(Array &P_old, Real denominator = metal_screening_length/epsilon_de*(coth - csch + FE_thickness/(2.*epsilonX_fe)); Real Qe = -numerator/denominator; - amrex::Print() << "average_P_r = " << average_P_r << "\n"; - amrex::Print() << "numerator = " << numerator << "\n"; - amrex::Print() << "denominator = " << denominator << "\n"; - amrex::Print() << "Qe = " << Qe << "\n"; + //amrex::Print() << "average_P_r = " << average_P_r << "\n"; + //amrex::Print() << "numerator = " << numerator << "\n"; + //amrex::Print() << "denominator = " << denominator << "\n"; + //amrex::Print() << "Qe = " << Qe << "\n"; // loop over boxes for (MFIter mfi(rho); mfi.isValid(); ++mfi) @@ -159,13 +159,20 @@ void InitializePandRho(Array &P_old, 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]); - } - //if(Qe != 0.) amrex::Print() << "initialization : Qe = " << Qe << "\n"; - charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); + //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]); + //} + ////if(Qe != 0.) amrex::Print() << "initialization : Qe = " << Qe << "\n"; + //charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); + + //Treat Metal as Semiconductor + hole_den_arr(i,j,k) = intrinsic_carrier_concentration; + e_den_arr(i,j,k) = intrinsic_carrier_concentration; + 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)); } else { if(use_Fermi_Dirac == 1){