diff --git a/Exec/inputs_newton_debug_mfis b/Exec/inputs_newton_debug_mfis new file mode 100644 index 0000000..03f503c --- /dev/null +++ b/Exec/inputs_newton_debug_mfis @@ -0,0 +1,121 @@ +################################# +###### PROBLEM DOMAIN ###### +################################# + +domain.prob_lo = -8.064e-9 -8.064e-9 0.e-9 +domain.prob_hi = 8.064e-9 8.064e-9 5.8e-9 + +domain.n_cell = 64 64 32 #dx = dy = 5.04A, size of HZO unit cell + +domain.max_grid_size = 64 64 32 +domain.blocking_factor = 64 64 32 + +domain.coord_sys = cartesian + +prob_type = 2 + +TimeIntegratorOrder = 1 + +nsteps = -1 #200 #200000 +plot_int = 500 + +dt = 0.5e-13 + +############################################ +###### COORDINATE TRANSFORMATION ########### +############################################ + +#Coordinate_Transformation = 0 +# +#tphase_geom.tphase_geom_function(x,y,z) = "1. - 1.*(x > -7.5e-9)*(x < 7.5e-9) *(y > -7.5e-9)*(y < 7.5e-9)" +# +##Set different Euler angles in each o-phase in degrees. Note that o-phase is located at (x,y,z) where tphase_geom_function(x,y,z) = 0.0 +##Copy the tphase_geom_function(x,y,z) above and specify the values where the coefficients of the logicals are zero to avoid mistakes. +##It doesn't matter what you set the angles in the t_phase. +#angle_alpha.alpha_function(x,y,z) = "0." +# +#angle_beta.beta_function(x,y,z) = "0." +# +#angle_theta.theta_function(x,y,z) = "0." +# +# +#epsilonX_fe_tphase = 40.0 +# +############################################ +###### POLARIZATION BOUNDARY CONDITIONS #### +############################################ + +P_BC_flag_lo = 3 3 0 +P_BC_flag_hi = 3 3 1 +lambda = 3.0e-10 + +############################################ +###### ELECTRICAL BOUNDARY CONDITIONS ###### +############################################ + +domain.is_periodic = 1 1 0 + +boundary.hi = per per dir(Zmax) +boundary.lo = per per dir(Zmin) + +boundary.Zmax_function = "0.9" +boundary.Zmin_function = "0.0" + +voltage_sweep = 1 +Phi_Bc_lo = 0.0 +Phi_Bc_hi = 0.9 + +Phi_Bc_inc = 0.05 +Phi_Bc_hi_max = 1. +phi_tolerance = 5.e-5 +num_Vapp_max = 20 +#mlmg_verbosity = 2 + +################################# +###### STACK GEOMETRY ########### +################################# + +SC_lo = -8.064e-9 -8.064e-9 0.0e-9 +SC_hi = 8.064e-9 8.064e-9 5.0e-9 + +DE_lo = -8.064e-9 -8.064e-9 5.0e-9 +DE_hi = 8.064e-9 8.064e-9 5.8e-9 + +FE_lo = -1. -1. -1. +FE_hi = -1. -1. -1. + +Channel_lo = -1. -1. -1. +Channel_hi = -1. -1. -1. +#FE_lo = -8.e-9 -8.e-9 4.0e-9 +#FE_hi = 8.e-9 8.e-9 9.0e-9 +##FE:0.0, DE:1.0, Source/Drain:2.0, Channel:3.0 + +#device_geom.device_geom_function(x,y,z) = "0.*(x > -8.e-9)*(x < 8.e-9) * (y > -8.e-9)*(y < 8.e-9) * (z > 4.0e-9)*(z < 9.e-9) +# + 1.*(x > -8.e-9)*(x < 8.e-9) * (y > -8.e-9)*(y < 8.e-9) * (z > 0.0)*(z < 4.e-9)" + +################################# +###### MATERIAL PROPERTIES ###### +################################# + +epsilon_0 = 8.85e-12 +epsilonX_fe = 25.0 +epsilonZ_fe = 25.0 +epsilon_de = 3.9 +epsilon_si = 11.7 +alpha = -2.5e9 +beta = 6.0e10 +gamma = 1.5e11 +BigGamma = 100 +g44 = 1.0e-10 #gy, gz +g11 = 1.0e-10 #-1.0e-12 #gx +#g11 = 1.0e-10 +g44_p = 0.0 +g12 = 0.0 +alpha_12 = 0.0 +alpha_112 = 0.0 +alpha_123 = 0.0 + +acceptor_doping = 1.e21 #9.696e+15 +acceptor_doping = 9.696e+15 +donor_doping = 0. + diff --git a/Source/FerroX.cpp b/Source/FerroX.cpp index 861ba14..886e47d 100644 --- a/Source/FerroX.cpp +++ b/Source/FerroX.cpp @@ -206,6 +206,8 @@ AMREX_GPU_MANAGED amrex::Real FerroX::kb; AMREX_GPU_MANAGED amrex::Real FerroX::T; AMREX_GPU_MANAGED amrex::Real FerroX::acceptor_doping; AMREX_GPU_MANAGED amrex::Real FerroX::donor_doping; +AMREX_GPU_MANAGED amrex::Real FerroX::intrinsic_carrier_concentration; +AMREX_GPU_MANAGED int FerroX::use_Fermi_Dirac; // P and Phi Bc AMREX_GPU_MANAGED amrex::Real FerroX::lambda; @@ -235,6 +237,7 @@ AMREX_GPU_MANAGED amrex::Real FerroX::Phi_Bc_inc; AMREX_GPU_MANAGED amrex::Real FerroX::Phi_Bc_hi_max; AMREX_GPU_MANAGED amrex::Real FerroX::phi_tolerance; AMREX_GPU_MANAGED int FerroX::random_seed; +AMREX_GPU_MANAGED int FerroX::num_Vapp_max; //Maximum number of applied voltage points to sweep void InitializeFerroXNamespace(const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi) { @@ -343,6 +346,9 @@ void InitializeFerroXNamespace(const amrex::GpuArray& hole_den_arr = p_den.array(mfi); const Array4& e_den_arr = e_den.array(mfi); @@ -30,36 +31,41 @@ void ComputeRho(MultiFab& PoissonPhi, if (mask(i,j,k) >= 2.0) { - //Maxwell-Boltzmann -// hole_den_arr(i,j,k) = Nv*exp(-(q*phi(i,j,k) - Ev*1.602e-19)/(kb*T)); -// e_den_arr(i,j,k) = Nc*exp(-(Ec*1.602e-19 - q*phi(i,j,k))/(kb*T)); -// charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k)); - - //Fermi-Dirac - Real eta_n = q*(phi(i,j,k) - Ec)/(kb*T); - Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2))); - Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8)); - Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0); - - e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; - - Real eta_p = q*(Ev - phi(i,j,k))/(kb*T); - Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2))); - Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8)); - Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0); - - hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; - - //If in channel, set acceptor doping, else (Source/Drain) set donor doping - if (mask(i,j,k) == 3.0) { - acceptor_den_arr(i,j,k) = acceptor_doping; - donor_den_arr(i,j,k) = 0.0; - } else { // Source / Drain - acceptor_den_arr(i,j,k) = 0.0; - donor_den_arr(i,j,k) = donor_doping; - } - - charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); + if(use_Fermi_Dirac == 1){ + //Fermi-Dirac + Real eta_n = q*(phi(i,j,k) - Ec)/(kb*T); + Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2))); + Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8)); + Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0); + + e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; + + Real eta_p = q*(Ev - phi(i,j,k))/(kb*T); + Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2))); + Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8)); + Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0); + + hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; + } else { + //Maxwell-Boltzmann + Real p_0 = acceptor_doping; + Real n_0 = intrinsic_carrier_concentration*intrinsic_carrier_concentration/p_0; + //Real n_0 = intrinsic_carrier_concentration; + //Real p_0 = intrinsic_carrier_concentration; + hole_den_arr(i,j,k) = p_0*exp(-(q*phi(i,j,k))/(kb*T)); + e_den_arr(i,j,k) = n_0*exp(q*phi(i,j,k)/(kb*T)); + } + + //If in channel, set acceptor doping, else (Source/Drain) set donor doping + if (mask(i,j,k) == 3.0) { + acceptor_den_arr(i,j,k) = acceptor_doping; + donor_den_arr(i,j,k) = 0.0; + } else { // Source / Drain + acceptor_den_arr(i,j,k) = 0.0; + donor_den_arr(i,j,k) = donor_doping; + } + + charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k) - acceptor_den_arr(i,j,k) + donor_den_arr(i,j,k)); } else { @@ -68,5 +74,5 @@ void ComputeRho(MultiFab& PoissonPhi, } }); } - } - + rho.FillBoundary(geom.periodicity()); + } diff --git a/Source/Solver/DerivativeAlgorithm.H b/Source/Solver/DerivativeAlgorithm.H index 8f3256d..0176ca8 100644 --- a/Source/Solver/DerivativeAlgorithm.H +++ b/Source/Solver/DerivativeAlgorithm.H @@ -12,13 +12,18 @@ using namespace FerroX; amrex::GpuArrayconst& prob_lo, amrex::GpuArrayconst& prob_hi) { - if (z_lo < prob_lo[2]){ // bottom metal - return (-4.*F(i,j,k-1) + 3.*F(i,j,k) + F(i,j,k+1)) / (3. * dx[2]); - } else if (z_hi > prob_hi[2]){ // top metal - return (4.*F(i,j,k+1) - 3.*F(i,j,k) - F(i,j,k-1)) / (3. * dx[2]); - } else { // inside stack - return (F(i,j,k+1) - F(i,j,k-1)) / (2. * dx[2]); - } + //if (z_lo < prob_lo[2]){ // bottom metal + // return (-4.*F(i,j,k-1) + 3.*F(i,j,k) + F(i,j,k+1)) / (3. * dx[2]); + //} else if (z_hi > prob_hi[2]){ // top metal + // return (4.*F(i,j,k+1) - 3.*F(i,j,k) - F(i,j,k-1)) / (3. * dx[2]); + //} else { // inside stack + // return (F(i,j,k+1) - F(i,j,k-1)) / (2. * dx[2]); + //} + return ( F(i,j, k+1) - F(i,j,k) + + F(i+1,j, k+1) - F(i+1,j,k) + + F(i,j+1, k+1) - F(i,j+1,k) + + F(i+1,j+1,k+1) - F(i+1,j+1,k) + )/4./dx[2]; } /** @@ -27,7 +32,11 @@ using namespace FerroX; static amrex::Real DFDx ( amrex::Array4 const& F, int const i, int const j, int const k, amrex::GpuArray dx) { - return (F(i+1,j,k) - F(i-1,j,k))/(2.*dx[0]); + return ( F(i+1,j, k) - F(i,j,k) + + F(i+1,j+1, k) - F(i,j+1,k) + + F(i+1,j, k+1) - F(i,j,k+1) + + F(i+1,j+1,k+1) - F(i,j+1,k+1) + )/4./dx[0]; } /** @@ -36,7 +45,82 @@ using namespace FerroX; static amrex::Real DFDy ( amrex::Array4 const& F, int const i, int const j, int const k, amrex::GpuArray dx) { - return (F(i,j+1,k) - F(i,j-1,k))/(2.*dx[1]); + return ( F(i, j+1, k) - F(i,j,k) + + F(i+1,j+1, k) - F(i+1,j,k) + + F(i ,j+1,k+1) - F(i,j,k+1) + + F(i+1,j+1,k+1) - F(i+1,j,k+1) + )/4./dx[1]; + } + +/** + * Perform first derivative dP/dx */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real NodalDPDx ( + amrex::Array4 const& F, + amrex::Array4 const& mask, + int const i, int const j, int const k, amrex::GpuArray dx +) { + // if (mask(i-1,j,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary + // + // if(P_BC_flag_lo[0] == 0){ + // Real F_lo = 0.0; + // return (-4.*F_lo + 3.*F(i,j,k) + F(i+1,j,k))/(3.*dx[0]);//2nd order using three point stencil using 0, pOld(i,j,k), and pOld(i+1,j,k) + + // } else if (P_BC_flag_lo[0] == 1){ + + // Real F_lo = F(i,j,k)/(1 + dx[0]/2/lambda); + // return (dx[0]*F_lo/lambda - F(i,j,k) + F(i+1,j,k))/(2.*dx[0]); // 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; + + // } else if (P_BC_flag_lo[0] == 2){ + // return ( - F(i,j,k) + F(i+1,j,k))/(2.*dx[0]); //dPdx = 0. + + // } else if (P_BC_flag_lo[0] == 3){ + // return (- F(i-1,j,k) + F(i+1,j,k))/(2.*dx[0]); //periodic + + // } else { + // amrex::Abort("Wrong flag of the lower x polarization boundary condition!!"); + // return 0.0; + // } + + // } else if (mask(i+1,j,k) != 0.0 && mask(i,j,k) == 0.0){ // FE higher boundary + + // if(P_BC_flag_hi[0] == 0){ + // Real F_hi = 0.0; + // return (4.*F_hi - 3.*F(i,j,k) - F(i-1,j,k))/(3.*dx[0]);//2nd order using three point stencil using 0, pOld(i,j,k), and pOld(i-1,j,k) + + // } else if (P_BC_flag_hi[0] == 1){ + // + // Real F_hi = F(i,j,k)/(1 - dx[0]/2/lambda); + // return (dx[0]*F_hi/lambda + F(i,j,k) - F(i-1,j,k))/(2.*dx[0]);//dPdz = P_hi/lambda; + + // // Real F_hi = (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_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; + + // } else if (P_BC_flag_hi[0] == 2){ + // return (F(i,j,k) - F(i-1,j,k))/(2.*dx[0]); //dPdx = 0. + + // } else if (P_BC_flag_hi[0] == 3){ + // return (F(i+1,j,k) - F(i-1,j,k))/(2.*dx[0]); //Periodic + + // } else { + // amrex::Abort("Wrong flag of the higher x polarization boundary condition!!"); + // return 0.0; + // } + // + // } else if (mask(i,j,k) == 0.0) { // inside FE + // return (F(i+1,j,k) - F(i-1,j,k))/(2.*dx[0]); + + // } else { + // return 0.0; + // } + return ( F(i,j-1,k-1) - F(i-1,j-1,k-1) + + F(i,j, k-1) - F(i-1,j ,k-1) + + F(i,j-1, k) - F(i-1,j-1,k) + + F(i,j, k) - F(i-1,j, k) + )/4./dx[0]; } /** @@ -105,6 +189,78 @@ using namespace FerroX; } } +/** + * Perform first derivative dP/dy */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real NodalDPDy ( + amrex::Array4 const& F, + amrex::Array4 const& mask, + int const i, int const j, int const k, amrex::GpuArray dx +) { + +// if (mask(i,j-1,k) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary +// +// if(P_BC_flag_lo[1] == 0){ +// Real F_lo = 0.0; +// return (-4.*F_lo + 3.*F(i,j,k) + F(i,j+1,k))/(3.*dx[1]);//2nd order using three point stencil using 0, pOld(i,j,k), and pOld(i,j,k+1) +// +// } else if (P_BC_flag_lo[1] == 1){ +// +// Real F_lo = F(i,j,k)/(1 + dx[1]/2/lambda); +// return (dx[1]*F_lo/lambda - F(i,j,k) + F(i,j+1,k))/(2.*dx[1]); // 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; +// +// } else if (P_BC_flag_lo[1] == 2){ +// return ( - F(i,j,k) + F(i,j+1,k))/(2.*dx[1]); //dPdy = 0. +// +// } else if (P_BC_flag_lo[1] == 3){ +// return ( - F(i,j-1,k) + F(i,j+1,k))/(2.*dx[1]); //Periodic +// +// } else { +// amrex::Abort("Wrong flag of the lower polarization boundary condition!!"); +// return 0.0; +// } +// +// } else if (mask(i,j+1,k) != 0.0 && mask(i,j,k) == 0.0){ // FE higher boundary +// +// if(P_BC_flag_hi[1] == 0){ +// Real F_hi = 0.0; +// return (4.*F_hi - 3.*F(i,j,k) - F(i,j-1,k))/(3.*dx[1]);//2nd order using three point stencil using 0, pOld(i,j,k), and pOld(i,j,k-1) +// +// } else if (P_BC_flag_hi[1] == 1){ +// +// Real F_hi = F(i,j,k)/(1 - dx[1]/2/lambda); +// return (dx[1]*F_hi/lambda + F(i,j,k) - F(i,j-1,k))/(2.*dx[1]);//dPdz = P_hi/lambda; +// +// // Real F_hi = (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_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; +// +// } else if (P_BC_flag_hi[1] == 2){ +// return (F(i,j,k) - F(i,j-1,k))/(2.*dx[1]); //dPdy = 0. +// +// } else if (P_BC_flag_hi[1] == 3){ +// return (F(i,j+1,k) - F(i,j-1,k))/(2.*dx[1]); //Periodic +// +// } else { +// amrex::Abort("Wrong flag of the higher y polarization boundary condition!!"); +// return 0.0; +// } +// +// } else if (mask(i,j,k) == 0.0) { // inside FE +// return (F(i,j+1,k) - F(i,j-1,k))/(2.*dx[1]); +// +// } else { +// return 0.0; +// } + return ( F(i-1,j,k-1) - F(i-1,j-1,k-1) + + F(i-1,j, k) - F(i-1,j-1,k) + + F(i, j,k-1) - F(i,j-1,k-1) + + F(i, j, k) - F(i,j-1, k) + )/4./dx[1]; + } + /** * Perform first derivative dP/dy */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -172,6 +328,77 @@ using namespace FerroX; } } +/** + * Perform first derivative dP/dz */ + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + static amrex::Real NodalDPDz ( + 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) != 0.0 && mask(i,j,k) == 0.0) { //FE lower boundary +// +// if(P_BC_flag_lo[2] == 0){ +// Real F_lo = 0.0; +// return (-4.*F_lo + 3.*F(i,j,k) + F(i,j,k+1))/(3.*dx[2]);//2nd order using three point stencil using 0, pOld(i,j,k), and pOld(i,j,k+1) +// +// } 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 = (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; +// +// } else if (P_BC_flag_lo[2] == 2){ +// return ( - F(i,j,k) + F(i,j,k+1))/(2.*dx[2]); //dPdz = 0. +// +// } else if (P_BC_flag_lo[2] == 3){ +// return ( - F(i,j,k-1) + F(i,j,k+1))/(2.*dx[2]); //No BC +// +// } else { +// amrex::Abort("Wrong flag of the lower polarization boundary condition!!"); +// return 0.0; +// } +// +// } else if ( mask(i,j,k+1) != 0.0 && mask(i,j,k) == 0.0 ){ // FE higher boundary +// +// if(P_BC_flag_hi[2] == 0){ +// Real F_hi = 0.0; +// return (4.*F_hi - 3.*F(i,j,k) - F(i,j,k-1))/(3.*dx[2]);//2nd order using three point stencil using 0, pOld(i,j,k), and pOld(i,j,k-1) +// +// } else if (P_BC_flag_hi[2] == 1){ +// +// Real F_hi = F(i,j,k)/(1 - dx[2]/2/lambda); +// return (dx[2]*F_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; +// +// // Real F_hi = (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_hi/lambda + F(i,j,k) - F(i,j,k-1))/(2.*dx[2]);//dPdz = P_hi/lambda; +// +// } else if (P_BC_flag_hi[2] == 2){ +// return (F(i,j,k) - F(i,j,k-1))/(2.*dx[2]); //dPdz = 0. +// +// } else if (P_BC_flag_hi[2] == 3){ +// return (0. - F(i,j,k-1))/(2.*dx[2]); //No BC +// +// } else { +// amrex::Abort("Wrong flag of the higher polarization boundary condition!!"); +// return 0.0; +// } +// +// } else if (mask(i,j,k) == 0.0) { // inside FE +// return (F(i,j,k+1) - F(i,j,k-1))/(2.*dx[2]); +// +// } else { +// return 0.0; +// } + return ( F(i-1,j-1,k) - F(i-1,j-1,k-1) + + F(i-1,j, k) - F(i-1,j, k-1) + + F(i, j-1,k) - F(i, j-1,k-1) + + F(i, j, k) - F(i, j, k-1) + )/4./dx[2]; + } /** * Perform first derivative dP/dz */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Solver/ElectrostaticSolver.H b/Source/Solver/ElectrostaticSolver.H index 592c99e..eb35447 100644 --- a/Source/Solver/ElectrostaticSolver.H +++ b/Source/Solver/ElectrostaticSolver.H @@ -53,10 +53,17 @@ void dF_dPhi(MultiFab& alpha_cc, void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& PoissonPhi, + const Geometry& geom, MultiFab& alpha_cc); +void ComputeAx_H(MultiFab& x_H, MultiFab& Ax_H, const amrex::GpuArray& n_cell); + void SetPhiBC_z(MultiFab& PossonPhi, const amrex::GpuArray& n_cell); +void SetPhiBC_z_after_solve(MultiFab& PossonPhi, const amrex::GpuArray& n_cell); + +void average_cc_to_nodes(MultiFab& mf_nodal, const MultiFab& mf_cc, 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); void Fill_Constant_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& PoissonPhi); diff --git a/Source/Solver/ElectrostaticSolver.cpp b/Source/Solver/ElectrostaticSolver.cpp index 486ca6d..2f40bad 100644 --- a/Source/Solver/ElectrostaticSolver.cpp +++ b/Source/Solver/ElectrostaticSolver.cpp @@ -32,8 +32,12 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, //Convert Euler angles from degrees to radians amrex::Real Pi = 3.14159265358979323846; - amrex::Real alpha_rad = Pi/180.*angle_alpha_arr(i,j,k); - amrex::Real beta_rad = Pi/180.*angle_beta_arr(i,j,k); + + // amrex::Real alpha_rad = 0.*Pi/180.;//Pi/180.*angle_alpha_arr(i,j,k); + // amrex::Real beta_rad = 45.*Pi/180.;//Pi/180.*angle_beta_arr(i,j,k); + // amrex::Real theta_rad = 0.*Pi/180.;//Pi/180.*angle_theta_arr(i,j,k); + amrex::Real alpha_rad = Pi/180.*angle_alpha_arr(i,j,k); + amrex::Real beta_rad = Pi/180.*angle_beta_arr(i,j,k); amrex::Real theta_rad = Pi/180.*angle_theta_arr(i,j,k); amrex::Real R_11, R_12, R_13, R_21, R_22, R_23, R_31, R_32, R_33; @@ -60,24 +64,43 @@ void ComputePoissonRHS(MultiFab& PoissonRHS, R_33 = cos(alpha_rad)*cos(beta_rad); } + //when coordinate transformation is OFF, + //alpha = beta = theta = 0. + //Therefore, R_11 = R_22 = R_33 = 1, R_12 = R_13 = R_21 = R_23 = R_31 = R_32 = 0. + + if (Coordinate_Transformation != 1){ + if (R_11 != 1.0 || R_12 != 0.0 || R_13 != 0.0 || + R_21 != 0.0 || R_22 != 1.0 || R_23 != 0.0 || + R_31 != 0.0 || R_32 != 0.0 || R_33 != 1.0 ){ + amrex::Print() << "alpha = " << alpha_rad << ", beta = " << beta_rad << ", theta = " << theta_rad << "\n"; + amrex::Abort("ComputePoissonRHS : Coordinate transformation is turned OFF, but rotation matrix is not an identity matrix!"); + } + } + if(mask(i,j,k) >= 2.0){ //SC region RHS(i,j,k) = charge_den_arr(i,j,k); - +// RHS(i,j,k) *= -1.; + //amrex::Print() << "RHS(i,j,k) = " << RHS(i,j,k) << "\n"; + //amrex::Print() << "charge_den_arr(i,j,k) = " << charge_den_arr(i,j,k) << "\n"; } else if(mask(i,j,k) == 1.0){ //DE region RHS(i,j,k) = 0.; } else { //mask(i,j,k) == 0.0 FE region + //RHS(i,j,k) = - (R_11*NodalDPDx(pOld_p, mask, i, j, k, dx) + R_12*NodalDPDy(pOld_p, mask, i, j, k, dx) + R_13*NodalDPDz(pOld_p, mask, i, j, k, dx)) + // - (R_21*NodalDPDx(pOld_q, mask, i, j, k, dx) + R_22*NodalDPDy(pOld_q, mask, i, j, k, dx) + R_23*NodalDPDz(pOld_q, mask, i, j, k, dx)) + // - (R_31*NodalDPDx(pOld_r, mask, i, j, k, dx) + R_32*NodalDPDy(pOld_r, mask, i, j, k, dx) + R_33*NodalDPDz(pOld_r, mask, i, j, k, dx)); RHS(i,j,k) = - (R_11*DPDx(pOld_p, mask, i, j, k, dx) + R_12*DPDy(pOld_p, mask, i, j, k, dx) + R_13*DPDz(pOld_p, mask, i, j, k, dx)) - (R_21*DPDx(pOld_q, mask, i, j, k, dx) + R_22*DPDy(pOld_q, mask, i, j, k, dx) + R_23*DPDz(pOld_q, mask, i, j, k, dx)) - (R_31*DPDx(pOld_r, mask, i, j, k, dx) + R_32*DPDy(pOld_r, mask, i, j, k, dx) + R_33*DPDz(pOld_r, mask, i, j, k, dx)); + /// RHS(i,j,k) *= -1.; } }); } - + PoissonRHS.FillBoundary(geom.periodicity()); } void dF_dPhi(MultiFab& alpha_cc, @@ -95,28 +118,47 @@ void dF_dPhi(MultiFab& alpha_cc, { - MultiFab PoissonPhi_plus_delta(PoissonPhi.boxArray(), PoissonPhi.DistributionMap(), 1, 0); - MultiFab PoissonRHS_phi_plus_delta(PoissonRHS.boxArray(), PoissonRHS.DistributionMap(), 1, 0); + MultiFab PoissonPhi_plus_delta(PoissonPhi.boxArray(), PoissonPhi.DistributionMap(), 1, 1); + MultiFab PoissonRHS_phi_plus_delta(PoissonRHS.boxArray(), PoissonRHS.DistributionMap(), 1, 1); - MultiFab::Copy(PoissonPhi_plus_delta, PoissonPhi, 0, 0, 1, 0); - PoissonPhi_plus_delta.plus(delta, 0, 1, 0); + MultiFab::Copy(PoissonPhi_plus_delta, PoissonPhi, 0, 0, 1, 1); + PoissonPhi_plus_delta.plus(delta, 0, 1, 1); // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, rho, e_den, p_den, MaterialMask); + ComputeRho(PoissonPhi_plus_delta, rho, e_den, p_den, geom, MaterialMask); //Compute RHS of Poisson equation ComputePoissonRHS(PoissonRHS_phi_plus_delta, P_old, rho, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); MultiFab::LinComb(alpha_cc, 1./delta, PoissonRHS_phi_plus_delta, 0, -1./delta, PoissonRHS, 0, 0, 1, 0); +// +// for ( MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi ) +// { +// const Box& bx = mfi.validbox(); +// //const Box& bx = mfi.growntilebox(1); +// +// const Array4& alpha = alpha_cc.array(mfi); +// const Array4& phi_p_delta = PoissonPhi_plus_delta.array(mfi); +// const Array4& rhs_phi_p_delta = PoissonRHS_phi_plus_delta.array(mfi); +// const Array4& rhs = PoissonRHS.array(mfi); +// +// amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) +// { +// if(std::isnan(alpha(i,j,k))) amrex::Print() <<"alpha(" << i << ", " << j << ", " << k << ") = " << alpha(i,j,k) << ", phi_p_delta = " << phi_p_delta(i,j,k) << ", rhs_phi_p_delta = " << rhs_phi_p_delta(i,j,k) << ", rhs = " << rhs(i,j,k) << "\n"; +// }); +// } + alpha_cc.FillBoundary(geom.periodicity()); } void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, MultiFab& PoissonPhi, + const Geometry& geom, MultiFab& alpha_cc) { for ( MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi ) { const Box& bx = mfi.validbox(); + //const Box& bx = mfi.growntilebox(1); const Array4& phi = PoissonPhi.array(mfi); const Array4& poissonRHS = PoissonRHS.array(mfi); @@ -124,9 +166,13 @@ void ComputePoissonRHS_Newton(MultiFab& PoissonRHS, amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + //if(std::isnan(alpha(i,j,k))) amrex::Print() << "poissonRHS(i,j,k) = " << poissonRHS(i,j,k) << ", alpha(i,j,k) = " << alpha(i,j,k) << ", phi(i,j,k) = " << phi(i,j,k) << ", phi_prev2(i,j,k) = " << phi_prev2(i,j,k) << "\n"; + //amrex::Print() << "poissonRHS(i,j,k) = " << poissonRHS(i,j,k) << ", alpha(i,j,k) = " << alpha(i,j,k) << ", phi(i,j,k) = " << phi(i,j,k) << ", phi_prev2(i,j,k) = " << phi_prev2(i,j,k) << "\n"; poissonRHS(i,j,k) = poissonRHS(i,j,k) - alpha(i,j,k)*phi(i,j,k) ; }); } + + PoissonRHS.FillBoundary(geom.periodicity()); } void ComputeEfromPhi(MultiFab& PoissonPhi, @@ -138,7 +184,7 @@ void ComputeEfromPhi(MultiFab& PoissonPhi, { // Calculate E from Phi - for ( MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi ) + for ( MFIter mfi(E[0]); mfi.isValid(); ++mfi ) { const Box& bx = mfi.validbox(); @@ -162,8 +208,11 @@ void ComputeEfromPhi(MultiFab& PoissonPhi, //Convert Euler angles from degrees to radians amrex::Real Pi = 3.14159265358979323846; +// amrex::Real alpha_rad = 0.*Pi/180.;//Pi/180.*angle_alpha_arr(i,j,k); +// amrex::Real beta_rad = 45.*Pi/180.;//Pi/180.*angle_beta_arr(i,j,k); +// amrex::Real theta_rad = 0.*Pi/180.;//Pi/180.*angle_theta_arr(i,j,k); amrex::Real alpha_rad = Pi/180.*angle_alpha_arr(i,j,k); - amrex::Real beta_rad = Pi/180.*angle_beta_arr(i,j,k); + amrex::Real beta_rad = Pi/180.*angle_beta_arr(i,j,k); amrex::Real theta_rad = Pi/180.*angle_theta_arr(i,j,k); amrex::Real R_11, R_12, R_13, R_21, R_22, R_23, R_31, R_32, R_33; @@ -190,6 +239,19 @@ void ComputeEfromPhi(MultiFab& PoissonPhi, R_33 = cos(alpha_rad)*cos(beta_rad); } + //when coordinate transformation is OFF, + //alpha = beta = theta = 0. + //Therefore, R_11 = R_22 = R_33 = 1, R_12 = R_13 = R_21 = R_23 = R_31 = R_32 = 0. + //So, Ep = Ex = -DFDx(phi), Eq = Ey = -DFDy(phi), Er = Ez = -DphiDz(phi) + + if (Coordinate_Transformation != 1){ + if (R_11 != 1.0 || R_12 != 0.0 || R_13 != 0.0 || + R_21 != 0.0 || R_22 != 1.0 || R_23 != 0.0 || + R_31 != 0.0 || R_32 != 0.0 || R_33 != 1.0 ){ + amrex::Abort("Coordinate transformation is turned OFF, but rotation matrix is not an identity matrix!"); + } + } + Ep_arr(i,j,k) = - (R_11*DFDx(phi, i, j, k, dx) + R_12*DFDy(phi, i, j, k, dx) + R_13*DphiDz(phi, z_hi, z_lo, i, j, k, dx, prob_lo, prob_hi)); Eq_arr(i,j,k) = - (R_21*DFDx(phi, i, j, k, dx) + R_22*DFDy(phi, i, j, k, dx) + R_23*DphiDz(phi, z_hi, z_lo, i, j, k, dx, prob_lo, prob_hi)); Er_arr(i,j,k) = - (R_31*DFDx(phi, i, j, k, dx) + R_32*DFDy(phi, i, j, k, dx) + R_33*DphiDz(phi, z_hi, z_lo, i, j, k, dx, prob_lo, prob_hi)); @@ -198,6 +260,11 @@ void ComputeEfromPhi(MultiFab& PoissonPhi, }); } + for (int i = 0; i < 3; i++){ + // fill periodic ghost cells + E[i].FillBoundary(geom.periodicity()); + } + } void InitializePermittivity(std::array,2>& LinOpBCType_2d, @@ -497,22 +564,83 @@ void Fill_FunctionBased_Inhomogeneous_Boundaries(c_FerroX& rFerroX, MultiFab& Po } } +//A multifab filled with zeros, but boundary cells filled to respect bc's void SetPhiBC_z(MultiFab& PoissonPhi, const amrex::GpuArray& n_cell) { +// PoissonPhi.setVal(0.); for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) { - const Box& bx = mfi.growntilebox(1); + const Box& bx = mfi.validbox(); + + const Array4& Phi = PoissonPhi.array(mfi); + + 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]){ + Phi(i,j,k) = Phi_Bc_hi; + } + }); + } +} + + +//A multifab filled with zeros, but boundary cells filled to respect bc's +void SetPhiBC_z_after_solve(MultiFab& PoissonPhi, const amrex::GpuArray& n_cell) +{ + for (MFIter mfi(PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); const Array4& Phi = PoissonPhi.array(mfi); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { - if(k < 0) { + if(k == 0) { Phi(i,j,k) = Phi_Bc_lo; - } else if(k >= n_cell[2]){ + } else if(k == n_cell[2]){ Phi(i,j,k) = Phi_Bc_hi; } }); } } +//Avergae cell-centered multifab to nodes +void average_cc_to_nodes(MultiFab& mf_nodal, const MultiFab& mf_cc, const Geometry& geom) +{ + for (MFIter mfi(mf_nodal); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + + const Array4& mf_cc_arr = mf_cc.array(mfi); + const Array4& mf_nodal_arr = mf_nodal.array(mfi); + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + mf_nodal_arr(i,j,k) = 1./8.*(mf_cc_arr(i-1, j-1, k-1) + + mf_cc_arr(i, j-1, k-1) + + mf_cc_arr(i-1, j, k-1) + + mf_cc_arr(i, j, k-1) + + mf_cc_arr(i-1, j-1, k ) + + mf_cc_arr(i, j-1, k ) + + mf_cc_arr(i-1, j, k ) + + mf_cc_arr(i, j, k )); + if(k == 0){ + mf_nodal_arr(i,j,k) = 1./4.*(mf_cc_arr(i-1, j-1, k ) + + mf_cc_arr(i, j-1, k ) + + mf_cc_arr(i-1, j, k ) + + mf_cc_arr(i, j, k )); + } + if(k == 32){ + mf_nodal_arr(i,j,k) = 1./4.*(mf_cc_arr(i-1, j-1,k-1 ) + + mf_cc_arr(i, j-1, k-1 ) + + mf_cc_arr(i-1, j, k-1 ) + + mf_cc_arr(i, j, k-1 )); + } + }); + + } + mf_nodal.FillBoundary(geom.periodicity()); +} + diff --git a/Source/Solver/Initialization.H b/Source/Solver/Initialization.H index 8a929d5..59e7ab9 100644 --- a/Source/Solver/Initialization.H +++ b/Source/Solver/Initialization.H @@ -24,7 +24,13 @@ void InitializeMaterialMask(MultiFab& MaterialMask, const amrex::GpuArray& prob_lo, const amrex::GpuArray& prob_hi); +void InitializeMaterialMask_Nodal(MultiFab& MaterialMask, + const Geometry& geom, + const amrex::GpuArray& prob_lo, + const amrex::GpuArray& prob_hi); + void InitializeMaterialMask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& MaterialMask); void Initialize_tphase_Mask(c_FerroX& rFerroX, const Geometry& geom, MultiFab& tphaseMask); void Initialize_Euler_angles(c_FerroX& rFerroX, const Geometry& geom, MultiFab& angle_alpha, MultiFab& angle_beta, MultiFab& angle_theta); +void Initialize_nodal_Euler_angles(c_FerroX& rFerroX, const Geometry& geom, MultiFab& Nodal_angle_alpha, MultiFab& Nodal_angle_beta, MultiFab& Nodal_angle_theta); diff --git a/Source/Solver/Initialization.cpp b/Source/Solver/Initialization.cpp index 2100462..e458a3b 100644 --- a/Source/Solver/Initialization.cpp +++ b/Source/Solver/Initialization.cpp @@ -59,8 +59,8 @@ void InitializePandRho(Array &P_old, rngs[i] = amrex::Random(); // uniform [0,1] option } - // loop over boxes - for (MFIter mfi(rho); mfi.isValid(); ++mfi) + // loop over cc boxes for P + for (MFIter mfi(P_old[0]); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); @@ -109,17 +109,21 @@ void InitializePandRho(Array &P_old, pOld_r(i,j,k) = 0.0; } + pOld_p(i,j,k) = 0.0; + pOld_q(i,j,k) = 0.0; + } else { + pOld_p(i,j,k) = 0.0; + pOld_q(i,j,k) = 0.0; pOld_r(i,j,k) = 0.0; Gam(i,j,k) = 0.0; } - pOld_p(i,j,k) = 0.0; - pOld_q(i,j,k) = 0.0; }); - // Calculate charge density from Phi, Nc, Nv, Ec, and Ev + + // 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); + MultiFab acceptor_den(rho.boxArray(), rho.DistributionMap(), 1, 1); + MultiFab donor_den(rho.boxArray(), rho.DistributionMap(), 1, 1); const Array4& hole_den_arr = p_den.array(mfi); const Array4& e_den_arr = e_den.array(mfi); @@ -134,49 +138,96 @@ void InitializePandRho(Array &P_old, //SC region if (mask(i,j,k) >= 2.0) { - Real Phi = 0.5*(Ec + Ev); //eV -// hole_den_arr(i,j,k) = Nv*exp(-(Phi - Ev)*1.602e-19/(kb*T)); -// e_den_arr(i,j,k) = Nc*exp(-(Ec - Phi)*1.602e-19/(kb*T)); -// charge_den_arr(i,j,k) = q*(hole_den_arr(i,j,k) - e_den_arr(i,j,k)); + if(use_Fermi_Dirac == 1){ + + //Approximate FD integral + Real Phi = 0.5*(Ec + Ev); //eV + Real eta_n = q*(Phi - Ec)/(kb*T); + Real nu_n = std::pow(eta_n, 4.0) + 50.0 + 33.6 * eta_n * (1 - 0.68 * exp(-0.17 * std::pow((eta_n + 1), 2.0))); + Real xi_n = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_n, 3/8)); + Real FD_half_n = std::pow(exp(-eta_n) + xi_n, -1.0); + + e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; + + Real eta_p = q*(Ev - Phi)/(kb*T); + Real nu_p = std::pow(eta_p, 4.0) + 50.0 + 33.6 * eta_p * (1 - 0.68 * exp(-0.17 * std::pow((eta_p + 1), 2.0))); + Real xi_p = 3.0 * sqrt(3.14)/(4.0 * std::pow(nu_p, 3/8)); + Real FD_half_p = std::pow(exp(-eta_p) + xi_p, -1.0); + + hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; + + } else { + + //hole_den_arr(i,j,k) = intrinsic_carrier_concentration; + //e_den_arr(i,j,k) = intrinsic_carrier_concentration; + hole_den_arr(i,j,k) = acceptor_doping; + e_den_arr(i,j,k) = intrinsic_carrier_concentration*intrinsic_carrier_concentration/acceptor_doping; + + } + } - //Approximate FD integral - 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); + //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)); - e_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nc*FD_half_n; + }); + } + for (int i = 0; i < 3; i++){ + // fill periodic ghost cells + P_old[i].FillBoundary(geom.periodicity()); + } - 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); + rho.FillBoundary(geom.periodicity()); + } - hole_den_arr(i,j,k) = 2.0/sqrt(3.14)*Nv*FD_half_p; +// create a mask filled with integers to idetify different material types +void InitializeMaterialMask_Nodal(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.validbox(); + const Box& bx = mfi.growntilebox(MaterialMask.nGrow()); - //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)); + // extract dx from the geometry object + GpuArray dx = geom.CellSizeArray(); - } else { + const Array4& mask = MaterialMask.array(mfi); - charge_den_arr(i,j,k) = 0.0; - } + amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real x = prob_lo[0] + i* dx[0]; + Real y = prob_lo[1] + j* dx[1]; + Real z = prob_lo[2] + k* 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 + } }); } - for (int i = 0; i < 3; i++){ - // fill periodic ghost cells - P_old[i].FillBoundary(geom.periodicity()); - } + MaterialMask.FillBoundary(geom.periodicity()); +} - } // create a mask filled with integers to idetify different material types void InitializeMaterialMask(MultiFab& MaterialMask, @@ -187,7 +238,8 @@ void InitializeMaterialMask(MultiFab& MaterialMask, // loop over boxes for (MFIter mfi(MaterialMask); mfi.isValid(); ++mfi) { - const Box& bx = mfi.validbox(); + //const Box& bx = mfi.validbox(); + const Box& bx = mfi.growntilebox(MaterialMask.nGrow()); // extract dx from the geometry object GpuArray dx = geom.CellSizeArray(); diff --git a/Source/Solver/TotalEnergyDensity.cpp b/Source/Solver/TotalEnergyDensity.cpp index a945f03..aed4696 100644 --- a/Source/Solver/TotalEnergyDensity.cpp +++ b/Source/Solver/TotalEnergyDensity.cpp @@ -73,95 +73,107 @@ void CalculateTDGL_RHS(Array &GL_rhs, R_33 = cos(alpha_rad)*cos(beta_rad); } - Real dFdPp_Landau = alpha*pOld_p(i,j,k) + beta*std::pow(pOld_p(i,j,k),3.) + FerroX::gamma*std::pow(pOld_p(i,j,k),5.) - + 2. * alpha_12 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) - + 2. * alpha_12 * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),2.) - + 4. * alpha_112 * std::pow(pOld_p(i,j,k),3.) * (std::pow(pOld_q(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) - + 2. * alpha_112 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),4.) - + 2. * alpha_112 * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),4.) - + 2. * alpha_123 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); - - Real dFdPq_Landau = alpha*pOld_q(i,j,k) + beta*std::pow(pOld_q(i,j,k),3.) + FerroX::gamma*std::pow(pOld_q(i,j,k),5.) - + 2. * alpha_12 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) - + 2. * alpha_12 * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),2.) - + 4. * alpha_112 * std::pow(pOld_q(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) - + 2. * alpha_112 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),4.) - + 2. * alpha_112 * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),4.) - + 2. * alpha_123 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); - - Real dFdPr_Landau = alpha*pOld_r(i,j,k) + beta*std::pow(pOld_r(i,j,k),3.) + FerroX::gamma*std::pow(pOld_r(i,j,k),5.) - + 2. * alpha_12 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) - + 2. * alpha_12 * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),2.) - + 4. * alpha_112 * std::pow(pOld_r(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_q(i,j,k),2.)) - + 2. * alpha_112 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),4.) - + 2. * alpha_112 * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),4.) - + 2. * alpha_123 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_q(i,j,k),2.); - - Real dFdPp_grad = - g11 * DoubleDPDx(pOld_p, mask, i, j, k, dx) - - (g44 + g44_p) * DoubleDPDy(pOld_p, mask, i, j, k, dx) - - (g44 + g44_p) * DoubleDPDz(pOld_p, mask, i, j, k, dx) - - (g12 + g44 - g44_p) * DoubleDPDxDy(pOld_q, mask, i, j, k, dx) // d2P/dxdy - - (g12 + g44 - g44_p) * DoubleDPDxDz(pOld_r, mask, i, j, k, dx); // d2P/dxdz - - Real dFdPq_grad = - g11 * DoubleDPDy(pOld_q, mask, i, j, k, dx) - - (g44 - g44_p) * DoubleDPDx(pOld_q, mask, i, j, k, dx) - - (g44 - g44_p) * DoubleDPDz(pOld_q, mask, i, j, k, dx) - - (g12 + g44 + g44_p) * DoubleDPDxDy(pOld_p, mask, i, j, k, dx) // d2P/dxdy - - (g12 + g44 - g44_p) * DoubleDPDyDz(pOld_r, mask, i, j, k, dx);// d2P/dydz - - //Switch g11 and g44 temporarily for multiphase simulations. This will be generalized later - Real dFdPr_grad = - g44 * ( R_31*R_31*DoubleDPDx(pOld_r, mask, i, j, k, dx) - +R_32*R_32*DoubleDPDy(pOld_r, mask, i, j, k, dx) - +R_33*R_33*DoubleDPDz(pOld_r, mask, i, j, k, dx) - +2.*R_31*R_32*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) - +2.*R_32*R_33*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) - +2.*R_33*R_31*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) - - - (g11 - g44_p) * ( R_11*R_11*DoubleDPDx(pOld_r, mask, i, j, k, dx) - +R_12*R_12*DoubleDPDy(pOld_r, mask, i, j, k, dx) - +R_13*R_13*DoubleDPDz(pOld_r, mask, i, j, k, dx) - +2.*R_11*R_12*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) - +2.*R_12*R_13*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) - +2.*R_13*R_11*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) - - - (g44 - g44_p) * ( R_21*R_21*DoubleDPDx(pOld_r, mask, i, j, k, dx) - +R_22*R_22*DoubleDPDy(pOld_r, mask, i, j, k, dx) - +R_23*R_23*DoubleDPDz(pOld_r, mask, i, j, k, dx) - +2.*R_21*R_22*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) - +2.*R_22*R_23*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) - +2.*R_23*R_21*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) - - - (g44 + g44_p + g12) * DoubleDPDyDz(pOld_q, mask, i, j, k, dx) // d2P/dydz - - (g44 + g44_p + g12) * DoubleDPDxDz(pOld_p, mask, i, j, k, dx); // d2P/dxdz - - GL_RHS_p(i,j,k) = -1.0 * Gam(i,j,k) * - ( dFdPp_Landau - + dFdPp_grad - - Ep(i,j,k) - //+ DFDx(phi, i, j, k, dx) - ); - - GL_RHS_q(i,j,k) = -1.0 * Gam(i,j,k) * - ( dFdPq_Landau - + dFdPq_grad - - Eq(i,j,k) - //+ DFDy(phi, i, j, k, dx) - ); - - GL_RHS_p(i,j,k) = 0.0; - GL_RHS_q(i,j,k) = 0.0; - GL_RHS_r(i,j,k) = -1.0 * Gam(i,j,k) * - ( dFdPr_Landau - + dFdPr_grad - - Er(i,j,k) - //+ DphiDz(phi, z_hi, z_lo, i, j, k, dx, prob_lo, prob_hi) - ); - - //set t_phase GL_RHS_r to zero so that it stays zero. It is initialized to zero in t-phase as well - //if(x <= t_phase_hi[0] && x >= t_phase_lo[0] && y <= t_phase_hi[1] && y >= t_phase_lo[1] && z <= t_phase_hi[2] && z >= t_phase_lo[2]){ - if(tphase(i,j,k) == 1.0){ - GL_RHS_r(i,j,k) = 0.0; - } + //when coordinate transformation is OFF, + //alpha = beta = theta = 0. + //Therefore, R_11 = R_22 = R_33 = 1, R_12 = R_13 = R_21 = R_23 = R_31 = R_32 = 0. + + if (Coordinate_Transformation != 1){ + if (R_11 != 1.0 || R_12 != 0.0 || R_13 != 0.0 || + R_21 != 0.0 || R_22 != 1.0 || R_23 != 0.0 || + R_31 != 0.0 || R_32 != 0.0 || R_33 != 1.0 ){ + amrex::Abort("Coordinate transformation is turned OFF, but rotation matrix is not an identity matrix!"); + } + } + + Real dFdPp_Landau = alpha*pOld_p(i,j,k) + beta*std::pow(pOld_p(i,j,k),3.) + FerroX::gamma*std::pow(pOld_p(i,j,k),5.) + + 2. * alpha_12 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) + + 2. * alpha_12 * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),2.) + + 4. * alpha_112 * std::pow(pOld_p(i,j,k),3.) * (std::pow(pOld_q(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) + + 2. * alpha_112 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),4.) + + 2. * alpha_112 * pOld_p(i,j,k) * std::pow(pOld_r(i,j,k),4.) + + 2. * alpha_123 * pOld_p(i,j,k) * std::pow(pOld_q(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); + + Real dFdPq_Landau = alpha*pOld_q(i,j,k) + beta*std::pow(pOld_q(i,j,k),3.) + FerroX::gamma*std::pow(pOld_q(i,j,k),5.) + + 2. * alpha_12 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) + + 2. * alpha_12 * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),2.) + + 4. * alpha_112 * std::pow(pOld_q(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_r(i,j,k),2.)) + + 2. * alpha_112 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),4.) + + 2. * alpha_112 * pOld_q(i,j,k) * std::pow(pOld_r(i,j,k),4.) + + 2. * alpha_123 * pOld_q(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_r(i,j,k),2.); + + Real dFdPr_Landau = alpha*pOld_r(i,j,k) + beta*std::pow(pOld_r(i,j,k),3.) + FerroX::gamma*std::pow(pOld_r(i,j,k),5.) + + 2. * alpha_12 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) + + 2. * alpha_12 * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),2.) + + 4. * alpha_112 * std::pow(pOld_r(i,j,k),3.) * (std::pow(pOld_p(i,j,k),2.) + std::pow(pOld_q(i,j,k),2.)) + + 2. * alpha_112 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),4.) + + 2. * alpha_112 * pOld_r(i,j,k) * std::pow(pOld_q(i,j,k),4.) + + 2. * alpha_123 * pOld_r(i,j,k) * std::pow(pOld_p(i,j,k),2.) * std::pow(pOld_q(i,j,k),2.); + + Real dFdPp_grad = - g11 * DoubleDPDx(pOld_p, mask, i, j, k, dx) + - (g44 + g44_p) * DoubleDPDy(pOld_p, mask, i, j, k, dx) + - (g44 + g44_p) * DoubleDPDz(pOld_p, mask, i, j, k, dx) + - (g12 + g44 - g44_p) * DoubleDPDxDy(pOld_q, mask, i, j, k, dx) // d2P/dxdy + - (g12 + g44 - g44_p) * DoubleDPDxDz(pOld_r, mask, i, j, k, dx); // d2P/dxdz + + Real dFdPq_grad = - g11 * DoubleDPDy(pOld_q, mask, i, j, k, dx) + - (g44 - g44_p) * DoubleDPDx(pOld_q, mask, i, j, k, dx) + - (g44 - g44_p) * DoubleDPDz(pOld_q, mask, i, j, k, dx) + - (g12 + g44 + g44_p) * DoubleDPDxDy(pOld_p, mask, i, j, k, dx) // d2P/dxdy + - (g12 + g44 - g44_p) * DoubleDPDyDz(pOld_r, mask, i, j, k, dx);// d2P/dydz + + //Switch g11 and g44 temporarily for multiphase simulations. This will be generalized later + Real dFdPr_grad = - g44 * ( R_31*R_31*DoubleDPDx(pOld_r, mask, i, j, k, dx) + +R_32*R_32*DoubleDPDy(pOld_r, mask, i, j, k, dx) + +R_33*R_33*DoubleDPDz(pOld_r, mask, i, j, k, dx) + +2.*R_31*R_32*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) + +2.*R_32*R_33*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) + +2.*R_33*R_31*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) + + - (g11 - g44_p) * ( R_11*R_11*DoubleDPDx(pOld_r, mask, i, j, k, dx) + +R_12*R_12*DoubleDPDy(pOld_r, mask, i, j, k, dx) + +R_13*R_13*DoubleDPDz(pOld_r, mask, i, j, k, dx) + +2.*R_11*R_12*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) + +2.*R_12*R_13*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) + +2.*R_13*R_11*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) + + - (g44 - g44_p) * ( R_21*R_21*DoubleDPDx(pOld_r, mask, i, j, k, dx) + +R_22*R_22*DoubleDPDy(pOld_r, mask, i, j, k, dx) + +R_23*R_23*DoubleDPDz(pOld_r, mask, i, j, k, dx) + +2.*R_21*R_22*DoubleDPDxDy(pOld_r, mask, i, j, k, dx) + +2.*R_22*R_23*DoubleDPDyDz(pOld_r, mask, i, j, k, dx) + +2.*R_23*R_21*DoubleDPDxDz(pOld_r, mask, i, j, k, dx)) + + - (g44 + g44_p + g12) * DoubleDPDyDz(pOld_q, mask, i, j, k, dx) // d2P/dydz + - (g44 + g44_p + g12) * DoubleDPDxDz(pOld_p, mask, i, j, k, dx); // d2P/dxdz + + GL_RHS_p(i,j,k) = -1.0 * Gam(i,j,k) * + ( dFdPp_Landau + + dFdPp_grad + - Ep(i,j,k) + //+ DFDx(phi, i, j, k, dx) + ); + + GL_RHS_q(i,j,k) = -1.0 * Gam(i,j,k) * + ( dFdPq_Landau + + dFdPq_grad + - Eq(i,j,k) + //+ DFDy(phi, i, j, k, dx) + ); + + GL_RHS_p(i,j,k) = 0.0; + GL_RHS_q(i,j,k) = 0.0; + GL_RHS_r(i,j,k) = -1.0 * Gam(i,j,k) * + ( dFdPr_Landau + + dFdPr_grad + - Er(i,j,k) + //+ DphiDz(phi, z_hi, z_lo, i, j, k, dx, prob_lo, prob_hi) + ); + + //set t_phase GL_RHS_r to zero so that it stays zero. It is initialized to zero in t-phase as well + //if(x <= t_phase_hi[0] && x >= t_phase_lo[0] && y <= t_phase_hi[1] && y >= t_phase_lo[1] && z <= t_phase_hi[2] && z >= t_phase_lo[2]){ + if(tphase(i,j,k) == 1.0){ + GL_RHS_r(i,j,k) = 0.0; + } }); } } diff --git a/Source/main.cpp b/Source/main.cpp index a9e6740..ad86a0f 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -1,6 +1,8 @@ #include #include #include +#include +#include #ifdef AMREX_USE_EB #include #endif @@ -47,6 +49,7 @@ void main_main (c_FerroX& rFerroX) auto& rGprop = rFerroX.get_GeometryProperties(); auto& geom = rGprop.geom; auto& ba = rGprop.ba; + const BoxArray& nba = amrex::convert(ba,IntVect::TheNodeVector()); auto& dm = rGprop.dm; auto& is_periodic = rGprop.is_periodic; auto& prob_lo = rGprop.prob_lo; @@ -106,31 +109,71 @@ void main_main (c_FerroX& rFerroX) E[dir].define(ba, dm, Ncomp, 0); } - MultiFab PoissonRHS(ba, dm, 1, 0); + MultiFab PoissonRHS(ba, dm, 1, 1); MultiFab PoissonPhi(ba, dm, 1, 1); MultiFab PoissonPhi_Old(ba, dm, 1, 1); MultiFab PoissonPhi_Prev(ba, dm, 1, 1); + MultiFab PoissonPhi_Prev2(ba, dm, 1, 1); MultiFab PhiErr(ba, dm, 1, 1); MultiFab Phidiff(ba, dm, 1, 1); - MultiFab Ex(ba, dm, 1, 0); - 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 charge_den(ba, dm, 1, 0); + //For nodal Poisson + MultiFab Nodal_PoissonRHS(nba, dm, 1, 0); + MultiFab Nodal_PoissonPhi(nba, dm, 1, 1); + MultiFab Nodal_PoissonPhi_Old(nba, dm, 1, 1); + MultiFab Nodal_PoissonPhi_Prev(nba, dm, 1, 1); + MultiFab Nodal_PoissonPhi_Prev2(nba, dm, 1, 1); + MultiFab Nodal_PoissonPhi_BC(nba, dm, 1, 0); + MultiFab APoissonPhi_BC(nba, dm, 1, 0); + MultiFab Nodal_PhiErr(nba, dm, 1, 1); + MultiFab Nodal_Phidiff(nba, dm, 1, 1); + + MultiFab hole_den(ba, dm, 1, 1); + MultiFab e_den(ba, dm, 1, 1); + MultiFab charge_den(ba, dm, 1, 1); + MultiFab MaterialMask(ba, dm, 1, 1); MultiFab tphaseMask(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); + MultiFab angle_alpha(ba, dm, 1, 1); + MultiFab angle_beta(ba, dm, 1, 1); + MultiFab angle_theta(ba, dm, 1, 1); + + for (int dir = 0; dir < AMREX_SPACEDIM; dir++) + { + P_old[dir].setVal(0.); + P_new[dir].setVal(0.); + P_new_pre[dir].setVal(0.); + GL_rhs[dir].setVal(0.); + GL_rhs_pre[dir].setVal(0.); + GL_rhs_avg[dir].setVal(0.); + E[dir].setVal(0.); + } + + PoissonPhi.setVal(0.); + PoissonPhi_Prev.setVal(0.); + PoissonPhi_Prev2.setVal(0.); + PoissonRHS.setVal(0.); + Phidiff.setVal(0.); + tphaseMask.setVal(0.); + angle_alpha.setVal(0.); + angle_beta.setVal(0.); + angle_theta.setVal(0.); + + Nodal_PoissonPhi.setVal(0.); + Nodal_PoissonPhi_Prev.setVal(0.); + Nodal_PoissonPhi_Prev2.setVal(0.); + Nodal_PoissonRHS.setVal(0.); + + hole_den.setVal(0.); + e_den.setVal(0.); + charge_den.setVal(0.); //Initialize material mask 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); + Initialize_Euler_angles(rFerroX, geom, angle_alpha, angle_beta, angle_theta); //cell-centered } else { tphaseMask.setVal(0.); } @@ -151,15 +194,20 @@ void main_main (c_FerroX& rFerroX) amrex::Print() << "contains_SC = " << contains_SC << "\n"; #ifdef AMREX_USE_EB - MultiFab Plt(ba, dm, 18, 0, MFInfo(), *rGprop.pEB->p_factory_union); + MultiFab Plt(ba, dm, 19, 0, MFInfo(), *rGprop.pEB->p_factory_union); #else - MultiFab Plt(ba, dm, 18, 0); + MultiFab Plt(ba, dm, 19, 0); #endif SetPoissonBC(rFerroX, LinOpBCType_2d, all_homogeneous_boundaries, some_functionbased_inhomogeneous_boundaries, some_constant_inhomogeneous_boundaries); // coefficients for solver - MultiFab alpha_cc(ba, dm, 1, 0); + MultiFab alpha_cc(ba, dm, 1, 1); + alpha_cc.setVal(0.); + + MultiFab alpha_nd(nba, dm, 1, 0); + alpha_nd.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);, @@ -220,74 +268,122 @@ void main_main (c_FerroX& rFerroX) pMLMG->setVerbose(mlmg_verbosity); #else - std::unique_ptr p_mlabec; - p_mlabec = std::make_unique(); - p_mlabec->define({geom}, {ba}, {dm}, info); +// std::unique_ptr p_mlabec; +// p_mlabec = std::make_unique(); +// p_mlabec->define({geom}, {ba}, {dm}, info); +// +// //Force singular system to be solvable +// p_mlabec->setEnforceSingularSolvable(false); +// +// p_mlabec->setMaxOrder(linop_maxorder); +// +// 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()); +// +// // set Dirichlet BC by reading in the ghost cell values +// p_mlabec->setLevelBC(amrlev, &PoissonPhi); +// +// // (A*alpha_cc - B * div beta grad) phi = rhs +// p_mlabec->setScalars(-1.0, 1.0); // A = -1.0, B = 1.0; solving (-alpha - div beta grad) phi = RHS +// p_mlabec->setBCoeffs(amrlev, amrex::GetArrOfConstPtrs(beta_face)); +// +// //Declare MLMG object +// pMLMG = std::make_unique(*p_mlabec); +// pMLMG->setVerbose(mlmg_verbosity); + + //Nodal Poisson + std::unique_ptr p_mlndabec; + p_mlndabec = std::make_unique(); + p_mlndabec->define({geom}, {ba}, {dm}, info); //Force singular system to be solvable - p_mlabec->setEnforceSingularSolvable(false); + p_mlndabec->setEnforceSingularSolvable(false); - p_mlabec->setMaxOrder(linop_maxorder); + p_mlndabec->setMaxOrder(linop_maxorder); - p_mlabec->setDomainBC(LinOpBCType_2d[0], LinOpBCType_2d[1]); + p_mlndabec->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, Nodal_PoissonPhi); + //} + //if(some_functionbased_inhomogeneous_boundaries) + //{ + // Fill_FunctionBased_Inhomogeneous_Boundaries(rFerroX, Nodal_PoissonPhi, time); + //} + //Nodal_PoissonPhi.FillBoundary(geom.periodicity()); + //SetPhiBC_z(Nodal_PoissonPhi, n_cell); + // set Dirichlet BC by reading in the ghost cell values - p_mlabec->setLevelBC(amrlev, &PoissonPhi); + //p_mlndabec->setLevelBC(amrlev, &Nodal_PoissonPhi); - // (A*alpha_cc - B * div beta grad) phi = rhs - p_mlabec->setScalars(-1.0, 1.0); // A = -1.0, B = 1.0; solving (-alpha - div beta grad) phi = RHS - p_mlabec->setBCoeffs(amrlev, amrex::GetArrOfConstPtrs(beta_face)); + // (div beta grad) phi = rhs + //p_mlnode->setSigma(amrlev, beta_cc); + + // (A*alpha_cc - B * div beta grad) phi = rhs where A, phi and rhs are nodal MultiFabs and B is cell-centered. + p_mlndabec->setScalars(-1.0, 1.0); // A = -1.0, B = 1.0; solving (-alpha - div beta grad) phi = RHS + //p_mlndabec->setBCoeffs(amrlev, amrex::GetArrOfConstPtrs(beta_face)); + p_mlndabec->setBCoeffs(amrlev, beta_cc); + //Declare MLMG object - pMLMG = std::make_unique(*p_mlabec); + pMLMG = std::make_unique(*p_mlndabec); pMLMG->setVerbose(mlmg_verbosity); + #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 - + //Obtain self consisten Phi and rho Real tol = 1.e-5; Real err = 1.0; int iter = 0; while(err > tol){ - - //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, geom, alpha_cc); + + average_cc_to_nodes(Nodal_PoissonRHS, PoissonRHS, geom); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + average_cc_to_nodes(alpha_nd, alpha_cc, geom); #ifdef AMREX_USE_EB p_mlebabec->setACoeffs(0, alpha_cc); #else - p_mlabec->setACoeffs(0, alpha_cc); + p_mlndabec->setACoeffs(0, alpha_nd); #endif + //Initial guess for phi - PoissonPhi.setVal(0.); + Nodal_PoissonPhi.setVal(0.); + SetPhiBC_z(Nodal_PoissonPhi, n_cell); //Poisson Solve - pMLMG->solve({&PoissonPhi}, {&PoissonRHS}, 1.e-10, -1); - PoissonPhi.FillBoundary(geom.periodicity()); + pMLMG->solve({&Nodal_PoissonPhi}, {&Nodal_PoissonRHS}, 1.e-10, -1); + + Nodal_PoissonPhi.FillBoundary(geom.periodicity()); - // Calculate rho from Phi in SC region - ComputeRho(PoissonPhi, charge_den, e_den, hole_den, MaterialMask); + amrex::average_node_to_cellcenter(PoissonPhi, 0, Nodal_PoissonPhi, 0, 1); + // Calculate rho from Phi in SC region + ComputeRho(PoissonPhi, charge_den, e_den, hole_den, geom, MaterialMask); + if (contains_SC == 0) { // no semiconductor region; set error to zero so the while loop terminates err = 0.; @@ -295,23 +391,24 @@ void main_main (c_FerroX& rFerroX) // 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()); + MultiFab::Copy(Nodal_PhiErr, Nodal_PoissonPhi, 0, 0, 1, 1); + MultiFab::Subtract(Nodal_PhiErr, Nodal_PoissonPhi_Prev, 0, 0, 1, 1); + err = Nodal_PhiErr.norm1(0, geom.periodicity())/Nodal_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); - + MultiFab::Copy(Nodal_PoissonPhi_Prev, Nodal_PoissonPhi, 0, 0, 1, 1); + iter = iter + 1; amrex::Print() << iter << " iterations :: err = " << err << std::endl; } - } + } amrex::Print() << "\n ========= Self-Consistent Initialization of P and Rho Done! ========== \n"<< iter << " iterations to obtain self consistent Phi with err = " << err << std::endl; - + // Calculate E from Phi - ComputeEfromPhi(PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + ComputeEfromPhi(Nodal_PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + + //amrex::average_node_to_cellcenter(PoissonPhi, 0, Nodal_PoissonPhi, 0, 1); // Write a plotfile of the initial data if plot_int > 0 if (plot_int > 0) @@ -337,10 +434,11 @@ void main_main (c_FerroX& rFerroX) MultiFab::Copy(Plt, angle_beta, 0, 15, 1, 0); MultiFab::Copy(Plt, angle_theta, 0, 16, 1, 0); MultiFab::Copy(Plt, Phidiff, 0, 17, 1, 0); + MultiFab::Copy(Plt, alpha_cc, 0, 18, 1, 0); #ifdef AMREX_USE_EB - amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "alpha_cc"}, geom, time, step); #else - amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "alpha_cc"}, geom, time, step); #endif } @@ -348,6 +446,10 @@ void main_main (c_FerroX& rFerroX) int steady_state_step = 1000000; //Initialize to a large number. It will be overwritten by the time step at which steady state condition is satidfied + int sign = 1; //change sign to -1*sign whenever abs(Phi_Bc_hi) == Phi_Bc_hi_max to do triangular wave sweep + int num_Vapp = 0; + Real tiny = 1.e-6; + for (int step = 1; step <= nsteps; ++step) { Real step_strt_time = ParallelDescriptor::second(); @@ -381,48 +483,54 @@ void main_main (c_FerroX& rFerroX) // iterate to compute Phi^{n+1,*} while(err > tol){ - - // Compute RHS of Poisson equation - ComputePoissonRHS(PoissonRHS, P_new_pre, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); + + 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, geom, alpha_cc); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + average_cc_to_nodes(Nodal_PoissonRHS, PoissonRHS, geom); + + average_cc_to_nodes(alpha_nd, alpha_cc, geom); #ifdef AMREX_USE_EB - p_mlebabec->setACoeffs(0, alpha_cc); -#else - p_mlabec->setACoeffs(0, alpha_cc); + p_mlebabec->setACoeffs(0, alpha_cc); +#else + p_mlndabec->setACoeffs(0, alpha_nd); #endif - //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, charge_den, e_den, hole_den, MaterialMask); - - if (contains_SC == 0) { - // no semiconductor region; set error to zero so the while loop terminates - err = 0.; - } else { - - // 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; - } + + //Initial guess for phi + Nodal_PoissonPhi.setVal(0.); + SetPhiBC_z(Nodal_PoissonPhi, n_cell); + + //Poisson Solve + pMLMG->solve({&Nodal_PoissonPhi}, {&Nodal_PoissonRHS}, 1.e-10, -1); + + Nodal_PoissonPhi.FillBoundary(geom.periodicity()); + + amrex::average_node_to_cellcenter(PoissonPhi, 0, Nodal_PoissonPhi, 0, 1); + + // Calculate rho from Phi in SC region + ComputeRho(PoissonPhi, charge_den, e_den, hole_den, geom, MaterialMask); + + if (contains_SC == 0) { + // no semiconductor region; set error to zero so the while loop terminates + err = 0.; + } else { + + // Calculate Error + if (iter > 0){ + MultiFab::Copy(Nodal_PhiErr, Nodal_PoissonPhi, 0, 0, 1, 1); + MultiFab::Subtract(Nodal_PhiErr, Nodal_PoissonPhi_Prev, 0, 0, 1, 1); + err = Nodal_PhiErr.norm1(0, geom.periodicity())/Nodal_PoissonPhi.norm1(0, geom.periodicity()); + } + + MultiFab::Copy(Nodal_PoissonPhi_Prev, Nodal_PoissonPhi, 0, 0, 1, 1); + + iter = iter + 1; + amrex::Print() << iter << " iterations :: err = " << err << std::endl; + } } if (TimeIntegratorOrder == 1) { @@ -452,47 +560,53 @@ void main_main (c_FerroX& rFerroX) // iterate to compute Phi^{n+1} while(err > tol){ - // Compute RHS of Poisson equation - ComputePoissonRHS(PoissonRHS, P_new, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); + 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, geom, alpha_cc); + + average_cc_to_nodes(Nodal_PoissonRHS, PoissonRHS, geom); + + average_cc_to_nodes(alpha_nd, alpha_cc, geom); #ifdef AMREX_USE_EB - p_mlebabec->setACoeffs(0, alpha_cc); -#else - p_mlabec->setACoeffs(0, alpha_cc); + p_mlebabec->setACoeffs(0, alpha_cc); +#else + p_mlndabec->setACoeffs(0, alpha_nd); #endif - - //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, charge_den, e_den, hole_den, MaterialMask); - - if (contains_SC == 0) { - // no semiconductor region; set error to zero so the while loop terminates - err = 0.; - } else { - - // 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; - } + + //Initial guess for phi + Nodal_PoissonPhi.setVal(0.); + SetPhiBC_z(Nodal_PoissonPhi, n_cell); + + //Poisson Solve + pMLMG->solve({&Nodal_PoissonPhi}, {&Nodal_PoissonRHS}, 1.e-10, -1); + + Nodal_PoissonPhi.FillBoundary(geom.periodicity()); + + amrex::average_node_to_cellcenter(PoissonPhi, 0, Nodal_PoissonPhi, 0, 1); + + // Calculate rho from Phi in SC region + ComputeRho(PoissonPhi, charge_den, e_den, hole_den, geom, MaterialMask); + + if (contains_SC == 0) { + // no semiconductor region; set error to zero so the while loop terminates + err = 0.; + } else { + + // Calculate Error + if (iter > 0){ + MultiFab::Copy(Nodal_PhiErr, Nodal_PoissonPhi, 0, 0, 1, 1); + MultiFab::Subtract(Nodal_PhiErr, Nodal_PoissonPhi_Prev, 0, 0, 1, 1); + err = Nodal_PhiErr.norm1(0, geom.periodicity())/Nodal_PoissonPhi.norm1(0, geom.periodicity()); + } + + MultiFab::Copy(Nodal_PoissonPhi_Prev, Nodal_PoissonPhi, 0, 0, 1, 1); + + iter = iter + 1; + amrex::Print() << iter << " iterations :: err = " << err << std::endl; + } } // copy new solution into old solution @@ -504,23 +618,43 @@ void main_main (c_FerroX& rFerroX) } // Check if steady state has reached - MultiFab::Copy(Phidiff, PoissonPhi, 0, 0, 1, 0); - MultiFab::Subtract(Phidiff, PoissonPhi_Old, 0, 0, 1, 0); - Real phi_err = Phidiff.norm0(); + //MultiFab::Copy(Phidiff, PoissonPhi, 0, 0, 1, 0); + //MultiFab::Subtract(Phidiff, PoissonPhi_Old, 0, 0, 1, 0); + + Real phi_max = Nodal_PoissonPhi_Old.norm0(); + + for (MFIter mfi(Nodal_PoissonPhi); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.growntilebox(1); + + const Array4& Phi = Nodal_PoissonPhi.array(mfi); + const Array4& PhiOld = Nodal_PoissonPhi_Old.array(mfi); + const Array4& Phi_err = Nodal_Phidiff.array(mfi); - if (phi_err < phi_tolerance) { + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + Phi_err(i,j,k) = amrex::Math::abs(Phi(i,j,k) - PhiOld(i,j,k)) / phi_max; + }); + } + + Real max_phi_err = Nodal_Phidiff.norm0(); + + if(step > 1){ + if (max_phi_err < phi_tolerance) { steady_state_step = step; inc_step = step; + } } //Copy PoissonPhi to PoissonPhi_Old to calculate difference at the next iteration - MultiFab::Copy(PoissonPhi_Old, PoissonPhi, 0, 0, 1, 0); + MultiFab::Copy(Nodal_PoissonPhi_Old, Nodal_PoissonPhi, 0, 0, 1, 1); - amrex::Print() << "Steady state check : (phi(t) - phi(t-1)).norm0() = " << phi_err << std::endl; + amrex::Print() << "Steady state check : (phi(t) - phi(t-1)).norm0() = " << max_phi_err << std::endl; // Calculate E from Phi - ComputeEfromPhi(PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); + ComputeEfromPhi(Nodal_PoissonPhi, E, angle_alpha, angle_beta, angle_theta, geom, prob_lo, prob_hi); Real step_stop_time = ParallelDescriptor::second() - step_strt_time; @@ -533,6 +667,7 @@ void main_main (c_FerroX& rFerroX) time = time + dt; + amrex::average_node_to_cellcenter(PoissonPhi, 0, Nodal_PoissonPhi, 0, 1); // Write a plotfile of the current data (plot_int was defined in the inputs file) if (plot_int > 0 && (step%plot_int == 0 || step == steady_state_step)) { @@ -556,10 +691,11 @@ void main_main (c_FerroX& rFerroX) MultiFab::Copy(Plt, angle_beta, 0, 15, 1, 0); MultiFab::Copy(Plt, angle_theta, 0, 16, 1, 0); MultiFab::Copy(Plt, Phidiff, 0, 17, 1, 0); + MultiFab::Copy(Plt, alpha_cc, 0, 18, 1, 0); #ifdef AMREX_USE_EB - amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::EB_WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "alpha_cc"}, geom, time, step); #else - amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff"}, geom, time, step); + amrex::WriteSingleLevelPlotfile(pltfile, Plt, {"Px","Py","Pz","Phi","PoissonRHS","Ex","Ey","Ez","holes","electrons","charge","epsilon", "mask", "tphase","alpha", "beta", "theta", "PhiDiff", "alpha_cc"}, geom, time, step); #endif } @@ -567,20 +703,13 @@ void main_main (c_FerroX& rFerroX) { //Update time-dependent Boundary Condition of Poisson's equation - amrex::Print() << "Applied voltage updated at time " << time << ", step = " << step << "\n"; - - Phi_Bc_hi += Phi_Bc_inc; - amrex::Print() << "step = " << step << ", Phi_Bc_hi = " << Phi_Bc_hi << std::endl; - - // Set Dirichlet BC for Phi in z - SetPhiBC_z(PoissonPhi, n_cell); - - // set Dirichlet BC by reading in the ghost cell values -#ifdef AMREX_USE_EB - p_mlebabec->setLevelBC(amrlev, &PoissonPhi); -#else - p_mlabec->setLevelBC(amrlev, &PoissonPhi); -#endif + Phi_Bc_hi += sign*Phi_Bc_inc; + num_Vapp += 1; + if(std::abs(std::abs(Phi_Bc_hi) - Phi_Bc_hi_max) <= tiny) { + sign *= -1; + amrex::Print() << "Direction of voltage sweep is reversed. Phi_Bc_hi = " << Phi_Bc_hi << ", and Phi_Bc_hi_max = " << Phi_Bc_hi_max << std::endl; + } + amrex::Print() << "step = " << step << ", Phi_Bc_hi = " << Phi_Bc_hi << ", num_Vapp = " << num_Vapp << ", sign = " << sign << std::endl; err = 1.0; iter = 0; @@ -588,54 +717,74 @@ void main_main (c_FerroX& rFerroX) // iterate to compute Phi^{n+1} with new Dirichlet value while(err > tol){ - // Compute RHS of Poisson equation - ComputePoissonRHS(PoissonRHS, P_old, charge_den, MaterialMask, angle_alpha, angle_beta, angle_theta, geom); + 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); + + ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, geom, alpha_cc); - 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); + average_cc_to_nodes(Nodal_PoissonRHS, PoissonRHS, geom); - ComputePoissonRHS_Newton(PoissonRHS, PoissonPhi, alpha_cc); + average_cc_to_nodes(alpha_nd, alpha_cc, geom); #ifdef AMREX_USE_EB - p_mlebabec->setACoeffs(0, alpha_cc); -#else - p_mlabec->setACoeffs(0, alpha_cc); + p_mlebabec->setACoeffs(0, alpha_cc); +#else + p_mlndabec->setACoeffs(0, alpha_nd); #endif - - //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, charge_den, e_den, hole_den, MaterialMask); - - if (contains_SC == 0) { - // no semiconductor region; set error to zero so the while loop terminates - err = 0.; - } else { - - // 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; - } + + //Initial guess for phi + Nodal_PoissonPhi.setVal(0.); + SetPhiBC_z(Nodal_PoissonPhi, n_cell); + + //Poisson Solve + pMLMG->solve({&Nodal_PoissonPhi}, {&Nodal_PoissonRHS}, 1.e-10, -1); + + Nodal_PoissonPhi.FillBoundary(geom.periodicity()); + + amrex::average_node_to_cellcenter(PoissonPhi, 0, Nodal_PoissonPhi, 0, 1); + + // Calculate rho from Phi in SC region + ComputeRho(PoissonPhi, charge_den, e_den, hole_den, geom, MaterialMask); + + if (contains_SC == 0) { + // no semiconductor region; set error to zero so the while loop terminates + err = 0.; + } else { + + // Calculate Error + if (iter > 0){ + MultiFab::Copy(Nodal_PhiErr, Nodal_PoissonPhi, 0, 0, 1, 1); + MultiFab::Subtract(Nodal_PhiErr, Nodal_PoissonPhi_Prev, 0, 0, 1, 1); + err = Nodal_PhiErr.norm1(0, geom.periodicity())/Nodal_PoissonPhi.norm1(0, geom.periodicity()); + } + + MultiFab::Copy(Nodal_PoissonPhi_Prev, Nodal_PoissonPhi, 0, 0, 1, 1); + + iter = iter + 1; + amrex::Print() << iter << " iterations :: err = " << err << std::endl; + } } }//end inc_step - - if (voltage_sweep == 0 && step == steady_state_step) break; - if (voltage_sweep == 1 && Phi_Bc_hi > Phi_Bc_hi_max) break; - + + if (voltage_sweep == 0 && step == steady_state_step) { + amrex::Print() << "voltage_sweep == 0 && step == steady_state_step!" << "\n"; + break; + } + if (voltage_sweep == 1 && Phi_Bc_hi > 0. && Phi_Bc_hi - Phi_Bc_hi_max > tiny) { + amrex::Print() << "voltage_sweep == 1 && Phi_Bc_hi > 0. && Phi_Bc_hi - Phi_Bc_hi_max > tiny!" << "\n"; + break; + } + if (voltage_sweep == 1 && Phi_Bc_hi < 0. && -Phi_Bc_hi - Phi_Bc_hi_max > tiny) { + amrex::Print() << "voltage_sweep == 1 && Phi_Bc_hi < 0. && -Phi_Bc_hi - Phi_Bc_hi_max > tiny!" << "\n"; + break; + } + if (voltage_sweep == 1 && num_Vapp == num_Vapp_max) { + amrex::Print() << "voltage_sweep == 1 && num_Vapp == num_Vapp_max!" << "\n"; + break; + } + } // end step // MultiFab memory usage