ForwardSolveEpetraModelEval2DSimMain.cpp

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

Generated on Tue Oct 20 12:47:19 2009 for EpetraExt/Thyra Adapters by doxygen 1.4.7