Albany: a Trilinos-based PDE code: QCAD

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_SOLVER_H
00008 #define QCAD_SOLVER_H
00009 
00010 #include <iostream>
00011 
00012 #include "LOCA.H"
00013 #include "LOCA_Epetra.H"
00014 #include "Epetra_Vector.h"
00015 #include "Epetra_LocalMap.h"
00016 #include "LOCA_Epetra_ModelEvaluatorInterface.H"
00017 #include <NOX_Epetra_MultiVector.H>
00018 
00019 #include "Albany_ModelEvaluator.hpp"
00020 #include "Albany_Utils.hpp"
00021 #include "Piro_Epetra_StokhosNOXObserver.hpp"
00022 
00023 #include "QCAD_MultiSolutionObserver.hpp"
00024 
00025 #ifdef ALBANY_CI
00026 #include "AlbanyCI_Types.hpp"
00027 #include "AlbanyCI_Tensor.hpp"
00028 #include "AlbanyCI_BlockTensor.hpp"
00029 #include "AlbanyCI_SingleParticleBasis.hpp"
00030 #include "AlbanyCI_BasisFactory.hpp"
00031 #include "AlbanyCI_ManyParticleBasis.hpp"
00032 #include "AlbanyCI_ManyParticleBasisBlock.hpp"
00033 #include "AlbanyCI_MatrixFactory.hpp"
00034 #include "AlbanyCI_ManyParticleMatrix.hpp"
00035 #include "AlbanyCI_Solver.hpp"
00036 #include "AlbanyCI_Solution.hpp"
00037 #include "AlbanyCI_qnumbers.hpp"
00038 #endif
00039 
00040 
00041 
00042 namespace QCAD {
00043   class SolverParamFn;
00044   class SolverResponseFn;
00045   class SolverSubSolver;
00046   class SolverSubSolverData;
00047 
00052   class Solver : public EpetraExt::ModelEvaluator {
00053   public:
00054 
00057 
00058       Solver(const Teuchos::RCP<Teuchos::ParameterList>& appParams,
00059        const Teuchos::RCP<const Epetra_Comm>& comm,
00060        const Teuchos::RCP<const Epetra_Vector>& initial_guess);
00062 
00063     ~Solver();
00064 
00065     Teuchos::RCP<const Epetra_Map> get_x_map() const;
00066     Teuchos::RCP<const Epetra_Map> get_f_map() const;
00067     Teuchos::RCP<const Epetra_Map> get_p_map(int l) const;
00068     Teuchos::RCP<const Epetra_Map> get_g_map(int j) const;
00069 
00070     Teuchos::RCP<const Epetra_Vector> get_x_init() const;
00071     Teuchos::RCP<const Epetra_Vector> get_p_init(int l) const;
00072 
00073     EpetraExt::ModelEvaluator::InArgs createInArgs() const;
00074     EpetraExt::ModelEvaluator::OutArgs createOutArgs() const;
00075 
00076     void evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const;    
00077 
00078 
00079   private:
00080     Teuchos::RCP<Teuchos::ParameterList> createPoissonInputFile(const Teuchos::RCP<Teuchos::ParameterList>& appParams,
00081                 int numDims, int nEigen, const std::string& specialProcessing,
00082                 const std::string& xmlOutputFile, const std::string& exoOutputFile) const;
00083     Teuchos::RCP<Teuchos::ParameterList> createSchrodingerInputFile(const Teuchos::RCP<Teuchos::ParameterList>& appParams,
00084                     int numDims, int nEigen, const std::string& specialProcessing,
00085                     const std::string& xmlOutputFile, const std::string& exoOutputFile) const;
00086     Teuchos::RCP<Teuchos::ParameterList> createPoissonSchrodingerInputFile(const Teuchos::RCP<Teuchos::ParameterList>& appParams,
00087                      int numDims, int nEigen, const std::string& xmlOutputFile,
00088                      const std::string& exoOutputFile) const;
00089 
00090     void evalPoissonSchrodingerModel(const InArgs& inArgs, const OutArgs& outArgs,
00091              std::vector<double>& eigenvalResponses, std::map<std::string, SolverSubSolver>& subSolvers) const;
00092     void evalPoissonCIModel(const InArgs& inArgs, const OutArgs& outArgs,
00093           std::vector<double>& eigenvalResponses, std::map<std::string, SolverSubSolver>& subSolvers) const;
00094     void evalCIModel(const InArgs& inArgs, const OutArgs& outArgs, 
00095          std::vector<double>& eigenvalResponses, std::map<std::string, SolverSubSolver>& subSolvers) const;
00096 
00097     bool doPSLoop(const std::string& mode, const InArgs& inArgs, std::map<std::string, SolverSubSolver>& subSolvers, 
00098       Teuchos::RCP<Albany::EigendataStruct>& eigenDataResult, bool bPrintNumOfQuantumElectrons) const;
00099 
00100     void setupParameterMapping(const Teuchos::ParameterList& list, const std::string& defaultSubSolver,
00101              const std::map<std::string, SolverSubSolverData>& subSolversData);
00102     void setupResponseMapping(const Teuchos::ParameterList& list, const std::string& defaultSubSolver, int nEigenvalues,
00103             const std::map<std::string, SolverSubSolverData>& subSolversData);
00104 
00105     void fillSingleSubSolverParams(const InArgs& inArgs, const std::string& name, 
00106            QCAD::SolverSubSolver& subSolver, int nLeaveOffEnd=0) const;
00107 
00108     SolverSubSolver CreateSubSolver(const Teuchos::RCP<Teuchos::ParameterList> appParams, const Epetra_Comm& comm,
00109             const Teuchos::RCP<const Epetra_Vector>& initial_guess  = Teuchos::null) const;
00110 
00111     SolverSubSolverData CreateSubSolverData(const QCAD::SolverSubSolver& sub) const;
00112 
00113 
00114     const Teuchos::RCP<Teuchos::ParameterList>& getSubSolverParams(const std::string& name) const;
00115     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00116 
00117     void printResponses(const QCAD::SolverSubSolver& solver, 
00118       const std::string& solverName, 
00119       Teuchos::RCP<Teuchos::FancyOStream> out) const;
00120     
00121   private:
00122     int numDims;
00123     std::string problemNameBase;
00124     std::string defaultSubSolver;
00125     Teuchos::RCP<Teuchos::ParameterList> mainAppParams;
00126     std::map<std::string, Teuchos::RCP<Teuchos::ParameterList> > subProblemAppParams;
00127 
00128     std::vector< std::vector<Teuchos::RCP<SolverParamFn> > > paramFnVecs;
00129     std::vector<Teuchos::RCP<SolverResponseFn> > responseFns;
00130 
00131     std::size_t maxIter;
00132     std::size_t nParameters;
00133     std::size_t nResponseDoubles;    
00134 
00135     std::string iterationMethod;
00136     int  nEigenvectors; //used in Poisson-CI coupling
00137 
00138     int num_p, num_g;
00139     Teuchos::RCP<Epetra_LocalMap> epetra_param_map;
00140     Teuchos::RCP<Epetra_LocalMap> epetra_response_map;
00141     Teuchos::RCP<Epetra_Map> epetra_x_map;
00142 
00143     Teuchos::RCP<Epetra_Vector> epetra_param_vec;
00144     DerivativeSupport deriv_support;
00145 
00146     Teuchos::RCP<const Epetra_Comm> solverComm;
00147     Teuchos::RCP<const Epetra_Vector> saved_initial_guess;
00148 
00149     bool bVerbose;
00150     bool bSupportDpDg;
00151     bool bRealEvecs;
00152 
00153     std::string eigensolverName;
00154     double ps_converge_tol;
00155     double shiftPercentBelowMin;  // for eigensolver shift-invert: shift point == minPotential * (1 + shiftPercent/100)
00156     int    minCIParticles;        // the minimum number of particles allowed to be used in CI calculation
00157     int    maxCIParticles;        // the maximum number of particles allowed to be used in CI calculation
00158     int    nCIParticles;          // the number of particles used in CI calculation
00159     int    nCIExcitations;        // the number of excitations used in CI calculation
00160     double fixedPSOcc;
00161     bool   bUseIntegratedPS;
00162     bool   bUseTotalSpinSymmetry; // use S2 symmetry in CI calculation
00163   };
00164 
00165 
00166   // helper classes - maybe nest inside Solver?
00167   class SolverParamFn {
00168   public:
00169     SolverParamFn(const std::string& fnString, 
00170       const std::map<std::string, SolverSubSolverData>& subSolversData);
00171     ~SolverParamFn() {};
00172 
00173     void fillSingleSubSolverParams(double parameterValue, const std::string& subSolverName,
00174            SolverSubSolver& subSolver) const;
00175 
00176     void fillSubSolverParams(double parameterValue, 
00177            const std::map<std::string, SolverSubSolver>& subSolvers) const;
00178 
00179     double getInitialParam(const std::map<std::string, SolverSubSolverData>& subSolversData) const;
00180 
00181     std::string getTargetName() const { return targetName; }
00182     std::vector<int> getTargetIndices() const { return targetIndices; }
00183     std::size_t getNumFilters() const { return filters.size(); }
00184     double getFilterScaling() const;
00185 
00186   protected:
00187     std::string targetName;
00188     std::vector<int> targetIndices;
00189     std::vector< std::vector<std::string> > filters;
00190   };
00191 
00192   class SolverResponseFn {
00193   public:
00194     SolverResponseFn(const std::string& fnString,
00195          const std::map<std::string, SolverSubSolverData>& subSolversData,
00196          int nEigenvalues);
00197     ~SolverResponseFn() {};
00198 
00199     void fillSolverResponses(Epetra_Vector& g, Teuchos::RCP<Epetra_MultiVector>& dgdp, int offset,
00200            const std::map<std::string, SolverSubSolver>& subSolvers,
00201            const std::vector<std::vector<Teuchos::RCP<SolverParamFn> > >& paramFnVecs,
00202            bool bSupportDpDg, const std::vector<double>& eigenvalueResponses) const;
00203 
00204     std::size_t getNumDoubles() const { return numDoubles; }
00205 
00206   protected:
00207     struct ArrayRef { 
00208       std::string name; 
00209       std::vector<int> indices;
00210     };
00211 
00212     std::string fnName;
00213     std::vector<ArrayRef> params;
00214     std::size_t numDoubles; //number of doubles produced by this response
00215   };
00216 
00217 
00218   class SolverSubSolver {
00219   public:
00220     Teuchos::RCP<Albany::Application> app;
00221     Teuchos::RCP<EpetraExt::ModelEvaluator> model;
00222     Teuchos::RCP<EpetraExt::ModelEvaluator::InArgs> params_in;
00223     Teuchos::RCP<EpetraExt::ModelEvaluator::OutArgs> responses_out;
00224     void freeUp() { app = Teuchos::null; model = Teuchos::null; }
00225   };
00226 
00227   class SolverSubSolverData {
00228   public:
00229     int Np, Ng;
00230     std::vector<int> pLength, gLength;
00231     Teuchos::RCP<const Epetra_Vector> p_init;
00232     EpetraExt::ModelEvaluator::DerivativeSupport deriv_support;
00233   };
00234 
00235 
00236 #ifdef ALBANY_CI
00237   class CISolver {
00238   public:
00239     CISolver(int n1PSpinlessStates, Teuchos::RCP<const Epetra_Comm> eComm, 
00240        Teuchos::RCP<Teuchos::FancyOStream> outStream);
00241 
00242     Teuchos::RCP<Teuchos::ParameterList> getDefaultParameterList() const;
00243 
00244     void fill1Pmx(const Teuchos::RCP<Albany::EigendataStruct>& eigenData1P);
00245     void fill1Pmx(const Teuchos::RCP<Albany::EigendataStruct>& eigenData1P,
00246       const Teuchos::RCP<Epetra_Vector>& g_noCharge,
00247       const Teuchos::RCP<Epetra_Vector>& g_delta,
00248       double deltaScale, bool bRealEvecs, bool bVerbose);
00249     void fill2Pmx(Teuchos::RCP<Albany::EigendataStruct> eigenData1P,
00250       const SolverSubSolver* coulombSolver, 
00251       const SolverSubSolver* coulombSolver_ImPart,
00252       const Teuchos::RCP<Epetra_Vector>& g_noCharge,
00253       bool bRealEvecs, bool bVerbose);
00254 
00255     Teuchos::RCP<AlbanyCI::Solution> Solve(Teuchos::RCP<Teuchos::ParameterList> AlbanyCIList) const;
00256 
00257     Teuchos::RCP<Epetra_MultiVector> ComputeStateDensities(Teuchos::RCP<Albany::EigendataStruct> eigenData1P,
00258                  Teuchos::RCP<AlbanyCI::Solution> soln);
00259 
00260 
00261   private:
00262     void SetCoulombParams(const Teuchos::RCP<EpetraExt::ModelEvaluator::InArgs> inArgs, int i2, int i4) const;
00263 
00264   private:
00265     // number of single particle states of each type of spin (up / down)
00266     int n1PperBlock;
00267 
00268     // 1P Blocks, accessed individually or as a vector
00269     Teuchos::RCP<AlbanyCI::Tensor<AlbanyCI::dcmplx> > blockU,blockD;
00270     std::vector<Teuchos::RCP<AlbanyCI::Tensor<AlbanyCI::dcmplx> > > blocks1P;
00271     Teuchos::RCP<AlbanyCI::BlockTensor<AlbanyCI::dcmplx> > mx1P;
00272 
00273     // 2P Blocks, accessed individually or as a vector
00274     Teuchos::RCP<AlbanyCI::Tensor<AlbanyCI::dcmplx> > blockUU, blockUD, blockDU, blockDD;
00275     std::vector<Teuchos::RCP<AlbanyCI::Tensor<AlbanyCI::dcmplx> > > blocks2P;
00276     Teuchos::RCP<AlbanyCI::BlockTensor<AlbanyCI::dcmplx> > mx2P;
00277 
00278     // 1P basis
00279     Teuchos::RCP<AlbanyCI::SingleParticleBasis> basis1P;
00280 
00281     // MPI Comm
00282     Teuchos::RCP<Teuchos::Comm<int> > comm;    
00283 
00284     // Output stream
00285     Teuchos::RCP<Teuchos::FancyOStream> out;
00286   };
00287 #endif
00288   
00289 }
00290 #endif
QCAD_Solver.hpp
