MOOCHO  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MultiPeriodNLPThyraEpetraAdvDiffReactOptMain.cpp
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
6 // Copyright (2003) 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 (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
45 #include "MoochoPack_MoochoThyraSolver.hpp"
46 #include "Thyra_EpetraModelEvaluator.hpp"
47 #include "Thyra_EpetraLinearOp.hpp"
48 #include "Thyra_DefaultClusteredSpmdProductVectorSpace.hpp"
49 #include "Thyra_DefaultMultiPeriodModelEvaluator.hpp"
50 #include "Thyra_VectorSpaceTester.hpp"
51 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
52 #include "Thyra_DefaultInverseModelEvaluator.hpp"
53 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
54 #include "Thyra_DefaultMultipliedLinearOp.hpp"
55 #include "Thyra_TestingTools.hpp"
56 #include "Teuchos_OpaqueWrapper.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
58 #include "Teuchos_CommandLineProcessor.hpp"
59 #include "Teuchos_StandardCatchMacros.hpp"
60 #include "Teuchos_VerboseObject.hpp"
61 #include "Teuchos_Tuple.hpp"
62 #include "Teuchos_Utils.hpp"
63 #include "Teuchos_DefaultComm.hpp"
64 #ifdef HAVE_MPI
65 # include "Teuchos_DefaultMpiComm.hpp"
66 # include "Epetra_MpiComm.h"
67 #else
68 # include "Teuchos_DefaultSerialComm.hpp"
69 # include "Epetra_SerialComm.h"
70 #endif
71 
72 namespace {
73 
74 typedef AbstractLinAlgPack::value_type Scalar;
75 
76 } // namespace
77 
78 int main( int argc, char* argv[] )
79 {
80 
81  using Teuchos::rcp;
82  using Teuchos::rcp_dynamic_cast;
83  using Teuchos::rcp_implicit_cast;
84  using Teuchos::null;
85  using Teuchos::RCP;
86  using Teuchos::outArg;
87  using Teuchos::Array;
88  using Teuchos::tuple;
91  using Teuchos::OSTab;
93  using Teuchos::toString;
94  using Thyra::VectorBase;
96  typedef Thyra::ModelEvaluatorBase MEB;
97  typedef Thyra::Ordinal Ordinal;
101 
102  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
103 
104  const int procRank = mpiSession.getRank();
105  const int numProcs = mpiSession.getNProc();
106 
107  Teuchos::Time timer("");
108 
109  bool result, success = true;
110 
111  RCP<Teuchos::FancyOStream>
113 
114  try {
115 
117  GLpApp::AdvDiffReactOptModelCreator epetraModelCreator;
118 
119  // Create the solver object
120  MoochoThyraSolver solver;
121 
122  //
123  // Get options from the command line
124  //
125 
126  CommandLineProcessor clp;
127  clp.throwExceptions(false);
128  clp.addOutputSetupOptions(true);
129 
130  epetraModelCreator.setupCLP(&clp);
131  lowsfCreator.setupCLP(&clp);
132  solver.setupCLP(&clp);
133 
134  int numProcsPerCluster = -1;
135  clp.setOption( "num-procs-per-cluster", &numProcsPerCluster,
136  "Number of processes in a cluster (<=0 means only one cluster)." );
137  int numPeriodsPerCluster = 1;
138  clp.setOption( "num-periods-per-cluster", &numPeriodsPerCluster,
139  "Number of periods in a cluster." );
140  bool dumpAll = false;
141  clp.setOption( "dump-all", "no-dump-all", &dumpAll,
142  "Set to true, then a bunch of debugging output will be created for the clustered vector tests." );
143  bool skipSolve = false;
144  clp.setOption( "skip-solve", "no-skip-solve", &skipSolve,
145  "Temporary flag for skip solve for testing." );
146  double perturbedParamScaling = 1.0;
147  clp.setOption( "p-perturb-scaling", &perturbedParamScaling,
148  "Scaling for perturbed paramters from the initial forward solve." );
149  bool doMultiPeriod = true;
150  clp.setOption( "multi-period", "no-multi-period", &doMultiPeriod,
151  "Do a mulit-period solve or not." );
152  bool useOuterInverse = true;
153  clp.setOption( "use-outer-inverse", "use-inner-inverse", &useOuterInverse,
154  "Determines if the outer inverse model will be used or the inner inverse." );
155  double periodParamScale = 1.0;
156  clp.setOption( "period-param-scale", &periodParamScale,
157  "Sets the scaling factor to scale z[i] from one period to the next." );
158  bool initialSolveContinuation = false;
159  clp.setOption( "init-solve-continuation", "init-solve-all-at-once", &initialSolveContinuation,
160  "Determines if the inital solve is done using continuation or all at once." );
161  bool useStatelessPeriodModel = false;
162  clp.setOption( "use-stateless-period-model", "use-statefull-period-model", &useStatelessPeriodModel,
163  "Determines if a stateless or a statefull period model should be used or not." );
164  double stateInvError = 1e-8;
165  clp.setOption( "state-inv-error", &stateInvError,
166  "The error in the l2 norm of the state inverse solution error." );
167  double paramInvError = 1e-8;
168  clp.setOption( "param-inv-error", &paramInvError,
169  "The error in the l2 norm of the parameter inverse solution error." );
170 
171  CommandLineProcessor::EParseCommandLineReturn
172  parse_return = clp.parse(argc,argv,&std::cerr);
173 
174  if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
175  return parse_return;
176 
177  lowsfCreator.readParameters(out.get());
178  solver.readParameters(out.get());
179 
180  *out
181  << "\n***"
182  << "\n*** NLPThyraEpetraAdvDiffReactOptMain, Global numProcs = "<<numProcs
183  << "\n***\n";
184 
185  int clusterRank = -1;
186  int numClusters = -1;
187 
188 #ifdef HAVE_MPI
189 
190  RCP<OpaqueWrapper<MPI_Comm> >
191  intraClusterMpiComm = Teuchos::opaqueWrapper<MPI_Comm>(MPI_COMM_WORLD),
192  interClusterMpiComm = Teuchos::null;
193 
194  {
195  if ( numProcsPerCluster <= 0 ) {
196  *out
197  << "\nnumProcsPerCluster = " << numProcsPerCluster
198  << " <= 0: Setting to " << numProcs << "...\n";
199  numProcsPerCluster = numProcs;
200  }
201  *out << "\nCreating communicator for local cluster of "<<numProcsPerCluster<<" processes ...\n";
202  numClusters = numProcs/numProcsPerCluster;
203  const int remainingProcs = numProcs%numProcsPerCluster;
205  remainingProcs!=0,std::logic_error
206  ,"Error, The number of processes per cluster numProcsPerCluster="<<numProcsPerCluster
207  << " does not divide into the global number of processes numProcs="<<numProcs
208  << " and instead has remainder="<<remainingProcs<<"!"
209  );
210  // Determine which cluster this process is part of and what the global
211  // process ranges are.
212  clusterRank = procRank / numProcsPerCluster; // Integer division!
213  *out << "\nclusterRank = " << clusterRank << "\n";
214  const int firstClusterProcRank = clusterRank * numProcsPerCluster;
215  const int lastClusterProcRank = firstClusterProcRank + numProcsPerCluster - 1;
216  *out << "\nclusterProcRange = ["<<firstClusterProcRank<<","<<lastClusterProcRank<<"]\n";
217  // Create the communicator for this cluster of processes
218  *out << "\nCreating intraClusterMpiComm ...";
219  MPI_Comm rawIntraClusterMpiComm = MPI_COMM_NULL;
220  MPI_Comm_split(
221  MPI_COMM_WORLD // comm
222  ,clusterRank // color (will all be put in the same output comm)
223  ,0 // key (not important here)
224  ,&rawIntraClusterMpiComm // newcomm
225  );
226  intraClusterMpiComm = Teuchos::opaqueWrapper(rawIntraClusterMpiComm,MPI_Comm_free);
227  {
228  *out << "\nintraClusterMpiComm:";
229  Teuchos::OSTab tab(out);
230  int rank, size;
231  MPI_Comm_size(*intraClusterMpiComm,&size);
232  MPI_Comm_rank(*intraClusterMpiComm,&rank);
233  *out << "\nsize="<<size;
234  *out << "\nrank="<<rank;
235  *out << "\n";
236  }
237  // Create the communicator for just the root process in each cluster
238  *out << "\nCreating interClusterMpiComm ...";
239  MPI_Comm rawInterClusterMpiComm = MPI_COMM_NULL;
240  MPI_Comm_split(
241  MPI_COMM_WORLD // comm
242  ,procRank==firstClusterProcRank?0:MPI_UNDEFINED // color
243  ,0 // key
244  ,&rawInterClusterMpiComm // newcomm
245  );
246  if(rawInterClusterMpiComm!=MPI_COMM_NULL)
247  interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm,MPI_Comm_free);
248  else
249  interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm);
250  {
251  *out << "\ninterClusterMpiComm:";
252  Teuchos::OSTab tab(out);
253  if(*interClusterMpiComm==MPI_COMM_NULL) {
254  *out << " NULL\n";
255  }
256  else {
257  int rank, size;
258  MPI_Comm_size(*interClusterMpiComm,&size);
259  MPI_Comm_rank(*interClusterMpiComm,&rank);
260  *out << "\nsize="<<size;
261  *out << "\nrank="<<rank;
262  *out << "\n";
263  }
264  }
265  }
266 
267 #endif
268 
269  RCP<Epetra_Comm> comm = Teuchos::null;
270 #ifdef HAVE_MPI
271  comm = Teuchos::rcp(new Epetra_MpiComm(*intraClusterMpiComm));
272  Teuchos::set_extra_data(intraClusterMpiComm, "mpiComm", outArg(comm));
273 #else
274  comm = Teuchos::rcp(new Epetra_SerialComm());
275 #endif
276 
277  //
278  // Create the Thyra::ModelEvaluator object
279  //
280 
281  *out << "\nCreate the GLpApp::AdvDiffReactOptModel wrapper object ...\n";
282 
283  RCP<GLpApp::AdvDiffReactOptModel>
284  epetraModel = epetraModelCreator.createModel(comm);
285 
286  *out << "\nCreate the Thyra::LinearOpWithSolveFactory object ...\n";
287 
288  RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
289  lowsFactory = lowsfCreator.createLinearSolveStrategy("");
290  // ToDo: Set the output stream before calling above!
291 
292  *out << "\nCreate the Thyra::EpetraModelEvaluator wrapper object ...\n";
293 
294  RCP<Thyra::EpetraModelEvaluator>
295  epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator()); // Sets default options!
296  epetraThyraModel->setOStream(out);
297  epetraThyraModel->initialize(epetraModel,lowsFactory);
298 
299  *out
300  << "\nnx = " << epetraThyraModel->get_x_space()->dim()
301  << "\nnp = " << epetraThyraModel->get_p_space(0)->dim() << "\n";
302 
303  //
304  // Create the parallel product spaces for x and f
305  //
306 
307  RCP<const Thyra::ProductVectorSpaceBase<Scalar> >
308  x_bar_space, f_bar_space;
309 
310 #ifdef HAVE_MPI
311 
312  // For now just build and test these vector spaces if we are not doing a
313  // solve! We have a lot more work to do to the "Clustered" support
314  // software before this will work for a solve.
315 
316  if (skipSolve) {
317 
318  *out << "\nCreate block parallel vector spaces for multi-period model.x and model.f ...\n";
319  RCP<const Teuchos::Comm<Ordinal> >
320  intraClusterComm = rcp(new Teuchos::MpiComm<Ordinal>(intraClusterMpiComm)),
321  interClusterComm = Teuchos::createMpiComm<Ordinal>(interClusterMpiComm);
322  x_bar_space = Teuchos::rcp(
324  intraClusterComm
325  ,0 // clusterRootRank
326  ,interClusterComm
327  ,1 // numBlocks
328  ,tuple<RCP<const Thyra::VectorSpaceBase<Scalar> > >(
329  epetraThyraModel->get_x_space()
330  ).getRawPtr()
331  )
332  );
333  f_bar_space = Teuchos::rcp(
335  intraClusterComm
336  ,0 // clusterRootRank
337  ,interClusterComm
338  ,1 // numBlocks
339  ,tuple<RCP<const Thyra::VectorSpaceBase<Scalar> > >(
340  epetraThyraModel->get_f_space()
341  ).getRawPtr()
342  )
343  );
344 
345  Thyra::VectorSpaceTester<Scalar> vectorSpaceTester;
346  vectorSpaceTester.show_all_tests(true);
347  vectorSpaceTester.dump_all(dumpAll);
348 
349 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
350  RTOpPack::show_mpi_apply_op_dump = dumpAll;
351 #endif
352 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
354 #endif
355 
356  *out << "\nTesting the vector space x_bar_space ...\n";
357  result = vectorSpaceTester.check(*x_bar_space,OSTab(out).get());
358  if(!result) success = false;
359 
360  *out << "\nTesting the vector space f_bar_space ...\n";
361  result = vectorSpaceTester.check(*f_bar_space,OSTab(out).get());
362  if(!result) success = false;
363 
364  RCP<const VectorBase<Scalar> >
365  x0 = epetraThyraModel->getNominalValues().get_x();
366 
367 /* 2008/02/21: rabartl: I have commented this out since it is causing an MPI error?
368 
369  double nrm_x0;
370 
371  *out << "\nTiming a global reduction across just this cluster: ||x0||_1 = ";
372  timer.start(true);
373  nrm_x0 = Thyra::norm_1(*x0);
374  *out << nrm_x0 << "\n";
375  timer.stop();
376  *out << "\n time = " << timer.totalElapsedTime() << " seconds\n";
377 
378  *out << "\nTiming a global reduction across the entire set of processes: ||x0||_1 = ";
379  timer.start(true);
380  RTOpPack::ROpNorm1<Scalar> norm_1_op;
381  RCP<RTOpPack::ReductTarget> norm_1_targ = norm_1_op.reduct_obj_create();
382  const Teuchos::RCP<const Teuchos::Comm<Teuchos_Ordinal> > comm
383  = Teuchos::DefaultComm<Ordinal>::getComm();
384  const Teuchos::ArrayView<const Teuchos::Ptr<const VectorBase<Scalar> > >
385  vecs = Teuchos::tuple(Teuchos::ptrInArg(*x0));
386  Teuchos::dyn_cast<const Thyra::SpmdVectorBase<Scalar> >(*x0).applyOpImplWithComm(
387  comm.ptr(),
388  norm_1_op, vecs, Teuchos::null, norm_1_targ.ptr(),
389  0, -1, 0
390  );
391  nrm_x0 = norm_1_op(*norm_1_targ);
392  *out << nrm_x0 << "\n";
393  timer.stop();
394  *out << "\n time = " << timer.totalElapsedTime() << " seconds\n";
395 
396 */
397 
398 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
399  RTOpPack::show_mpi_apply_op_dump = false;
400 #endif
401 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
403 #endif
404 
405  }
406 
407 #endif // HAVE_MPI
408 
409  if(skipSolve) {
410 
411  if(success)
412  *out << "\nEnd Result: TEST PASSED" << std::endl;
413  else
414  *out << "\nEnd Result: TEST FAILED" << std::endl;
415 
416  return ( success ? 0 : 1 );
417 
418  }
419 
420  const int N = numPeriodsPerCluster;
421 
422  Array<RCP<Thyra::ModelEvaluator<Scalar> > >
423  inverseThyraModels(N);
424  if (useOuterInverse) {
425  *out << "\nUsing Thyra::DefaultInverseModelEvaluator for the objective function ...\n";
426  if ( useStatelessPeriodModel ) {
427  *out << "\nBuilding a single Thyra::DefaultInverseModelEvaluator object where the matching vector will be maintained externally ...\n";
428  }
429  else {
430  *out << "\nBuilding multiple Thyra::DefaultInverseModelEvaluator objects where the matching vector is held internally ...\n";
431  }
432  RCP<ParameterList> invMEPL = Teuchos::parameterList();
433  invMEPL->set( "Observation Multiplier", 1.0 );
434  invMEPL->set( "Parameter Multiplier", epetraModel->getDataPool()->getbeta() );
435  if ( useStatelessPeriodModel )
436  invMEPL->set( "Observation Target as Parameter", true );
437  RCP<const Thyra::EpetraLinearOp>
438  H = Thyra::epetraLinearOp(epetraModel->getDataPool()->getH(),"H"),
439  R = Thyra::epetraLinearOp(epetraModel->getDataPool()->getR(),"R");
440  RCP<const Thyra::MultiVectorBase<Scalar> >
442  epetraModel->get_B_bar(),
443  R->spmdRange()
444  );
445  RCP<const Thyra::LinearOpBase<Scalar> >
446  R_bar = Thyra::multiply<Scalar>(Thyra::adjoint<Scalar>(B_bar),R,B_bar);
447  for ( int i = 0; i < N; ++i ) {
448  if ( ( useStatelessPeriodModel && i==0 ) || !useStatelessPeriodModel ) {
449  RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
450  _inverseThyraModel = Thyra::defaultInverseModelEvaluator<Scalar>(
451  epetraThyraModel );
452  _inverseThyraModel->setParameterList(invMEPL);
453  _inverseThyraModel->set_observationMatchWeightingOp(H);
454  _inverseThyraModel->set_parameterRegularizationWeightingOp(R_bar);
455  inverseThyraModels[i] = _inverseThyraModel;
456  }
457  else {
458 #ifdef TEUCHOS_DEBUG
459  TEUCHOS_TEST_FOR_EXCEPT( ! ( useStatelessPeriodModel && i > 0 ) );
460 #endif
461  inverseThyraModels[i] = inverseThyraModels[0];
462  }
463  }
464  }
465  else {
466  *out << "\nUsing built-in inverse objective function ...\n";
468  N != 1, std::logic_error,
469  "Error, you can't have N = "<<N<<" > 1\n"
470  "and be using an internal inverse objective!" );
471  inverseThyraModels[0] = epetraThyraModel;
472  }
473 
474  const int p_index = 0;
475  const int z_index = 1;
476  const int z_p_index = 1; // Ordinal of the reaction rate parameter parameter subvector
477  const int z_x_index = 2; // Ordinal of the state matching subvector parameter
478  Array<int> z_indexes;
479  if (useStatelessPeriodModel)
480  z_indexes = tuple<int>(z_p_index, z_x_index);
481  else
482  z_indexes = tuple<int>(z_p_index);
483  Array<Array<RCP<const VectorBase<Scalar> > > > z;
484  const int g_index = ( useOuterInverse ? 1 : 0 );
485  Array<Scalar> weights;
486  RCP<VectorBase<Scalar> >
487  z_base = inverseThyraModels[0]->getNominalValues().get_p(z_index)->clone_v();
488  *out << "\nz_base =\n" << Teuchos::describe(*z_base,Teuchos::VERB_EXTREME);
489  Scalar scale_z_i = 1.0;
490  for ( int i = 0; i < N; ++i ) {
491  weights.push_back(1.0);
492  const RCP<VectorBase<Scalar> > z_i = z_base->clone_v();
493  Vt_S( z_i.ptr(), scale_z_i );
494  *out << "\nz["<<i<<"] =\n" << Teuchos::describe(*z_i,Teuchos::VERB_EXTREME);
495  if ( useStatelessPeriodModel ) {
496  z.push_back(
497  tuple<RCP<const VectorBase<Scalar> > >(
498  z_i,
499  null // We will set this again later!
500  )
501  );
502  }
503  else {
504  z.push_back(
505  tuple<RCP<const VectorBase<Scalar> > >(z_i)
506  );
507  }
508  scale_z_i *= periodParamScale;
509  }
510 
511  RCP<Thyra::ModelEvaluator<Scalar> >
512  thyraModel = inverseThyraModels[0];
513 
514  if (doMultiPeriod) {
515  thyraModel =
516  rcp(
518  N, inverseThyraModels, z_indexes, z, g_index, weights,
519  x_bar_space, f_bar_space
520  )
521  );
522  }
523 
524  MoochoSolver::ESolutionStatus solution_status;
525 
526  //
527  *out << "\n***\n*** Solving the initial forward problem\n***\n";
528  //
529 
530  // Set the solve mode to solve the forward problem
531  solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_FORWARD);
532 
533  // Save the solution for model.x and model.p to be used later
534  RCP<const VectorBase<Scalar> >
535  x_opt, // Will be set below
536  x_init, // Will be set below
537  p_opt = inverseThyraModels[0]->getNominalValues().get_p(0)->clone_v();
538 
539  *out << "\np_opt =\n" << Teuchos::describe(*p_opt,Teuchos::VERB_EXTREME);
540 
541  if ( initialSolveContinuation ) {
542 
543  *out << "\nSolving individual period problems one at time using continuation ...\n";
544 
545  RCP<ProductVectorBase<Scalar> >
546  x_opt_prod = rcp_dynamic_cast<ProductVectorBase<Scalar> >(
547  createMember( thyraModel->get_x_space() ), true );
548 
549  RCP<const VectorBase<Scalar> > period_x;
550 
551  for ( int i = 0; i < N; ++i ) {
552 
553  *out << "\nSolving period i = " << i << " using guess from last period ...\n";
554 
555  // Set the deliminator for the output files!
556  solver.getSolver().set_output_context("fwd-init-"+toString(i));
557 
558  // Set the period model
559  solver.setModel(inverseThyraModels[i]);
560 
561  // Set the initial guess and the parameter values
562  MEB::InArgs<Scalar> initialGuess = inverseThyraModels[i]->createInArgs();
563  initialGuess.set_p(z_index,z[i][0]->clone_v());
564  initialGuess.set_p(p_index,p_opt->clone_v());
565  if ( i == 0 ) {
566  // For the first period just use whatever initial guess is built
567  // into the model
568  }
569  else {
570  // Set the final solution for x from the last period!
571  initialGuess.set_x(period_x);
572  }
573  solver.setInitialGuess(initialGuess);
574 
575  // Solve the period model
576  solution_status = solver.solve();
577  TEUCHOS_TEST_FOR_EXCEPT( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED );
578 
579  // Save the final solution for the next period!
580  period_x = solver.getFinalPoint().get_x()->clone_v();
581  assign( x_opt_prod->getNonconstVectorBlock(i).ptr(), *period_x );
582  if ( useStatelessPeriodModel )
583  z[i][1] = period_x->clone_v(); // This is our matching vector!
584 
585  }
586 
587  x_opt = x_opt_prod;
588  x_init = x_opt->clone_v();
589 
590  if ( useStatelessPeriodModel ) {
592  thyraModel
593  )->reset_z(z);
594  }
595 
596  }
597  else {
598 
599  *out << "\nSolving all periods simultaniously ...\n";
600 
601  // Set the deliminator for the output files!
602  solver.getSolver().set_output_context("fwd-init");
603 
604  // Set the model
605  solver.setModel(thyraModel);
606 
607  // Set the initial guess from files (if specified on commandline)
608  solver.readInitialGuess(out.get());
609 
610  // Solve the initial forward problem
611  solution_status = solver.solve();
612  TEUCHOS_TEST_FOR_EXCEPT( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED );
613 
614  // Save the solution for model.x and model.p to be used later
615  x_opt = solver.getFinalPoint().get_x()->clone_v();
616  x_init = solver.getFinalPoint().get_x()->clone_v();
617 
618  }
619 
620  //
621  *out << "\n***\n*** Solving the perturbed forward problem\n***\n";
622  //
623 
624  // Set the deliminator for the output files!
625  solver.getSolver().set_output_context("fwd");
626 
627  // Set the solve mode to solve the forward problem
628  solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_FORWARD);
629 
630  // Set the model
631  solver.setModel(thyraModel);
632 
633  // Set the initial guess and the perturbed parameters
634  RCP<VectorBase<Scalar> >
635  p_init = p_opt->clone_v();
636  {
637  MEB::InArgs<Scalar> initialGuess = thyraModel->createInArgs();
638  initialGuess.setArgs(thyraModel->getNominalValues());
639  initialGuess.set_x(x_init);
640  Thyra::Vt_S(p_init.ptr(), perturbedParamScaling);
641  initialGuess.set_p(0,p_init);
642  //*out << "\nInitial Guess:\n" << Teuchos::describe(initialGuess,Teuchos::VERB_EXTREME);
643  solver.setInitialGuess(initialGuess);
644  }
645 
646  // Solve the perturbed forward problem
647  solution_status = solver.solve();
648 
649  // Save the solution for model.x and model.p to be used later
650  x_init = solver.getFinalPoint().get_x()->clone_v();
651  p_init = solver.getFinalPoint().get_p(0)->clone_v();
652 
653  *out
654  << "\nrelVectorErr(x_perturb,x_opt) = " << Thyra::relVectorErr(*x_init,*x_opt)
655  << "\nrelVectorErr(p_perturb,p_opt) = " << Thyra::relVectorErr(*p_init,*p_opt)
656  << "\n";
657 
658  //
659  *out << "\n***\n*** Solving the perturbed inverse problem\n***\n";
660  //
661 
662  // Set the deliminator for the output files!
663  solver.getSolver().set_output_context("inv");
664 
665  //TEUCHOS_TEST_FOR_EXCEPT("ToDo: We need to use the DefaultInverseModelEvaluator to set the matching vector correctly!");
666 
667  // Set the matching vector
668  if ( N > 1 ) {
670  !useOuterInverse, std::logic_error,
671  "Error, if N > 1, you have to use the outer inverse objective function\n"
672  "since each target vector will be different!" );
673  RCP<const ProductVectorBase<Scalar> >
674  x_opt_prod = rcp_dynamic_cast<const ProductVectorBase<Scalar> >(
675  rcp_implicit_cast<const VectorBase<Scalar> >(x_opt), true
676  ); // This cast can *not* fail!
677  for ( int i = 0; i < N; ++i ) {
679  inverseThyraModels[i], true
680  )->set_observationTarget(x_opt_prod->getVectorBlock(i));
681  }
682  }
683  else if ( 1 == N ) {
684  RCP<const ProductVectorBase<Scalar> >
685  x_opt_prod = rcp_dynamic_cast<const ProductVectorBase<Scalar> >(
686  rcp_implicit_cast<const VectorBase<Scalar> >(x_opt)
687  ); // This cast can fail!
688  RCP<const VectorBase<Scalar> >
689  x_opt_i =
690  ( !is_null(x_opt_prod)
691  ? x_opt_prod->getVectorBlock(0)
692  : rcp_implicit_cast<const VectorBase<Scalar> >(x_opt)
693  );
694  if (useOuterInverse) {
696  inverseThyraModels[0], true
697  )->set_observationTarget(x_opt_i);
698  }
699  else {
700  epetraModel->set_q(
701  Thyra::get_Epetra_Vector(*epetraModel->get_x_map(), x_opt_i)
702  );
703  }
704  }
705  else {
706  TEUCHOS_TEST_FOR_EXCEPT("Error, should not get here!");
707  }
708 
709  // Set the solve mode to solve the inverse problem
710  solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_OPTIMIZE);
711 
712  // Set the model
713  solver.setModel(thyraModel);
714 
715  // Set the initial guess for model.x and model.p
716  {
717  MEB::InArgs<Scalar> initialGuess = thyraModel->createInArgs();
718  initialGuess.setArgs(thyraModel->getNominalValues());
719  initialGuess.set_x(x_init);
720  initialGuess.set_p(0,p_init);
721  //*out << "\nInitial Guess:\n" << Teuchos::describe(initialGuess,Teuchos::VERB_EXTREME);
722  solver.setInitialGuess(initialGuess);
723  }
724 
725  // Solve the inverse problem
726  solution_status = solver.solve();
727  TEUCHOS_TEST_FOR_EXCEPT( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED );
728 
729  //
730  *out << "\n***\n*** Testing the error in the inversion\n***\n";
731  //
732 
733  // Get the inverted for solution and compare it to the optimal solution
734 
735  RCP<const VectorBase<Scalar> >
736  x_inv = solver.getFinalPoint().get_x(),
737  p_inv = solver.getFinalPoint().get_p(0);
738 
739  *out << "\np_opt =\n" << Teuchos::describe(*p_opt,Teuchos::VERB_EXTREME);
740  *out << "\np_inv =\n" << Teuchos::describe(*p_inv,Teuchos::VERB_EXTREME);
741 
742  const Scalar
743  x_err = Thyra::relVectorErr( *x_inv, *x_opt ),
744  p_err = Thyra::relVectorErr( *p_inv, *p_opt );
745 
746  const bool
747  x_test_passed = ( x_err <= stateInvError ),
748  p_test_passed = ( p_err <= paramInvError );
749 
750  *out
751  << "\nrelVectorErr(x_inv,x_opt) = " << x_err << " <= " << stateInvError
752  << " : " << Thyra::passfail(x_test_passed)
753  << "\nrelVectorErr(p_inv,p_opt) = " << p_err << " <= " << paramInvError
754  << " : " << Thyra::passfail(p_test_passed)
755  << "\n";
756 
757  // Write the final solution
758  solver.writeFinalSolution(out.get());
759 
760  // Write the final parameters to file
761  lowsfCreator.writeParamsFile(*lowsFactory);
762  solver.writeParamsFile();
763 
765  solution_status != MoochoSolver::SOLVE_RETURN_SOLVED
766  );
767 
768  }
769  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
770 
771  if(success)
772  *out << "\nEnd Result: TEST PASSED" << std::endl;
773  else
774  *out << "\nEnd Result: TEST FAILED" << std::endl;
775 
776  return ( success ? 0 : 1 );
777 
778 }
bool is_null(const boost::shared_ptr< T > &p)
basic_OSTab< char > OSTab
void setupCLP(Teuchos::CommandLineProcessor *clp)
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
const std::string passfail(const bool result)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
void readParameters(std::ostream *out)
void assign(MatrixOp *M_lhs, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Ordinal Ordinal
RCP< Thyra::LinearOpWithSolveFactoryBase< double > > createLinearSolveStrategy(const std::string &linearSolveStrategyName) const
static RCP< FancyOStream > getDefaultOStream()
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
void writeParamsFile(const Thyra::LinearOpWithSolveFactoryBase< double > &lowsFactory, const std::string &outputXmlFileName="") const
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
void setupCLP(Teuchos::CommandLineProcessor *clp)
void dump_all(const bool dump_all)
Teuchos::ScalarTraits< Scalar >::magnitudeType relVectorErr(const VectorBase< Scalar > &v1, const VectorBase< Scalar > &v2)
bool check(const VectorSpaceBase< Scalar > &vs, Teuchos::FancyOStream *out) const
void show_all_tests(const bool show_all_tests)
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
Teuchos::RCP< AdvDiffReactOptModel > createModel(const Teuchos::RCP< const Epetra_Comm > &comm, std::ostream *out=NULL) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
std::string toString(const T &t)