38CDR_Model(
const Teuchos::RCP<const Epetra_Comm>& comm,
39 const int num_global_elements,
45 num_global_elements_(num_global_elements),
50 showGetInvalidArg_(false)
54 using ::Thyra::VectorBase;
55 typedef ::Thyra::ModelEvaluatorBase MEB;
56 typedef Teuchos::ScalarTraits<Scalar> ST;
58 TEUCHOS_ASSERT(nonnull(
comm_));
67 if (
comm_->NumProc() == 1) {
71 int OverlapNumMyElements;
73 OverlapNumMyElements =
x_owned_map_->NumMyElements() + 2;
74 if ( (
comm_->MyPID() == 0) || (
comm_->MyPID() == (
comm_->NumProc() - 1)) )
75 OverlapNumMyElements --;
77 if (
comm_->MyPID() == 0)
82 int* OverlapMyGlobalElements =
new int[OverlapNumMyElements];
84 for (
int i = 0; i < OverlapNumMyElements; i ++)
85 OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
88 Teuchos::rcp(
new Epetra_Map(-1,
90 OverlapMyGlobalElements,
94 delete [] OverlapMyGlobalElements;
105 V_S(
x0_.ptr(), ST::zero());
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 );
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);
150 auto x_dot_init = Thyra::createMember(this->
get_x_space());
151 Thyra::put_scalar(Scalar(0.0),x_dot_init.ptr());
318 const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
319 const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
324 using Teuchos::rcp_dynamic_cast;
326 TEUCHOS_ASSERT(nonnull(inArgs.get_x()));
327 TEUCHOS_ASSERT(nonnull(inArgs.get_x_dot()));
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();
336 if ( nonnull(f_out) || nonnull(W_out) || nonnull(W_prec_out) ) {
342 RCP<Epetra_Vector> f;
343 if (nonnull(f_out)) {
344 f = Thyra::get_Epetra_Vector(*f_owned_map_,outArgs.get_f());
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));
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_));
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;
373 u_ptr = Teuchos::rcp(
new Epetra_Vector(*x_ghosted_map_));
375 u_ptr->Import(*(Thyra::get_Epetra_Vector(*x_owned_map_,inArgs.get_x())), *importer_, Insert);
377 if (is_null(u_dot_ptr))
378 u_dot_ptr = Teuchos::rcp(
new Epetra_Vector(*x_ghosted_map_));
380 u_dot_ptr->Import(*(Thyra::get_Epetra_Vector(*x_owned_map_,inArgs.get_x_dot())), *importer_, Insert);
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);
387 Epetra_Vector& u = *u_ptr;
388 Epetra_Vector& u_dot = *u_dot_ptr;
389 Epetra_Vector& x = *x_ptr;
392 int OverlapNumMyElements = x_ghosted_map_->NumMyElements();
398 const auto alpha = inArgs.get_alpha();
399 const auto beta = inArgs.get_beta();
407 M_inv->PutScalar(0.0);
412 for (
int ne=0; ne < OverlapNumMyElements-1; ne++) {
415 for(
int gp=0; gp < 2; gp++) {
422 uu_dot[1]=u_dot[ne+1];
424 basis.computeBasis(gp, xx, uu, uu_dot);
427 for (
int i=0; i< 2; i++) {
428 int row=x_ghosted_map_->GID(ne+i);
431 if (x_owned_map_->MyGID(row)) {
433 (*f)[x_owned_map_->LID(x_ghosted_map_->GID(ne+i))]+=
435 *( basis.uu_dot*basis.phi[i]
436 +(a_/basis.dz*basis.duu*basis.phi[i]
437 +1.0/(basis.dz*basis.dz))*basis.duu*basis.dphide[i]
438 +k_*basis.uu*basis.uu*basis.phi[i] );
443 for(
int j=0;j < 2; j++) {
444 if (x_owned_map_->MyGID(row)) {
445 int column=x_ghosted_map_->GID(ne+j);
448 alpha * basis.phi[i] * basis.phi[j]
450 +a_/basis.dz*basis.dphide[j]*basis.phi[i]
451 +(1.0/(basis.dz*basis.dz))*basis.dphide[j]*basis.dphide[i]
452 +2.0*k_*basis.uu*basis.phi[j]*basis.phi[i]
455 ierr=J->SumIntoGlobalValues(row, 1, &jac, &column);
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);
467 alpha * basis.phi[i] * basis.phi[j]
469 +a_/basis.dz*basis.dphide[j]*basis.phi[i]
470 +(1.0/(basis.dz*basis.dz))*basis.dphide[j]*basis.dphide[i]
471 +2.0*k_*basis.uu*basis.phi[j]*basis.phi[i]
474 ierr = M_inv->SumIntoGlobalValues(row, 1, &jac, &column);
485 if (comm_->MyPID() == 0) {
492 ierr = J->ReplaceGlobalValues(0, 1, &jac, &column);
495 ierr = J->ReplaceGlobalValues(0, 1, &jac, &column);
497 if (nonnull(M_inv)) {
500 ierr = M_inv->ReplaceGlobalValues(0, 1, &jac, &column);
503 ierr = M_inv->ReplaceGlobalValues(0, 1, &jac, &column);
510 if (nonnull(M_inv)) {
512 M_inv->ExtractDiagonalCopy(*J_diagonal_);
514 for (
int i=0; i < J_diagonal_->MyLength(); ++i)
515 (*J_diagonal_)[i] = 1.0 / ((*J_diagonal_)[i]);
517 M_inv->PutScalar(0.0);
518 M_inv->ReplaceDiagonalValues(*J_diagonal_);
519 M_inv->FillComplete();
522 TEUCHOS_ASSERT(ierr > -1);