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