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::RefCountPtr;
00046 
00047 
00048 const std::string Implicit_name = "Implicit";
00049 const bool Implicit_default = true;
00050 
00051 const std::string Lambda_min_name = "Lambda_min";
00052 const double Lambda_min_default = -0.9;
00053 
00054 const std::string Lambda_max_name = "Lambda_max";
00055 const double Lambda_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 Lambda_fit_name = "Lambda_fit";
00061 const std::string Lambda_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& lambda, const double& s )
00075 {
00076   return (exp(lambda*t)*sin(s*t));
00077 }
00078 
00079 
00080 inline
00081 double d_evalR_d_s( const double& t, const double& lambda, const double& s )
00082 {
00083   return (exp(lambda*t)*cos(s*t)*t);
00084 }
00085 
00086 
00087 inline
00088 double f_func(
00089   const double& x, const double& t, const double& lambda, const double& s
00090   )
00091 {
00092   return ( lambda*x + evalR(t,lambda,s) );
00093 }
00094 
00095 
00096 inline
00097 double x_exact(
00098   const double& x0, const double& t, const double& lambda, const double& s
00099   )
00100 {
00101   if ( s == 0.0 )
00102     return ( x0 * exp(lambda*t) );
00103   return ( exp(lambda*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& lambda, const double& s
00114   )
00115 {
00116   if ( s == 0.0 )
00117     return 0.0;
00118   return ( -exp(lambda*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 RefCountPtr<RefCountPtr<const Epetra_Vector> > &vec
00131     )
00132     {
00133       setVector(vec);
00134     }
00135   void setVector( const RefCountPtr<RefCountPtr<const Epetra_Vector> > &vec )
00136     {
00137       vec_ = vec;
00138     }
00139 private:
00140   RefCountPtr<RefCountPtr<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::RefCountPtr<Epetra_Comm> const& epetra_comm
00155   )
00156   : epetra_comm_(epetra_comm.assert_not_null()),
00157     implicit_(Implicit_default),
00158     numElements_(NumElements_default),
00159     lambda_min_(Lambda_min_default),
00160     lambda_max_(Lambda_max_default),
00161     coeff_s_(Teuchos::fromStringToArray<double>(Coeff_s_default)),
00162     lambda_fit_(LAMBDA_FIT_LINEAR), // Must be the same as Lambda_fit_default!
00163     x0_(x0_default),
00164     exactSolutionAsResponse_(ExactSolutionAsResponse_default),
00165     isIntialized_(false)
00166 {
00167   initialize();
00168 }
00169 
00170 
00171 Teuchos::RefCountPtr<const Epetra_Vector>
00172 DiagonalTransientModel::getExactSolution(
00173   const double t, const Epetra_Vector *coeff_s_p
00174   ) const
00175 {
00176   set_coeff_s_p(Teuchos::rcp(coeff_s_p,false));
00177   UnsetParameterVector unsetParameterVector(Teuchos::rcp(&coeff_s_p_,false));
00178   Teuchos::RefCountPtr<Epetra_Vector>
00179     x_star_ptr = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
00180   Epetra_Vector& x_star = *x_star_ptr;
00181   Epetra_Vector& lambda = *lambda_;
00182   int myN = x_star.MyLength();
00183   for ( int i=0 ; i<myN ; ++i ) {
00184     x_star[i] = x_exact( x0_, t, lambda[i], coeff_s(i) );
00185   }
00186   return(x_star_ptr);
00187 }
00188 
00189 
00190 Teuchos::RefCountPtr<const Epetra_MultiVector>
00191 DiagonalTransientModel::getExactSensSolution(
00192   const double t, const Epetra_Vector *coeff_s_p
00193   ) const
00194 {
00195   set_coeff_s_p(Teuchos::rcp(coeff_s_p,false));
00196   UnsetParameterVector unsetParameterVector(Teuchos::rcp(&coeff_s_p_,false));
00197   Teuchos::RefCountPtr<Epetra_MultiVector>
00198     dxds_star_ptr = Teuchos::rcp(new Epetra_MultiVector(*epetra_map_,np_,false));
00199   Epetra_MultiVector& dxds_star = *dxds_star_ptr;
00200   dxds_star.PutScalar(0.0);
00201   Epetra_Vector& lambda = *lambda_;
00202   int myN = dxds_star.MyLength();
00203   for ( int i=0 ; i<myN ; ++i ) {
00204     const int coeff_s_idx = this->coeff_s_idx(i);
00205     (*dxds_star(coeff_s_idx))[i] = dxds_exact( t, lambda[i], coeff_s(i) );
00206     // Note: Above, at least the column access will be validated in debug mode
00207     // but the row index i will not ever be.  Perhaps we can augment Epetra to
00208     // fix this?
00209   }
00210   return (dxds_star_ptr);
00211 }
00212 
00213 
00214 // Overridden from ParameterListAcceptor
00215 
00216 
00217 void DiagonalTransientModel::setParameterList(
00218   Teuchos::RefCountPtr<Teuchos::ParameterList> const& paramList
00219   )
00220 {
00221   using Teuchos::get; using Teuchos::getIntegralValue;
00222   using Teuchos::getArrayFromStringParameter;
00223   TEST_FOR_EXCEPT( is_null(paramList) );
00224   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00225   isIntialized_ = false;
00226   paramList_ = paramList;
00227   implicit_ = get<bool>(*paramList_,Implicit_name);
00228   numElements_ = get<int>(*paramList_,NumElements_name);
00229   lambda_min_ = get<double>(*paramList_,Lambda_min_name);
00230   lambda_max_ = get<double>(*paramList_,Lambda_max_name);
00231   coeff_s_ = getArrayFromStringParameter<double>(*paramList_,Coeff_s_name);
00232   lambda_fit_ = getIntegralValue<ELambdaFit>(*paramList_,Lambda_fit_name);
00233   x0_ = get<double>(*paramList_,x0_name);
00234   exactSolutionAsResponse_ = get<bool>(*paramList_,ExactSolutionAsResponse_name);
00235   initialize();
00236 }
00237 
00238 
00239 Teuchos::RefCountPtr<Teuchos::ParameterList> 
00240 DiagonalTransientModel::getParameterList()
00241 {
00242   return paramList_;
00243 }
00244 
00245 
00246 Teuchos::RefCountPtr<Teuchos::ParameterList>
00247 DiagonalTransientModel::unsetParameterList()
00248 {
00249   Teuchos::RefCountPtr<Teuchos::ParameterList> _paramList = paramList_;
00250   paramList_ = Teuchos::null;
00251   return _paramList;
00252 }
00253 
00254 
00255 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00256 DiagonalTransientModel::getParameterList() const
00257 {
00258   return paramList_;
00259 }
00260 
00261 
00262 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00263 DiagonalTransientModel::getValidParameters() const
00264 {
00265   using Teuchos::RefCountPtr; using Teuchos::ParameterList;
00266   using Teuchos::tuple;
00267   using Teuchos::setIntParameter; using Teuchos::setDoubleParameter;
00268   using Teuchos::setStringToIntegralParameter;
00269   static RefCountPtr<const ParameterList> validPL;
00270   if (is_null(validPL)) {
00271     RefCountPtr<ParameterList> pl = Teuchos::parameterList();
00272     pl->set(Implicit_name,true);
00273     setDoubleParameter(
00274       Lambda_min_name, Lambda_min_default, "",
00275       &*pl
00276       );
00277     setDoubleParameter(
00278       Lambda_max_name, Lambda_max_default, "",
00279       &*pl
00280       );
00281     setStringToIntegralParameter(
00282       Lambda_fit_name, Lambda_fit_default, "",
00283       tuple<std::string>("Linear","Random"),
00284       tuple<ELambdaFit>(LAMBDA_FIT_LINEAR,LAMBDA_FIT_RANDOM),
00285       &*pl
00286       );
00287     setIntParameter(
00288       NumElements_name, NumElements_default, "",
00289       &*pl
00290       );
00291     setDoubleParameter(
00292       x0_name, x0_default, "",
00293       &*pl
00294       );
00295     pl->set( Coeff_s_name, Coeff_s_default );
00296     pl->set( ExactSolutionAsResponse_name, ExactSolutionAsResponse_default );
00297     validPL = pl;
00298   }
00299   return validPL;
00300 }
00301 
00302 
00303 // Overridden from EpetraExt::ModelEvaluator
00304 
00305 
00306 Teuchos::RefCountPtr<const Epetra_Map>
00307 DiagonalTransientModel::get_x_map() const
00308 {
00309   return epetra_map_;
00310 }
00311 
00312 
00313 Teuchos::RefCountPtr<const Epetra_Map>
00314 DiagonalTransientModel::get_f_map() const
00315 {
00316   return epetra_map_;
00317 }
00318 
00319 
00320 Teuchos::RefCountPtr<const Epetra_Map>
00321 DiagonalTransientModel::get_p_map(int l) const
00322 {
00323 #ifdef TEUCHOS_DEBUG
00324   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
00325 #endif
00326   return map_p_[l];
00327 }
00328 
00329 
00330 Teuchos::RefCountPtr<const Teuchos::Array<std::string> >
00331 DiagonalTransientModel::get_p_names(int l) const
00332 {
00333 #ifdef TEUCHOS_DEBUG
00334   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( l, 0, Np_ );
00335 #endif
00336   return names_p_[l];
00337 }
00338 
00339 
00340 Teuchos::RefCountPtr<const Epetra_Map>
00341 DiagonalTransientModel::get_g_map(int j) const
00342 {
00343 #ifdef TEUCHOS_DEBUG
00344   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( j, 0, Ng_ );
00345 #endif
00346   return map_g_[j];
00347 }
00348 
00349 
00350 Teuchos::RefCountPtr<const Epetra_Vector>
00351 DiagonalTransientModel::get_x_init() const
00352 {
00353   return x_init_;
00354 }
00355 
00356 
00357 Teuchos::RefCountPtr<const Epetra_Vector>
00358 DiagonalTransientModel::get_x_dot_init() const
00359 {
00360   return x_dot_init_;
00361 }
00362 
00363 
00364 Teuchos::RefCountPtr<const Epetra_Vector>
00365 DiagonalTransientModel::get_p_init(int l) const
00366 {
00367 #ifdef TEUCHOS_DEBUG
00368   TEUCHOS_ASSERT( l == 0 );
00369 #endif
00370   return p_init_[l];
00371 }
00372 
00373 
00374 Teuchos::RefCountPtr<Epetra_Operator>
00375 DiagonalTransientModel::create_W() const
00376 {
00377   if(implicit_)
00378     return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00379   return Teuchos::null;
00380 }
00381 
00382 
00383 EpetraExt::ModelEvaluator::InArgs
00384 DiagonalTransientModel::createInArgs() const
00385 {
00386   InArgsSetup inArgs;
00387   inArgs.set_Np(Np_);
00388   inArgs.setSupports(IN_ARG_t,true);
00389   inArgs.setSupports(IN_ARG_x,true);
00390   if(implicit_) {
00391     inArgs.setSupports(IN_ARG_x_dot,true);
00392     inArgs.setSupports(IN_ARG_alpha,true);
00393     inArgs.setSupports(IN_ARG_beta,true);
00394   }
00395   return inArgs;
00396 }
00397 
00398 
00399 EpetraExt::ModelEvaluator::OutArgs
00400 DiagonalTransientModel::createOutArgs() const
00401 {
00402   OutArgsSetup outArgs;
00403   outArgs.set_Np_Ng(Np_,Ng_);
00404   outArgs.setSupports(OUT_ARG_f,true);
00405   if(implicit_) {
00406     outArgs.setSupports(OUT_ARG_W,true);
00407     outArgs.set_W_properties(
00408       DerivativeProperties(
00409         DERIV_LINEARITY_NONCONST
00410         ,DERIV_RANK_FULL
00411         ,true // supportsAdjoint
00412         )
00413       );
00414   }
00415   outArgs.setSupports(OUT_ARG_DfDp,0,DERIV_MV_BY_COL);
00416   outArgs.set_DfDp_properties(
00417     0,DerivativeProperties(
00418       DERIV_LINEARITY_NONCONST,
00419       DERIV_RANK_DEFICIENT,
00420       true // supportsAdjoint
00421       )
00422     );
00423   if (exactSolutionAsResponse_) {
00424     outArgs.setSupports(OUT_ARG_DgDp,0,0,DERIV_MV_BY_COL);
00425     outArgs.set_DgDp_properties(
00426       0,0,DerivativeProperties(
00427         DERIV_LINEARITY_NONCONST,
00428         DERIV_RANK_DEFICIENT,
00429         true // supportsAdjoint
00430         )
00431       );
00432   }
00433   return outArgs;
00434 }
00435 
00436 
00437 void DiagonalTransientModel::evalModel(
00438   const InArgs& inArgs, const OutArgs& outArgs
00439   ) const
00440 {
00441 
00442   using Teuchos::rcp;
00443   using Teuchos::RefCountPtr;
00444   using Teuchos::null;
00445   using Teuchos::dyn_cast;
00446 
00447   const Epetra_Vector &x = *(inArgs.get_x());
00448   const double t = inArgs.get_t();
00449   if (Np_) set_coeff_s_p(inArgs.get_p(0)); // Sets for coeff_s(...) function!
00450   UnsetParameterVector unsetParameterVector(rcp(&coeff_s_p_,false));
00451   // Note: Above, the destructor for unsetParameterVector will ensure that the
00452   // RCP to the parameter vector will be unset no matter if an exception is
00453   // thrown or not.
00454 
00455   Epetra_Operator *W_out = ( implicit_ ? outArgs.get_W().get() : 0 );
00456   Epetra_Vector *f_out = outArgs.get_f().get();
00457   Epetra_MultiVector *DfDp_out = 0;
00458   if (Np_) DfDp_out = get_DfDp_mv(0,outArgs).get();
00459   Epetra_Vector *g_out = 0;
00460   Epetra_MultiVector *DgDp_out = 0;
00461   if (exactSolutionAsResponse_) {
00462     g_out = outArgs.get_g(0).get();
00463     DgDp_out = get_DgDp_mv(0,0,outArgs,DERIV_MV_BY_COL).get();
00464   }
00465 
00466   const Epetra_Vector &lambda = *lambda_;
00467 
00468   int localNumElements = x.MyLength();
00469 
00470   if (f_out) {
00471     Epetra_Vector &f = *f_out;
00472     if (implicit_) {
00473       const Epetra_Vector *x_dot_in = inArgs.get_x_dot().get();
00474       if (x_dot_in) {
00475         const Epetra_Vector &x_dot = *x_dot_in;
00476         for ( int i=0 ; i<localNumElements ; ++i )
00477           f[i] = x_dot[i] - f_func(x[i],t,lambda[i],coeff_s(i));
00478       }
00479       else {
00480         for ( int i=0 ; i<localNumElements ; ++i )
00481           f[i] = - f_func(x[i],t,lambda[i],coeff_s(i));
00482       }
00483     }
00484     else {
00485       for ( int i=0 ; i<localNumElements ; ++i )
00486         f[i] = f_func(x[i],t,lambda[i],coeff_s(i));
00487     }
00488   }
00489 
00490   if ( W_out ) {
00491     // If we get here then we are in implicit mode!
00492     const double alpha = inArgs.get_alpha();
00493     const double beta = inArgs.get_beta();
00494     Epetra_CrsMatrix &crsW = dyn_cast<Epetra_CrsMatrix>(*W_out);
00495     double values[1];
00496     int indices[1];
00497     const int offset
00498       = epetra_comm_->MyPID()*localNumElements + epetra_map_->IndexBase(); 
00499     for( int i = 0; i < localNumElements; ++i ) {
00500       values[0] = alpha - beta*lambda[i];
00501       indices[0] = i + offset;  // global column
00502       crsW.ReplaceGlobalValues(
00503         i + offset // GlobalRow
00504         ,1 // NumEntries
00505         ,values // Values
00506         ,indices // Indices
00507         );
00508     }
00509   }
00510 
00511   if (DfDp_out) {
00512     Epetra_MultiVector &DfDp = *DfDp_out;
00513     DfDp.PutScalar(0.0);
00514     if (implicit_) {
00515       for( int i = 0; i < localNumElements; ++i ) {
00516         DfDp[coeff_s_idx(i)][i]
00517           = - d_evalR_d_s(t,lambda[i],coeff_s(i));
00518       }
00519     }
00520     else {
00521       for( int i = 0; i < localNumElements; ++i ) {
00522         DfDp[coeff_s_idx(i)][i]
00523           = + d_evalR_d_s(t,lambda[i],coeff_s(i));
00524       }
00525     }
00526   }
00527 
00528   if (g_out) {
00529     *g_out = *getExactSolution(t,&*coeff_s_p_);
00530     // Note: Above will wipe out coeff_s_p_ as a side effect!
00531   }
00532 
00533   if (DgDp_out) {
00534     *DgDp_out = *getExactSensSolution(t,&*coeff_s_p_);
00535     // Note: Above will wipe out coeff_s_p_ as a side effect!
00536   }
00537 
00538   // Warning: From here on out coeff_s_p_ is unset!
00539   
00540 }
00541 
00542 
00543 // private
00544 
00545 
00546 void DiagonalTransientModel::initialize()
00547 {
00548 
00549   using Teuchos::rcp;
00550   using Teuchos::as;
00551 
00552   //
00553   // Setup map
00554   //
00555 
00556   const int numProcs = epetra_comm_->NumProc();
00557   const int procRank = epetra_comm_->MyPID();
00558   epetra_map_ = rcp( new Epetra_Map( numElements_ * numProcs, 0, *epetra_comm_ ) );
00559 
00560   //
00561   // Create lambda
00562   //
00563 
00564   lambda_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_));
00565   Epetra_Vector &lambda = *lambda_;
00566   switch(lambda_fit_) {
00567     case LAMBDA_FIT_LINEAR: {
00568       const int N = lambda.GlobalLength();
00569       const double slope = (lambda_max_ - lambda_min_)/(N-1);
00570       const int MyLength = lambda.MyLength();
00571       if (1==MyLength) {
00572         lambda[0] = lambda_min_;
00573       }
00574       else {
00575         for ( int i=0 ; i<MyLength ; ++i )
00576         {
00577           lambda[i] = slope*(procRank*MyLength+i)+lambda_min_;
00578         }
00579       }
00580       break;
00581     }
00582     case LAMBDA_FIT_RANDOM: {
00583       unsigned int seed = Teuchos::ScalarTraits<int>::random()+10*procRank; 
00584       seed *= seed;
00585       lambda.SetSeed(seed);
00586       lambda.Random(); // fill with random numbers in (-1,1)
00587       // Scale random numbers to (lambda_min_,lambda_max_)
00588       const double slope = (lambda_min_ - lambda_max_)/2.0;
00589       const double tmp = (lambda_max_ + lambda_min_)/2.0;
00590       int MyLength = lambda.MyLength();
00591       for (int i=0 ; i<MyLength ; ++i)
00592       {
00593         lambda[i] = slope*lambda[i] + tmp;
00594       }
00595       break;
00596     }
00597     default:
00598       TEST_FOR_EXCEPT("Should never get here!");
00599   }
00600 
00601   //
00602   // Setup for parameters
00603   //
00604 
00605   Np_ = 1;
00606   np_ = coeff_s_.size();
00607   map_p_.clear();
00608   map_p_.push_back(
00609     rcp( new Epetra_LocalMap(np_,0,*epetra_comm_) )
00610     );
00611   names_p_.clear();
00612   {
00613     Teuchos::RefCountPtr<Teuchos::Array<std::string> >
00614       names_p_0 = Teuchos::rcp(new Teuchos::Array<std::string>());
00615     for ( int i = 0; i < np_; ++i )
00616       names_p_0->push_back("coeff_s("+Teuchos::toString(i)+")");
00617     names_p_.push_back(names_p_0);
00618   }
00619 
00620   coeff_s_idx_.clear();
00621   const int num_func_per_coeff_s = numElements_ / np_;
00622   const int num_func_per_coeff_s_rem = numElements_ % np_;
00623   for ( int coeff_s_idx_i = 0; coeff_s_idx_i < np_; ++coeff_s_idx_i ) {
00624     const int num_func_per_coeff_s_idx_i
00625       = num_func_per_coeff_s
00626       + ( coeff_s_idx_i < num_func_per_coeff_s_rem ? 1 : 0 );
00627     for (
00628       int coeff_s_idx_i_j = 0;
00629       coeff_s_idx_i_j < num_func_per_coeff_s_idx_i; 
00630       ++coeff_s_idx_i_j
00631       )
00632     {
00633       coeff_s_idx_.push_back(coeff_s_idx_i);
00634     }
00635   }
00636 #ifdef TEUCHOS_DEBUG
00637   TEST_FOR_EXCEPT(
00638     ( as<int>(coeff_s_idx_.size()) != numElements_ ) && "Internal programming error!" );
00639 #endif
00640 
00641   //
00642   // Setup exact solution as response function
00643   //
00644 
00645   if (exactSolutionAsResponse_) {
00646     Ng_ = 1;
00647     map_g_.clear();
00648     map_g_.push_back(
00649       rcp( new Epetra_LocalMap(1,0,*epetra_comm_) )
00650       );
00651   }
00652   else {
00653     Ng_ = 0;
00654   }
00655 
00656   //
00657   // Setup graph for W
00658   //
00659 
00660   if(implicit_) {
00661 
00662     W_graph_ = rcp(
00663       new Epetra_CrsGraph(
00664         ::Copy,*epetra_map_,
00665         1 // elements per row
00666         )
00667       );
00668 
00669     int indices[1];
00670     const int offset = procRank*numElements_ + epetra_map_->IndexBase(); 
00671 
00672     for( int i = 0; i < numElements_; ++i ) {
00673       indices[0] = i + offset;  // global column
00674       W_graph_->InsertGlobalIndices(
00675         i + offset // GlobalRow
00676         ,1 // NumEntries
00677         ,indices // Indices
00678         );
00679     }
00680 
00681     W_graph_->FillComplete();
00682 
00683   }
00684 
00685   //
00686   // Setup initial guess
00687   //
00688 
00689   // Set x_init
00690   x_init_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
00691   x_init_->PutScalar(x0_);
00692 
00693   // Set x_dot_init to provide for a consistent inital guess for implicit mode
00694   // such that f(x_dot,x) = 0!
00695   if (implicit_) {
00696     x_dot_init_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_,false));
00697     InArgs inArgs = this->createInArgs();
00698     inArgs.set_x(x_init_);
00699     inArgs.set_t(0.0);
00700     OutArgs outArgs = this->createOutArgs();
00701     outArgs.set_f(x_dot_init_);
00702     this->evalModel(inArgs,outArgs);
00703     x_dot_init_->Scale(-1.0);
00704   }
00705 
00706   // Set p_init
00707   p_init_.clear();
00708   p_init_.push_back(
00709     rcp( new Epetra_Vector( ::Copy, *map_p_[0], &coeff_s_[0] ) )
00710     );
00711 
00712 }
00713 
00714 
00715 void DiagonalTransientModel::set_coeff_s_p( 
00716   const Teuchos::RefCountPtr<const Epetra_Vector> &coeff_s_p
00717   ) const
00718 {
00719   if (!is_null(coeff_s_p))
00720     coeff_s_p_ = coeff_s_p;
00721   else
00722     unset_coeff_s_p();
00723 }
00724 
00725 
00726 void DiagonalTransientModel::unset_coeff_s_p() const
00727 {
00728   using Teuchos::as;
00729 #ifdef TEUCHOS_DEBUG
00730   TEUCHOS_ASSERT(
00731     as<int>(get_p_map(0)->NumGlobalElements()) == as<int>(coeff_s_.size()) );
00732 #endif
00733   coeff_s_p_ = Teuchos::rcp(
00734     new Epetra_Vector(
00735       ::View,
00736       *get_p_map(0),
00737       const_cast<double*>(&coeff_s_[0])
00738       )
00739     );
00740   // Note: The above const cast is okay since the coeff_s_p_ RCP is to a const
00741   // Epetr_Vector!
00742 }
00743 
00744 
00745 
00746 } // namespace EpetraExt
00747 
00748 
00749 // Nonmembers
00750 
00751 
00752 Teuchos::RefCountPtr<EpetraExt::DiagonalTransientModel>
00753 EpetraExt::diagonalTransientModel(
00754   Teuchos::RefCountPtr<Epetra_Comm> const& epetra_comm,
00755   Teuchos::RefCountPtr<Teuchos::ParameterList> const& paramList
00756   )
00757 {
00758   Teuchos::RefCountPtr<DiagonalTransientModel>
00759     diagonalTransientModel
00760     = Teuchos::rcp(new DiagonalTransientModel(epetra_comm));
00761   if (!is_null(paramList))
00762     diagonalTransientModel->setParameterList(paramList);
00763   return diagonalTransientModel;
00764 }
00765 

Generated on Tue Oct 20 12:45:29 2009 for EpetraExt by doxygen 1.4.7