Thyra Package Browser (Single Doxygen Collection)  Version of the Day
ForwardSolveEpetraModelEval2DSimMain.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
6 // Copyright (2004) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 #include "EpetraModelEval2DSim.hpp"
45 #include "EpetraModelEval4DOpt.hpp"
47 #include "Thyra_DefaultModelEvaluatorWithSolveFactory.hpp"
48 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
49 #include "Thyra_DampenedNewtonNonlinearSolver.hpp"
50 #include "Teuchos_VerboseObject.hpp"
51 #include "Teuchos_CommandLineProcessor.hpp"
52 #include "Teuchos_StandardCatchMacros.hpp"
53 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
54 
55 
56 namespace {
57 
58 
59 const Teuchos::RCP<const Epetra_Vector>
60 createScalingVec(const double &scale, const Epetra_Map &map)
61 {
62  Teuchos::RCP<Epetra_Vector> scalingVec = Teuchos::rcp(new Epetra_Vector(map));
63  scalingVec->PutScalar(scale);
64  return scalingVec;
65 }
66 
67 
68 void scaleEpetraModelEvaluator( const double &s_x, const double &s_f,
69  const Teuchos::Ptr<Thyra::EpetraModelEvaluator> &model
70  )
71 {
72  if (s_x != 1.0) {
73  model->setStateVariableScalingVec(
74  createScalingVec(s_x, *model->getEpetraModel()->get_x_map())
75  );
76  }
77  if (s_f != 1.0) {
78  model->setStateFunctionScalingVec(
79  createScalingVec(s_f, *model->getEpetraModel()->get_f_map())
80  );
81  }
82 }
83 
84 
85 } // namespace
86 
87 
88 int main( int argc, char* argv[] )
89 {
90 
91  using Teuchos::RCP;
92  using Teuchos::rcp;
93  using Teuchos::CommandLineProcessor;
94  using Teuchos::outArg;
95  typedef RCP<Thyra::VectorBase<double> > VectorPtr;
96 
97  bool success = true;
98 
99  try {
100 
101  //
102  // Get options from the command line
103  //
104 
105  CommandLineProcessor clp;
106  clp.throwExceptions(false);
107  clp.addOutputSetupOptions(true);
108 
109  clp.setDocString(
110  "This example program solves a simple 2 x 2 set of nonlinear equations using a simple\n"
111  "dampened Newton method.\n\n"
112 
113  "The equations that are solved are:\n\n"
114 
115  " f[0] = x[0] + x[1]*x[1] - p[0];\n"
116  " f[1] = d * ( x[0]*x[0] - x[1] - p[1] );\n\n"
117 
118  "The Jacobian for these equations is nonsingular for every point except x=(-0.5,0.5)\n"
119  "and x=(0.5,-0.5) You can cause the Jacobian to be singular at the solution by setting\n"
120  "p[0]=x[0]+x[1]*x[1] and p[1] = x[0]*x[0]-x[1] for these values of x.\n\n"
121 
122  "The equations are solved using a simple dampended Newton method that uses a Armijo\n"
123  "line search which is implemented in the general class Thyra::DampenedNewtonNonlinearsolver\n"
124  "You can get different levels of detail about the Newton method by adjustingthe command-line\n"
125  "option \"verb-level\" (see above)\n"
126  );
127 
128  double d = 10.0;
129  clp.setOption( "d", &d, "Model constant d" );
130  double p0 = 2.0;
131  clp.setOption( "p0", &p0, "Model constant p[0]" );
132  double p1 = 0.0;
133  clp.setOption( "p1", &p1, "Model constant p[1]" );
134  double x00 = 0.0;
135  clp.setOption( "x00", &x00, "Initial guess for x[0]" );
136  double x01 = 1.0;
137  clp.setOption( "x01", &x01, "Initial guess for x[1]" );
138  Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
139  setVerbosityLevelOption( "verb-level", &verbLevel, "Verbosity level.", &clp );
140  double tol = 1e-10;
141  clp.setOption( "tol", &tol, "Nonlinear solve tolerance" );
142  int maxIters = 100;
143  clp.setOption( "max-iters", &maxIters, "Maximum number of nonlinear iterations" );
144  bool use4DOpt = false;
145  clp.setOption( "use-4D-opt", "use-2D-sim", &use4DOpt,
146  "Determines if the EpetraModelEval4DOpt or EpetraModelEval2DSim subclasses are used" );
147  bool externalFactory = false;
148  clp.setOption( "external-lowsf", "internal-lowsf", &externalFactory,
149  "Determines of the Thyra::LinearOpWithSolveFactory is used externally or internally to the Thyra::EpetraModelEvaluator object" );
150  bool showSetInvalidArg = false;
151  clp.setOption( "show-set-invalid-arg", "no-show-set-invalid-arg", &showSetInvalidArg,
152  "Determines if an attempt is made to set an invalid/unsupported ModelEvaluator input argument" );
153  bool showGetInvalidArg = false;
154  clp.setOption( "show-get-invalid-arg", "no-show-get-invalid-arg", &showGetInvalidArg,
155  "Determines if an attempt is made to get an invalid/unsupported ModelEvaluator output argument (2DSim only)" );
156  double s_x = 1.0;
157  clp.setOption( "x-scale", &s_x, "State variables scaling." );
158  double s_f = 1.0;
159  clp.setOption( "f-scale", &s_f, "State function scaling." );
160 
161  CommandLineProcessor::EParseCommandLineReturn
162  parse_return = clp.parse(argc,argv,&std::cerr);
163 
164  if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
165  return parse_return;
166 
167  RCP<Teuchos::FancyOStream>
168  out = Teuchos::VerboseObjectBase::getDefaultOStream();
169 
170  *out << "\nCreating the nonlinear equations object ...\n";
171 
172  RCP<EpetraExt::ModelEvaluator> epetraModel;
173  if(use4DOpt) {
174  epetraModel = rcp(new EpetraModelEval4DOpt(0.0,0.0,p0,p1,d,x00,x01,p0,p1));
175  }
176  else {
177  epetraModel = rcp(new EpetraModelEval2DSim(d,p0,p1,x00,x01,showGetInvalidArg));
178  }
179 
180  *out << "\nCreating linear solver strategy ...\n";
181 
182  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
183  linearSolverBuilder.setParameterList(Teuchos::parameterList());
184  RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
185  lowsFactory = linearSolverBuilder.createLinearSolveStrategy("Amesos");
186  // Above, we are just using the stratimkikos class
187  // DefaultLinearSolverBuilder to create a default Amesos solver
188  // (typically KLU) with a default set of options. By setting a parameter
189  // list on linearSolverBuilder, you build from a number of solvers.
190 
191  RCP<Thyra::EpetraModelEvaluator>
192  epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator());
193 
194  RCP<Thyra::ModelEvaluator<double> > thyraModel;
195  if(externalFactory) {
196  epetraThyraModel->initialize(epetraModel, Teuchos::null);
197  thyraModel = rcp(
198  new Thyra::DefaultModelEvaluatorWithSolveFactory<double>(
199  epetraThyraModel, lowsFactory
200  )
201  );
202  }
203  else {
204  epetraThyraModel->initialize(epetraModel, lowsFactory);
205  thyraModel = epetraThyraModel;
206  }
207 
208  scaleEpetraModelEvaluator( s_x, s_f, epetraThyraModel.ptr() );
209 
210  if( showSetInvalidArg ) {
211  *out << "\nAttempting to set an invalid input argument that throws an exception ...\n\n";
212  Thyra::ModelEvaluatorBase::InArgs<double> inArgs = thyraModel->createInArgs();
213  inArgs.set_x_dot(createMember(thyraModel->get_x_space()));
214  }
215 
216  *out << "\nCreating the nonlinear solver and solving the equations ...\n\n";
217 
218  Thyra::DampenedNewtonNonlinearSolver<double> newtonSolver; // Set defaults
219  newtonSolver.setVerbLevel(verbLevel);
220 
221  VectorPtr x = createMember(thyraModel->get_x_space());
222  V_V( &*x, *thyraModel->getNominalValues().get_x() );
223 
224  Thyra::SolveCriteria<double> solveCriteria; // Sets defaults
225  solveCriteria.solveMeasureType.set(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_NORM_RHS);
226  solveCriteria.requestedTol = tol;
227  solveCriteria.extraParameters = Teuchos::parameterList("Nonlinear Solve");
228  solveCriteria.extraParameters->set("Max Iters",int(maxIters));
229 
230  newtonSolver.setModel(thyraModel);
231  Thyra::SolveStatus<double>
232  solveStatus = Thyra::solve( newtonSolver, &*x, &solveCriteria );
233 
234  *out << "\nNonlinear solver return status:\n";
235  {
236  Teuchos::OSTab tab(out);
237  *out << solveStatus;
238  }
239  *out << "\nFinal solution for x=\n" << *x;
240 
241  }
242  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,std::cerr,success)
243 
244  return success ? 0 : 1;
245 }
Concrete Adapter subclass that takes an EpetraExt::ModelEvaluator object and wraps it as a Thyra::Mod...
int main(int argc, char *argv[])