ForwardSolveEpetraModelEval2DSimMain.cpp

00001 #include "EpetraModelEval2DSim.hpp"
00002 #include "EpetraModelEval4DOpt.hpp"
00003 #include "Thyra_EpetraModelEvaluator.hpp"
00004 #include "Thyra_DefaultModelEvaluatorWithSolveFactory.hpp"
00005 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00006 #include "Thyra_DampenedNewtonNonlinearSolver.hpp"
00007 #include "Teuchos_VerboseObject.hpp"
00008 #include "Teuchos_CommandLineProcessor.hpp"
00009 #include "Teuchos_StandardCatchMacros.hpp"
00010 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00011 
00012 
00013 namespace {
00014   
00015 
00016 const Teuchos::RCP<const Epetra_Vector>
00017 createScalingVec(const double &scale, const Epetra_Map &map)
00018 {
00019   Teuchos::RCP<Epetra_Vector> scalingVec = Teuchos::rcp(new Epetra_Vector(map));
00020   scalingVec->PutScalar(scale);
00021   return scalingVec;
00022 }
00023 
00024 
00025 void scaleEpetraModelEvaluator( const double &s_x, const double &s_f,
00026   const Teuchos::Ptr<Thyra::EpetraModelEvaluator> &model
00027   )
00028 {
00029   if (s_x != 1.0) {
00030     model->setStateVariableScalingVec(
00031       createScalingVec(s_x, *model->getEpetraModel()->get_x_map())
00032       );
00033   }
00034   if (s_f != 1.0) {
00035     model->setStateFunctionScalingVec(
00036       createScalingVec(s_f, *model->getEpetraModel()->get_f_map())
00037       );
00038   }
00039 }
00040 
00041 
00042 } // namespace
00043 
00044 
00045 int main( int argc, char* argv[] )
00046 {
00047 
00048   using Teuchos::RCP;
00049   using Teuchos::rcp;
00050   using Teuchos::CommandLineProcessor;
00051   using Teuchos::outArg;
00052   typedef RCP<Thyra::VectorBase<double> > VectorPtr;
00053 
00054   bool success = true;
00055 
00056   try {
00057   
00058     //
00059     // Get options from the command line
00060     //
00061 
00062     CommandLineProcessor  clp;
00063     clp.throwExceptions(false);
00064     clp.addOutputSetupOptions(true);
00065 
00066     clp.setDocString(
00067       "This example program solves a simple 2 x 2 set of nonlinear equations using a simple\n"
00068       "dampened Newton method.\n\n"
00069 
00070       "The equations that are solved are:\n\n"
00071 
00072       "  f[0] =       x[0]      + x[1]*x[1] - p[0];\n"
00073       "  f[1] = d * ( x[0]*x[0] - x[1]      - p[1] );\n\n"
00074 
00075       "The Jacobian for these equations is nonsingular for every point except x=(-0.5,0.5)\n"
00076       "and x=(0.5,-0.5)  You can cause the Jacobian to be singular at the solution by setting\n"
00077       "p[0]=x[0]+x[1]*x[1] and p[1] = x[0]*x[0]-x[1] for these values of x.\n\n"
00078 
00079       "The equations are solved using a simple dampended Newton method that uses a Armijo\n"
00080       "line search which is implemented in the general class Thyra::DampenedNewtonNonlinearsolver\n"
00081       "You can get different levels of detail about the Newton method by adjustingthe command-line\n"
00082       "option \"verb-level\" (see above)\n"
00083       );
00084 
00085     double d = 10.0;
00086     clp.setOption( "d", &d, "Model constant d" );
00087     double p0 = 2.0;
00088     clp.setOption( "p0", &p0, "Model constant p[0]" );
00089     double p1 = 0.0;
00090     clp.setOption( "p1", &p1, "Model constant p[1]" );
00091     double x00 = 0.0;
00092     clp.setOption( "x00", &x00, "Initial guess for x[0]" );
00093     double x01 = 1.0;
00094     clp.setOption( "x01", &x01, "Initial guess for x[1]" );
00095     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
00096     setVerbosityLevelOption( "verb-level", &verbLevel, "Verbosity level.", &clp );
00097     double tol = 1e-10;
00098     clp.setOption( "tol", &tol, "Nonlinear solve tolerance" );
00099     int maxIters = 100;
00100     clp.setOption( "max-iters", &maxIters, "Maximum number of nonlinear iterations" );
00101     bool use4DOpt = false;
00102     clp.setOption( "use-4D-opt", "use-2D-sim", &use4DOpt,
00103       "Determines if the EpetraModelEval4DOpt or EpetraModelEval2DSim subclasses are used" );
00104     bool externalFactory = false;
00105     clp.setOption( "external-lowsf", "internal-lowsf", &externalFactory,
00106       "Determines of the Thyra::LinearOpWithSolveFactory is used externally or internally to the Thyra::EpetraModelEvaluator object" );
00107     bool showSetInvalidArg = false;
00108     clp.setOption( "show-set-invalid-arg", "no-show-set-invalid-arg", &showSetInvalidArg,
00109       "Determines if an attempt is made to set an invalid/unsupported ModelEvaluator input argument" );
00110     bool showGetInvalidArg = false;
00111     clp.setOption( "show-get-invalid-arg", "no-show-get-invalid-arg", &showGetInvalidArg,
00112       "Determines if an attempt is made to get an invalid/unsupported ModelEvaluator output argument (2DSim only)" );
00113     double s_x = 1.0;
00114     clp.setOption( "x-scale", &s_x, "State variables scaling." );
00115     double s_f = 1.0;
00116     clp.setOption( "f-scale", &s_f, "State function scaling." );
00117   
00118     CommandLineProcessor::EParseCommandLineReturn
00119       parse_return = clp.parse(argc,argv,&std::cerr);
00120 
00121     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
00122       return parse_return;
00123 
00124     RCP<Teuchos::FancyOStream>
00125       out = Teuchos::VerboseObjectBase::getDefaultOStream();
00126 
00127     *out << "\nCreating the nonlinear equations object ...\n";
00128     
00129     RCP<EpetraExt::ModelEvaluator> epetraModel;
00130     if(use4DOpt) {
00131       epetraModel = rcp(new EpetraModelEval4DOpt(0.0,0.0,p0,p1,d,x00,x01,p0,p1));
00132     }
00133     else {
00134       epetraModel = rcp(new EpetraModelEval2DSim(d,p0,p1,x00,x01,showGetInvalidArg));
00135     }
00136 
00137     *out << "\nCreating linear solver strategy ...\n";
00138 
00139     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00140     linearSolverBuilder.setParameterList(Teuchos::parameterList());
00141     RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
00142       lowsFactory = linearSolverBuilder.createLinearSolveStrategy("Amesos");
00143     // Above, we are just using the stratimkikos class
00144     // DefaultLinearSolverBuilder to create a default Amesos solver
00145     // (typically KLU) with a default set of options.  By setting a parameter
00146     // list on linearSolverBuilder, you build from a number of solvers.
00147 
00148     RCP<Thyra::EpetraModelEvaluator>
00149       epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator());
00150     
00151     RCP<Thyra::ModelEvaluator<double> > thyraModel;
00152     if(externalFactory) {
00153       epetraThyraModel->initialize(epetraModel, Teuchos::null);
00154       thyraModel = rcp(
00155         new Thyra::DefaultModelEvaluatorWithSolveFactory<double>(
00156           epetraThyraModel, lowsFactory
00157           )
00158         );
00159     }
00160     else {
00161       epetraThyraModel->initialize(epetraModel, lowsFactory);
00162       thyraModel = epetraThyraModel;
00163     }
00164 
00165     scaleEpetraModelEvaluator( s_x, s_f, epetraThyraModel.ptr() );
00166     
00167     if( showSetInvalidArg ) {
00168       *out << "\nAttempting to set an invalid input argument that throws an exception ...\n\n";
00169       Thyra::ModelEvaluatorBase::InArgs<double> inArgs = thyraModel->createInArgs();
00170       inArgs.set_x_dot(createMember(thyraModel->get_x_space()));
00171     }
00172     
00173     *out << "\nCreating the nonlinear solver and solving the equations ...\n\n";
00174 
00175     Thyra::DampenedNewtonNonlinearSolver<double> newtonSolver; // Set defaults
00176     newtonSolver.setVerbLevel(verbLevel);
00177     
00178     VectorPtr x = createMember(thyraModel->get_x_space());
00179     V_V( &*x, *thyraModel->getNominalValues().get_x() );
00180 
00181     Thyra::SolveCriteria<double> solveCriteria; // Sets defaults
00182     solveCriteria.solveMeasureType.set(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_NORM_RHS);
00183     solveCriteria.requestedTol = tol;
00184     solveCriteria.extraParameters = Teuchos::parameterList("Nonlinear Solve");
00185     solveCriteria.extraParameters->set("Max Iters",int(maxIters));
00186 
00187     newtonSolver.setModel(thyraModel);
00188     Thyra::SolveStatus<double>
00189       solveStatus = Thyra::solve( newtonSolver, &*x, &solveCriteria );
00190     
00191     *out << "\nNonlinear solver return status:\n";
00192     {
00193       Teuchos::OSTab tab(out);
00194       *out << solveStatus;
00195     }
00196     *out << "\nFinal solution for x=\n" << *x;
00197     
00198   }
00199   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,std::cerr,success)
00200     
00201   return  success ? 0 : 1;
00202 }

Generated on Thu Feb 11 16:28:43 2010 for EpetraExt/Thyra Adapters by  doxygen 1.4.7