VertexCFD  0.0-dev
CDR_Model_impl.hpp
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 VERTEXCFD_CDR_MODEL_IMPL_HPP
10 #define VERTEXCFD_CDR_MODEL_IMPL_HPP
11 
12 // Thyra support
13 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
14 #include "Thyra_DefaultSpmdVectorSpace.hpp"
15 #include "Thyra_DetachedMultiVectorView.hpp"
16 #include "Thyra_DetachedVectorView.hpp"
17 #include "Thyra_MultiVectorStdOps.hpp"
18 #include "Thyra_PreconditionerBase.hpp"
19 #include "Thyra_VectorStdOps.hpp"
20 
21 // Epetra support
22 #include "Epetra_Comm.h"
23 #include "Epetra_CrsGraph.h"
24 #include "Epetra_CrsMatrix.h"
25 #include "Epetra_Import.h"
26 #include "Epetra_Map.h"
27 #include "Epetra_Vector.h"
28 #include "Thyra_EpetraLinearOp.hpp"
29 #include "Thyra_EpetraThyraWrappers.hpp"
30 #include "Thyra_get_Epetra_Operator.hpp"
31 
32 namespace VertexCFD
33 {
34 namespace Test
35 {
36 // Constructor
37 
38 template<class Scalar>
39 CDR_Model<Scalar>::CDR_Model(const Teuchos::RCP<const Epetra_Comm>& comm,
40  const int num_global_elements,
41  const Scalar z_min,
42  const Scalar z_max,
43  const Scalar a,
44  const Scalar k)
45  : comm_(comm)
46  , num_global_elements_(num_global_elements)
47  , z_min_(z_min)
48  , z_max_(z_max)
49  , a_(a)
50  , k_(k)
51  , showGetInvalidArg_(false)
52 {
53  using Teuchos::RCP;
54  using Teuchos::rcp;
55  using ::Thyra::VectorBase;
56  typedef ::Thyra::ModelEvaluatorBase MEB;
57  typedef Teuchos::ScalarTraits<Scalar> ST;
58 
59  TEUCHOS_ASSERT(nonnull(comm_));
60 
61  const int num_nodes = num_global_elements_ + 1;
62 
63  // owned space
64  x_owned_map_ = rcp(new Epetra_Map(num_nodes, 0, *comm_));
65  x_space_ = ::Thyra::create_VectorSpace(x_owned_map_);
66 
67  // ghosted space
68  if (comm_->NumProc() == 1)
69  {
70  x_ghosted_map_ = x_owned_map_;
71  }
72  else
73  {
74  int OverlapNumMyElements;
75  int OverlapMinMyGID;
76  OverlapNumMyElements = x_owned_map_->NumMyElements() + 2;
77  if ((comm_->MyPID() == 0) || (comm_->MyPID() == (comm_->NumProc() - 1)))
78  OverlapNumMyElements--;
79 
80  if (comm_->MyPID() == 0)
81  OverlapMinMyGID = x_owned_map_->MinMyGID();
82  else
83  OverlapMinMyGID = x_owned_map_->MinMyGID() - 1;
84 
85  int* OverlapMyGlobalElements = new int[OverlapNumMyElements];
86 
87  for (int i = 0; i < OverlapNumMyElements; i++)
88  OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
89 
90  x_ghosted_map_ = Teuchos::rcp(new Epetra_Map(
91  -1, OverlapNumMyElements, OverlapMyGlobalElements, 0, *comm_));
92 
93  delete[] OverlapMyGlobalElements;
94  }
95 
96  importer_ = Teuchos::rcp(new Epetra_Import(*x_ghosted_map_, *x_owned_map_));
97 
98  // residual space
99  f_owned_map_ = x_owned_map_;
100  f_space_ = x_space_;
101 
102  x0_ = ::Thyra::createMember(x_space_);
103  V_S(x0_.ptr(), ST::zero());
104 
105  // set_p(Teuchos::tuple<Scalar>(p0, p1)());
106  // set_x0(Teuchos::tuple<Scalar>(x0, x1)());
107 
108  // Initialize the graph for W CrsMatrix object
109  W_graph_ = createGraph();
110 
111  // Create the nodal coordinates
112  {
113  node_coordinates_ = Teuchos::rcp(new Epetra_Vector(*x_owned_map_));
114  Scalar length = z_max_ - z_min_;
115  Scalar dx = length / ((double)num_global_elements_ - 1);
116  for (int i = 0; i < x_owned_map_->NumMyElements(); i++)
117  {
118  (*node_coordinates_)[i]
119  = z_min_ + dx * ((double)x_owned_map_->MinMyGID() + i);
120  }
121  }
122 
123  MEB::InArgsSetup<Scalar> inArgs;
124  inArgs.setModelEvalDescription(this->description());
125  inArgs.setSupports(MEB::IN_ARG_t);
126  inArgs.setSupports(MEB::IN_ARG_x);
127  inArgs.setSupports(MEB::IN_ARG_beta);
128  inArgs.setSupports(MEB::IN_ARG_x_dot);
129  inArgs.setSupports(MEB::IN_ARG_alpha);
130  prototypeInArgs_ = inArgs;
131 
132  MEB::OutArgsSetup<Scalar> outArgs;
133  outArgs.setModelEvalDescription(this->description());
134  outArgs.setSupports(MEB::OUT_ARG_f);
135  outArgs.setSupports(MEB::OUT_ARG_W);
136  outArgs.setSupports(MEB::OUT_ARG_W_op);
137  outArgs.setSupports(MEB::OUT_ARG_W_prec);
138  // outArgs.set_W_properties(DerivativeProperties(
139  // DERIV_LINEARITY_NONCONST
140  // ,DERIV_RANK_FULL
141  // ,true // supportsAdjoint
142  // ));
143  prototypeOutArgs_ = outArgs;
144 
145  // Setup nominal values
146  nominalValues_ = inArgs;
147  nominalValues_.set_x(x0_);
148  auto x_dot_init = Thyra::createMember(this->get_x_space());
149  Thyra::put_scalar(Scalar(0.0), x_dot_init.ptr());
150  nominalValues_.set_x_dot(x_dot_init);
151 }
152 
153 // Initializers/Accessors
154 
155 template<class Scalar>
156 Teuchos::RCP<Epetra_CrsGraph> CDR_Model<Scalar>::createGraph()
157 {
158  Teuchos::RCP<Epetra_CrsGraph> W_graph;
159 
160  // Create the shell for the
161  W_graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, *x_owned_map_, 5));
162 
163  // Declare required variables
164  int row;
165  int column;
166  const int OverlapNumMyElements = x_ghosted_map_->NumMyElements();
167 
168  // Loop Over # of Finite Elements on Processor
169  for (int ne = 0; ne < OverlapNumMyElements - 1; ne++)
170  {
171  // Loop over Nodes in Element
172  for (int i = 0; i < 2; i++)
173  {
174  row = x_ghosted_map_->GID(ne + i);
175 
176  // Loop over Trial Functions
177  for (int j = 0; j < 2; j++)
178  {
179  // If this row is owned by current processor, add the index
180  if (x_owned_map_->MyGID(row))
181  {
182  column = x_ghosted_map_->GID(ne + j);
183  W_graph->InsertGlobalIndices(row, 1, &column);
184  }
185  }
186  }
187  }
188  W_graph->FillComplete();
189  return W_graph;
190 }
191 
192 template<class Scalar>
193 void CDR_Model<Scalar>::set_x0(const Teuchos::ArrayView<const Scalar>& x0_in)
194 {
195 #ifdef TEUCHOS_DEBUG
196  TEUCHOS_ASSERT_EQUALITY(x_space_->dim(), x0_in.size());
197 #endif
198  Thyra::DetachedVectorView<Scalar> x0(x0_);
199  x0.sv().values()().assign(x0_in);
200 }
201 
202 template<class Scalar>
203 void CDR_Model<Scalar>::setShowGetInvalidArgs(bool showGetInvalidArg)
204 {
205  showGetInvalidArg_ = showGetInvalidArg;
206 }
207 
208 template<class Scalar>
209 void CDR_Model<Scalar>::set_W_factory(
210  const Teuchos::RCP<const ::Thyra::LinearOpWithSolveFactoryBase<Scalar>>&
211  W_factory)
212 {
213  W_factory_ = W_factory;
214 }
215 
216 // Public functions overridden from ModelEvaluator
217 
218 template<class Scalar>
219 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
220 CDR_Model<Scalar>::get_x_space() const
221 {
222  return x_space_;
223 }
224 
225 template<class Scalar>
226 Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar>>
227 CDR_Model<Scalar>::get_f_space() const
228 {
229  return f_space_;
230 }
231 
232 template<class Scalar>
233 Thyra::ModelEvaluatorBase::InArgs<Scalar>
234 CDR_Model<Scalar>::getNominalValues() const
235 {
236  return nominalValues_;
237 }
238 
239 template<class Scalar>
240 Teuchos::RCP<Thyra::LinearOpWithSolveBase<double>>
241 CDR_Model<Scalar>::create_W() const
242 {
243  const Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<double>> W_factory
244  = this->get_W_factory();
245 
246  TEUCHOS_TEST_FOR_EXCEPTION(is_null(W_factory),
247  std::runtime_error,
248  "W_factory in CDR_Model has a null W_factory!");
249 
250  const Teuchos::RCP<Thyra::LinearOpBase<double>> matrix
251  = this->create_W_op();
252 
253  Teuchos::RCP<Thyra::LinearOpWithSolveBase<double>> W
254  = Thyra::linearOpWithSolve<double>(*W_factory, matrix);
255 
256  return W;
257 }
258 
259 template<class Scalar>
260 Teuchos::RCP<Thyra::LinearOpBase<Scalar>> CDR_Model<Scalar>::create_W_op() const
261 {
262  const Teuchos::RCP<Epetra_CrsMatrix> W_epetra
263  = Teuchos::rcp(new Epetra_CrsMatrix(::Copy, *W_graph_));
264 
265  return Thyra::nonconstEpetraLinearOp(W_epetra);
266 }
267 
268 template<class Scalar>
269 Teuchos::RCP<::Thyra::PreconditionerBase<Scalar>>
270 CDR_Model<Scalar>::create_W_prec() const
271 {
272  const Teuchos::RCP<Epetra_CrsMatrix> W_epetra
273  = Teuchos::rcp(new Epetra_CrsMatrix(::Copy, *W_graph_));
274 
275  const Teuchos::RCP<Thyra::LinearOpBase<Scalar>> W_op
276  = Thyra::nonconstEpetraLinearOp(W_epetra);
277 
278  const Teuchos::RCP<Thyra::DefaultPreconditioner<Scalar>> prec
279  = Teuchos::rcp(new Thyra::DefaultPreconditioner<Scalar>);
280 
281  prec->initializeRight(W_op);
282 
283  return prec;
284 }
285 
286 template<class Scalar>
287 Teuchos::RCP<const Thyra::LinearOpWithSolveFactoryBase<Scalar>>
288 CDR_Model<Scalar>::get_W_factory() const
289 {
290  return W_factory_;
291 }
292 
293 template<class Scalar>
294 Thyra::ModelEvaluatorBase::InArgs<Scalar>
295 CDR_Model<Scalar>::createInArgs() const
296 {
297  return prototypeInArgs_;
298 }
299 
300 // Private functions overridden from ModelEvaluatorDefaultBase
301 
302 template<class Scalar>
303 Thyra::ModelEvaluatorBase::OutArgs<Scalar>
304 CDR_Model<Scalar>::createOutArgsImpl() const
305 {
306  return prototypeOutArgs_;
307 }
308 
309 template<class Scalar>
310 void CDR_Model<Scalar>::evalModelImpl(
311  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
312  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
313 {
314  using Teuchos::RCP;
315  using Teuchos::rcp;
316  using Teuchos::rcp_dynamic_cast;
317 
318  TEUCHOS_ASSERT(nonnull(inArgs.get_x()));
319  TEUCHOS_ASSERT(nonnull(inArgs.get_x_dot()));
320 
321  // const Thyra::ConstDetachedVectorView<Scalar> x(inArgs.get_x());
322 
323  const RCP<Thyra::VectorBase<Scalar>> f_out = outArgs.get_f();
324  const RCP<Thyra::LinearOpBase<Scalar>> W_out = outArgs.get_W_op();
325  const RCP<Thyra::PreconditionerBase<Scalar>> W_prec_out
326  = outArgs.get_W_prec();
327 
328  if (nonnull(f_out) || nonnull(W_out) || nonnull(W_prec_out))
329  {
330  // ****************
331  // Get the underlying epetra objects
332  // ****************
333 
334  RCP<Epetra_Vector> f;
335  if (nonnull(f_out))
336  {
337  f = Thyra::get_Epetra_Vector(*f_owned_map_, outArgs.get_f());
338  }
339 
340  RCP<Epetra_CrsMatrix> J;
341  if (nonnull(W_out))
342  {
343  const RCP<Epetra_Operator> W_epetra
344  = Thyra::get_Epetra_Operator(*W_out);
345  J = rcp_dynamic_cast<Epetra_CrsMatrix>(W_epetra);
346  TEUCHOS_ASSERT(nonnull(J));
347  }
348 
349  RCP<Epetra_CrsMatrix> M_inv;
350  if (nonnull(W_prec_out))
351  {
352  const RCP<Epetra_Operator> M_epetra = Thyra::get_Epetra_Operator(
353  *(W_prec_out->getNonconstRightPrecOp()));
354  M_inv = rcp_dynamic_cast<Epetra_CrsMatrix>(M_epetra);
355  TEUCHOS_ASSERT(nonnull(M_inv));
356  J_diagonal_ = Teuchos::rcp(new Epetra_Vector(*x_owned_map_));
357  }
358 
359  // ****************
360  // Create ghosted objects
361  // ****************
362 
363  // Set the boundary condition directly. Works for both x and xDot
364  // solves.
365  if (comm_->MyPID() == 0)
366  {
367  const RCP<Thyra::VectorBase<Scalar>> x
368  = Teuchos::rcp_const_cast<Thyra::VectorBase<Scalar>>(
369  inArgs.get_x());
370  (*Thyra::get_Epetra_Vector(*x_owned_map_, x))[0] = 1.0;
371  }
372 
373  if (is_null(u_ptr))
374  u_ptr = Teuchos::rcp(new Epetra_Vector(*x_ghosted_map_));
375 
376  u_ptr->Import(
377  *(Thyra::get_Epetra_Vector(*x_owned_map_, inArgs.get_x())),
378  *importer_,
379  Insert);
380 
381  if (is_null(u_dot_ptr))
382  u_dot_ptr = Teuchos::rcp(new Epetra_Vector(*x_ghosted_map_));
383 
384  u_dot_ptr->Import(
385  *(Thyra::get_Epetra_Vector(*x_owned_map_, inArgs.get_x_dot())),
386  *importer_,
387  Insert);
388 
389  if (is_null(x_ptr))
390  {
391  x_ptr = Teuchos::rcp(new Epetra_Vector(*x_ghosted_map_));
392  x_ptr->Import(*node_coordinates_, *importer_, Insert);
393  }
394 
395  Epetra_Vector& u = *u_ptr;
396  Epetra_Vector& u_dot = *u_dot_ptr;
397  Epetra_Vector& x = *x_ptr;
398 
399  int ierr = 0;
400  const int OverlapNumMyElements = x_ghosted_map_->NumMyElements();
401 
402  double xx[2];
403  double uu[2];
404  double uu_dot[2];
405  Basis basis;
406  const auto alpha = inArgs.get_alpha();
407  const auto beta = inArgs.get_beta();
408 
409  // Zero out the objects that will be filled
410  if (nonnull(f))
411  f->PutScalar(0.0);
412  if (nonnull(J))
413  J->PutScalar(0.0);
414  if (nonnull(M_inv))
415  M_inv->PutScalar(0.0);
416 
417  // Loop Over # of Finite Elements on Processor
418  for (int ne = 0; ne < OverlapNumMyElements - 1; ne++)
419  {
420  // Loop Over Gauss Points
421  for (int gp = 0; gp < 2; gp++)
422  {
423  // Get the solution and coordinates at the nodes
424  xx[0] = x[ne];
425  xx[1] = x[ne + 1];
426  uu[0] = u[ne];
427  uu[1] = u[ne + 1];
428  uu_dot[0] = u_dot[ne];
429  uu_dot[1] = u_dot[ne + 1];
430  // Calculate the basis function at the gauss point
431  basis.computeBasis(gp, xx, uu, uu_dot);
432 
433  // Loop over Nodes in Element
434  for (int i = 0; i < 2; i++)
435  {
436  const int row = x_ghosted_map_->GID(ne + i);
437  // printf("Proc=%d GlobalRow=%d LocalRow=%d Owned=%d\n",
438  // MyPID, row, ne+i,x_owned_map_.MyGID(row));
439  if (x_owned_map_->MyGID(row))
440  {
441  if (nonnull(f))
442  {
443  (*f)[x_owned_map_->LID(x_ghosted_map_->GID(ne + i))]
444  += +basis.wt * basis.dz
445  * (basis.uu_dot * basis.phi[i] // transient
446  + (a_ / basis.dz * basis.duu
447  * basis.phi[i] // convection
448  + 1.0 / (basis.dz * basis.dz))
449  * basis.duu
450  * basis.dphide[i] // diffusion
451  + k_ * basis.uu * basis.uu
452  * basis.phi[i]); // source
453  }
454  }
455  // Loop over Trial Functions
456  if (nonnull(J))
457  {
458  for (int j = 0; j < 2; j++)
459  {
460  if (x_owned_map_->MyGID(row))
461  {
462  const int column = x_ghosted_map_->GID(ne + j);
463  const double jac
464  = basis.wt * basis.dz
465  * (alpha * basis.phi[i]
466  * basis.phi[j] // transient
467  + beta
468  * (+a_ / basis.dz
469  * basis.dphide[j]
470  * basis.phi[i] // convection
471  + (1.0 / (basis.dz * basis.dz))
472  * basis.dphide[j]
473  * basis.dphide[i] // diffusion
474  + 2.0 * k_ * basis.uu
475  * basis.phi[j]
476  * basis.phi[i] // source
477  ));
478  ierr = J->SumIntoGlobalValues(
479  row, 1, &jac, &column);
480  }
481  }
482  }
483  if (nonnull(M_inv))
484  {
485  for (int j = 0; j < 2; j++)
486  {
487  if (x_owned_map_->MyGID(row))
488  {
489  const int column = x_ghosted_map_->GID(ne + j);
490  // The prec will be the diagonal of J. No need
491  // to assemble the other entries
492  if (row == column)
493  {
494  const double jac
495  = basis.wt * basis.dz
496  * (alpha * basis.phi[i]
497  * basis.phi[j] // transient
498  + beta
499  * (+a_ / basis.dz
500  * basis.dphide[j]
501  * basis.phi[i] // convection
502  + (1.0
503  / (basis.dz * basis.dz))
504  * basis.dphide[j]
505  * basis.dphide[i] // diffusion
506  + 2.0 * k_ * basis.uu
507  * basis.phi[j]
508  * basis.phi[i] // source
509  ));
510  ierr = M_inv->SumIntoGlobalValues(
511  row, 1, &jac, &column);
512  }
513  }
514  }
515  }
516  }
517  }
518  }
519 
520  // Insert Boundary Conditions and modify Jacobian and function (F)
521  // U(0)=1
522  if (comm_->MyPID() == 0)
523  {
524  if (nonnull(f))
525  (*f)[0] = 0.0; // Setting BC above and zero residual here works
526  // for x and xDot solves.
527  //(*f)[0]= u[0] - 1.0; // BC equation works for x solves.
528  if (nonnull(J))
529  {
530  int column = 0;
531  double jac = 1.0;
532  ierr = J->ReplaceGlobalValues(0, 1, &jac, &column);
533  column = 1;
534  jac = 0.0;
535  ierr = J->ReplaceGlobalValues(0, 1, &jac, &column);
536  }
537  if (nonnull(M_inv))
538  {
539  int column = 0;
540  double jac = 1.0;
541  ierr = M_inv->ReplaceGlobalValues(0, 1, &jac, &column);
542  column = 1;
543  jac = 0.0;
544  ierr = M_inv->ReplaceGlobalValues(0, 1, &jac, &column);
545  }
546  }
547 
548  if (nonnull(J))
549  J->FillComplete();
550 
551  if (nonnull(M_inv))
552  {
553  // invert the Jacobian diagonal for the preconditioner
554  M_inv->ExtractDiagonalCopy(*J_diagonal_);
555 
556  for (int i = 0; i < J_diagonal_->MyLength(); ++i)
557  (*J_diagonal_)[i] = 1.0 / ((*J_diagonal_)[i]);
558 
559  M_inv->PutScalar(0.0);
560  M_inv->ReplaceDiagonalValues(*J_diagonal_);
561  M_inv->FillComplete();
562  }
563 
564  TEUCHOS_ASSERT(ierr > -1);
565  }
566 }
567 
568 //====================================================================
569 // Basis vector
570 
571 // Constructor
572 Basis::Basis()
573  : uu(0.0)
574  , zz(0.0)
575  , duu(0.0)
576  , eta(0.0)
577  , wt(0.0)
578  , dz(0.0)
579  , uu_dot(0.0)
580  , duu_dot(0.0)
581 {
582  phi = new double[2];
583  dphide = new double[2];
584 }
585 
586 // Destructor
587 Basis::~Basis()
588 {
589  delete[] phi;
590  delete[] dphide;
591 }
592 
593 // Calculates a linear 1D basis
594 void Basis::computeBasis(int gp, double* z, double* u, double* u_dot)
595 {
596  const int N = 2;
597  if (gp == 0)
598  {
599  eta = -1.0 / sqrt(3.0);
600  wt = 1.0;
601  }
602  if (gp == 1)
603  {
604  eta = 1.0 / sqrt(3.0);
605  wt = 1.0;
606  }
607 
608  // Calculate basis function and derivatives at nodel pts
609  phi[0] = (1.0 - eta) / 2.0;
610  phi[1] = (1.0 + eta) / 2.0;
611  dphide[0] = -0.5;
612  dphide[1] = 0.5;
613 
614  // Caculate basis function and derivative at GP.
615  dz = 0.5 * (z[1] - z[0]);
616  zz = 0.0;
617  uu = 0.0;
618  duu = 0.0;
619  uu_dot = 0.0;
620  duu_dot = 0.0;
621  for (int i = 0; i < N; i++)
622  {
623  zz += z[i] * phi[i];
624  uu += u[i] * phi[i];
625  duu += u[i] * dphide[i];
626  if (u_dot)
627  {
628  uu_dot += u_dot[i] * phi[i];
629  duu_dot += u_dot[i] * dphide[i];
630  }
631  }
632 
633  return;
634 }
635 
636 } // namespace Test
637 } // namespace VertexCFD
638 
639 #endif
VertexCFD
Definition: tstMethodManufacturedSolutionBC.cpp:23
Generated by doxygen 1.8.20