EpetraExt_DiagonalTransientModel.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 
00030 #include "EpetraExt_DiagonalTransientModel.hpp"
00031 #include "Epetra_Comm.h"
00032 #include "Epetra_Map.h"
00033 #include "Epetra_CrsGraph.h"
00034 #include "Epetra_CrsMatrix.h"
00035 #include "Epetra_LocalMap.h"
00036 #include "Teuchos_ParameterList.hpp"
00037 #include "Teuchos_StandardParameterEntryValidators.hpp"
00038 #include "Teuchos_Assert.hpp"
00039 #include "Teuchos_as.hpp"
00040 
00041 
00042 namespace {
00043 
00044 
00045 using Teuchos::RCP;
00046 
00047 
00048 const std::string Implicit_name = "Implicit";
00049 const bool Implicit_default = true;
00050 
00051 const std::string Gamma_min_name = "Gamma_min";
00052 const double Gamma_min_default = -0.9;
00053 
00054 const std::string Gamma_max_name = "Gamma_max";
00055 const double Gamma_max_default = -0.01;
00056 
00057 const std::string Coeff_s_name = "Coeff_s";
00058 const string Coeff_s_default = "{ 0.0 }";
00059 
00060 const std::string Gamma_fit_name = "Gamma_fit";
00061 const std::string Gamma_fit_default = "Linear"; // Will be validated at runtime!
00062 
00063 const std::string NumElements_name = "NumElements";
00064 const int NumElements_default = 1;
00065 
00066 const std::string x0_name = "x0";
00067 const double x0_default = 10.0;
00068 
00069 const std::string ExactSolutionAsResponse_name = "Exact Solution as Response";
00070 const bool ExactSolutionAsResponse_default = false;
00071 
00072 
00073 inline
00074 double evalR( const double& t, const double& gamma, const double& s )
00075 {
00076   return (exp(gamma*t)*sin(s*t));
00077 }
00078 
00079 
00080 inline
00081 double d_evalR_d_s( const double& t, const double& gamma, const double& s )
00082 {
00083   return (exp(gamma*t)*cos(s*t)*t);
00084 }
00085 
00086 
00087 inline
00088 double f_func(
00089   const double& x, const double& t, const double& gamma, const double& s
00090   )
00091 {
00092   return ( gamma*x + evalR(t,gamma,s) );
00093 }
00094 
00095 
00096 inline
00097 double x_exact(
00098   const double& x0, const double& t, const double& gamma, const double& s
00099   )
00100 {
00101   if ( s == 0.0 )
00102     return ( x0 * exp(gamma*t) );
00103   return ( exp(gamma*t) * (x0 + (1.0/s) * ( 1.0 - cos(s*t) ) ) );
00104   // Note that the limit of (1.0/s) * ( 1.0 - cos(s*t) ) as s goes to zero is
00105   // zero.  This limit is neeed to avoid the 0/0 that would occur if floating
00106   // point was used to evaluate this term.  This means that cos(t*s) goes to
00107   // one at the same rate as s goes to zero giving 1-1=0..
00108 }
00109 
00110 
00111 inline
00112 double dxds_exact(
00113   const double& t, const double& gamma, const double& s
00114   )
00115 {
00116   if ( s == 0.0 )
00117     return 0.0;
00118   return ( -exp(gamma*t)/(s*s) * ( 1.0 - sin(s*t)*(s*t) - cos(s*t) ) );
00119 }
00120 
00121 
00122 class UnsetParameterVector {
00123 public:
00124   ~UnsetParameterVector()
00125     {
00126       if (!is_null(vec_))
00127         *vec_ = Teuchos::null;
00128     }
00129   UnsetParameterVector(
00130     const RCP<RCP<const Epetra_Vector> > &vec
00131     )
00132     {
00133       setVector(vec);
00134     }
00135   void setVector( const RCP<RCP<const Epetra_Vector> > &vec )
00136     {
00137       vec_ = vec;
00138     }
00139 private:
00140   RCP<RCP<const Epetra_Vector> > vec_;
00141 };
00142 
00143 
00144 } // namespace
00145 
00146 
00147 namespace EpetraExt {
00148 
00149 
00150 // Constructors
00151 
00152 
00153 DiagonalTransientModel::DiagonalTransientModel(
00154   Teuchos::RCP<Epetra_Comm> const& epetra_comm
00155   )
00156   : epetra_comm_(epetra_comm.assert_not_null()),
00157     implicit_(Implicit_default),
00158     numElements_(NumElements_default),
00159     gamma_min_(Gamma_min_default),
00160     gamma_max_(Gamma_max_default),
00161     coeff_s_(Teuchos::fromStringToArray<double>(Coeff_s_default)),
00162     gamma_fit_(GAMMA_FIT_LINEAR), // Must be the same as Gamma_fit_default!
00163     x0_(x0_default),
00164     exactSolutionAsResponse_(ExactSolutionAsResponse_default),
00165     isIntialized_(false)
00166 {
00167   initialize();
00168 }
00169 
00170 
00171 Teuchos::RCP<const Epetra_Vector>
00172 DiagonalTransientModel::get_gamma() const
00173 {
00174   return gamma_;
00175 }
00176 
00177 
00178 Teuchos::RCP<const Epetra_Vector>
00179 DiagonalTransientModel::getExactSolution(
00180   const double t, const Epetra_Vector *coeff_s_p
00181   ) const
00182 {
00183   set_coeff_s_p(Teuchos::rcp(coeff_s_p,false));
00184   UnsetParameterVector unsetParameterVector(Teuchos::rcp(&coeff_s_p_,false));
00185   Teuchos::RCP<Epetra_Vector>
00186     x_star_ptr = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
00187   Epetra_Vector& x_star = *x_star_ptr;
00188   Epetra_Vector& gamma = *gamma_;
00189   int myN = x_star.MyLength();
00190   for ( int i=0 ; i<myN ; ++i ) {
00191     x_star[i] = x_exact( x0_, t, gamma[i], coeff_s(i) );
00192   }
00193   return(x_star_ptr);
00194 }
00195 
00196 
00197 Teuchos::RCP<const Epetra_MultiVector>
00198 DiagonalTransientModel::getExactSensSolution(
00199   const double t, const Epetra_Vector *coeff_s_p
00200   ) const
00201 {
00202   set_coeff_s_p(Teuchos::rcp(coeff_s_p,false));
00203   UnsetParameterVector unsetParameterVector(Teuchos::rcp(&coeff_s_p_,false));
00204   Teuchos::RCP<Epetra_MultiVector>
00205     dxds_star_ptr = Teuchos::rcp(new Epetra_MultiVector(*epetra_map_,np_,false));
00206   Epetra_MultiVector& dxds_star = *dxds_star_ptr;
00207   dxds_star.PutScalar(0.0);
00208   Epetra_Vector& gamma = *gamma_;
00209   int myN = dxds_star.MyLength();
00210   for ( int i=0 ; i<myN ; ++i ) {
00211     const int coeff_s_idx = this->coeff_s_idx(i);
00212     (*dxds_star(coeff_s_idx))[i] = dxds_exact( t, gamma[i], coeff_s(i) );
00213     // Note: Above, at least the column access will be validated in debug mode
00214     // but the row index i will not ever be.  Perhaps we can augment Epetra to
00215     // fix this?
00216   }
00217   return (dxds_star_ptr);
00218 }
00219 
00220 
00221 // Overridden from ParameterListAcceptor
00222 
00223 
00224 void DiagonalTransientModel::setParameterList(
00225   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00226   )
00227 {
00228   using Teuchos::get; using Teuchos::getIntegralValue;
00229   using Teuchos::getArrayFromStringParameter;
00230   TEST_FOR_EXCEPT( is_null(paramList) );
00231   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00232   isIntialized_ = false;
00233   paramList_ = paramList;
00234   implicit_ = get<bool>(*paramList_,Implicit_name);
00235   numElements_ = get<int>(*paramList_,NumElements_name);
00236   gamma_min_ = get<double>(*paramList_,Gamma_min_name);
00237   gamma_max_ = get<double>(*paramList_,Gamma_max_name);
00238   coeff_s_ = getArrayFromStringParameter<double>(*paramList_,Coeff_s_name);
00239   gamma_fit_ = getIntegralValue<EGammaFit>(*paramList_,Gamma_fit_name);
00240   x0_ = get<double>(*paramList_,x0_name);
00241   exactSolutionAsResponse_ = get<bool>(*paramList_,ExactSolutionAsResponse_name);
00242   initialize();
00243 }
00244 
00245 
00246 Teuchos::RCP<Teuchos::ParameterList> 
00247 DiagonalTransientModel::getNonconstParameterList()
00248 {
00249   return paramList_;
00250 }
00251 
00252 
00253 Teuchos::RCP<Teuchos::ParameterList>
00254 DiagonalTransientModel::unsetParameterList()
00255 {
00256   Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
00257   paramList_ = Teuchos::null;
00258   return _paramList;
00259 }
00260 
00261 
00262 Teuchos::RCP<const Teuchos::ParameterList>
00263 DiagonalTransientModel::getParameterList() const
00264 {
00265   return paramList_;
00266 }
00267 
00268 
00269 Teuchos::RCP<const Teuchos::ParameterList>
00270 DiagonalTransientModel::getValidParameters() const
00271 {
00272   using Teuchos::RCP; using Teuchos::ParameterList;
00273   using Teuchos::tuple;
00274   using Teuchos::setIntParameter; using Teuchos::setDoubleParameter;
00275   using Teuchos::setStringToIntegralParameter;
00276   static RCP<const ParameterList> validPL;
00277   if (is_null(validPL)) {
00278     RCP<ParameterList> pl = Teuchos::parameterList();
00279     pl->set(Implicit_name,true);
00280     setDoubleParameter(
00281       Gamma_min_name, Gamma_min_default, "",
00282       &*pl
00283       );
00284     setDoubleParameter(
00285       Gamma_max_name, Gamma_max_default, "",
00286       &*pl
00287       );
00288     setStringToIntegralParameter<EGammaFit>(
00289       Gamma_fit_name, Gamma_fit_default, "",
00290       tuple<std::string>("Linear","Random"),
00291       tuple<EGammaFit>(GAMMA_FIT_LINEAR,GAMMA_FIT_RANDOM),
00292       &*pl
00293       );
00294     setIntParameter(
00295       NumElements_name, NumElements_default, "",
00296       &*pl
00297       );
00298     setDoubleParameter(
00299       x0_name, x0_default, "",
00300       &*pl
00301       );
00302     pl->set( Coeff_s_name, Coeff_s_default );
00303     pl->set( ExactSolutionAsResponse_name, ExactSolutionAsResponse_default );
00304     validPL = pl;
00305   }
00306   return validPL;
00307 }
00308 
00309 
00310 // Overridden from EpetraExt::ModelEvaluator
00311 
00312 
00313 Teuchos::RCP<const Epetra_Map>
00314 DiagonalTransientModel::get_x_map() const
00315 {
00316   return epetra_map_;
00317 }
00318 
00319 
00320 Teuchos::RCP<const Epetra_Map>
00321 DiagonalTransientModel::get_f_map() const
00322 {
00323   return epetra_map_;
00324 }
00325 
00326 
00327 Teuchos::RCP<const Epetra_Map>
00328 DiagonalTransientModel::get_p_map(int l) const
00329 {
00330 #ifdef TEUCHOS_DEBUG
00331   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
00332 #endif
00333   return map_p_[l];
00334 }
00335 
00336 
00337 Teuchos::RCP<const Teuchos::Array<std::string> >
00338 DiagonalTransientModel::get_p_names(int l) const
00339 {
00340 #ifdef TEUCHOS_DEBUG
00341   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
00342 #endif
00343   return names_p_[l];
00344 }
00345 
00346 
00347 Teuchos::RCP<const Epetra_Map>
00348 DiagonalTransientModel::get_g_map(int j) const
00349 {
00350 #ifdef TEUCHOS_DEBUG
00351   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
00352 #endif
00353   return map_g_[j];
00354 }
00355 
00356 
00357 Teuchos::RCP<const Epetra_Vector>
00358 DiagonalTransientModel::get_x_init() const
00359 {
00360   return x_init_;
00361 }
00362 
00363 
00364 Teuchos::RCP<const Epetra_Vector>
00365 DiagonalTransientModel::get_x_dot_init() const
00366 {
00367   return x_dot_init_;
00368 }
00369 
00370 
00371 Teuchos::RCP<const Epetra_Vector>
00372 DiagonalTransientModel::get_p_init(int l) const
00373 {
00374 #ifdef TEUCHOS_DEBUG
00375   TEUCHOS_ASSERT( l == 0 );
00376 #endif
00377   return p_init_[l];
00378 }
00379 
00380 
00381 Teuchos::RCP<Epetra_Operator>
00382 DiagonalTransientModel::create_W() const
00383 {
00384   if(implicit_)
00385     return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00386   return Teuchos::null;
00387 }
00388 
00389 
00390 EpetraExt::ModelEvaluator::InArgs
00391 DiagonalTransientModel::createInArgs() const
00392 {
00393   InArgsSetup inArgs;
00394   inArgs.set_Np(Np_);
00395   inArgs.setSupports(IN_ARG_t,true);
00396   inArgs.setSupports(IN_ARG_x,true);
00397   if(implicit_) {
00398     inArgs.setSupports(IN_ARG_x_dot,true);
00399     inArgs.setSupports(IN_ARG_alpha,true);
00400     inArgs.setSupports(IN_ARG_beta,true);
00401   }
00402   return inArgs;
00403 }
00404 
00405 
00406 EpetraExt::ModelEvaluator::OutArgs
00407 DiagonalTransientModel::createOutArgs() const
00408 {
00409   OutArgsSetup outArgs;
00410   outArgs.set_Np_Ng(Np_,Ng_);
00411   outArgs.setSupports(OUT_ARG_f,true);
00412   if(implicit_) {
00413     outArgs.setSupports(OUT_ARG_W,true);
00414     outArgs.set_W_properties(
00415       DerivativeProperties(
00416         DERIV_LINEARITY_NONCONST
00417         ,DERIV_RANK_FULL
00418         ,true // supportsAdjoint
00419         )
00420       );
00421   }
00422   outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
00423   outArgs.set_DfDp_properties(
00424     0,DerivativeProperties(
00425       DERIV_LINEARITY_NONCONST,
00426       DERIV_RANK_DEFICIENT,
00427       true // supportsAdjoint
00428       )
00429     );
00430   if (exactSolutionAsResponse_) {
00431     outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_MV_BY_COL);
00432     outArgs.set_DgDp_properties(
00433       0,0,DerivativeProperties(
00434         DERIV_LINEARITY_NONCONST,
00435         DERIV_RANK_DEFICIENT,
00436         true // supportsAdjoint
00437         )
00438       );
00439   }
00440   return outArgs;
00441 }
00442 
00443 
00444 void DiagonalTransientModel::evalModel(
00445   const InArgs& inArgs, const OutArgs& outArgs
00446   ) const
00447 {
00448 
00449   using Teuchos::rcp;
00450   using Teuchos::RCP;
00451   using Teuchos::null;
00452   using Teuchos::dyn_cast;
00453 
00454   const Epetra_Vector &x = *(inArgs.get_x());
00455   const double t = inArgs.get_t();
00456   if (Np_) set_coeff_s_p(inArgs.get_p(0)); // Sets for coeff_s(...) function!
00457   UnsetParameterVector unsetParameterVector(rcp(&coeff_s_p_,false));
00458   // Note: Above, the destructor for unsetParameterVector will ensure that the
00459   // RCP to the parameter vector will be unset no matter if an exception is
00460   // thrown or not.
00461 
00462   Epetra_Operator *W_out = ( implicit_ ? outArgs.get_W().get() : 0 );
00463   Epetra_Vector *f_out = outArgs.get_f().get();
00464   Epetra_MultiVector *DfDp_out = 0;
00465   if (Np_) DfDp_out = get_DfDp_mv(0,outArgs).get();
00466   Epetra_Vector *g_out = 0;
00467   Epetra_MultiVector *DgDp_out = 0;
00468   if (exactSolutionAsResponse_) {
00469     g_out = outArgs.get_g(0).get();
00470     DgDp_out = get_DgDp_mv(0,0,outArgs,DERIV_MV_BY_COL).get();
00471   }
00472 
00473   const Epetra_Vector &gamma = *gamma_;
00474 
00475   int localNumElements = x.MyLength();
00476 
00477   if (f_out) {
00478     Epetra_Vector &f = *f_out;
00479     if (implicit_) {
00480       const Epetra_Vector *x_dot_in = inArgs.get_x_dot().get();
00481       if (x_dot_in) {
00482         const Epetra_Vector &x_dot = *x_dot_in;
00483         for ( int i=0 ; i<localNumElements ; ++i )
00484           f[i] = x_dot[i] - f_func(x[i],t,gamma[i],coeff_s(i));
00485       }
00486       else {
00487         for ( int i=0 ; i<localNumElements ; ++i )
00488           f[i] = - f_func(x[i],t,gamma[i],coeff_s(i));
00489       }
00490     }
00491     else {
00492       for ( int i=0 ; i<localNumElements ; ++i )
00493         f[i] = f_func(x[i],t,gamma[i],coeff_s(i));
00494     }
00495   }
00496 
00497   if ( W_out ) {
00498     // If we get here then we are in implicit mode!
00499     const double alpha = inArgs.get_alpha();
00500     const double beta = inArgs.get_beta();
00501     Epetra_CrsMatrix &crsW = dyn_cast<Epetra_CrsMatrix>(*W_out);
00502     double values[1];
00503     int indices[1];
00504     const int offset
00505       = epetra_comm_->MyPID()*localNumElements + epetra_map_->IndexBase(); 
00506     for( int i = 0; i < localNumElements; ++i ) {
00507       values[0] = alpha - beta*gamma[i];
00508       indices[0] = i + offset;  // global column
00509       crsW.ReplaceGlobalValues(
00510         i + offset // GlobalRow
00511         ,1 // NumEntries
00512         ,values // Values
00513         ,indices // Indices
00514         );
00515     }
00516   }
00517 
00518   if (DfDp_out) {
00519     Epetra_MultiVector &DfDp = *DfDp_out;
00520     DfDp.PutScalar(0.0);
00521     if (implicit_) {
00522       for( int i = 0; i < localNumElements; ++i ) {
00523         DfDp[coeff_s_idx(i)][i]
00524           = - d_evalR_d_s(t,gamma[i],coeff_s(i));
00525       }
00526     }
00527     else {
00528       for( int i = 0; i < localNumElements; ++i ) {
00529         DfDp[coeff_s_idx(i)][i]
00530           = + d_evalR_d_s(t,gamma[i],coeff_s(i));
00531       }
00532     }
00533   }
00534 
00535   if (g_out) {
00536     *g_out = *getExactSolution(t,&*coeff_s_p_);
00537     // Note: Above will wipe out coeff_s_p_ as a side effect!
00538   }
00539 
00540   if (DgDp_out) {
00541     *DgDp_out = *getExactSensSolution(t,&*coeff_s_p_);
00542     // Note: Above will wipe out coeff_s_p_ as a side effect!
00543   }
00544 
00545   // Warning: From here on out coeff_s_p_ is unset!
00546   
00547 }
00548 
00549 
00550 // private
00551 
00552 
00553 void DiagonalTransientModel::initialize()
00554 {
00555 
00556   using Teuchos::rcp;
00557   using Teuchos::as;
00558 
00559   //
00560   // Setup map
00561   //
00562 
00563   const int numProcs = epetra_comm_->NumProc();
00564   const int procRank = epetra_comm_->MyPID();
00565   epetra_map_ = rcp( new Epetra_Map( numElements_ * numProcs, 0, *epetra_comm_ ) );
00566 
00567   //
00568   // Create gamma
00569   //
00570 
00571   gamma_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_));
00572   Epetra_Vector &gamma = *gamma_;
00573   switch(gamma_fit_) {
00574     case GAMMA_FIT_LINEAR: {
00575       const int N = gamma.GlobalLength();
00576       const double slope = (gamma_max_ - gamma_min_)/(N-1);
00577       const int MyLength = gamma.MyLength();
00578       if (1==MyLength) {
00579         gamma[0] = gamma_min_;
00580       }
00581       else {
00582         for ( int i=0 ; i<MyLength ; ++i )
00583         {
00584           gamma[i] = slope*(procRank*MyLength+i)+gamma_min_;
00585         }
00586       }
00587       break;
00588     }
00589     case GAMMA_FIT_RANDOM: {
00590       unsigned int seed = Teuchos::ScalarTraits<int>::random()+10*procRank; 
00591       seed *= seed;
00592       gamma.SetSeed(seed);
00593       gamma.Random(); // fill with random numbers in (-1,1)
00594       // Scale random numbers to (gamma_min_,gamma_max_)
00595       const double slope = (gamma_min_ - gamma_max_)/2.0;
00596       const double tmp = (gamma_max_ + gamma_min_)/2.0;
00597       int MyLength = gamma.MyLength();
00598       for (int i=0 ; i<MyLength ; ++i)
00599       {
00600         gamma[i] = slope*gamma[i] + tmp;
00601       }
00602       break;
00603     }
00604     default:
00605       TEST_FOR_EXCEPT("Should never get here!");
00606   }
00607 
00608   //
00609   // Setup for parameters
00610   //
00611 
00612   Np_ = 1;
00613   np_ = coeff_s_.size();
00614   map_p_.clear();
00615   map_p_.push_back(
00616     rcp( new Epetra_LocalMap(np_,0,*epetra_comm_) )
00617     );
00618   names_p_.clear();
00619   {
00620     Teuchos::RCP<Teuchos::Array<std::string> >
00621       names_p_0 = Teuchos::rcp(new Teuchos::Array<std::string>());
00622     for ( int i = 0; i < np_; ++i )
00623       names_p_0->push_back("coeff_s("+Teuchos::toString(i)+")");
00624     names_p_.push_back(names_p_0);
00625   }
00626 
00627   coeff_s_idx_.clear();
00628   const int num_func_per_coeff_s = numElements_ / np_;
00629   const int num_func_per_coeff_s_rem = numElements_ % np_;
00630   for ( int coeff_s_idx_i = 0; coeff_s_idx_i < np_; ++coeff_s_idx_i ) {
00631     const int num_func_per_coeff_s_idx_i
00632       = num_func_per_coeff_s
00633       + ( coeff_s_idx_i < num_func_per_coeff_s_rem ? 1 : 0 );
00634     for (
00635       int coeff_s_idx_i_j = 0;
00636       coeff_s_idx_i_j < num_func_per_coeff_s_idx_i; 
00637       ++coeff_s_idx_i_j
00638       )
00639     {
00640       coeff_s_idx_.push_back(coeff_s_idx_i);
00641     }
00642   }
00643 #ifdef TEUCHOS_DEBUG
00644   TEST_FOR_EXCEPT(
00645     ( as<int>(coeff_s_idx_.size()) != numElements_ ) && "Internal programming error!" );
00646 #endif
00647 
00648   //
00649   // Setup exact solution as response function
00650   //
00651 
00652   if (exactSolutionAsResponse_) {
00653     Ng_ = 1;
00654     map_g_.clear();
00655     map_g_.push_back(
00656       rcp( new Epetra_LocalMap(1,0,*epetra_comm_) )
00657       );
00658   }
00659   else {
00660     Ng_ = 0;
00661   }
00662 
00663   //
00664   // Setup graph for W
00665   //
00666 
00667   if(implicit_) {
00668 
00669     W_graph_ = rcp(
00670       new Epetra_CrsGraph(
00671         ::Copy,*epetra_map_,
00672         1 // elements per row
00673         )
00674       );
00675 
00676     int indices[1];
00677     const int offset = procRank*numElements_ + epetra_map_->IndexBase(); 
00678 
00679     for( int i = 0; i < numElements_; ++i ) {
00680       indices[0] = i + offset;  // global column
00681       W_graph_->InsertGlobalIndices(
00682         i + offset // GlobalRow
00683         ,1 // NumEntries
00684         ,indices // Indices
00685         );
00686     }
00687 
00688     W_graph_->FillComplete();
00689 
00690   }
00691 
00692   //
00693   // Setup initial guess
00694   //
00695 
00696   // Set x_init
00697   x_init_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
00698   x_init_->PutScalar(x0_);
00699 
00700   // Set x_dot_init to provide for a consistent inital guess for implicit mode
00701   // such that f(x_dot,x) = 0!
00702   if (implicit_) {
00703     x_dot_init_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
00704     InArgs inArgs = this->createInArgs();
00705     inArgs.set_x(x_init_);
00706     inArgs.set_t(0.0);
00707     OutArgs outArgs = this->createOutArgs();
00708     outArgs.set_f(x_dot_init_);
00709     this->evalModel(inArgs,outArgs);
00710     x_dot_init_->Scale(-1.0);
00711   }
00712 
00713   // Set p_init
00714   p_init_.clear();
00715   p_init_.push_back(
00716     rcp( new Epetra_Vector( ::Copy, *map_p_[0], &coeff_s_[0] ) )
00717     );
00718 
00719 }
00720 
00721 
00722 void DiagonalTransientModel::set_coeff_s_p( 
00723   const Teuchos::RCP<const Epetra_Vector> &coeff_s_p
00724   ) const
00725 {
00726   if (!is_null(coeff_s_p))
00727     coeff_s_p_ = coeff_s_p;
00728   else
00729     unset_coeff_s_p();
00730 }
00731 
00732 
00733 void DiagonalTransientModel::unset_coeff_s_p() const
00734 {
00735   using Teuchos::as;
00736 #ifdef TEUCHOS_DEBUG
00737   TEUCHOS_ASSERT(
00738     as<int>(get_p_map(0)->NumGlobalElements()) == as<int>(coeff_s_.size()) );
00739 #endif
00740   coeff_s_p_ = Teuchos::rcp(
00741     new Epetra_Vector(
00742       ::View,
00743       *get_p_map(0),
00744       const_cast<double*>(&coeff_s_[0])
00745       )
00746     );
00747   // Note: The above const cast is okay since the coeff_s_p_ RCP is to a const
00748   // Epetr_Vector!
00749 }
00750 
00751 
00752 
00753 } // namespace EpetraExt
00754 
00755 
00756 // Nonmembers
00757 
00758 
00759 Teuchos::RCP<EpetraExt::DiagonalTransientModel>
00760 EpetraExt::diagonalTransientModel(
00761   Teuchos::RCP<Epetra_Comm> const& epetra_comm,
00762   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00763   )
00764 {
00765   Teuchos::RCP<DiagonalTransientModel> diagonalTransientModel =
00766     Teuchos::rcp(new DiagonalTransientModel(epetra_comm));
00767   if (!is_null(paramList))
00768     diagonalTransientModel->setParameterList(paramList);
00769   return diagonalTransientModel;
00770 }

Generated on Wed May 12 21:40:37 2010 for EpetraExt by  doxygen 1.4.7