diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 680074f..d8c435d 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -176,6 +176,9 @@ AMREX_GPU_MANAGED amrex::GpuArray FerroX::FE_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::SC_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Channel_hi; AMREX_GPU_MANAGED amrex::GpuArray FerroX::Channel_lo; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::Metal_hi; +AMREX_GPU_MANAGED amrex::GpuArray FerroX::Metal_lo; + // material parameters AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_0; @@ -184,6 +187,9 @@ AMREX_GPU_MANAGED amrex::Real FerroX::epsilonX_fe_tphase; AMREX_GPU_MANAGED amrex::Real FerroX::epsilonZ_fe; AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_de; AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_si; +AMREX_GPU_MANAGED amrex::Real FerroX::epsilon_metal; +AMREX_GPU_MANAGED amrex::Real FerroX::metal_screening_length; +AMREX_GPU_MANAGED amrex::Real FerroX::metal_thickness; AMREX_GPU_MANAGED amrex::Real FerroX::alpha; // alpha = 2*alpha_1 AMREX_GPU_MANAGED amrex::Real FerroX::beta; // beta = 4*alpha_11 AMREX_GPU_MANAGED amrex::Real FerroX::gamma; // gamma = 6*alpha_111 @@ -284,6 +290,16 @@ void InitializeFerroXNamespace(const amrex::GpuArray SC_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Channel_hi; extern AMREX_GPU_MANAGED amrex::GpuArray Channel_lo; + extern AMREX_GPU_MANAGED amrex::GpuArray Metal_hi; + extern AMREX_GPU_MANAGED amrex::GpuArray Metal_lo; + // material parameters extern AMREX_GPU_MANAGED amrex::Real epsilon_0; @@ -23,6 +26,9 @@ namespace FerroX { extern AMREX_GPU_MANAGED amrex::Real epsilonZ_fe; extern AMREX_GPU_MANAGED amrex::Real epsilon_de; extern AMREX_GPU_MANAGED amrex::Real epsilon_si; + extern AMREX_GPU_MANAGED amrex::Real epsilon_metal; + extern AMREX_GPU_MANAGED amrex::Real metal_screening_length; + extern AMREX_GPU_MANAGED amrex::Real metal_thickness; extern AMREX_GPU_MANAGED amrex::Real alpha; // alpha = 2*alpha_1 extern AMREX_GPU_MANAGED amrex::Real beta; // beta = 4*alpha_11 extern AMREX_GPU_MANAGED amrex::Real gamma; // gamma = 6*alpha_111 diff --git a/Source/Solver/ChargeDensity.H b/Source/Solver/ChargeDensity.H index af92798..ab29aec 100644 --- a/Source/Solver/ChargeDensity.H +++ b/Source/Solver/ChargeDensity.H @@ -6,9 +6,19 @@ using namespace amrex; using namespace FerroX; + void ComputeRho(MultiFab& PoissonPhi, + Array &P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - const MultiFab& MaterialMask); + const MultiFab& MaterialMask, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi); + +void Compute_P_Sum(const std::array& P, Real& sum); + +void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count); +void Compute_P_av(const std::array& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real &P_av_z); diff --git a/Source/Solver/ChargeDensity.cpp b/Source/Solver/ChargeDensity.cpp index 7c93c24..f7777a1 100644 --- a/Source/Solver/ChargeDensity.cpp +++ b/Source/Solver/ChargeDensity.cpp @@ -2,20 +2,51 @@ // Compute rho in SC region for given phi void ComputeRho(MultiFab& PoissonPhi, + Array &P_old, MultiFab& rho, MultiFab& e_den, MultiFab& p_den, - const MultiFab& MaterialMask) + const MultiFab& MaterialMask, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi) { + // Calculate average Pr. We take the average only over the FE region + Real average_P_r = 0.; + Real total_P_r = 0.; + Real FE_index_counter = 0.; + + Compute_P_av(P_old, total_P_r, MaterialMask, FE_index_counter, average_P_r); + + //Calculate integrated electrode charge (Qe) based on eq 13 of https://pubs.aip.org/aip/jap/article/44/8/3379/6486/Depolarization-fields-in-thin-ferroelectric-films + Real FE_thickness = FE_hi[2] - FE_lo[2]; + Real coth = (exp(2.*metal_thickness/metal_screening_length) + 1.0) / (exp(2.*metal_thickness/metal_screening_length) - 1.0); + Real csch = (2.*exp(metal_thickness/metal_screening_length)) / (exp(2.*metal_thickness/metal_screening_length) - 1.0); + Real numerator = 0.5 * FE_thickness * average_P_r / epsilonX_fe; + Real denominator = metal_screening_length/epsilon_de*(coth - csch) + FE_thickness/(2.*epsilonX_fe); + Real Qe = -numerator/denominator; + + amrex::Print() << "ChargeDen average_P_r = " << average_P_r << "\n"; + amrex::Print() << "ChargeDen numerator = " << numerator << "\n"; + amrex::Print() << "ChargeDen denominator = " << denominator << "\n"; + amrex::Print() << "ChargeDen Qe = " << Qe << "\n"; + // loop over boxes for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); + // extract dx from the geometry object + GpuArray dx = geom.CellSizeArray(); + // Calculate charge density from Phi, Nc, Nv, Ec, and Ev MultiFab acceptor_den(rho.boxArray(), rho.DistributionMap(), 1, 0); MultiFab donor_den(rho.boxArray(), rho.DistributionMap(), 1, 0); + const Array4 &pOld_p = P_old[0].array(mfi); + const Array4 &pOld_q = P_old[1].array(mfi); + const Array4 &pOld_r = P_old[2].array(mfi); + const Array4& hole_den_arr = p_den.array(mfi); const Array4& e_den_arr = e_den.array(mfi); const Array4& charge_den_arr = rho.array(mfi); @@ -24,45 +55,70 @@ void ComputeRho(MultiFab& PoissonPhi, const Array4& donor_den_arr = donor_den.array(mfi); const Array4& mask = MaterialMask.array(mfi); - + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - if (mask(i,j,k) >= 2.0) { + Real z = prob_lo[2] + (k+0.5) * dx[2]; + Real z_metal = 0.; - if(use_Fermi_Dirac == 1){ - //Fermi-Dirac - Real eta_n = q*(phi(i,j,k) - Ec)/(kb*T); - Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2))); - Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8)); - Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0); - - e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; - - Real eta_p = q*(Ev - phi(i,j,k))/(kb*T); - Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2))); - Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8)); - Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0); + if (mask(i,j,k) >= 2.0) { - hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; - } else { - //Maxwell-Boltzmann - Real n_0 = intrinsic_carrier_concentration; - Real p_0 = intrinsic_carrier_concentration; - hole_den_arr(i,j,k) = n_0*exp(-(q*phi(i,j,k))/(kb*T)); - e_den_arr(i,j,k) = p_0*exp(q*phi(i,j,k)/(kb*T)); - } - - //If in channel, set acceptor doping, else (Source/Drain) set donor doping - if (mask(i,j,k) == 3.0) { + if (mask(i,j,k) == 4.0) { //Metal + + //if(z <= FE_lo[2]){ + // z_metal = std::abs(FE_lo[2] - (k + 0.5) * dx[2]); + //} else if (z >= FE_hi[2]){ + // z_metal = std::abs((k + 0.5) * dx[2] - FE_hi[2]); + //} + // // amrex::Print() << "Qe = " << Qe << "\n"; + //charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); + + //Treat Metal as SC + Real n_0 = intrinsic_carrier_concentration; + Real p_0 = intrinsic_carrier_concentration; + hole_den_arr(i,j,k) = n_0*exp(-(q*phi(i,j,k))/(kb*T)); + e_den_arr(i,j,k) = p_0*exp(q*phi(i,j,k)/(kb*T)); acceptor_den_arr(i,j,k) = acceptor_doping; - donor_den_arr(i,j,k) = 0.0; - } else { // Source / Drain - acceptor_den_arr(i,j,k) = 0.0; donor_den_arr(i,j,k) = donor_doping; - } + charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); + + } else { - charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); + if(use_Fermi_Dirac == 1){ + //Fermi-Dirac + Real eta_n = q*(phi(i,j,k) - Ec)/(kb*T); + Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2))); + Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8)); + Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0); + + e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; + + Real eta_p = q*(Ev - phi(i,j,k))/(kb*T); + Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2))); + Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8)); + Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0); + + hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; + } else { + //Maxwell-Boltzmann + Real n_0 = intrinsic_carrier_concentration; + Real p_0 = intrinsic_carrier_concentration; + hole_den_arr(i,j,k) = n_0*exp(-(q*phi(i,j,k))/(kb*T)); + e_den_arr(i,j,k) = p_0*exp(q*phi(i,j,k)/(kb*T)); + } + + //If in channel, set acceptor doping, else (Source/Drain) set donor doping + if (mask(i,j,k) == 3.0) { + acceptor_den_arr(i,j,k) = acceptor_doping; + donor_den_arr(i,j,k) = 0.0; + } else if (mask(i,j,k) == 2.0){ // Source / Drain + acceptor_den_arr(i,j,k) = 0.0; + donor_den_arr(i,j,k) = donor_doping; + } + + charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); + } } else { @@ -73,3 +129,74 @@ void ComputeRho(MultiFab& PoissonPhi, } } +void Compute_P_Sum(const std::array& P, Real& sum) +{ + + // Initialize to zero + sum = 0.; + + ReduceOps reduce_op; + + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + for (MFIter mfi(P[2],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Box& bx_grid = mfi.validbox(); + + auto const& fab = P[2].array(mfi); + + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple + { + return {fab(i,j,k)}; + }); + } + + sum = amrex::get<0>(reduce_data.value()); + ParallelDescriptor::ReduceRealSum(sum); +} + + +void Compute_P_index_Sum(const MultiFab& MaterialMask, Real& count) +{ + + // Initialize to zero + count = 0.; + + ReduceOps reduce_op; + + ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + for (MFIter mfi(MaterialMask, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Box& bx_grid = mfi.validbox(); + + auto const& fab = MaterialMask.array(mfi); + + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple + { + if(fab(i,j,k) == 0.) { + return {1.}; + } else { + return {0.}; + } + + }); + } + + count = amrex::get<0>(reduce_data.value()); + ParallelDescriptor::ReduceRealSum(count); +} + +void Compute_P_av(const std::array& P, Real& sum, const MultiFab& MaterialMask, Real& count, Real& P_av_z) +{ + Compute_P_Sum(P, sum); + Compute_P_index_Sum(MaterialMask, count); + P_av_z = sum/count; +} + diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index a1cac53..f903d9f 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -102,8 +102,8 @@ void dF_dPhi(MultiFab& alpha_cc, PoissonPhi_plus_delta.plus(delta, 0, 1, 0); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi_plus_delta, rho, e_den, p_den, MaterialMask); - + ComputeRho(PoissonPhi_plus_delta, P_old, rho, e_den, p_den, MaterialMask, geom, prob_lo, prob_hi); + //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); @@ -231,7 +231,7 @@ void InitializePermittivity(std::array= 2.0){ + } else if (mask(i,j,k) == 2.0 || mask(i,j,k) == 3.0){ beta(i,j,k) = epsilon_si * epsilon_0; //SC layer + } else if (mask(i,j,k) == 4.0){ + beta(i,j,k) = epsilon_metal * epsilon_0; //Metal layer } else { beta(i,j,k) = epsilon_de * epsilon_0; //Spacer is same as DE } diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 82cf1a0..5e1b8de 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -1,4 +1,5 @@ #include "Initialization.H" +#include "ChargeDensity.H" #include "Utils/eXstaticUtils/eXstaticUtil.H" #include "../../Utils/SelectWarpXUtils/WarpXUtil.H" @@ -59,6 +60,25 @@ void InitializePandRho(Array &P_old, rngs[i] = amrex::Random(); // uniform [0,1] option } + Real average_P_r = 1.e-6; + //Real total_P_r = 0.; + //Real FE_index_counter = 0.; + + //Compute_P_av(P_old, total_P_r, MaterialMask, FE_index_counter, average_P_r); + + //Calculate integrated electrode charge (Qe) based on eq 13 of https://pubs.aip.org/aip/jap/article/44/8/3379/6486/Depolarization-fields-in-thin-ferroelectric-films + Real FE_thickness = FE_hi[2] - FE_lo[2]; + Real coth = (exp(2.*metal_thickness/metal_screening_length) + 1.0) / (exp(2.*metal_thickness/metal_screening_length) - 1.0); + Real csch = (2.*exp(metal_thickness/metal_screening_length)) / (exp(2.*metal_thickness/metal_screening_length) - 1.0); + Real numerator = 0.5 * FE_thickness * average_P_r / epsilonX_fe; + Real denominator = metal_screening_length/epsilon_de*(coth - csch + FE_thickness/(2.*epsilonX_fe)); + Real Qe = -numerator/denominator; + + //amrex::Print() << "average_P_r = " << average_P_r << "\n"; + //amrex::Print() << "numerator = " << numerator << "\n"; + //amrex::Print() << "denominator = " << denominator << "\n"; + //amrex::Print() << "Qe = " << Qe << "\n"; + // loop over boxes for (MFIter mfi(rho); mfi.isValid(); ++mfi) { @@ -131,45 +151,67 @@ void InitializePandRho(Array &P_old, amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real z = prob_lo[2] + (k+0.5) * dx[2]; + Real z_metal = 0.; + //SC region if (mask(i,j,k) >= 2.0) { - if(use_Fermi_Dirac == 1){ + if (mask(i,j,k) == 4.0) { //Metal + + //if(z <= FE_lo[2]){ + // z_metal = std::abs(FE_lo[2] - (k + 0.5) * dx[2]); + //} else if (z >= FE_hi[2]){ + // z_metal = std::abs((k + 0.5) * dx[2] - FE_hi[2]); + //} + ////if(Qe != 0.) amrex::Print() << "initialization : Qe = " << Qe << "\n"; + //charge_den_arr(i,j,k) = Qe/metal_screening_length*exp(-z_metal/metal_screening_length); + + //Treat Metal as Semiconductor + hole_den_arr(i,j,k) = intrinsic_carrier_concentration; + e_den_arr(i,j,k) = intrinsic_carrier_concentration; + acceptor_den_arr(i,j,k) = acceptor_doping; + donor_den_arr(i,j,k) = donor_doping; + charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); + + } else { + if(use_Fermi_Dirac == 1){ - //Approximate FD integral - Real Phi = 0.5*(Ec + Ev); //eV - Real eta_n = q*(Phi - Ec)/(kb*T); - Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2.0))); - Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8)); - Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0); + //Approximate FD integral + Real Phi = 0.5*(Ec + Ev); //eV + Real eta_n = q*(Phi - Ec)/(kb*T); + Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2.0))); + Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8)); + Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0); - e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; + e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; - Real eta_p = q*(Ev - Phi)/(kb*T); - Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2.0))); - Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8)); - Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0); + Real eta_p = q*(Ev - Phi)/(kb*T); + Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2.0))); + Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8)); + Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0); - hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; + hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; - } else { - - hole_den_arr(i,j,k) = intrinsic_carrier_concentration; - e_den_arr(i,j,k) = intrinsic_carrier_concentration; - + } else { + + hole_den_arr(i,j,k) = intrinsic_carrier_concentration; + e_den_arr(i,j,k) = intrinsic_carrier_concentration; + + } + + //If in channel, set acceptor doping, else (Source/Drain) set donor doping + if (mask(i,j,k) == 3.0) { + acceptor_den_arr(i,j,k) = acceptor_doping; + donor_den_arr(i,j,k) = 0.0; + } else { // Source / Drain + acceptor_den_arr(i,j,k) = 0.0; + donor_den_arr(i,j,k) = donor_doping; + } + charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); } } - - //If in channel, set acceptor doping, else (Source/Drain) set donor doping - if (mask(i,j,k) == 3.0) { - acceptor_den_arr(i,j,k) = acceptor_doping; - donor_den_arr(i,j,k) = 0.0; - } else { // Source / Drain - acceptor_den_arr(i,j,k) = 0.0; - donor_den_arr(i,j,k) = donor_doping; - } - charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); - + }); } for (int i = 0; i < 3; i++){ @@ -211,6 +253,8 @@ void InitializeMaterialMask(MultiFab& MaterialMask, if (x <= Channel_hi[0] && x >= Channel_lo[0] && y <= Channel_hi[1] && y >= Channel_lo[1] && z <= Channel_hi[2] && z >= Channel_lo[2]){ mask(i,j,k) = 3.; } + } else if (x <= Metal_hi[0] && x >= Metal_lo[0] && y <= Metal_hi[1] && y >= Metal_lo[1] && z <= Metal_hi[2] && z >= Metal_lo[2]) { + mask(i,j,k) = 4.; } else { mask(i,j,k) = 1.; //spacer is DE } diff --git a/Source/Utils/FerroXUtils/FerroXUtil.cpp b/Source/Utils/FerroXUtils/FerroXUtil.cpp index 0237195..32e614c 100644 --- a/Source/Utils/FerroXUtils/FerroXUtil.cpp +++ b/Source/Utils/FerroXUtils/FerroXUtil.cpp @@ -25,7 +25,7 @@ void FerroX_Util::Contains_sc(MultiFab& MaterialMask, bool& contains_SC) for (auto k = lo.z; k <= hi.z; ++k) { for (auto j = lo.y; j <= hi.y; ++j) { for (auto i = lo.x; i <= hi.x; ++i) { - if (mask(i,j,k) >= 2.0) { + if (mask(i,j,k) == 2.0 || mask(i,j,k) == 3.0) { has_SC = 1; } } diff --git a/Source/main.cpp b/Source/main.cpp index 6cbcf9d..e84a4de 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -144,8 +144,8 @@ void main_main (c_FerroX& rFerroX) angle_theta.setVal(0.); //Initialize material mask - InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); - //InitializeMaterialMask(rFerroX, geom, MaterialMask); + //InitializeMaterialMask(MaterialMask, geom, prob_lo, prob_hi); + InitializeMaterialMask(rFerroX, geom, MaterialMask); if(Coordinate_Transformation == 1){ Initialize_tphase_Mask(rFerroX, geom, tphaseMask); Initialize_Euler_angles(rFerroX, geom, angle_alpha, angle_beta, angle_theta); @@ -285,9 +285,9 @@ void main_main (c_FerroX& rFerroX) //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); @@ -302,7 +302,7 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho(PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -405,9 +405,9 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_new_pre, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); @@ -423,7 +423,7 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho(PoissonPhi, P_new_pre, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -475,9 +475,9 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_new, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); @@ -493,7 +493,7 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho(PoissonPhi, P_new, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates @@ -634,9 +634,9 @@ void main_main (c_FerroX& rFerroX) // Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); - dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + //dF_dPhi(alpha_cc, PoissonRHS, PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + //ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); @@ -652,7 +652,7 @@ void main_main (c_FerroX& rFerroX) PoissonPhi.FillBoundary(geom.periodicity()); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + ComputeRho(PoissonPhi, P_old, charge_den, e_den, hole_den, MaterialMask, geom, prob_lo, prob_hi); if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates