Albany: a Trilinos-based PDE code: QCAD_SchrodingerPotential

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 <fstream>
00008 #include "Teuchos_TestForException.hpp"
00009 #include "Phalanx_DataLayout.hpp"
00010 #include "Sacado_ParameterRegistration.hpp"
00011 #include "QCAD_StringFormulaEvaluator.hpp"
00012 
00013 template<typename EvalT, typename Traits>
00014 QCAD::SchrodingerPotential<EvalT, Traits>::
00015 SchrodingerPotential(Teuchos::ParameterList& p,
00016                  const Teuchos::RCP<Albany::Layouts>& dl) :
00017   coordVec(p.get<std::string>("QP Coordinate Vector Name"), dl->qp_gradient),
00018   psi(p.get<std::string>("QP Variable Name"), dl->qp_scalar),
00019   V(p.get<std::string>("QP Potential Name"), dl->qp_scalar)
00020 {
00021   // psList == "Potential" sublist of main input file
00022   Teuchos::ParameterList* psList = p.get<Teuchos::ParameterList*>("Parameter List");
00023 
00024   Teuchos::RCP<const Teuchos::ParameterList> reflist = 
00025       this->getValidSchrodingerPotentialParameters();
00026   psList->validateParameters(*reflist,0);
00027 
00028   std::vector<PHX::DataLayout::size_type> dims;
00029   dl->qp_gradient->dimensions(dims);
00030   numQPs  = dims[1];
00031   numDims = dims[2];
00032 
00033   energy_unit_in_eV = p.get<double>("Energy unit in eV");
00034   length_unit_in_m = p.get<double>("Length unit in m");
00035   
00036   potentialType = psList->get("Type", "defaultType");
00037   E0 = psList->get("E0", 1.0);
00038   scalingFactor = psList->get("Scaling Factor", 1.0);
00039   
00040   // Parameters for String Formula
00041   stringFormula = psList->get("Formula", "0");
00042 
00043   // Parameters for Finite Wall 
00044   barrEffMass = psList->get<double>("Barrier Effective Mass", 0.0);
00045   barrWidth = psList->get<double>("Barrier Width", 0.0);
00046   wellEffMass = psList->get<double>("Well Effective Mass", 0.0);
00047   wellWidth = psList->get<double>("Well Width", 0.0);
00048   
00049   potentialStateName = psList->get<std::string>("State Name","Not Specified");
00050 
00051   // Add E0 as a Sacado-ized parameter
00052   Teuchos::RCP<ParamLib> paramLib =
00053       p.get< Teuchos::RCP<ParamLib> >("Parameter Library", Teuchos::null);
00054   new Sacado::ParameterRegistration<EvalT, SPL_Traits>(
00055       "Schrodinger Potential E0", this, paramLib);
00056   new Sacado::ParameterRegistration<EvalT, SPL_Traits> (
00057       "Schrodinger Potential Scaling Factor", this, paramLib);
00058 
00059   this->addDependentField(psi);
00060   this->addDependentField(coordVec);
00061 
00062   this->addEvaluatedField(V);
00063   this->setName("Schrodinger Potential"+PHX::TypeString<EvalT>::value);
00064 }
00065 
00066 // **********************************************************************
00067 template<typename EvalT, typename Traits>
00068 void QCAD::SchrodingerPotential<EvalT, Traits>::
00069 postRegistrationSetup(typename Traits::SetupData d,
00070                       PHX::FieldManager<Traits>& fm)
00071 {
00072   this->utils.setFieldData(V,fm);
00073   this->utils.setFieldData(psi,fm);
00074   this->utils.setFieldData(coordVec,fm);
00075 }
00076 
00077 // **********************************************************************
00078 template<typename EvalT, typename Traits>
00079 void QCAD::SchrodingerPotential<EvalT, Traits>::
00080 evaluateFields(typename Traits::EvalData workset)
00081 {
00082 
00083   // Parabolic potential (test case)
00084   if (potentialType == "Parabolic") 
00085   {
00086     for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
00087       for (std::size_t qp = 0; qp < numQPs; ++qp) {
00088         V(cell, qp) = parabolicPotentialValue(numDims, &coordVec(cell,qp,0));
00089       }
00090     }
00091   }
00092   
00093   // Infinite barrier wall
00094   else if (potentialType == "Infinite Wall")
00095   {
00096     for (std::size_t cell = 0; cell < workset.numCells; ++cell)
00097       for (std::size_t qp = 0; qp < numQPs; ++qp)
00098         V(cell, qp) = 0.0; 
00099   } 
00100 
00101   // Finite barrier wall
00102   else if (potentialType == "Finite Wall")
00103   {
00104     for (std::size_t cell = 0; cell < workset.numCells; ++cell)
00105     {
00106       for (std::size_t qp = 0; qp < numQPs; ++qp)
00107         V(cell, qp) = finiteWallPotential(numDims, &coordVec(cell,qp,0));
00108     }    
00109   }
00110 
00111   // Specifed by formula string
00112   else if (potentialType == "String Formula")
00113   {
00114     for (std::size_t cell = 0; cell < workset.numCells; ++cell)
00115     {
00116       for (std::size_t qp = 0; qp < numQPs; ++qp)
00117         V(cell, qp) = stringFormulaPotential(numDims, &coordVec(cell,qp,0));
00118     }    
00119   }
00120 
00121   
00122   // Potential energy taken from a StateManager state
00123   else if (potentialType == "From State") 
00124   {
00125     Albany::StateArray& states = *workset.stateArrayPtr;
00126     Albany::MDArray& potentialState = states[potentialStateName];
00127 
00128     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00129       for (std::size_t qp=0; qp < numQPs; ++qp) {
00130         double d =  potentialState(cell, qp);
00131         V(cell, qp) = d; //scalingFactor * d;
00132 
00133   //ANDY: behavior I don't understand - scalingFactor gets set to 1000?
00134   //if( scalingFactor != 1.0 ) 
00135   //  std::cout << "DEBUG: scaling factor = " << scalingFactor << std::endl;
00136   //HACK to help anasazi solve
00137   //if(workset.EBName == "silicon" || scalingFactor < 0) {
00138   //  V(cell, qp) =  d;
00139   //}
00140   //else V(cell, qp) = scalingFactor;
00141       }
00142     }
00143   }    
00144 
00145   /***** otherwise error? ******/
00146   else 
00147   {
00148     TEUCHOS_TEST_FOR_EXCEPT(true);
00149   }
00150 }
00151 
00152 // **********************************************************************
00153 template<typename EvalT,typename Traits>
00154 typename QCAD::SchrodingerPotential<EvalT,Traits>::ScalarT& 
00155 QCAD::SchrodingerPotential<EvalT,Traits>::getValue(const std::string &n)
00156 {
00157   if(n == "Schrodinger Potential Scaling Factor") return scalingFactor;
00158   else if(n == "Schrodinger Potential E0") return E0;
00159   else TEUCHOS_TEST_FOR_EXCEPT(true); return E0; //dummy so all control paths return
00160 }
00161 
00162 // **********************************************************************
00163 template<typename EvalT,typename Traits>
00164 Teuchos::RCP<const Teuchos::ParameterList>
00165 QCAD::SchrodingerPotential<EvalT,Traits>::getValidSchrodingerPotentialParameters() const
00166 {
00167   Teuchos::RCP<Teuchos::ParameterList> validPL =
00168        rcp(new Teuchos::ParameterList("Valid Schrodinger Potential Params"));;
00169 
00170   validPL->set<std::string>("Type", "defaultType", "Switch between different potential types");
00171   validPL->set<double>("E0", 1.0, "Energy scale - dependent on type");
00172   validPL->set<double>("Scaling Factor", 1.0, "Constant scaling factor");
00173 
00174   // For string formula potential
00175   validPL->set<std::string>("Formula", "0", "Mathematical expression containing x,y,z that specifies the potential at coordinates (x,y,z).  In 1D and 2D use just x and x,y respectively.");
00176   
00177   // For Finite Wall potential 
00178   validPL->set<double>("Barrier Effective Mass", 0.0, "Barrier effective mass in [m0]");
00179   validPL->set<double>("Barrier Width", 0.0, "Barrier width in length_unit_in_m");
00180   validPL->set<double>("Well Effective Mass", 0.0, "Well effective mass in [m0]");
00181   validPL->set<double>("Well Width", 0.0, "Well width in length_unit_in_m");
00182 
00183   // For 1D MOSCapacitor to test the 1D P-S iteration
00184   /* specific parameters for 1D MOSCapacitor to set correct effective mass 
00185      for oxide and silicon regions (could be given in the Poisson Coupling
00186      section, putting here for less perturbation to the main code). */
00187 
00188   validPL->set<double>("Oxide Width", 0.0, "Oxide width in length_unit_in_m");
00189   validPL->set<double>("Silicon Width", 0.0, "Silicon width in length_unit_in_m");
00190 
00191   // For potential imported from a state
00192   validPL->set<std::string>("State Name", "", "Name of state to import potential from");
00193 
00194   return validPL;
00195 }
00196 
00197 // **********************************************************************
00198 
00199 //Return potential in energy_unit_in_eV * eV units
00200 template<typename EvalT,typename Traits>
00201 typename QCAD::SchrodingerPotential<EvalT,Traits>::ScalarT
00202 QCAD::SchrodingerPotential<EvalT,Traits>::parabolicPotentialValue(
00203     const int numDim, const MeshScalarT* coord)
00204 {
00205   ScalarT val;
00206   MeshScalarT r2;
00207   int i;
00208 
00209   //std::cout << "x = " << coord[0] << endl; //in 1D, x-coords run from zero to 1
00210 
00211   /***** define universal constants as double constants *****/
00212   const double hbar = 1.0546e-34;  // Planck constant [J s]
00213   const double emass = 9.1094e-31; // Electron mass [kg]
00214   const double evPerJ = 6.2415e18; // eV per Joule (J/eV)
00215 
00216   const double parabolicFctr = 0.5 * (emass / (evPerJ * hbar*hbar) );
00217 
00218   // prefactor from constant, including scaling due to units
00219   ScalarT prefactor;  
00220   
00221   prefactor = parabolicFctr * E0*E0 * (energy_unit_in_eV * pow(length_unit_in_m,2));  
00222   for(i=0, r2=0.0; i<numDim; i++)
00223     r2 += (coord[i]-0.5)*(coord[i]-0.5);
00224   val = prefactor * r2;  
00225   
00226   return scalingFactor * val;
00227 }
00228 
00229 // **********************************************************************
00230 
00231 //Return potential in energy_unit_in_eV * eV units
00232 template<typename EvalT,typename Traits>
00233 typename QCAD::SchrodingerPotential<EvalT,Traits>::ScalarT
00234 QCAD::SchrodingerPotential<EvalT,Traits>::finiteWallPotential(
00235     const int numDim, const MeshScalarT* coord)
00236 {
00237   ScalarT val;
00238   
00239   switch (numDim)
00240   {
00241     case 1: // 1D: total width = 2*barrWidth + wellWidth
00242     {
00243       if ( (coord[0] >= 0) && (coord[0] < barrWidth) )
00244         val = E0;  // barrier
00245       else if ( (coord[0] >= barrWidth) && (coord[0] <= (barrWidth+wellWidth)) )
00246         val = 0;   // well
00247       else if ( (coord[0] > (barrWidth+wellWidth)) && (coord[0] <= (2*barrWidth+wellWidth)) )
00248         val = E0;  // barrier
00249       else 
00250         TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00251           << "Error! x coordinate is outside [0, 2*barrWidth+wellWidth] range,"
00252           << " make sure 1D Scale in Discretization equal to 2*barrWidth+wellWidth !" << std::endl);
00253       break;
00254     }
00255     
00256     case 2: // 2D
00257     {
00258       if ( (coord[0] >= barrWidth) && (coord[0] <= (barrWidth+wellWidth)) &&
00259            (coord[1] >= barrWidth) && (coord[1] <= (barrWidth+wellWidth)) )
00260         val = 0.0;  // well
00261       else
00262         val = E0;   // barrier
00263       break;
00264     }
00265 
00266     case 3: // 3D
00267     {
00268       if ( (coord[0] >= barrWidth) && (coord[0] <= (barrWidth+wellWidth)) &&
00269            (coord[1] >= barrWidth) && (coord[1] <= (barrWidth+wellWidth)) && 
00270            (coord[2] >= barrWidth) && (coord[2] <= (barrWidth+wellWidth)) )
00271         val = 0.0;  // well
00272       else
00273         val = E0;   // barrier
00274       break;
00275     }
00276     
00277     default: 
00278     {
00279       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00280         << "Error! Invalid numDim = " << numDim << ", must be 1 or 2 or 3 !" << std::endl);
00281       break;  
00282     }
00283     
00284   }  // end of switch (numDim)
00285   
00286   return scalingFactor * val;
00287 }
00288 
00289 
00290 // **********************************************************************
00291 
00292 //Return potential in energy_unit_in_eV * eV units
00293 template<typename EvalT,typename Traits>
00294 typename QCAD::SchrodingerPotential<EvalT,Traits>::ScalarT
00295 QCAD::SchrodingerPotential<EvalT,Traits>::stringFormulaPotential(
00296     const int numDim, const MeshScalarT* coord)
00297 {
00298   ScalarT val;
00299   MeshScalarT x, y, z, result;
00300   
00301   switch (numDim)
00302   {
00303     case 1: // 1D: total width = 2*barrWidth + wellWidth
00304     {
00305       x = coord[0];
00306       y = 0; z=0;
00307       break;
00308     }
00309     
00310     case 2: // 2D
00311     {
00312       x = coord[0];
00313       y = coord[1];
00314       z = 0;
00315       break;
00316     }
00317 
00318     case 3: // 3D
00319     {
00320       x = coord[0];
00321       y = coord[1];
00322       z = coord[2];
00323       break;
00324     }
00325     
00326     default: 
00327     {
00328       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00329         << "Error! Invalid numDim = " << numDim << ", must be 1 or 2 or 3 !" << std::endl);
00330       break;  
00331     }    
00332   }  // end of switch (numDim)
00333   
00334   result = Evaluate(stringFormula,x,y,z); //Parse and evaluate formula string 
00335   val = result;                                // (could speed this up in future by saving intermediate parse info)
00336   return scalingFactor * val;
00337 }
00338 
00339 
00340 // *****************************************************************************
00341
QCAD_SchrodingerPotential_Def.hpp
