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