Tempus  Version of the Day
Time Integration
Tempus_WrapperModelEvaluatorSecondOrder_impl.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef Tempus_WrapperModelEvaluatorSecondOrder_impl_hpp
10 #define Tempus_WrapperModelEvaluatorSecondOrder_impl_hpp
11 
12 namespace Tempus {
13 
14 
15 template <typename Scalar>
16 Thyra::ModelEvaluatorBase::InArgs<Scalar>
17 WrapperModelEvaluatorSecondOrder<Scalar>::
18 createInArgs() const
19 {
20 #ifdef VERBOSE_DEBUG_OUTPUT
21  *out_ << "DEBUG: " << __PRETTY_FUNCTION__ << "\n";
22 #endif
23  typedef Thyra::ModelEvaluatorBase MEB;
24 
25  MEB::InArgsSetup<Scalar> inArgs;
26  inArgs.setModelEvalDescription(this->description());
27  inArgs.set_Np(appModel_->Np());
28  inArgs.setSupports(MEB::IN_ARG_x);
29 
30  return inArgs;
31 }
32 
33 
34 template <typename Scalar>
35 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
36 WrapperModelEvaluatorSecondOrder<Scalar>::
37 createOutArgsImpl() const
38 {
39 #ifdef VERBOSE_DEBUG_OUTPUT
40  *out_ << "DEBUG: " << __PRETTY_FUNCTION__ << "\n";
41 #endif
42  typedef Thyra::ModelEvaluatorBase MEB;
43 
44  MEB::OutArgsSetup<Scalar> outArgs;
45  outArgs.setModelEvalDescription(this->description());
46  outArgs.set_Np_Ng(appModel_->Np(),0);
47  outArgs.setSupports(MEB::OUT_ARG_f);
48  outArgs.setSupports(MEB::OUT_ARG_W_op);
49 
50  return outArgs;
51 }
52 
53 
54 template <typename Scalar>
55 void
56 WrapperModelEvaluatorSecondOrder<Scalar>::
57 evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
58  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs) const
59 {
60 #ifdef VERBOSE_DEBUG_OUTPUT
61  *out_ << "DEBUG: " << __PRETTY_FUNCTION__ << "\n";
62 #endif
63  typedef Thyra::ModelEvaluatorBase MEB;
64  using Teuchos::RCP;
65 
66  //Setup initial condition
67  //Create and populate inArgs
68  MEB::InArgs<Scalar> appInArgs = appModel_->createInArgs();
69  MEB::OutArgs<Scalar> appOutArgs = appModel_->createOutArgs();
70 
71  switch (schemeType_)
72  {
73  case NEWMARK_IMPLICIT_AFORM: {
74  //Specific for the Newmark Implicit a-Form stepper. May want to redesign this for a generic
75  //second order scheme to not have an if statement here...
76  //IKT, 3/14/17: this is effectively the same as the Piro::NewmarkDecorator::evalModel function.
77  //the solution variable in NOX is the acceleration, a_{n+1}
78  appInArgs.set_x_dot_dot(inArgs.get_x());
79  RCP<Thyra::VectorBase<Scalar> > velocity = Thyra::createMember(inArgs.get_x()->space());
80  //compute the velocity, v_{n+1}(a_{n+1}) = velocity_{pred} + \gamma dt a_{n+1}
81  Thyra::V_StVpStV(Teuchos::ptrFromRef(*velocity), 1.0, *v_pred_, delta_t_*gamma_, *inArgs.get_x());
82  appInArgs.set_x_dot(velocity);
83  RCP<Thyra::VectorBase<Scalar> > displacement = Thyra::createMember(inArgs.get_x()->space());
84  //compute the displacement, d_{n+1}(a_{n+1}) = displacement_{pred} + \beta dt^2 a_{n+1}
85  Thyra::V_StVpStV(Teuchos::ptrFromRef(*displacement), 1.0, *d_pred_, beta_*delta_t_*delta_t_, *inArgs.get_x());
86  appInArgs.set_x(displacement);
87  appInArgs.set_W_x_dot_dot_coeff(1.0); // da/da
88  appInArgs.set_alpha(gamma_*delta_t_); // dv/da
89  appInArgs.set_beta(beta_*delta_t_*delta_t_); // dd/da
90 
91  appInArgs.set_t(t_);
92  for (int i=0; i<appModel_->Np(); ++i) {
93  if (inArgs.get_p(i) != Teuchos::null)
94  appInArgs.set_p(i, inArgs.get_p(i));
95  }
96 
97  //Setup output condition
98  //Create and populate outArgs
99  appOutArgs.set_f(outArgs.get_f());
100  appOutArgs.set_W_op(outArgs.get_W_op());
101 
102  // build residual and jacobian
103  appModel_->evalModel(appInArgs,appOutArgs);
104  break;
105  }
106 
107  case NEWMARK_IMPLICIT_DFORM: {
108  // Setup initial condition
109  // Populate inArgs
110  RCP<Thyra::VectorBase<Scalar> const>
111  d = inArgs.get_x();
112 
113  RCP<Thyra::VectorBase<Scalar>>
114  v = Thyra::createMember(inArgs.get_x()->space());
115 
116  RCP<Thyra::VectorBase<Scalar>>
117  a = Thyra::createMember(inArgs.get_x()->space());
118 
119 #ifdef DEBUG_OUTPUT
120  Teuchos::Range1D range;
121 
122  *out_ << "\n*** d_bef ***\n";
123  RTOpPack::ConstSubVectorView<Scalar> dov;
124  d->acquireDetachedView(range, &dov);
125  auto doa = dov.values();
126  for (auto i = 0; i < doa.size(); ++i) *out_ << doa[i] << " ";
127  *out_ << "\n*** d_bef ***\n";
128 
129  *out_ << "\n*** v_bef ***\n";
130  RTOpPack::ConstSubVectorView<Scalar> vov;
131  v->acquireDetachedView(range, &vov);
132  auto voa = vov.values();
133  for (auto i = 0; i < voa.size(); ++i) *out_ << voa[i] << " ";
134  *out_ << "\n*** v_bef ***\n";
135 
136  *out_ << "\n*** a_bef ***\n";
137  RTOpPack::ConstSubVectorView<Scalar> aov;
138  a->acquireDetachedView(range, &aov);
139  auto aoa = aov.values();
140  for (auto i = 0; i < aoa.size(); ++i) *out_ << aoa[i] << " ";
141  *out_ << "\n*** a_bef ***\n";
142 #endif
143 
144  Scalar const
145  c = 1.0 / beta_ / delta_t_ / delta_t_;
146 
147  // compute acceleration
148  // a_{n+1} = (d_{n+1} - d_pred) / dt / dt / beta
149  Thyra::V_StVpStV(Teuchos::ptrFromRef(*a), c, *d, -c, *d_pred_);
150 
151  // compute velocity
152  // v_{n+1} = v_pred + \gamma dt a_{n+1}
153  Thyra::V_StVpStV(
154  Teuchos::ptrFromRef(*v), 1.0, *v_pred_, delta_t_ * gamma_, *a);
155 
156  appInArgs.set_x(d);
157  appInArgs.set_x_dot(v);
158  appInArgs.set_x_dot_dot(a);
159 
160  appInArgs.set_W_x_dot_dot_coeff(c); // da/dd
161  appInArgs.set_alpha(gamma_ / delta_t_ / beta_); // dv/dd
162  appInArgs.set_beta(1.0); // dd/dd
163 
164  appInArgs.set_t(t_);
165  for (int i = 0; i < appModel_->Np(); ++i) {
166  if (inArgs.get_p(i) != Teuchos::null)
167  appInArgs.set_p(i, inArgs.get_p(i));
168  }
169 
170  // Setup output condition
171  // Create and populate outArgs
172  appOutArgs.set_f(outArgs.get_f());
173  appOutArgs.set_W_op(outArgs.get_W_op());
174 
175  // build residual and jacobian
176  appModel_->evalModel(appInArgs, appOutArgs);
177 
178  // compute acceleration
179  // a_{n+1} = (d_{n+1} - d_pred) / dt / dt / beta
180  Thyra::V_StVpStV(Teuchos::ptrFromRef(*a), c, *d, -c, *d_pred_);
181 
182  // compute velocity
183  // v_{n+1} = v_pred + \gamma dt a_{n+1}
184  Thyra::V_StVpStV(
185  Teuchos::ptrFromRef(*v), 1.0, *v_pred_, delta_t_ * gamma_, *a);
186 
187  appInArgs.set_x(d);
188  appInArgs.set_x_dot(v);
189  appInArgs.set_x_dot_dot(a);
190 
191 #ifdef DEBUG_OUTPUT
192  *out_ << "\n*** d_aft ***\n";
193  RTOpPack::ConstSubVectorView<Scalar> dnv;
194  d->acquireDetachedView(range, &dnv);
195  auto dna = dnv.values();
196  for (auto i = 0; i < dna.size(); ++i) *out_ << dna[i] << " ";
197  *out_ << "\n*** d_aft ***\n";
198 
199  *out_ << "\n*** v_aft ***\n";
200  RTOpPack::ConstSubVectorView<Scalar> vnv;
201  v->acquireDetachedView(range, &vnv);
202  auto vna = vnv.values();
203  for (auto i = 0; i < vna.size(); ++i) *out_ << vna[i] << " ";
204  *out_ << "\n*** v_aft ***\n";
205 
206  *out_ << "\n*** a_aft ***\n";
207  RTOpPack::ConstSubVectorView<Scalar> anv;
208  a->acquireDetachedView(range, &anv);
209  auto ana = anv.values();
210  for (auto i = 0; i < ana.size(); ++i) *out_ << ana[i] << " ";
211  *out_ << "\n*** a_aft ***\n";
212 #endif
213  break;
214  }
215  }
216 }
217 
218 
219 } // namespace Tempus
220 
221 #endif // Tempus_WrapperModelEvaluatorSecondOrder_impl_hpp
Tempus::WrapperModelEvaluatorSecondOrder::evalModelImpl
void evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs< Scalar > &inArgs, const Thyra::ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const
Definition: Tempus_WrapperModelEvaluatorSecondOrder_impl.hpp:57
Tempus
Definition: Tempus_AdjointAuxSensitivityModelEvaluator_decl.hpp:20
Tempus::WrapperModelEvaluatorSecondOrder::createOutArgsImpl
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Definition: Tempus_WrapperModelEvaluatorSecondOrder_impl.hpp:37
Tempus::WrapperModelEvaluatorSecondOrder::createInArgs
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Definition: Tempus_WrapperModelEvaluatorSecondOrder_impl.hpp:18