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"
65 # include "Teuchos_DefaultMpiComm.hpp"
66 # include "Epetra_MpiComm.h"
68 # include "Teuchos_DefaultSerialComm.hpp"
69 # include "Epetra_SerialComm.h"
74 typedef AbstractLinAlgPack::value_type Scalar;
78 int main(
int argc,
char* argv[] )
82 using Teuchos::rcp_dynamic_cast;
83 using Teuchos::rcp_implicit_cast;
86 using Teuchos::outArg;
104 const int procRank = mpiSession.getRank();
105 const int numProcs = mpiSession.getNProc();
109 bool result, success =
true;
111 RCP<Teuchos::FancyOStream>
120 MoochoThyraSolver solver;
126 CommandLineProcessor clp;
127 clp.throwExceptions(
false);
128 clp.addOutputSetupOptions(
true);
132 solver.setupCLP(&clp);
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", ¶mInvError,
169 "The error in the l2 norm of the parameter inverse solution error." );
171 CommandLineProcessor::EParseCommandLineReturn
172 parse_return = clp.parse(argc,argv,&std::cerr);
174 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
178 solver.readParameters(out.get());
182 <<
"\n*** NLPThyraEpetraAdvDiffReactOptMain, Global numProcs = "<<numProcs
185 int clusterRank = -1;
186 int numClusters = -1;
190 RCP<OpaqueWrapper<MPI_Comm> >
191 intraClusterMpiComm = Teuchos::opaqueWrapper<MPI_Comm>(MPI_COMM_WORLD),
192 interClusterMpiComm = Teuchos::null;
195 if ( numProcsPerCluster <= 0 ) {
197 <<
"\nnumProcsPerCluster = " << numProcsPerCluster
198 <<
" <= 0: Setting to " << numProcs <<
"...\n";
199 numProcsPerCluster = numProcs;
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<<
"!"
212 clusterRank = procRank / numProcsPerCluster;
213 *out <<
"\nclusterRank = " << clusterRank <<
"\n";
214 const int firstClusterProcRank = clusterRank * numProcsPerCluster;
215 const int lastClusterProcRank = firstClusterProcRank + numProcsPerCluster - 1;
216 *out <<
"\nclusterProcRange = ["<<firstClusterProcRank<<
","<<lastClusterProcRank<<
"]\n";
218 *out <<
"\nCreating intraClusterMpiComm ...";
219 MPI_Comm rawIntraClusterMpiComm = MPI_COMM_NULL;
224 ,&rawIntraClusterMpiComm
226 intraClusterMpiComm = Teuchos::opaqueWrapper(rawIntraClusterMpiComm,MPI_Comm_free);
228 *out <<
"\nintraClusterMpiComm:";
231 MPI_Comm_size(*intraClusterMpiComm,&size);
232 MPI_Comm_rank(*intraClusterMpiComm,&rank);
233 *out <<
"\nsize="<<size;
234 *out <<
"\nrank="<<rank;
238 *out <<
"\nCreating interClusterMpiComm ...";
239 MPI_Comm rawInterClusterMpiComm = MPI_COMM_NULL;
242 ,procRank==firstClusterProcRank?0:MPI_UNDEFINED
244 ,&rawInterClusterMpiComm
246 if(rawInterClusterMpiComm!=MPI_COMM_NULL)
247 interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm,MPI_Comm_free);
249 interClusterMpiComm = Teuchos::opaqueWrapper(rawInterClusterMpiComm);
251 *out <<
"\ninterClusterMpiComm:";
253 if(*interClusterMpiComm==MPI_COMM_NULL) {
258 MPI_Comm_size(*interClusterMpiComm,&size);
259 MPI_Comm_rank(*interClusterMpiComm,&rank);
260 *out <<
"\nsize="<<size;
261 *out <<
"\nrank="<<rank;
269 RCP<Epetra_Comm> comm = Teuchos::null;
272 Teuchos::set_extra_data(intraClusterMpiComm,
"mpiComm", outArg(comm));
281 *out <<
"\nCreate the GLpApp::AdvDiffReactOptModel wrapper object ...\n";
283 RCP<GLpApp::AdvDiffReactOptModel>
284 epetraModel = epetraModelCreator.
createModel(comm);
286 *out <<
"\nCreate the Thyra::LinearOpWithSolveFactory object ...\n";
288 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
292 *out <<
"\nCreate the Thyra::EpetraModelEvaluator wrapper object ...\n";
294 RCP<Thyra::EpetraModelEvaluator>
296 epetraThyraModel->setOStream(out);
297 epetraThyraModel->initialize(epetraModel,lowsFactory);
300 <<
"\nnx = " << epetraThyraModel->get_x_space()->dim()
301 <<
"\nnp = " << epetraThyraModel->get_p_space(0)->dim() <<
"\n";
307 RCP<const Thyra::ProductVectorSpaceBase<Scalar> >
308 x_bar_space, f_bar_space;
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);
329 epetraThyraModel->get_x_space()
340 epetraThyraModel->get_f_space()
347 vectorSpaceTester.
dump_all(dumpAll);
349 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
350 RTOpPack::show_mpi_apply_op_dump = dumpAll;
352 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
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;
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;
364 RCP<const VectorBase<Scalar> >
365 x0 = epetraThyraModel->getNominalValues().get_x();
398 #ifdef RTOPPACK_SPMD_APPLY_OP_DUMP
399 RTOpPack::show_mpi_apply_op_dump =
false;
401 #ifdef THYRA_SPMD_VECTOR_BASE_DUMP
412 *out <<
"\nEnd Result: TEST PASSED" << std::endl;
414 *out <<
"\nEnd Result: TEST FAILED" << std::endl;
416 return ( success ? 0 : 1 );
420 const int N = numPeriodsPerCluster;
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";
430 *out <<
"\nBuilding multiple Thyra::DefaultInverseModelEvaluator objects where the matching vector is held internally ...\n";
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(),
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>(
452 _inverseThyraModel->setParameterList(invMEPL);
453 _inverseThyraModel->set_observationMatchWeightingOp(H);
454 _inverseThyraModel->set_parameterRegularizationWeightingOp(R_bar);
455 inverseThyraModels[i] = _inverseThyraModel;
461 inverseThyraModels[i] = inverseThyraModels[0];
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;
474 const int p_index = 0;
475 const int z_index = 1;
476 const int z_p_index = 1;
477 const int z_x_index = 2;
478 Array<int> z_indexes;
479 if (useStatelessPeriodModel)
480 z_indexes = tuple<int>(z_p_index, z_x_index);
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();
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 );
495 if ( useStatelessPeriodModel ) {
497 tuple<RCP<
const VectorBase<Scalar> > >(
505 tuple<RCP<
const VectorBase<Scalar> > >(z_i)
508 scale_z_i *= periodParamScale;
511 RCP<Thyra::ModelEvaluator<Scalar> >
512 thyraModel = inverseThyraModels[0];
518 N, inverseThyraModels, z_indexes, z, g_index, weights,
519 x_bar_space, f_bar_space
524 MoochoSolver::ESolutionStatus solution_status;
527 *out <<
"\n***\n*** Solving the initial forward problem\n***\n";
531 solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_FORWARD);
534 RCP<const VectorBase<Scalar> >
537 p_opt = inverseThyraModels[0]->getNominalValues().get_p(0)->clone_v();
541 if ( initialSolveContinuation ) {
543 *out <<
"\nSolving individual period problems one at time using continuation ...\n";
545 RCP<ProductVectorBase<Scalar> >
546 x_opt_prod = rcp_dynamic_cast<ProductVectorBase<Scalar> >(
547 createMember( thyraModel->get_x_space() ),
true );
549 RCP<const VectorBase<Scalar> > period_x;
551 for (
int i = 0; i < N; ++i ) {
553 *out <<
"\nSolving period i = " << i <<
" using guess from last period ...\n";
556 solver.getSolver().set_output_context(
"fwd-init-"+
toString(i));
559 solver.setModel(inverseThyraModels[i]);
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());
571 initialGuess.set_x(period_x);
573 solver.setInitialGuess(initialGuess);
576 solution_status = solver.solve();
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();
588 x_init = x_opt->clone_v();
590 if ( useStatelessPeriodModel ) {
599 *out <<
"\nSolving all periods simultaniously ...\n";
602 solver.getSolver().set_output_context(
"fwd-init");
605 solver.setModel(thyraModel);
608 solver.readInitialGuess(out.get());
611 solution_status = solver.solve();
615 x_opt = solver.getFinalPoint().get_x()->clone_v();
616 x_init = solver.getFinalPoint().get_x()->clone_v();
621 *out <<
"\n***\n*** Solving the perturbed forward problem\n***\n";
625 solver.getSolver().set_output_context(
"fwd");
628 solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_FORWARD);
631 solver.setModel(thyraModel);
634 RCP<VectorBase<Scalar> >
635 p_init = p_opt->clone_v();
637 MEB::InArgs<Scalar> initialGuess = thyraModel->createInArgs();
638 initialGuess.setArgs(thyraModel->getNominalValues());
639 initialGuess.set_x(x_init);
641 initialGuess.set_p(0,p_init);
643 solver.setInitialGuess(initialGuess);
647 solution_status = solver.solve();
650 x_init = solver.getFinalPoint().get_x()->clone_v();
651 p_init = solver.getFinalPoint().get_p(0)->clone_v();
659 *out <<
"\n***\n*** Solving the perturbed inverse problem\n***\n";
663 solver.getSolver().set_output_context(
"inv");
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
677 for (
int i = 0; i < N; ++i ) {
679 inverseThyraModels[i], true
680 )->set_observationTarget(x_opt_prod->getVectorBlock(i));
684 RCP<const ProductVectorBase<Scalar> >
685 x_opt_prod = rcp_dynamic_cast<
const ProductVectorBase<Scalar> >(
686 rcp_implicit_cast<
const VectorBase<Scalar> >(x_opt)
688 RCP<const VectorBase<Scalar> >
691 ? x_opt_prod->getVectorBlock(0)
692 : rcp_implicit_cast<
const VectorBase<Scalar> >(x_opt)
694 if (useOuterInverse) {
696 inverseThyraModels[0], true
697 )->set_observationTarget(x_opt_i);
710 solver.setSolveMode(MoochoThyraSolver::SOLVE_MODE_OPTIMIZE);
713 solver.setModel(thyraModel);
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);
722 solver.setInitialGuess(initialGuess);
726 solution_status = solver.solve();
730 *out <<
"\n***\n*** Testing the error in the inversion\n***\n";
735 RCP<const VectorBase<Scalar> >
736 x_inv = solver.getFinalPoint().get_x(),
737 p_inv = solver.getFinalPoint().get_p(0);
747 x_test_passed = ( x_err <= stateInvError ),
748 p_test_passed = ( p_err <= paramInvError );
751 <<
"\nrelVectorErr(x_inv,x_opt) = " << x_err <<
" <= " << stateInvError
753 <<
"\nrelVectorErr(p_inv,p_opt) = " << p_err <<
" <= " << paramInvError
758 solver.writeFinalSolution(out.get());
762 solver.writeParamsFile();
765 solution_status != MoochoSolver::SOLVE_RETURN_SOLVED
772 *out << "\nEnd Result: TEST PASSED" << std::endl;
774 *out << "\nEnd Result: TEST FAILED" << std::endl;
776 return ( success ? 0 : 1 );
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)
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)