From 7b51e247de665a4b80e045fb975fa38c3afe6d8a Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Fri, 29 Mar 2024 11:43:39 -0700 Subject: [PATCH 1/6] add nucleation parser and sro params --- Source/FerroX.cpp | 31 +++++-- Source/FerroX_namespace.H | 2 + Source/Solver/ChargeDensity.cpp | 17 +++- Source/Solver/ElectrostaticSolver.H | 2 +- Source/Solver/ElectrostaticSolver.cpp | 24 ++++++ Source/Solver/Initialization.H | 1 + Source/Solver/Initialization.cpp | 119 ++++++++++++++++++++++++-- Source/Solver/TotalEnergyDensity.cpp | 1 - Source/main.cpp | 12 ++- 9 files changed, 186 insertions(+), 23 deletions(-) diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 91544d5..feed3f3 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -212,6 +212,8 @@ 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; @@ -447,16 +449,23 @@ void InitializeFerroXNamespace(const amrex::GpuArray &P_old, MultiFab& MaterialMask); diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index bda6a91..1e837cd 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -521,3 +521,27 @@ void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray PoissonPhi.FillBoundary(geom.periodicity()); } + +void SetNucleation(Array &P_old, MultiFab& NucleationMask) +{ + for (MFIter mfi(P_old[0]); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + + 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& mask = NucleationMask.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + if (mask(i,j,k) == 1.) { + pOld_r(i,j,k) = Remnant_P[2]; + } + if (mask(i,j,k) == -1.) { + pOld_r(i,j,k) = -Remnant_P[2]; + } + }); + } + +} diff --git a/Source/Solver/Initialization.H b/Source/Solver/Initialization.H index 8a929d5..a5c7c51 100644 --- a/Source/Solver/Initialization.H +++ b/Source/Solver/Initialization.H @@ -26,5 +26,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); diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 447adf4..d9ecf4c 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -2,6 +2,17 @@ #include "Utils/eXstaticUtils/eXstaticUtil.H" #include "../../Utils/SelectWarpXUtils/WarpXUtil.H" +// Approximation to the Fermi-Dirac Integral of Order 1/2 +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real FD_half(const amrex::Real eta) +{ + amrex::Real nu = std::pow(eta, 4.0) + 50.0 + 33.6 * eta * (1.0 - 0.68 * exp(-0.17 * std::pow((eta + 1.0), 2.0))); + amrex::Real xi = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu, 3./8.)); + amrex::Real integral = std::pow(exp(-eta) + xi, -1.0); + return integral; +} + + // INITIALIZE rho in SC region void InitializePandRho(Array &P_old, MultiFab& Gamma, @@ -148,14 +159,65 @@ void InitializePandRho(Array &P_old, //SC region if (mask(i,j,k) >= 2.0) { - 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; - } + //Following: http://dx.doi.org/10.1063/1.4825209 + + amrex::Real Ef; + if (mask(i,j,k) == 2.0){ + Ef = 0.0; + } else { + Ef = -q*Phi_Bc_hi; + } + + amrex::Real Phi = 0.0; + amrex::Real Eg = bandgap; + amrex::Real Chi = affinity; + amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc/Nv)/q; + amrex::Real Ec_corr = -q*(Phi - phi_ref) - Chi*q; + amrex::Real Ev_corr = Ec_corr - q*Eg; + + //g_A is the acceptor ground state degeneracy factor and is equal to 4 + //because in most semiconductors each acceptor level can accept one hole of either spin + //and the impurity level is doubly degenerate as a result of the two degenerate valence bands + //(heavy hole and light hole bands) at the \Gamma point. + + //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 Ea = acceptor_ionization_energy; + amrex::Real Ed = donor_ionization_energy; + + if(use_Fermi_Dirac == 1){ + //Fermi-Dirac + + Real eta_n = -(Ec_corr - q*Ef)/(kb*T); + Real eta_p = -(q*Ef - Ev_corr)/(kb*T); + e_den_arr(i,j,k) = Nc*FD_half(eta_n); + hole_den_arr(i,j,k) = Nv*FD_half(eta_p); - 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)); + 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)/(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) / (kb*T) )); + } 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)/(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) / (kb*T) )); + + } + + 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; + + } }); } for (int i = 0; i < 3; i++){ @@ -382,3 +444,48 @@ void Initialize_Euler_angles(c_FerroX& rFerroX, const Geometry& geom, MultiFab& angle_theta.FillBoundary(geom.periodicity()); } +// initialization of nucleation site mask with parser +void Initialize_nucleation_Mask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& NucleationMask) +{ + auto& rGprop = rFerroX.get_GeometryProperties(); + Box const& domain = rGprop.geom.Domain(); + + const auto dx = rGprop.geom.CellSizeArray(); + const auto& real_box = rGprop.geom.ProbDomain(); + const auto iv = NucleationMask.ixType().toIntVect(); + + for (MFIter mfi(NucleationMask, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const auto& mask_arr = NucleationMask.array(mfi); + const auto& bx = mfi.tilebox(); + + std::string nucleation_mask_s; + std::unique_ptr nucleation_mask_parser; + std::string m_str_nucleation_geom_function; + + ParmParse pp_mask("nucleation_geom"); + + + if (pp_mask.query("nucleation_geom_function(x,y,z)", m_str_nucleation_geom_function) ) { + nucleation_mask_s = "parse_nucleation_geom_function"; + } + + if (nucleation_mask_s == "parse_nucleation_geom_function") { + Store_parserString(pp_mask, "nucleation_geom_function(x,y,z)", m_str_nucleation_geom_function); + nucleation_mask_parser = std::make_unique( + makeParser(m_str_nucleation_geom_function,{"x","y","z"})); + } + + const auto& macro_parser = nucleation_mask_parser->compile<3>(); + + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + eXstatic_MFab_Util::ConvertParserIntoMultiFab_3vars(i,j,k,dx,real_box,iv,macro_parser,mask_arr); + }); + + } + NucleationMask.FillBoundary(geom.periodicity()); +} + + diff --git a/Source/Solver/TotalEnergyDensity.cpp b/Source/Solver/TotalEnergyDensity.cpp index 53ef0bc..387c3a9 100644 --- a/Source/Solver/TotalEnergyDensity.cpp +++ b/Source/Solver/TotalEnergyDensity.cpp @@ -167,4 +167,3 @@ void CalculateTDGL_RHS(Array &GL_rhs, } } - diff --git a/Source/main.cpp b/Source/main.cpp index 86b708d..f738089 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -121,6 +121,7 @@ void main_main (c_FerroX& rFerroX) MultiFab charge_den(ba, dm, 1, 0); MultiFab MaterialMask(ba, dm, 1, 1); MultiFab tphaseMask(ba, dm, 1, 1); + MultiFab NucleationMask(ba, dm, 1, 1); MultiFab angle_alpha(ba, dm, 1, 0); MultiFab angle_beta(ba, dm, 1, 0); MultiFab angle_theta(ba, dm, 1, 0); @@ -141,13 +142,14 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.setVal(0.); PoissonRHS.setVal(0.); tphaseMask.setVal(0.); + NucleationMask.setVal(0.); angle_alpha.setVal(0.); angle_beta.setVal(0.); 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); @@ -277,6 +279,9 @@ void main_main (c_FerroX& rFerroX) //InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, geom, prob_lo, prob_hi);//old InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi);//mask based + + Initialize_nucleation_Mask(rFerroX, geom, NucleationMask); + SetNucleation(P_old, NucleationMask); //Obtain self consisten Phi and rho Real tol = 1.e-5; @@ -352,7 +357,7 @@ void main_main (c_FerroX& rFerroX) MultiFab::Copy(Plt, beta_cc, 0, 11, 1, 0); MultiFab::Copy(Plt, MaterialMask, 0, 12, 1, 0); - MultiFab::Copy(Plt, tphaseMask, 0, 13, 1, 0); + MultiFab::Copy(Plt, NucleationMask, 0, 13, 1, 0); MultiFab::Copy(Plt, angle_alpha, 0, 14, 1, 0); MultiFab::Copy(Plt, angle_beta, 0, 15, 1, 0); MultiFab::Copy(Plt, angle_theta, 0, 16, 1, 0); @@ -385,6 +390,7 @@ void main_main (c_FerroX& rFerroX) P_new_pre[i].FillBoundary(geom.periodicity()); } + SetNucleation(P_new_pre, NucleationMask); /** * \brief dst = a*x + b*y */ From 491809e3d2f6fd6f4a5269d9447fa6c311ecc20f Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Fri, 29 Mar 2024 12:51:44 -0700 Subject: [PATCH 2/6] option for both complete and partial ionization --- Source/Solver/ChargeDensity.H | 5 ++- Source/Solver/ChargeDensity.cpp | 47 +++++++++++++++------------ Source/Solver/ElectrostaticSolver.H | 2 ++ Source/Solver/ElectrostaticSolver.cpp | 4 ++- Source/Solver/Initialization.H | 2 ++ Source/Solver/Initialization.cpp | 25 +++++++++----- Source/main.cpp | 23 ++++++++----- 7 files changed, 68 insertions(+), 40 deletions(-) diff --git a/Source/Solver/ChargeDensity.H b/Source/Solver/ChargeDensity.H index af92798..f54a333 100644 --- a/Source/Solver/ChargeDensity.H +++ b/Source/Solver/ChargeDensity.H @@ -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); diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index c21f03b..8291fb4 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -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) { @@ -67,14 +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, g_D; - if (complete_ionization == 0){ - g_A = 4.0; - g_D = 2.0; - } else { - g_A = 0.0; //acceptor_den_arr(i,j,k) = acceptor_doping - g_D = 0.0; //donor_den_arr(i,j,k) = donor_doping - } + amrex::Real g_A = 4.0; + amrex::Real g_D = 2.0; amrex::Real Ea = acceptor_ionization_energy; amrex::Real Ed = donor_ionization_energy; @@ -87,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)); @@ -110,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()); } diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index ed1e174..7b20bee 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -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, diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 1e837cd..179270b 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -87,6 +87,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, @@ -102,7 +104,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); diff --git a/Source/Solver/Initialization.H b/Source/Solver/Initialization.H index a5c7c51..8e771ab 100644 --- a/Source/Solver/Initialization.H +++ b/Source/Solver/Initialization.H @@ -12,6 +12,8 @@ void InitializePandRho(Array &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& n_cell, diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index d9ecf4c..e22a6ec 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -19,6 +19,8 @@ void InitializePandRho(Array &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& n_cell, @@ -143,9 +145,6 @@ void InitializePandRho(Array &P_old, }); // 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& hole_den_arr = p_den.array(mfi); const Array4& e_den_arr = e_den.array(mfi); const Array4& charge_den_arr = rho.array(mfi); @@ -197,17 +196,27 @@ void InitializePandRho(Array &P_old, 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)/(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) / (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)/(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) / (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)/(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) / (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)/(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) / (kb*T) )); + } else { + acceptor_den_arr(i,j,k) = acceptor_doping; + donor_den_arr(i,j,k) = donor_doping; + } } diff --git a/Source/main.cpp b/Source/main.cpp index f738089..6959689 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -119,6 +119,8 @@ void main_main (c_FerroX& rFerroX) MultiFab hole_den(ba, dm, 1, 0); MultiFab e_den(ba, dm, 1, 0); MultiFab charge_den(ba, dm, 1, 0); + MultiFab acceptor_den(ba, dm, 1, 0); + MultiFab donor_den(ba, dm, 1, 0); MultiFab MaterialMask(ba, dm, 1, 1); MultiFab tphaseMask(ba, dm, 1, 1); MultiFab NucleationMask(ba, dm, 1, 1); @@ -139,6 +141,9 @@ void main_main (c_FerroX& rFerroX) e_den.setVal(0.); hole_den.setVal(0.); + charge_den.setVal(0.); + acceptor_den.setVal(0.); + donor_den.setVal(0.); PoissonPhi.setVal(0.); PoissonRHS.setVal(0.); tphaseMask.setVal(0.); @@ -278,7 +283,7 @@ void main_main (c_FerroX& rFerroX) // INITIALIZE P in FE and rho in SC regions //InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, geom, prob_lo, prob_hi);//old - InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi);//mask based + InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, acceptor_den, donor_den, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi);//mask based Initialize_nucleation_Mask(rFerroX, geom, NucleationMask); SetNucleation(P_old, NucleationMask); @@ -293,7 +298,7 @@ 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, acceptor_den, donor_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -310,7 +315,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, charge_den, e_den, hole_den, acceptor_den, donor_den, MaterialMask, geom); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -415,7 +420,7 @@ 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, acceptor_den, donor_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -433,7 +438,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, charge_den, e_den, hole_den, acceptor_den, donor_den, MaterialMask, geom); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -486,7 +491,7 @@ 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, acceptor_den, donor_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -504,7 +509,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, charge_den, e_den, hole_den, acceptor_den, donor_den, MaterialMask, geom); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -646,7 +651,7 @@ 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, acceptor_den, donor_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -664,7 +669,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, charge_den, e_den, hole_den, acceptor_den, donor_den, MaterialMask, geom); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates From c4417f59e0ed99b40d4a23d72b0ffaf912ac6313 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Tue, 30 Apr 2024 11:20:51 -0700 Subject: [PATCH 3/6] uniform initialization of P --- Source/Solver/Initialization.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index e22a6ec..0cdbfc0 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -117,6 +117,12 @@ void InitializePandRho(Array &P_old, pOld_q(i,j,k) = Remnant_P[1]*exp(-(x*x/(2.0*5.e-9*5.e-9) + y*y/(2.0*5.e-9*5.e-9) + (z-1.5*DE_hi[2])*(z - 1.5*DE_hi[2])/(2.0*2.0e-9*2.0e-9))); pOld_r(i,j,k) = Remnant_P[2]*exp(-(x*x/(2.0*5.e-9*5.e-9) + y*y/(2.0*5.e-9*5.e-9) + (z-1.5*DE_hi[2])*(z - 1.5*DE_hi[2])/(2.0*2.0e-9*2.0e-9))); + } else if (prob_type == 4) { // uniform P + + pOld_p(i,j,k) = Remnant_P[0]; + pOld_q(i,j,k) = Remnant_P[1]; + pOld_r(i,j,k) = Remnant_P[2]; + } else { Abort("Invalid prob_type"); From 667ccb8c22eee33ec8a27cd3123c48b7255d2ae1 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Tue, 30 Apr 2024 11:26:40 -0700 Subject: [PATCH 4/6] add warning message for the new initialization --- Source/Solver/Initialization.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 0cdbfc0..934ad1d 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -47,11 +47,18 @@ void InitializePandRho(Array &P_old, "P is initialized for convergence test." "\n" "==================================""\n" ; + } else if (prob_type == 4) { + + amrex::Print() << "==================================""\n" + "Uniform P is initialized." "\n" + "==================================""\n" ; + } else { amrex::Print() << "Undefine problem type!! Set prob_type in input script." "\n" "prob_type = 1 for 2D problems" "\n" "prob_type = 2 for 3D problems" "\n" - "prob_type = 3 for convergence tests." "\n"; + "prob_type = 3 for convergence tests." "\n" + "prob_type = 4 uniform P." "\n"; amrex::Abort(); } From 47caecc8ede274798a17ff8f30dde530726ceedf Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Fri, 10 May 2024 10:15:24 -0700 Subject: [PATCH 5/6] noise --- Source/Solver/ElectrostaticSolver.H | 2 +- Source/Solver/ElectrostaticSolver.cpp | 44 ++++++++++++++++++++++----- Source/main.cpp | 4 +-- 3 files changed, 40 insertions(+), 10 deletions(-) diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index 7b20bee..6deba04 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -63,4 +63,4 @@ void SetPoissonBC(c_FerroX& rFerroX, std::array &P_old, MultiFab& MaterialMask); +void SetNucleation(Array &P_old, MultiFab& MaterialMask, const amrex::GpuArray& n_cell); diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 179270b..265c094 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -1,3 +1,4 @@ +#include #include "ElectrostaticSolver.H" #include "DerivativeAlgorithm.H" #include "ChargeDensity.H" @@ -524,8 +525,27 @@ void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray } -void SetNucleation(Array &P_old, MultiFab& NucleationMask) +void SetNucleation(Array &P_old, MultiFab& NucleationMask, const amrex::GpuArray& 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 rngs(nrand, 0.0); + + // generate random numbers on the host + for (int i=0; i &P_old, MultiFab& NucleationM const Array4 &pOld_r = P_old[2].array(mfi); const Array4& mask = NucleationMask.array(mfi); - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + 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.) { - pOld_r(i,j,k) = Remnant_P[2]; - } - if (mask(i,j,k) == -1.) { - pOld_r(i,j,k) = -Remnant_P[2]; - } + if (prob_type == 1) { //2D + + pOld_p(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[0]*0.001; + pOld_q(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[1]*0.001; + pOld_r(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[2]*0.001;; + + } else if (prob_type == 2) { //3D + + pOld_p(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[0]*0.001; + pOld_q(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[1]*0.001; + pOld_r(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[2]*0.001; + } + } }); } diff --git a/Source/main.cpp b/Source/main.cpp index 6959689..52ad80e 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -286,7 +286,7 @@ void main_main (c_FerroX& rFerroX) InitializePandRho(P_old, Gamma, charge_den, e_den, hole_den, acceptor_den, donor_den, MaterialMask, tphaseMask, n_cell, geom, prob_lo, prob_hi);//mask based Initialize_nucleation_Mask(rFerroX, geom, NucleationMask); - SetNucleation(P_old, NucleationMask); + SetNucleation(P_old, NucleationMask, n_cell); //Obtain self consisten Phi and rho Real tol = 1.e-5; @@ -395,7 +395,7 @@ void main_main (c_FerroX& rFerroX) P_new_pre[i].FillBoundary(geom.periodicity()); } - SetNucleation(P_new_pre, NucleationMask); + SetNucleation(P_new_pre, NucleationMask, n_cell); /** * \brief dst = a*x + b*y */ From d50c9346eb37c0021418e5db2e709230816ccbf0 Mon Sep 17 00:00:00 2001 From: Prabhat Kumar Date: Wed, 26 Jun 2024 12:08:33 -0700 Subject: [PATCH 6/6] noise --- Source/FerroX.cpp | 4 ++++ Source/FerroX_namespace.H | 1 + Source/Solver/DerivativeAlgorithm.H | 8 ++++---- Source/Solver/ElectrostaticSolver.cpp | 12 ++++++------ 4 files changed, 15 insertions(+), 10 deletions(-) diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index feed3f3..257beb9 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -220,6 +220,7 @@ AMREX_GPU_MANAGED amrex::Real FerroX::lambda; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_hi; AMREX_GPU_MANAGED amrex::GpuArray 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; @@ -360,6 +361,9 @@ void InitializeFerroXNamespace(const amrex::GpuArray P_BC_flag_lo; extern AMREX_GPU_MANAGED amrex::GpuArray P_BC_flag_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Remnant_P; + extern AMREX_GPU_MANAGED amrex::Real noise_amplitude; //problem type : initialization of P for 2D/3D/convergence problems diff --git a/Source/Solver/DerivativeAlgorithm.H b/Source/Solver/DerivativeAlgorithm.H index 80ee376..d8e3a6c 100644 --- a/Source/Solver/DerivativeAlgorithm.H +++ b/Source/Solver/DerivativeAlgorithm.H @@ -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; @@ -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; diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 265c094..269e403 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -562,15 +562,15 @@ void SetNucleation(Array &P_old, MultiFab& NucleationM 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]*0.001; - pOld_q(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[1]*0.001; - pOld_r(i,j,k) += (-1.0 + 2.0*rng[i + k*n_cell[2]])*Remnant_P[2]*0.001;; + 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]*0.001; - pOld_q(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[1]*0.001; - pOld_r(i,j,k) += (-1.0 + 2.0*Random(engine))*Remnant_P[2]*0.001; + 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; } } });