Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 std::move(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 std::move(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
75  // redesign this for a generic second-order scheme to not have an
76  // if statement here...
77  // IKT, 3/14/17: this is effectively the same as the
78  // Piro::NewmarkDecorator::evalModel function.
79 
80  // The solution variable in NOX is the acceleration, a_{n+1}
81  appInArgs.set_x_dot_dot(inArgs.get_x());
82 
83  // Compute the velocity.
84  // v_n = velocity_{pred} + \gamma dt a_n
85  auto velocity = Thyra::createMember(inArgs.get_x()->space());
86  Thyra::V_StVpStV(velocity.ptr(),
87  Scalar(1.0), *v_pred_, delta_t_*gamma_, *inArgs.get_x());
88  appInArgs.set_x_dot(velocity);
89 
90  // Compute the displacement.
91  // d_n = displacement_{pred} + \beta dt^2 a_n
92  auto displacement = Thyra::createMember(inArgs.get_x()->space());
93  Thyra::V_StVpStV(displacement.ptr(),
94  Scalar(1.0), *d_pred_, beta_*delta_t_*delta_t_, *inArgs.get_x());
95  appInArgs.set_x(displacement);
96 
97  appInArgs.set_W_x_dot_dot_coeff(Scalar(1.0)); // da/da
98  appInArgs.set_alpha(gamma_*delta_t_); // dv/da
99  appInArgs.set_beta(beta_*delta_t_*delta_t_); // dd/da
100 
101  appInArgs.set_t(t_);
102  for (int i=0; i<appModel_->Np(); ++i) {
103  if (inArgs.get_p(i) != Teuchos::null)
104  appInArgs.set_p(i, inArgs.get_p(i));
105  }
106 
107  //Setup output condition
108  //Create and populate outArgs
109  appOutArgs.set_f(outArgs.get_f());
110  appOutArgs.set_W_op(outArgs.get_W_op());
111 
112  // build residual and jacobian
113  appModel_->evalModel(appInArgs,appOutArgs);
114  break;
115  }
116 
117  case NEWMARK_IMPLICIT_DFORM: {
118  // Setup initial condition
119  // Populate inArgs
120  RCP<Thyra::VectorBase<Scalar> const>
121  d = inArgs.get_x();
122 
123  RCP<Thyra::VectorBase<Scalar>>
124  v = Thyra::createMember(inArgs.get_x()->space());
125 
126  RCP<Thyra::VectorBase<Scalar>>
127  a = Thyra::createMember(inArgs.get_x()->space());
128 
129 #ifdef DEBUG_OUTPUT
130  Teuchos::Range1D range;
131 
132  *out_ << "\n*** d_bef ***\n";
133  RTOpPack::ConstSubVectorView<Scalar> dov;
134  d->acquireDetachedView(range, &dov);
135  auto doa = dov.values();
136  for (auto i = 0; i < doa.size(); ++i) *out_ << doa[i] << " ";
137  *out_ << "\n*** d_bef ***\n";
138 
139  *out_ << "\n*** v_bef ***\n";
140  RTOpPack::ConstSubVectorView<Scalar> vov;
141  v->acquireDetachedView(range, &vov);
142  auto voa = vov.values();
143  for (auto i = 0; i < voa.size(); ++i) *out_ << voa[i] << " ";
144  *out_ << "\n*** v_bef ***\n";
145 
146  *out_ << "\n*** a_bef ***\n";
147  RTOpPack::ConstSubVectorView<Scalar> aov;
148  a->acquireDetachedView(range, &aov);
149  auto aoa = aov.values();
150  for (auto i = 0; i < aoa.size(); ++i) *out_ << aoa[i] << " ";
151  *out_ << "\n*** a_bef ***\n";
152 #endif
153 
154  Scalar const
155  c = 1.0 / beta_ / delta_t_ / delta_t_;
156 
157  // compute acceleration
158  // a_{n+1} = (d_{n+1} - d_pred) / dt / dt / beta
159  Thyra::V_StVpStV(Teuchos::ptrFromRef(*a), c, *d, -c, *d_pred_);
160 
161  // compute velocity
162  // v_{n+1} = v_pred + \gamma dt a_{n+1}
163  Thyra::V_StVpStV(
164  Teuchos::ptrFromRef(*v), 1.0, *v_pred_, delta_t_ * gamma_, *a);
165 
166  appInArgs.set_x(d);
167  appInArgs.set_x_dot(v);
168  appInArgs.set_x_dot_dot(a);
169 
170  appInArgs.set_W_x_dot_dot_coeff(c); // da/dd
171  appInArgs.set_alpha(gamma_ / delta_t_ / beta_); // dv/dd
172  appInArgs.set_beta(1.0); // dd/dd
173 
174  appInArgs.set_t(t_);
175  for (int i = 0; i < appModel_->Np(); ++i) {
176  if (inArgs.get_p(i) != Teuchos::null)
177  appInArgs.set_p(i, inArgs.get_p(i));
178  }
179 
180  // Setup output condition
181  // Create and populate outArgs
182  appOutArgs.set_f(outArgs.get_f());
183  appOutArgs.set_W_op(outArgs.get_W_op());
184 
185  // build residual and jacobian
186  appModel_->evalModel(appInArgs, appOutArgs);
187 
188  // compute acceleration
189  // a_{n+1} = (d_{n+1} - d_pred) / dt / dt / beta
190  Thyra::V_StVpStV(Teuchos::ptrFromRef(*a), c, *d, -c, *d_pred_);
191 
192  // compute velocity
193  // v_{n+1} = v_pred + \gamma dt a_{n+1}
194  Thyra::V_StVpStV(
195  Teuchos::ptrFromRef(*v), 1.0, *v_pred_, delta_t_ * gamma_, *a);
196 
197  appInArgs.set_x(d);
198  appInArgs.set_x_dot(v);
199  appInArgs.set_x_dot_dot(a);
200 
201 #ifdef DEBUG_OUTPUT
202  *out_ << "\n*** d_aft ***\n";
203  RTOpPack::ConstSubVectorView<Scalar> dnv;
204  d->acquireDetachedView(range, &dnv);
205  auto dna = dnv.values();
206  for (auto i = 0; i < dna.size(); ++i) *out_ << dna[i] << " ";
207  *out_ << "\n*** d_aft ***\n";
208 
209  *out_ << "\n*** v_aft ***\n";
210  RTOpPack::ConstSubVectorView<Scalar> vnv;
211  v->acquireDetachedView(range, &vnv);
212  auto vna = vnv.values();
213  for (auto i = 0; i < vna.size(); ++i) *out_ << vna[i] << " ";
214  *out_ << "\n*** v_aft ***\n";
215 
216  *out_ << "\n*** a_aft ***\n";
217  RTOpPack::ConstSubVectorView<Scalar> anv;
218  a->acquireDetachedView(range, &anv);
219  auto ana = anv.values();
220  for (auto i = 0; i < ana.size(); ++i) *out_ << ana[i] << " ";
221  *out_ << "\n*** a_aft ***\n";
222 #endif
223  break;
224  }
225  }
226 }
227 
228 
229 } // namespace Tempus
230 
231 #endif // Tempus_WrapperModelEvaluatorSecondOrder_impl_hpp
Tempus::WrapperModelEvaluatorSecondOrder::createOutArgsImpl
Thyra::ModelEvaluatorBase::OutArgs< Scalar > createOutArgsImpl() const
Definition: Tempus_WrapperModelEvaluatorSecondOrder_impl.hpp:37
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::WrapperModelEvaluatorSecondOrder::createInArgs
Thyra::ModelEvaluatorBase::InArgs< Scalar > createInArgs() const
Definition: Tempus_WrapperModelEvaluatorSecondOrder_impl.hpp:18