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

QCAD_SchrodingerProblem.hpp

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 #ifndef QCAD_SCHRODINGERPROBLEM_HPP
00008 #define QCAD_SCHRODINGERPROBLEM_HPP
00009 
00010 #include "Teuchos_RCP.hpp"
00011 #include "Teuchos_ParameterList.hpp"
00012 
00013 #include "Albany_AbstractProblem.hpp"
00014 #include "Albany_ProblemUtils.hpp"
00015 #include "Albany_ResponseUtilities.hpp"
00016 
00017 #include "Phalanx.hpp"
00018 #include "PHAL_Workset.hpp"
00019 #include "PHAL_Dimension.hpp"
00020 
00021 namespace QCAD {
00022 
00027   class SchrodingerProblem : public Albany::AbstractProblem {
00028   public:
00029   
00031     SchrodingerProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00032            const Teuchos::RCP<ParamLib>& paramLib,
00033            const int numDim_,
00034            const Teuchos::RCP<const Epetra_Comm>& comm_);
00035 
00037     ~SchrodingerProblem();
00038 
00040     virtual int spatialDimension() const { return numDim; }
00041 
00043     virtual void buildProblem(
00044       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00045       Albany::StateManager& stateMgr);
00046 
00047     // Build evaluators
00048     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00049     buildEvaluators(
00050       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00051       const Albany::MeshSpecsStruct& meshSpecs,
00052       Albany::StateManager& stateMgr,
00053       Albany::FieldManagerChoice fmchoice,
00054       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00055 
00057     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00058 
00059   private:
00060 
00062     SchrodingerProblem(const SchrodingerProblem&);
00063     
00065     SchrodingerProblem& operator=(const SchrodingerProblem&);
00066 
00067   public:
00068 
00070     template <typename EvalT>
00071     Teuchos::RCP<const PHX::FieldTag>
00072     constructEvaluators(
00073       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00074       const Albany::MeshSpecsStruct& meshSpecs,
00075       Albany::StateManager& stateMgr,
00076       Albany::FieldManagerChoice fmchoice,
00077       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00078 
00079     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00080 
00081   protected:
00082     Teuchos::RCP<const Epetra_Comm> comm;
00083     bool havePotential;
00084     double energy_unit_in_eV, length_unit_in_m;
00085     std::string potentialFieldName;
00086     int potentialAuxIndex;
00087     std::string mtrlDbFilename;
00088 
00089     int numDim;
00090     int nEigenvectorsToOuputAsStates;
00091     bool bOnlySolveInQuantumBlocks;
00092   };
00093 
00094 }
00095 
00096 #include "QCAD_MaterialDatabase.hpp"
00097 
00098 #include "Intrepid_FieldContainer.hpp"
00099 #include "Intrepid_DefaultCubatureFactory.hpp"
00100 #include "Shards_CellTopology.hpp"
00101 
00102 #include "Albany_Utils.hpp"
00103 #include "Albany_EvaluatorUtils.hpp"
00104 
00105 #include "PHAL_DOFInterpolation.hpp"
00106 
00107 #include "QCAD_SchrodingerPotential.hpp"
00108 #include "QCAD_SchrodingerResid.hpp"
00109 #include "QCAD_ResponseSaddleValue.hpp"
00110 
00111 
00112 template <typename EvalT>
00113 Teuchos::RCP<const PHX::FieldTag> 
00114 QCAD::SchrodingerProblem::constructEvaluators(
00115   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00116   const Albany::MeshSpecsStruct& meshSpecs,
00117   Albany::StateManager& stateMgr,
00118   Albany::FieldManagerChoice fieldManagerChoice,
00119   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00120 {
00121    using Teuchos::RCP;
00122    using Teuchos::rcp;
00123    using Teuchos::ParameterList;
00124    using PHX::DataLayout;
00125    using std::vector;
00126    using std::string;
00127    using PHAL::AlbanyTraits;
00128 
00129    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00130    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00131      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00132 
00133    const int numNodes = intrepidBasis->getCardinality();
00134    const int worksetSize = meshSpecs.worksetSize;
00135 
00136    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00137    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00138 
00139    const int numQPts = cubature->getNumPoints();
00140    const int numVertices = cellType->getNodeCount();
00141 
00142    *out << "Field Dimensions: Workset=" << worksetSize 
00143         << ", Vertices= " << numVertices
00144         << ", Nodes= " << numNodes
00145         << ", QuadPts= " << numQPts
00146         << ", Dim= " << numDim << std::endl;
00147 
00148    RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00149    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00150    bool supportsTransient=true;
00151 
00152    // Temporary variable used numerous times below
00153    Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00154 
00155    // Define Field Names
00156 
00157    Teuchos::ArrayRCP<string> dof_names(neq);
00158      dof_names[0] = "psi";
00159 
00160    Teuchos::ArrayRCP<string> dof_names_dot(neq);
00161    if (supportsTransient) {
00162      for (unsigned int i=0; i<neq; i++) dof_names_dot[i] = dof_names[i]+"_dot";
00163    }
00164 
00165    Teuchos::ArrayRCP<string> resid_names(neq);
00166      for (unsigned int i=0; i<neq; i++) resid_names[i] = dof_names[i]+" Residual";
00167 
00168    if (supportsTransient) fm0.template registerEvaluator<EvalT>
00169        (evalUtils.constructGatherSolutionEvaluator(false, dof_names, dof_names_dot));
00170    else fm0.template registerEvaluator<EvalT>
00171        (evalUtils.constructGatherSolutionEvaluator_noTransient(false, dof_names));
00172 
00173    fm0.template registerEvaluator<EvalT>
00174      (evalUtils.constructScatterResidualEvaluator(false, resid_names));
00175 
00176    fm0.template registerEvaluator<EvalT>
00177      (evalUtils.constructGatherCoordinateVectorEvaluator());
00178 
00179    fm0.template registerEvaluator<EvalT>
00180      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00181 
00182    fm0.template registerEvaluator<EvalT>
00183      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00184 
00185    for (unsigned int i=0; i<neq; i++) {
00186      fm0.template registerEvaluator<EvalT>
00187        (evalUtils.constructDOFInterpolationEvaluator(dof_names[i]));
00188 
00189      if (supportsTransient)
00190      fm0.template registerEvaluator<EvalT>
00191          (evalUtils.constructDOFInterpolationEvaluator(dof_names_dot[i]));
00192 
00193      fm0.template registerEvaluator<EvalT>
00194        (evalUtils.constructDOFGradInterpolationEvaluator(dof_names[i]));
00195   }
00196 
00197    // Create Material Database
00198    RCP<QCAD::MaterialDatabase> materialDB = rcp(new QCAD::MaterialDatabase(mtrlDbFilename, comm));
00199 
00200   if (havePotential) { // If a "Potential" sublist is specified in the input, add a potential energy term
00201 
00202     if(potentialAuxIndex < 0) { 
00203 
00204       // Case when potential is given using the "Potential" parameter sublist
00205 
00206       RCP<ParameterList> p = rcp(new ParameterList);
00207 
00208       p->set<string>("QP Variable Name", "psi");
00209       p->set<string>("QP Potential Name", potentialFieldName);
00210       p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00211 
00212       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00213       Teuchos::ParameterList& paramList = params->sublist("Potential");
00214       p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00215 
00216       //Global Problem Parameters
00217       p->set<double>("Energy unit in eV", energy_unit_in_eV);
00218       p->set<double>("Length unit in m", length_unit_in_m);
00219 
00220       ev = rcp(new QCAD::SchrodingerPotential<EvalT,AlbanyTraits>(*p,dl));
00221       fm0.template registerEvaluator<EvalT>(ev);
00222     }
00223     else {
00224       
00225       //Case when we load the potential from an aux data vector (on the nodes)
00226       // to the potential field (on the quad points).  Note this requires the
00227       // "Type" parameter of the "Potential" input file sublist to be "From Aux Data Vector"
00228       // directives found in the "Potential" input file sublist.  This is used
00229       // when coupling the Poisson and Schrodinger problems (see QCAD::CoupledPoissonSchrodinger)
00230 
00231       Teuchos::ParameterList& paramList = params->sublist("Potential");
00232       std::string potentialType = paramList.get("Type", "not-given");
00233       TEUCHOS_TEST_FOR_EXCEPTION (potentialType != "From Aux Data Vector", Teuchos::Exceptions::InvalidParameter, std::endl 
00234     << "Error! Schrodinger potential type must be \"From Aux Data Vector\" when an aux vector index is specified!" << std::endl);
00235 
00236       // Gather the aux data vector into potentialFieldName
00237       RCP<ParameterList> p = rcp(new ParameterList);
00238       p->set<string>("Field Name", potentialFieldName);
00239       p->set<int>("Aux Data Vector Index", potentialAuxIndex);
00240 
00241       ev = rcp(new PHAL::GatherAuxData<EvalT,AlbanyTraits>(*p,dl));
00242       fm0.template registerEvaluator<EvalT>(ev);
00243 
00244       // Interpolate potential to quad points (use DOFInterpolation)
00245       p = rcp(new ParameterList("Interpolate potential to quad points"));
00246       p->set<string>("Variable Name", potentialFieldName); // assumes same Name for output as for input 
00247       p->set<string>("BF Name", "BF");
00248       p->set<int>("Offset of First DOF", 0);
00249 
00250       ev = rcp(new PHAL::DOFInterpolation<EvalT,AlbanyTraits>(*p,dl));
00251       fm0.template registerEvaluator<EvalT>(ev);
00252     }
00253   }
00254 
00255   { // Wavefunction (psi) Resid
00256     RCP<ParameterList> p = rcp(new ParameterList("Wavefunction Resid"));
00257 
00258     //Input
00259     p->set<string>("Weighted BF Name", "wBF");
00260     p->set<string>("QP Variable Name", "psi");
00261     p->set<string>("QP Time Derivative Variable Name", "psi_dot");
00262     p->set<string>("QP Coordinate Vector Name", "Coord Vec");
00263 
00264     p->set<bool>("Have Potential", havePotential);
00265     p->set<string>("Potential Name", potentialFieldName); // was "V"
00266 
00267     p->set<string>("Gradient QP Variable Name", "psi Gradient");
00268 
00269     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00270 
00271     //Output
00272     p->set<string>("Residual Name", "psi Residual");
00273 
00274     //Global Problem Parameters
00275     p->set<double>("Energy unit in eV", energy_unit_in_eV);
00276     p->set<double>("Length unit in m", length_unit_in_m);
00277     p->set<bool>("Only solve in quantum blocks", bOnlySolveInQuantumBlocks);
00278     p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00279 
00280     //Pass the Potential parameter list to test Finite Wall with different effective mass
00281     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00282     Teuchos::ParameterList& paramList = params->sublist("Potential");
00283     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00284 
00285     ev = rcp(new QCAD::SchrodingerResid<EvalT,AlbanyTraits>(*p,dl));
00286     fm0.template registerEvaluator<EvalT>(ev);
00287   }
00288 
00289 
00290   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00291     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00292     fm0.requireField<EvalT>(res_tag);
00293     return res_tag.clone();
00294   }
00295 
00296   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00297 
00298     // Parameters to be sent to all response constructors (whether they use them or not).
00299     RCP<ParameterList> pFromProb = rcp(new ParameterList("Response Parameters from Problem"));
00300     pFromProb->set<double>("Length unit in m", length_unit_in_m);
00301     pFromProb->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00302 
00303     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00304     return respUtils.constructResponses(fm0, *responseList, pFromProb, stateMgr);
00305   }
00306 
00307   return Teuchos::null;
00308 }
00309 #endif

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