ExampleApplication.cpp

00001 //@HEADER
00002 
00003 // ***********************************************************************
00004 //
00005 //                     Rythmos Package
00006 //                 Copyright (2006) Sandia Corporation
00007 //
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 //
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00026 //
00027 // ***********************************************************************
00028 //@HEADER
00029 
00030 #include "ExampleApplication.hpp"
00031 #include "Teuchos_ScalarTraits.hpp"
00032 
00033 #ifdef HAVE_MPI
00034 #include "Epetra_MpiComm.h"
00035 #include "mpi.h"
00036 #else
00037 #include "Epetra_SerialComm.h"
00038 #endif // HAVE_MPI
00039 
00040 #include "Epetra_CrsMatrix.h"
00041 
00042 #ifdef EXAMPLEAPPLICATION_DEBUG
00043 #include <iostream>
00044 #endif
00045 
00046 ExampleApplication::ExampleApplication(Teuchos::RCP<Epetra_Comm> &epetra_comm_ptr_, Teuchos::ParameterList &params)
00047 {
00048   implicit_ = params.get<bool>( "Implicit" );
00049   lambda_min_ = params.get<double>( "Lambda_min" );
00050   lambda_max_ = params.get<double>( "Lambda_max" );
00051   lambda_fit_ = params.get<std::string>( "Lambda_fit" );
00052   numElements_ = params.get<int>( "NumElements" );
00053   x0_ = params.get<double>( "x0" );
00054   coeff_s_ = params.get<double>( "Coeff_s" );
00055 
00056   // Construct a Map with NumElements and index base of 0
00057   epetra_map_ptr_ = Teuchos::rcp( new Epetra_Map(numElements_, 0, *epetra_comm_ptr_) );
00058 
00059   lambda_ptr_ = Teuchos::rcp(new Epetra_Vector(*epetra_map_ptr_));
00060   Epetra_Vector &lambda = *lambda_ptr_;
00061   if ( lambda_fit_ == "linear" )
00062   {
00063     int N = lambda.GlobalLength();
00064     double tmp;
00065     if (N == 1)
00066       tmp = 0; // lambda[0] = lambda_min_
00067     else
00068       tmp = (lambda_max_ - lambda_min_)/(N-1);
00069     int MyLength = lambda.MyLength();
00070     int MyPID = lambda.Comm().MyPID();
00071     for (int i=0 ; i<MyLength ; ++i)
00072     {
00073       lambda[i] = tmp*(MyPID*MyLength+i)+lambda_min_;
00074     }
00075   }
00076   else // if ( lambda_fit_ == "random" )
00077   {
00078     int MyPID = lambda.Comm().MyPID();
00079     unsigned int seed = Teuchos::ScalarTraits<int>::random()+10*MyPID; 
00080     seed *= seed;
00081     lambda.SetSeed(seed);
00082     lambda.Random(); // fill with random numbers in (-1,1)
00083     // Scale random numbers to (lambda_min_,lambda_max_)
00084     double slope = (lambda_min_ - lambda_max_)/2.0;
00085     double tmp = (lambda_max_ + lambda_min_)/2.0;
00086     int MyLength = lambda.MyLength();
00087     for (int i=0 ; i<MyLength ; ++i)
00088     {
00089       lambda[i] = slope*lambda[i] + tmp;
00090     }
00091   }
00092 
00093   if(implicit_) {
00094     int localNumElements = lambda.MyLength();
00095     W_graph_ = Teuchos::rcp(new Epetra_CrsGraph(::Copy,*epetra_map_ptr_,1));
00096     int indices[1];
00097     const int IB = epetra_map_ptr_->IndexBase();
00098     for( int i = 0; i < localNumElements; ++i ) {
00099       indices[0] = i + IB;  // global column
00100       W_graph_->InsertGlobalIndices(
00101         i + IB              // GlobalRow
00102         ,1                  // NumEntries
00103         ,indices            // Indices
00104         );
00105     }
00106     W_graph_->FillComplete();
00107   }
00108  
00109 }
00110 
00111 Teuchos::RCP<const Epetra_Vector> ExampleApplication::get_coeff() const
00112 {
00113   return(lambda_ptr_);
00114 }
00115 
00116 Teuchos::RCP<const Epetra_Vector> ExampleApplication::getExactSolution(double t) const
00117 {
00118   Teuchos::RCP<Epetra_Vector> x_star_ptr = Teuchos::rcp(new Epetra_Vector(*epetra_map_ptr_));
00119   Epetra_Vector& x_star = *x_star_ptr;
00120   Epetra_Vector& lambda = *lambda_ptr_;
00121   x_star.PutScalar(0.0);
00122   x_star.Scale(t,lambda);
00123   int myN = x_star.MyLength();
00124   if (coeff_s_ == 0.0)
00125   {
00126     for ( int i=0 ; i<myN ; ++i )
00127       x_star[i] = x0_*exp(x_star[i]);
00128   }
00129   else
00130   {
00131     for ( int i=0 ; i<myN ; ++i )
00132       x_star[i] = (x0_+(1.0/coeff_s_))*exp(x_star[i]) - exp(x_star[i])*cos(t*coeff_s_)/coeff_s_;
00133   }
00134   return(x_star_ptr);
00135 }
00136 
00137 // Overridden from EpetraExt::ModelEvaluator
00138 
00139 Teuchos::RCP<const Epetra_Map>
00140 ExampleApplication::get_x_map() const
00141 {
00142   return epetra_map_ptr_;
00143 }
00144 
00145 Teuchos::RCP<const Epetra_Map>
00146 ExampleApplication::get_f_map() const
00147 {
00148   return epetra_map_ptr_;
00149 }
00150 
00151 Teuchos::RCP<const Epetra_Vector>
00152 ExampleApplication::get_x_init() const
00153 {
00154   Teuchos::RCP<Epetra_Vector> x_init = Teuchos::rcp(new Epetra_Vector(*epetra_map_ptr_));
00155   x_init->PutScalar(x0_);
00156   return x_init;
00157 }
00158 
00159 Teuchos::RCP<const Epetra_Vector>
00160 ExampleApplication::get_x_dot_init() const
00161 {
00162   Teuchos::RCP<Epetra_Vector> x_dot_init = Teuchos::rcp(new Epetra_Vector(*epetra_map_ptr_));
00163   x_dot_init->PutScalar(0.0);
00164   return x_dot_init;
00165 }
00166 
00167 Teuchos::RCP<Epetra_Operator>
00168 ExampleApplication::create_W() const
00169 {
00170   if(implicit_)
00171     return Teuchos::rcp(new Epetra_CrsMatrix(::Copy,*W_graph_));
00172   return Teuchos::null;
00173 }
00174 
00175 EpetraExt::ModelEvaluator::InArgs
00176 ExampleApplication::createInArgs() const
00177 {
00178   InArgsSetup inArgs;
00179   inArgs.setSupports(IN_ARG_t,true);
00180   inArgs.setSupports(IN_ARG_x,true);
00181   if(implicit_) {
00182     inArgs.setSupports(IN_ARG_x_dot,true);
00183     inArgs.setSupports(IN_ARG_alpha,true);
00184     inArgs.setSupports(IN_ARG_beta,true);
00185   }
00186   return inArgs;
00187 }
00188 
00189 EpetraExt::ModelEvaluator::OutArgs
00190 ExampleApplication::createOutArgs() const
00191 {
00192   OutArgsSetup outArgs;
00193   outArgs.setSupports(OUT_ARG_f,true);
00194   if(implicit_) {
00195     outArgs.setSupports(OUT_ARG_W,true);
00196   }
00197   return outArgs;
00198 }
00199 
00200 void ExampleApplication::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
00201 {
00202   const Epetra_Vector &x = *(inArgs.get_x());
00203   const double t = inArgs.get_t();
00204   const Epetra_Vector &lambda = *lambda_ptr_;
00205 #ifdef EXAMPLEAPPLICATION_DEBUG
00206       std::cout << "----------------------------------------------------------------------" << std::endl;
00207       std::cout << "ExampleApplication::evalModel x = " << std::endl;
00208       x.Print(std::cout);
00209       std::cout << "ExampleApplication::evalModel lambda = " << std::endl;
00210       lambda.Print(std::cout);
00211 #endif
00212   int localNumElements = x.MyLength();
00213   if(implicit_) 
00214   {
00215     const Epetra_Vector &x_dot = *inArgs.get_x_dot();
00216     if(outArgs.get_f().get()) 
00217     {
00218       Epetra_Vector &f = *outArgs.get_f();
00219       for (int i=0 ; i<localNumElements ; ++i)
00220       {
00221         f[i] = x_dot[i] - lambda[i]*x[i] - evalR(t,lambda[i],coeff_s_);
00222       }
00223 #ifdef EXAMPLEAPPLICATION_DEBUG
00224       std::cout << "ExampleApplication::evalModel (implicit) x_dot = " << std::endl;
00225       x_dot.Print(std::cout);
00226       std::cout << "ExampleApplication::evalModel (implicit) f = " << std::endl;
00227       f.Print(std::cout);
00228 #endif
00229     }
00230     Teuchos::RCP<Epetra_Operator> W;
00231     if( (W = outArgs.get_W()).get() ) {
00232       const double alpha = inArgs.get_alpha();
00233       const double beta = inArgs.get_beta();
00234       Epetra_CrsMatrix &crsW = Teuchos::dyn_cast<Epetra_CrsMatrix>(*W);
00235       double values[1];
00236       int indices[1];
00237       const int IB = epetra_map_ptr_->IndexBase();
00238       for( int i = 0; i < localNumElements; ++i ) 
00239       {
00240         values[0] = alpha - beta*lambda[i];
00241         indices[0] = i + IB;  // global column
00242         crsW.ReplaceGlobalValues(i + IB              // GlobalRow
00243                                  ,1                  // NumEntries
00244                                  ,values             // Values
00245                                  ,indices            // Indices
00246                                           );
00247       }
00248 #ifdef EXAMPLEAPPLICATION_DEBUG
00249       std::cout << "ExampleApplication::evalModel (implicit) alpha, beta = " << std::endl;
00250       std::cout << "alpha = " << alpha << std::endl;
00251       std::cout << "beta = "  << beta  << std::endl;
00252       std::cout << "ExampleApplication::evalModel (implicit) W = " << std::endl;
00253       crsW.Print(std::cout);
00254 #endif
00255     }
00256   }
00257   else 
00258   {
00259     Epetra_Vector &f = *outArgs.get_f();
00260     int localNumElements = x.MyLength();
00261     for (int i=0 ; i<localNumElements ; ++i)
00262     {
00263       f[i] = lambda[i]*x[i]+evalR(t,lambda[i],coeff_s_);
00264     }
00265 #ifdef EXAMPLEAPPLICATION_DEBUG
00266     std::cout << "ExampleApplication::evalModel (explicit) f = " << std::endl;
00267     f.Print(std::cout);
00268 #endif
00269   }
00270 #ifdef EXAMPLEAPPLICATION_DEBUG
00271   std::cout << "----------------------------------------------------------------------" << std::endl;
00272 #endif
00273 }
00274 
00275 double ExampleApplication::evalR(const double& t, const double& lambda, const double& s) const
00276 {
00277   return(exp(lambda*t)*sin(s*t));
00278 }
00279 

Generated on Tue Oct 20 12:46:07 2009 for Rythmos - Transient Integration for Differential Equations by doxygen 1.4.7