diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 91544d5..257beb9 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -212,12 +212,15 @@ 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; 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; @@ -358,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/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 f19f18b..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) { @@ -47,7 +44,12 @@ void ComputeRho(MultiFab& PoissonPhi, //Following: http://dx.doi.org/10.1063/1.4825209 - amrex::Real Ef = 0.0; + amrex::Real Ef; + if (mask(i,j,k) == 2.0){ + Ef = 0.0; + } else { + Ef = -q*Phi_Bc_hi; + } amrex::Real Eg = bandgap; amrex::Real Chi = affinity; amrex::Real phi_ref = Chi + 0.5*Eg + 0.5*kb*T*log(Nc/Nv)/q; @@ -62,8 +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 = 4.0; - amrex::Real g_D = 2.0; + amrex::Real g_A = 4.0; + amrex::Real g_D = 2.0; amrex::Real Ea = acceptor_ionization_energy; amrex::Real Ed = donor_ionization_energy; @@ -76,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)); @@ -99,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/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.H b/Source/Solver/ElectrostaticSolver.H index c8da3ce..6deba04 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, @@ -61,4 +63,4 @@ void SetPoissonBC(c_FerroX& rFerroX, std::array &P_old, MultiFab& MaterialMask, const amrex::GpuArray& n_cell); diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index bda6a91..269e403 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" @@ -87,6 +88,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 +105,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); @@ -521,3 +524,56 @@ void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray PoissonPhi.FillBoundary(geom.periodicity()); } + +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 &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); + + 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.) { + if (prob_type == 1) { //2D + + 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]*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; + } + } + }); + } + +} diff --git a/Source/Solver/Initialization.H b/Source/Solver/Initialization.H index 8a929d5..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, @@ -26,5 +28,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..934ad1d 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -2,12 +2,25 @@ #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, 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, @@ -34,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(); } @@ -104,6 +124,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"); @@ -132,9 +158,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); @@ -148,14 +171,75 @@ 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 - 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::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); + + 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) ); + + 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; + } + + } + + 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 +466,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..52ad80e 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -119,8 +119,11 @@ 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); MultiFab angle_alpha(ba, dm, 1, 0); MultiFab angle_beta(ba, dm, 1, 0); MultiFab angle_theta(ba, dm, 1, 0); @@ -138,16 +141,20 @@ 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.); + 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); @@ -276,7 +283,10 @@ 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, n_cell); //Obtain self consisten Phi and rho Real tol = 1.e-5; @@ -288,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); @@ -305,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 @@ -352,7 +362,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 +395,7 @@ void main_main (c_FerroX& rFerroX) P_new_pre[i].FillBoundary(geom.periodicity()); } + SetNucleation(P_new_pre, NucleationMask, n_cell); /** * \brief dst = a*x + b*y */ @@ -409,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); @@ -427,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 @@ -480,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); @@ -498,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 @@ -640,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); @@ -658,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