256 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
260 int numStates = solutionHistory->getNumStates();
262 TEUCHOS_TEST_FOR_EXCEPTION(numStates < 1, std::logic_error,
263 "Error - setInitialConditions() needs at least one SolutionState\n"
264 " to set the initial condition. Number of States = " << numStates);
267 RCP<Teuchos::FancyOStream> out = this->getOStream();
268 out->setOutputToRootOnly(0);
269 Teuchos::OSTab ostab(out,1,
"StepperNewmarkImplicitAForm::setInitialConditions()");
270 *out <<
"Warning -- SolutionHistory has more than one state!\n"
271 <<
"Setting the initial conditions on the currentState.\n"<<std::endl;
274 RCP<SolutionState<Scalar> > initialState = solutionHistory->getCurrentState();
275 RCP<Thyra::VectorBase<Scalar> > x = initialState->getX();
276 RCP<Thyra::VectorBase<Scalar> > xDot = initialState->getXDot();
278 auto inArgs = this->wrapperModel_->getNominalValues();
279 TEUCHOS_TEST_FOR_EXCEPTION(
280 !((x != Teuchos::null && xDot != Teuchos::null) ||
281 (inArgs.get_x() != Teuchos::null &&
282 inArgs.get_x_dot() != Teuchos::null)), std::logic_error,
283 "Error - We need to set the initial conditions for x and xDot from\n"
284 " either initialState or appModel_->getNominalValues::InArgs\n"
285 " (but not from a mixture of the two).\n");
288 if ( x == Teuchos::null || xDot == Teuchos::null ) {
289 using Teuchos::rcp_const_cast;
290 TEUCHOS_TEST_FOR_EXCEPTION( (inArgs.get_x() == Teuchos::null) ||
291 (inArgs.get_x_dot() == Teuchos::null), std::logic_error,
292 "Error - setInitialConditions() needs the ICs from the initialState\n"
293 " or getNominalValues()!\n");
294 x = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x());
295 initialState->setX(x);
296 xDot = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
297 initialState->setXDot(xDot);
301 if (initialState->getXDotDot() == Teuchos::null)
302 initialState->setXDotDot(initialState->getX()->clone_v());
304 this->setStepperXDotDot(initialState->getXDotDot());
307 std::string icConsistency = this->getICConsistency();
308 if (icConsistency ==
"None") {
309 if (initialState->getXDotDot() == Teuchos::null) {
310 RCP<Teuchos::FancyOStream> out = this->getOStream();
311 out->setOutputToRootOnly(0);
312 Teuchos::OSTab ostab(out,1,
313 "StepperNewmarkImplicitAForm::setInitialConditions()");
314 *out <<
"Warning -- Requested IC consistency of 'None' but\n"
315 <<
" initialState does not have an xDot.\n"
316 <<
" Setting a 'Zero' xDot!\n" << std::endl;
318 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
321 else if (icConsistency ==
"Zero")
322 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), Scalar(0.0));
323 else if (icConsistency ==
"App") {
324 auto xDotDot = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar> >(
325 inArgs.get_x_dot_dot());
326 TEUCHOS_TEST_FOR_EXCEPTION(xDotDot == Teuchos::null, std::logic_error,
327 "Error - setInitialConditions() requested 'App' for IC consistency,\n"
328 " but 'App' returned a null pointer for xDotDot!\n");
329 Thyra::assign(this->getStepperXDotDot(initialState).ptr(), *xDotDot);
331 else if (icConsistency ==
"Consistent") {
333 const Scalar time = initialState->getTime();
334 auto xDotDot = this->getStepperXDotDot(initialState);
338 if (this->initialGuess_ != Teuchos::null) {
339 TEUCHOS_TEST_FOR_EXCEPTION(
340 !((xDotDot->space())->isCompatible(*this->initialGuess_->space())),
342 "Error - User-provided initial guess for Newton is not compatible\n"
343 " with solution vector!\n");
344 Thyra::copy(*this->initialGuess_, xDotDot.ptr());
347 Thyra::put_scalar(0.0, xDotDot.ptr());
351 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
352 this->wrapperModel_);
354 wrapperModel->initializeNewmark(xDot, x, 0.0, time, beta_, gamma_);
355 const Thyra::SolveStatus<Scalar> sStatus = (*(this->solver_)).solve(&*xDotDot);
357 TEUCHOS_TEST_FOR_EXCEPTION(
358 sStatus.solveStatus != Thyra::SOLVE_STATUS_CONVERGED, std::logic_error,
359 "Error - Solver failed while determining the initial conditions.\n"
360 " Solver status is "<<Thyra::toString(sStatus.solveStatus)<<
".\n");
363 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
364 "Error - setInitialConditions() invalid IC consistency, "
365 << icConsistency <<
".\n");
370 initialState->setIsSynced(
true);
373 if (this->getICConsistencyCheck()) {
374 auto f = initialState->getX()->clone_v();
375 auto xDotDot = this->getStepperXDotDot(initialState);
377 typedef Thyra::ModelEvaluatorBase MEB;
378 MEB::InArgs<Scalar> appInArgs =
379 this->wrapperModel_->getAppModel()->createInArgs();
380 MEB::OutArgs<Scalar> appOutArgs =
381 this->wrapperModel_->getAppModel()->createOutArgs();
383 appInArgs.set_x (x );
384 appInArgs.set_x_dot (xDot );
385 appInArgs.set_x_dot_dot(xDotDot);
387 appOutArgs.set_f(appOutArgs.get_f());
389 appInArgs.set_W_x_dot_dot_coeff(Scalar(0.0));
390 appInArgs.set_alpha (Scalar(0.0));
391 appInArgs.set_beta (Scalar(0.0));
393 appInArgs.set_t (initialState->getTime() );
395 this->wrapperModel_->getAppModel()->evalModel(appInArgs, appOutArgs);
397 Scalar reldiff = Thyra::norm(*f);
398 Scalar normx = Thyra::norm(*x);
399 Scalar eps = Scalar(100.0)*std::abs(Teuchos::ScalarTraits<Scalar>::eps());
400 if (normx > eps*reldiff) reldiff /= normx;
403 RCP<Teuchos::FancyOStream> out = this->getOStream();
404 out->setOutputToRootOnly(0);
405 Teuchos::OSTab ostab(out,1,
406 "StepperNewmarkImplicitAForm::setInitialConditions()");
407 *out <<
"Warning -- Failed consistency check but continuing!\n"
408 <<
" ||f(x,xDot,xDotDot,t)||/||x|| > eps" << std::endl
409 <<
" ||f(x,xDot,xDotDot,t)|| = " << Thyra::norm(*f)<< std::endl
410 <<
" ||x|| = " << Thyra::norm(*x)<< std::endl
411 <<
" ||f(x,xDot,xDotDot,t)||/||x|| = " << reldiff << std::endl
412 <<
" eps = " << eps << std::endl;
416 if (!(this->getUseFSAL())) {
417 RCP<Teuchos::FancyOStream> out = this->getOStream();
418 out->setOutputToRootOnly(0);
419 Teuchos::OSTab ostab(out,1,
420 "StepperNewmarkImplicitAForm::setInitialConditions()");
421 *out <<
"\nWarning -- The First-Same-As-Last (FSAL) principle is "
422 <<
"part of the Newmark Implicit A-Form. The default is to "
423 <<
"set useFSAL=true, and useFSAL=false will be ignored." << std::endl;
430 const Teuchos::RCP<SolutionHistory<Scalar> >& solutionHistory)
432#ifdef VERBOSE_DEBUG_OUTPUT
433 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
435 this->checkInitialized();
439 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperNewmarkImplicitAForm::takeStep()");
441 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
443 "Error - StepperNewmarkImplicitAForm<Scalar>::takeStep(...)\n"
444 "Need at least two SolutionStates for NewmarkImplicitAForm.\n"
445 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
446 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
447 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
449 auto thisStepper = Teuchos::rcpFromRef(*
this);
450 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,
453 RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
454 RCP<SolutionState<Scalar> > currentState=solutionHistory->getCurrentState();
457 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
458 this->wrapperModel_);
461 RCP<const Thyra::VectorBase<Scalar> > d_old = currentState->getX();
462 RCP<const Thyra::VectorBase<Scalar> > v_old = currentState->getXDot();
463 RCP< Thyra::VectorBase<Scalar> > a_old = currentState->getXDotDot();
466 RCP<Thyra::VectorBase<Scalar> > d_new = workingState->getX();
467 RCP<Thyra::VectorBase<Scalar> > v_new = workingState->getXDot();
468 RCP<Thyra::VectorBase<Scalar> > a_new = workingState->getXDotDot();
471 const Scalar time = currentState->getTime();
472 const Scalar dt = workingState->getTimeStep();
476 predictDisplacement(*d_new, *d_old, *v_old, *a_old, dt);
477 predictVelocity(*v_new, *v_old, *a_old, dt);
480 wrapperModel->initializeNewmark(v_new,d_new,dt,t,beta_,gamma_);
482 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,
485 if (this->getZeroInitialGuess())
486 Thyra::assign(a_new.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
489 const Thyra::SolveStatus<Scalar> sStatus = (*(this->solver_)).solve(&*a_new);
491 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,
495 correctVelocity(*v_new, *v_new, *a_new, dt);
496 correctDisplacement(*d_new, *d_new, *a_new, dt);
498 workingState->setSolutionStatus(sStatus);
499 workingState->setOrder(this->getOrder());
500 workingState->computeNorms(currentState);
502 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,