263 const Teuchos::RCP<SolutionHistory<Scalar>>& solutionHistory) {
264#ifdef VERBOSE_DEBUG_OUTPUT
265 *out_ <<
"DEBUG: " << __PRETTY_FUNCTION__ <<
"\n";
267 this->checkInitialized();
271 TEMPUS_FUNC_TIME_MONITOR(
"Tempus::StepperNewmarkImplicitDForm::takeStep()");
273 TEUCHOS_TEST_FOR_EXCEPTION(solutionHistory->getNumStates() < 2,
275 "Error - StepperNewmarkImplicitDForm<Scalar>::takeStep(...)\n"
276 "Need at least two SolutionStates for NewmarkImplicitDForm.\n"
277 " Number of States = " << solutionHistory->getNumStates() <<
"\n"
278 "Try setting in \"Solution History\" \"Storage Type\" = \"Undo\"\n"
279 " or \"Storage Type\" = \"Static\" and \"Storage Limit\" = \"2\"\n");
281 auto thisStepper = Teuchos::rcpFromRef(*
this);
282 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,
285 RCP<SolutionState<Scalar>> workingState =solutionHistory->getWorkingState();
286 RCP<SolutionState<Scalar>> currentState =solutionHistory->getCurrentState();
288 Teuchos::RCP<WrapperModelEvaluatorSecondOrder<Scalar> > wrapperModel =
289 Teuchos::rcp_dynamic_cast<WrapperModelEvaluatorSecondOrder<Scalar> >(
290 this->wrapperModel_);
293 RCP<const Thyra::VectorBase<Scalar>> d_old = currentState->getX();
294 RCP<Thyra::VectorBase<Scalar>> v_old = currentState->getXDot();
295 RCP<Thyra::VectorBase<Scalar>> a_old = currentState->getXDotDot();
299 RCP<Thyra::VectorBase<Scalar>> d_new = workingState->getX();
300 RCP<Thyra::VectorBase<Scalar>> v_new = workingState->getXDot();
301 RCP<Thyra::VectorBase<Scalar>> a_new = workingState->getXDotDot();
304 const Scalar time = currentState->getTime();
305 const Scalar dt = workingState->getTimeStep();
307 Scalar t = time + dt;
311 Teuchos::Range1D range;
313 *out_ <<
"\n*** d_old ***\n";
314 RTOpPack::ConstSubVectorView<Scalar> dov;
315 d_old->acquireDetachedView(range, &dov);
316 auto doa = dov.values();
317 for (
auto i = 0; i < doa.size(); ++i) *out_ << doa[i] <<
" ";
318 *out_ <<
"\n*** d_old ***\n";
320 *out_ <<
"\n*** v_old ***\n";
321 RTOpPack::ConstSubVectorView<Scalar> vov;
322 v_old->acquireDetachedView(range, &vov);
323 auto voa = vov.values();
324 for (
auto i = 0; i < voa.size(); ++i) *out_ << voa[i] <<
" ";
325 *out_ <<
"\n*** v_old ***\n";
327 *out_ <<
"\n*** a_old ***\n";
328 RTOpPack::ConstSubVectorView<Scalar> aov;
329 a_old->acquireDetachedView(range, &aov);
330 auto aoa = aov.values();
331 for (
auto i = 0; i < aoa.size(); ++i) *out_ << aoa[i] <<
" ";
332 *out_ <<
"\n*** a_old ***\n";
336 RCP<Thyra::VectorBase<Scalar>> d_pred = Thyra::createMember(d_old->space());
337 RCP<Thyra::VectorBase<Scalar>> v_pred = Thyra::createMember(v_old->space());
340 predictDisplacement(*d_pred, *d_old, *v_old, *a_old, dt);
341 predictVelocity(*v_pred, *v_old, *a_old, dt);
344 *out_ <<
"\n*** d_pred ***\n";
345 RTOpPack::ConstSubVectorView<Scalar> dpv;
346 d_pred->acquireDetachedView(range, &dpv);
347 auto dpa = dpv.values();
348 for (
auto i = 0; i < dpa.size(); ++i) *out_ << dpa[i] <<
" ";
349 *out_ <<
"\n*** d_pred ***\n";
351 *out_ <<
"\n*** v_pred ***\n";
352 RTOpPack::ConstSubVectorView<Scalar> vpv;
353 v_pred->acquireDetachedView(range, &vpv);
354 auto vpa = vpv.values();
355 for (
auto i = 0; i < vpa.size(); ++i) *out_ << vpa[i] <<
" ";
356 *out_ <<
"\n*** v_pred ***\n";
360 wrapperModel->initializeNewmark(v_pred, d_pred, dt, t, beta_, gamma_);
363 RCP<Thyra::VectorBase<Scalar>> initial_guess = Thyra::createMember(d_pred->space());
364 if ((time == solutionHistory->minTime()) && (this->initialGuess_ != Teuchos::null)) {
367 bool is_compatible = (initial_guess->space())->isCompatible(*this->initialGuess_->space());
368 TEUCHOS_TEST_FOR_EXCEPTION(
369 is_compatible !=
true, std::logic_error,
370 "Error in Tempus::NemwarkImplicitDForm takeStep(): user-provided initial guess'!\n"
371 <<
"for Newton is not compatible with solution vector!\n");
372 Thyra::copy(*this->initialGuess_, initial_guess.ptr());
376 Thyra::copy(*d_pred, initial_guess.ptr());
379 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,
384 const Thyra::SolveStatus<Scalar> sStatus = (*(this->solver_)).solve(&*initial_guess);
386 workingState->setSolutionStatus(sStatus);
388 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,
393 Thyra::copy(*initial_guess, d_new.ptr());
396 correctAcceleration(*a_new, *d_pred, *d_new, dt);
397 correctVelocity(*v_new, *v_pred, *a_new, dt);
400 *out_ <<
"\n*** d_new ***\n";
401 RTOpPack::ConstSubVectorView<Scalar> dnv;
402 d_new->acquireDetachedView(range, &dnv);
403 auto dna = dnv.values();
404 for (
auto i = 0; i < dna.size(); ++i) *out_ << dna[i] <<
" ";
405 *out_ <<
"\n*** d_new ***\n";
407 *out_ <<
"\n*** v_new ***\n";
408 RTOpPack::ConstSubVectorView<Scalar> vnv;
409 v_new->acquireDetachedView(range, &vnv);
410 auto vna = vnv.values();
411 for (
auto i = 0; i < vna.size(); ++i) *out_ << vna[i] <<
" ";
412 *out_ <<
"\n*** v_new ***\n";
414 *out_ <<
"\n*** a_new ***\n";
415 RTOpPack::ConstSubVectorView<Scalar> anv;
416 a_new->acquireDetachedView(range, &anv);
417 auto ana = anv.values();
418 for (
auto i = 0; i < ana.size(); ++i) *out_ << ana[i] <<
" ";
419 *out_ <<
"\n*** a_new ***\n";
422 workingState->setOrder(this->getOrder());
423 workingState->computeNorms(currentState);
425 stepperNewmarkImpAppAction_->execute(solutionHistory, thisStepper,