Albany: a Trilinos-based PDE code: Albany_SolverFactory.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 "Albany_SolverFactory.hpp"
00008 #include "Albany_ObserverFactory.hpp"
00009 #include "Albany_PiroObserver.hpp"
00010 #include "Albany_SaveEigenData.hpp"
00011 #include "Albany_ModelFactory.hpp"
00012 
00013 #include "Piro_Epetra_SolverFactory.hpp"
00014 #include "Piro_ProviderBase.hpp"
00015 
00016 #include "Piro_NOXSolver.hpp"
00017 #include "Piro_NullSpaceUtils.hpp"
00018 #include "Piro_StratimikosUtils.hpp"
00019 
00020 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00021 
00022 #ifdef ALBANY_QCAD
00023   #include "QCAD_Solver.hpp"
00024   #include "QCAD_CoupledPoissonSchrodinger.hpp"
00025   #include "QCAD_CoupledPSObserver.hpp"
00026   #include "QCAD_GenEigensolver.hpp"
00027 #endif
00028 
00029 //#include "Thyra_EpetraModelEvaluator.hpp"
00030 //#include "AAdapt_AdaptiveModelFactory.hpp"
00031 
00032 #include "Thyra_DetachedVectorView.hpp"
00033 
00034 #include "Teuchos_XMLParameterListHelpers.hpp"
00035 #include "Teuchos_TestForException.hpp"
00036 
00037 #include "NOX_Epetra_Observer.H"
00038 #include "Rythmos_IntegrationObserverBase.hpp"
00039 
00040 #include "Albany_Application.hpp"
00041 #include "Albany_Utils.hpp"
00042 
00043 //#include <stdexcept>
00044 
00045 
00046 namespace Albany {
00047 
00048 class NOXObserverConstructor : public Piro::ProviderBase<NOX::Epetra::Observer> {
00049 public:
00050   explicit NOXObserverConstructor(const Teuchos::RCP<Application> &app) :
00051     factory_(app),
00052     instance_(Teuchos::null)
00053   {}
00054 
00055   virtual Teuchos::RCP<NOX::Epetra::Observer> getInstance(
00056       const Teuchos::RCP<Teuchos::ParameterList> &params);
00057 
00058 private:
00059   NOXObserverFactory factory_;
00060   Teuchos::RCP<NOX::Epetra::Observer> instance_;
00061 };
00062 
00063 Teuchos::RCP<NOX::Epetra::Observer>
00064 NOXObserverConstructor::getInstance(const Teuchos::RCP<Teuchos::ParameterList> &/*params*/)
00065 {
00066   if (Teuchos::is_null(instance_)) {
00067     instance_ = factory_.createInstance();
00068   }
00069   return instance_;
00070 }
00071 
00072 class RythmosObserverConstructor : public Piro::ProviderBase<Rythmos::IntegrationObserverBase<double> > {
00073 public:
00074   explicit RythmosObserverConstructor(const Teuchos::RCP<Application> &app) :
00075     factory_(app),
00076     instance_(Teuchos::null)
00077   {}
00078 
00079   virtual Teuchos::RCP<Rythmos::IntegrationObserverBase<double> > getInstance(
00080       const Teuchos::RCP<Teuchos::ParameterList> &params);
00081 
00082 private:
00083   RythmosObserverFactory factory_;
00084   Teuchos::RCP<Rythmos::IntegrationObserverBase<double> > instance_;
00085 };
00086 
00087 Teuchos::RCP<Rythmos::IntegrationObserverBase<double> >
00088 RythmosObserverConstructor::getInstance(const Teuchos::RCP<Teuchos::ParameterList> &/*params*/)
00089 {
00090   if (Teuchos::is_null(instance_)) {
00091     instance_ = factory_.createInstance();
00092   }
00093   return instance_;
00094 }
00095 
00096 class SaveEigenDataConstructor : public Piro::ProviderBase<LOCA::SaveEigenData::AbstractStrategy> {
00097 public:
00098   SaveEigenDataConstructor(
00099       Teuchos::ParameterList &locaParams,
00100       StateManager* pStateMgr,
00101       const Teuchos::RCP<Piro::ProviderBase<NOX::Epetra::Observer> > &observerProvider) :
00102     locaParams_(locaParams),
00103     pStateMgr_(pStateMgr),
00104     observerProvider_(observerProvider)
00105   {}
00106 
00107   virtual Teuchos::RCP<LOCA::SaveEigenData::AbstractStrategy> getInstance(
00108       const Teuchos::RCP<Teuchos::ParameterList> &params);
00109 
00110 private:
00111   Teuchos::ParameterList &locaParams_;
00112   StateManager* pStateMgr_;
00113 
00114   Teuchos::RCP<Piro::ProviderBase<NOX::Epetra::Observer> > observerProvider_;
00115 };
00116 
00117 Teuchos::RCP<LOCA::SaveEigenData::AbstractStrategy>
00118 SaveEigenDataConstructor::getInstance(const Teuchos::RCP<Teuchos::ParameterList> &params)
00119 {
00120   const Teuchos::RCP<NOX::Epetra::Observer> noxObserver = observerProvider_->getInstance(params);
00121   return Teuchos::rcp(new SaveEigenData(locaParams_, noxObserver, pStateMgr_));
00122 }
00123 
00124 } // namespace Albany
00125 
00126 using Teuchos::RCP;
00127 using Teuchos::rcp;
00128 using Teuchos::ParameterList;
00129 
00130 
00131 
00132 Albany::SolverFactory::SolverFactory(
00133         const std::string& inputFile,
00134         const Albany_MPI_Comm& mcomm)
00135   : out(Teuchos::VerboseObjectBase::getDefaultOStream())
00136 {
00137   RCP<Teuchos::Comm<int> > tcomm = Albany::createTeuchosCommFromMpiComm(mcomm);
00138 
00139   // Set up application parameters: read and broadcast XML file, and set defaults
00140   //RCP<ParameterList> input_
00141   appParams = Teuchos::createParameterList("Albany Parameters");
00142   Teuchos::updateParametersFromXmlFileAndBroadcast(inputFile, appParams.ptr(), *tcomm);
00143 
00144   //do not set default solver parameters for QCAD::Solver problems, as it handles this itself
00145   if (appParams->sublist("Problem").get("Solution Method", "Steady") != "QCAD Multi-Problem") {  
00146     RCP<ParameterList> defaultSolverParams = rcp(new ParameterList());
00147     setSolverParamDefaults(defaultSolverParams.get(), tcomm->getRank());
00148     appParams->setParametersNotAlreadySet(*defaultSolverParams);
00149   }
00150 
00151   appParams->validateParametersAndSetDefaults(*getValidAppParameters(),0);
00152 }
00153   
00154 
00155 Albany::SolverFactory::SolverFactory(
00156                   const Teuchos::RCP<Teuchos::ParameterList>& input_appParams,
00157         const Albany_MPI_Comm& mcomm)
00158   : appParams(input_appParams), out(Teuchos::VerboseObjectBase::getDefaultOStream())
00159 {
00160   RCP<Teuchos::Comm<int> > tcomm = Albany::createTeuchosCommFromMpiComm(mcomm);
00161 
00162   //do not set default solver parameters for QCAD::Solver problems, as it handles this itself
00163   if (appParams->sublist("Problem").get("Solution Method", "Steady") != "QCAD Multi-Problem") {  
00164     RCP<ParameterList> defaultSolverParams = rcp(new ParameterList());
00165     setSolverParamDefaults(defaultSolverParams.get(), tcomm->getRank());
00166     appParams->setParametersNotAlreadySet(*defaultSolverParams);
00167   }
00168 
00169   appParams->validateParametersAndSetDefaults(*getValidAppParameters(),0);
00170 }
00171 
00172 Albany::SolverFactory::~SolverFactory(){
00173   
00174   // Release the model to eliminate RCP circular reference
00175   if(Teuchos::nonnull(thyraModelFactory))
00176     thyraModelFactory->releaseModel();
00177 
00178 #ifdef ALBANY_DEBUG
00179   *out << "Calling destructor for Albany_SolverFactory" << std::endl;
00180 #endif
00181 }
00182 
00183 
00184 Teuchos::RCP<EpetraExt::ModelEvaluator>
00185 Albany::SolverFactory::create(
00186   const Teuchos::RCP<const Epetra_Comm>& appComm,
00187   const Teuchos::RCP<const Epetra_Comm>& solverComm,
00188   const Teuchos::RCP<const Epetra_Vector>& initial_guess)
00189 {
00190   Teuchos::RCP<Albany::Application> dummyAlbanyApp;
00191   return createAndGetAlbanyApp(dummyAlbanyApp, appComm, solverComm, initial_guess);
00192 }
00193 
00194 Teuchos::RCP<EpetraExt::ModelEvaluator>
00195 Albany::SolverFactory::createAndGetAlbanyApp(
00196   Teuchos::RCP<Albany::Application>& albanyApp,
00197   const Teuchos::RCP<const Epetra_Comm>& appComm,
00198   const Teuchos::RCP<const Epetra_Comm>& solverComm,
00199   const Teuchos::RCP<const Epetra_Vector>& initial_guess)
00200 {
00201     const RCP<ParameterList> problemParams = Teuchos::sublist(appParams, "Problem");
00202     const std::string solutionMethod = problemParams->get("Solution Method", "Steady");
00203 
00204     if (solutionMethod == "QCAD Multi-Problem") {
00205 #ifdef ALBANY_QCAD
00206       return rcp(new QCAD::Solver(appParams, solverComm, initial_guess));
00207 #else /* ALBANY_QCAD */
00208       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Must activate QCAD\n");
00209 #endif /* ALBANY_QCAD */
00210     }
00211 
00212     if (solutionMethod == "QCAD Poisson-Schrodinger") {
00213 #ifdef ALBANY_QCAD
00214       const RCP<QCAD::CoupledPoissonSchrodinger> ps_model = rcp(new QCAD::CoupledPoissonSchrodinger(appParams, solverComm, initial_guess));
00215       const RCP<ParameterList> piroParams = Teuchos::sublist(appParams, "Piro");
00216 
00217       // Create and setup the Piro solver factory
00218       Piro::Epetra::SolverFactory piroFactory;
00219       {
00220         // Do we need: Observers for output from time-stepper ??
00221   const RCP<Piro::ProviderBase<NOX::Epetra::Observer> > noxObserverProvider =
00222     rcp(new QCAD::CoupledPS_NOXObserverConstructor(ps_model));
00223     //  rcp(new NOXObserverConstructor(poisson_app));
00224   piroFactory.setSource<NOX::Epetra::Observer>(noxObserverProvider);
00225 
00226   // LOCA auxiliary objects -- needed?
00227       }
00228       return piroFactory.createSolver(piroParams, ps_model);
00229 
00230 #else /* ALBANY_QCAD */
00231       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Must activate QCAD\n");
00232 #endif /* ALBANY_QCAD */
00233     }
00234 
00235     if (solutionMethod == "Eigensolve") {
00236 #ifdef ALBANY_QCAD
00237 
00238       RCP<Albany::Application> app;
00239       const RCP<EpetraExt::ModelEvaluator> model = createAlbanyAppAndModel(app, appComm, initial_guess);
00240       albanyApp = app;
00241       
00242       //QCAD::GenEigensolver uses a state manager as an observer (for now)
00243       RCP<Albany::StateManager> observer = rcp( &(app->getStateMgr()), false);
00244 
00245       // Currently, QCAD eigensolver just uses LOCA's eigensolver list under Piro -- maybe give it it's own list
00246       //   outside of Piro?
00247       const RCP<ParameterList> eigensolveParams = rcp(&(appParams->sublist("Piro").sublist("LOCA").sublist("Stepper").sublist("Eigensolver")), false);
00248       const RCP<QCAD::GenEigensolver> es_model = rcp(new QCAD::GenEigensolver(eigensolveParams, model, observer, solverComm));
00249       return es_model;
00250 
00251 #else /* ALBANY_QCAD */
00252       TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Must activate QCAD\n");
00253 #endif /* ALBANY_QCAD */
00254     }
00255 
00256     // Solver uses a single app, create it here along with observer
00257     RCP<Albany::Application> app;
00258     const RCP<EpetraExt::ModelEvaluator> model = createAlbanyAppAndModel(app, appComm, initial_guess);
00259 
00260     //Pass back albany app so that interface beyond ModelEvaluator can be used.
00261     // This is essentially a hack to allow additional in/out arguments beyond
00262     //  what ModelEvaluator specifies.
00263     albanyApp = app;
00264 
00265     const RCP<ParameterList> piroParams = Teuchos::sublist(appParams, "Piro");
00266 
00267     if (solutionMethod == "Continuation") {
00268       ParameterList& locaParams = piroParams->sublist("LOCA");
00269       if (app->getDiscretization()->hasRestartSolution()) {
00270         // Pick up problem time from restart file
00271         locaParams.sublist("Stepper").set("Initial Value", app->getDiscretization()->restartDataTime());
00272       }
00273     }
00274 
00275     // Create and setup the Piro solver factory
00276     Piro::Epetra::SolverFactory piroFactory;
00277     {
00278       // Observers for output from time-stepper
00279       const RCP<Piro::ProviderBase<Rythmos::IntegrationObserverBase<double> > > rythmosObserverProvider =
00280         rcp(new RythmosObserverConstructor(app));
00281       piroFactory.setSource<Rythmos::IntegrationObserverBase<double> >(rythmosObserverProvider);
00282 
00283       const RCP<Piro::ProviderBase<NOX::Epetra::Observer> > noxObserverProvider =
00284         rcp(new NOXObserverConstructor(app));
00285       piroFactory.setSource<NOX::Epetra::Observer>(noxObserverProvider);
00286 
00287       // LOCA auxiliary objects
00288       {
00289         const RCP<AAdapt::AdaptiveSolutionManager> adaptMgr = app->getAdaptSolMgr();
00290         piroFactory.setSource<Piro::Epetra::AdaptiveSolutionManager>(adaptMgr);
00291 
00292         const RCP<Piro::ProviderBase<LOCA::SaveEigenData::AbstractStrategy> > saveEigenDataProvider =
00293           rcp(new SaveEigenDataConstructor(piroParams->sublist("LOCA"), &app->getStateMgr(), noxObserverProvider));
00294         piroFactory.setSource<LOCA::SaveEigenData::AbstractStrategy>(saveEigenDataProvider);
00295       }
00296     }
00297 
00298     return piroFactory.createSolver(piroParams, model);
00299 }
00300 
00301 Teuchos::RCP<Thyra::ModelEvaluator<double> >
00302 Albany::SolverFactory::createThyraSolverAndGetAlbanyApp(
00303     Teuchos::RCP<Application>& albanyApp,
00304     const Teuchos::RCP<const Epetra_Comm>& appComm,
00305     const Teuchos::RCP<const Epetra_Comm>& solverComm,
00306     const Teuchos::RCP<const Epetra_Vector>& initial_guess)
00307 {
00308   const RCP<ParameterList> piroParams = Teuchos::sublist(appParams, "Piro");
00309   const Teuchos::Ptr<const std::string> solverToken(piroParams->getPtr<std::string>("Solver Type"));
00310 
00311   const RCP<ParameterList> problemParams = Teuchos::sublist(appParams, "Problem");
00312   const std::string solutionMethod = problemParams->get("Solution Method", "Steady");
00313 
00314   if (Teuchos::nonnull(solverToken) && *solverToken == "ThyraNOX") {
00315     piroParams->set("Solver Type", "NOX");
00316 
00317     RCP<Albany::Application> app;
00318 
00319     // Creates the Albany::ModelEvaluator
00320     const RCP<EpetraExt::ModelEvaluator> model = createAlbanyAppAndModel(app, appComm, initial_guess);
00321 
00322     // Pass back albany app so that interface beyond ModelEvaluator can be used.
00323     // This is essentially a hack to allow additional in/out arguments beyond
00324     // what ModelEvaluator specifies.
00325     albanyApp = app;
00326 
00327     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00328     linearSolverBuilder.setParameterList(Piro::extractStratimikosParams(piroParams));
00329 
00330     const RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory =
00331       createLinearSolveStrategy(linearSolverBuilder);
00332 
00333     if (solutionMethod == "QCAD Multi-Problem" || solutionMethod == "QCAD Poisson-Schrodinger") {
00334        // These QCAD solvers do not contain a primary Albany::Application instance and so albanyApp is null.
00335        // For now, do not resize the response vectors. FIXME sort out this issue.
00336        const RCP<Thyra::ModelEvaluator<double> > thyraModel = Thyra::epetraModelEvaluator(model, lowsFactory);
00337        const RCP<Piro::ObserverBase<double> > observer = rcp(new PiroObserver(app));
00338        return rcp(new Piro::NOXSolver<double>(piroParams, thyraModel, observer));
00339     }
00340     else {
00341 //      const RCP<AAdapt::AdaptiveModelFactory> thyraModelFactory = albanyApp->getAdaptSolMgr()->modelFactory();
00342       thyraModelFactory = albanyApp->getAdaptSolMgr()->modelFactory();
00343       const RCP<Thyra::ModelEvaluator<double> > thyraModel = thyraModelFactory->create(model, lowsFactory);
00344       const RCP<Piro::ObserverBase<double> > observer = rcp(new PiroObserver(app));
00345       return rcp(new Piro::NOXSolver<double>(piroParams, thyraModel, observer));
00346     }
00347   }
00348 
00349   const Teuchos::RCP<EpetraExt::ModelEvaluator> epetraSolver =
00350     this->createAndGetAlbanyApp(albanyApp, appComm, solverComm, initial_guess);
00351 
00352   if (solutionMethod == "QCAD Multi-Problem" || solutionMethod == "QCAD Poisson-Schrodinger") {
00353     return Thyra::epetraModelEvaluator(epetraSolver, Teuchos::null);
00354   }
00355   else {
00356 //    const RCP<AAdapt::AdaptiveModelFactory> thyraModelFactory = albanyApp->getAdaptSolMgr()->modelFactory();
00357     thyraModelFactory = albanyApp->getAdaptSolMgr()->modelFactory();
00358     return thyraModelFactory->create(epetraSolver, Teuchos::null);
00359   }
00360 }
00361 
00362 Teuchos::RCP<EpetraExt::ModelEvaluator>
00363 Albany::SolverFactory::createAlbanyAppAndModel(
00364   Teuchos::RCP<Albany::Application>& albanyApp,
00365   const Teuchos::RCP<const Epetra_Comm>& appComm,
00366   const Teuchos::RCP<const Epetra_Vector>& initial_guess)
00367 {
00368   // Create application
00369   albanyApp = rcp(new Albany::Application(appComm, appParams, initial_guess));
00370 
00371   // Validate Response list: may move inside individual Problem class
00372   const RCP<ParameterList> problemParams = Teuchos::sublist(appParams, "Problem");
00373   problemParams->sublist("Response Functions").
00374     validateParameters(*getValidResponseParameters(),0);
00375 
00376   // Create model evaluator
00377   Albany::ModelFactory modelFactory(appParams, albanyApp);
00378 
00379   return modelFactory.create();
00380 
00381 }
00382 
00383 int Albany::SolverFactory::checkSolveTestResults(
00384   int response_index,
00385   int parameter_index,
00386   const Epetra_Vector* g,
00387   const Epetra_MultiVector* dgdp) const
00388 {
00389   ParameterList* testParams = getTestParameters(response_index);
00390 
00391   int failures = 0;
00392   int comparisons = 0;
00393   const double relTol = testParams->get<double>("Relative Tolerance");
00394   const double absTol = testParams->get<double>("Absolute Tolerance");
00395 
00396   // Get number of responses (g) to test
00397   const int numResponseTests = testParams->get<int>("Number of Comparisons");
00398   if (numResponseTests > 0) {
00399     if (g == NULL || numResponseTests > g->MyLength()) failures += 1000;
00400     else { // do comparisons
00401       Teuchos::Array<double> testValues =
00402         testParams->get<Teuchos::Array<double> >("Test Values");
00403 
00404       TEUCHOS_TEST_FOR_EXCEPT(numResponseTests != testValues.size());
00405       for (int i=0; i<testValues.size(); i++) {
00406         failures += scaledCompare((*g)[i], testValues[i], relTol, absTol);
00407         comparisons++;
00408       }
00409     }
00410   }
00411 
00412   // Repeat comparisons for sensitivities
00413   Teuchos::ParameterList *sensitivityParams;
00414   std::string sensitivity_sublist_name =
00415     Albany::strint("Sensitivity Comparisons", parameter_index);
00416   if (parameter_index == 0 && !testParams->isSublist(sensitivity_sublist_name))
00417     sensitivityParams = testParams;
00418   else
00419     sensitivityParams = &(testParams->sublist(sensitivity_sublist_name));
00420   const int numSensTests =
00421     sensitivityParams->get<int>("Number of Sensitivity Comparisons");
00422   if (numSensTests > 0) {
00423     if (dgdp == NULL || numSensTests > dgdp->MyLength()) failures += 10000;
00424     else {
00425       for (int i=0; i<numSensTests; i++) {
00426         Teuchos::Array<double> testSensValues =
00427           sensitivityParams->get<Teuchos::Array<double> >(Albany::strint("Sensitivity Test Values",i));
00428         TEUCHOS_TEST_FOR_EXCEPT(dgdp->NumVectors() != testSensValues.size());
00429         for (int j=0; j<dgdp->NumVectors(); j++) {
00430           failures += scaledCompare((*dgdp)[j][i], testSensValues[j], relTol, absTol);
00431           comparisons++;
00432         }
00433       }
00434     }
00435   }
00436 
00437   storeTestResults(testParams, failures, comparisons);
00438 
00439   return failures;
00440 }
00441 
00442 int Albany::SolverFactory::checkDakotaTestResults(
00443   int response_index,
00444   const Teuchos::SerialDenseVector<int,double>* drdv) const
00445 {
00446   ParameterList* testParams = getTestParameters(response_index);
00447 
00448   int failures = 0;
00449   int comparisons = 0;
00450   const double relTol = testParams->get<double>("Relative Tolerance");
00451   const double absTol = testParams->get<double>("Absolute Tolerance");
00452 
00453   const int numDakotaTests = testParams->get<int>("Number of Dakota Comparisons");
00454   if (numDakotaTests > 0 && drdv != NULL) {
00455 
00456     if (numDakotaTests > drdv->length()) {
00457       failures += 100000;
00458     } else {
00459       // Read accepted test results
00460       Teuchos::Array<double> testValues =
00461         testParams->get<Teuchos::Array<double> >("Dakota Test Values");
00462 
00463       TEUCHOS_TEST_FOR_EXCEPT(numDakotaTests != testValues.size());
00464       for (int i=0; i<numDakotaTests; i++) {
00465         failures += scaledCompare((*drdv)[i], testValues[i], relTol, absTol);
00466         comparisons++;
00467       }
00468     }
00469   }
00470 
00471   storeTestResults(testParams, failures, comparisons);
00472 
00473   return failures;
00474 }
00475 
00476 int Albany::SolverFactory::checkAnalysisTestResults(
00477   int response_index,
00478   const Teuchos::RCP<Thyra::VectorBase<double> >& tvec) const
00479 {
00480   ParameterList* testParams = getTestParameters(response_index);
00481 
00482   int failures = 0;
00483   int comparisons = 0;
00484   const double relTol = testParams->get<double>("Relative Tolerance");
00485   const double absTol = testParams->get<double>("Absolute Tolerance");
00486 
00487   int numPiroTests = testParams->get<int>("Number of Piro Analysis Comparisons");
00488   if (numPiroTests > 0 && tvec != Teuchos::null) {
00489 
00490      // Create indexable thyra vector
00491       ::Thyra::DetachedVectorView<double> p(tvec);
00492 
00493     if (numPiroTests > p.subDim()) failures += 300000;
00494     else { // do comparisons
00495       // Read accepted test results
00496       Teuchos::Array<double> testValues =
00497         testParams->get<Teuchos::Array<double> >("Piro Analysis Test Values");
00498 
00499       TEUCHOS_TEST_FOR_EXCEPT(numPiroTests != testValues.size());
00500       for (int i=0; i<numPiroTests; i++) {
00501         failures += scaledCompare(p[i], testValues[i], relTol, absTol);
00502         comparisons++;
00503       }
00504     }
00505   }
00506 
00507   storeTestResults(testParams, failures, comparisons);
00508 
00509   return failures;
00510 }
00511 
00512 int Albany::SolverFactory::checkSGTestResults(
00513   int response_index,
00514   const Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>& g_sg,
00515   const Epetra_Vector* g_mean,
00516   const Epetra_Vector* g_std_dev) const
00517 {
00518   ParameterList* testParams = getTestParameters(response_index);
00519 
00520   int failures = 0;
00521   int comparisons = 0;
00522   const double relTol = testParams->get<double>("Relative Tolerance");
00523   const double absTol = testParams->get<double>("Absolute Tolerance");
00524 
00525   int numSGTests = testParams->get<int>("Number of Stochastic Galerkin Comparisons", 0);
00526   if (numSGTests > 0 && g_sg != Teuchos::null) {
00527     if (numSGTests > (*g_sg)[0].MyLength()) failures += 10000;
00528     else {
00529       for (int i=0; i<numSGTests; i++) {
00530         Teuchos::Array<double> testSGValues =
00531           testParams->get<Teuchos::Array<double> >
00532             (Albany::strint("Stochastic Galerkin Expansion Test Values",i));
00533         TEUCHOS_TEST_FOR_EXCEPT(g_sg->size() != testSGValues.size());
00534   for (int j=0; j<g_sg->size(); j++) {
00535     failures +=
00536       scaledCompare((*g_sg)[j][i], testSGValues[j], relTol, absTol);
00537           comparisons++;
00538         }
00539       }
00540     }
00541   }
00542 
00543   // Repeat comparisons for SG mean statistics
00544   int numMeanResponseTests = testParams->get<int>("Number of Stochastic Galerkin Mean Comparisons", 0);
00545   if (numMeanResponseTests > 0) {
00546 
00547     if (g_mean == NULL || numMeanResponseTests > g_mean->MyLength()) failures += 30000;
00548     else { // do comparisons
00549       Teuchos::Array<double> testValues =
00550         testParams->get<Teuchos::Array<double> >("Stochastic Galerkin Mean Test Values");
00551 
00552       TEUCHOS_TEST_FOR_EXCEPT(numMeanResponseTests != testValues.size());
00553       for (int i=0; i<testValues.size(); i++) {
00554         failures += scaledCompare((*g_mean)[i], testValues[i], relTol, absTol);
00555         comparisons++;
00556       }
00557     }
00558   }
00559 
00560   // Repeat comparisons for SG standard deviation statistics
00561   int numSDResponseTests = testParams->get<int>("Number of Stochastic Galerkin Standard Deviation Comparisons", 0);
00562   if (numSDResponseTests > 0) {
00563 
00564     if (g_std_dev == NULL || numSDResponseTests > g_std_dev->MyLength()) failures +=50000;
00565     else { // do comparisons
00566       Teuchos::Array<double> testValues =
00567         testParams->get<Teuchos::Array<double> >("Stochastic Galerkin Standard Deviation Test Values");
00568 
00569       TEUCHOS_TEST_FOR_EXCEPT(numSDResponseTests != testValues.size());
00570       for (int i=0; i<testValues.size(); i++) {
00571         failures += scaledCompare((*g_std_dev)[i], testValues[i], relTol, absTol);
00572         comparisons++;
00573       }
00574     }
00575   }
00576 
00577   storeTestResults(testParams, failures, comparisons);
00578 
00579   return failures;
00580 }
00581 
00582 ParameterList* Albany::SolverFactory::getTestParameters(int response_index) const
00583 {
00584   ParameterList* result;
00585 
00586   if (response_index == 0 && appParams->isSublist("Regression Results")) {
00587     result = &(appParams->sublist("Regression Results"));
00588   } else {
00589     result = &(appParams->sublist(Albany::strint("Regression Results", response_index)));
00590   }
00591 
00592   TEUCHOS_TEST_FOR_EXCEPTION(result->isType<std::string>("Test Values"), std::logic_error,
00593     "Array information in XML file must now be of type Array(double)\n");
00594   result->validateParametersAndSetDefaults(*getValidRegressionResultsParameters(),0);
00595 
00596   return result;
00597 }
00598 
00599 void Albany::SolverFactory::storeTestResults(
00600   ParameterList* testParams,
00601   int failures,
00602   int comparisons) const
00603 {
00604   // Store failures in param list (this requires mutable appParams!)
00605   testParams->set("Number of Failures", failures);
00606   testParams->set("Number of Comparisons Attempted", comparisons);
00607   *out << "\nCheckTestResults: Number of Comparisons Attempted = "
00608        << comparisons << std::endl;
00609 }
00610 
00611 int Albany::SolverFactory::scaledCompare(double x1, double x2, double relTol, double absTol) const
00612 {
00613   double diff = fabs(x1 - x2) / (0.5*fabs(x1) + 0.5*fabs(x2) + fabs(absTol));
00614 
00615   if (diff < relTol) return 0; //pass
00616   else               return 1; //fail
00617 }
00618 
00619 
00620 void Albany::SolverFactory::setSolverParamDefaults(
00621               ParameterList* appParams_, int myRank)
00622 {
00623     // Set the nonlinear solver method
00624     ParameterList& piroParams = appParams_->sublist("Piro");
00625     ParameterList& noxParams = piroParams.sublist("NOX");
00626     noxParams.set("Nonlinear Solver", "Line Search Based");
00627 
00628     // Set the printing parameters in the "Printing" sublist
00629     ParameterList& printParams = noxParams.sublist("Printing");
00630     printParams.set("MyPID", myRank);
00631     printParams.set("Output Precision", 3);
00632     printParams.set("Output Processor", 0);
00633     printParams.set("Output Information",
00634         NOX::Utils::OuterIteration +
00635         NOX::Utils::OuterIterationStatusTest +
00636         NOX::Utils::InnerIteration +
00637         NOX::Utils::Parameters +
00638         NOX::Utils::Details +
00639         NOX::Utils::LinearSolverDetails +
00640         NOX::Utils::Warning +
00641         NOX::Utils::Error);
00642 
00643     // Sublist for line search
00644     ParameterList& searchParams = noxParams.sublist("Line Search");
00645     searchParams.set("Method", "Full Step");
00646 
00647     // Sublist for direction
00648     ParameterList& dirParams = noxParams.sublist("Direction");
00649     dirParams.set("Method", "Newton");
00650     ParameterList& newtonParams = dirParams.sublist("Newton");
00651     newtonParams.set("Forcing Term Method", "Constant");
00652 
00653     // Sublist for linear solver for the Newton method
00654     ParameterList& lsParams = newtonParams.sublist("Linear Solver");
00655     lsParams.set("Aztec Solver", "GMRES");
00656     lsParams.set("Max Iterations", 43);
00657     lsParams.set("Tolerance", 1e-4);
00658     lsParams.set("Output Frequency", 20);
00659     lsParams.set("Preconditioner", "Ifpack");
00660 
00661     // Sublist for status tests
00662     ParameterList& statusParams = noxParams.sublist("Status Tests");
00663     statusParams.set("Test Type", "Combo");
00664     statusParams.set("Combo Type", "OR");
00665     statusParams.set("Number of Tests", 2);
00666     ParameterList& normF = statusParams.sublist("Test 0");
00667     normF.set("Test Type", "NormF");
00668     normF.set("Tolerance", 1.0e-8);
00669     normF.set("Norm Type", "Two Norm");
00670     normF.set("Scale Type", "Unscaled");
00671     ParameterList& maxiters = statusParams.sublist("Test 1");
00672     maxiters.set("Test Type", "MaxIters");
00673     maxiters.set("Maximum Iterations", 10);
00674 
00675 }
00676 
00677 RCP<const ParameterList>
00678 Albany::SolverFactory::getValidAppParameters() const
00679 {
00680   RCP<ParameterList> validPL = rcp(new ParameterList("ValidAppParams"));;
00681   validPL->sublist("Problem",            false, "Problem sublist");
00682   validPL->sublist("Debug Output",       false, "Debug Output sublist");
00683   validPL->sublist("Discretization",     false, "Discretization sublist");
00684   validPL->sublist("Quadrature",         false, "Quadrature sublist");
00685   validPL->sublist("Regression Results", false, "Regression Results sublist");
00686   validPL->sublist("VTK",                false, "DEPRECATED  VTK sublist");
00687   validPL->sublist("Piro",               false, "Piro sublist");
00688   validPL->sublist("Coupled System",     false, "Coupled system sublist");
00689 
00690   // validPL->set<std::string>("Jacobian Operator", "Have Jacobian", "Flag to allow Matrix-Free specification in Piro");
00691   // validPL->set<double>("Matrix-Free Perturbation", 3.0e-7, "delta in matrix-free formula");
00692 
00693   return validPL;
00694 }
00695 
00696 RCP<const ParameterList>
00697 Albany::SolverFactory::getValidRegressionResultsParameters() const
00698 {
00699   using Teuchos::Array;
00700   RCP<ParameterList> validPL = rcp(new ParameterList("ValidRegressionParams"));;
00701   Array<double> ta;; // std::string to be converted to teuchos array
00702 
00703   validPL->set<double>("Relative Tolerance", 1.0e-4,
00704           "Relative Tolerance used in regression testing");
00705   validPL->set<double>("Absolute Tolerance", 1.0e-8,
00706           "Absolute Tolerance used in regression testing");
00707 
00708   validPL->set<int>("Number of Comparisons", 0,
00709           "Number of responses to regress against");
00710   validPL->set<Array<double> >("Test Values", ta,
00711           "Array of regression values for responses");
00712 
00713   validPL->set<int>("Number of Sensitivity Comparisons", 0,
00714           "Number of sensitivity vectors to regress against");
00715 
00716   const int maxSensTests=10;
00717   for (int i=0; i<maxSensTests; i++) {
00718     validPL->set<Array<double> >(
00719        Albany::strint("Sensitivity Test Values",i), ta,
00720        Albany::strint("Array of regression values for Sensitivities w.r.t parameter",i));
00721   }
00722 
00723   validPL->set<int>("Number of Dakota Comparisons", 0,
00724           "Number of paramters from Dakota runs to regress against");
00725   validPL->set<Array<double> >("Dakota Test Values", ta,
00726           "Array of regression values for final parameters from Dakota runs");
00727 
00728   validPL->set<int>("Number of Piro Analysis Comparisons", 0,
00729           "Number of paramters from Analysis to regress against");
00730   validPL->set<Array<double> >("Piro Analysis Test Values", ta,
00731           "Array of regression values for final parameters from Analysis runs");
00732 
00733   validPL->set<int>("Number of Stochastic Galerkin Comparisons", 0,
00734           "Number of stochastic Galerkin expansions to regress against");
00735 
00736   const int maxSGTests=10;
00737   for (int i=0; i<maxSGTests; i++) {
00738     validPL->set<Array<double> >(
00739        Albany::strint("Stochastic Galerkin Expansion Test Values",i), ta,
00740        Albany::strint("Array of regression values for stochastic Galerkin expansions",i));
00741   }
00742 
00743   validPL->set<int>("Number of Stochastic Galerkin Mean Comparisons", 0,
00744           "Number of SG mean responses to regress against");
00745   validPL->set<Array<double> >("Stochastic Galerkin Mean Test Values", ta,
00746           "Array of regression values for SG mean responses");
00747   validPL->set<int>(
00748     "Number of Stochastic Galerkin Standard Deviation Comparisons", 0,
00749     "Number of SG standard deviation responses to regress against");
00750   validPL->set<Array<double> >(
00751     "Stochastic Galerkin Standard Deviation Test Values", ta,
00752     "Array of regression values for SG standard deviation responses");
00753 
00754   // These two are typically not set on input, just output.
00755   validPL->set<int>("Number of Failures", 0,
00756      "Output information from regression tests reporting number of failed tests");
00757   validPL->set<int>("Number of Comparisons Attempted", 0,
00758      "Output information from regression tests reporting number of comparisons attempted");
00759 
00760   return validPL;
00761 }
00762 
00763 RCP<const ParameterList>
00764 Albany::SolverFactory::getValidParameterParameters() const
00765 {
00766   RCP<ParameterList> validPL = rcp(new ParameterList("ValidParameterParams"));;
00767 
00768   validPL->set<int>("Number", 0);
00769   const int maxParameters = 100;
00770   for (int i=0; i<maxParameters; i++) {
00771     validPL->set<std::string>(Albany::strint("Parameter",i), "");
00772   }
00773   return validPL;
00774 }
00775 
00776 RCP<const ParameterList>
00777 Albany::SolverFactory::getValidResponseParameters() const
00778 {
00779   RCP<ParameterList> validPL = rcp(new ParameterList("ValidResponseParams"));;
00780 
00781   validPL->set<int>("Number of Response Vectors", 0);
00782   validPL->set<int>("Number", 0);
00783   validPL->set<int>("Equation", 0);
00784   const int maxParameters = 500;
00785   for (int i=0; i<maxParameters; i++) {
00786     validPL->set<std::string>(Albany::strint("Response",i), "");
00787     validPL->sublist(Albany::strint("ResponseParams",i));
00788     validPL->sublist(Albany::strint("Response Vector",i));
00789   }
00790   return validPL;
00791 }