244 int numStates = solutionHistory->getNumStates();
246 TEUCHOS_TEST_FOR_EXCEPTION(
247 numStates < 1, std::logic_error,
248 "Error - setInitialConditions() needs at least one SolutionState\n"
249 " to set the initial condition. Number of States = "
253 RCP<Teuchos::FancyOStream> out = this->getOStream();
254 out->setOutputToRootOnly(0);
255 Teuchos::OSTab ostab(out, 1,
256 "StepperNewmarkImplicitAForm::setInitialConditions()");
257 *out <<
"Warning -- SolutionHistory has more than one state!\n"
258 <<
"Setting the initial conditions on the currentState.\n"
262 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
263 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
264 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
266 auto inArgs = this->wrapperModel_->getNominalValues();
267 TEUCHOS_TEST_FOR_EXCEPTION(
268 !((x != Teuchos::null && xDot != Teuchos::null) ||
269 (inArgs.get_x() != Teuchos::null &&
270 inArgs.get_x_dot() != Teuchos::null)),
272 "Error - We need to set the initial conditions for x and xDot from\n"
273 " either initialState or appModel_->getNominalValues::InArgs\n"
274 " (but not from a mixture of the two).\n");
277 if (x == Teuchos::null || xDot == Teuchos::null) {
278 using Teuchos::rcp_const_cast;
279 TEUCHOS_TEST_FOR_EXCEPTION(
280 (inArgs.get_x() == Teuchos::null) ||
281 (inArgs.get_x_dot() == Teuchos::null),
283 "Error - setInitialConditions() needs the ICs from the initialState\n"
284 " or getNominalValues()!\n");
285 x = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
286 initialState->setX(x);
287 xDot = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
288 initialState->setXDot(xDot);
292 if (initialState->getXDotDot() == Teuchos::null)
293 initialState->setXDotDot(initialState->getX()->clone_v());
295 this->setStepperXDotDot(initialState->getXDotDot());
298 std::string icConsistency = this->getICConsistency();
299 if (icConsistency ==
"None") {
300 if (initialState->getXDotDot() == Teuchos::null) {
301 RCP<Teuchos::FancyOStream> out = this->getOStream();
302 out->setOutputToRootOnly(0);
303 Teuchos::OSTab ostab(
304 out, 1,
"StepperNewmarkImplicitAForm::setInitialConditions()");
305 *out <<
"Warning -- Requested IC consistency of 'None' but\n"
306 <<
" initialState does not have an xDot.\n"
307 <<
" Setting a 'Zero' xDot!\n"
310 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
313 else if (icConsistency ==
"Zero")
314 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
315 else if (icConsistency ==
"App") {
316 auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
317 inArgs.get_x_dot_dot());
318 TEUCHOS_TEST_FOR_EXCEPTION(
319 xDotDot == Teuchos::null, std::logic_error,
320 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
321 " but 'App' returned a null pointer for xDotDot!\n");
322 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
324 else if (icConsistency ==
"Consistent") {
326 const Scalar time = initialState->getTime();
327 auto xDotDot = this->getStepperXDotDot(initialState);
331 if (this->initialGuess_ != Teuchos::null) {
332 TEUCHOS_TEST_FOR_EXCEPTION(
333 !((xDotDot->space())->isCompatible(*this->initialGuess_->space())),
335 "Error - User-provided initial guess for Newton is not compatible\n"
336 " with solution vector!\n");
337 Thyra::copy(*this->initialGuess_, xDotDot.ptr());
340 Thyra::put_scalar(0.0, xDotDot.ptr());
344 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
345 this->wrapperModel_);
347 wrapperModel->initializeNewmark(xDot, x, 0.0, time, beta_, gamma_);
348 const Thyra::SolveStatus<Scalar> sStatus =
349 (*(this->solver_)).solve(&*xDotDot);
351 TEUCHOS_TEST_FOR_EXCEPTION(
352 sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED, std::logic_error,
353 "Error - Solver failed while determining the initial conditions.\n"
355 << Thyra::toString(sStatus.solveStatus) <<
".\n");
358 TEUCHOS_TEST_FOR_EXCEPTION(
359 true, std::logic_error,
360 "Error - setInitialConditions() invalid IC consistency, "
361 << icConsistency <<
".\n");
366 initialState->setIsSynced(
true);
369 if (this->getICConsistencyCheck()) {
370 auto f = initialState->getX()->clone_v();
371 auto xDotDot = this->getStepperXDotDot(initialState);
373 typedef Thyra::ModelEvaluatorBase MEB;
374 MEB::InArgs<Scalar> appInArgs =
375 this->wrapperModel_->getAppModel()->createInArgs();
376 MEB::OutArgs<Scalar> appOutArgs =
377 this->wrapperModel_->getAppModel()->createOutArgs();
380 appInArgs.set_x_dot(xDot);
381 appInArgs.set_x_dot_dot(xDotDot);
383 appOutArgs.set_f(appOutArgs.get_f());
385 appInArgs.set_W_x_dot_dot_coeff(Scalar(0.0));
386 appInArgs.set_alpha(Scalar(0.0));
387 appInArgs.set_beta(Scalar(0.0));
389 appInArgs.set_t(initialState->getTime());
391 this->wrapperModel_->getAppModel()->evalModel(appInArgs, appOutArgs);
393 Scalar reldiff = Thyra::norm(*f);
394 Scalar normx = Thyra::norm(*x);
395 Scalar eps = Scalar(100.0) * std::abs(Teuchos::ScalarTraits<Scalar>::eps());
396 if (normx > eps * reldiff) reldiff /= normx;
399 RCP<Teuchos::FancyOStream> out = this->getOStream();
400 out->setOutputToRootOnly(0);
401 Teuchos::OSTab ostab(
402 out, 1,
"StepperNewmarkImplicitAForm::setInitialConditions()");
403 *out <<
"Warning -- Failed consistency check but continuing!\n"
404 <<
" ||f(x,xDot,xDotDot,t)||/||x|| > eps" << std::endl
405 <<
" ||f(x,xDot,xDotDot,t)|| = " << Thyra::norm(*f)
407 <<
" ||x|| = " << Thyra::norm(*x)
409 <<
" ||f(x,xDot,xDotDot,t)||/||x|| = " << reldiff << std::endl
410 <<
" eps = " << eps << std::endl;
414 if (!(this->getUseFSAL())) {
415 RCP<Teuchos::FancyOStream> out = this->getOStream();
416 out->setOutputToRootOnly(0);
417 Teuchos::OSTab ostab(out, 1,
418 "StepperNewmarkImplicitAForm::setInitialConditions()");
419 *out <<
"\nWarning -- The First-Same-As-Last (FSAL) principle is "
420 <<
"part of the Newmark Implicit A-Form. The default is to "
421 <<
"set useFSAL=true, and useFSAL=false will be ignored." << std::endl;
429#ifdef VERBOSE_DEBUG_OUTPUT
430 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
432 this->checkInitialized();
436 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperNewmarkImplicitAForm::takeStep()");
438 TEUCHOS_TEST_FOR_EXCEPTION(
439 solutionHistory->getNumStates() < 2, std::logic_error,
440 "Error - StepperNewmarkImplicitAForm<Scalar>::takeStep(...)\n"
441 <<
"Need at least two SolutionStates for NewmarkImplicitAForm.\n"
442 <<
" Number of States = " << solutionHistory->getNumStates()
443 <<
"\nTry setting in \"Solution History\" \"Storage Type\" = "
445 <<
" or \"Storage Type\" = \"Static\" and \"Storage Limit\" = "
448 auto thisStepper = Teuchos::rcpFromRef(*
this);
449 stepperNewmarkImpAppAction_->execute(
450 solutionHistory, thisStepper,
452 Scalar>::ACTION_LOCATION::BEGIN_STEP);
454 RCP<SolutionState<Scalar> > workingState =
455 solutionHistory->getWorkingState();
456 RCP<SolutionState<Scalar> > currentState =
457 solutionHistory->getCurrentState();
460 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
461 this->wrapperModel_);
464 RCP<const Thyra::VectorBase<Scalar> > d_old = currentState->getX();
465 RCP<const Thyra::VectorBase<Scalar> > v_old = currentState->getXDot();
466 RCP<Thyra::VectorBase<Scalar> > a_old = currentState->getXDotDot();
469 RCP<Thyra::VectorBase<Scalar> > d_new = workingState->getX();
470 RCP<Thyra::VectorBase<Scalar> > v_new = workingState->getXDot();
471 RCP<Thyra::VectorBase<Scalar> > a_new = workingState->getXDotDot();
474 const Scalar time = currentState->getTime();
475 const Scalar dt = workingState->getTimeStep();
476 Scalar t = time + dt;
479 predictDisplacement(*d_new, *d_old, *v_old, *a_old, dt);
480 predictVelocity(*v_new, *v_old, *a_old, dt);
483 wrapperModel->initializeNewmark(v_new, d_new, dt, t, beta_, gamma_);
485 stepperNewmarkImpAppAction_->execute(
486 solutionHistory, thisStepper,
488 Scalar>::ACTION_LOCATION::BEFORE_SOLVE);
490 if (this->getZeroInitialGuess())
491 Thyra::assign(a_new.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
494 const Thyra::SolveStatus<Scalar> sStatus =
495 (*(this->solver_)).solve(&*a_new);
497 stepperNewmarkImpAppAction_->execute(
498 solutionHistory, thisStepper,
500 Scalar>::ACTION_LOCATION::AFTER_SOLVE);
503 correctVelocity(*v_new, *v_new, *a_new, dt);
504 correctDisplacement(*d_new, *d_new, *a_new, dt);
506 workingState->setSolutionStatus(sStatus);
507 workingState->setOrder(this->getOrder());
508 workingState->computeNorms(currentState);
510 stepperNewmarkImpAppAction_->execute(
511 solutionHistory, thisStepper,
513 Scalar>::ACTION_LOCATION::END_STEP);
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...