Albany: a Trilinos-based PDE code: QCAD_CoupledPoissonSchrodinger.cpp Source File

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_CoupledPoissonSchrodinger.hpp"
00008 #include "QCAD_CoupledPSJacobian.hpp"
00009 #include "Piro_Epetra_LOCASolver.hpp"
00010 #include "Stokhos.hpp"
00011 #include "Stokhos_Epetra.hpp"
00012 #include "Sacado_PCE_OrthogPoly.hpp"
00013 #include "Teuchos_XMLParameterListHelpers.hpp"
00014 #include "Teuchos_TestForException.hpp"
00015 
00016 #include "Teuchos_ParameterList.hpp"
00017 
00018 #include "Albany_ModelFactory.hpp"
00019 #include "Albany_Utils.hpp"
00020 #include "Albany_SolverFactory.hpp"
00021 #include "Albany_StateInfoStruct.hpp"
00022 #include "Albany_EigendataInfoStruct.hpp"
00023 
00024 #include "QCAD_CoupledPSJacobian.hpp"
00025 #include "QCAD_CoupledPSPreconditioner.hpp"
00026 #include "QCAD_MultiSolutionObserver.hpp" //for utility functions
00027 
00028 //For creating discretiation object without a problem object
00029 #include "Albany_DiscretizationFactory.hpp"
00030 #include "Albany_StateInfoStruct.hpp"
00031 #include "Albany_AbstractFieldContainer.hpp"
00032 #include "Piro_NullSpaceUtils.hpp"
00033 
00034 //Ifpack includes
00035 #include "Ifpack_ConfigDefs.h"
00036 #include "Ifpack.h"
00037 
00038 
00039 std::string QCAD::strdim(const std::string s, const int dim) {
00040   std::ostringstream ss;
00041   ss << s << " " << dim << "D";
00042   return ss.str();
00043 }
00044 
00045 
00046 QCAD::CoupledPoissonSchrodinger::
00047 CoupledPoissonSchrodinger(const Teuchos::RCP<Teuchos::ParameterList>& appParams_,
00048         const Teuchos::RCP<const Epetra_Comm>& comm,
00049         const Teuchos::RCP<const Epetra_Vector>& initial_guess)
00050 {
00051   using std::string;
00052   
00053   // make a copy of the appParams, since we modify them below (e.g. discretization list)
00054   Teuchos::RCP<Teuchos::ParameterList> appParams = Teuchos::rcp( new Teuchos::ParameterList(*appParams_) );
00055 
00056   const Albany_MPI_Comm& mcomm = Albany::getMpiCommFromEpetraComm(*comm);
00057   Teuchos::RCP<Teuchos::Comm<int> > tcomm = Albany::createTeuchosCommFromMpiComm(mcomm);
00058 
00059   // Get sub-problem input xml files from problem parameters
00060   Teuchos::ParameterList& problemParams = appParams->sublist("Problem");
00061 
00062   // Validate Problem parameters against list for this specific problem
00063   problemParams.validateParameters(*getValidProblemParameters(),0);
00064 
00065   // Get the number of dimensions
00066   numDims = 0;
00067   string name = problemParams.get<std::string>("Name");
00068   if(name == "Poisson Schrodinger 1D") numDims = 1;
00069   else if(name == "Poisson Schrodinger 2D") numDims = 2;
00070   else if(name == "Poisson Schrodinger 3D") numDims = 3;
00071   else TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00072         << "Error!  Invalid problem name " << name << std::endl);
00073 
00074   // Get parameters from Problem sublist used to generate poisson and schrodinger app lists
00075   int vizDetail         = problemParams.get<int>("Phalanx Graph Visualization Detail");
00076   double Temp           = problemParams.get<double>("Temperature");
00077   double lenUnit        = problemParams.get<double>("Length Unit In Meters", 1e-6);
00078   double energyUnit     = problemParams.get<double>("Energy Unit In Electron Volts", 1.0);
00079   std::string matrlFile = problemParams.get<std::string>("MaterialDB Filename", "materials.xml");
00080   bool   bXCPot         = problemParams.get<bool>("Include exchange-correlation potential",false);
00081   bool   bQBOnly        = problemParams.get<bool>("Only solve schrodinger in quantum blocks",true);
00082   
00083   nEigenvals   = problemParams.get<int>("Number of Eigenvalues");
00084   Teuchos::ParameterList& discList = appParams->sublist("Discretization");
00085   Teuchos::ParameterList& poisson_subList = problemParams.sublist("Poisson Problem", false);
00086   Teuchos::ParameterList& schro_subList = problemParams.sublist("Schrodinger Problem", false);
00087 
00088   // Process debug options to write poisson and schrodinger app params to files
00089   Teuchos::ParameterList& debugList = appParams->sublist("Debug Output", true);
00090   std::string poissonXmlFile       = debugList.get("Poisson XML Input", "");
00091   std::string schrodingerXmlFile   = debugList.get("Schrodinger XML Input", "");
00092   std::string poissonExoOutput     = debugList.get("Poisson Exodus Output", "");
00093   std::string schrodingerExoOutput = debugList.get("Schrodinger Exodus Output", "");
00094 
00095 
00096   // Create input parameter list for poission app which mimics a separate input file
00097   Teuchos::RCP<Teuchos::ParameterList> poisson_appParams = 
00098     Teuchos::createParameterList("Poisson Application Parameters");
00099   Teuchos::ParameterList& poisson_probParams = poisson_appParams->sublist("Problem",false);
00100   
00101   poisson_probParams.set("Name", QCAD::strdim("Poisson",numDims));
00102   poisson_probParams.set("Phalanx Graph Visualization Detail", vizDetail);
00103   poisson_probParams.set("Length Unit In Meters",lenUnit);
00104   poisson_probParams.set("Energy Unit In Electron Volts",energyUnit);
00105   poisson_probParams.set("Temperature",Temp);
00106   poisson_probParams.set("MaterialDB Filename", matrlFile);
00107 
00108   {
00109     Teuchos::ParameterList auto_sourceList;
00110     auto_sourceList.set("Factor",1.0);
00111     auto_sourceList.set("Device","elementblocks");
00112     auto_sourceList.set("Quantum Region Source", "schrodinger");
00113     auto_sourceList.set("Non Quantum Region Source", bQBOnly ? "semiclassical" : "schrodinger");
00114     auto_sourceList.set("Eigenvectors to Import", nEigenvals);
00115     auto_sourceList.set("Use predictor-corrector method", false);
00116     auto_sourceList.set("Include exchange-correlation potential", bXCPot);
00117 
00118     Teuchos::ParameterList& sourceList = poisson_probParams.sublist("Poisson Source", false);
00119     if(poisson_subList.isSublist("Poisson Source"))
00120       sourceList.setParameters( poisson_subList.sublist("Poisson Source") );
00121     sourceList.setParametersNotAlreadySet( auto_sourceList );
00122   }
00123 
00124   {
00125     Teuchos::ParameterList auto_permList;
00126     auto_permList.set("Permittivity Type","Block Dependent");
00127 
00128     Teuchos::ParameterList& permList = poisson_probParams.sublist("Permittivity", false);
00129     if(!poisson_subList.isSublist("Permittivity"))
00130       permList.setParameters( poisson_subList.sublist("Permittivity") );
00131     permList.setParametersNotAlreadySet( auto_permList );
00132   }  
00133 
00134             
00135   if(poisson_subList.isSublist("Dirichlet BCs")) {
00136     Teuchos::ParameterList& tmp = poisson_probParams.sublist("Dirichlet BCs", false);
00137     tmp.setParameters(poisson_subList.sublist("Dirichlet BCs"));
00138   }
00139 
00140   // Copy Parameter list over, or create an empty list if it was omitted (so validation passes)
00141   Teuchos::ParameterList& poisson_params = poisson_probParams.sublist("Parameters", false);
00142   if(poisson_subList.isSublist("Parameters"))
00143     poisson_params.setParameters(poisson_subList.sublist("Parameters"));
00144   else poisson_params.set("Number",0);
00145 
00146   // Copy Response Functions list over, or create an empty list if it was omitted (so validation passes)
00147   Teuchos::ParameterList& poisson_resps = poisson_probParams.sublist("Response Functions", false);
00148   if(poisson_subList.isSublist("Response Functions"))
00149     poisson_resps.setParameters(poisson_subList.sublist("Response Functions"));
00150   else poisson_resps.set("Number",0);  
00151   
00152 
00153   Teuchos::ParameterList& poisson_discList = poisson_appParams->sublist("Discretization", false);
00154   poisson_discList.setParameters(discList);
00155   if(poissonExoOutput.length() > 0) 
00156     poisson_discList.set("Exodus Output File Name",poissonExoOutput);
00157   else poisson_discList.remove("Exodus Output File Name",false); 
00158 
00159   if(poissonXmlFile.length() > 0 and tcomm->getRank() == 0)
00160     Teuchos::writeParameterListToXmlFile(*poisson_appParams, poissonXmlFile);
00161 
00162 
00163 
00164   // Create input parameter list for schrodinger app which mimics a separate input file
00165   Teuchos::RCP<Teuchos::ParameterList> schro_appParams = 
00166     Teuchos::createParameterList("Schrodinger Application Parameters");  
00167   Teuchos::ParameterList& schro_probParams = schro_appParams->sublist("Problem",false);
00168 
00169   schro_probParams.set("Name", QCAD::strdim("Schrodinger",numDims));
00170   schro_probParams.set("Solution Method", "Continuation");
00171   schro_probParams.set("Phalanx Graph Visualization Detail", vizDetail);
00172   schro_probParams.set("Energy Unit In Electron Volts",energyUnit);
00173   schro_probParams.set("Length Unit In Meters",lenUnit);
00174   schro_probParams.set("MaterialDB Filename", matrlFile);
00175   schro_probParams.set("Only solve in quantum blocks", bQBOnly);
00176 
00177   {
00178     Teuchos::ParameterList auto_potList;
00179     auto_potList.set("Type","From Aux Data Vector");
00180     auto_potList.set("Aux Index", 0); //we import potential using aux vector 0
00181     
00182     Teuchos::ParameterList& potList = schro_probParams.sublist("Potential", false);
00183     if(schro_subList.isSublist("Potential"))
00184       potList.setParameters(schro_subList.sublist("Potential"));
00185     potList.setParametersNotAlreadySet(auto_potList);
00186   }
00187 
00188   if(!schro_subList.isSublist("Dirichlet BCs") && poisson_subList.isSublist("Dirichlet BCs")) {
00189     Teuchos::ParameterList& schro_dbcList = schro_probParams.sublist("Dirichlet BCs", false);
00190     const Teuchos::ParameterList& poisson_dbcList = poisson_subList.sublist("Dirichlet BCs");
00191     Teuchos::ParameterList::ConstIterator it;
00192     for(it = poisson_dbcList.begin(); it != poisson_dbcList.end(); ++it) {
00193       std::string dbcName = poisson_dbcList.name(it);
00194       std::size_t k = dbcName.find("Phi");
00195       if( k != std::string::npos ) {
00196   dbcName.replace(k, 3 /* len("Phi") */, "psi");  // replace Phi -> psi
00197   schro_dbcList.set( dbcName, 0.0 ); //copy all poisson DBCs but set to zero
00198       }
00199     }
00200   }
00201 
00202   // Copy Parameter list over, or create an empty list if it was omitted (so validation passes)
00203   Teuchos::ParameterList& schro_params = schro_probParams.sublist("Parameters", false);
00204   if(schro_subList.isSublist("Parameters"))
00205     schro_params.setParameters(schro_subList.sublist("Parameters"));
00206   else schro_params.set("Number",0);
00207 
00208   // Copy Response Functions list over, or create an empty list if it was omitted (so validation passes)
00209   Teuchos::ParameterList& schro_resps = schro_probParams.sublist("Response Functions", false);
00210   if(schro_subList.isSublist("Response Functions"))
00211     schro_resps.setParameters(schro_subList.sublist("Response Functions"));
00212   else schro_resps.set("Number",0);
00213 
00214   Teuchos::ParameterList& schro_discList = schro_appParams->sublist("Discretization", false);
00215   schro_discList.setParameters(discList);
00216   if(schrodingerExoOutput.length() > 0) 
00217     schro_discList.set("Exodus Output File Name",schrodingerExoOutput);
00218   else schro_discList.remove("Exodus Output File Name",false); 
00219 
00220   if(schrodingerXmlFile.length() > 0 and tcomm->getRank() == 0)
00221     Teuchos::writeParameterListToXmlFile(*schro_appParams, schrodingerXmlFile);
00222 
00223     
00224   //TODO: need to add meshmover initialization, as in Albany::Application constructor??
00225 
00226   //Create a dummy solverFactory for validating application parameter lists
00227   Albany::SolverFactory validFactory( Teuchos::createParameterList("Empty dummy for Validation"), mcomm );
00228   Teuchos::RCP<const Teuchos::ParameterList> validAppParams = validFactory.getValidAppParameters();
00229   Teuchos::RCP<const Teuchos::ParameterList> validParameterParams = validFactory.getValidParameterParameters();
00230   Teuchos::RCP<const Teuchos::ParameterList> validResponseParams = validFactory.getValidResponseParameters();
00231 
00235 
00236   saved_initial_guess = initial_guess;
00237 
00238   // Validate common parts of poisson app param list: may move inside individual Problem class
00239   poisson_appParams->validateParametersAndSetDefaults(*validAppParams,0);
00240   poisson_appParams->sublist("Problem").sublist("Parameters").validateParameters(*validParameterParams, 0);
00241   poisson_appParams->sublist("Problem").sublist("Response Functions").validateParameters(*validResponseParams, 0);
00242   poissonApp = Teuchos::rcp(new Albany::Application(comm, poisson_appParams, Teuchos::null)); //validates problem params
00243 
00244   // Create model evaluator
00245   Albany::ModelFactory poissonModelFactory(poisson_appParams, poissonApp);
00246   poissonModel = poissonModelFactory.create();
00247 
00248 
00252 
00253   // Validate common parts of schrodinger app param list: may move inside individual Problem class
00254   schro_appParams->validateParametersAndSetDefaults(*validAppParams,0);
00255   schro_appParams->sublist("Problem").sublist("Parameters").validateParameters(*validParameterParams, 0);
00256   schro_appParams->sublist("Problem").sublist("Response Functions").validateParameters(*validResponseParams, 0);
00257   schrodingerApp = Teuchos::rcp(new Albany::Application(comm, schro_appParams, Teuchos::null)); //validates problem params
00258 
00259   // Create model evaluator
00260   Albany::ModelFactory schrodingerModelFactory(schro_appParams, schrodingerApp);
00261   schrodingerModel = schrodingerModelFactory.create();
00262 
00263   //Save the discretization's maps for convenience (should be the same for Poisson and Schrodinger apps)
00264   disc_map = poissonApp->getMap();
00265   disc_overlap_map =  poissonApp->getStateMgr().getDiscretization()->getOverlapMap();
00266   
00267   //Create map for the entire coupled S-P application from the maps from the individual Poisson and Schrodinger applications
00268   //  We need to create a map which is the product of 1 disc_map (for P), N disc_maps (for S's), +N extra (for norm. eqns)
00269   //  in such a way that the elements for each disc_map are contiguous in index space (so that we can easily get Epetra vector views
00270   //  to them separately)
00271   combined_SP_map = QCAD::CreateCombinedMap(disc_map, 1+nEigenvals, nEigenvals, comm);
00272 
00273   // Parameter vectors:  Parameter vectors of coupled PS model evaluator are just the parameter vectors
00274   //   of the Poisson then Schrodinger model evaluators (in order).
00275 
00276   //Get the number of parameter vectors of the Poisson model evaluator
00277   EpetraExt::ModelEvaluator::InArgs poisson_inArgs = poissonModel->createInArgs();
00278   num_poisson_param_vecs = poisson_inArgs.Np();
00279 
00280   //Get the number of parameter vectors of the Schrodginer model evaluator
00281   EpetraExt::ModelEvaluator::InArgs schrodinger_inArgs = schrodingerModel->createInArgs();
00282   num_schrodinger_param_vecs = schrodinger_inArgs.Np();
00283 
00284   num_param_vecs = num_poisson_param_vecs + num_schrodinger_param_vecs;
00285 
00286   // Create sacado parameter vectors of appropriate size for use in evalModel
00287   poisson_sacado_param_vec.resize(num_poisson_param_vecs);
00288   schrodinger_sacado_param_vec.resize(num_schrodinger_param_vecs);
00289 
00290   // Response vectors:  Response vectors of coupled PS model evaluator are just the response vectors
00291   //   of the Poisson then Schrodinger model evaluators (in order).
00292   num_response_vecs = poissonApp->getNumResponses() + schrodingerApp->getNumResponses();
00293 
00294 
00295   // Set member variables based on parameters from the main list
00296   temperature = Temp;
00297   length_unit_in_m  = lenUnit;
00298   energy_unit_in_eV = energyUnit;
00299 
00300 
00301   // Get conduction band offset from reference level (solution to poisson problem), as this
00302   //   is needed to convert the poisson solution vector to conduction band values expected by the schrodinger problem
00303 
00304     // Material database
00305   std::string mtrlDbFilename = "materials.xml";
00306   if(problemParams.isType<std::string>("MaterialDB Filename"))
00307     mtrlDbFilename = problemParams.get<std::string>("MaterialDB Filename");
00308   materialDB = Teuchos::rcp(new QCAD::MaterialDatabase(mtrlDbFilename, comm));
00309   
00310   std::string refMtrlName = materialDB->getParam<std::string>("Reference Material");
00311   double refmatChi = materialDB->getMaterialParam<double>(refMtrlName,"Electron Affinity");
00312 
00313   // compute energy reference
00314   double qPhiRef;
00315   {
00316     const double kB = 8.617332e-5;  // eV/K
00317     std::string category = materialDB->getMaterialParam<std::string>(refMtrlName,"Category");
00318     if (category == "Semiconductor") 
00319     {
00320       // Same qPhiRef needs to be used for the entire structure
00321       double mdn = materialDB->getMaterialParam<double>(refMtrlName,"Electron DOS Effective Mass");
00322       double mdp = materialDB->getMaterialParam<double>(refMtrlName,"Hole DOS Effective Mass");
00323       double Chi = materialDB->getMaterialParam<double>(refMtrlName,"Electron Affinity");
00324       double Eg0 = materialDB->getMaterialParam<double>(refMtrlName,"Zero Temperature Band Gap");
00325       double alpha = materialDB->getMaterialParam<double>(refMtrlName,"Band Gap Alpha Coefficient");
00326       double beta = materialDB->getMaterialParam<double>(refMtrlName,"Band Gap Beta Coefficient");
00327       
00328       double Eg = Eg0 - alpha*pow(temperature,2.0)/(beta+temperature); // in [eV]
00329       double kbT = kB * temperature;      // in [eV]
00330       double Eic = -Eg/2. + 3./4.*kbT*log(mdp/mdn);  // (Ei-Ec) in [eV]
00331       qPhiRef = Chi - Eic;  // (Evac-Ei) in [eV] where Evac = vacuum level
00332     }
00333     else 
00334       TEUCHOS_TEST_FOR_EXCEPTION (true, Teuchos::Exceptions::InvalidParameter, std::endl 
00335         << "Error!  Invalid category " << category << " for reference material !" << std::endl);
00336   } 
00337 
00338   // NOTE: this works for element blocks of the reference material, but really needs to have refmatChi replaced by the
00339   //     chi (electron affinity) for the element block that owns each node...
00340   this->offset_to_CB = qPhiRef - refmatChi; // Conduction Band = offset - poisson_solution
00341 
00342 
00343   // Add discretization parameters to appParams to describe to discretization object how to 
00344   //   interpret the "combined" solution vector used by the coupled P-S solver.
00345   std::vector<std::string> solnVecComps( 2*(1+nEigenvals) ), residVecComps( 2*(1+nEigenvals) );
00346   
00347   solnVecComps[0] = "solution"; solnVecComps[1] = "S"; //was "Potential" but keep as "solution" for backward compatibility
00348   residVecComps[0] = "PoissonRes"; residVecComps[1] = "S";
00349   for(int i=0; i<nEigenvals; i++) {
00350     std::ostringstream ss1; ss1 << "Eigenvector" << i;
00351     solnVecComps[ 2*(i+1) ] = ss1.str(); solnVecComps[ 2*(i+1)+1 ] = "S";
00352     std::ostringstream ss2; ss2 << "SchroRes" << i;
00353     residVecComps[ 2*(i+1) ] = ss2.str(); ; residVecComps[ 2*(i+1)+1 ] = "S";
00354   }
00355   
00356   discList.set("Solution Vector Components", Teuchos::Array<std::string>(solnVecComps));
00357   discList.set("Residual Vector Components", Teuchos::Array<std::string>(residVecComps));
00358   discList.set("Interleaved Ordering", false); //combined vector is concatenated, not "interleaved"
00359 
00360   /* -- Example XML this would generate --
00361      <Parameter name="Solution Vector Components" type="Array(string)" value="{Potential, S, Eigenvector0, S, Eigenvector1, S}"/>
00362      <Parameter name="Residual Vector Components" type="Array(string)" value="{PoissonRes, S, SchroRes0, S, SchroRes1, S}"/>
00363      <Parameter name="Interleaved Ordering" type="bool" value="false"/>
00364   */
00365 
00366   // Create discretization object solely for producing collected output
00367   Albany::DiscretizationFactory discFactory(appParams, comm);
00368 
00369   // Get mesh specification object: worksetSize, cell topology, etc
00370   Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> > meshSpecs =
00371     discFactory.createMeshSpecs();
00372 
00373   Albany::AbstractFieldContainer::FieldContainerRequirements requirements; //empty?
00374   Teuchos::RCP<Albany::StateInfoStruct> stateInfo = poissonApp->getStateMgr().getStateInfoStruct(); //for now, just use Poisson app's states
00375                             //Teuchos::rcp(new Albany::StateInfoStruct); //empty
00376 
00377   int neq = 1+nEigenvals; // number of mesh-equations
00378   Teuchos::RCP<Piro::MLRigidBodyModes> rigidBodyModes(Teuchos::rcp(new Piro::MLRigidBodyModes(neq)));
00379   disc = discFactory.createDiscretization(neq, stateInfo,requirements,rigidBodyModes);
00380 
00381   myComm = comm;
00382 }
00383 
00384 QCAD::CoupledPoissonSchrodinger::~CoupledPoissonSchrodinger()
00385 {
00386 }
00387 
00388 
00389 Teuchos::RCP<const Epetra_Map> QCAD::CoupledPoissonSchrodinger::get_x_map() const
00390 {
00391   return combined_SP_map;
00392 }
00393 
00394 Teuchos::RCP<const Epetra_Map> QCAD::CoupledPoissonSchrodinger::get_f_map() const
00395 {
00396   return combined_SP_map;
00397 }
00398 
00399 Teuchos::RCP<const Epetra_Map> QCAD::CoupledPoissonSchrodinger::get_p_map(int l) const
00400 {
00401   TEUCHOS_TEST_FOR_EXCEPTION(l >= num_param_vecs || l < 0, Teuchos::Exceptions::InvalidParameter,
00402                      std::endl <<
00403                      "Error in QCAD::CoupledPoissonSchrodinger::get_p_map():  " <<
00404                      "Invalid parameter index l = " << l << std::endl);
00405   if(l < num_poisson_param_vecs)
00406     return poissonModel->get_p_map(l);
00407   else
00408     return schrodingerModel->get_p_map(l - num_poisson_param_vecs);
00409 }
00410 
00411 Teuchos::RCP<const Epetra_Map> QCAD::CoupledPoissonSchrodinger::get_g_map(int j) const
00412 {
00413   TEUCHOS_TEST_FOR_EXCEPTION(j > num_response_vecs || j < 0, Teuchos::Exceptions::InvalidParameter,
00414                      std::endl <<
00415                      "Error in QCAD::CoupledPoissonSchrodinger::get_g_map():  " <<
00416                      "Invalid response index j = " << j << std::endl);
00417   
00418   if(j < poissonApp->getNumResponses())
00419     return poissonModel->get_g_map(j);
00420   else
00421     return schrodingerModel->get_g_map(j - poissonApp->getNumResponses());
00422 }
00423 
00424 Teuchos::RCP<const Teuchos::Array<std::string> > QCAD::CoupledPoissonSchrodinger::get_p_names(int l) const
00425 {
00426   TEUCHOS_TEST_FOR_EXCEPTION(l >= num_param_vecs || l < 0, 
00427          Teuchos::Exceptions::InvalidParameter,
00428                      std::endl << 
00429                      "Error!  Albany::ModelEvaluator::get_p_names():  " <<
00430                      "Invalid parameter index l = " << l << std::endl);
00431   if(l < num_poisson_param_vecs)
00432     return poissonModel->get_p_names(l);
00433   else
00434     return schrodingerModel->get_p_names(l - num_poisson_param_vecs);
00435 }
00436 
00437 
00438 Teuchos::RCP<const Epetra_Vector> QCAD::CoupledPoissonSchrodinger::get_x_init() const
00439 {
00440   if(saved_initial_guess != Teuchos::null) {
00441     std::cout << "DEBUG CPS: returning saved initial guess!" << std::endl;
00442     return saved_initial_guess;
00443   }
00444 
00445   //Put together x_init's from Poisson and Schrodinger for now (but does this make sense for eigenvectors?) -- TODO: discuss
00446   Teuchos::RCP<const Epetra_Vector> poisson_x_init = poissonModel->get_x_init(); // should have disc_map
00447   Teuchos::RCP<const Epetra_Vector> schrodinger_x_init = schrodingerModel->get_x_init(); // should have disc_map
00448   
00449   Teuchos::RCP<Epetra_Vector> x_init = Teuchos::rcp(new Epetra_Vector(*combined_SP_map));
00450   Teuchos::RCP<Epetra_Vector> x_init_poisson;
00451   Teuchos::RCP<Epetra_MultiVector> x_init_schrodinger;
00452 
00453   separateCombinedVector(x_init, x_init_poisson, x_init_schrodinger);
00454 
00455   std::vector<int> localInds( poisson_x_init->MyLength() );
00456   for(int i=0; i < poisson_x_init->MyLength(); i++) localInds[i] = i;
00457 
00458   x_init_poisson->ReplaceMyValues( poisson_x_init->MyLength(), &(*poisson_x_init)[0], &localInds[0] );
00459   for(int k=0; k < nEigenvals; k++)
00460     (*x_init_schrodinger)(k)->ReplaceMyValues( schrodinger_x_init->MyLength(), &(*schrodinger_x_init)[0], &localInds[0] ); //localInds are the same
00461   
00462   return x_init;
00463 }
00464 
00465 Teuchos::RCP<const Epetra_Vector> QCAD::CoupledPoissonSchrodinger::get_x_dot_init() const
00466 {
00467   //Put together x_dot_init's from Poisson and Schrodinger for now (but does this make sense for eigenvectors?) -- TODO: discuss
00468   Teuchos::RCP<const Epetra_Vector> poisson_x_dot_init = poissonModel->get_x_dot_init(); // should have disc_map
00469   Teuchos::RCP<const Epetra_Vector> schrodinger_x_dot_init = schrodingerModel->get_x_dot_init(); // should have disc_map
00470   
00471   Teuchos::RCP<Epetra_Vector> x_dot_init = Teuchos::rcp(new Epetra_Vector(*combined_SP_map));
00472   Teuchos::RCP<Epetra_Vector> x_dot_init_poisson;
00473   Teuchos::RCP<Epetra_MultiVector> x_dot_init_schrodinger;
00474 
00475   separateCombinedVector(x_dot_init, x_dot_init_poisson, x_dot_init_schrodinger);
00476 
00477   std::vector<int> localInds( poisson_x_dot_init->MyLength() );
00478   for(int i=0; i < poisson_x_dot_init->MyLength(); i++) localInds[i] = i;
00479 
00480   x_dot_init_poisson->ReplaceMyValues( poisson_x_dot_init->MyLength(), &(*poisson_x_dot_init)[0], &localInds[0] );
00481   for(int k=0; k < nEigenvals; k++)
00482     (*x_dot_init_schrodinger)(k)->ReplaceMyValues( schrodinger_x_dot_init->MyLength(), &(*schrodinger_x_dot_init)[0], &localInds[0] ); //same localInds are the same
00483   
00484   //Teuchos::RCP<const Epetra_Vector> const_x_dot_init = Teuchos::rcp(new const Epetra_Vector(*x_dot_init));
00485   return x_dot_init;
00486 }
00487 
00488 
00489 Teuchos::RCP<const Epetra_Vector> QCAD::CoupledPoissonSchrodinger::get_p_init(int l) const
00490 {
00491   TEUCHOS_TEST_FOR_EXCEPTION(l >= num_param_vecs || l < 0, Teuchos::Exceptions::InvalidParameter,
00492                      std::endl <<
00493                      "Error in QCAD::CoupledPoissonSchrodinger::get_p_init():  " <<
00494                      "Invalid parameter index l = " << l << std::endl);
00495 
00496   if(l < num_poisson_param_vecs)
00497     return poissonModel->get_p_init(l);
00498   else
00499     return schrodingerModel->get_p_init(l - num_poisson_param_vecs);
00500 }
00501 
00502 
00503 Teuchos::RCP<Epetra_Operator>
00504 QCAD::CoupledPoissonSchrodinger::create_W() const
00505 {
00506   // Get material parameters for quantum region, used in computing quantum density
00507   std::string quantumMtrlName = materialDB->getParam<std::string>("Quantum Material");
00508   int valleyDegeneracyFactor = materialDB->getMaterialParam<int>(quantumMtrlName,"Number of conduction band min",2);
00509   double effMass = materialDB->getMaterialParam<double>(quantumMtrlName,"Transverse Electron Effective Mass");
00510 
00511   std::cout << "DEBUG:  CPS create_W called!!" << std::endl;
00512   return Teuchos::rcp( new QCAD::CoupledPSJacobian(nEigenvals, disc_map, combined_SP_map, myComm, 
00513                numDims, valleyDegeneracyFactor, temperature,
00514                length_unit_in_m, energy_unit_in_eV, effMass, offset_to_CB) );
00515 }
00516 
00517 Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
00518 QCAD::CoupledPoissonSchrodinger::create_WPrec() const
00519 {
00520   //std::cout << "DEBUG:  CPS create_WPrec called!!" << std::endl;
00521   Teuchos::RCP<Epetra_Operator> precOp = Teuchos::rcp( new QCAD::CoupledPSPreconditioner(nEigenvals, disc_map, combined_SP_map, myComm) );
00522 
00523   // bool is answer to: "Prec is already inverted?"
00524   return Teuchos::rcp(new EpetraExt::ModelEvaluator::Preconditioner(precOp,true));
00525 }
00526 
00527 Teuchos::RCP<Epetra_Operator>
00528 QCAD::CoupledPoissonSchrodinger::create_DgDx_op(int j) const
00529 {
00530   TEUCHOS_TEST_FOR_EXCEPTION(
00531     j >= num_response_vecs || j < 0, 
00532     Teuchos::Exceptions::InvalidParameter,
00533     std::endl << 
00534     "Error!  Albany::ModelEvaluator::create_DgDx_op():  " << 
00535     "Invalid response index j = " << j << std::endl);
00536   
00537   if(j < poissonApp->getNumResponses())
00538     return poissonApp->getResponse(j)->createGradientOp();
00539   else
00540     return schrodingerApp->getResponse(j - poissonApp->getNumResponses())->createGradientOp();
00541 }
00542 
00543 Teuchos::RCP<Epetra_Operator>
00544 QCAD::CoupledPoissonSchrodinger::create_DgDx_dot_op(int j) const
00545 {
00546   TEUCHOS_TEST_FOR_EXCEPTION(
00547     j >= num_response_vecs || j < 0, 
00548     Teuchos::Exceptions::InvalidParameter,
00549     std::endl << 
00550     "Error!  Albany::ModelEvaluator::create_DgDx_dot_op():  " << 
00551     "Invalid response index j = " << j << std::endl);
00552   
00553   if(j < poissonApp->getNumResponses())
00554     return poissonApp->getResponse(j)->createGradientOp();
00555   else
00556     return schrodingerApp->getResponse(j - poissonApp->getNumResponses())->createGradientOp();
00557 }
00558 
00559 
00560 EpetraExt::ModelEvaluator::InArgs QCAD::CoupledPoissonSchrodinger::createInArgs() const
00561 {
00562   InArgsSetup inArgs;
00563   inArgs.setModelEvalDescription("QCAD Coupled Poisson-Schrodinger Model Evaluator");
00564 
00565   inArgs.setSupports(IN_ARG_t,true);
00566   inArgs.setSupports(IN_ARG_x,true);
00567   inArgs.setSupports(IN_ARG_x_dot,true);
00568   inArgs.setSupports(IN_ARG_alpha,true);
00569   inArgs.setSupports(IN_ARG_beta,true);
00570   inArgs.set_Np(num_param_vecs);
00571 
00572   // Note: no SG support yet...
00573 
00574   return inArgs;
00575 }
00576 
00577 EpetraExt::ModelEvaluator::OutArgs QCAD::CoupledPoissonSchrodinger::createOutArgs() const
00578 {
00579   OutArgsSetup outArgs;
00580   outArgs.setModelEvalDescription("QCAD Coupled Poisson-Schrodinger Model Evaluator");
00581 
00582   int n_g = num_response_vecs;
00583   bool bScalarResponse;
00584 
00585   // Deterministic
00586   outArgs.setSupports(OUT_ARG_f,true);
00587   outArgs.setSupports(OUT_ARG_W,true);
00588   outArgs.set_W_properties(
00589     DerivativeProperties(DERIV_LINEARITY_UNKNOWN, DERIV_RANK_FULL, true));
00590   outArgs.setSupports(OUT_ARG_WPrec, true);
00591   outArgs.set_Np_Ng(num_param_vecs, n_g);
00592 
00593   for (int i=0; i<num_param_vecs; i++)
00594     outArgs.setSupports(OUT_ARG_DfDp, i, DerivativeSupport(DERIV_MV_BY_COL));
00595   for (int i=0; i<n_g; i++) {
00596 
00597     if(i < poissonApp->getNumResponses())
00598       bScalarResponse = poissonApp->getResponse(i)->isScalarResponse();
00599     else
00600       bScalarResponse = schrodingerApp->getResponse(i - poissonApp->getNumResponses())->isScalarResponse();
00601 
00602     if(bScalarResponse) {
00603       outArgs.setSupports(OUT_ARG_DgDx, i,
00604                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00605       outArgs.setSupports(OUT_ARG_DgDx_dot, i,
00606                           DerivativeSupport(DERIV_TRANS_MV_BY_ROW));
00607     }
00608     else {
00609       outArgs.setSupports(OUT_ARG_DgDx, i,
00610                           DerivativeSupport(DERIV_LINEAR_OP));
00611       outArgs.setSupports(OUT_ARG_DgDx_dot, i,
00612                           DerivativeSupport(DERIV_LINEAR_OP));
00613     }
00614 
00615     for (int j=0; j<num_param_vecs; j++)
00616       outArgs.setSupports(OUT_ARG_DgDp, i, j,
00617                           DerivativeSupport(DERIV_MV_BY_COL));
00618   }
00619 
00620   //Note: no SG support yet...
00621 
00622   return outArgs;
00623 }
00624 
00625 
00626 void 
00627 QCAD::CoupledPoissonSchrodinger::evalModel(const InArgs& inArgs,
00628       const OutArgs& outArgs ) const
00629 {
00630   //?? Teuchos::TimeMonitor Timer(*timer); //start timer
00631 
00632   //
00633   // Get the input arguments
00634   //
00635   Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
00636   Teuchos::RCP<const Epetra_Vector> x_dot;
00637   double alpha     = 0.0;  // M coeff
00638   double beta      = 1.0;  // J coeff
00639   double curr_time = 0.0;
00640   x_dot = inArgs.get_x_dot();
00641   if (x_dot != Teuchos::null) {
00642     alpha = inArgs.get_alpha();
00643     beta = inArgs.get_beta();
00644     curr_time  = inArgs.get_t();
00645     std::cout << "DEBUG: WARNING: x_dot given to CoupledPoissonSchrodinger evalModel!!" << std::endl;
00646   }
00647   for (int i=0; i<inArgs.Np(); i++) {
00648     Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i);
00649     if (p != Teuchos::null) {
00650       if(i < num_poisson_param_vecs) {
00651   for (unsigned int j=0; j<poisson_sacado_param_vec[i].size(); j++)
00652     poisson_sacado_param_vec[i][j].baseValue = (*p)[j];
00653       }
00654       else {
00655   for (unsigned int j=0; j<schrodinger_sacado_param_vec[i-num_poisson_param_vecs].size(); j++)
00656     schrodinger_sacado_param_vec[i-num_poisson_param_vecs][j].baseValue = (*p)[j];
00657       }
00658     }
00659   }
00660 
00661 
00662   //
00663   // Get the output arguments
00664   //
00665   EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out = outArgs.get_f();
00666   Teuchos::RCP<Epetra_Operator> W_out = outArgs.get_W();
00667 
00668   // Get preconditioner operator, if requested
00669   Teuchos::RCP<Epetra_Operator> WPrec_out;
00670   if (outArgs.supports(OUT_ARG_WPrec)) WPrec_out = outArgs.get_WPrec();
00671 
00672 
00673   //
00674   // Get views into 'x' (and 'xdot'?) vectors to use for separate poisson and schrodinger application object calls
00675   //
00676   int disc_nMyElements = disc_map->NumMyElements();
00677 
00678   Teuchos::RCP<const Epetra_Vector> x_poisson, xdot_poisson, eigenvals_dist;
00679   Teuchos::RCP<const Epetra_MultiVector> x_schrodinger, xdot_schrodinger;
00680   std::vector<const Epetra_Vector*> xdot_schrodinger_vec(nEigenvals);
00681   separateCombinedVector(x, x_poisson, x_schrodinger, eigenvals_dist);
00682     
00683   if (x_dot != Teuchos::null) {  //maybe unnecessary - it seems that the coupled PS model evaluator shouldn't support x_dot ...
00684     separateCombinedVector(x_dot, xdot_poisson, xdot_schrodinger);
00685     for(int i=0; i<nEigenvals; i++) xdot_schrodinger_vec[i] = (*xdot_schrodinger)(i);
00686   }
00687   else {
00688     xdot_poisson = Teuchos::null;
00689     for(int i=0; i<nEigenvals; i++) 
00690       xdot_schrodinger_vec[i] = NULL;
00691   }
00692 
00693   //
00694   // Communicate all the eigenvalues to every processor, since all parts of the mesh need them
00695   //
00696   int my_nEigenvals = combined_SP_map->NumMyElements() - disc_nMyElements * (1+nEigenvals);
00697   Epetra_Map dist_eigenval_map(nEigenvals, my_nEigenvals, 0, *myComm);
00698   Epetra_LocalMap local_eigenval_map(nEigenvals, 0, *myComm);
00699   Epetra_Import eigenval_importer(local_eigenval_map, dist_eigenval_map);
00700 
00701   Teuchos::RCP<Epetra_Vector> eigenvals =  Teuchos::rcp(new Epetra_Vector(local_eigenval_map));
00702   eigenvals->Import(*eigenvals_dist, eigenval_importer, Insert);
00703   Teuchos::RCP<std::vector<double> > stdvec_eigenvals = Teuchos::rcp(new std::vector<double>(&(*eigenvals)[0], &(*eigenvals)[0] + nEigenvals));
00704 
00705   //
00706   // Get views into 'f' residual vector to use for separate poisson and schrodinger application object calls
00707   //
00708   Teuchos::RCP<Epetra_Vector>   f_poisson, f_norm_local, f_norm_dist;
00709   Teuchos::RCP<Epetra_MultiVector> f_schrodinger;
00710   std::vector<Epetra_Vector*> f_schrodinger_vec(nEigenvals);
00711 
00712   if(f_out != Teuchos::null) {
00713     separateCombinedVector(f_out, f_poisson, f_schrodinger, f_norm_dist);
00714     for(int i=0; i<nEigenvals; i++) f_schrodinger_vec[i] = (*f_schrodinger)(i);
00715 
00716     // Create local vector for holding the residual of the normalization equations on each proc.
00717     //   (later we sum all procs contributions together and copy into distributed f_norm_dist vector)
00718     f_norm_local = Teuchos::rcp(new Epetra_Vector(local_eigenval_map));
00719   }
00720   else {
00721     f_poisson = Teuchos::null;
00722     for(int i=0; i<nEigenvals; i++) f_schrodinger_vec[i] = NULL;
00723     f_norm_local = f_norm_dist = Teuchos::null;
00724   }
00725 
00726 
00727   // Create an eigendata struct for passing the eigenvectors to the poisson app
00728   //  -- note that this requires the *overlapped* eigenvectors
00729   Teuchos::RCP<Albany::EigendataStruct> eigenData = Teuchos::rcp( new Albany::EigendataStruct );
00730 
00731   eigenData->eigenvalueRe = stdvec_eigenvals;
00732   eigenData->eigenvectorRe = 
00733     Teuchos::rcp(new Epetra_MultiVector(*disc_overlap_map, nEigenvals));
00734   eigenData->eigenvectorIm = Teuchos::null; // no imaginary eigenvalue data... 
00735 
00736     // Importer for overlapped data
00737   Teuchos::RCP<Epetra_Import> overlap_importer =
00738     Teuchos::rcp(new Epetra_Import(*disc_overlap_map, *disc_map));
00739 
00740     // Overlapped eigenstate vectors
00741   for(int i=0; i<nEigenvals; i++) {
00742     (*(eigenData->eigenvectorRe))(i)->Import( *((*x_schrodinger)(i)), *overlap_importer, Insert );
00743     //(*(eigenData->eigenvectorRe))(i)->PutScalar(0.0); //DEBUG - zero out eigenvectors passed to Poisson
00744   }
00745 
00746     // set eigenvalues / eigenvectors for use in poisson problem:
00747   poissonApp->getStateMgr().setEigenData(eigenData);
00748 
00749 
00750   // Get overlapped version of potential (x_poisson) for passing as auxData to schrodinger app
00751   Teuchos::RCP<Epetra_MultiVector> overlapped_V = Teuchos::rcp(new Epetra_MultiVector(*disc_overlap_map, 1));
00752   Teuchos::RCP<Epetra_Vector> ones_vec = Teuchos::rcp(new Epetra_Vector(*disc_overlap_map));
00753   ones_vec->PutScalar(1.0);
00754   (*overlapped_V)(0)->Import( *x_poisson, *overlap_importer, Insert );
00755   (*overlapped_V)(0)->Update(offset_to_CB, *ones_vec, -1.0);
00756   //std::cout << "DEBUG: Offset to conduction band = " << offset_to_CB << std::endl;
00757 
00758   // set potential for use in schrodinger problem
00759   schrodingerApp->getStateMgr().setAuxData(overlapped_V);
00760 
00761 
00762   
00763   //
00764   // Compute the functions
00765   //
00766   bool f_poisson_already_computed = false;
00767   std::vector<bool> f_schrodinger_already_computed(nEigenvals, false);
00768 
00769   Teuchos::RCP<Epetra_CrsMatrix> W_out_poisson_crs; //possibly used by preconditioner, so declare here
00770   Teuchos::RCP<Epetra_CrsMatrix> W_out_schrodinger_crs; //possibly used by preconditioner, so declare here
00771 
00772   // Mass Matrix -- needed even if we don't need to compute the Jacobian, since it enters into the normalization equations
00773   //   --> Compute mass matrix using schrodinger equation -- independent of eigenvector so can just use 0th
00774   //       Note: to compute this, we need to evaluate the schrodinger problem as a transient problem, so create a dummy xdot...
00775   Teuchos::RCP<Epetra_Operator> M_out_schrodinger = schrodingerModel->create_W(); //maybe re-use this and not create it every time?
00776   Teuchos::RCP<Epetra_CrsMatrix> M_out_schrodinger_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(M_out_schrodinger, true);
00777   Teuchos::RCP<const Epetra_Vector> dummy_xdot = schrodingerModel->get_x_dot_init(); // I think this would work as well: Teuchos::rcp(new Epetra_Vector(*disc_map)) 
00778   schrodingerApp->computeGlobalJacobian(1.0, 0.0, 0.0, curr_time, dummy_xdot.get(), NULL, *((*x_schrodinger)(0)), 
00779               schrodinger_sacado_param_vec, f_schrodinger_vec[0], *M_out_schrodinger_crs);
00780 
00781 
00782   // Hamiltionan Matrix -- needed even if we don't need to compute the Jacobian, since this is how we compute the schrodinger residuals
00783   //   --> Computed as jacobian matrix of schrodinger equation -- independent of eigenvector so can just use 0th
00784   Teuchos::RCP<Epetra_Operator> J_out_schrodinger = schrodingerModel->create_W(); //maybe re-use this and not create it every time?
00785   Teuchos::RCP<Epetra_CrsMatrix> J_out_schrodinger_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(J_out_schrodinger, true);
00786   schrodingerApp->computeGlobalJacobian(0.0, 1.0, 0.0, curr_time, dummy_xdot.get(), NULL, *((*x_schrodinger)(0)), 
00787               schrodinger_sacado_param_vec, f_schrodinger_vec[0], *J_out_schrodinger_crs);
00788 
00789   f_schrodinger_already_computed[0] = true; //residual is not affected by alpha & beta, so both of the above calls compute it.
00790 
00791 
00792   // W 
00793   if (W_out != Teuchos::null) { 
00794     // W = alpha*M + beta*J where M is mass mx and J is jacobian
00795 
00796     //if we need to compute the jacobian, get the jacobians of the poisson and schrodinger
00797     //  applications (as crs matrices), as well as the mass matrix (from the schrodinger problem,
00798     //  since it includes xdot - maybe need to fabricate this??) and from these construct a CoupledPoissonSchrodingerJacobian object (an Epetra_Operator)
00799     
00800     //TODO - how to allow general alpha and beta?  This won't work given current logic, so we should test that alpha=0, beta=1 and throw an error otherwise...
00801 
00802     // Compute poisson Jacobian
00803     Teuchos::RCP<Epetra_Operator> W_out_poisson = poissonModel->create_W(); //maybe re-use this and not create it every time?
00804     W_out_poisson_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out_poisson, true);
00805 
00806     poissonApp->computeGlobalJacobian(alpha, beta, 0.0, curr_time, xdot_poisson.get(), NULL, *x_poisson, 
00807               poisson_sacado_param_vec, f_poisson.get(), *W_out_poisson_crs);
00808     f_poisson_already_computed = true;
00809 
00810     
00811     TEUCHOS_TEST_FOR_EXCEPTION(nEigenvals <= 0, Teuchos::Exceptions::InvalidParameter,"Error! The number of eigenvalues must be greater than zero.");
00812       
00813     //Compute schrodinger Jacobian using first eigenvector -- independent of eigenvector since Schro. eqn is linear
00814     // (test for cases we've already computed above)
00815     if(alpha == 1.0 && beta == 0.0)
00816       W_out_schrodinger_crs = M_out_schrodinger_crs;
00817     else if(alpha == 0.0 && beta == 1.0)
00818       W_out_schrodinger_crs = J_out_schrodinger_crs;
00819     else {
00820       Teuchos::RCP<Epetra_Operator> W_out_schrodinger = schrodingerModel->create_W(); //maybe re-use this and not create it every time?
00821       W_out_schrodinger_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out_schrodinger, true);
00822       schrodingerApp->computeGlobalJacobian(alpha, beta, 0.0, curr_time, xdot_schrodinger_vec[0], NULL, *((*x_schrodinger)(0)), 
00823             schrodinger_sacado_param_vec, f_schrodinger_vec[0], *W_out_schrodinger_crs);
00824       f_schrodinger_already_computed[0] = true;
00825     }
00826     
00827     Teuchos::RCP<QCAD::CoupledPSJacobian> W_out_psj = Teuchos::rcp_dynamic_cast<QCAD::CoupledPSJacobian>(W_out, true);
00828     W_out_psj->initialize(W_out_poisson_crs, W_out_schrodinger_crs, M_out_schrodinger_crs, eigenvals, x_schrodinger);
00829 
00830 
00831 
00832     /*
00833     // DEBUG --- JACOBIAN TEST -----------------------------------------------------------------------------------------
00834      
00835     Teuchos::RCP<Epetra_Vector> x_plus_dx = Teuchos::rcp( new Epetra_Vector(*combined_SP_map) );
00836     Teuchos::RCP<Epetra_Vector> dx_test = Teuchos::rcp( new Epetra_Vector(*combined_SP_map) );
00837     double eps;
00838     
00839     //Init dx_test
00840     Teuchos::RCP<Epetra_Vector> dx_test_poisson, dx_test_eigenvals;
00841     Teuchos::RCP<Epetra_MultiVector> dx_test_schrodinger;
00842     separateCombinedVector(dx_test, dx_test_poisson, dx_test_schrodinger, dx_test_eigenvals);
00843 
00844     eps = 1e-7;
00845     dx_test->PutScalar(1.0);
00846     //dx_test_poisson->PutScalar(1.0);
00847     //dx_test_schrodinger->PutScalar(1.0);
00848     //dx_test_eigenvals->PutScalar(1.0);
00849 
00850     Teuchos::RCP<Epetra_Vector> f_test_manual = Teuchos::rcp( new Epetra_Vector(*combined_SP_map) );
00851     Teuchos::RCP<Epetra_Vector> f_test_tmp = Teuchos::rcp( new Epetra_Vector(*combined_SP_map) );
00852     Teuchos::RCP<Epetra_Vector> f_test_jac = Teuchos::rcp( new Epetra_Vector(*combined_SP_map) );
00853 
00854     computeResidual(x, f_test_manual, M_out_schrodinger_crs);
00855     //f_test_manual->Print( std::cout << "JACOBIAN TEST:  MANUAL 1" << std::endl);
00856     x_plus_dx->Update(1.0, *x, eps, *dx_test, 0.0); // x + eps*dx
00857     computeResidual(x_plus_dx, f_test_tmp, M_out_schrodinger_crs);
00858     //f_test_tmp->Print( std::cout << "JACOBIAN TEST:  MANUAL 2" << std::endl);
00859     f_test_manual->Update(1.0/eps, *f_test_tmp, -1.0/eps); // f_test_manual = (resid(x+eps*dx) - resid(x)) / eps ~ jacobian * dx
00860 
00861     W_out_psj->Apply(*dx_test, *f_test_jac);
00862 
00863 
00864     //x_schrodinger->Print( std::cout << "JACOBIAN TEST:  X_SCHRODINGER" << std::endl);
00865     f_test_manual->Print( std::cout << "JACOBIAN TEST:  MANUAL DIFF" << std::endl);
00866     f_test_jac->Print( std::cout << "JACOBIAN TEST:  JACOBIAN DIFF" << std::endl);
00867 
00868     f_test_tmp->Update(1.0, *f_test_manual, -1.0, *f_test_jac, 0.0);
00869     f_test_tmp->Print( std::cout << "JACOBIAN TEST:  COMPARISON VECTOR" << std::endl);
00870     
00871     double test_norm;
00872     f_test_tmp->Norm2(&test_norm);
00873     std::cout << "JACOBIAN TEST: COMPARISON VECTOR 2-NORM = " << test_norm << std::endl;
00874 
00875     // DEBUG --- JACOBIAN TEST -----------------------------------------------------------------------------------------
00876     */
00877   }
00878 
00879 
00880   if (WPrec_out != Teuchos::null) {
00881      // Get Poisson Preconditioner
00882      Teuchos::RCP<Epetra_Operator> WPrec_poisson;
00883      
00884      // Get the Poisson Jacobian -- (just copy the it if we already computed it)
00885      Teuchos::RCP<Epetra_CrsMatrix> Extra_W_crs_poisson;
00886      if(W_out != Teuchos::null)
00887        Extra_W_crs_poisson = Teuchos::rcp( new Epetra_CrsMatrix(*W_out_poisson_crs) ); //Check: does this need to copy?
00888      else {
00889        Extra_W_crs_poisson = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(poissonModel->create_W(), true);
00890        poissonApp->computeGlobalJacobian(alpha, beta, 0.0, curr_time, xdot_poisson.get(), NULL, *x_poisson, 
00891            poisson_sacado_param_vec, f_poisson.get(), *Extra_W_crs_poisson);
00892        f_poisson_already_computed = true;
00893      }
00894 
00895      bool poisson_supports_teko_prec = false;  // TODO: I think this should = whether poisson outargs supports OUT_ARG_WPrec
00896      if( poisson_supports_teko_prec ) {
00897        Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner> WPrec_poisson_pre = poissonModel->create_WPrec(); //maybe re-use this and not create it every time?
00898        WPrec_poisson = WPrec_poisson_pre->PrecOp;       
00899        poissonApp->computeGlobalPreconditioner(Extra_W_crs_poisson, WPrec_poisson);
00900      }
00901      else {
00902        // Use Ifpack to get a pseudo inverse of Extra_W_crs_poisson
00903        Teuchos::ParameterList Ifpack_list;
00904        Ifpack Ifpack_factory; // allocate an IFPACK factory.
00905 
00906        // create the preconditioner. -- maybe pull this info from input file in FUTURE
00907        std::string PrecType = "ILU"; // incomplete LU
00908        int OverlapLevel = 1; // must be >= 0. If Comm.NumProc() == 1, it is ignored.
00909 
00910        Teuchos::RCP<Ifpack_Preconditioner> WPrec_poisson_pre = Teuchos::rcp( Ifpack_factory.Create(PrecType, &*Extra_W_crs_poisson, OverlapLevel) );
00911        assert(WPrec_poisson_pre != Teuchos::null);
00912 
00913        // specify parameters for ILU -- maybe pull this info from input file in FUTURE
00914        Ifpack_list.set("fact: drop tolerance", 1e-9);
00915        Ifpack_list.set("fact: level-of-fill", 1);
00916        // the combine mode is one of  the following:
00917        // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax" (Their meaning is as defined in file Epetra_CombineMode.h)
00918        Ifpack_list.set("schwarz: combine mode", "Add");
00919 
00920 
00921        if( WPrec_poisson_pre->SetParameters(Ifpack_list) != 0 ) // sets the parameters
00922    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,"Error! Invalid IFPACK Parameters.");  
00923        if( WPrec_poisson_pre->Initialize() != 0)                // initialize preconditioner (must fillComplete matrix by now)
00924    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,"Error Inializing Ifpack preconditioner");   
00925        if( WPrec_poisson_pre->Compute() != 0)                   // compute preconditioner
00926    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,"Error Computing Ifpack preconditioner");  
00927 
00928        WPrec_poisson = Teuchos::rcp_dynamic_cast<Epetra_Operator>(WPrec_poisson_pre, true);
00929      }
00930 
00931      // Get Schrodinger Preconditioner
00932      Teuchos::RCP<Epetra_Operator> WPrec_schrodinger;
00933 
00934        // Get the Schrodinger Jacobian
00935      Teuchos::RCP<Epetra_CrsMatrix> Extra_W_crs_schrodinger;
00936      if(W_out != Teuchos::null)
00937        Extra_W_crs_schrodinger = Teuchos::rcp( new Epetra_CrsMatrix(*W_out_schrodinger_crs) ); //Check: does this need to copy?
00938      else {
00939        Extra_W_crs_schrodinger = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(schrodingerModel->create_W(), true);
00940        schrodingerApp->computeGlobalJacobian(alpha, beta, 0.0, curr_time, xdot_schrodinger_vec[0], NULL, *((*x_schrodinger)(0)), 
00941                schrodinger_sacado_param_vec, f_schrodinger_vec[0], *Extra_W_crs_schrodinger);
00942        f_schrodinger_already_computed[0] = true;
00943      }
00944 
00945      bool schrodinger_supports_teko_prec = false;  // I think this should = whether poisson outargs supports OUT_ARG_WPrec
00946      if( schrodinger_supports_teko_prec ) {
00947        Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner> WPrec_schrodinger_pre = schrodingerModel->create_WPrec(); //maybe re-use this and not create every time?
00948        WPrec_schrodinger = WPrec_schrodinger_pre->PrecOp;
00949        schrodingerApp->computeGlobalPreconditioner(Extra_W_crs_schrodinger, WPrec_schrodinger);
00950      }
00951      else {
00952        // Use Ifpack to get a pseudo inverse of Extra_W_crs_schrodinger
00953        Teuchos::ParameterList Ifpack_list;
00954        Ifpack Ifpack_factory; // allocate an IFPACK factory.
00955 
00956        // create the preconditioner. -- maybe pull this info from input file in FUTURE
00957        std::string PrecType = "ILU"; // incomplete LU
00958        int OverlapLevel = 1; // must be >= 0. If Comm.NumProc() == 1, it is ignored.
00959 
00960        Teuchos::RCP<Ifpack_Preconditioner> WPrec_schrodinger_pre = Teuchos::rcp( Ifpack_factory.Create(PrecType, &*Extra_W_crs_schrodinger, OverlapLevel) );
00961        assert(WPrec_schrodinger_pre != Teuchos::null);
00962 
00963        // specify parameters for ILU -- maybe pull this info from input file in FUTURE
00964        Ifpack_list.set("fact: drop tolerance", 1e-9);
00965        Ifpack_list.set("fact: level-of-fill", 1);
00966        // the combine mode is one of the following:
00967        // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax"   (Their meaning is as defined in file Epetra_CombineMode.h)
00968        Ifpack_list.set("schwarz: combine mode", "Add");
00969 
00970 
00971        if( WPrec_schrodinger_pre->SetParameters(Ifpack_list) != 0 ) // sets the parameters
00972    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,"Error! Invalid IFPACK Parameters.");  
00973        if( WPrec_schrodinger_pre->Initialize() != 0)                // initialize preconditioner (must fillComplete matrix by now)
00974    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,"Error Inializing Ifpack preconditioner");   
00975        if( WPrec_schrodinger_pre->Compute() != 0)                   // compute preconditioner
00976    TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,"Error Computing Ifpack preconditioner");  
00977 
00978        WPrec_schrodinger = Teuchos::rcp_dynamic_cast<Epetra_Operator>(WPrec_schrodinger_pre, true);
00979      }
00980 
00981      Teuchos::RCP<QCAD::CoupledPSPreconditioner> WPrec_out_psp = Teuchos::rcp_dynamic_cast<QCAD::CoupledPSPreconditioner>(WPrec_out, true);
00982      WPrec_out_psp->initialize(WPrec_poisson, WPrec_schrodinger);
00983   }
00984 
00985   // df/dp
00986   for (int i=0; i<outArgs.Np(); i++) {
00987     Teuchos::RCP<Epetra_MultiVector> dfdp_out = 
00988       outArgs.get_DfDp(i).getMultiVector();
00989     if (dfdp_out != Teuchos::null) {
00990 
00991 
00992       // Get views into dfdp_out vectors for poisson and schrodinger parts
00993       //  Note that df/dp will be zero for parts of f corresponding to an app
00994       //    different from the one owning the p vector.  E.g. if i==0 corresponds
00995       //    to p being a poisson parameter vector then df/dp == 0 for all the schrodinger
00996       //    parts of f.
00997 
00998       int nParamComponents = dfdp_out->NumVectors();
00999       Teuchos::RCP<Epetra_MultiVector> dfdp_poisson;
01000       std::vector< Teuchos::RCP<Epetra_MultiVector> > dfdp_schrodinger(nEigenvals);
01001 
01002       double *dfdp_data; int myLDA;
01003       if(dfdp_out->ExtractView(&dfdp_data, &myLDA) != 0) 
01004   TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
01005            "Error!  QCAD::CoupledPoissonSchrodinger -- cannot extract dgdp vector view");
01006 
01007       dfdp_poisson = Teuchos::rcp(new Epetra_MultiVector(::View, *disc_map, &(dfdp_data[0]), myLDA, nParamComponents));
01008       for(int k=0; k<nEigenvals; k++)
01009   dfdp_schrodinger[k] = Teuchos::rcp(new Epetra_MultiVector(::View, *disc_map, &(dfdp_data[(k+1)*disc_nMyElements]), myLDA, nParamComponents));
01010 
01011       // Assemble p_vec
01012       Teuchos::Array<int> p_indexes = 
01013         outArgs.get_DfDp(i).getDerivativeMultiVector().getParamIndexes();
01014       Teuchos::RCP<ParamVec> p_vec;
01015 
01016       Teuchos::Array<ParamVec>& sacado_param_vec = 
01017   (i < num_poisson_param_vecs) ? poisson_sacado_param_vec : schrodinger_sacado_param_vec;
01018       int offset = (i < num_poisson_param_vecs) ? 0 : num_poisson_param_vecs;
01019 
01020       if (p_indexes.size() == 0)
01021         p_vec = Teuchos::rcp(&sacado_param_vec[i-offset],false);
01022       else {
01023         p_vec = Teuchos::rcp(new ParamVec);
01024         for (int j=0; j<p_indexes.size(); j++)
01025           p_vec->addParam(sacado_param_vec[i-offset][p_indexes[j]].family, 
01026                           sacado_param_vec[i-offset][p_indexes[j]].baseValue);
01027       }
01028 
01029       dfdp_out->PutScalar(0.0);
01030 
01031       // Compute full dfdp by computing non-zero parts and leaving zeros in others
01032       if (i < num_poisson_param_vecs) {
01033   // "Poisson-owned" param vector, so only poisson part of dfdp vector can be nonzero
01034   poissonApp->computeGlobalTangent(0.0, 0.0, 0.0, curr_time, false, xdot_poisson.get(), NULL, *x_poisson, 
01035           poisson_sacado_param_vec, p_vec.get(),
01036           NULL, NULL, NULL, NULL, f_poisson.get(), NULL, 
01037           dfdp_poisson.get());
01038 
01039   f_poisson_already_computed=true;
01040       }
01041       else {
01042   // "Schrodinger-owned" param vector, so only schrodinger parts of dfdp vector can be nonzero
01043   for(int k=0; k<nEigenvals; k++) {
01044     schrodingerApp->computeGlobalTangent(0.0, 0.0, 0.0, curr_time, false, xdot_schrodinger_vec[k], NULL, *((*x_schrodinger)(k)),
01045                  schrodinger_sacado_param_vec, p_vec.get(),
01046                  NULL, NULL, NULL, NULL, f_schrodinger_vec[k], NULL, 
01047                  dfdp_schrodinger[k].get());  
01048     f_schrodinger_already_computed[k]=true;
01049   }
01050       }
01051     }
01052   }
01053 
01054   // f
01055   /*if (app->is_adjoint) {  //TODO: support Adjoints?
01056     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
01057            "Error!  QCAD::CoupledPoissonSchrodinger -- adjoints not implemented yet");
01058     Derivative f_deriv(f_out, DERIV_TRANS_MV_BY_ROW);
01059     int response_index = 0; // need to add capability for sending this in
01060     app->evaluateResponseDerivative(response_index, curr_time, x_dot.get(), *x, 
01061             sacado_param_vec, NULL, 
01062             NULL, f_deriv, Derivative(), Derivative());
01063   }
01064   else {  */
01065     if (f_out != Teuchos::null) { 
01066       Epetra_Vector M_vec(*disc_map);  //temp storage for mass matrix times vec -- maybe don't allocate this on the stack??
01067 
01068       if(!f_poisson_already_computed) {
01069   poissonApp->computeGlobalResidual(curr_time, xdot_poisson.get(), NULL, *x_poisson, 
01070             poisson_sacado_param_vec, *f_poisson);
01071       }
01072       
01073       for(int i=0; i<nEigenvals; i++) {
01074 
01075   // Compute Mass_matrix * eigenvector[i]
01076   const Epetra_Vector& vec = *((*x_schrodinger)(i));
01077   M_out_schrodinger_crs->Multiply(false, vec, M_vec);  
01078 
01079 
01080   // Compute the schrodinger residual f_schrodinger_vec[i]: H*eigenvector[i] - eigenvalue[i] * M * eigenvector[i]
01081 
01082   if(!f_schrodinger_already_computed[i]) {
01083 
01084       //Could call this, but multiply below is faster
01085     /*schrodingerApp->computeGlobalResidual(curr_time, xdot_schrodinger_vec[i], *((*x_schrodinger)(i)), 
01086                   schrodinger_sacado_param_vec, *(f_schrodinger_vec[i]) );  // H*evec[i] */
01087 
01088     // H * Psi - E * M * Psi
01089     const Epetra_CrsMatrix& Hamiltonian_crs =  *J_out_schrodinger_crs;
01090     Hamiltonian_crs.Multiply(false, vec, *(f_schrodinger_vec[i]));
01091   }
01092 
01093   /* ---- DEBUG ----
01094      double He_norm, Me_norm, H_expect;
01095      f_schrodinger_vec[i]->Norm2(&He_norm);
01096      M_vec.Norm2(&Me_norm);
01097      std::cout << "DEBUG " << i << ": norm(-H*evec) = " << He_norm << ", norm(M*evec) = " << Me_norm 
01098      << ", eval = " << (*stdvec_eigenvals)[i] << std::endl;
01099      f_schrodinger_vec[i]->Dot(vec, &H_expect);
01100   */
01101 
01102   // add -eval[i]*M*evec[i] to H*evec[i] (recall evals are really negative_evals)
01103   f_schrodinger_vec[i]->Update( (*stdvec_eigenvals)[i], M_vec, 1.0); 
01104 
01105 
01106         // Compute normalization equation residuals:  f_norm[i] = abs(1 - evec[i] . M . evec[i])
01107   double vec_M_vec;
01108   vec.Dot( M_vec, &vec_M_vec );
01109   (*f_norm_local)[i] = 1.0 - vec_M_vec;
01110       }
01111 
01112       // Fill elements of f_norm_dist that belong to this processor, i.e. loop over
01113       // eigenvalue indices "owned" by the current proc in the combined distributed map
01114       std::vector<int> eval_global_elements(my_nEigenvals);
01115       dist_eigenval_map.MyGlobalElements(&eval_global_elements[0]);
01116       for(int i=0; i<my_nEigenvals; i++)
01117   (*f_norm_dist)[i] = (*f_norm_local)[eval_global_elements[i]];
01118 
01119       
01120 
01121       //DEBUG -- print residual in gory detail for debugging
01122       if(1) {
01123   if(myComm->MyPID() == 0) std::cout << "DEBUG: ----------------- Coupled Schrodinger Poisson Info Dump ---------------------" << std::endl;
01124   double norm, mean;
01125 
01126   /*std::cout << "x map has " << x->Map().NumGlobalElements() << " global els" << std::endl;
01127     std::cout << "x_poisson map has " << x_poisson->Map().NumGlobalElements() << " global els" << std::endl;
01128     std::cout << "x_schrodinger map has " << x_schrodinger->Map().NumGlobalElements() << " global els (each vec)" << std::endl;
01129     std::cout << "dist_eval_map has " << dist_eigenval_map.NumGlobalElements() << " global els" << std::endl;
01130   */
01131 
01132   x->Norm2(&norm); x->MeanValue(&mean);
01133   std::cout << std::setprecision(10);
01134   if(myComm->MyPID() == 0) std::cout << "X Norm & Mean = " << norm << " , " << mean << std::endl;
01135   
01136   x_poisson->Norm2(&norm); x_poisson->MeanValue(&mean);
01137   if(myComm->MyPID() == 0) std::cout << "Poisson-part X Norm & Mean = " << norm << " , " << mean << std::endl;
01138   for(int i=0; i<nEigenvals; i++) {
01139     (*x_schrodinger)(i)->Norm2(&norm);
01140     if(myComm->MyPID() == 0) std::cout << "Schrodinger[" << i << "]-part X Norm = " << norm << std::endl;
01141   }
01142   for(int i=0; i<nEigenvals; i++) 
01143     if(myComm->MyPID() == 0) std::cout << "Eigenvalue[" << i << "] = " << (*stdvec_eigenvals)[i] << std::endl;
01144   
01145   f_poisson->Norm2(&norm);
01146   if(myComm->MyPID() == 0) std::cout << "Poisson-part Residual Norm = " << norm << std::endl; //f_poisson->Print(std::cout);
01147   for(int i=0; i<nEigenvals; i++) {
01148     if(f_schrodinger_vec[i] != NULL) {
01149       f_schrodinger_vec[i]->Norm2(&norm);
01150       if(myComm->MyPID() == 0) std::cout << "Schrodinger[" << i << "]-part Residual Norm = " << norm << std::endl; //f_schrodinger_vec[i]->Print(std::cout);
01151     }
01152   }
01153   if(myComm->MyPID() == 0) std::cout << "Eigenvalue-part Residual: " << std::endl;
01154   f_norm_dist->Print(std::cout); // only rank 0 prints
01155       }
01156     }
01157 
01158   // Response functions
01159   for (int i=0; i<outArgs.Ng(); i++) {
01160     Teuchos::RCP<Epetra_Vector> g_out = outArgs.get_g(i);
01161    
01162     bool g_computed = false;
01163 
01164     Derivative dgdx_out = outArgs.get_DgDx(i);
01165     Derivative dgdxdot_out = outArgs.get_DgDx_dot(i);
01166 
01167     // dg/dx, dg/dxdot
01168     if (!dgdx_out.isEmpty() || !dgdxdot_out.isEmpty()) {
01169       if(i < poissonApp->getNumResponses()) {
01170   poissonApp->evaluateResponseDerivative(i, curr_time, xdot_poisson.get(), NULL, *(x_poisson),
01171                                       poisson_sacado_param_vec, NULL,
01172                                       g_out.get(), dgdx_out,
01173                                       dgdxdot_out, Derivative(), Derivative());
01174       }
01175       else {
01176   // take response derivatives using lowest eigenstate only (is there something better??)
01177   schrodingerApp->evaluateResponseDerivative(i - poissonApp->getNumResponses(), curr_time, xdot_schrodinger_vec[0], NULL,
01178                                       *((*x_schrodinger)(0)),
01179                                       schrodinger_sacado_param_vec, NULL,
01180                                       g_out.get(), dgdx_out,
01181                                       dgdxdot_out, Derivative(), Derivative());
01182       }
01183       g_computed = true;
01184     }
01185 
01186     // dg/dp
01187     for (int j=0; j<outArgs.Np(); j++) {
01188       Teuchos::RCP<Epetra_MultiVector> dgdp_out =
01189         outArgs.get_DgDp(i,j).getMultiVector();
01190       if (dgdp_out != Teuchos::null) {
01191         Teuchos::Array<int> p_indexes =
01192           outArgs.get_DgDp(i,j).getDerivativeMultiVector().getParamIndexes();
01193 
01194         Teuchos::RCP<ParamVec> p_vec;
01195 
01196   Teuchos::Array<ParamVec>& sacado_param_vec = 
01197     (j < num_poisson_param_vecs) ? poisson_sacado_param_vec : schrodinger_sacado_param_vec;
01198   int offset = (j < num_poisson_param_vecs) ? 0 : num_poisson_param_vecs;
01199 
01200         if (p_indexes.size() == 0)
01201           p_vec = Teuchos::rcp(&sacado_param_vec[j-offset],false);
01202         else {
01203           p_vec = Teuchos::rcp(new ParamVec);
01204           for (int k=0; k<p_indexes.size(); k++)
01205             p_vec->addParam(sacado_param_vec[j-offset][p_indexes[k]].family,
01206                             sacado_param_vec[j-offset][p_indexes[k]].baseValue);
01207         }
01208 
01209   if(i < poissonApp->getNumResponses() && j < num_poisson_param_vecs) {
01210     //both response and param vectors belong to poisson problem
01211     poissonApp->evaluateResponseTangent(i, alpha, beta, 0.0, curr_time, false,
01212                 xdot_poisson.get(), NULL, *x_poisson,
01213                 poisson_sacado_param_vec, p_vec.get(),
01214                 NULL, NULL, NULL, NULL, g_out.get(), NULL,
01215                 dgdp_out.get());
01216   }
01217   else if(i >= poissonApp->getNumResponses() && j >= num_poisson_param_vecs) {
01218     //both response and param vectors belong to schrodinger problem -- evaluate dg/dp using first eigenvector
01219     schrodingerApp->evaluateResponseTangent(i - poissonApp->getNumResponses(), alpha, beta, 0.0, curr_time, false,
01220               xdot_schrodinger_vec[0], NULL, *((*x_schrodinger)(0)),
01221               schrodinger_sacado_param_vec, p_vec.get(),
01222               NULL, NULL, NULL, NULL, g_out.get(), NULL,
01223               dgdp_out.get());
01224   }
01225   else {
01226     // response and param vectors belong to different sub-problems (Poisson or Schrodinger)
01227     dgdp_out->PutScalar(0.0);
01228   }
01229 
01230         g_computed = true;
01231       }
01232     }
01233 
01234     if (g_out != Teuchos::null && !g_computed) {
01235       if(i < poissonApp->getNumResponses()) {
01236   poissonApp->evaluateResponse(i, curr_time, xdot_poisson.get(), NULL, *x_poisson,
01237              poisson_sacado_param_vec, *g_out);
01238       }
01239       else {
01240   schrodingerApp->evaluateResponse(i - poissonApp->getNumResponses(), curr_time, xdot_schrodinger_vec[0], NULL, *((*x_schrodinger)(0)), 
01241            schrodinger_sacado_param_vec, *g_out);
01242       }
01243     }
01244 
01245   }
01246 
01247 }
01248 
01249 Teuchos::RCP<Albany::Application>
01250 QCAD::CoupledPoissonSchrodinger::getPoissonApp() const
01251 {
01252   return poissonApp;
01253 }
01254 
01255 Teuchos::RCP<Albany::Application>
01256 QCAD::CoupledPoissonSchrodinger::getSchrodingerApp() const
01257 {
01258   return schrodingerApp;
01259 }
01260 
01261 
01262 
01263 // Note: we could use QCAD::separateCombinedVector(...) to implement this function and those below
01264 void QCAD::CoupledPoissonSchrodinger::separateCombinedVector(const Teuchos::RCP<Epetra_Vector>& combinedVector,
01265                    Teuchos::RCP<Epetra_Vector>& poisson_part,
01266                    Teuchos::RCP<Epetra_MultiVector>& schrodinger_part) const
01267 {
01268   double* data;
01269   int disc_nMyElements = disc_map->NumMyElements();
01270 
01271   if(combinedVector->ExtractView(&data) != 0)
01272     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
01273              "Error!  QCAD::CoupledPoissonSchrodinger cannot extract vector views");
01274 
01275   poisson_part = Teuchos::rcp(new Epetra_Vector(::View, *disc_map, &data[0]));
01276   schrodinger_part = Teuchos::rcp(new Epetra_MultiVector(::View, *disc_map, &data[disc_nMyElements], disc_nMyElements, nEigenvals));
01277 }
01278 
01279 
01280 void QCAD::CoupledPoissonSchrodinger::separateCombinedVector(const Teuchos::RCP<Epetra_Vector>& combinedVector,
01281                    Teuchos::RCP<Epetra_Vector>& poisson_part,
01282                    Teuchos::RCP<Epetra_MultiVector>& schrodinger_part,
01283                    Teuchos::RCP<Epetra_Vector>& eigenvalue_part) const
01284 {
01285   this->separateCombinedVector(combinedVector, poisson_part, schrodinger_part);
01286 
01287   double* data;
01288   int disc_nMyElements = disc_map->NumMyElements();
01289   int my_nEigenvals = combined_SP_map->NumMyElements() - disc_nMyElements * (1+nEigenvals);
01290   Epetra_Map dist_eigenval_map(nEigenvals, my_nEigenvals, 0, *myComm);
01291 
01292   combinedVector->ExtractView(&data); //above call tests for failure
01293   eigenvalue_part = Teuchos::rcp(new Epetra_Vector(::View, dist_eigenval_map, &data[(1+nEigenvals)*disc_nMyElements]));
01294 }
01295 
01296 
01297 
01298 void QCAD::CoupledPoissonSchrodinger::separateCombinedVector(const Teuchos::RCP<const Epetra_Vector>& combinedVector,
01299                    Teuchos::RCP<const Epetra_Vector>& poisson_part,
01300                    Teuchos::RCP<const Epetra_MultiVector>& schrodinger_part) const
01301 {
01302   double* data;
01303   int disc_nMyElements = disc_map->NumMyElements();
01304 
01305   if(combinedVector->ExtractView(&data) != 0)
01306     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
01307              "Error!  QCAD::CoupledPoissonSchrodinger cannot extract vector views");
01308 
01309   poisson_part = Teuchos::rcp(new const Epetra_Vector(::View, *disc_map, &data[0]));
01310   schrodinger_part = Teuchos::rcp(new const Epetra_MultiVector(::View, *disc_map, &data[disc_nMyElements], disc_nMyElements, nEigenvals));
01311 }
01312 
01313 
01314 void QCAD::CoupledPoissonSchrodinger::separateCombinedVector(const Teuchos::RCP<const Epetra_Vector>& combinedVector,
01315                    Teuchos::RCP<const Epetra_Vector>& poisson_part,
01316                    Teuchos::RCP<const Epetra_MultiVector>& schrodinger_part,
01317                    Teuchos::RCP<const Epetra_Vector>& eigenvalue_part) const
01318 {
01319   this->separateCombinedVector(combinedVector, poisson_part, schrodinger_part);
01320 
01321   double* data;
01322   int disc_nMyElements = disc_map->NumMyElements();
01323   int my_nEigenvals = combined_SP_map->NumMyElements() - disc_nMyElements * (1+nEigenvals);
01324   Epetra_Map dist_eigenval_map(nEigenvals, my_nEigenvals, 0, *myComm);
01325 
01326   combinedVector->ExtractView(&data); //above call tests for failure
01327   eigenvalue_part = Teuchos::rcp(new const Epetra_Vector(::View, dist_eigenval_map, &data[(1+nEigenvals)*disc_nMyElements]));
01328 }
01329 
01330 
01331 
01332 //Copied from Albany::SolverFactory -- used to validate applicaton parameters of applications not created via a SolverFactory
01333 Teuchos::RCP<const Teuchos::ParameterList>
01334 QCAD::CoupledPoissonSchrodinger::getValidAppParameters() const
01335 {  
01336   Teuchos::RCP<Teuchos::ParameterList> validPL = Teuchos::rcp(new Teuchos::ParameterList("ValidAppParams"));;
01337   validPL->sublist("Problem",            false, "Problem sublist");
01338   validPL->sublist("Debug Output",       false, "Debug Output sublist");
01339   validPL->sublist("Discretization",     false, "Discretization sublist");
01340   validPL->sublist("Quadrature",         false, "Quadrature sublist");
01341   validPL->sublist("Regression Results", false, "Regression Results sublist");
01342   validPL->sublist("VTK",                false, "DEPRECATED  VTK sublist");
01343   validPL->sublist("Piro",               false, "Piro sublist");
01344   validPL->sublist("Coupled System",     false, "Coupled system sublist");
01345 
01346   return validPL;
01347 }
01348 
01349 Teuchos::RCP<const Teuchos::ParameterList>
01350 QCAD::CoupledPoissonSchrodinger::getValidProblemParameters() const
01351 {
01352   Teuchos::RCP<Teuchos::ParameterList> validPL = Teuchos::createParameterList("ValidPoissonSchrodingerProblemParams");
01353 
01354   validPL->set<std::string>("Name", "", "String to designate Problem Class");
01355   validPL->set<int>("Phalanx Graph Visualization Detail", 0,
01356                     "Flag to select output of Phalanx Graph and level of detail");
01357 
01358   validPL->set<double>("Length Unit In Meters",1e-6,"Length unit in meters");
01359   validPL->set<double>("Energy Unit In Electron Volts",1.0,"Energy (voltage) unit in electron volts");
01360   validPL->set<double>("Temperature",300,"Temperature in Kelvin");
01361   validPL->set<std::string>("MaterialDB Filename","materials.xml","Filename of material database xml file");
01362   validPL->set<int>("Number of Eigenvalues",0,"The number of eigenvalue-eigenvector pairs");
01363   validPL->set<bool>("Verbose Output",false,"Enable detailed output mode");
01364 
01365   validPL->set<bool>("Include exchange-correlation potential",false,"Include exchange-correlation potential in poisson source term");
01366   validPL->set<bool>("Only solve schrodinger in quantum blocks",true,"Limit schrodinger solution to elements blocks labeled as quantum in the materials DB");
01367 
01368   validPL->sublist("Poisson Problem", false, "");
01369   validPL->sublist("Schrodinger Problem", false, "");
01370 
01371   // Candidates for deprecation. Pertain to the solution rather than the problem definition.
01372   validPL->set<std::string>("Solution Method", "Steady", "Flag for Steady, Transient, or Continuation");
01373   
01374   return validPL;
01375 }
01376 
01377 
01378 //This function is used solely for Jacobian debugging
01379 /*void QCAD::CoupledPoissonSchrodinger::computeResidual(const Teuchos::RCP<const Epetra_Vector>& x,
01380                   Teuchos::RCP<Epetra_Vector>& f,
01381                   Teuchos::RCP<Epetra_CrsMatrix>& massMx) const
01382 {
01383   double curr_time = 0.0;
01384   Epetra_Vector M_vec(*disc_map);  //temp storage for mass matrix times vec -- maybe don't allocate this on the stack??
01385   Epetra_LocalMap local_eigenval_map(nEigenvals, 0, *myComm);
01386 
01387   int disc_nMyElements = disc_map->NumMyElements();
01388   int my_nEigenvals = combined_SP_map->NumMyElements() - disc_nMyElements * (1+nEigenvals);
01389 
01390   Teuchos::RCP<const Epetra_Vector> x_poisson, eigenvals_dist;
01391   Teuchos::RCP<const Epetra_MultiVector> x_schrodinger;
01392   separateCombinedVector(x, x_poisson, x_schrodinger, eigenvals_dist);
01393 
01394   Teuchos::RCP<Epetra_Vector>   f_poisson, f_norm_local, f_norm_dist;
01395   Teuchos::RCP<Epetra_MultiVector> f_schrodinger;
01396   std::vector<Epetra_Vector*> f_schrodinger_vec(nEigenvals);
01397   separateCombinedVector(f, f_poisson, f_schrodinger, f_norm_dist);
01398   for(int i=0; i<nEigenvals; i++) f_schrodinger_vec[i] = (*f_schrodinger)(i);
01399   f_norm_local = Teuchos::rcp(new Epetra_Vector(local_eigenval_map));
01400 
01401 
01402   //update schrodinger wavefunctions for poisson
01403 
01404   Epetra_Import eigenval_importer(local_eigenval_map, eigenvals_dist->Map() );
01405   Teuchos::RCP<Epetra_Vector> eigenvals =  Teuchos::rcp(new Epetra_Vector(local_eigenval_map));
01406   eigenvals->Import(*eigenvals_dist, eigenval_importer, Insert);
01407   Teuchos::RCP<std::vector<double> > stdvec_eigenvals = Teuchos::rcp(new std::vector<double>(&(*eigenvals)[0], &(*eigenvals)[0] + nEigenvals));
01408 
01409   Teuchos::RCP<Albany::EigendataStruct> eigenData = Teuchos::rcp( new Albany::EigendataStruct );
01410   eigenData->eigenvalueRe = stdvec_eigenvals;
01411   eigenData->eigenvectorRe = Teuchos::rcp(new Epetra_MultiVector(*disc_overlap_map, nEigenvals));
01412   eigenData->eigenvectorIm = Teuchos::null;
01413   Teuchos::RCP<Epetra_Import> overlap_importer = Teuchos::rcp(new Epetra_Import(*disc_overlap_map, *disc_map));
01414 
01415   for(int i=0; i<nEigenvals; i++)
01416     (*(eigenData->eigenvectorRe))(i)->Import( *((*x_schrodinger)(i)), *overlap_importer, Insert );
01417 
01418     // set eigenvalues / eigenvectors for use in poisson problem:
01419   poissonApp->getStateMgr().setEigenData(eigenData);
01420   poissonApp->computeGlobalResidual(curr_time, NULL, *x_poisson, 
01421             poisson_sacado_param_vec, *f_poisson);
01422       
01423 
01424     // Get overlapped version of potential (x_poisson) for passing as auxData to schrodinger app
01425   Teuchos::RCP<Epetra_MultiVector> overlapped_V = Teuchos::rcp(new Epetra_MultiVector(*disc_overlap_map, 1));
01426   Teuchos::RCP<Epetra_Vector> ones_vec = Teuchos::rcp(new Epetra_Vector(*disc_overlap_map));
01427   ones_vec->PutScalar(1.0);
01428   (*overlapped_V)(0)->Import( *x_poisson, *overlap_importer, Insert );
01429   (*overlapped_V)(0)->Update(offset_to_CB, *ones_vec, -1.0);
01430   schrodingerApp->getStateMgr().setAuxData(overlapped_V);
01431 
01432     // compute schrodinger Hamiltonian
01433   Teuchos::RCP<Epetra_Operator> hamMx = schrodingerModel->create_W(); //maybe re-use this and not create it every time?
01434   Teuchos::RCP<Epetra_CrsMatrix> hamMx_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(hamMx, true);
01435   Teuchos::RCP<const Epetra_Vector> dummy_xdot = schrodingerModel->get_x_dot_init();
01436   schrodingerApp->computeGlobalJacobian(0.0, 1.0, 0.0, curr_time, dummy_xdot.get(), NULL, *((*x_schrodinger)(0)), 
01437           schrodinger_sacado_param_vec, f_schrodinger_vec[0], *hamMx_crs);
01438 
01439   for(int i=0; i<nEigenvals; i++) {
01440     const Epetra_Vector& vec = *((*x_schrodinger)(i));
01441     massMx->Multiply(false, vec, M_vec);  
01442     hamMx_crs->Multiply(false, vec, *(f_schrodinger_vec[i]));  
01443     f_schrodinger_vec[i]->Update( (*stdvec_eigenvals)[i], M_vec, 1.0); 
01444 
01445     // Compute normalization equation residuals:  f_norm[i] = abs(1 - evec[i] . M . evec[i])
01446     double vec_M_vec;
01447     vec.Dot( M_vec, &vec_M_vec );
01448     (*f_norm_local)[i] = 1.0 - vec_M_vec;
01449   }
01450 
01451   // Fill elements of f_norm_dist that belong to this processor, i.e. loop over
01452   // eigenvalue indices "owned" by the current proc in the combined distributed map
01453   std::vector<int> eval_global_elements(my_nEigenvals);
01454   eigenvals_dist->Map().MyGlobalElements(&eval_global_elements[0]);
01455   for(int i=0; i<my_nEigenvals; i++)
01456     (*f_norm_dist)[i] = (*f_norm_local)[eval_global_elements[i]];
01457 }*/