• Main Page
  • Related Pages
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

QCAD_CoupledPSJacobian.cpp

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 "QCAD_CoupledPSJacobian.hpp"
00008 #include "Teuchos_ParameterListExceptions.hpp"
00009 #include "Teuchos_TestForException.hpp"
00010 #include "Epetra_LocalMap.h"
00011 
00012 //Forward Prototypes for utility functions
00013 double n_prefactor(int numDims, int valleyDegeneracyFactor, double T, double length_unit_in_m, double energy_unit_in_eV, double effmass);
00014 double n_weight_factor(double eigenvalue, int numDims, double T, double energy_unit_in_eV);
00015 double dn_weight_factor(double eigenvalue, int numDims, double T, double energy_unit_in_eV);
00016 double compute_FDIntOneHalf(const double x);
00017 double compute_dFDIntOneHalf(const double x);
00018 double compute_FDIntMinusOneHalf(const double x);
00019 double compute_dFDIntMinusOneHalf(const double x);
00020 
00021 QCAD::CoupledPSJacobian::CoupledPSJacobian(int nEigenvals, 
00022              const Teuchos::RCP<const Epetra_Map>& discretizationMap, 
00023              const Teuchos::RCP<const Epetra_Map>& fullPSMap, 
00024              const Teuchos::RCP<const Epetra_Comm>& comm,
00025              int dim, int valleyDegen, double temp, 
00026              double lengthUnitInMeters, double energyUnitInElectronVolts,
00027              double effMass, double conductionBandOffset)
00028 {
00029   discMap = discretizationMap;
00030   domainMap = rangeMap = fullPSMap;
00031   myComm = comm;
00032   bUseTranspose = false;
00033   bInitialized = false;
00034 
00035   numDims = dim;
00036   valleyDegenFactor = valleyDegen;
00037   temperature = temp;
00038   length_unit_in_m = lengthUnitInMeters;
00039   energy_unit_in_eV = energyUnitInElectronVolts;
00040   effmass = effMass;
00041   offset_to_CB = conductionBandOffset;
00042 }
00043 
00044 QCAD::CoupledPSJacobian::~CoupledPSJacobian()
00045 {
00046 }
00047 
00048 
00050 void QCAD::CoupledPSJacobian::initialize(const Teuchos::RCP<Epetra_CrsMatrix>& poissonJac, const Teuchos::RCP<Epetra_CrsMatrix>& schrodingerJac, 
00051            const Teuchos::RCP<Epetra_CrsMatrix>& massMx,
00052            const Teuchos::RCP<Epetra_Vector>& neg_eigenvals, const Teuchos::RCP<const Epetra_MultiVector>& eigenvecs)
00053 {
00054   // Set member variables
00055   poissonJacobian = poissonJac;
00056   schrodingerJacobian = schrodingerJac;
00057   massMatrix = massMx;
00058   neg_eigenvalues = neg_eigenvals;
00059   psiVectors = eigenvecs;
00060 
00061   // Fill vectors that will be needed in Apply, but do so here so 
00062   //   it is only done once per evalModel call, not each time Apply is called
00063 
00064   int nEigenvalues = neg_eigenvalues->GlobalLength();
00065   int num_discMap_myEls = discMap->NumMyElements();
00066 
00067   // dn_dPsi : vectors of dn/dPsi[i] values
00068   dn_dPsi = Teuchos::rcp(new Epetra_MultiVector( *psiVectors ));
00069   for(int i=0; i<nEigenvalues; i++) {
00070     (*dn_dPsi)(i)->Scale( n_prefactor(numDims, valleyDegenFactor, temperature, length_unit_in_m, energy_unit_in_eV, effmass) 
00071         * 2 * n_weight_factor( -(*neg_eigenvalues)[i], numDims, temperature, energy_unit_in_eV), *((*psiVectors)(i)));
00072   }
00073 
00074   // dn_dEval : vectors of dn/dEval[i]
00075   dn_dEval = Teuchos::rcp(new Epetra_MultiVector( *discMap, nEigenvalues ));
00076   for(int i=0; i<nEigenvalues; i++) {
00077     double prefactor = n_prefactor(numDims, valleyDegenFactor, temperature, length_unit_in_m, energy_unit_in_eV, effmass);
00078     double dweight = dn_weight_factor(-(*neg_eigenvalues)[i], numDims, temperature, energy_unit_in_eV);
00079 
00080     //DEBUG
00081     //double eps = 1e-7;
00082     //double dweight = (n_weight_factor( -(*neg_eigenvalues)[i] + eps, numDims, temperature, energy_unit_in_eV) - 
00083     //          n_weight_factor( -(*neg_eigenvalues)[i], numDims, temperature, energy_unit_in_eV)) / eps; 
00084     //const double kbBoltz = 8.617343e-05;
00085     //std::cout << "DEBUG: dn_dEval["<<i<<"] dweight arg = " <<  (*neg_eigenvalues)[i]/(kbBoltz*temperature) << std::endl;  // in [eV]
00086     //std::cout << "DEBUG: dn_dEval["<<i<<"] factor = " <<  prefactor * dweight << std::endl;  // in [eV]
00087     //DEBUG
00088 
00089     for(int k=0; k<num_discMap_myEls; k++)
00090       ( *((*dn_dEval)(i)) )[k] =  prefactor * pow( (*((*psiVectors)(i)))[k], 2.0 ) * dweight;
00091     //(*dn_dEval)(i)->Print( std::cout << "DEBUG: dn_dEval["<<i<<"]:" << std::endl );
00092   }
00093 
00094   // mass matrix multiplied by Psi: M*Psi
00095   M_Psi = Teuchos::rcp(new Epetra_MultiVector( *discMap, nEigenvalues ));
00096   massMatrix->Multiply(false, *psiVectors, *M_Psi);
00097 
00098   // transpose(mass matrix) multiplied by Psi: MT*Psi
00099   MT_Psi = Teuchos::rcp(new Epetra_MultiVector( *discMap, nEigenvalues ));
00100   massMatrix->Multiply(true, *psiVectors, *MT_Psi);
00101 
00102   // Prepare the distributed and local eigenvalue maps
00103   int my_nEigenvals = domainMap->NumMyElements() - num_discMap_myEls * (1+nEigenvalues);
00104   dist_evalMap  = Teuchos::rcp(new Epetra_Map(nEigenvalues, my_nEigenvals, 0, *myComm));
00105   local_evalMap = Teuchos::rcp(new Epetra_LocalMap(nEigenvalues, 0, *myComm));
00106   eval_importer = Teuchos::rcp(new Epetra_Import(*local_evalMap, *dist_evalMap));
00107              
00108   // Allocate temporary space for local storage of the eigenvalue part of the x vectors passed to Apply()
00109   x_neg_evals_local = Teuchos::rcp(new Epetra_Vector(*local_evalMap));
00110 
00111 }
00112 
00113 
00115 int QCAD::CoupledPSJacobian::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00116 { 
00117   //int disc_nGlobalElements = discMap->NumGlobalElements();
00118   int disc_nMyElements = discMap->NumMyElements();
00119   int nEigenvals = neg_eigenvalues->GlobalLength();
00120 
00121   Teuchos::RCP<Epetra_Vector> x_poisson, y_poisson;
00122   Teuchos::RCP<Epetra_MultiVector> x_schrodinger, y_schrodinger;
00123   Teuchos::RCP<Epetra_Vector> x_neg_evals, y_neg_evals;
00124   double *x_data, *y_data;
00125 
00126   Epetra_Vector tempVec(*discMap); //we could allocate this in initialize instead of on stack...
00127   Epetra_Vector tempVec2(*discMap); //we could allocate this in initialize instead of on stack...
00128 
00129   for(int i=0; i < X.NumVectors(); i++) {
00130 
00131     // Split X(i) and Y(i) into vector views for Poisson and Schrodinger parts
00132     if(X(i)->ExtractView(&x_data) != 0 || Y(i)->ExtractView(&y_data) != 0) 
00133       TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00134          "Error!  QCAD::CoupledPSJacobian -- cannot extract vector views");
00135 
00136     x_poisson = Teuchos::rcp(new Epetra_Vector(::View, *discMap, &x_data[0]));
00137     x_schrodinger = Teuchos::rcp(new Epetra_MultiVector(::View, *discMap, &x_data[disc_nMyElements], disc_nMyElements, nEigenvals));
00138     x_neg_evals = Teuchos::rcp(new Epetra_Vector(::View, *dist_evalMap, &x_data[(1+nEigenvals)*disc_nMyElements]));
00139 
00140     y_poisson = Teuchos::rcp(new Epetra_Vector(::View, *discMap, &y_data[0]));
00141     y_schrodinger = Teuchos::rcp(new Epetra_MultiVector(::View, *discMap, &y_data[disc_nMyElements], disc_nMyElements, nEigenvals));
00142     y_neg_evals = Teuchos::rcp(new Epetra_Vector(::View, *dist_evalMap, &y_data[(1+nEigenvals)*disc_nMyElements]));
00143 
00144     // Communicate all the x_evals to every processor, since all parts of the mesh need them
00145     x_neg_evals_local->Import(*x_neg_evals, *eval_importer, Insert);
00146 
00147 
00148     // Jacobian Matrix is:
00149     //
00150     //                   Phi                    Psi[i]                            -Eval[i]
00151     //          | ------------------------------------------------------------------------------------------|
00152     //          |                      |                             |                                      |
00153     // Poisson  |    Jac_poisson       |   M*diag(dn/d{Psi[i](x)})   |        -M*col(dn/dEval[i])           |
00154     //          |                      |                             |                                      |
00155     //          | ------------------------------------------------------------------------------------------|
00156     //          |                      |                             |                                      |
00157     // Schro[j] |  M*diag(-Psi[j](x))  | delta(i,j)*[ H-Eval[i]*M ]  |        delta(i,j)*M*Psi[i](x)        |    
00158     //          |                      |                             |                                      |
00159     //          | ------------------------------------------------------------------------------------------|
00160     //          |                      |                             |                                      |
00161     // Norm[j]  |    0                 | -delta(i,j)*(M+M^T)*Psi[i]  |                   0                  |
00162     //          |                      |                             |                                      |
00163     //          | ------------------------------------------------------------------------------------------|
00164     //
00165     //
00166     //   Where:
00167     //       n = quantum density function which depends on dimension
00168     
00169     // Scratch:  val = sum_ij x_i * M_ij * x_j
00170     //         dval/dx_k = sum_j!=k M_kj * x_j + sum_i!=k x_i * M_ik + 2* M_kk * x_k
00171     //                   = sum_i (M_ki + M_ik) * x_i
00172     //                   = sum_i (M + M^T)_ki * x_i  == k-th el of (M + M^T)*x
00173     //   So d(x*M*x)/dx = (M+M^T)*x in matrix form
00174 
00175     // Do multiplication block-wise
00176 
00177     // y_poisson =     
00178     poissonJacobian->Multiply(false, *x_poisson, *y_poisson); // y_poisson = Jac_poisson * x_poisson
00179     for(int i=0; i<nEigenvals; i++) {
00180       tempVec.Multiply(1.0, *((*dn_dPsi)(i)), *((*x_schrodinger)(i)), 0.0);  //tempVec = dn_dPsi[i] @ x_schrodinger[i]   (@ = el-wise product)
00181       massMatrix->Multiply(false, tempVec, tempVec2);  
00182       y_poisson->Update(1.0, tempVec2, 1.0);                                 // y_poisson += M * (dn_dPsi[i] @ x_schrodinger[i])
00183 
00184       tempVec.Update( -(*x_neg_evals_local)[i], *((*dn_dEval)(i)), 0.0); // tempVec = -dn_dEval[i] * scalar(x_neg_eval[i])
00185       massMatrix->Multiply(false, tempVec, tempVec2);
00186       y_poisson->Update( 1.0, tempVec2, 1.0); // y_poisson += M * (-dn_dEval[i] * scalar(x_neg_eval[i]))
00187 
00188       //OLD: y_poisson->Multiply(1.0, *((*dn_dPsi)(i)), *((*x_schrodinger)(i)), 1.0); // y_poisson += dn_dPsi[i] @ x_schrodinger[i]   (@ = el-wise product)
00189       //OLD: y_poisson->Update( -(*x_neg_evals_local)[i], *((*dn_dEval)(i)), 1.0); // y_poisson += -dn_dEval[i] * scalar(x_neg_eval[i])
00190     }
00191 
00192     // y_schrodinger =
00193     for(int j=0; j<nEigenvals; j++) {
00194 
00195       schrodingerJacobian->Multiply(false, *((*x_schrodinger)(j)), *((*y_schrodinger)(j)) ); // y_schrodinger[j] = H * x_schrodinger[j]
00196       massMatrix->Multiply(false, *((*x_schrodinger)(j)), tempVec );
00197       (*y_schrodinger)(j)->Update( (*neg_eigenvalues)[j], tempVec, 1.0); // y_schrodinger[j] += -eval[j] * M * x_schrodinger[j] 
00198 
00199       //tempVec.PutScalar(offset_to_CB);
00200       //tempVec.Update( -1.0, *((*psiVectors)(j)), 1.0);
00201       tempVec.Multiply( -1.0, *((*psiVectors)(j)), *x_poisson, 0.0); // tempVec = (-Psi[j]) @ x_poisson
00202       massMatrix->Multiply(false, tempVec, tempVec2);  
00203       (*y_schrodinger)(j)->Update( 1.0, tempVec2, 1.0); // y_schrodinger[j] += M * ( (-Psi[j]) @ x_poisson )
00204 
00205       //OLD: (*y_schrodinger)(j)->Multiply( -1.0, *((*psiVectors)(j)), *x_poisson, 1.0); // y_schrodinger[j] += (-Psi[j]) @ x_poisson
00206       (*y_schrodinger)(j)->Update( (*x_neg_evals_local)[j], *((*M_Psi)(j)), 1.0); // y_schrodinger[j] += M*Psi[j] * scalar(x_neg_eval[j])
00207     }
00208       
00209     // y_neg_evals =
00210     std::vector<double> y_neg_evals_local(nEigenvals);
00211     for(int j=0; j<nEigenvals; j++) {
00212       tempVec.Update(-1.0, *((*M_Psi)(j)), -1.0, *((*MT_Psi)(j)), 0.0); //tempVec = -(M*Psi[j] + MT*Psi[j]) = -(M+MT)*Psi[j]
00213       tempVec.Dot( *((*x_schrodinger)(j)), &y_neg_evals_local[j] );
00214     }
00215       // Fill elements of y_neg_evals that belong to this processor
00216     int my_nEigenvals = dist_evalMap->NumMyElements();
00217     std::vector<int> eval_global_elements(my_nEigenvals);
00218     dist_evalMap->MyGlobalElements(&eval_global_elements[0]);
00219     for(int i=0; i<my_nEigenvals; i++)
00220       (*y_neg_evals)[i] = y_neg_evals_local[eval_global_elements[i]];
00221 
00222   }
00223   return 0;
00224 }
00225 
00227 int QCAD::CoupledPSJacobian::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00228 {
00229     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00230              std::endl << "QCAD::CoupledPSJacobian::ApplyInverse is not implemented." << std::endl);
00231     return 0;
00232 }
00233 
00234 
00235 
00236 
00237 // *****************************************************************************
00238 // Quantum Density function
00239 //
00240 //  In all dimensions, it has the form:
00241 //  n( x, Eval[*], Psi[*]) = prefactor * sum_k( |Psi[k](x)|^2 * weight_factor(Eval[k]) )
00242 //
00243 //  Thus (assuming Psi's are real):
00244 //     dn/d{Psi[i](x)} = prefactor * 2 * Psi[i] * weight_factor(Eval[i])
00245 //     dn/dEval[i]     = prefactor * Psi[i]^2 * d{weight_factor(Eval[i])}/dEval[i]
00246 //
00247 //  Below we define functions for prefactor, weight_factor, and d{weight_factor(E)}/dE
00248 // *****************************************************************************
00249 
00250 // CONSTANTS -- copied from QCAD_PoissonSource.cpp
00251 
00252 // Boltzmann constant in [eV/K]
00253 const double kbBoltz = 8.617343e-05;
00254 
00255 // vacuum permittivity in [C/(V.cm)]
00256 const double eps0 = 8.854187817e-12*0.01;
00257 
00258 // electron elemental charge in [C]
00259 const double eleQ = 1.602176487e-19; 
00260 
00261 // vacuum electron mass in [kg]
00262 const double m0 = 9.10938215e-31; 
00263 
00264 // reduced planck constant in [J.s]
00265 const double hbar = 1.054571628e-34; 
00266 
00267 // pi constant (unitless)
00268 const double pi = 3.141592654; 
00269 
00270 // unit conversion factors
00271 const double eVPerJ = 1.0/eleQ; 
00272 const double cm2Perm2 = 1.0e4; 
00273 
00274 // maximum allowed exponent in an exponential function (unitless)
00275 const int MAX_EXPONENT = 100.0;
00276 
00277 
00278 double n_prefactor(int numDims, int valleyDegeneracyFactor, double T, double length_unit_in_m, double energy_unit_in_eV, double effmass)
00279 {
00280   // Scaling factors
00281   double X0 = length_unit_in_m/1e-2; // length scaling to get to [cm] (structure dimension in [um])
00282   double kbT = kbBoltz*T;  // in [eV]
00283   double eDenPrefactor;
00284 
00285   // unit factor applied to poisson source rhs1
00286   double poissonSourcePrefactor = (eleQ*X0*X0)/eps0 * 1.0/energy_unit_in_eV;
00287 
00288   // compute quantum electron density prefactor according to dimensionality
00289   switch (numDims)
00290   {
00291     case 1: // 1D wavefunction (1D confinement)
00292     {  
00293       // effmass = in-plane effective mass (y-z plane when the 1D wavefunc. is along x)
00294       // For Delta2-band (or valley), it is the transverse effective mass (0.19). 
00295         
00296       // 2D density of states in [#/(eV.cm^2)] where 2D is the unconfined y-z plane
00297       // dos2D below includes the spin degeneracy of 2
00298       double dos2D = effmass*m0/(pi*hbar*hbar*eVPerJ*cm2Perm2); 
00299         
00300       // subband-independent prefactor in calculating electron density
00301       // X0 is used to scale wavefunc. squared from [um^-1] or [nm^-1] to [cm^-1]
00302       eDenPrefactor = valleyDegeneracyFactor*dos2D*kbT/X0;
00303       break;
00304     }
00305     case 2: // 2D wavefunction (2D confinement)
00306     {
00307       // mUnconfined = effective mass in the unconfined direction (z dir. when the 2D wavefunc. is in x-y plane)
00308       // For Delta2-band and assume SiO2/Si interface parallel to [100] plane, mUnconfined=0.19. 
00309       double mUnconfined = effmass;
00310         
00311       // n1D below is a factor that is part of the line electron density 
00312       // in the unconfined dir. and includes spin degeneracy and in unit of [cm^-1]
00313       double n1D = sqrt(2.0*mUnconfined*m0*kbT/(pi*hbar*hbar*eVPerJ*cm2Perm2));
00314         
00315       // subband-independent prefactor in calculating electron density
00316       // X0^2 is used to scale wavefunc. squared from [um^-2] or [nm^-2] to [cm^-2]
00317       eDenPrefactor = valleyDegeneracyFactor*n1D/pow(X0,2.);
00318       break;
00319     }
00320     case 3: // 3D wavefunction (3D confinement)
00321     { 
00322       //degeneracy factor
00323       int spinDegeneracyFactor = 2;
00324       double degeneracyFactor = spinDegeneracyFactor * valleyDegeneracyFactor;
00325         
00326       // subband-independent prefactor in calculating electron density
00327       // X0^3 is used to scale wavefunc. squared from [um^-3] or [nm^-3] to [cm^-3]
00328       eDenPrefactor = degeneracyFactor/pow(X0,3.);
00329       break;
00330     }
00331     default:
00332       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter,
00333         std::endl << "Error!  Invalid number of dimensions " << numDims << "!"<< std::endl);
00334       eDenPrefactor = 0.0;
00335       break; 
00336   }
00337   return eDenPrefactor * poissonSourcePrefactor;
00338 }
00339 
00340 
00341 //Note: assumes fermi level Ef == 0, so this is assumed to be true in the quantum region... perhaps pass in as a param later?
00342 double n_weight_factor(double eigenvalue, int numDims, double T, double energy_unit_in_eV)
00343 {
00344   double Ef = 0.0; //Fermi level of the quantum region
00345   double kbT = kbBoltz*T / energy_unit_in_eV;  // in [myV]
00346 
00347   switch (numDims)
00348   {
00349     case 1: // 1D wavefunction (1D confinement)
00350     {
00351       double tmpArg = (Ef-eigenvalue)/kbT;
00352       double logFunc; 
00353       if (tmpArg > MAX_EXPONENT)
00354   logFunc = tmpArg;  // exp(tmpArg) blows up for large positive tmpArg, leading to bad derivative
00355       else
00356   logFunc = log(1.0 + exp(tmpArg));
00357       return logFunc;
00358     }
00359     case 2: // 2D wavefunction (2D confinement)
00360     {
00361       double inArg = (Ef-eigenvalue)/kbT;
00362       return compute_FDIntMinusOneHalf(inArg);
00363     }
00364     case 3: // 3D wavefunction (3D confinement)
00365     { 
00366       double tmpArg = (eigenvalue-Ef)/kbT;
00367       double fermiFactor; 
00368         
00369       if (tmpArg > MAX_EXPONENT) 
00370   fermiFactor = exp(-tmpArg);  // use Boltzmann statistics for large positive tmpArg,
00371       else                           // as the Fermi statistics leads to bad derivative        
00372   fermiFactor = 1.0/( exp(tmpArg) + 1.0 ); 
00373       return fermiFactor;
00374     }      
00375     default:
00376       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter,
00377         std::endl << "Error!  Invalid number of dimensions " << numDims << "!"<< std::endl);
00378       break; 
00379   }
00380   return 0.0;
00381 }
00382 
00383 double dn_weight_factor(double eigenvalue, int numDims, double T, double energy_unit_in_eV)
00384 {
00385   double Ef = 0.0; //Fermi level of the quantum region
00386   double kbT = kbBoltz*T / energy_unit_in_eV;  // in [myV]
00387 
00388   switch (numDims)
00389   {
00390     case 1: // 1D wavefunction (1D confinement)
00391     {
00392       double tmpArg = (Ef-eigenvalue)/kbT;
00393       double dlogFunc; 
00394       if (tmpArg > MAX_EXPONENT)
00395   dlogFunc = -1.0/kbT;
00396       else
00397   dlogFunc = -1.0/(1.0 + exp(tmpArg)) * exp(tmpArg)/kbT;
00398       return dlogFunc;
00399     }
00400     case 2: // 2D wavefunction (2D confinement)
00401     {
00402       double inArg = (Ef-eigenvalue)/kbT;
00403       return compute_dFDIntMinusOneHalf(inArg) * -1.0/kbT;
00404     }
00405     case 3: // 3D wavefunction (3D confinement)
00406     { 
00407       double tmpArg = (eigenvalue-Ef)/kbT;
00408       double dfermiFactor; 
00409         
00410       if (tmpArg > MAX_EXPONENT) 
00411   dfermiFactor = exp(-tmpArg) * -1.0/kbT;  // use Boltzmann statistics for large positive tmpArg,
00412       else                                       // as the Fermi statistics leads to bad derivative        
00413   dfermiFactor = -1.0/pow( exp(tmpArg) + 1.0, 2) * exp(tmpArg)/kbT ; 
00414       return dfermiFactor;
00415     }      
00416     default:
00417       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter,
00418         std::endl << "Error!  Invalid number of dimensions " << numDims << "!"<< std::endl);
00419       break; 
00420   }
00421   return 0.0;
00422 }
00423 
00424 
00425 double compute_FDIntOneHalf(const double x)
00426 {
00427    // Use the approximate 1/2 FD integral by D. Bednarczyk and J. Bednarczyk, 
00428    // "The approximation of the Fermi-Dirac integral F_{1/2}(x),"
00429    // Physics Letters A, vol.64, no.4, pp.409-410, 1978. The approximation 
00430    // has error < 4e-3 in the entire x range.  
00431    
00432    double fdInt; 
00433    if (x >= -50.0)
00434    {
00435      fdInt = pow(x,4.) + 50. + 33.6*x*(1.-0.68*exp(-0.17*pow((x+1.),2.0)));
00436      fdInt = pow((exp(-x) + (3./4.*sqrt(pi)) * pow(fdInt, -3./8.)),-1.0);
00437    }      
00438    else
00439      fdInt = exp(x); // for x<-50, the 1/2 FD integral is well approximated by exp(x)
00440      
00441    return fdInt;
00442 }
00443 
00444 double compute_dFDIntOneHalf(const double x)
00445 {
00446    // Use the approximate 1/2 FD integral by D. Bednarczyk and J. Bednarczyk, 
00447    // "The approximation of the Fermi-Dirac integral F_{1/2}(x),"
00448    // Physics Letters A, vol.64, no.4, pp.409-410, 1978. The approximation 
00449    // has error < 4e-3 in the entire x range.  
00450    
00451    double dfdInt, v,dv,u,du; 
00452    if (x >= -50.0)
00453    {
00454      v = pow(x,4.) + 50. + 33.6*x*(1.-0.68*exp(-0.17*pow((x+1.),2.0)));
00455      dv = 4.*pow(x,3.) + 33.6 * ( x*(-0.68*-0.17*2*(x+1.)*exp(-0.17*pow((x+1.),2.0))) + (1.-0.68*exp(-0.17*pow((x+1.),2.0))) );
00456      u = exp(-x) + (3./4.*sqrt(pi)) * pow(v, -3./8.);
00457      du = -exp(-x) + (3./4.*sqrt(pi)) * -3./8. * pow(v, -11./8.) * dv;
00458      dfdInt = -1.0 * pow( u ,-2.0) * du;
00459    }      
00460    else
00461      dfdInt = exp(x); // for x<-50, the 1/2 FD integral is well approximated by exp(x)
00462      
00463    return dfdInt;
00464 }
00465 
00466 
00467 
00468 double compute_FDIntMinusOneHalf(const double x)
00469 {
00470    // Use the approximate -1/2 FD integral by P. Van Halen and D. L. Pulfrey, 
00471    // "Accurate, short series approximations to Fermi-Dirac integrals of order 
00472    // -/2, 1/2, 1, 3/2, 2, 5/2, 3, and 7/2" J. Appl. Phys. 57, 5271 (1985) 
00473    // and its Erratum in J. Appl. Phys. 59, 2264 (1986). The approximation
00474    // has error < 1e-5 in the entire x range.  
00475    
00476    double fdInt; 
00477    double a1, a2, a3, a4, a5, a6, a7; 
00478    if (x <= 0.)  // eqn.(4) in the reference
00479    {
00480      a1 = 0.999909;  // Table I in Erratum
00481      a2 = 0.706781;
00482      a3 = 0.572752;
00483      a4 = 0.466318;
00484      a5 = 0.324511;
00485      a6 = 0.152889;
00486      a7 = 0.033673;
00487      fdInt = a1*exp(x)-a2*exp(2.*x)+a3*exp(3.*x)-a4*exp(4.*x)+a5*exp(5.*x)-a6*exp(6.*x)+a7*exp(7.*x);
00488    }
00489    else if (x >= 5.)  // eqn.(6) in Erratum
00490    {
00491      a1 = 1.12837;  // Table II in Erratum
00492      a2 = -0.470698;
00493      a3 = -0.453108;
00494      a4 = -228.975;
00495      a5 = 8303.50;
00496      a6 = -118124;
00497      a7 = 632895;
00498      fdInt = sqrt(x)*(a1+ a2/pow(x,2.)+ a3/pow(x,4.)+ a4/pow(x,6.)+ a5/pow(x,8.)+ a6/pow(x,10.)+ a7/pow(x,12.));
00499    }
00500    else if ((x > 0.) && (x <= 2.5))  // eqn.(7) in Erratum
00501    {
00502      double a8, a9;
00503      a1 = 0.604856;  // Table III in Erratum
00504      a2 = 0.380080;
00505      a3 = 0.059320;
00506      a4 = -0.014526;
00507      a5 = -0.004222;
00508      a6 = 0.001335;
00509      a7 = 0.000291;
00510      a8 = -0.000159;
00511      a9 = 0.000018;
00512      fdInt = a1+ a2*x+ a3*pow(x,2.)+ a4*pow(x,3.)+ a5*pow(x,4.)+ a6*pow(x,5.)+ a7*pow(x,6.)+ a8*pow(x,7.)+ a9*pow(x,8.);
00513    }
00514    else  // 2.5<x<5, eqn.(7) in Erratum
00515    {
00516      double a8;
00517      a1 = 0.638086;  // Table III in Erratum
00518      a2 = 0.292266;
00519      a3 = 0.159486;
00520      a4 = -0.077691;
00521      a5 = 0.018650;
00522      a6 = -0.002736;
00523      a7 = 0.000249;
00524      a8 = -0.000013;
00525      fdInt = a1+ a2*x+ a3*pow(x,2.)+ a4*pow(x,3.)+ a5*pow(x,4.)+ a6*pow(x,5.)+ a7*pow(x,6.)+ a8*pow(x,7.);
00526    }
00527    
00528    return fdInt;
00529 }
00530 
00531 
00532 
00533 double compute_dFDIntMinusOneHalf(const double x)
00534 {
00535    // Use the approximate -1/2 FD integral by P. Van Halen and D. L. Pulfrey, 
00536    // "Accurate, short series approximations to Fermi-Dirac integrals of order 
00537    // -/2, 1/2, 1, 3/2, 2, 5/2, 3, and 7/2" J. Appl. Phys. 57, 5271 (1985) 
00538    // and its Erratum in J. Appl. Phys. 59, 2264 (1986). The approximation
00539    // has error < 1e-5 in the entire x range.  
00540    
00541    double dfdInt; 
00542    double u,v,du,dv;
00543    double a1, a2, a3, a4, a5, a6, a7; 
00544    if (x <= 0.)  // eqn.(4) in the reference
00545    {
00546      a1 = 0.999909;  // Table I in Erratum
00547      a2 = 0.706781;
00548      a3 = 0.572752;
00549      a4 = 0.466318;
00550      a5 = 0.324511;
00551      a6 = 0.152889;
00552      a7 = 0.033673;
00553      dfdInt = a1*exp(x)-2*a2*exp(2.*x)+3*a3*exp(3.*x)-4*a4*exp(4.*x)+5*a5*exp(5.*x)-6*a6*exp(6.*x)+7*a7*exp(7.*x);
00554    }
00555    else if (x >= 5.)  // eqn.(6) in Erratum
00556    {
00557      a1 = 1.12837;  // Table II in Erratum
00558      a2 = -0.470698;
00559      a3 = -0.453108;
00560      a4 = -228.975;
00561      a5 = 8303.50;
00562      a6 = -118124;
00563      a7 = 632895;
00564 
00565      u = sqrt(x);
00566      du = 0.5/sqrt(x);
00567      v =  a1 + a2/pow(x,2.)+ a3/pow(x,4.)+ a4/pow(x,6.)+ a5/pow(x,8.)+ a6/pow(x,10.)+ a7/pow(x,12.);
00568      dv = -2*a2/pow(x,3.) -4*a3/pow(x,5.) -6*a4/pow(x,7.) -8*a5/pow(x,9.) -10*a6/pow(x,11) -12*a7/pow(x,13);
00569      dfdInt = u*dv + du*v;
00570    }
00571    else if ((x > 0.) && (x <= 2.5))  // eqn.(7) in Erratum
00572    {
00573      double a8, a9;
00574      a1 = 0.604856;  // Table III in Erratum
00575      a2 = 0.380080;
00576      a3 = 0.059320;
00577      a4 = -0.014526;
00578      a5 = -0.004222;
00579      a6 = 0.001335;
00580      a7 = 0.000291;
00581      a8 = -0.000159;
00582      a9 = 0.000018;
00583      dfdInt = a2+ 2*a3*x + 3*a4*pow(x,2.)+ 4*a5*pow(x,3.)+ 5*a6*pow(x,4.)+ 6*a7*pow(x,5.)+ 7*a8*pow(x,6.)+ 8*a9*pow(x,7.);
00584    }
00585    else  // 2.5<x<5, eqn.(7) in Erratum
00586    {
00587      double a8;
00588      a1 = 0.638086;  // Table III in Erratum
00589      a2 = 0.292266;
00590      a3 = 0.159486;
00591      a4 = -0.077691;
00592      a5 = 0.018650;
00593      a6 = -0.002736;
00594      a7 = 0.000249;
00595      a8 = -0.000013;
00596      dfdInt = a2 + 2*a3*x + 3*a4*pow(x,2.)+ 4*a5*pow(x,3.)+ 5*a6*pow(x,4.)+ 6*a7*pow(x,5.)+ 7*a8*pow(x,6.);
00597    }
00598    
00599    return dfdInt;
00600 }

Generated on Wed Mar 26 2014 18:36:43 for Albany: a Trilinos-based PDE code by  doxygen 1.7.1