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