diff --git a/Source/FerroX.H b/Source/FerroX.H index 8541df5..3b854ab 100644 --- a/Source/FerroX.H +++ b/Source/FerroX.H @@ -128,8 +128,12 @@ void WritePlotfile(c_FerroX& rFerroX, MultiFab& PoissonRHS, Array< MultiFab, AMREX_SPACEDIM>& P_old, Array< MultiFab, AMREX_SPACEDIM>& E, + Array< MultiFab, AMREX_SPACEDIM>& Jn, + Array< MultiFab, AMREX_SPACEDIM>& Jp, MultiFab& hole_den, MultiFab& e_den, + MultiFab& acceptor_den, + MultiFab& donor_den, MultiFab& charge_den, MultiFab& beta_cc, MultiFab& MaterialMask, diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index bf8ad0d..6063239 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -170,8 +170,12 @@ AMREX_GPU_MANAGED amrex::Real FerroX::dt; int FerroX::plot_Phi; int FerroX::plot_PoissonRHS; int FerroX::plot_E; +int FerroX::plot_Jn; +int FerroX::plot_Jp; int FerroX::plot_holes; int FerroX::plot_electrons; +int FerroX::plot_acceptors; +int FerroX::plot_donors; int FerroX::plot_charge; int FerroX::plot_epsilon; int FerroX::plot_mask; @@ -189,8 +193,10 @@ AMREX_GPU_MANAGED amrex::GpuArray FerroX::SC_lo; AMREX_GPU_MANAGED amrex::GpuArray FerroX::DE_hi; 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::p_type_hi; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::p_type_lo; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::n_type_hi; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::n_type_lo; // material parameters AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_0; @@ -228,6 +234,13 @@ 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; + +//Transport Model +AMREX_GPU_MANAGED amrex::Real FerroX::electron_mobility; +AMREX_GPU_MANAGED amrex::Real FerroX::hole_mobility; +AMREX_GPU_MANAGED amrex::Real FerroX::electron_diffusion_coefficient; +AMREX_GPU_MANAGED amrex::Real FerroX::hole_diffusion_coefficient; + // P and Phi Bc AMREX_GPU_MANAGED amrex::Real FerroX::lambda; AMREX_GPU_MANAGED amrex::GpuArray FerroX::P_BC_flag_lo; @@ -303,10 +316,18 @@ void InitializeFerroXNamespace(const amrex::GpuArray DE_hi; extern AMREX_GPU_MANAGED amrex::GpuArray FE_hi; extern AMREX_GPU_MANAGED 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 p_type_hi; + extern AMREX_GPU_MANAGED amrex::GpuArray p_type_lo; + extern AMREX_GPU_MANAGED amrex::GpuArray n_type_hi; + extern AMREX_GPU_MANAGED amrex::GpuArray n_type_lo; // material parameters extern AMREX_GPU_MANAGED amrex::Real epsilon_0; @@ -65,7 +71,12 @@ namespace FerroX { extern AMREX_GPU_MANAGED int use_Fermi_Dirac; extern AMREX_GPU_MANAGED int use_work_function; extern AMREX_GPU_MANAGED amrex::Real metal_work_function; - + + //Transport model + extern AMREX_GPU_MANAGED amrex::Real electron_mobility; + extern AMREX_GPU_MANAGED amrex::Real hole_mobility; + extern AMREX_GPU_MANAGED amrex::Real electron_diffusion_coefficient; + extern AMREX_GPU_MANAGED amrex::Real hole_diffusion_coefficient; // P and Phi Bc extern AMREX_GPU_MANAGED amrex::Real lambda; diff --git a/Source/Input/GeometryProperties/EmbeddedBoundaries.H b/Source/Input/GeometryProperties/EmbeddedBoundaries.H index 974c489..3d56598 100644 --- a/Source/Input/GeometryProperties/EmbeddedBoundaries.H +++ b/Source/Input/GeometryProperties/EmbeddedBoundaries.H @@ -77,7 +77,15 @@ public: amrex::MultiFab& get_soln_mf (int i) {return (*m_p_soln_mf[i]);} //amrex::MultiFab& get_beta_mf (int i) {return (*m_p_beta_mf[i]);} - void clear_soln_mf() { m_p_soln_mf.clear(); }; + //void clear_soln_mf() { m_p_soln_mf.clear(); }; + void clear_soln_mf() + { + for (int i = 0; i < m_p_soln_mf.size(); ++i) + { + m_p_soln_mf[i].reset(); + } + m_p_soln_mf.clear(); + }; //void clear_beta_mf() { m_p_beta_mf.clear(); }; void BuildGeometry(const amrex::Geometry* geom, const amrex::BoxArray* ba, const amrex::DistributionMapping* dm); diff --git a/Source/Plotfile.cpp b/Source/Plotfile.cpp index 0887f47..841d646 100644 --- a/Source/Plotfile.cpp +++ b/Source/Plotfile.cpp @@ -7,8 +7,12 @@ void WritePlotfile(c_FerroX& rFerroX, MultiFab& PoissonRHS, Array< MultiFab, AMREX_SPACEDIM>& P_old, Array< MultiFab, AMREX_SPACEDIM>& E, + Array< MultiFab, AMREX_SPACEDIM>& Jn, + Array< MultiFab, AMREX_SPACEDIM>& Jp, MultiFab& hole_den, MultiFab& e_den, + MultiFab& acceptor_den, + MultiFab& donor_den, MultiFab& charge_den, MultiFab& beta_cc, MultiFab& MaterialMask, @@ -55,6 +59,20 @@ void WritePlotfile(c_FerroX& rFerroX, var_names.push_back("Ez"); } + if (plot_Jn) { + nvar += 3; + var_names.push_back("Jnx"); + var_names.push_back("Jny"); + var_names.push_back("Jnz"); + } + + if (plot_Jp) { + nvar += 3; + var_names.push_back("Jpx"); + var_names.push_back("Jpy"); + var_names.push_back("Jpz"); + } + if (plot_holes) { ++nvar; var_names.push_back("holes"); @@ -65,6 +83,16 @@ void WritePlotfile(c_FerroX& rFerroX, var_names.push_back("electrons"); } + if (plot_acceptors) { + ++nvar; + var_names.push_back("acceptors"); + } + + if (plot_donors) { + ++nvar; + var_names.push_back("donors"); + } + if (plot_charge) { ++nvar; var_names.push_back("charge"); @@ -134,6 +162,18 @@ void WritePlotfile(c_FerroX& rFerroX, MultiFab::Copy(Plt, E[2], 0, counter++, 1, 0); } + if (plot_Jn) { + MultiFab::Copy(Plt, Jn[0], 0, counter++, 1, 0); + MultiFab::Copy(Plt, Jn[1], 0, counter++, 1, 0); + MultiFab::Copy(Plt, Jn[2], 0, counter++, 1, 0); + } + + if (plot_Jp) { + MultiFab::Copy(Plt, Jp[0], 0, counter++, 1, 0); + MultiFab::Copy(Plt, Jp[1], 0, counter++, 1, 0); + MultiFab::Copy(Plt, Jp[2], 0, counter++, 1, 0); + } + if (plot_holes) { MultiFab::Copy(Plt, hole_den, 0, counter++, 1, 0); } @@ -142,6 +182,14 @@ void WritePlotfile(c_FerroX& rFerroX, MultiFab::Copy(Plt, e_den, 0, counter++, 1, 0); } + if (plot_acceptors) { + MultiFab::Copy(Plt, acceptor_den, 0, counter++, 1, 0); + } + + if (plot_donors) { + MultiFab::Copy(Plt, donor_den, 0, counter++, 1, 0); + } + if (plot_charge) { MultiFab::Copy(Plt, charge_den, 0, counter++, 1, 0); } @@ -175,4 +223,7 @@ void WritePlotfile(c_FerroX& rFerroX, } WriteSingleLevelPlotfile(pltfile, Plt, var_names, geom, time, plt_step); -} \ No newline at end of file +#ifdef AMREX_USE_EB + EB_WriteSingleLevelPlotfile(pltfile, Plt, var_names, geom, time, plt_step); +#endif +} diff --git a/Source/Solver/ChargeDensity.H b/Source/Solver/ChargeDensity.H index af92798..411000e 100644 --- a/Source/Solver/ChargeDensity.H +++ b/Source/Solver/ChargeDensity.H @@ -12,3 +12,35 @@ void ComputeRho(MultiFab& PoissonPhi, MultiFab& p_den, const MultiFab& MaterialMask); +void CalculateDriftDiffusionCurrents( + amrex::Array& Jn, // OUT: Electron current density components (x,y,z) + amrex::Array& Jp, // OUT: Hole current density components (x,y,z) + const amrex::MultiFab& e_den, // IN: Electron density + const amrex::MultiFab& p_den, // IN: Hole density + const amrex::MultiFab& PoissonPhi, // IN: Electric potential + const amrex::Geometry& geom // IN: Simulation geometry +); +/* +void Compute_Effective_Potentials(const MultiFab& PoissonPhi, + MultiFab& e_potential, + MultiFab& p_potential, + const Geometry& geom); +*/ + +void Compute_Effective_Potentials(const MultiFab& PoissonPhi, + const MultiFab& e_den, + const MultiFab& p_den, + MultiFab& e_potential, + MultiFab& p_potential, + const Geometry& geom); + +void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi, + MultiFab& rho, + Array &Jn, + Array &Jp, + MultiFab& e_den, + MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, + MultiFab& MaterialMask, + const Geometry& geom); diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index f19f18b..597a176 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -1,4 +1,5 @@ #include "ChargeDensity.H" +#include "DerivativeAlgorithm.H" // Approximation to the Fermi-Dirac Integral of Order 1/2 AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -10,7 +11,6 @@ amrex::Real FD_half(const amrex::Real eta) return integral; } - // Compute rho in SC region for given phi void ComputeRho(MultiFab& PoissonPhi, MultiFab& rho, @@ -18,6 +18,7 @@ void ComputeRho(MultiFab& PoissonPhi, MultiFab& p_den, const MultiFab& MaterialMask) { + amrex::Print() << "Calculating steady-state carrier distribution." << "\n"; //Define acceptor and donor multifabs for doping and fill them with zero. MultiFab acceptor_den(rho.boxArray(), rho.DistributionMap(), 1, 0); @@ -67,7 +68,20 @@ void ComputeRho(MultiFab& PoissonPhi, amrex::Real Ea = acceptor_ionization_energy; amrex::Real Ed = donor_ionization_energy; - + + amrex::Real Na, Nd; + + if (mask(i,j,k) == 2.0) {//intrinsic + Na = 0.0; + Nd = 0.0; + } else if (mask(i,j,k) == 3.0) { // p-type + Na = acceptor_doping; + Nd = 0.0; + } else if (mask(i,j,k) == 4.0) { // n-type + Na = 0.0; + Nd = donor_doping; + } + if(use_Fermi_Dirac == 1){ //Fermi-Dirac @@ -76,8 +90,8 @@ 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) )); + acceptor_den_arr(i,j,k) = Na/(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) = Nd/(1.0 + g_D*exp( (q*Ef + q*Ed - q*phi_ref + q*Chi + q*phi(i,j,k)) / (kb*T) )); } else { @@ -85,8 +99,8 @@ void ComputeRho(MultiFab& PoissonPhi, 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) )); + acceptor_den_arr(i,j,k) = Na/(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) = Nd/(1.0 + g_D*exp( (q*Ef + q*Ed - q*phi_ref + q*Chi + q*phi(i,j,k)) / (kb*T) )); } @@ -101,3 +115,489 @@ void ComputeRho(MultiFab& PoissonPhi, } } + +// --- Bernoulli Function Implementation --- +// This is a crucial part of the Sharfetter-Gummel scheme. +// It's defined as Bern(x) = x / (exp(x) - 1). +// Special care is needed for x close to 0 to avoid division by zero (use Taylor expansion). +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real Bern(amrex::Real x) +{ + // Use a small epsilon for robustness around x=0 + if (amrex::Math::abs(x) < 1.0e-6) { + // Taylor expansion for small x: 1 - x/2 + x^2/12 - x^4/720 + ... + return 1.0 - x/2.0 + x*x/12.0; + } else { + return x / (exp(x) - 1.0); + } +} + +// --- CalculateDriftDiffusionCurrents --- +void CalculateDriftDiffusionCurrents( + amrex::Array& Jn, // OUT: Electron current density components (x,y,z) + amrex::Array& Jp, // OUT: Hole current density components (x,y,z) + const amrex::MultiFab& e_den, // IN: Electron density + const amrex::MultiFab& p_den, // IN: Hole density + const amrex::MultiFab& MaterialMask, // IN: Hole density + const amrex::MultiFab& PoissonPhi, // IN: Electric potential + const amrex::Geometry& geom // IN: Simulation geometry +) +{ + + const amrex::Real kBT_over_q = kb * T / q; + const amrex::GpuArray dx = geom.CellSizeArray(); + + MultiFab e_potential(PoissonPhi.boxArray(), PoissonPhi.DistributionMap(), 1, 1); + MultiFab p_potential(PoissonPhi.boxArray(), PoissonPhi.DistributionMap(), 1, 1); + e_potential.setVal(0.); + p_potential.setVal(0.); + + Compute_Effective_Potentials(PoissonPhi, e_den, p_den, e_potential, p_potential, geom); + + for (amrex::MFIter mfi(e_den, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& bx = mfi.validbox(); + + // Get Array4 views for electron density, hole density, potential + amrex::Array4 const& e_den_arr = e_den.const_array(mfi); + amrex::Array4 const& p_den_arr = p_den.const_array(mfi); + amrex::Array4 const& phi_arr = PoissonPhi.const_array(mfi); + // amrex::Array4const& phi_n_arr = e_potential.array(mfi); + // amrex::Array4const& phi_p_arr = p_potential.array(mfi); + amrex::Array4 const& phi_n_arr = e_potential.const_array(mfi); + amrex::Array4 const& phi_p_arr = p_potential.const_array(mfi); + amrex::Array4 const& mask = MaterialMask.const_array(mfi); + + // Get Array4 views for current components (output) + amrex::Array4 const& Jnx_arr = Jn[0].array(mfi); + amrex::Array4 const& Jny_arr = Jn[1].array(mfi); + amrex::Array4 const& Jnz_arr = Jn[2].array(mfi); + amrex::Array4 const& Jpx_arr = Jp[0].array(mfi); + amrex::Array4 const& Jpy_arr = Jp[1].array(mfi); + amrex::Array4 const& Jpz_arr = Jp[2].array(mfi); + + amrex::Real mu_n = electron_mobility; + amrex::Real mu_p = hole_mobility; + + amrex::Real D_n = electron_diffusion_coefficient; + amrex::Real D_p = hole_diffusion_coefficient; + + + amrex::ParallelFor(bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + // Initialize current components for this cell (will be overwritten) + 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; + + //phi_n_arr(i,j,k) = 1.0*phi_arr(i,j,k); + //phi_p_arr(i,j,k) = 1.0*phi_arr(i,j,k); + + if (mask(i,j,k) >= 2.0) { + // --- Calculate J_x (current across faces normal to x-axis) --- + // For a cell (i,j,k), we're interested in the current *through* its faces. + // Jnx_arr(i,j,k) will store the current through the FACE at (i+1/2, j, k). + // Jpx_arr(i,j,k) will store the current through the FACE at (i+1/2, j, k). + // Similarly for y and z. + + // Only compute if we are not at the very right boundary of the box for the current's central difference. + // If the box is not periodic, the currents at the domain boundaries will be handled by specific boundary conditions + // or implicitly be zero if not computed beyond the physical domain. + if (i <= bx.bigEnd(0)) { // This calculates J_x at (i+1/2, j, k) + amrex::Real dPhi_n = phi_n_arr(i+1, j, k) - phi_n_arr(i, j, k); + amrex::Real arg_n = dPhi_n / kBT_over_q; + + amrex::Real dPhi_p = phi_p_arr(i+1, j, k) - phi_p_arr(i, j, k); + amrex::Real arg_p = dPhi_p / kBT_over_q; + + Jnx_arr(i, j, k) = q * D_n / dx[0] * (e_den_arr(i+1,j,k) * Bern(arg_n) - e_den_arr(i,j,k) * Bern(-arg_n)); + Jpx_arr(i, j, k) = q * D_p / dx[0] * (p_den_arr(i,j,k) * Bern(-arg_p) - p_den_arr(i+1,j,k) * Bern(arg_p)); + } + + // --- Calculate J_y (current across faces normal to y-axis) --- + if (j <= bx.bigEnd(1)) { // This calculates J_y at (i, j+1/2, k) + amrex::Real dPhi_n = phi_n_arr(i, j+1, k) - phi_n_arr(i, j, k); + amrex::Real arg_n = dPhi_n / kBT_over_q; + + amrex::Real dPhi_p = phi_p_arr(i, j+1, k) - phi_p_arr(i, j, k); + amrex::Real arg_p = dPhi_p / kBT_over_q; + + Jny_arr(i, j, k) = q * D_n / dx[1] * (e_den_arr(i,j+1,k) * Bern(arg_n) - e_den_arr(i,j,k) * Bern(-arg_n)); + Jpy_arr(i, j, k) = q * D_p / dx[1] * (p_den_arr(i,j,k) * Bern(-arg_p) - p_den_arr(i,j+1,k) * Bern(arg_p)); + } + + // --- Calculate J_z (current across faces normal to z-axis) --- + if (k <= bx.bigEnd(2)) { // This calculates J_z at (i, j, k+1/2) + amrex::Real dPhi_n = phi_n_arr(i, j, k+1) - phi_n_arr(i, j, k); + amrex::Real arg_n = dPhi_n / kBT_over_q; + + amrex::Real dPhi_p = phi_p_arr(i, j, k+1) - phi_p_arr(i, j, k); + amrex::Real arg_p = dPhi_p / kBT_over_q; + + Jnz_arr(i, j, k) = q * D_n / dx[2] * (e_den_arr(i,j,k+1) * Bern(arg_n) - e_den_arr(i,j,k) * Bern(-arg_n)); + Jpz_arr(i, j, k) = q * D_p / dx[2] * (p_den_arr(i,j,k) * Bern(-arg_p) - p_den_arr(i,j,k+1) * Bern(arg_p)); + } + } + }); + } + + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + Jn[d].FillBoundary(geom.periodicity()); + Jp[d].FillBoundary(geom.periodicity()); + } +} + +// Compute rho in SC region for given phi +void ComputeRho_DriftDiffusion(MultiFab& PoissonPhi, + MultiFab& rho, + Array &Jn, + Array &Jp, + MultiFab& e_den, + MultiFab& p_den, + MultiFab& e_den_old, + MultiFab& p_den_old, + MultiFab& MaterialMask, + const Geometry& geom) +{ + + amrex::Print() << "Calculating carrier transport using Drift-Diffusion model." << "\n"; + + //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.); + + // First, calculate the current components and store them in Jn and Jp + CalculateDriftDiffusionCurrents(Jn, Jp, e_den, p_den, MaterialMask, PoissonPhi, geom); + + // Get cell spacing from geometry + const amrex::GpuArray dx = geom.CellSizeArray(); + + // Loop over grids (boxes) in the MultiFab for updating densities + for (amrex::MFIter mfi(e_den, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& bx = mfi.tilebox(); + + // Get Array4 views for densities (for update) + amrex::Array4 const& e_den_arr = e_den.array(mfi); + amrex::Array4 const& p_den_arr = p_den.array(mfi); + amrex::Array4 const& charge_den_arr = rho.array(mfi); + amrex::Array4 const& phi = PoissonPhi.array(mfi); + const Array4& acceptor_den_arr = acceptor_den.array(mfi); + const Array4& donor_den_arr = donor_den.array(mfi); + const Array4& mask = MaterialMask.array(mfi); + + // Get Array4 views for current components (from `CalculateDriftDiffusionCurrents`) + amrex::Array4 const& Jnx_arr = Jn[0].const_array(mfi); + amrex::Array4 const& Jny_arr = Jn[1].const_array(mfi); + amrex::Array4 const& Jnz_arr = Jn[2].const_array(mfi); + amrex::Array4 const& Jpx_arr = Jp[0].const_array(mfi); + amrex::Array4 const& Jpy_arr = Jp[1].const_array(mfi); + amrex::Array4 const& Jpz_arr = Jp[2].const_array(mfi); + + amrex::Real ni_sq_val = Nc * Nv * exp(-q*bandgap / (kb * T)); // Assuming bandgap is in eV, kb*T in eV + + amrex::Real ni_val = std::sqrt(ni_sq_val); + + amrex::Real tau_n_val = 1.0e-4; //taun_const; // Example: Pass as captured variable or global + amrex::Real tau_p_val = 1.0e-4; //taup_const; // Example: Pass as captured variable or global + + amrex::ParallelFor(bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + + amrex::Real Ef = 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; + + if (mask(i,j,k) >= 2.0) { + + // --- Calculate Divergence --- + // div(J) = (J_x_R - J_x_L)/dx + (J_y_R - J_y_L)/dy + (J_z_R - J_z_L)/dz + // Note: J_x_R for cell (i,j,k) is Jnx_arr(i,j,k) (current at i+1/2 face) + // J_x_L for cell (i,j,k) is Jnx_arr(i-1,j,k) (current at i-1/2 face) + + amrex::Real div_Jn = (Jnx_arr(i, j, k) - Jnx_arr(i-1, j, k)) / dx[0] + + (Jny_arr(i, j, k) - Jny_arr(i, j-1, k)) / dx[1] + + (Jnz_arr(i, j, k) - Jnz_arr(i, j, k-1)) / dx[2]; + + amrex::Real div_Jp = (Jpx_arr(i, j, k) - Jpx_arr(i-1, j, k)) / dx[0] + + (Jpy_arr(i, j, k) - Jpy_arr(i, j-1, k)) / dx[1] + + (Jpz_arr(i, j, k) - Jpz_arr(i, j, k-1)) / dx[2]; + + // --- Calculate SRH Net Recombination Rate (R_SRH) --- + amrex::Real SRH_numerator = (e_den_arr(i, j, k) * p_den_arr(i, j, k)) - ni_sq_val; + amrex::Real SRH_denominator = tau_p_val * (e_den_arr(i, j, k) + ni_val) + tau_n_val * (p_den_arr(i, j, k) + ni_val); + + // Handle potential division by zero if both carrier densities and ni are very small + // For practical device simulations, carrier densities are rarely zero in active regions. + // If they could be, you might need a small epsilon or check. + amrex::Real R_SRH = 0.0; + if (SRH_denominator > 1.0e-30) { // Add a small epsilon to avoid division by zero + R_SRH = SRH_numerator / SRH_denominator; + } + +// if(i == 0 && j == 0 && k == 32) amrex::Print() << "R_SRH = " << R_SRH << ", ni_val = " << ni_val << "\n"; + + // --- Update Densities (including recombination term) --- + // For electrons: q dn/dt = div(Jn) - qR + // dn/dt = (1/q) * div(Jn) - R + e_den_arr(i, j, k) += dt * ((1.0/q * div_Jn) - R_SRH); + + // For holes: q dp/dt = -div(Jp) - qR + // dp/dt = (-1/q) * div(Jp) - R + p_den_arr(i, j, k) += dt * ((-1.0/q * div_Jp) - R_SRH); + + // --- Update Densities (ignoring recombination) --- + //e_den_arr(i, j, k) += dt * (1.0/q * div_Jn); + //p_den_arr(i, j, k) += dt * (-1.0/q * div_Jp); + + e_den_arr(i, j, -1) = 0.5*donor_doping + std::sqrt(std::pow(0.5*donor_doping,2.0) + intrinsic_carrier_concentration * intrinsic_carrier_concentration); + p_den_arr(i, j, -1) = intrinsic_carrier_concentration * intrinsic_carrier_concentration / e_den_arr(i, j, 0); + + e_den_arr(i, j, 0) = 0.5*donor_doping + std::sqrt(std::pow(0.5*donor_doping,2.0) + intrinsic_carrier_concentration * intrinsic_carrier_concentration); + p_den_arr(i, j, 0) = intrinsic_carrier_concentration * intrinsic_carrier_concentration / e_den_arr(i, j, 0); + + p_den_arr(i, j, 63) = 0.5*acceptor_doping + std::sqrt(std::pow(0.5*acceptor_doping,2.0) + intrinsic_carrier_concentration * intrinsic_carrier_concentration); + e_den_arr(i, j, 63) = intrinsic_carrier_concentration * intrinsic_carrier_concentration / p_den_arr(i, j, 63); + + p_den_arr(i, j, 64) = 0.5*acceptor_doping + std::sqrt(std::pow(0.5*acceptor_doping,2.0) + intrinsic_carrier_concentration * intrinsic_carrier_concentration); + e_den_arr(i, j, 64) = intrinsic_carrier_concentration * intrinsic_carrier_concentration / p_den_arr(i, j, 64); + + + //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. + + //setting it to zero assuming complete ionization + amrex::Real g_A = 0.0; //4.0; + amrex::Real g_D = 0.0; //2.0; + + amrex::Real Ea = acceptor_ionization_energy; + amrex::Real Ed = donor_ionization_energy; + + amrex::Real Na, Nd; + + if (mask(i,j,k) == 2.0) {//intrinsic + Na = 0.0; + Nd = 0.0; + } else if (mask(i,j,k) == 3.0) { // p-type + Na = acceptor_doping; + Nd = 0.0; + } else if (mask(i,j,k) == 4.0) { // n-type + Na = 0.0; + Nd = donor_doping; + } + + acceptor_den_arr(i,j,k) = Na/(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) = Nd/(1.0 + g_D*exp( (q*Ef + q*Ed - q*phi_ref + q*Chi + q*phi(i,j,k)) / (kb*T) )); + + charge_den_arr(i,j,k) = q*(p_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); + } + }); + } + + // After updating, fill ghost cells for density multifabs if needed + e_den.FillBoundary(geom.periodicity()); + p_den.FillBoundary(geom.periodicity()); + rho.FillBoundary(geom.periodicity()); + } +/* +void Compute_Effective_Potentials(const MultiFab& PoissonPhi, + MultiFab& e_potential, + MultiFab& p_potential, + const Geometry& geom) +{ + + // loop over boxes + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + + const Array4& phi = PoissonPhi.array(mfi); + const Array4& e_phi = e_potential.array(mfi); + const Array4& p_phi = p_potential.array(mfi); + + amrex::Real Ef = 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 Delta_Eg = 0.; // Delta_Eg is an energy. Assume it's in eV. Convert to Joules if needed. + + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + + amrex::Real Ec = -q*(phi(i,j,k) - phi_ref) - Chi*q; + amrex::Real Ev = Ec - q*Eg; + + // eta_n and eta_p are dimensionless (Energy/Energy) + amrex::Real eta_n = (q*Ef - Ec)/(kb*T); + amrex::Real eta_p = (Ev - q*Ef)/(kb*T); + + amrex::Real gamma_n = 1.0;// FD_half(eta_n)/exp(eta_n);// 1.; // Default for Maxwell-Boltzmann + amrex::Real gamma_p = 1.0;// FD_half(eta_p)/exp(eta_p); //1.; // Default for Maxwell-Boltzmann + + amrex::Real E_i_J = q*phi_ref + - q*Chi + - q*phi(i,j,k) // phi is in Volts, q*phi is Joules + - 0.5*q*Eg + - 0.5*kb*T*log( (Nc*gamma_n) / (Nv*gamma_p) ); + + // Calculate the effective potentials in Joules + amrex::Real E_n_eff_J = E_i_J - 0.5*Delta_Eg - 0.5*kb*T*log(gamma_n*gamma_p); + amrex::Real E_p_eff_J = E_i_J + 0.5*Delta_Eg + 0.5*kb*T*log(gamma_n*gamma_p); + + e_phi(i,j,k) = 1./q*E_n_eff_J; // J + p_phi(i,j,k) = 1./q*E_p_eff_J; // J + }); + } + + e_potential.FillBoundary(geom.periodicity()); + p_potential.FillBoundary(geom.periodicity()); +} +*/ +// Approximation to the inverse of the Fermi-Dirac Integral of Order 1/2 +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real Inverse_FD_half(amrex::Real u) +{ + + amrex::Real sqrt_pi = std::sqrt(3.14); + amrex::Real nu = std::pow( (3.0 * sqrt_pi * u / 4.0), 2.0 / 3.0 ); + + amrex::Real log_term = -std::log(u) / (u*u - 1.0); + amrex::Real denom = 1.0 + std::pow(0.24 + 1.08 * nu, -2.0); + amrex::Real eta = log_term + nu / denom; + + return eta; +} + + +void Compute_Effective_Potentials(const MultiFab& PoissonPhi, + const MultiFab& e_den, + const MultiFab& p_den, + MultiFab& e_potential, + MultiFab& p_potential, + const Geometry& geom) + +{ + // Need to get Nc and Nv (effective density of states) + // Make sure these are consistent with your material parameters + // Example values for Silicon at 300K, typically in m^-3 + amrex::Real Nc_val = 2.8e25; // Example: 2.8e19 cm^-3 = 2.8e25 m^-3 + amrex::Real Nv_val = 1.04e25; // Example: 1.04e19 cm^-3 = 1.04e25 m^-3 + + // Constants from your previous code + amrex::Real Eg = bandgap; // Bandgap in eV + amrex::Real Chi = affinity; // Electron affinity in eV + amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc_val/Nv_val)/q; // Assuming consistent with your band definition + + // loop over boxes + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + + const Array4& phi = PoissonPhi.array(mfi); + const Array4& e_den_arr = e_den.const_array(mfi); + const Array4& p_den_arr = p_den.const_array(mfi); + + const Array4& e_phi = e_potential.array(mfi); + const Array4& p_phi = p_potential.array(mfi); + + amrex::ParallelFor( bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) // Add AMREX_GPU_HOST_DEVICE + { + // Calculate band edges (Ec, Ev) in Joules + amrex::Real Ec_J = -q*(phi(i,j,k) - phi_ref) - Chi*q; + amrex::Real Ev_J = Ec_J - q*Eg; + + // --- Electron Quasi-Fermi Level (E_Fn) --- + amrex::Real current_e_den = e_den_arr(i,j,k); + if (current_e_den <= 0.0) { + current_e_den = intrinsic_carrier_concentration; // Or a very small positive number + } + + amrex::Real E_Fn_J; + if (use_Fermi_Dirac == 1) { + // Fermi-Dirac statistics for electrons + amrex::Real ratio_n = current_e_den / Nc_val; + if (ratio_n < 1.0e-10) ratio_n = 1.0e-10; // Clamp for numerical stability + E_Fn_J = Ec_J + Inverse_FD_half(ratio_n) * kb*T; + // Make sure inverse_FD_half is accessible in GPU code if you use this branch. + // Placeholder for now: + //E_Fn_J = Ec_J + kb*T * log(current_e_den / Nc_val); // Fallback to MB if inverse_FD_half is not implemented + } else { + // Maxwell-Boltzmann statistics for electrons (non-degenerate) + E_Fn_J = Ec_J + kb*T * log(current_e_den / Nc_val); + } + + // --- Hole Quasi-Fermi Level (E_Fp) --- + amrex::Real current_p_den = p_den_arr(i,j,k); + if (current_p_den <= 0.0) { + current_p_den = intrinsic_carrier_concentration; // Or a very small positive number + } + + amrex::Real E_Fp_J; + if (use_Fermi_Dirac == 1) { + // Fermi-Dirac statistics for holes + amrex::Real ratio_p = current_p_den / Nv_val; + if (ratio_p < 1.0e-10) ratio_p = 1.0e-10; // Clamp for numerical stability + E_Fp_J = Ev_J - Inverse_FD_half(ratio_p) * kb*T; // Note the sign difference for holes + // Placeholder for now: + //E_Fp_J = Ev_J - kb*T * log(current_p_den / Nv_val); // Fallback to MB if inverse_FD_half is not implemented + } else { + // Maxwell-Boltzmann statistics for holes (non-degenerate) + E_Fp_J = Ev_J - kb*T * log(current_p_den / Nv_val); + } + + // Assign quasi-Fermi levels (in Volts) to the output MultiFabs + e_phi(i,j,k) = E_Fn_J / q; + p_phi(i,j,k) = E_Fp_J / q; + }); + + const int lo_z = bx.smallEnd(2); + const int hi_z = bx.bigEnd(2); + + // Left Boundary (Z = 0, assuming smallEnd(2) == 0 for the domain boundary) + if (lo_z == 0) { // Check if this box is at the global domain's left Z boundary + Box z_minus_1_ghost_slice = bx; // Start with the full box + z_minus_1_ghost_slice.setSmall(2, -1); // Set Z-low to -1 + z_minus_1_ghost_slice.setBig(2, -1); // Set Z-high to -1 (for a single slice) + + amrex::ParallelFor(z_minus_1_ghost_slice, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + // The '0' here refers to the actual physical cell index, not a relative index. + e_phi(i,j,k) = e_phi(i,j,0); // e_phi(i,j,-1) = e_phi(i,j,0) + p_phi(i,j,k) = p_phi(i,j,0); // p_phi(i,j,-1) = p_phi(i,j,0) + }); + } + + // Right Boundary (Z = 1000 nm, assuming bigEnd(2) == 63 for the domain boundary with a 0-indexed system up to 63) + if (hi_z == 63) { // Check if this box is at the global domain's right Z boundary + Box z_plus_1_ghost_slice = bx; // Start with the full box + z_plus_1_ghost_slice.setSmall(2, 64); // Set Z-low to 64 + z_plus_1_ghost_slice.setBig(2, 64); // Set Z-high to 64 (for a single slice) + + amrex::ParallelFor(z_plus_1_ghost_slice, + [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept + { + // The '63' here refers to the actual physical cell index. + e_phi(i,j,k) = e_phi(i,j,63); // e_phi(i,j,64) = e_phi(i,j,63) + p_phi(i,j,k) = p_phi(i,j,63); // p_phi(i,j,64) = p_phi(i,j,63) + }); + } + // --- End of corrected section for QFL BCs --- + } + + e_potential.FillBoundary(geom.periodicity()); + p_potential.FillBoundary(geom.periodicity()); +} diff --git a/Source/Solver/DerivativeAlgorithm.H b/Source/Solver/DerivativeAlgorithm.H index 80ee376..87d3f6b 100644 --- a/Source/Solver/DerivativeAlgorithm.H +++ b/Source/Solver/DerivativeAlgorithm.H @@ -39,6 +39,22 @@ 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 const& F, + amrex::Array4 const& mask, + int const i, int const j, int const k, amrex::GpuArray dx) { + if (mask(i, j, k+1) < 2.0 && mask(i,j,k) >= 2.0){ // SC higher boundary + return (F(i,j,k) - F(i,j,k-1))/(dx[2]); + } else if (mask(i, j, k-1) < 2.0 && mask(i,j,k) >= 2.0){ // SC lower boundary + return (F(i,j,k+1) - F(i,j,k))/(dx[2]); + } else { + 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 diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index 128e6cb..66eff08 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -49,9 +49,13 @@ void dF_dPhi(MultiFab& alpha_cc, MultiFab& PoissonRHS, MultiFab& PoissonPhi, Array& P_old, + Array& Jn, + Array& 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, @@ -62,7 +66,8 @@ void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& PoissonPhi, MultiFab& alpha_cc); -void SetPhiBC_z(MultiFab& PossonPhi, const amrex::GpuArray& n_cell, const Geometry& geom); +//void SetPhiBC_z(MultiFab& PossonPhi, const amrex::GpuArray& n_cell, const Geometry& geom); +void SetPhiBC_z(MultiFab& PoissonPhi, MultiFab& MaterialMask, const amrex::GpuArray& n_cell, const Geometry& geom); void SetPoissonBC(c_FerroX& rFerroX, std::array,2>& LinOpBCType_2d, bool& all_homogeneous_boundaries, bool& some_functionbased_inhomogeneous_boundaries, bool& some_constant_inhomogeneous_boundaries); @@ -75,6 +80,7 @@ void SetupMLMG(std::unique_ptr& pMLMG, std::array,2>& LinOpBCType_2d, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + MultiFab& MaterialMask, c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info); #ifdef AMREX_USE_EB @@ -83,6 +89,7 @@ void SetupMLMG_EB(std::unique_ptr& pMLMG, std::array,2>& LinOpBCType_2d, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + MultiFab& MaterialMask, MultiFab& beta_cc, c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info); #endif @@ -94,14 +101,39 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, MultiFab& PoissonPhi, MultiFab& PoissonPhi_Prev, MultiFab& PhiErr, - Array& P_old, + Array& P_old, MultiFab& rho, + Array& Jn, + Array& Jp, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, + 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& prob_lo, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi); + +void ComputePhi_Rho_Equilibrium(std::unique_ptr& pMLMG, + std::unique_ptr& p_mlabec, + MultiFab& alpha_cc, + MultiFab& PoissonRHS, + MultiFab& PoissonPhi, + MultiFab& PoissonPhi_Prev, + MultiFab& PhiErr, + Array& P_old, + MultiFab& rho, + Array& Jn, + Array& Jp, + 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& prob_lo, const amrex::GpuArray& prob_hi); #ifdef AMREX_USE_EB @@ -112,13 +144,17 @@ void ComputePhi_Rho_EB(std::unique_ptr& pMLMG, MultiFab& PoissonPhi, MultiFab& PoissonPhi_Prev, MultiFab& PhiErr, - Array& P_old, + Array& P_old, MultiFab& rho, + Array& Jn, + Array& Jp, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, + 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& prob_lo, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); #endif diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index feb80c7..d848e60 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -86,9 +86,13 @@ void dF_dPhi(MultiFab& alpha_cc, MultiFab& PoissonRHS, MultiFab& PoissonPhi, Array& P_old, + Array& Jn, + Array& 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, @@ -105,6 +109,7 @@ void dF_dPhi(MultiFab& alpha_cc, // Calculate rho from Phi in SC region ComputeRho(PoissonPhi_plus_delta, rho, e_den, p_den, MaterialMask); + //ComputeRho_DriftDiffusion(PoissonPhi_plus_delta, rho, Jn, Jp, e_den, p_den, e_den_old, p_den_old, MaterialMask, geom); //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); @@ -499,25 +504,152 @@ void Fill_FunctionBased_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& Po } } -void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray& n_cell, const Geometry& geom) +// Approximation to the inverse of the Fermi-Dirac Integral of Order 1/2 +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real Inverse_FD_half(amrex::Real u) +{ + + amrex::Real sqrt_pi = std::sqrt(3.14); + amrex::Real nu = std::pow( (3.0 * sqrt_pi * u / 4.0), 2.0 / 3.0 ); + + amrex::Real log_term = -std::log(u) / (u*u - 1.0); + amrex::Real denom = 1.0 + std::pow(0.24 + 1.08 * nu, -2.0); + amrex::Real eta = log_term + nu / denom; + + return eta; +} + +// --- Main Calculation Function for Poisson BC values --- +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::GpuArray CalculatePoissonBoundaryPotentials( + amrex::Real Nc, amrex::Real Nv, amrex::Real bandgap_eV, amrex::Real affinity_eV, + amrex::Real q_coulombs, amrex::Real kb_joules_per_k, amrex::Real T_kelvin, + amrex::Real donor_doping_val, amrex::Real acceptor_doping_val, + amrex::Real Phi_Bc_lo = 0.0, amrex::Real Phi_Bc_hi = 0.0) +{ + /** + * @brief Calculates the Dirichlet boundary conditions for the Poisson equation at + * the top and bottom boundaries of a p-n diode at zero bias. + * + * @param Nc Effective density of states in conduction band. + * @param Nv Effective density of states in valence band. + * @param bandgap_eV Bandgap energy in electron-Volts (eV). + * @param affinity_eV Electron affinity in electron-Volts (eV). + * @param q_coulombs Elementary charge in Coulombs. + * @param kb_joules_per_k Boltzmann constant in J/K. + * @param T_kelvin Temperature in Kelvin. + * @param donor_doping_val Donor concentration for n-type material. + * @param acceptor_doping_val Acceptor concentration for p-type material. + * @param Phi_Bc_lo Applied voltage at the bottom boundary (z=0). + * @param Phi_Bc_hi Applied voltage at the top boundary (z=1). + * @return A GpuArray containing the calculated BCs: [bc_at_bottom_ntype, bc_at_top_ptype]. + * Returns NaNs if Inverse_FD_half returns NaN. + */ + + amrex::Real kbT_over_q_V = (kb_joules_per_k * T_kelvin) / q_coulombs; + + //Calculate phi_ref_V (intrinsic Fermi level relative to vacuum, as a potential) + amrex::Real log_Nc_over_Nv = log(Nc / Nv); + amrex::Real phi_ref_V = affinity_eV + (0.5 * bandgap_eV) + (0.5 * kbT_over_q_V * log_Nc_over_Nv); + + // --- Calculate BC at n-type contact --- + amrex::Real u_ntype = donor_doping_val / Nc; + amrex::Real eta_ntype = Inverse_FD_half(u_ntype); + + if (amrex::Gpu::isnan(eta_ntype)) { + return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; + } + + // Formula for n-type: phi_ohm^n = phi_ref_V - Chi_eV + (kbT/q)*eta + Va + // affinity_eV is used directly as a potential in Volts. + amrex::Real bc_at_bottom_ntype = phi_ref_V - affinity_eV + (kbT_over_q_V * eta_ntype) + Phi_Bc_lo; + + // --- Calculate BC at p-type contact --- + amrex::Real u_ptype = acceptor_doping_val / Nv; + amrex::Real eta_ptype = Inverse_FD_half(u_ptype); + + if (amrex::Gpu::isnan(eta_ptype)) { + return {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}; + } + + // Formula for p-type: phi_ohm^p = phi_ref_V - Chi_eV - Eg_eV - (kbT/q)*eta + Va + // affinity_eV and bandgap_eV are used directly as potentials in Volts. + amrex::Real bc_at_top_ptype = phi_ref_V - affinity_eV - bandgap_eV - (kbT_over_q_V * eta_ptype) + Phi_Bc_hi; + + // Return the calculated boundary conditions + return {bc_at_bottom_ntype, bc_at_top_ptype}; +} + +void SetPhiBC_z(MultiFab& PoissonPhi, MultiFab& MaterialMask, const amrex::GpuArray& n_cell, const Geometry& geom) { for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(1); const Array4& Phi = PoissonPhi.array(mfi); + const Array4& mask = MaterialMask.array(mfi); + + amrex::Real Eg_eV = bandgap; // Bandgap energy in eV + amrex::Real Chi_eV = affinity; // Electron affinity energy in eV + + //Calculate phi_ref_V (intrinsic Fermi level relative to vacuum, as a potential) + amrex::Real log_Nc_over_Nv = log(Nc / Nv); + amrex::Real phi_ref_V = affinity + (0.5 * bandgap) + (0.5 * kb * T / q * log_Nc_over_Nv); + + amrex::Real phi_m_V = use_work_function ? metal_work_function : phi_ref_V; + + // Calculate the boundary conditions + amrex::GpuArray bc_values = CalculatePoissonBoundaryPotentials( + Nc, Nv, bandgap, affinity, q, kb, T, + donor_doping, acceptor_doping, + Phi_Bc_lo, + Phi_Bc_hi + ); + + //bc_values = {0.0, 0.0}; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - if(k < 0) { - Phi(i,j,k) = Phi_Bc_lo; - } else if(k >= n_cell[2]){ - 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 phi_m = use_work_function ? metal_work_function : phi_ref; //in eV When not used, applied voltgae is set as the potential on the metal interface - Phi(i,j,k) = Phi_Bc_hi - (phi_m - phi_ref); - } + if(i == 0 && j == 0 && k == 0)amrex::Print() << "bc_values = " << bc_values[0] << ", " << bc_values[1] << "\n"; + + + // Boundary condition for the lower z-face (k=0 in problem domain, so k < 0 in ghost cells) + if (k < 0) { + if (mask(i,j,0) == 3.0) { // lo_z touches p-type + + Phi(i,j,k) = bc_values[1]; + if(i == 0 && j == 0)amrex::Print() << "lo z : setting up dirichlet BC at p-type contact with Phi_Bc_lo = " << Phi_Bc_lo << ", Phi = " << Phi(i,j,k) << "\n"; + + } else if (mask(i,j,0) == 4.0) { // lo_z touches n-type + + Phi(i,j,k) = bc_values[0]; + if(i == 0 && j == 0)amrex::Print() << "lo z : setting up dirichlet BC at n-type contact with Phi_Bc_lo = " << Phi_Bc_lo << ", Phi = " << Phi(i,j,k) << "\n"; + + } else { // lo_z touches insulator or intrinsic SC or metal + + if(i == 0 && j == 0)amrex::Print() << "lo z : setting up dirichlet BC at insulator or intrinsic SC contact " << "\n"; + Phi(i,j,k) = Phi_Bc_lo - (phi_m_V - phi_ref_V); + } + } + + // Boundary condition for the upper z-face (k=n_cell[2] in problem domain, so k >= n_cell[2] in ghost cells) + if (k >= n_cell[2]) { + if (mask(i,j,n_cell[2]-1) == 3.0) { // hi_z touches p-type + + Phi(i,j,k) = bc_values[1]; + if(i == 0 && j == 0)amrex::Print() << "hi z : setting up dirichlet BC at p-type contact with Phi_Bc_hi = " << Phi_Bc_hi << ", Phi = " << Phi(i,j,k) << "\n"; + + } else if (mask(i,j,n_cell[2]-1) == 4.0) { // hi_z touches n-type + + Phi(i,j,k) = bc_values[0]; + if(i == 0 && j == 0)amrex::Print() << "hi z : setting up dirichlet BC at n-type contact with Phi_Bc_hi = " << Phi_Bc_hi << ", Phi = " << Phi(i,j,k) << "\n"; + + } else { // hi_z touches insulator or intrinsic SC or metal + + if(i == 0 && j == 0)amrex::Print() << "hi z : setting up dirichlet BC at insulator or intrinsic SC contact " << "\n"; + Phi(i,j,k) = Phi_Bc_hi - (phi_m_V - phi_ref_V); + } + } }); } PoissonPhi.FillBoundary(geom.periodicity()); @@ -564,6 +696,7 @@ void SetupMLMG(std::unique_ptr& pMLMG, std::array,2>& LinOpBCType_2d, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + MultiFab& MaterialMask, c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info) { auto& rGprop = rFerroX.get_GeometryProperties(); @@ -586,18 +719,24 @@ void SetupMLMG(std::unique_ptr& pMLMG, p_mlabec->setDomainBC(LinOpBCType_2d[0], LinOpBCType_2d[1]); - if(some_constant_inhomogeneous_boundaries) - { - Fill_Constant_Inhomogeneous_Boundaries(rFerroX, PoissonPhi); - } - if(some_functionbased_inhomogeneous_boundaries) - { - Fill_FunctionBased_Inhomogeneous_Boundaries(rFerroX, PoissonPhi, time); - } - PoissonPhi.FillBoundary(geom.periodicity()); + SetPoissonBC(rFerroX, LinOpBCType_2d, all_homogeneous_boundaries, some_functionbased_inhomogeneous_boundaries, some_constant_inhomogeneous_boundaries); + + //if(some_constant_inhomogeneous_boundaries) + //{ + // Fill_Constant_Inhomogeneous_Boundaries(rFerroX, PoissonPhi); + //} + //if(some_functionbased_inhomogeneous_boundaries) + //{ + // Fill_FunctionBased_Inhomogeneous_Boundaries(rFerroX, PoissonPhi, time); + //} + //PoissonPhi.FillBoundary(geom.periodicity()); + // For now only this option implements Ohmic contacts and metal work function // set Dirichlet BC by reading in the ghost cell values - SetPhiBC_z(PoissonPhi, n_cell, geom); + //if(some_constant_inhomogeneous_boundaries){ + SetPhiBC_z(PoissonPhi, MaterialMask, n_cell, geom); + //} + p_mlabec->setLevelBC(amrlev, &PoissonPhi); // (A*alpha_cc - B * div beta grad) phi = rhs @@ -616,6 +755,7 @@ void SetupMLMG(std::unique_ptr& pMLMG, std::array,2>& LinOpBCType_2d, const amrex::GpuArray& n_cell, std::array< MultiFab, AMREX_SPACEDIM >& beta_face, + MultiFab& MaterialMask, MultiFab& beta_cc, c_FerroX& rFerroX, MultiFab& PoissonPhi, amrex::Real& time, amrex::LPInfo& info) { @@ -641,18 +781,18 @@ void SetupMLMG(std::unique_ptr& pMLMG, // assign domain boundary conditions to the solver p_mlebabec->setDomainBC(LinOpBCType_2d[0], LinOpBCType_2d[1]); - if(some_constant_inhomogeneous_boundaries) - { - Fill_Constant_Inhomogeneous_Boundaries(rFerroX, PoissonPhi); - } - if(some_functionbased_inhomogeneous_boundaries) - { - Fill_FunctionBased_Inhomogeneous_Boundaries(rFerroX, PoissonPhi, time); - } - PoissonPhi.FillBoundary(geom.periodicity()); + //if(some_constant_inhomogeneous_boundaries) + //{ + // Fill_Constant_Inhomogeneous_Boundaries(rFerroX, PoissonPhi); + //} + //if(some_functionbased_inhomogeneous_boundaries) + //{ + // Fill_FunctionBased_Inhomogeneous_Boundaries(rFerroX, PoissonPhi, time); + //} + //PoissonPhi.FillBoundary(geom.periodicity()); // Set Dirichlet BC for Phi in z - SetPhiBC_z(PoissonPhi, n_cell, geom); + SetPhiBC_z(PoissonPhi, MaterialMask, n_cell, geom); p_mlebabec->setLevelBC(amrlev, &PoissonPhi); // (A*alpha_cc - B * div beta grad) phi = rhs @@ -675,26 +815,30 @@ void SetupMLMG(std::unique_ptr& pMLMG, } #endif -void ComputePhi_Rho(std::unique_ptr& pMLMG, +void ComputePhi_Rho_Equilibrium(std::unique_ptr& pMLMG, std::unique_ptr& p_mlabec, MultiFab& alpha_cc, MultiFab& PoissonRHS, MultiFab& PoissonPhi, MultiFab& PoissonPhi_Prev, MultiFab& PhiErr, - Array& P_old, + Array& P_old, MultiFab& rho, + Array& Jn, + Array& Jp, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, + 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& prob_lo, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { //Obtain self consisten Phi and rho - Real tol = 1.e-5; + Real tol = 1.e-3; Real err = 1.0; int iter = 0; bool contains_SC = false; @@ -705,7 +849,7 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, rho, e_den, p_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, Jn, Jp, rho, e_den, p_den, e_den_old, p_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); @@ -717,15 +861,95 @@ void ComputePhi_Rho(std::unique_ptr& pMLMG, //Poisson Solve pMLMG->solve({&PoissonPhi}, {&PoissonRHS}, 1.e-10, -1); - PoissonPhi.FillBoundary(geom.periodicity()); + + PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region ComputeRho(PoissonPhi, rho, e_den, p_den, MaterialMask); - + if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates err = 0.; } else { + // err = 0.; + + // Calculate Error + if (iter > 0){ + MultiFab::Copy(PhiErr, PoissonPhi, 0, 0, 1, 0); + MultiFab::Subtract(PhiErr, PoissonPhi_Prev, 0, 0, 1, 0); + err = PhiErr.norm1(0, geom.periodicity())/PoissonPhi.norm1(0, geom.periodicity()); + } + + //Copy PoissonPhi to PoissonPhi_Prev to calculate error at the next iteration + MultiFab::Copy(PoissonPhi_Prev, PoissonPhi, 0, 0, 1, 0); + + iter = iter + 1; + amrex::Print() << iter << " iterations :: err = " << err << std::endl; + if( iter > 20 ) amrex::Print() << "Failed to reach self consistency between Phi and Rho in 20 iterations!! " << std::endl; + } + } + + // amrex::Print() << "\n ========= Self-Consistent Initialization of Phi and Rho Done! ========== \n"<< iter << " iterations to obtain self consistent Phi with err = " << err << std::endl; +} + +void ComputePhi_Rho(std::unique_ptr& pMLMG, + std::unique_ptr& p_mlabec, + MultiFab& alpha_cc, + MultiFab& PoissonRHS, + MultiFab& PoissonPhi, + MultiFab& PoissonPhi_Prev, + MultiFab& PhiErr, + Array& P_old, + MultiFab& rho, + Array& Jn, + Array& Jp, + 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& prob_lo, + const amrex::GpuArray& prob_hi) + +{ +//Obtain self consisten Phi and rho + Real tol = 1.e-3; + Real err = 1.0; + int iter = 0; + bool contains_SC = false; + FerroX_Util::Contains_sc(MaterialMask, contains_SC); + + while(err > tol){ + + //Compute RHS of Poisson equation + ComputePoissonRHS(PoissonRHS, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); + + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, Jn, Jp, rho, e_den, p_den, e_den_old, p_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + + + p_mlabec->setACoeffs(0, alpha_cc); + + //Initial guess for phi + PoissonPhi.setVal(0.); + + //Poisson Solve + pMLMG->solve({&PoissonPhi}, {&PoissonRHS}, 1.e-10, -1); + + PoissonPhi.FillBoundary(geom.periodicity()); + + // Calculate rho from Phi in SC region + //ComputeRho(PoissonPhi, rho, e_den, p_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, rho, Jn, Jp, e_den, p_den, e_den_old, p_den_old, MaterialMask, geom); + + if (contains_SC == 0) { + // no semiconductor region; set error to zero so the while loop terminates + err = 0.; + } else { + // err = 0.; // Calculate Error if (iter > 0){ @@ -754,14 +978,18 @@ void ComputePhi_Rho_EB(std::unique_ptr& pMLMG, MultiFab& PoissonPhi, MultiFab& PoissonPhi_Prev, MultiFab& PhiErr, - Array& P_old, + Array& P_old, MultiFab& rho, + Array& Jn, + Array& Jp, MultiFab& e_den, MultiFab& p_den, - MultiFab& MaterialMask, + 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& prob_lo, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { @@ -777,9 +1005,9 @@ void ComputePhi_Rho_EB(std::unique_ptr& pMLMG, //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, rho, e_den, p_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, Jn, Jp, rho, e_den, p_den, e_den_old, p_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); p_mlebabec->setACoeffs(0, alpha_cc); @@ -792,7 +1020,8 @@ void ComputePhi_Rho_EB(std::unique_ptr& pMLMG, PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, rho, e_den, p_den, MaterialMask); + //ComputeRho(PoissonPhi, rho, e_den, p_den, MaterialMask); + ComputeRho_DriftDiffusion(PoissonPhi, rho, Jn, Jp, e_den, p_den, e_den_old, p_den_old, MaterialMask, geom); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates diff --git a/Source/Solver/Initialization.H b/Source/Solver/Initialization.H index 8a929d5..722259a 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 447adf4..33f1e0c 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -1,13 +1,16 @@ #include "Initialization.H" #include "Utils/eXstaticUtils/eXstaticUtil.H" #include "../../Utils/SelectWarpXUtils/WarpXUtil.H" - +#define MATERIAL_UNKNOWN -1.0 // A unique value for undefined/overlapping regions + // INITIALIZE rho in SC region void InitializePandRho(Array &P_old, MultiFab& Gamma, 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, @@ -130,76 +133,446 @@ void InitializePandRho(Array &P_old, pOld_q(i,j,k) = 0.0; } }); - // Calculate charge density from Phi, Nc, Nv, Ec, and Ev + } + + for (int i = 0; i < 3; i++){ + // fill periodic ghost cells + P_old[i].FillBoundary(geom.periodicity()); + } + + // 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); + // loop over boxes + for (MFIter mfi(rho); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); 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); + //const Array4& charge_den_arr = rho.array(mfi); const Array4& acceptor_den_arr = acceptor_den.array(mfi); const Array4& donor_den_arr = donor_den.array(mfi); + const Array4& mask = MaterialMask.array(mfi); + + // extract dx from the geometry object + GpuArray dx = geom.CellSizeArray(); + /* + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + amrex::Real Na_val, Nd_val; // Use temporary values for N_A and N_D for this cell + amrex::Real initial_n, initial_p; + + Real cell_z = prob_lo[2] + (k+0.5) * dx[2]; + + const amrex::Real z_junction = 0.5e-6; // Center of the p-n junction + const amrex::Real junction_width = 0.05e-6; // Characteristic width of the transition (adjust as needed) + + // SC region (mask >= 2.0 indicates semiconductor) + if (mask(i,j,k) >= 2.0) { + + // Smooth doping profile for p-n junction + Nd_val = 0.5 * donor_doping * (1.0 - std::tanh((cell_z - z_junction) / junction_width)); + Na_val = 0.5 * acceptor_doping * (1.0 + std::tanh((cell_z - z_junction) / junction_width)); + + // If explicitly an intrinsic region, override doping (optional, depends on model) + if (mask(i,j,k) == 2.0) { + Na_val = 0.0; + Nd_val = 0.0; + } + + // --- MODIFIED INITIALIZATION FOR CARRIER CONCENTRATIONS --- + amrex::Real ni = intrinsic_carrier_concentration; + amrex::Real net_doping = Nd_val - Na_val; + + // Simple approach: Assume majority carriers are equal to net doping, + // minority carriers follow mass action. This works for neutral regions. + // For the depletion region, we will use a more explicit model. + + // This is a simplified depletion approximation for initial guess. + // It's still approximate and the solver will refine it. + if (net_doping > 0.0) { // Effectively n-type + initial_n = std::max(Nd_val - Na_val, ni); // Assume n ~ Nd - Na, but not less than ni + initial_p = ni * ni / initial_n; + } else if (net_doping < 0.0) { // Effectively p-type + initial_p = std::max(Na_val - Nd_val, ni); // Assume p ~ Na - Nd, but not less than ni + initial_n = ni * ni / initial_p; + } else { // Intrinsic or perfectly compensated + initial_n = ni; + initial_p = ni; + } + + // Now, the crucial part: if we are in the *anticipated* depletion region, + // we might want to override carrier concentrations to reflect charge. + // This is still an approximation for initial guess. + // A better way is to solve Poisson with initial potential guess. + + // Let's refine the initial_n and initial_p calculation to explicitly consider the + // charge density directly in the formula, rather than assuming local neutrality. + // This is typically done by solving the full drift-diffusion-Poisson system. + // For a *pure initialization of charge for Poisson*, you can often simplify. + + // Let's go back to the idea that in neutral regions, the charge is zero, + // and in the depletion region, it's roughly the ionized dopants. + // The tanh function *already* defines the doping profile. + // The most robust way is to provide a *potential* profile guess, + // then derive n and p from it, and then calculate charge. + + // For a good initial guess for Poisson, you often calculate: + // rho = q * (p - n + Nd+ - Na-) + // and then you solve Poisson for phi. + // + // The issue you're facing is that your *initial guess for n and p* already + // enforces p - n - Na + Nd = 0, which makes rho vanish. + // + // If you want to *see* charge at initialization, you must provide an initial + // guess for n and p that *does not* locally enforce neutrality in the depletion region. + + // REVISED INITIALIZATION FOR CHARGE_DEN_ARR: + // Let's try to directly define the charge density for a simplified + // depletion region, primarily driven by the dopants. + + // Assign initial carrier concentrations to MultiFabs + // These are just *guesses* to kick off the solver. + // The solver will find the true equilibrium n and p. + hole_den_arr(i,j,k) = initial_p; // Based on local neutrality in bulk + e_den_arr(i,j,k) = initial_n; // Based on local neutrality in bulk + + // Assign doping concentrations to MultiFabs + acceptor_den_arr(i,j,k) = Na_val; + donor_den_arr(i,j,k) = Nd_val; + + // Calculate the charge density for Poisson's RHS + // THIS IS THE CRITICAL PART FOR SEEING INITIAL CHARGE + // The charge density *is* primarily due to the fixed dopants in depletion. + // Mobile carriers will be very low in the depletion region. + + // For an initial guess, one common way is to assume full ionization of dopants + // and negligible mobile carriers *in the regions that will become depleted*. + // However, defining these regions precisely without a potential guess is tricky. + + // Let's stick with the definition: charge_den_arr = q * (p - n + Nd+ - Na-) + // And understand *why* it's small. + // Your current `initial_n` and `initial_p` are *designed* to make (p - n + Nd - Na) small. + + // A common strategy for initialization of a *depletion region* is to + // assume that n and p are very low in the "middle" region and then + // let the dopants define the charge. + + // Let's try a different approach to defining initial_n and initial_p + // that allows for charge in the depletion region. + // This still relies on a *conceptual* depletion region. + + amrex::Real phi_bi = (kb * T / q) * std::log((donor_doping * acceptor_doping) / (ni * ni)); // Built-in potential approximation + + // A linear or tanh-smoothed potential profile across the junction + // This is a rough guess for the potential. + amrex::Real initial_potential; + if (cell_z < z_junction) { // N-side + initial_potential = 0.5 * phi_bi * (1.0 - std::tanh((z_junction - cell_z) / (junction_width))); + } else { // P-side + initial_potential = -0.5 * phi_bi * (1.0 - std::tanh((cell_z - z_junction) / (junction_width))); + } + // Adjust potential to be relative to some reference (e.g., intrinsic Fermi level) + // For simplicity, let's assume the N-side is higher potential, P-side lower. + // Potential on N side should be higher than P side by phi_bi. + // So, potential goes from +phi_bi/2 to -phi_bi/2 across the junction. + + // Let's refine the initial potential guess based on typically assumed potential for an ideal junction + // N-side (far from junction) potential ~ 0 (or some reference) + // P-side (far from junction) potential ~ -phi_bi + initial_potential = -0.5 * phi_bi * (1.0 + std::tanh((cell_z - z_junction) / junction_width)); + // This will go from near 0 (n-side) to near -phi_bi (p-side) + // The sign convention for potential can vary; often e.g. n-side is +phi_bi/2 and p-side -phi_bi/2 + + // Let's try another common convention: potential difference from intrinsic level + // Assuming EF_n - E_i = kT ln(Nd/ni) and EF_p - E_i = -kT ln(Na/ni) + // and for a p-n junction at equilibrium, EF is flat. + // So, the potential should track the difference between Ei and EF. + // E_i - E_F = -q*potential + // Thus, potential = (E_F - E_i) / q + + // Let V_T = k_B * T / q; + // Potential on N-side (far from junction) should be approx V_T * ln(donor_doping / ni) + // Potential on P-side (far from junction) should be approx -V_T * ln(acceptor_doping / ni) + amrex::Real V_T = kb * T / q; + amrex::Real phi_n_bulk = V_T * std::log(donor_doping / ni); + amrex::Real phi_p_bulk = -V_T * std::log(acceptor_doping / ni); + + // This combines the potential profile with the smooth doping. + // The tanh term will make potential transition from phi_n_bulk to phi_p_bulk. + initial_potential = 0.5 * (phi_n_bulk + phi_p_bulk) - 0.5 * (phi_n_bulk - phi_p_bulk) * std::tanh((cell_z - z_junction) / junction_width); + + // Store this initial potential (if you have a potential MultiFab) + // potential_arr(i,j,k) = initial_potential; // You'd need a potential MultiFab + + // Now calculate n and p from this initial potential + initial_n = ni * std::exp(initial_potential / V_T); + initial_p = ni * std::exp(-initial_potential / V_T); + + // Assign carrier concentrations to MultiFabs + hole_den_arr(i,j,k) = initial_p; + e_den_arr(i,j,k) = initial_n; + + // Assign doping concentrations to MultiFabs + acceptor_den_arr(i,j,k) = Na_val; + donor_den_arr(i,j,k) = Nd_val; + + } else { // Non-semiconductor regions + hole_den_arr(i,j,k) = 0.0; + e_den_arr(i,j,k) = 0.0; + acceptor_den_arr(i,j,k) = 0.0; + donor_den_arr(i,j,k) = 0.0; + } + + // Calculate the charge density for Poisson's RHS + // q is the elementary charge (e.g., 1.602e-19 C) + // 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)); + }); + /* + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + amrex::Real Na_val, Nd_val; // Use temporary values for N_A and N_D for this cell + amrex::Real initial_n, initial_p; + + // Get the physical coordinate of the cell center in z + Real cell_z = prob_lo[2] + (k+0.5) * dx[2]; + + // Define junction parameters (these should be passed in or be global constants) + const amrex::Real z_junction = 0.5e-6; // Center of the p-n junction + const amrex::Real junction_width = 0.05e-6; // Characteristic width of the transition (adjust as needed) + + // SC region (mask >= 2.0 indicates semiconductor) + if (mask(i,j,k) >= 2.0) { + + // Apply smooth doping profile for p-n junction + // Using tanh for a smooth transition from n-type to p-type + // At z_junction, Na_val and Nd_val should be roughly equal (and half of max doping) + // For z < z_junction (n-type side), Nd_val should be high, Na_val low + // For z > z_junction (p-type side), Na_val should be high, Nd_val low + + // Nd_val will be high for z < z_junction + Nd_val = 0.5 * donor_doping * (1.0 - std::tanh((cell_z - z_junction) / junction_width)); + // Na_val will be high for z > z_junction + Na_val = 0.5 * acceptor_doping * (1.0 + std::tanh((cell_z - z_junction) / junction_width)); + + if (mask(i,j,k) == 2.0) { // If it's explicitly marked as intrinsic in the mask, keep it intrinsic + Na_val = 0.0; + Nd_val = 0.0; + } + // Calculate initial carrier concentrations based on local doping. + // This assumes charge neutrality and full ionization at equilibrium. + amrex::Real net_doping = Nd_val - Na_val; + amrex::Real ni = intrinsic_carrier_concentration; + + if (net_doping > 0) { // n-type effective region + initial_n = 0.5 * (net_doping + std::sqrt(net_doping * net_doping + 4.0 * ni * ni)); + initial_p = ni * ni / initial_n; + } else if (net_doping < 0) { // p-type effective region + initial_p = 0.5 * (-net_doping + std::sqrt(net_doping * net_doping + 4.0 * ni * ni)); + initial_n = ni * ni / initial_p; + } else { // intrinsic (or very close to intrinsic) + initial_n = ni; + initial_p = ni; + } + + // Assign initial carrier concentrations to MultiFabs + hole_den_arr(i,j,k) = initial_p; + e_den_arr(i,j,k) = initial_n; + + // Assign doping concentrations to MultiFabs + acceptor_den_arr(i,j,k) = Na_val; + donor_den_arr(i,j,k) = Nd_val; + + } else { // Non-semiconductor regions (e.g., oxide, metal contacts if you have them) + // Set carrier and doping densities to zero outside SC region + hole_den_arr(i,j,k) = 0.0; + e_den_arr(i,j,k) = 0.0; + acceptor_den_arr(i,j,k) = 0.0; + donor_den_arr(i,j,k) = 0.0; + // charge_den_arr(i,j,k) = 0.0; + } + + // Calculate the charge density for Poisson's RHS + // q is the elementary charge (e.g., 1.602e-19 C) + //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(i == 32 && j == 32 && k == 32) amrex::Print() << "hole_den_arr = " << hole_den_arr(i,j,k) << "\n" << "e_den_arr = " << e_den_arr(i,j,k) << "\n" << "acceptor_den_arr = " << acceptor_den_arr(i,j,k) << "\n" << "donor_den_arr = " << donor_den_arr(i,j,k) << "\n" << "charge_den_arr = " << charge_den_arr(i,j,k) << "\n"; + }); +*/ amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + amrex::Real Na_val, Nd_val; // Use temporary values for N_A and N_D for this cell + amrex::Real initial_n, initial_p; + + // SC region (mask >= 2.0 indicates semiconductor) + if (mask(i,j,k) >= 2.0) { + + if (mask(i,j,k) == 2.0) { // intrinsic + Na_val = 0.0; + Nd_val = 0.0; + initial_n = intrinsic_carrier_concentration; + initial_p = intrinsic_carrier_concentration; + } else if (mask(i,j,k) == 3.0) { // p-type + Na_val = acceptor_doping; + Nd_val = 0.0; + // In p-type, p is majority, n is minority + initial_p = acceptor_doping; // Assume full ionization and charge neutrality + initial_n = intrinsic_carrier_concentration * intrinsic_carrier_concentration / initial_p; + } else if (mask(i,j,k) == 4.0) { // n-type + Na_val = 0.0; + Nd_val = donor_doping; + // In n-type, n is majority, p is minority + initial_n = donor_doping; // Assume full ionization and charge neutrality + initial_p = intrinsic_carrier_concentration * intrinsic_carrier_concentration / initial_n; + } - //SC region - if (mask(i,j,k) >= 2.0) { + // Assign initial carrier concentrations to MultiFabs + hole_den_arr(i,j,k) = initial_p; + e_den_arr(i,j,k) = initial_n; - 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; - } + // Assign doping concentrations to MultiFabs + acceptor_den_arr(i,j,k) = Na_val; + donor_den_arr(i,j,k) = Nd_val; - 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 { // Non-semiconductor regions (e.g., oxide, metal contacts if you have them) + // Set carrier and doping densities to zero outside SC region + hole_den_arr(i,j,k) = 0.0; + e_den_arr(i,j,k) = 0.0; + acceptor_den_arr(i,j,k) = 0.0; + donor_den_arr(i,j,k) = 0.0; + } + // Calculate the charge density for Poisson's RHS + // q is the elementary charge (e.g., 1.602e-19 C) + // 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(i == 32 && j == 32 && k == 32) amrex::Print() << "hole_den_arr = " << hole_den_arr(i,j,k) << "\n" << "e_den_arr = " << e_den_arr(i,j,k) << "\n" << "acceptor_den_arr = " << acceptor_den_arr(i,j,k) << "\n" << "donor_den_arr = " << donor_den_arr(i,j,k) << "\n" << "charge_den_arr = " << charge_den_arr(i,j,k) << "\n"; }); +// */ } - for (int i = 0; i < 3; i++){ - // fill periodic ghost cells - P_old[i].FillBoundary(geom.periodicity()); - } + // Fill boundaries for all MultiFabs + e_den.FillBoundary(geom.periodicity()); + p_den.FillBoundary(geom.periodicity()); + acceptor_den.FillBoundary(geom.periodicity()); + donor_den.FillBoundary(geom.periodicity()); + + for (MFIter mfi(rho); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); // Or mfi.tilebox() + //const Box& bx = mfi.growntilebox(1); + + 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); + const Array4& acceptor_den_arr = acceptor_den.array(mfi); + const Array4& donor_den_arr = donor_den.array(mfi); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + 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)); + //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)); + }); + } + rho.FillBoundary(geom.periodicity()); } // create a mask filled with integers to idetify different material types -void InitializeMaterialMask(MultiFab& MaterialMask, - const Geometry& geom, - const amrex::GpuArray& prob_lo, +//void InitializeMaterialMask(MultiFab& MaterialMask, +// const Geometry& geom, +// const amrex::GpuArray& prob_lo, +// const amrex::GpuArray& prob_hi) +//{ +// // loop over boxes +// for (MFIter mfi(MaterialMask); mfi.isValid(); ++mfi) +// { +// const Box& bx = mfi.growntilebox(MaterialMask.nGrow()); +// // extract dx from the geometry object +// GpuArray dx = geom.CellSizeArray(); +// +// const Array4& mask = MaterialMask.array(mfi); +// +// +// amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) +// { +// Real x = prob_lo[0] + (i+0.5) * dx[0]; +// Real y = prob_lo[1] + (j+0.5) * dx[1]; +// Real z = prob_lo[2] + (k+0.5) * dx[2]; +// +// //FE:0, DE:1, Source/Drain:2, p_type:3, n_type:4 +// if (x <= FE_hi[0] && x >= FE_lo[0] && y <= FE_hi[1] && y >= FE_lo[1] && z <= FE_hi[2] && z >= FE_lo[2]) { +// mask(i,j,k) = 0.; +// } else if (x <= DE_hi[0] && x >= DE_lo[0] && y <= DE_hi[1] && y >= DE_lo[1] && z <= DE_hi[2] && z >= DE_lo[2]) { +// mask(i,j,k) = 1.; +// //intrinsic +// } else if (x <= SC_hi[0] && x >= SC_lo[0] && y <= SC_hi[1] && y >= SC_lo[1] && z <= SC_hi[2] && z >= SC_lo[2]) { +// mask(i,j,k) = 2.; +// //p_type +// } else if (x <= p_type_hi[0] && x >= p_type_lo[0] && y <= p_type_hi[1] && y >= p_type_lo[1] && z <= p_type_hi[2] && z >= p_type_lo[2]){ +// mask(i,j,k) = 3.; +// //n_type +// } else if (x <= n_type_hi[0] && x >= n_type_lo[0] && y <= n_type_hi[1] && y >= n_type_lo[1] && z <= n_type_hi[2] && z >= n_type_lo[2]){ +// mask(i,j,k) = 4.; +// } else { +// mask(i,j,k) = 1.; //spacer is DE +// } +// }); +// } +// MaterialMask.FillBoundary(geom.periodicity()); +//} +// +void InitializeMaterialMask(MultiFab& MaterialMask, + const Geometry& geom, + const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { - // loop over boxes + for (MFIter mfi(MaterialMask); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(MaterialMask.nGrow()); - // extract dx from the geometry object GpuArray dx = geom.CellSizeArray(); - const Array4& mask = MaterialMask.array(mfi); - amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Real x = prob_lo[0] + (i+0.5) * dx[0]; - Real y = prob_lo[1] + (j+0.5) * dx[1]; - Real z = prob_lo[2] + (k+0.5) * dx[2]; - - //FE:0, DE:1, Source/Drain:2, Channel:3 - if (x <= FE_hi[0] && x >= FE_lo[0] && y <= FE_hi[1] && y >= FE_lo[1] && z <= FE_hi[2] && z >= FE_lo[2]) { - mask(i,j,k) = 0.; - } else if (x <= DE_hi[0] && x >= DE_lo[0] && y <= DE_hi[1] && y >= DE_lo[1] && z <= DE_hi[2] && z >= DE_lo[2]) { - mask(i,j,k) = 1.; - } else if (x <= SC_hi[0] && x >= SC_lo[0] && y <= SC_hi[1] && y >= SC_lo[1] && z <= SC_hi[2] && z >= SC_lo[2]) { - mask(i,j,k) = 2.; - 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 { - mask(i,j,k) = 1.; //spacer is DE - } + Real x = prob_lo[0] + (i+0.5) * dx[0]; + Real y = prob_lo[1] + (j+0.5) * dx[1]; + Real z = prob_lo[2] + (k+0.5) * dx[2]; + + if (x >= n_type_lo[0] && x <= n_type_hi[0] && + y >= n_type_lo[1] && y <= n_type_hi[1] && + z >= n_type_lo[2] && z <= n_type_hi[2]) + { + mask(i,j,k) = 4.; // n_type + } + else if (x >= p_type_lo[0] && x <= p_type_hi[0] && + y >= p_type_lo[1] && y <= p_type_hi[1] && + z >= p_type_lo[2] && z <= p_type_hi[2]) + { + mask(i,j,k) = 3.; // p_type + } + else if (x >= FE_lo[0] && x <= FE_hi[0] && + y >= FE_lo[1] && y <= FE_hi[1] && + z >= FE_lo[2] && z <= FE_hi[2]) + { + mask(i,j,k) = 0.; // FE + } + else if (x >= DE_lo[0] && x <= DE_hi[0] && + y >= DE_lo[1] && y <= DE_hi[1] && + z >= DE_lo[2] && z <= DE_hi[2]) + { + mask(i,j,k) = 1.; // DE + } + else if (x >= SC_lo[0] && x <= SC_hi[0] && + y >= SC_lo[1] && y <= SC_hi[1] && + z >= SC_lo[2] && z <= SC_hi[2]) + { + mask(i,j,k) = 2.; // SC (Source/Drain) + } + else + { + mask(i,j,k) = 1.; // Default: Assign to DE (Dielectric) for any region not explicitly defined above + } + }); } MaterialMask.FillBoundary(geom.periodicity()); diff --git a/Source/main.cpp b/Source/main.cpp index a1d8daa..71345d3 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -132,6 +132,15 @@ void main_main (c_FerroX& rFerroX) E[dir].define(ba, dm, Ncomp, 0); } + //Drift-Diffusion + Array Jn; + Array Jp; + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) + { + Jn[dir].define(ba, dm, Ncomp, 1); + Jp[dir].define(ba, dm, Ncomp, 1); + } + MultiFab PoissonRHS(ba, dm, 1, 0); MultiFab PoissonPhi(ba, dm, 1, 1); MultiFab PoissonPhi_Old(ba, dm, 1, 1); @@ -143,8 +152,12 @@ void main_main (c_FerroX& rFerroX) MultiFab Ey(ba, dm, 1, 0); MultiFab Ez(ba, dm, 1, 0); - MultiFab hole_den(ba, dm, 1, 0); - MultiFab e_den(ba, dm, 1, 0); + MultiFab hole_den(ba, dm, 1, 1); + MultiFab e_den(ba, dm, 1, 1); + MultiFab acceptor_den(ba, dm, 1, 1); + MultiFab donor_den(ba, dm, 1, 1); + MultiFab hole_den_old(ba, dm, 1, 1); + MultiFab e_den_old(ba, dm, 1, 1); MultiFab charge_den(ba, dm, 1, 0); MultiFab MaterialMask(ba, dm, 1, 1); MultiFab tphaseMask(ba, dm, 1, 1); @@ -161,10 +174,16 @@ void main_main (c_FerroX& rFerroX) GL_rhs_pre[dir].setVal(0.); GL_rhs_avg[dir].setVal(0.); E[dir].setVal(0.); + Jn[dir].setVal(0.); + Jp[dir].setVal(0.); } e_den.setVal(0.); hole_den.setVal(0.); + e_den_old.setVal(0.); + hole_den_old.setVal(0.); + acceptor_den.setVal(0.); + donor_den.setVal(0.); PoissonPhi.setVal(0.); PoissonRHS.setVal(0.); tphaseMask.setVal(0.); @@ -194,6 +213,7 @@ void main_main (c_FerroX& rFerroX) // coefficients for solver MultiFab alpha_cc(ba, dm, 1, 0); + alpha_cc.setVal(0.); MultiFab beta_cc(ba, dm, 1, 1); std::array< MultiFab, AMREX_SPACEDIM > beta_face; AMREX_D_TERM(beta_face[0].define(convert(ba,IntVect(AMREX_D_DECL(1,0,0))), dm, 1, 0);, @@ -213,23 +233,27 @@ void main_main (c_FerroX& rFerroX) int linop_maxorder = 2; int amrlev = 0; //refers to the setcoarsest level of the solve - SetupMLMG(pMLMG, p_mlabec, LinOpBCType_2d, n_cell, beta_face, rFerroX, PoissonPhi, time, info); + SetupMLMG(pMLMG, p_mlabec, LinOpBCType_2d, n_cell, beta_face, MaterialMask, rFerroX, PoissonPhi, time, info); #ifdef AMREX_USE_EB std::unique_ptr p_mlebabec; - SetupMLMG_EB(pMLMG, p_mlebabec, LinOpBCType_2d, n_cell, beta_face, beta_cc, rFerroX, PoissonPhi, time, info); + SetupMLMG_EB(pMLMG, p_mlebabec, LinOpBCType_2d, n_cell, beta_face, MaterialMask, beta_cc, rFerroX, PoissonPhi, time, info); #endif // 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 + ComputePhi_Rho_Equilibrium(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, + P_old, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, + angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + // Write a plotfile of the initial data if plot_int > 0 if (plot_int > 0) { int plt_step = 0; - WritePlotfile(rFerroX, PoissonPhi, PoissonRHS, P_old, E, hole_den, e_den, charge_den, beta_cc, + WritePlotfile(rFerroX, PoissonPhi, PoissonRHS, P_old, E, Jn, Jp, hole_den, e_den, acceptor_den, donor_den, charge_den, beta_cc, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, Phidiff, geom, time, plt_step); } @@ -266,11 +290,11 @@ void main_main (c_FerroX& rFerroX) #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - P_old, charge_den, e_den, hole_den, MaterialMask, + P_old, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - P_old, charge_den, e_den, hole_den, MaterialMask, + P_old, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #endif @@ -309,11 +333,11 @@ void main_main (c_FerroX& rFerroX) #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - P_new_pre, charge_den, e_den, hole_den, MaterialMask, + P_new_pre, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - P_new_pre, charge_den, e_den, hole_den, MaterialMask, + P_new_pre, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #endif @@ -380,11 +404,11 @@ void main_main (c_FerroX& rFerroX) #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - ar_state, charge_den, e_den, hole_den, MaterialMask, + ar_state, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - ar_state, charge_den, e_den, hole_den, MaterialMask, + ar_state, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #endif @@ -453,11 +477,11 @@ void main_main (c_FerroX& rFerroX) #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - ar_state, charge_den, e_den, hole_den, MaterialMask, + ar_state, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - ar_state, charge_den, e_den, hole_den, MaterialMask, + ar_state, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #endif @@ -564,7 +588,7 @@ void main_main (c_FerroX& rFerroX) if (plot_int > 0 && (step%plot_int == 0 || step == steady_state_step)) { int plt_step = step; - WritePlotfile(rFerroX, PoissonPhi, PoissonRHS, P_old, E, hole_den, e_den, charge_den, beta_cc, + WritePlotfile(rFerroX, PoissonPhi, PoissonRHS, P_old, E, Jn, Jp, hole_den, e_den, acceptor_den, donor_den, charge_den, beta_cc, MaterialMask, tphaseMask, angle_alpha, angle_beta, angle_theta, Phidiff, geom, time, plt_step); } @@ -582,25 +606,33 @@ void main_main (c_FerroX& rFerroX) amrex::Print() << "step = " << step << ", Phi_Bc_hi = " << Phi_Bc_hi << ", num_Vapp = " << num_Vapp << ", sign = " << sign << std::endl; // Set Dirichlet BC for Phi in z - SetPhiBC_z(PoissonPhi, n_cell, geom); + SetPhiBC_z(PoissonPhi, MaterialMask, n_cell, geom); // set Dirichlet BC by reading in the ghost cell values #ifdef AMREX_USE_EB p_mlebabec->setLevelBC(amrlev, &PoissonPhi); + + if(rGprop.pEB->specify_inhomogeneous_dirichlet == 0) + { + p_mlebabec->setEBHomogDirichlet(amrlev, beta_cc); + } + else + { + p_mlebabec->setEBDirichlet(amrlev, *rGprop.pEB->p_surf_soln_union, beta_cc); + } #else p_mlabec->setLevelBC(amrlev, &PoissonPhi); #endif #ifdef AMREX_USE_EB ComputePhi_Rho_EB(pMLMG, p_mlebabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - P_old, charge_den, e_den, hole_den, MaterialMask, + P_old, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #else ComputePhi_Rho(pMLMG, p_mlabec, alpha_cc, PoissonRHS, PoissonPhi, PoissonPhi_Prev, PhiErr, - P_old, charge_den, e_den, hole_den, MaterialMask, + P_old, charge_den, Jn, Jp, e_den, hole_den, e_den_old, hole_den_old, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); #endif - }//end inc_step if (voltage_sweep == 0 && step == steady_state_step) {