Thyra Package Browser (Single Doxygen Collection) Version of the Day
ForwardSolveEpetraModelEval2DSimMain.cpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 #include "EpetraModelEval2DSim.hpp"
00045 #include "EpetraModelEval4DOpt.hpp"
00046 #include "Thyra_EpetraModelEvaluator.hpp"
00047 #include "Thyra_DefaultModelEvaluatorWithSolveFactory.hpp"
00048 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00049 #include "Thyra_DampenedNewtonNonlinearSolver.hpp"
00050 #include "Teuchos_VerboseObject.hpp"
00051 #include "Teuchos_CommandLineProcessor.hpp"
00052 #include "Teuchos_StandardCatchMacros.hpp"
00053 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00054 
00055 
00056 namespace {
00057   
00058 
00059 const Teuchos::RCP<const Epetra_Vector>
00060 createScalingVec(const double &scale, const Epetra_Map &map)
00061 {
00062   Teuchos::RCP<Epetra_Vector> scalingVec = Teuchos::rcp(new Epetra_Vector(map));
00063   scalingVec->PutScalar(scale);
00064   return scalingVec;
00065 }
00066 
00067 
00068 void scaleEpetraModelEvaluator( const double &s_x, const double &s_f,
00069   const Teuchos::Ptr<Thyra::EpetraModelEvaluator> &model
00070   )
00071 {
00072   if (s_x != 1.0) {
00073     model->setStateVariableScalingVec(
00074       createScalingVec(s_x, *model->getEpetraModel()->get_x_map())
00075       );
00076   }
00077   if (s_f != 1.0) {
00078     model->setStateFunctionScalingVec(
00079       createScalingVec(s_f, *model->getEpetraModel()->get_f_map())
00080       );
00081   }
00082 }
00083 
00084 
00085 } // namespace
00086 
00087 
00088 int main( int argc, char* argv[] )
00089 {
00090 
00091   using Teuchos::RCP;
00092   using Teuchos::rcp;
00093   using Teuchos::CommandLineProcessor;
00094   using Teuchos::outArg;
00095   typedef RCP<Thyra::VectorBase<double> > VectorPtr;
00096 
00097   bool success = true;
00098 
00099   try {
00100   
00101     //
00102     // Get options from the command line
00103     //
00104 
00105     CommandLineProcessor  clp;
00106     clp.throwExceptions(false);
00107     clp.addOutputSetupOptions(true);
00108 
00109     clp.setDocString(
00110       "This example program solves a simple 2 x 2 set of nonlinear equations using a simple\n"
00111       "dampened Newton method.\n\n"
00112 
00113       "The equations that are solved are:\n\n"
00114 
00115       "  f[0] =       x[0]      + x[1]*x[1] - p[0];\n"
00116       "  f[1] = d * ( x[0]*x[0] - x[1]      - p[1] );\n\n"
00117 
00118       "The Jacobian for these equations is nonsingular for every point except x=(-0.5,0.5)\n"
00119       "and x=(0.5,-0.5)  You can cause the Jacobian to be singular at the solution by setting\n"
00120       "p[0]=x[0]+x[1]*x[1] and p[1] = x[0]*x[0]-x[1] for these values of x.\n\n"
00121 
00122       "The equations are solved using a simple dampended Newton method that uses a Armijo\n"
00123       "line search which is implemented in the general class Thyra::DampenedNewtonNonlinearsolver\n"
00124       "You can get different levels of detail about the Newton method by adjustingthe command-line\n"
00125       "option \"verb-level\" (see above)\n"
00126       );
00127 
00128     double d = 10.0;
00129     clp.setOption( "d", &d, "Model constant d" );
00130     double p0 = 2.0;
00131     clp.setOption( "p0", &p0, "Model constant p[0]" );
00132     double p1 = 0.0;
00133     clp.setOption( "p1", &p1, "Model constant p[1]" );
00134     double x00 = 0.0;
00135     clp.setOption( "x00", &x00, "Initial guess for x[0]" );
00136     double x01 = 1.0;
00137     clp.setOption( "x01", &x01, "Initial guess for x[1]" );
00138     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
00139     setVerbosityLevelOption( "verb-level", &verbLevel, "Verbosity level.", &clp );
00140     double tol = 1e-10;
00141     clp.setOption( "tol", &tol, "Nonlinear solve tolerance" );
00142     int maxIters = 100;
00143     clp.setOption( "max-iters", &maxIters, "Maximum number of nonlinear iterations" );
00144     bool use4DOpt = false;
00145     clp.setOption( "use-4D-opt", "use-2D-sim", &use4DOpt,
00146       "Determines if the EpetraModelEval4DOpt or EpetraModelEval2DSim subclasses are used" );
00147     bool externalFactory = false;
00148     clp.setOption( "external-lowsf", "internal-lowsf", &externalFactory,
00149       "Determines of the Thyra::LinearOpWithSolveFactory is used externally or internally to the Thyra::EpetraModelEvaluator object" );
00150     bool showSetInvalidArg = false;
00151     clp.setOption( "show-set-invalid-arg", "no-show-set-invalid-arg", &showSetInvalidArg,
00152       "Determines if an attempt is made to set an invalid/unsupported ModelEvaluator input argument" );
00153     bool showGetInvalidArg = false;
00154     clp.setOption( "show-get-invalid-arg", "no-show-get-invalid-arg", &showGetInvalidArg,
00155       "Determines if an attempt is made to get an invalid/unsupported ModelEvaluator output argument (2DSim only)" );
00156     double s_x = 1.0;
00157     clp.setOption( "x-scale", &s_x, "State variables scaling." );
00158     double s_f = 1.0;
00159     clp.setOption( "f-scale", &s_f, "State function scaling." );
00160   
00161     CommandLineProcessor::EParseCommandLineReturn
00162       parse_return = clp.parse(argc,argv,&std::cerr);
00163 
00164     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
00165       return parse_return;
00166 
00167     RCP<Teuchos::FancyOStream>
00168       out = Teuchos::VerboseObjectBase::getDefaultOStream();
00169 
00170     *out << "\nCreating the nonlinear equations object ...\n";
00171     
00172     RCP<EpetraExt::ModelEvaluator> epetraModel;
00173     if(use4DOpt) {
00174       epetraModel = rcp(new EpetraModelEval4DOpt(0.0,0.0,p0,p1,d,x00,x01,p0,p1));
00175     }
00176     else {
00177       epetraModel = rcp(new EpetraModelEval2DSim(d,p0,p1,x00,x01,showGetInvalidArg));
00178     }
00179 
00180     *out << "\nCreating linear solver strategy ...\n";
00181 
00182     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00183     linearSolverBuilder.setParameterList(Teuchos::parameterList());
00184     RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
00185       lowsFactory = linearSolverBuilder.createLinearSolveStrategy("Amesos");
00186     // Above, we are just using the stratimkikos class
00187     // DefaultLinearSolverBuilder to create a default Amesos solver
00188     // (typically KLU) with a default set of options.  By setting a parameter
00189     // list on linearSolverBuilder, you build from a number of solvers.
00190 
00191     RCP<Thyra::EpetraModelEvaluator>
00192       epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator());
00193     
00194     RCP<Thyra::ModelEvaluator<double> > thyraModel;
00195     if(externalFactory) {
00196       epetraThyraModel->initialize(epetraModel, Teuchos::null);
00197       thyraModel = rcp(
00198         new Thyra::DefaultModelEvaluatorWithSolveFactory<double>(
00199           epetraThyraModel, lowsFactory
00200           )
00201         );
00202     }
00203     else {
00204       epetraThyraModel->initialize(epetraModel, lowsFactory);
00205       thyraModel = epetraThyraModel;
00206     }
00207 
00208     scaleEpetraModelEvaluator( s_x, s_f, epetraThyraModel.ptr() );
00209     
00210     if( showSetInvalidArg ) {
00211       *out << "\nAttempting to set an invalid input argument that throws an exception ...\n\n";
00212       Thyra::ModelEvaluatorBase::InArgs<double> inArgs = thyraModel->createInArgs();
00213       inArgs.set_x_dot(createMember(thyraModel->get_x_space()));
00214     }
00215     
00216     *out << "\nCreating the nonlinear solver and solving the equations ...\n\n";
00217 
00218     Thyra::DampenedNewtonNonlinearSolver<double> newtonSolver; // Set defaults
00219     newtonSolver.setVerbLevel(verbLevel);
00220     
00221     VectorPtr x = createMember(thyraModel->get_x_space());
00222     V_V( &*x, *thyraModel->getNominalValues().get_x() );
00223 
00224     Thyra::SolveCriteria<double> solveCriteria; // Sets defaults
00225     solveCriteria.solveMeasureType.set(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_NORM_RHS);
00226     solveCriteria.requestedTol = tol;
00227     solveCriteria.extraParameters = Teuchos::parameterList("Nonlinear Solve");
00228     solveCriteria.extraParameters->set("Max Iters",int(maxIters));
00229 
00230     newtonSolver.setModel(thyraModel);
00231     Thyra::SolveStatus<double>
00232       solveStatus = Thyra::solve( newtonSolver, &*x, &solveCriteria );
00233     
00234     *out << "\nNonlinear solver return status:\n";
00235     {
00236       Teuchos::OSTab tab(out);
00237       *out << solveStatus;
00238     }
00239     *out << "\nFinal solution for x=\n" << *x;
00240     
00241   }
00242   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,std::cerr,success)
00243     
00244   return  success ? 0 : 1;
00245 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines