Rythmos - Transient Integration for Differential Equations Version of the Day
Loading...
Searching...
No Matches
diagonalTransientMain.cpp
1//@HEADER
2
3// ***********************************************************************
4//
5// Rythmos Package
6// Copyright (2006) 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// This library is free software; you can redistribute it and/or modify
12// it under the terms of the GNU Lesser General Public License as
13// published by the Free Software Foundation; either version 2.1 of the
14// License, or (at your option) any later version.
15//
16// This library is distributed in the hope that it will be useful, but
17// WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19// Lesser General Public License for more details.
20//
21// You should have received a copy of the GNU Lesser General Public
22// License along with this library; if not, write to the Free Software
23// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
24// USA
25// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
26//
27// ***********************************************************************
28//@HEADER
29
30#include "EpetraExt_DiagonalTransientModel.hpp"
31#include "Rythmos_BackwardEulerStepper.hpp"
32#include "Rythmos_ImplicitBDFStepper.hpp"
33#include "Rythmos_ImplicitRKStepper.hpp"
34#include "Rythmos_ForwardSensitivityStepper.hpp"
35#include "Rythmos_TimeStepNonlinearSolver.hpp"
36#include "Rythmos_DefaultIntegrator.hpp"
37#include "Rythmos_SimpleIntegrationControlStrategy.hpp"
38#include "Rythmos_StepperAsModelEvaluator.hpp"
39#include "Rythmos_RKButcherTableauBuilder.hpp"
40#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
41#include "Thyra_EpetraModelEvaluator.hpp"
42#include "Thyra_DirectionalFiniteDiffCalculator.hpp"
43#include "Thyra_TestingTools.hpp"
44#include "Teuchos_StandardCatchMacros.hpp"
45#include "Teuchos_GlobalMPISession.hpp"
46#include "Teuchos_DefaultComm.hpp"
47#include "Teuchos_VerboseObject.hpp"
48#include "Teuchos_XMLParameterListHelpers.hpp"
49#include "Teuchos_CommandLineProcessor.hpp"
50#include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
51#include "Teuchos_as.hpp"
52
53#ifdef HAVE_MPI
54# include "Epetra_MpiComm.h"
55#else
56# include "Epetra_SerialComm.h"
57#endif // HAVE_MPI
58
59
60namespace {
61
62
63const std::string TimeStepNonlinearSolver_name = "TimeStepNonlinearSolver";
64
65const std::string Stratimikos_name = "Stratimikos";
66
67const std::string DiagonalTransientModel_name = "DiagonalTransientModel";
68
69const std::string RythmosStepper_name = "Rythmos Stepper";
70
71const std::string RythmosIntegrator_name = "Rythmos Integrator";
72
73const std::string FdCalc_name = "FD Calc";
74
75
76Teuchos::RCP<const Teuchos::ParameterList>
77getValidParameters()
78{
79 using Teuchos::RCP; using Teuchos::ParameterList;
80 static RCP<const ParameterList> validPL;
81 if (is_null(validPL)) {
82 RCP<ParameterList> pl = Teuchos::parameterList();
83 pl->sublist(TimeStepNonlinearSolver_name);
84 pl->sublist(Stratimikos_name);
85 pl->sublist(DiagonalTransientModel_name);
86 pl->sublist(RythmosStepper_name);
87 pl->sublist(RythmosIntegrator_name);
88 pl->sublist(FdCalc_name);
89 validPL = pl;
90 }
91 return validPL;
92}
93
94
95} // namespace
96
97
98int main(int argc, char *argv[])
99{
100
101 using std::endl;
102 typedef double Scalar;
103 // typedef double ScalarMag; // unused
104 using Teuchos::describe;
105 using Teuchos::RCP;
106 using Teuchos::rcp;
107 using Teuchos::rcp_implicit_cast;
108 using Teuchos::rcp_dynamic_cast;
109 using Teuchos::as;
110 using Teuchos::ParameterList;
111 using Teuchos::CommandLineProcessor;
112 typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
113 typedef Thyra::ModelEvaluatorBase MEB;
114 // typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS; // unused
115 using Thyra::productVectorBase;
116
117 bool result, success = true;
118
119 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
120
121 RCP<Epetra_Comm> epetra_comm;
122#ifdef HAVE_MPI
123 epetra_comm = rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
124#else
125 epetra_comm = rcp( new Epetra_SerialComm );
126#endif // HAVE_MPI
127
128 RCP<Teuchos::FancyOStream>
129 out = Teuchos::VerboseObjectBase::getDefaultOStream();
130
131 try {
132
133 //
134 // Read commandline options
135 //
136
137 CommandLineProcessor clp;
138 clp.throwExceptions(false);
139 clp.addOutputSetupOptions(true);
140
141 std::string paramsFileName = "";
142 clp.setOption( "params-file", &paramsFileName,
143 "File name for XML parameters" );
144
145 std::string extraParamsString = "";
146 clp.setOption( "extra-params", &extraParamsString,
147 "Extra XML parameters" );
148
149 std::string extraParamsFile = "";
150 clp.setOption( "extra-params-file", &extraParamsFile, "File containing extra parameters in XML format.");
151
152 double maxStateError = 1e-6;
153 clp.setOption( "max-state-error", &maxStateError,
154 "The maximum allowed error in the integrated state in relation to the exact state solution" );
155
156 double finalTime = 1e-3;
157 clp.setOption( "final-time", &finalTime,
158 "Final integration time (initial time is 0.0)" );
159
160 int numTimeSteps = 10;
161 clp.setOption( "num-time-steps", &numTimeSteps,
162 "Number of (fixed) time steps. If <= 0.0, then variable time steps are taken" );
163
164 bool useBDF = false;
165 clp.setOption( "use-BDF", "use-BE", &useBDF,
166 "Use BDF or Backward Euler (BE)" );
167
168 bool useIRK = false;
169 clp.setOption( "use-IRK", "use-other", &useIRK,
170 "Use IRK or something" );
171
172 bool doFwdSensSolve = false;
173 clp.setOption( "fwd-sens-solve", "state-solve", &doFwdSensSolve,
174 "Do the forward sensitivity solve or just the state solve" );
175
176 bool doFwdSensErrorControl = false;
177 clp.setOption( "fwd-sens-err-cntrl", "no-fwd-sens-err-cntrl", &doFwdSensErrorControl,
178 "Do error control on the forward sensitivity solve or not" );
179
180 double maxRestateError = 0.0;
181 clp.setOption( "max-restate-error", &maxRestateError,
182 "The maximum allowed error between the state integrated by itself verses integrated along with DxDp" );
183
184 double maxSensError = 1e-4;
185 clp.setOption( "max-sens-error", &maxSensError,
186 "The maximum allowed error in the integrated sensitivity in relation to"
187 " the finite-difference sensitivity" );
188
189 Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
190 setVerbosityLevelOption( "verb-level", &verbLevel,
191 "Top-level verbosity level. By default, this gets deincremented as you go deeper into numerical objects.",
192 &clp );
193
194 bool testExactSensitivity = false;
195 clp.setOption( "test-exact-sens", "no-test-exact-sens", &testExactSensitivity,
196 "Test the exact sensitivity with finite differences or not." );
197
198 bool dumpFinalSolutions = false;
199 clp.setOption(
200 "dump-final-solutions", "no-dump-final-solutions", &dumpFinalSolutions,
201 "Determine if the final solutions are dumpped or not." );
202
203 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
204 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
205
206 if ( Teuchos::VERB_DEFAULT == verbLevel )
207 verbLevel = Teuchos::VERB_LOW;
208
209 const Teuchos::EVerbosityLevel
210 solnVerbLevel = ( dumpFinalSolutions ? Teuchos::VERB_EXTREME : verbLevel );
211
212 //
213 // Get the base parameter list that all other parameter lists will be read
214 // from.
215 //
216
217 RCP<ParameterList>
218 paramList = Teuchos::parameterList();
219 if (paramsFileName.length())
220 updateParametersFromXmlFile( paramsFileName, paramList.ptr() );
221 if(extraParamsFile.length())
222 Teuchos::updateParametersFromXmlFile( "./"+extraParamsFile, paramList.ptr() );
223 if (extraParamsString.length())
224 updateParametersFromXmlString( extraParamsString, paramList.ptr() );
225
226 if (testExactSensitivity) {
227 paramList->sublist(DiagonalTransientModel_name).set("Exact Solution as Response",true);
228 }
229
230 paramList->validateParameters(*getValidParameters(),0); // Only validate top level lists!
231
232 //
233 // Create the Stratimikos linear solver factory.
234 //
235 // This is the linear solve strategy that will be used to solve for the
236 // linear system with the W.
237 //
238
239 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
240 linearSolverBuilder.setParameterList(sublist(paramList,Stratimikos_name));
241 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
242 W_factory = createLinearSolveStrategy(linearSolverBuilder);
243
244 //
245 // Create the underlying EpetraExt::ModelEvaluator
246 //
247
248 RCP<EpetraExt::DiagonalTransientModel>
249 epetraStateModel = EpetraExt::diagonalTransientModel(
250 epetra_comm,
251 sublist(paramList,DiagonalTransientModel_name)
252 );
253
254 *out <<"\nepetraStateModel valid options:\n";
255 epetraStateModel->getValidParameters()->print(
256 *out, PLPrintOptions().indent(2).showTypes(true).showDoc(true)
257 );
258
259 //
260 // Create the Thyra-wrapped ModelEvaluator
261 //
262
263 RCP<Thyra::ModelEvaluator<double> >
264 stateModel = epetraModelEvaluator(epetraStateModel,W_factory);
265
266 *out << "\nParameter names = " << *stateModel->get_p_names(0) << "\n";
267
268 //
269 // Create the Rythmos stateStepper
270 //
271
272 RCP<Rythmos::TimeStepNonlinearSolver<double> >
273 nonlinearSolver = Rythmos::timeStepNonlinearSolver<double>();
274 RCP<ParameterList>
275 nonlinearSolverPL = sublist(paramList,TimeStepNonlinearSolver_name);
276 nonlinearSolverPL->get("Default Tol",1e-3*maxStateError); // Set default if not set
277 nonlinearSolver->setParameterList(nonlinearSolverPL);
278
279 RCP<Rythmos::StepperBase<Scalar> > stateStepper;
280
281 if (useBDF) {
282 stateStepper = rcp(
283 new Rythmos::ImplicitBDFStepper<double>(
284 stateModel, nonlinearSolver
285 )
286 );
287 }
288 else if (useIRK) {
289 // We need a separate LOWSFB object for the IRK stepper
290 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
291 irk_W_factory = createLinearSolveStrategy(linearSolverBuilder);
292 RCP<Rythmos::RKButcherTableauBase<double> > irkbt = Rythmos::createRKBT<double>("Backward Euler");
293 stateStepper = Rythmos::implicitRKStepper<double>(
294 stateModel, nonlinearSolver, irk_W_factory, irkbt
295 );
296 }
297 else {
298 stateStepper = rcp(
299 new Rythmos::BackwardEulerStepper<double>(
300 stateModel, nonlinearSolver
301 )
302 );
303 }
304
305 *out <<"\nstateStepper:\n" << describe(*stateStepper,verbLevel);
306 *out <<"\nstateStepper valid options:\n";
307 stateStepper->getValidParameters()->print(
308 *out, PLPrintOptions().indent(2).showTypes(true).showDoc(true)
309 );
310
311 stateStepper->setParameterList(sublist(paramList,RythmosStepper_name));
312
313 //
314 // Setup finite difference objects that will be used for tests
315 //
316
317 Thyra::DirectionalFiniteDiffCalculator<Scalar> fdCalc;
318 fdCalc.setParameterList(sublist(paramList,FdCalc_name));
319 fdCalc.setOStream(out);
320 fdCalc.setVerbLevel(verbLevel);
321
322 //
323 // Use a StepperAsModelEvaluator to integrate the state
324 //
325
326 const MEB::InArgs<Scalar>
327 state_ic = stateModel->getNominalValues();
328 *out << "\nstate_ic:\n" << describe(state_ic,verbLevel);
329
330 RCP<Rythmos::IntegratorBase<Scalar> > integrator;
331 {
332 RCP<ParameterList>
333 integratorPL = sublist(paramList,RythmosIntegrator_name);
334 integratorPL->set( "Take Variable Steps", as<bool>(numTimeSteps < 0) );
335 integratorPL->set( "Fixed dt", as<double>((finalTime - state_ic.get_t())/numTimeSteps) );
336 RCP<Rythmos::IntegratorBase<Scalar> >
337 defaultIntegrator = Rythmos::controlledDefaultIntegrator<Scalar>(
338 Rythmos::simpleIntegrationControlStrategy<Scalar>(integratorPL)
339 );
340 integrator = defaultIntegrator;
341 }
342
343 RCP<Rythmos::StepperAsModelEvaluator<Scalar> >
344 stateIntegratorAsModel = Rythmos::stepperAsModelEvaluator(
345 stateStepper, integrator, state_ic
346 );
347 stateIntegratorAsModel->setVerbLevel(verbLevel);
348
349 *out << "\nUse the StepperAsModelEvaluator to integrate state x(p,finalTime) ... \n";
350
351 RCP<Thyra::VectorBase<Scalar> > x_final;
352
353 {
354
355 Teuchos::OSTab tab(out);
356
357 x_final = createMember(stateIntegratorAsModel->get_g_space(0));
358
359 eval_g(
360 *stateIntegratorAsModel,
361 0, *state_ic.get_p(0),
362 finalTime,
363 0, &*x_final
364 );
365
366 *out
367 << "\nx_final = x(p,finalTime) evaluated using stateIntegratorAsModel:\n"
368 << describe(*x_final,solnVerbLevel);
369
370 }
371
372 //
373 // Test the integrated state against the exact analytical state solution
374 //
375
376 RCP<const Thyra::VectorBase<Scalar> >
377 exact_x_final = create_Vector(
378 epetraStateModel->getExactSolution(finalTime),
379 stateModel->get_x_space()
380 );
381
382 result = Thyra::testRelNormDiffErr(
383 "exact_x_final", *exact_x_final, "x_final", *x_final,
384 "maxStateError", maxStateError, "warningTol", 1.0, // Don't warn
385 &*out, solnVerbLevel
386 );
387 if (!result) success = false;
388
389 //
390 // Solve and test the forward sensitivity computation
391 //
392
393 if (doFwdSensSolve) {
394
395 //
396 // Create the forward sensitivity stepper
397 //
398
399 RCP<Rythmos::ForwardSensitivityStepper<Scalar> > stateAndSensStepper =
400 Rythmos::forwardSensitivityStepper<Scalar>();
401 if (doFwdSensErrorControl) {
402 stateAndSensStepper->initializeDecoupledSteppers(
403 stateModel, 0, stateModel->getNominalValues(),
404 stateStepper, nonlinearSolver,
405 integrator->cloneIntegrator(), finalTime
406 );
407 }
408 else {
409 stateAndSensStepper->initializeSyncedSteppers(
410 stateModel, 0, stateModel->getNominalValues(),
411 stateStepper, nonlinearSolver
412 );
413 // The above call will result in stateStepper and nonlinearSolver being
414 // cloned. This helps to ensure consistency between the state and
415 // sensitivity computations!
416 }
417
418 //
419 // Set the initial condition for the state and forward sensitivities
420 //
421
422 RCP<Thyra::VectorBase<Scalar> > s_bar_init
423 = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
424 assign( s_bar_init.ptr(), 0.0 );
425 RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init
426 = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
427 assign( s_bar_dot_init.ptr(), 0.0 );
428 // Above, I believe that these are the correct initial conditions for
429 // s_bar and s_bar_dot given how the EpetraExt::DiagonalTransientModel
430 // is currently implemented!
431
432 RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
433 stateAndSensModel = stateAndSensStepper->getStateAndFwdSensModel();
434
435 MEB::InArgs<Scalar>
436 state_and_sens_ic = stateAndSensStepper->getModel()->createInArgs();
437
438 // Copy time, parameters etc.
439 state_and_sens_ic.setArgs(state_ic);
440 // Set initial condition for x_bar = [ x; s_bar ]
441 state_and_sens_ic.set_x(
442 stateAndSensModel->create_x_bar_vec(state_ic.get_x(),s_bar_init)
443 );
444 // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
445 state_and_sens_ic.set_x_dot(
446 stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(),s_bar_dot_init)
447 );
448
449 *out << "\nstate_and_sens_ic:\n" << describe(state_and_sens_ic,verbLevel);
450
451 stateAndSensStepper->setInitialCondition(state_and_sens_ic);
452
453 //
454 // Use a StepperAsModelEvaluator to integrate the state+sens
455 //
456
457 RCP<Rythmos::StepperAsModelEvaluator<Scalar> >
458 stateAndSensIntegratorAsModel = Rythmos::stepperAsModelEvaluator(
459 rcp_implicit_cast<Rythmos::StepperBase<Scalar> >(stateAndSensStepper),
460 integrator, state_and_sens_ic
461 );
462 stateAndSensIntegratorAsModel->setVerbLevel(verbLevel);
463
464 *out << "\nUse the StepperAsModelEvaluator to integrate state + sens x_bar(p,finalTime) ... \n";
465
466 RCP<Thyra::VectorBase<Scalar> > x_bar_final;
467
468 {
469
470 Teuchos::OSTab tab(out);
471
472 x_bar_final = createMember(stateAndSensIntegratorAsModel->get_g_space(0));
473
474 eval_g(
475 *stateAndSensIntegratorAsModel,
476 0, *state_ic.get_p(0),
477 finalTime,
478 0, &*x_bar_final
479 );
480
481 *out
482 << "\nx_bar_final = x_bar(p,finalTime) evaluated using stateAndSensIntegratorAsModel:\n"
483 << describe(*x_bar_final,solnVerbLevel);
484
485 }
486
487 //
488 // Test that the state computed above is same as computed initially!
489 //
490
491 *out << "\nChecking that x(p,finalTime) computed as part of x_bar above is the same ...\n";
492
493 {
494
495 Teuchos::OSTab tab(out);
496
497 RCP<const Thyra::VectorBase<Scalar> >
498 x_in_x_bar_final = productVectorBase<Scalar>(x_bar_final)->getVectorBlock(0);
499
500 result = Thyra::testRelNormDiffErr<Scalar>(
501 "x_final", *x_final,
502 "x_in_x_bar_final", *x_in_x_bar_final,
503 "maxRestateError", maxRestateError,
504 "warningTol", 1.0, // Don't warn
505 &*out, solnVerbLevel
506 );
507 if (!result) success = false;
508
509 }
510
511 //
512 // Compute DxDp using finite differences
513 //
514
515 *out << "\nApproximating DxDp(p,t) using directional finite differences of integrator for x(p,t) ...\n";
516
517 RCP<Thyra::MultiVectorBase<Scalar> > DxDp_fd_final;
518
519 {
520
521 Teuchos::OSTab tab(out);
522
523
524 MEB::InArgs<Scalar>
525 fdBasePoint = stateIntegratorAsModel->createInArgs();
526
527 fdBasePoint.set_t(finalTime);
528 fdBasePoint.set_p(0,stateModel->getNominalValues().get_p(0));
529
530 DxDp_fd_final = createMembers(
531 stateIntegratorAsModel->get_g_space(0),
532 stateIntegratorAsModel->get_p_space(0)->dim()
533 );
534
535 typedef Thyra::DirectionalFiniteDiffCalculatorTypes::SelectedDerivatives
536 SelectedDerivatives;
537
538 MEB::OutArgs<Scalar> fdOutArgs =
539 fdCalc.createOutArgs(
540 *stateIntegratorAsModel,
541 SelectedDerivatives().supports(MEB::OUT_ARG_DgDp,0,0)
542 );
543 fdOutArgs.set_DgDp(0,0,DxDp_fd_final);
544
545 // Silence the model evaluators that are called. The fdCal object
546 // will show all of the inputs and outputs for each call.
547 stateStepper->setVerbLevel(Teuchos::VERB_NONE);
548 stateIntegratorAsModel->setVerbLevel(Teuchos::VERB_NONE);
549
550 fdCalc.calcDerivatives(
551 *stateIntegratorAsModel, fdBasePoint,
552 stateIntegratorAsModel->createOutArgs(), // Don't bother with function value
553 fdOutArgs
554 );
555
556 *out
557 << "\nFinite difference DxDp_fd_final = DxDp(p,finalTime): "
558 << describe(*DxDp_fd_final,solnVerbLevel);
559
560 }
561
562 //
563 // Test that the integrated sens and the F.D. sens are similar
564 //
565
566 *out << "\nChecking that integrated DxDp(p,finalTime) and finite-diff DxDp(p,finalTime) are similar ...\n";
567
568 {
569
570 Teuchos::OSTab tab(out);
571
572 RCP<const Thyra::VectorBase<Scalar> >
573 DxDp_vec_final = Thyra::productVectorBase<Scalar>(x_bar_final)->getVectorBlock(1);
574
575 RCP<const Thyra::VectorBase<Scalar> >
576 DxDp_fd_vec_final = Thyra::multiVectorProductVector(
577 rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> >(
578 DxDp_vec_final->range()
579 ),
580 DxDp_fd_final
581 );
582
583 result = Thyra::testRelNormDiffErr(
584 "DxDp_vec_final", *DxDp_vec_final,
585 "DxDp_fd_vec_final", *DxDp_fd_vec_final,
586 "maxSensError", maxSensError,
587 "warningTol", 1.0, // Don't warn
588 &*out, solnVerbLevel
589 );
590 if (!result) success = false;
591
592 }
593
594 }
595
596 }
597 TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
598
599 if(success)
600 *out << "\nEnd Result: TEST PASSED" << endl;
601 else
602 *out << "\nEnd Result: TEST FAILED" << endl;
603
604 return ( success ? 0 : 1 );
605
606} // end main() [Doxygen looks for this!]
607