Nodal poisson

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.
+

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<amrex::Real, AMREX_SPACEDIM>& prob_lo,
                                const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& prob_hi) {
@@ -343,6 +346,9 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
      random_seed = 1;
      pp.query("random_seed",random_seed);
 
+     num_Vapp_max = 1;
+     pp.query("num_Vapp_max",num_Vapp_max);
+
      //stack dimensions in 3D. This is an alternate way of initializing the device geometry, which works in simpler scenarios.
      //A more general way of initializing device geometry is accomplished through masks which use function parsers
 
@@ -440,6 +446,12 @@ void InitializeFerroXNamespace(const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM
      donor_doping = 0.0;
      pp.query("donor_doping",donor_doping);
 
+     intrinsic_carrier_concentration = 9.696e+15;
+     pp.query("intrinsic_carrier_concentration",intrinsic_carrier_concentration);
+
+     use_Fermi_Dirac = 0;
+     pp.query("use_Fermi_Dirac",use_Fermi_Dirac);
+
      Coordinate_Transformation = 0;
      pp.query("Coordinate_Transformation",Coordinate_Transformation);
 

Source/FerroX_namespace.H

@@ -45,6 +45,8 @@ namespace FerroX {
     extern AMREX_GPU_MANAGED amrex::Real T;
     extern AMREX_GPU_MANAGED amrex::Real acceptor_doping;
     extern AMREX_GPU_MANAGED amrex::Real donor_doping;
+    extern AMREX_GPU_MANAGED amrex::Real intrinsic_carrier_concentration;
+    extern AMREX_GPU_MANAGED int use_Fermi_Dirac;
 
     // P and Phi Bc
     extern AMREX_GPU_MANAGED amrex::Real lambda;
@@ -75,6 +77,7 @@ namespace FerroX {
     extern AMREX_GPU_MANAGED amrex::Real Phi_Bc_inc;
     extern AMREX_GPU_MANAGED amrex::Real Phi_Bc_hi_max;
     extern AMREX_GPU_MANAGED amrex::Real phi_tolerance;
+    extern AMREX_GPU_MANAGED int num_Vapp_max;
 
     extern AMREX_GPU_MANAGED int random_seed;
 }

Source/Solver/ChargeDensity.H

@@ -10,5 +10,6 @@ void ComputeRho(MultiFab&      PoissonPhi,
                 MultiFab&      rho,
                 MultiFab&      e_den,
                 MultiFab&      p_den,
+                const       Geometry& geom,
                 const MultiFab& MaterialMask);
 

Source/Solver/ChargeDensity.cpp

@@ -5,6 +5,7 @@ void ComputeRho(MultiFab&      PoissonPhi,
                 MultiFab&      rho,
                 MultiFab&      e_den,
                 MultiFab&      p_den,
+                const       Geometry& geom,
 		const MultiFab& MaterialMask)
 {
     // loop over boxes
@@ -13,8 +14,8 @@ void ComputeRho(MultiFab&      PoissonPhi,
         const Box& bx = mfi.validbox();
 
         // 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<Real>& hole_den_arr = p_den.array(mfi);
         const Array4<Real>& 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());
+ }