Albany: a Trilinos-based PDE code: QCAD_PoissonDirichlet

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009 #include "Sacado_ParameterRegistration.hpp"
00010 #include "Albany_AbstractDiscretization.hpp"
00011 
00012 
00013 namespace QCAD {
00014 
00015 template<typename EvalT,typename Traits>
00016 PoissonDirichlet<EvalT, Traits>::
00017 PoissonDirichlet(Teuchos::ParameterList& p) :
00018   PHAL::Dirichlet<EvalT,Traits>(p)
00019 {
00020   // get parameters from ParameterList
00021   user_value = p.get<RealType>("Dirichlet Value");
00022   materialDB = p.get< Teuchos::RCP<QCAD::MaterialDatabase> >("MaterialDB");
00023   energy_unit_in_eV = p.get<double>("Energy unit in eV");
00024 
00025   Teuchos::ParameterList* psList = p.get<Teuchos::ParameterList*>("Poisson Source Parameter List");
00026 
00027   carrierStatistics = psList->get("Carrier Statistics", "Boltzmann Statistics");
00028   incompIonization = psList->get("Incomplete Ionization", "False");
00029   //dopingDonor = psList->get("Donor Doping", 1e14);
00030   //dopingAcceptor = psList->get("Acceptor Doping", 1e14);
00031   //donorActE = psList->get("Donor Activation Energy", 0.040);
00032   //acceptorActE = psList->get("Acceptor Activation Energy", 0.045);
00033  
00034   temperature = p.get<double>("Temperature"); //To be replaced by SharedParameter evaluator access
00035 
00036   // obtain material or eb name for a given nodeset 
00037   std::string nodeSetName = PHAL::DirichletBase<EvalT,Traits>::nodeSetID;
00038   material = materialDB->getNodeSetParam<std::string>(nodeSetName,"material","");
00039   if (material.length() == 0) {
00040     ebName = materialDB->getNodeSetParam<std::string>(nodeSetName,"elementBlock","");
00041     if (ebName.length() > 0)
00042       material = materialDB->getElementBlockParam<std::string>(ebName,"material","");
00043   }
00044 
00045   // private scaling parameter (note: kbT is used in calculating qPhiRef below)
00046   const double kbBoltz = 8.617343e-05; //[eV/K]
00047   kbT = kbBoltz*temperature;
00048   V0 = kbBoltz*temperature/1.0; // kb*T/q in [V], scaling for potential
00049 
00050   // obtain qPhiRef (energy reference for heterogeneous structures)
00051   if ( (material.length() > 0) || (ebName.length() > 0) )
00052   {
00053     std::string refMtrlName, category;
00054     refMtrlName = materialDB->getParam<std::string>("Reference Material");
00055     category = materialDB->getMaterialParam<std::string>(refMtrlName,"Category");
00056     if (category == "Semiconductor") {
00057       double mdn = materialDB->getMaterialParam<double>(refMtrlName,"Electron DOS Effective Mass");
00058       double mdp = materialDB->getMaterialParam<double>(refMtrlName,"Hole DOS Effective Mass");
00059       double Chi = materialDB->getMaterialParam<double>(refMtrlName,"Electron Affinity");
00060       double Eg0 = materialDB->getMaterialParam<double>(refMtrlName,"Zero Temperature Band Gap");
00061       double alpha = materialDB->getMaterialParam<double>(refMtrlName,"Band Gap Alpha Coefficient");
00062       double beta = materialDB->getMaterialParam<double>(refMtrlName,"Band Gap Beta Coefficient");
00063       
00064       ScalarT Eg = Eg0-alpha*pow(temperature,2.0)/(beta+temperature); // in [eV]
00065       ScalarT Eic = -Eg/2. + 3./4.*kbT*log(mdp/mdn);  // (Ei-Ec) in [eV]
00066       qPhiRef = Chi - Eic;  // (Evac-Ei) in [eV] where Evac = vacuum level
00067     }
00068     else if (category == "Insulator") {
00069       double Chi = materialDB->getMaterialParam<double>(refMtrlName,"Electron Affinity");
00070       qPhiRef = Chi;
00071     }
00072     else if (category == "Metal") {
00073       double workFn = materialDB->getMaterialParam<double>(refMtrlName,"Work Function");
00074       qPhiRef = workFn;
00075     }
00076     else {
00077       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00078         << "Error!  Invalid category " << category 
00079         << " for reference material !" << std::endl);
00080     }
00081   }
00082 }
00083 
00084 
00085 // ****************************************************************************
00086 template<typename EvalT,typename Traits>
00087 void PoissonDirichlet<EvalT, Traits>::
00088 evaluateFields(typename Traits::EvalData dirichletWorkset)
00089 {
00090   ScalarT bcValue = 0.0; 
00091 
00092 
00093   if (material.length() > 0)  
00094   {
00095     std::string category = materialDB->getMaterialParam<std::string>(material,"Category");
00096 
00098     if (category != "Semiconductor") {
00099       double metalWorkFunc = materialDB->getMaterialParam<double>(material,"Work Function");
00100 
00101       // TODO - get electron affinity instead for an insulator??
00102       ScalarT offsetDueToWorkFunc = (metalWorkFunc-qPhiRef)/1.0;  // 1.0 converts from [eV] to [V]
00103 
00104       bcValue = (user_value - offsetDueToWorkFunc) / energy_unit_in_eV;  //[myV]
00105     }
00106   
00108     else {
00109 
00110       // Universal constants
00111       const double kbBoltz = 8.617343e-05;  // [eV/K]
00112       const double eleQ = 1.602176487e-19;  // [C]
00113       const double m0 = 9.10938215e-31;     // [kg]
00114       const double hbar = 1.054571628e-34;  // [J.s]
00115       const double pi = 3.141592654;
00116 
00117   // Temperature-independent material parameters
00118       double mdn, mdp, Tref, Eg0, alpha, beta, Chi;
00119       if (ebName.length() > 0)  {
00120   mdn = materialDB->getElementBlockParam<double>(ebName,"Electron DOS Effective Mass");
00121   mdp = materialDB->getElementBlockParam<double>(ebName,"Hole DOS Effective Mass");
00122   Tref = materialDB->getElementBlockParam<double>(ebName,"Reference Temperature");
00123   
00124   Eg0 = materialDB->getElementBlockParam<double>(ebName,"Zero Temperature Band Gap");
00125   alpha = materialDB->getElementBlockParam<double>(ebName,"Band Gap Alpha Coefficient");
00126   beta = materialDB->getElementBlockParam<double>(ebName,"Band Gap Beta Coefficient");
00127   Chi = materialDB->getElementBlockParam<double>(ebName,"Electron Affinity");
00128       }
00129       else {
00130   mdn = materialDB->getMaterialParam<double>(material,"Electron DOS Effective Mass");
00131   mdp = materialDB->getMaterialParam<double>(material,"Hole DOS Effective Mass");
00132   Tref = materialDB->getMaterialParam<double>(material,"Reference Temperature");
00133   
00134   Eg0 = materialDB->getMaterialParam<double>(material,"Zero Temperature Band Gap");
00135   alpha = materialDB->getMaterialParam<double>(material,"Band Gap Alpha Coefficient");
00136   beta = materialDB->getMaterialParam<double>(material,"Band Gap Beta Coefficient");
00137   Chi = materialDB->getMaterialParam<double>(material,"Electron Affinity");
00138       }
00139 
00140       // Constant prefactor in calculating Nc and Nv in [cm-3]
00141       double NcvFactor = 2.0*pow((kbBoltz*eleQ*m0*Tref)/(2*pi*pow(hbar,2)),3./2.)*1.e-6;
00142           // eleQ converts kbBoltz in [eV/K] to [J/K], 1e-6 converts [m-3] to [cm-3]
00143     
00144       // Strong temperature-dependent material parameters
00145       ScalarT Nc;  // conduction band effective DOS in [cm-3]
00146       ScalarT Nv;  // valence band effective DOS in [cm-3]
00147       ScalarT Eg;  // band gap at T [K] in [eV]
00148     
00149       Nc = NcvFactor*pow(mdn,1.5)*pow(temperature/Tref,1.5);  // in [cm-3]
00150       Nv = NcvFactor*pow(mdp,1.5)*pow(temperature/Tref,1.5); 
00151       Eg = Eg0-alpha*pow(temperature,2.0)/(beta+temperature); // in [eV]
00152     
00153       ScalarT builtinPotential = 0.0; 
00154       std::string dopantType;
00155       if (ebName.length() > 0)
00156   dopantType = materialDB->getElementBlockParam<std::string>(ebName,"Dopant Type","None");
00157       else
00158   dopantType = materialDB->getMaterialParam<std::string>(material,"Dopant Type","None");
00159 
00160       // Intrinsic semiconductor (no doping)
00161       if (dopantType == "None")  
00162       {
00163   // apply charge neutrality (p=n) and MB statistics
00164   builtinPotential = (qPhiRef-Chi-0.5*Eg)/1.0 + 0.5*kbT*log(Nv/Nc)/1.0;
00165   bcValue = (user_value + builtinPotential) / energy_unit_in_eV;  // [myV]
00166       }
00167     
00168       else // Extrinsic semiconductor (doped)
00169       {
00170   double dopingConc, dopantActE;
00171   if (ebName.length() > 0) {
00172     dopingConc = materialDB->getElementBlockParam<double>(ebName,"Doping Value");
00173     dopantActE = materialDB->getElementBlockParam<double>(ebName,"Dopant Activation Energy", 0.045);
00174   }
00175   else {
00176     dopingConc = materialDB->getMaterialParam<double>(material,"Doping Value");
00177     dopantActE = materialDB->getMaterialParam<double>(material,"Dopant Activation Energy", 0.045);
00178   }
00179 
00180   if ((carrierStatistics=="Boltzmann Statistics") && (incompIonization=="False"))
00181     builtinPotential = potentialForMBComplIon(Nc,Nv,Eg,Chi,dopantType,dopingConc);
00182     
00183   else if ((carrierStatistics=="Boltzmann Statistics") && (incompIonization=="True"))
00184     builtinPotential = potentialForMBIncomplIon(Nc,Nv,Eg,Chi,dopantType,dopingConc,dopantActE);
00185 
00186   else if ((carrierStatistics=="Fermi-Dirac Statistics") && (incompIonization=="False"))
00187     builtinPotential = potentialForFDComplIon(Nc,Nv,Eg,Chi,dopantType,dopingConc);
00188   
00189   else if ((carrierStatistics=="0-K Fermi-Dirac Statistics") && (incompIonization=="False"))
00190     builtinPotential = potentialForZeroKFDComplIon(Nc,Nv,Eg,Chi,dopantType,dopingConc);
00191     
00192   // For cases of FD and 0-K FD with incompIonization==True, one needs to 
00193   // numerically solve a non-trivial equation. Since when incompIonization==True
00194   // is enabled, the MB almost always holds, so use potentialForMBIncomplIon.
00195   else  
00196     builtinPotential = potentialForMBIncomplIon(Nc,Nv,Eg,Chi,dopantType,dopingConc,dopantActE);
00197         
00198   bcValue = (user_value + builtinPotential) / energy_unit_in_eV;  // [myV]
00199       }
00200     }
00201   } // end of if (material.length() > 0)
00202     
00204   //    (DBCs are always specified in volts by the user)
00205   else
00206     bcValue = user_value / energy_unit_in_eV;
00207 
00209   PHAL::DirichletBase<EvalT,Traits>::value = bcValue;
00210     
00212   PHAL::Dirichlet<EvalT,Traits>::evaluateFields(dirichletWorkset);
00213 
00214   // std::string nodeSetName = PHAL::DirichletBase<EvalT,Traits>::nodeSetID;
00215   // std::cout << "nodeSetName = " << nodeSetName << ", bcValue = " << bcValue << std::endl;   
00216 
00217 }
00218 
00219 
00220 // *****************************************************************************
00221 template<typename EvalT,typename Traits>
00222 inline typename QCAD::PoissonDirichlet<EvalT,Traits>::ScalarT
00223 QCAD::PoissonDirichlet<EvalT,Traits>::inverseFDIntOneHalf(const ScalarT x)
00224 {
00225   // Use the approximate inverse of 1/2 FD integral by N. G. Nillsson, 
00226   // "An Accurate Approximation of the Generalized Einstein Relation for 
00227   // Degenerate Semiconductors," physica status solidi (a) 19, K75 (1973).
00228   
00229   const double pi = 3.1415926536;
00230   const double a = 0.24;
00231   const double b = 1.08; 
00232 
00233   ScalarT nu = pow(3./4.*sqrt(pi)*x, 2./3.);
00234   ScalarT invFDInt = 0.0; 
00235   
00236   if (x > 0.)
00237     invFDInt = log(x)/(1.-pow(x,2.)) + nu/(1.+pow(a+b*nu,-2.));
00238   else
00239     TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00240       << "Error! x in inverseFDIntOneHalf(x) must be greater than 0 " << std::endl);
00241 
00242   return invFDInt;
00243 }
00244 
00245 
00246 // *****************************************************************************
00247 template<typename EvalT,typename Traits>
00248 typename QCAD::PoissonDirichlet<EvalT,Traits>::ScalarT
00249 QCAD::PoissonDirichlet<EvalT,Traits>::potentialForMBComplIon(const ScalarT &Nc, 
00250    const ScalarT &Nv, const ScalarT &Eg, const double &Chi, const std::string &dopType, const double &dopingConc)
00251 {
00252   ScalarT Cn = Nc*exp((-qPhiRef+Chi)/kbT); 
00253   ScalarT Cp = Nv*exp((qPhiRef-Chi-Eg)/kbT);
00254   ScalarT builtinPotential = 0.0; 
00255   
00256   // for high-T, include n and p in charge neutrality: p=n+Na or p+Nd=n
00257   if ((Cn > 1.e-5) && (Cp > 1.e-5))
00258   {
00259     double signedDopingConc;
00260     if(dopType == "Donor") 
00261       signedDopingConc = dopingConc;
00262     else if(dopType == "Acceptor") 
00263       signedDopingConc = -dopingConc;
00264     else 
00265     {
00266       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00267         << "Error!  Unknown dopant type " << dopType << "!"<< std::endl);
00268     }
00269     ScalarT tmp1 = signedDopingConc/(2.0*Cn);
00270     ScalarT tmp2 = tmp1 + sqrt(pow(tmp1,2.0) + Cp/Cn);
00271     if (tmp2 <= 0.0)
00272     {
00273       if(dopType == "Donor") 
00274         builtinPotential = (qPhiRef-Chi)/1.0 + V0*log(dopingConc/Nc);  
00275       else if(dopType == "Acceptor") 
00276         builtinPotential = (qPhiRef-Chi-Eg)/1.0 - V0*log(dopingConc/Nv);  
00277     }
00278     else
00279       builtinPotential = V0*log(tmp2); 
00280   }
00281   
00282   // for low-T (where Cn=0, Cp=0), consider only p=Na or n=Nd
00283   else
00284   {
00285     if(dopType == "Donor") 
00286       builtinPotential = (qPhiRef-Chi)/1.0 + V0*log(dopingConc/Nc);  
00287     else if(dopType == "Acceptor") 
00288       builtinPotential = (qPhiRef-Chi-Eg)/1.0 - V0*log(dopingConc/Nv);  
00289     else 
00290     {
00291       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00292         << "Error!  Unknown dopant type " << dopType << "!"<< std::endl);
00293     }
00294   }
00295   
00296   return builtinPotential;
00297 }
00298 
00299 
00300 // *****************************************************************************
00301 template<typename EvalT,typename Traits>
00302 typename QCAD::PoissonDirichlet<EvalT,Traits>::ScalarT
00303 QCAD::PoissonDirichlet<EvalT,Traits>::potentialForMBIncomplIon(const ScalarT &Nc, 
00304       const ScalarT &Nv, const ScalarT &Eg, const double &Chi, const std::string &dopType, 
00305       const double &dopingConc, const double &dopantActE )
00306 {
00307   // maximum allowed exponent in an exponential function (unitless)
00308   const double MAX_EXPONENT = 100.0; 
00309   ScalarT builtinPotential;
00310   
00311   // assume n = Nd+ to have an analytical expression (neglect p)  
00312   if(dopType == "Donor") 
00313   {
00314     ScalarT offset = (-dopantActE +qPhiRef -Chi) /1.0;
00315     if (dopantActE/kbT > MAX_EXPONENT)  // exp(dopantActE/kbT) -> +infinity for very small T
00316     {
00317       builtinPotential = offset +V0 *(0.5*log(dopingConc/(2.*Nc)) + dopantActE/(2.*kbT) );
00318     }
00319     else
00320     { 
00321       ScalarT tmp = -1./4.+1./4.*sqrt(1.+8.*dopingConc/Nc*exp(dopantActE/kbT));
00322       builtinPotential = offset + V0*log(tmp);
00323     }  
00324   }
00325   
00326   // assume p = Na- to have an analytical expression (neglect n)
00327   else if(dopType == "Acceptor") 
00328   {
00329     ScalarT offset = (dopantActE +qPhiRef -Chi -Eg) /1.0; 
00330     if (dopantActE/kbT > MAX_EXPONENT)  // exp(dopantActE/kbT) -> +infinity for very small T
00331     {
00332       builtinPotential = offset -V0 *(0.5*log(dopingConc/(4.*Nv)) + dopantActE/(2.*kbT) );
00333     }
00334     else
00335     {
00336       ScalarT tmp = -1./8.+1./8.*sqrt(1.+16.*dopingConc/Nv*exp(dopantActE/kbT));
00337       builtinPotential = offset - V0*log(tmp);
00338     }
00339   }
00340 
00341   else 
00342   {
00343     TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00344       << "Error!  Unknown dopant type " << dopType << "!"<< std::endl);
00345   }
00346 
00347   return builtinPotential;
00348 }
00349 
00350 
00351 // *****************************************************************************
00352 template<typename EvalT,typename Traits>
00353 typename QCAD::PoissonDirichlet<EvalT,Traits>::ScalarT
00354 QCAD::PoissonDirichlet<EvalT,Traits>::potentialForFDComplIon(const ScalarT &Nc, 
00355     const ScalarT &Nv, const ScalarT &Eg, const double &Chi, const std::string &dopType, const double &dopingConc)
00356 {
00357   ScalarT builtinPotential;
00358     
00359   // assume n = Nd to have an analytical expression (neglect p)
00360   if(dopType == "Donor") 
00361   {
00362     ScalarT invFDInt = inverseFDIntOneHalf(dopingConc/Nc);
00363     builtinPotential = (qPhiRef-Chi)/1.0 + V0*invFDInt;
00364   }
00365   
00366   // assume p = Na to have an analytical expression (neglect n)
00367   else if(dopType == "Acceptor") 
00368   {
00369     ScalarT invFDInt = inverseFDIntOneHalf(dopingConc/Nv);
00370     builtinPotential = (qPhiRef-Chi-Eg)/1.0 - V0*invFDInt;
00371   }
00372   else {
00373     TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00374       << "Error!  Unknown dopant type " << dopType << "!"<< std::endl);
00375   }
00376     
00377   return builtinPotential;
00378 }
00379 
00380 
00381 // *****************************************************************************
00382 template<typename EvalT,typename Traits>
00383 typename QCAD::PoissonDirichlet<EvalT,Traits>::ScalarT
00384 QCAD::PoissonDirichlet<EvalT,Traits>::potentialForZeroKFDComplIon(const ScalarT &Nc, 
00385     const ScalarT &Nv, const ScalarT &Eg, const double &Chi, const std::string &dopType, const double &dopingConc)
00386 {
00387   const double pi = 3.1415926536;
00388 
00389   ScalarT builtinPotential;
00390     
00391   // assume n = Nd to have an analytical expression (neglect p)
00392   if(dopType == "Donor") 
00393   {
00394     if (dopingConc < Nc)  // Fermi level is below conduction band (due to doping)
00395       builtinPotential = potentialForFDComplIon(Nc,Nv,Eg,Chi,dopType,dopingConc);
00396     else  // Fermi level is in conduction band (degenerate)
00397     { 
00398       ScalarT invFDInt = pow(3./4.*sqrt(pi)*(dopingConc/Nc),2./3.);
00399       builtinPotential = (qPhiRef-Chi)/1.0+ V0*invFDInt;
00400     }
00401   }
00402   
00403   // assume p = Na to have an analytical expression (neglect n)
00404   else if(dopType == "Acceptor") 
00405   {
00406     if (dopingConc < Nv)  // Fermi level is above valence band (due to doping) 
00407       builtinPotential = potentialForFDComplIon(Nc,Nv,Eg,Chi,dopType,dopingConc); 
00408     else  // Fermi level is in valence band (degenerate)
00409     {  
00410       ScalarT invFDInt = pow(3./4.*sqrt(pi)*(dopingConc/Nv),2./3.);
00411       builtinPotential = (qPhiRef-Chi-Eg)/1.0- V0*invFDInt;
00412     }
00413   }
00414   
00415   else 
00416   {
00417     TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00418       << "Error!  Unknown dopant type " << dopType << "!"<< std::endl);
00419   }
00420     
00421   return builtinPotential;
00422 }
00423 
00424 }  //end of QCAD namespace
QCAD_PoissonDirichlet_Def.hpp
