Rythmos - Transient Integration for Differential Equations  Version of the Day
 All Classes Functions Variables Typedefs Pages
EpetraExt_DiagonalTransientModel.cpp
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 
42 
43 #include "EpetraExt_DiagonalTransientModel.hpp"
44 #include "Epetra_Comm.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_CrsGraph.h"
47 #include "Epetra_CrsMatrix.h"
48 #include "Epetra_LocalMap.h"
49 #include "Teuchos_ParameterList.hpp"
50 #include "Teuchos_StandardParameterEntryValidators.hpp"
51 #include "Teuchos_Assert.hpp"
52 #include "Teuchos_as.hpp"
53 
54 
55 namespace {
56 
57 
58 using Teuchos::RCP;
59 
60 
61 const std::string Implicit_name = "Implicit";
62 const bool Implicit_default = true;
63 
64 const std::string Gamma_min_name = "Gamma_min";
65 const double Gamma_min_default = -0.9;
66 
67 const std::string Gamma_max_name = "Gamma_max";
68 const double Gamma_max_default = -0.01;
69 
70 const std::string Coeff_s_name = "Coeff_s";
71 const std::string Coeff_s_default = "{ 0.0 }";
72 
73 const std::string Gamma_fit_name = "Gamma_fit";
74 const std::string Gamma_fit_default = "Linear"; // Will be validated at runtime!
75 
76 const std::string NumElements_name = "NumElements";
77 const int NumElements_default = 1;
78 
79 const std::string x0_name = "x0";
80 const double x0_default = 10.0;
81 
82 const std::string ExactSolutionAsResponse_name = "Exact Solution as Response";
83 const bool ExactSolutionAsResponse_default = false;
84 
85 
86 inline
87 double evalR( const double& t, const double& gamma, const double& s )
88 {
89  return (exp(gamma*t)*sin(s*t));
90 }
91 
92 
93 inline
94 double d_evalR_d_s( const double& t, const double& gamma, const double& s )
95 {
96  return (exp(gamma*t)*cos(s*t)*t);
97 }
98 
99 
100 inline
101 double f_func(
102  const double& x, const double& t, const double& gamma, const double& s
103  )
104 {
105  return ( gamma*x + evalR(t,gamma,s) );
106 }
107 
108 
109 inline
110 double x_exact(
111  const double& x0, const double& t, const double& gamma, const double& s
112  )
113 {
114  if ( s == 0.0 )
115  return ( x0 * exp(gamma*t) );
116  return ( exp(gamma*t) * (x0 + (1.0/s) * ( 1.0 - cos(s*t) ) ) );
117  // Note that the limit of (1.0/s) * ( 1.0 - cos(s*t) ) as s goes to zero is
118  // zero. This limit is neeed to avoid the 0/0 that would occur if floating
119  // point was used to evaluate this term. This means that cos(t*s) goes to
120  // one at the same rate as s goes to zero giving 1-1=0..
121 }
122 
123 
124 inline
125 double dxds_exact(
126  const double& t, const double& gamma, const double& s
127  )
128 {
129  if ( s == 0.0 )
130  return 0.0;
131  return ( -exp(gamma*t)/(s*s) * ( 1.0 - sin(s*t)*(s*t) - cos(s*t) ) );
132 }
133 
134 
135 class UnsetParameterVector {
136 public:
137  ~UnsetParameterVector()
138  {
139  if (!is_null(vec_))
140  *vec_ = Teuchos::null;
141  }
142  UnsetParameterVector(
143  const RCP<RCP<const Epetra_Vector> > &vec
144  )
145  {
146  setVector(vec);
147  }
148  void setVector( const RCP<RCP<const Epetra_Vector> > &vec )
149  {
150  vec_ = vec;
151  }
152 private:
153  RCP<RCP<const Epetra_Vector> > vec_;
154 };
155 
156 
157 } // namespace
158 
159 
160 namespace EpetraExt {
161 
162 
163 // Constructors
164 
165 
167  Teuchos::RCP<Epetra_Comm> const& epetra_comm
168  )
169  : epetra_comm_(epetra_comm.assert_not_null()),
170  implicit_(Implicit_default),
171  numElements_(NumElements_default),
172  gamma_min_(Gamma_min_default),
173  gamma_max_(Gamma_max_default),
174  coeff_s_(Teuchos::fromStringToArray<double>(Coeff_s_default)),
175  gamma_fit_(GAMMA_FIT_LINEAR), // Must be the same as Gamma_fit_default!
176  x0_(x0_default),
177  exactSolutionAsResponse_(ExactSolutionAsResponse_default),
178  isIntialized_(false)
179 {
180  initialize();
181 }
182 
183 
184 Teuchos::RCP<const Epetra_Vector>
186 {
187  return gamma_;
188 }
189 
190 
191 Teuchos::RCP<const Epetra_Vector>
193  const double t, const Epetra_Vector *coeff_s_p
194  ) const
195 {
196  set_coeff_s_p(Teuchos::rcp(coeff_s_p,false));
197  UnsetParameterVector unsetParameterVector(Teuchos::rcp(&coeff_s_p_,false));
198  Teuchos::RCP<Epetra_Vector>
199  x_star_ptr = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
200  Epetra_Vector& x_star = *x_star_ptr;
201  Epetra_Vector& gamma = *gamma_;
202  int myN = x_star.MyLength();
203  for ( int i=0 ; i<myN ; ++i ) {
204  x_star[i] = x_exact( x0_, t, gamma[i], coeff_s(i) );
205  }
206  return(x_star_ptr);
207 }
208 
209 
210 Teuchos::RCP<const Epetra_MultiVector>
212  const double t, const Epetra_Vector *coeff_s_p
213  ) const
214 {
215  set_coeff_s_p(Teuchos::rcp(coeff_s_p,false));
216  UnsetParameterVector unsetParameterVector(Teuchos::rcp(&coeff_s_p_,false));
217  Teuchos::RCP<Epetra_MultiVector>
218  dxds_star_ptr = Teuchos::rcp(new Epetra_MultiVector(*epetra_map_,np_,false));
219  Epetra_MultiVector& dxds_star = *dxds_star_ptr;
220  dxds_star.PutScalar(0.0);
221  Epetra_Vector& gamma = *gamma_;
222  int myN = dxds_star.MyLength();
223  for ( int i=0 ; i<myN ; ++i ) {
224  const int coeff_s_idx_i = this->coeff_s_idx(i);
225  (*dxds_star(coeff_s_idx_i))[i] = dxds_exact( t, gamma[i], coeff_s(i) );
226  // Note: Above, at least the column access will be validated in debug mode
227  // but the row index i will not ever be. Perhaps we can augment Epetra to
228  // fix this?
229  }
230  return (dxds_star_ptr);
231 }
232 
233 
234 // Overridden from ParameterListAcceptor
235 
236 
238  Teuchos::RCP<Teuchos::ParameterList> const& paramList
239  )
240 {
241  using Teuchos::get; using Teuchos::getIntegralValue;
242  using Teuchos::getArrayFromStringParameter;
243  TEUCHOS_TEST_FOR_EXCEPT( is_null(paramList) );
244  paramList->validateParametersAndSetDefaults(*this->getValidParameters());
245  isIntialized_ = false;
246  paramList_ = paramList;
247  implicit_ = get<bool>(*paramList_,Implicit_name);
248  numElements_ = get<int>(*paramList_,NumElements_name);
249  gamma_min_ = get<double>(*paramList_,Gamma_min_name);
250  gamma_max_ = get<double>(*paramList_,Gamma_max_name);
251  coeff_s_ = getArrayFromStringParameter<double>(*paramList_,Coeff_s_name);
252  gamma_fit_ = getIntegralValue<EGammaFit>(*paramList_,Gamma_fit_name);
253  x0_ = get<double>(*paramList_,x0_name);
254  exactSolutionAsResponse_ = get<bool>(*paramList_,ExactSolutionAsResponse_name);
255  initialize();
256 }
257 
258 
259 Teuchos::RCP<Teuchos::ParameterList>
261 {
262  return paramList_;
263 }
264 
265 
266 Teuchos::RCP<Teuchos::ParameterList>
268 {
269  Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
270  paramList_ = Teuchos::null;
271  return _paramList;
272 }
273 
274 
275 Teuchos::RCP<const Teuchos::ParameterList>
277 {
278  return paramList_;
279 }
280 
281 
282 Teuchos::RCP<const Teuchos::ParameterList>
284 {
285  using Teuchos::RCP; using Teuchos::ParameterList;
286  using Teuchos::tuple;
287  using Teuchos::setIntParameter; using Teuchos::setDoubleParameter;
288  using Teuchos::setStringToIntegralParameter;
289  static RCP<const ParameterList> validPL;
290  if (is_null(validPL)) {
291  RCP<ParameterList> pl = Teuchos::parameterList();
292  pl->set(Implicit_name,true);
293  setDoubleParameter(
294  Gamma_min_name, Gamma_min_default, "",
295  &*pl
296  );
297  setDoubleParameter(
298  Gamma_max_name, Gamma_max_default, "",
299  &*pl
300  );
301  setStringToIntegralParameter<EGammaFit>(
302  Gamma_fit_name, Gamma_fit_default, "",
303  tuple<std::string>("Linear","Random"),
304  tuple<EGammaFit>(GAMMA_FIT_LINEAR,GAMMA_FIT_RANDOM),
305  &*pl
306  );
307  setIntParameter(
308  NumElements_name, NumElements_default, "",
309  &*pl
310  );
311  setDoubleParameter(
312  x0_name, x0_default, "",
313  &*pl
314  );
315  pl->set( Coeff_s_name, Coeff_s_default );
316  pl->set( ExactSolutionAsResponse_name, ExactSolutionAsResponse_default );
317  validPL = pl;
318  }
319  return validPL;
320 }
321 
322 
323 // Overridden from EpetraExt::ModelEvaluator
324 
325 
326 Teuchos::RCP<const Epetra_Map>
328 {
329  return epetra_map_;
330 }
331 
332 
333 Teuchos::RCP<const Epetra_Map>
335 {
336  return epetra_map_;
337 }
338 
339 
340 Teuchos::RCP<const Epetra_Map>
342 {
343 #ifdef TEUCHOS_DEBUG
344  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
345 #endif
346  return map_p_[l];
347 }
348 
349 
350 Teuchos::RCP<const Teuchos::Array<std::string> >
352 {
353 #ifdef TEUCHOS_DEBUG
354  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
355 #endif
356  return names_p_[l];
357 }
358 
359 
360 Teuchos::RCP<const Epetra_Map>
362 {
363 #ifdef TEUCHOS_DEBUG
364  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
365 #endif
366  return map_g_[j];
367 }
368 
369 
370 Teuchos::RCP<const Epetra_Vector>
372 {
373  return x_init_;
374 }
375 
376 
377 Teuchos::RCP<const Epetra_Vector>
379 {
380  return x_dot_init_;
381 }
382 
383 
384 Teuchos::RCP<const Epetra_Vector>
386 {
387 #ifdef TEUCHOS_DEBUG
388  TEUCHOS_ASSERT( l == 0 );
389 #endif
390  return p_init_[l];
391 }
392 
393 
394 Teuchos::RCP<Epetra_Operator>
396 {
397  if(implicit_)
398  return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
399  return Teuchos::null;
400 }
401 
402 
403 EpetraExt::ModelEvaluator::InArgs
405 {
406  InArgsSetup inArgs;
407  inArgs.set_Np(Np_);
408  inArgs.setSupports(IN_ARG_t,true);
409  inArgs.setSupports(IN_ARG_x,true);
410  if(implicit_) {
411  inArgs.setSupports(IN_ARG_x_dot,true);
412  inArgs.setSupports(IN_ARG_alpha,true);
413  inArgs.setSupports(IN_ARG_beta,true);
414  }
415  return inArgs;
416 }
417 
418 
419 EpetraExt::ModelEvaluator::OutArgs
421 {
422  OutArgsSetup outArgs;
423  outArgs.set_Np_Ng(Np_,Ng_);
424  outArgs.setSupports(OUT_ARG_f,true);
425  if(implicit_) {
426  outArgs.setSupports(OUT_ARG_W,true);
427  outArgs.set_W_properties(
428  DerivativeProperties(
429  DERIV_LINEARITY_NONCONST
430  ,DERIV_RANK_FULL
431  ,true // supportsAdjoint
432  )
433  );
434  }
435  outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
436  outArgs.set_DfDp_properties(
437  0,DerivativeProperties(
438  DERIV_LINEARITY_NONCONST,
439  DERIV_RANK_DEFICIENT,
440  true // supportsAdjoint
441  )
442  );
443  if (exactSolutionAsResponse_) {
444  outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_MV_BY_COL);
445  outArgs.set_DgDp_properties(
446  0,0,DerivativeProperties(
447  DERIV_LINEARITY_NONCONST,
448  DERIV_RANK_DEFICIENT,
449  true // supportsAdjoint
450  )
451  );
452  }
453  return outArgs;
454 }
455 
456 
458  const InArgs& inArgs, const OutArgs& outArgs
459  ) const
460 {
461 
462  using Teuchos::rcp;
463  using Teuchos::RCP;
464  using Teuchos::null;
465  using Teuchos::dyn_cast;
466 
467  const Epetra_Vector &x = *(inArgs.get_x());
468  const double t = inArgs.get_t();
469  if (Np_) set_coeff_s_p(inArgs.get_p(0)); // Sets for coeff_s(...) function!
470  UnsetParameterVector unsetParameterVector(rcp(&coeff_s_p_,false));
471  // Note: Above, the destructor for unsetParameterVector will ensure that the
472  // RCP to the parameter vector will be unset no matter if an exception is
473  // thrown or not.
474 
475  Epetra_Operator *W_out = ( implicit_ ? outArgs.get_W().get() : 0 );
476  Epetra_Vector *f_out = outArgs.get_f().get();
477  Epetra_MultiVector *DfDp_out = 0;
478  if (Np_) DfDp_out = get_DfDp_mv(0,outArgs).get();
479  Epetra_Vector *g_out = 0;
480  Epetra_MultiVector *DgDp_out = 0;
481  if (exactSolutionAsResponse_) {
482  g_out = outArgs.get_g(0).get();
483  DgDp_out = get_DgDp_mv(0,0,outArgs,DERIV_MV_BY_COL).get();
484  }
485 
486  const Epetra_Vector &gamma = *gamma_;
487 
488  int localNumElements = x.MyLength();
489 
490  if (f_out) {
491  Epetra_Vector &f = *f_out;
492  if (implicit_) {
493  const Epetra_Vector *x_dot_in = inArgs.get_x_dot().get();
494  if (x_dot_in) {
495  const Epetra_Vector &x_dot = *x_dot_in;
496  for ( int i=0 ; i<localNumElements ; ++i )
497  f[i] = x_dot[i] - f_func(x[i],t,gamma[i],coeff_s(i));
498  }
499  else {
500  for ( int i=0 ; i<localNumElements ; ++i )
501  f[i] = - f_func(x[i],t,gamma[i],coeff_s(i));
502  }
503  }
504  else {
505  for ( int i=0 ; i<localNumElements ; ++i )
506  f[i] = f_func(x[i],t,gamma[i],coeff_s(i));
507  }
508  }
509 
510  if ( W_out ) {
511  // If we get here then we are in implicit mode!
512  const double alpha = inArgs.get_alpha();
513  const double beta = inArgs.get_beta();
514  Epetra_CrsMatrix &crsW = dyn_cast<Epetra_CrsMatrix>(*W_out);
515  double values[1];
516  int indices[1];
517  const int offset
518  = epetra_comm_->MyPID()*localNumElements + epetra_map_->IndexBase();
519  for( int i = 0; i < localNumElements; ++i ) {
520  values[0] = alpha - beta*gamma[i];
521  indices[0] = i + offset; // global column
522  crsW.ReplaceGlobalValues(
523  i + offset // GlobalRow
524  ,1 // NumEntries
525  ,values // Values
526  ,indices // Indices
527  );
528  }
529  }
530 
531  if (DfDp_out) {
532  Epetra_MultiVector &DfDp = *DfDp_out;
533  DfDp.PutScalar(0.0);
534  if (implicit_) {
535  for( int i = 0; i < localNumElements; ++i ) {
536  DfDp[coeff_s_idx(i)][i]
537  = - d_evalR_d_s(t,gamma[i],coeff_s(i));
538  }
539  }
540  else {
541  for( int i = 0; i < localNumElements; ++i ) {
542  DfDp[coeff_s_idx(i)][i]
543  = + d_evalR_d_s(t,gamma[i],coeff_s(i));
544  }
545  }
546  }
547 
548  if (g_out) {
549  *g_out = *getExactSolution(t,&*coeff_s_p_);
550  // Note: Above will wipe out coeff_s_p_ as a side effect!
551  }
552 
553  if (DgDp_out) {
554  *DgDp_out = *getExactSensSolution(t,&*coeff_s_p_);
555  // Note: Above will wipe out coeff_s_p_ as a side effect!
556  }
557 
558  // Warning: From here on out coeff_s_p_ is unset!
559 
560 }
561 
562 
563 // private
564 
565 
566 void DiagonalTransientModel::initialize()
567 {
568 
569  using Teuchos::rcp;
570  using Teuchos::as;
571 
572  //
573  // Setup map
574  //
575 
576  const int numProcs = epetra_comm_->NumProc();
577  const int procRank = epetra_comm_->MyPID();
578  epetra_map_ = rcp( new Epetra_Map( numElements_ * numProcs, 0, *epetra_comm_ ) );
579 
580  //
581  // Create gamma
582  //
583 
584  gamma_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_));
585  Epetra_Vector &gamma = *gamma_;
586  switch(gamma_fit_) {
587  case GAMMA_FIT_LINEAR: {
588  const int N = gamma.GlobalLength();
589  const double slope = (gamma_max_ - gamma_min_)/(N-1);
590  const int MyLength = gamma.MyLength();
591  if (1==MyLength) {
592  gamma[0] = gamma_min_;
593  }
594  else {
595  for ( int i=0 ; i<MyLength ; ++i )
596  {
597  gamma[i] = slope*(procRank*MyLength+i)+gamma_min_;
598  }
599  }
600  break;
601  }
602  case GAMMA_FIT_RANDOM: {
603  unsigned int seed = Teuchos::ScalarTraits<int>::random()+10*procRank;
604  seed *= seed;
605  gamma.SetSeed(seed);
606  gamma.Random(); // fill with random numbers in (-1,1)
607  // Scale random numbers to (gamma_min_,gamma_max_)
608  const double slope = (gamma_min_ - gamma_max_)/2.0;
609  const double tmp = (gamma_max_ + gamma_min_)/2.0;
610  int MyLength = gamma.MyLength();
611  for (int i=0 ; i<MyLength ; ++i)
612  {
613  gamma[i] = slope*gamma[i] + tmp;
614  }
615  break;
616  }
617  default:
618  TEUCHOS_TEST_FOR_EXCEPT("Should never get here!");
619  }
620 
621  //
622  // Setup for parameters
623  //
624 
625  Np_ = 1;
626  np_ = coeff_s_.size();
627  map_p_.clear();
628  map_p_.push_back(
629  rcp( new Epetra_LocalMap(np_,0,*epetra_comm_) )
630  );
631  names_p_.clear();
632  {
633  Teuchos::RCP<Teuchos::Array<std::string> >
634  names_p_0 = Teuchos::rcp(new Teuchos::Array<std::string>());
635  for ( int i = 0; i < np_; ++i )
636  names_p_0->push_back("coeff_s("+Teuchos::toString(i)+")");
637  names_p_.push_back(names_p_0);
638  }
639 
640  coeff_s_idx_.clear();
641  const int num_func_per_coeff_s = numElements_ / np_;
642  const int num_func_per_coeff_s_rem = numElements_ % np_;
643  for ( int coeff_s_idx_i = 0; coeff_s_idx_i < np_; ++coeff_s_idx_i ) {
644  const int num_func_per_coeff_s_idx_i
645  = num_func_per_coeff_s
646  + ( coeff_s_idx_i < num_func_per_coeff_s_rem ? 1 : 0 );
647  for (
648  int coeff_s_idx_i_j = 0;
649  coeff_s_idx_i_j < num_func_per_coeff_s_idx_i;
650  ++coeff_s_idx_i_j
651  )
652  {
653  coeff_s_idx_.push_back(coeff_s_idx_i);
654  }
655  }
656 #ifdef TEUCHOS_DEBUG
657  TEUCHOS_TEST_FOR_EXCEPT(
658  ( as<int>(coeff_s_idx_.size()) != numElements_ ) && "Internal programming error!" );
659 #endif
660 
661  //
662  // Setup exact solution as response function
663  //
664 
665  if (exactSolutionAsResponse_) {
666  Ng_ = 1;
667  map_g_.clear();
668  map_g_.push_back(
669  rcp( new Epetra_LocalMap(1,0,*epetra_comm_) )
670  );
671  }
672  else {
673  Ng_ = 0;
674  }
675 
676  //
677  // Setup graph for W
678  //
679 
680  if(implicit_) {
681 
682  W_graph_ = rcp(
683  new Epetra_CrsGraph(
684  ::Copy,*epetra_map_,
685  1 // elements per row
686  )
687  );
688 
689  int indices[1];
690  const int offset = procRank*numElements_ + epetra_map_->IndexBase();
691 
692  for( int i = 0; i < numElements_; ++i ) {
693  indices[0] = i + offset; // global column
694  W_graph_->InsertGlobalIndices(
695  i + offset // GlobalRow
696  ,1 // NumEntries
697  ,indices // Indices
698  );
699  }
700 
701  W_graph_->FillComplete();
702 
703  }
704 
705  //
706  // Setup initial guess
707  //
708 
709  // Set x_init
710  x_init_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
711  x_init_->PutScalar(x0_);
712 
713  // Set x_dot_init to provide for a consistent inital guess for implicit mode
714  // such that f(x_dot,x) = 0!
715  if (implicit_) {
716  x_dot_init_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
717  InArgs inArgs = this->createInArgs();
718  inArgs.set_x(x_init_);
719  inArgs.set_t(0.0);
720  OutArgs outArgs = this->createOutArgs();
721  outArgs.set_f(x_dot_init_);
722  this->evalModel(inArgs,outArgs);
723  x_dot_init_->Scale(-1.0);
724  }
725 
726  // Set p_init
727  p_init_.clear();
728  p_init_.push_back(
729  rcp( new Epetra_Vector( ::Copy, *map_p_[0], &coeff_s_[0] ) )
730  );
731 
732 }
733 
734 
735 void DiagonalTransientModel::set_coeff_s_p(
736  const Teuchos::RCP<const Epetra_Vector> &coeff_s_p
737  ) const
738 {
739  if (!is_null(coeff_s_p))
740  coeff_s_p_ = coeff_s_p;
741  else
742  unset_coeff_s_p();
743 }
744 
745 
746 void DiagonalTransientModel::unset_coeff_s_p() const
747 {
748  using Teuchos::as;
749 #ifdef TEUCHOS_DEBUG
750  TEUCHOS_ASSERT(
751  as<int>(get_p_map(0)->NumGlobalElements()) == as<int>(coeff_s_.size()) );
752 #endif
753  coeff_s_p_ = Teuchos::rcp(
754  new Epetra_Vector(
755  ::View,
756  *get_p_map(0),
757  const_cast<double*>(&coeff_s_[0])
758  )
759  );
760  // Note: The above const cast is okay since the coeff_s_p_ RCP is to a const
761  // Epetr_Vector!
762 }
763 
764 
765 
766 } // namespace EpetraExt
767 
768 
769 // Nonmembers
770 
771 
772 Teuchos::RCP<EpetraExt::DiagonalTransientModel>
773 EpetraExt::diagonalTransientModel(
774  Teuchos::RCP<Epetra_Comm> const& epetra_comm,
775  Teuchos::RCP<Teuchos::ParameterList> const& paramList
776  )
777 {
778  Teuchos::RCP<DiagonalTransientModel> diagonalTransientModel =
779  Teuchos::rcp(new DiagonalTransientModel(epetra_comm));
780  if (!is_null(paramList))
781  diagonalTransientModel->setParameterList(paramList);
782  return diagonalTransientModel;
783 }
Teuchos::RCP< const Epetra_Vector > getExactSolution(const double t, const Epetra_Vector *coeff_s_p=0) const
Return the exact solution as a function of time.
Teuchos::RCP< const Epetra_Vector > get_gamma() const
Return the model vector gamma,.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
DiagonalTransientModel(Teuchos::RCP< Epetra_Comm > const &epetra_comm)
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Teuchos::RCP< const Epetra_Map > get_f_map() const
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< DiagonalTransientModel > diagonalTransientModel(Teuchos::RCP< Epetra_Comm > const &epetra_comm, Teuchos::RCP< Teuchos::ParameterList > const &paramList=Teuchos::null)
Nonmember constructor.
Teuchos::RCP< const Epetra_MultiVector > getExactSensSolution(const double t, const Epetra_Vector *coeff_s_p=0) const
Return the exact sensitivity of x as a function of time.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Teuchos::RCP< Epetra_Operator > create_W() const
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
Teuchos::RCP< const Epetra_Vector > get_x_dot_init() const
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< const Epetra_Map > get_g_map(int j) const