PlayaNOXSolver.cpp

00001 // $Id$ 
00002 // $Source$ 
00003 
00004 //@HEADER
00005 //   
00006 //@HEADER
00007 
00008 #include "PlayaNOXSolver.hpp"         
00009 #include "NOX_StatusTest_SafeCombo.hpp"         
00010 #include "NOX.H"         
00011 //#include "NOX_Parameter_Teuchos2NOX.H"         
00012 #include "PlayaLinearSolverBuilder.hpp"
00013 #include "Teuchos_Time.hpp"
00014 #include "Teuchos_TimeMonitor.hpp"
00015 #include "PlayaOut.hpp"
00016 #include "PlayaTabs.hpp"
00017 #include "PlayaExceptions.hpp"
00018 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00019 #include "PlayaVectorImpl.hpp"
00020 #include "PlayaLinearOperatorImpl.hpp"
00021 #endif
00022 
00023 
00024 using namespace NOX;
00025 using namespace NOX::NOXPlaya;
00026 using namespace Teuchos;
00027 using namespace Playa;
00028 
00029 
00030 static Time& noxSolverTimer() 
00031 {
00032   static RCP<Time> rtn 
00033     = TimeMonitor::getNewTimer("NOX solve"); 
00034   return *rtn;
00035 }
00036 
00037 NOXSolver::NOXSolver(const ParameterList& params)
00038   : linSolver_(),
00039     statusTest_(),
00040     params_(),
00041     printParams_()
00042 {
00043   TEUCHOS_TEST_FOR_EXCEPTION(!params.isSublist("NOX Solver"), runtime_error,
00044                      "did not find NOX Solver sublist in " << params);
00045   
00046   params_ = params.sublist("NOX Solver");
00047   /* NOX wants to have the process ID in a parameter list???? */
00048   params_.sublist("Printing").set("MyPID", MPIComm::world().getRank());
00049 
00050   if (params_.isSublist("Status Test"))
00051     {
00052       statusTest_ = StatusTestBuilder::makeStatusTest(params_);
00053     }
00054   else
00055     {
00056       RCP<StatusTest::Generic> A = rcp(new StatusTest::NormF(1.0e-12));
00057       RCP<StatusTest::Generic> B = rcp(new StatusTest::MaxIters(20));
00058       statusTest_ = 
00059         rcp(new StatusTest::SafeCombo(StatusTest::SafeCombo::OR, A, B));
00060     }
00061   
00062   if (params_.isSublist("Linear Solver"))
00063     {
00064       linSolver_ = LinearSolverBuilder::createSolver(params_);
00065     }
00066   else
00067   {
00068     TEUCHOS_TEST_FOR_EXCEPTION(!params_.isSublist("Linear Solver"),
00069       RuntimeError, "no linear solver specified in NOX parameters");
00070   }
00071   
00072   if (params_.isSublist("Printing"))
00073     {
00074       printParams_ = params_.sublist("Printing");
00075     }
00076   
00077   TEUCHOS_TEST_FOR_EXCEPTION(linSolver_.ptr().get()==0, runtime_error,
00078                      "null linear solver object in NOXSolver ctor");
00079 
00080   TEUCHOS_TEST_FOR_EXCEPTION(statusTest_.get()==0, runtime_error,
00081                      "null status test object in NOXSolver ctor");
00082 
00083 }
00084 
00085 NOXSolver::NOXSolver(const ParameterList& nonlinParams,
00086       const LinearSolver<double>& linSolver)
00087   : linSolver_(linSolver),
00088     statusTest_(),
00089     params_(),
00090     printParams_()
00091 {
00092   Tabs tab(0);
00093   TEUCHOS_TEST_FOR_EXCEPTION(!nonlinParams.isSublist("NOX Solver"), runtime_error,
00094                      "did not find NOX Solver sublist in " << nonlinParams);
00095   
00096   params_ = nonlinParams.sublist("NOX Solver");
00097   /* NOX wants to have the process ID in a parameter list???? */
00098   params_.sublist("Printing").set("MyPID", MPIComm::world().getRank());
00099 
00100   if (params_.isSublist("Status Test"))
00101     {
00102       statusTest_ = StatusTestBuilder::makeStatusTest(params_);
00103     }
00104   else
00105     {
00106       RCP<StatusTest::Generic> A = rcp(new StatusTest::NormF(1.0e-12));
00107       RCP<StatusTest::Generic> B = rcp(new StatusTest::MaxIters(20));
00108       statusTest_ = 
00109         rcp(new StatusTest::SafeCombo(StatusTest::SafeCombo::OR, A, B));
00110     }
00111   
00112   if (params_.isSublist("Linear Solver"))
00113     {
00114       Out::root() << tab << "WARNING: linear solver in NOX parameter list "
00115         "will be overridden by alternate solver" << std::endl;
00116     }
00117   
00118   if (params_.isSublist("Printing"))
00119     {
00120       printParams_ = params_.sublist("Printing");
00121     }
00122   
00123   TEUCHOS_TEST_FOR_EXCEPTION(linSolver_.ptr().get()==0, runtime_error,
00124                      "null linear solver object in NOXSolver ctor");
00125 
00126   TEUCHOS_TEST_FOR_EXCEPTION(statusTest_.get()==0, runtime_error,
00127                      "null status test object in NOXSolver ctor");
00128 
00129 }
00130 
00131 
00132 
00133 
00134 SolverState<double>
00135 NOXSolver::solve(const NonlinearOperator<double>& F, 
00136                  Playa::Vector<double>& solnVec) const 
00137 {
00138   TimeMonitor timer(noxSolverTimer());
00139 
00140   Vector<double> x0 = F.getInitialGuess();
00141   RCP<NOX::NOXPlaya::Group> grp = rcp(new NOX::NOXPlaya::Group(x0, F, linSolver_));
00142   RCP<Teuchos::ParameterList> noxParams 
00143     = Teuchos::rcp(&params_, false);
00144   RCP<NOX::Solver::Generic> solver 
00145     = NOX::Solver::buildSolver(grp, statusTest_, noxParams);
00146 
00147   NOX::StatusTest::StatusType rtn = solver->solve();
00148 
00149 
00150   const NOX::NOXPlaya::Group* solnGrp 
00151     = dynamic_cast<const NOX::NOXPlaya::Group*>(&(solver->getSolutionGroup()));
00152 
00153   TEUCHOS_TEST_FOR_EXCEPTION(solnGrp==0, runtime_error,
00154                      "Solution group could not be cast to NOX::NOXPlaya::Group");
00155 
00156   double resid = solnGrp->getNormF();
00157   int itersUsed = solver->getNumIterations();
00158 
00159 
00160 
00161   const NOX::NOXPlaya::Vector* x 
00162     = dynamic_cast<const NOX::NOXPlaya::Vector*>(&(solnGrp->getX()));
00163 
00164   TEUCHOS_TEST_FOR_EXCEPTION(x==0, runtime_error,
00165     "Solution vector could not be cast to NOX::NOXPlaya::Vector");
00166   
00167   solnVec = x->getPlayaVector();
00168 
00169   if (rtn==NOX::StatusTest::Converged)
00170   {
00171     return SolverState<double>(SolveConverged, "Solve converged", itersUsed, resid);
00172   }
00173   else if (rtn==NOX::StatusTest::Unconverged)
00174   {
00175     return SolverState<double>(SolveFailedToConverge, "Solve failed to converge", itersUsed, resid);
00176   }
00177   else
00178   {
00179     return SolverState<double>(SolveCrashed, "Solve crashed", itersUsed, resid);
00180   }
00181 
00182 }

doxygen