Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/test/Discretization/Integration/test_17.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) 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 // 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 Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00044 
00051 #include "Intrepid_AdaptiveSparseGrid.hpp"
00052 //#include "Intrepid_CubatureLineSorted.hpp"
00053 #include "Intrepid_Utils.hpp"
00054 #include "Teuchos_oblackholestream.hpp"
00055 #include "Teuchos_RCP.hpp"
00056 #include "Teuchos_RefCountPtr.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 
00059 using namespace Intrepid;
00060 
00061 
00062 template<class Scalar>
00063 class StdVector {
00064 private:
00065   Teuchos::RefCountPtr<std::vector<Scalar> >  std_vec_;
00066 
00067 public:
00068 
00069   StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec ) 
00070   : std_vec_(std_vec) {}
00071 
00072   Teuchos::RefCountPtr<StdVector<Scalar> > Create() const {
00073     return Teuchos::rcp( new StdVector<Scalar>(
00074            Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0))));
00075   }
00076 
00077   void Update( StdVector<Scalar> & s ) {
00078     int dimension  = (int)(std_vec_->size());
00079     for (int i=0; i<dimension; i++)
00080       (*std_vec_)[i] += s[i]; 
00081   }
00082 
00083   void Update( Scalar alpha, StdVector<Scalar> & s )  {
00084     int dimension  = (int)(std_vec_->size());
00085     for (int i=0; i<dimension; i++)
00086       (*std_vec_)[i] += alpha*s[i];
00087   }
00088 
00089   Scalar operator[](int i) {
00090     return (*std_vec_)[i];
00091   }
00092 
00093   void clear() {
00094     std_vec_->clear();
00095   }
00096 
00097   void resize(int n, Scalar p) {
00098     std_vec_->resize(n,p);
00099   }
00100 
00101   int size() {
00102     return (int)std_vec_->size();
00103   }
00104 
00105   void Set( Scalar alpha )  {
00106     int dimension  = (int)(std_vec_->size());
00107     for (int i=0; i<dimension; i++)
00108       (*std_vec_)[i] = alpha;
00109   }
00110 };
00111 
00112 template<class Scalar, class UserVector>
00113 class ASGdata : public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector>
00114 {  
00115 public:  
00116   ~ASGdata() {}
00117 
00118   ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D,
00119           std::vector<EIntrepidGrowth> growth1D, int maxLevel,
00120           bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>(
00121           dimension,rule1D,growth1D,maxLevel,isNormalized) {}
00122 
00123   void eval_integrand(
00124          UserVector & output, 
00125          std::vector<Scalar> & input) {
00126  
00127     output.clear(); output.resize(1,0.0);
00128     int dimension = (int)input.size();
00129     std::vector<Scalar> point = input;
00130     for (int i=0; i<dimension; i++) {
00131       point[i] = 0.5*point[i]+0.5;
00132     }
00133     Teuchos::RCP<std::vector<long double> > etmp = 
00134       Teuchos::rcp(new std::vector<long double>(1,0.0));
00135     StdVector<long double> tmp(etmp); 
00136     Scalar gamma     = 0.5;
00137     Scalar x         = 0.0;
00138     
00139     Scalar prod1     = gamma*(1.0-x);
00140     Scalar prod2     = (1.0-x)*point[0];
00141     
00142     for (int i=1; i<dimension; i++) {
00143       prod1 = powl(gamma*(1.0-x),(long double)i); prod2 = 1.0-x;
00144       for (int j=0; j<i; j++) {
00145         prod2 *= point[j];
00146         if (j<i-1) {
00147           prod1 *= powl(point[j],(long double)(i-(j+1)));
00148         }
00149       }
00150       (*etmp)[0] =  prod1*(1.0-prod2);
00151         //output[0] += prod1*(1.0-prod2);
00152       output.Update(tmp); tmp.Set(0.0);
00153     }
00154   }
00155   Scalar error_indicator(UserVector & input) {
00156     int dimension = (int)input.size();
00157     Scalar norm2  = 0.0;
00158     for (int i=0; i<dimension; i++)
00159       norm2 += input[i]*input[i];
00160     
00161     Scalar ID = AdaptiveSparseGridInterface<Scalar,UserVector>::
00162       getInitialDiff();
00163     norm2 = std::sqrt(norm2)/ID;
00164     return norm2;
00165   }
00166 };
00167 
00168 long double adaptSG(StdVector<long double> & iv,
00169     AdaptiveSparseGridInterface<long double,StdVector<long double> > & 
00170     problem_data, long double TOL) {
00171 
00172   // Construct a Container for the adapted rule
00173   int dimension = problem_data.getDimension();
00174   std::vector<int> index(dimension,1);
00175   
00176   // Initialize global error indicator
00177   long double eta = 1.0;
00178   
00179   // Initialize the Active index set
00180   std::multimap<long double,std::vector<int> > activeIndex;  
00181   activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
00182 
00183   // Initialize the old index set
00184   std::set<std::vector<int> > oldIndex;
00185   // Perform Adaptation
00186   while (eta > TOL) {
00187     eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(
00188             activeIndex,oldIndex,iv,eta,problem_data);
00189   }
00190   return eta;
00191 }
00192 
00193 int main(int argc, char *argv[]) {
00194   
00195   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00196 
00197   // This little trick lets us print to std::cout only if
00198   // a (dummy) command-line argument is provided.
00199   int iprint     = argc - 1;
00200   Teuchos::RCP<std::ostream> outStream;
00201   Teuchos::oblackholestream bhs; // outputs nothing
00202   if (iprint > 0)
00203     outStream = Teuchos::rcp(&std::cout, false);
00204   else
00205     outStream = Teuchos::rcp(&bhs, false);
00206 
00207   // Save the format state of the original std::cout.
00208   Teuchos::oblackholestream oldFormatState;
00209   oldFormatState.copyfmt(std::cout);
00210  
00211   *outStream \
00212   << "===============================================================================\n" \
00213   << "|                                                                             |\n" \
00214   << "|                         Unit Test (AdaptiveSparseGrid)                      |\n" \
00215   << "|                                                                             |\n" \
00216   << "|     1) Particle traveling through 1D slab of length 1.                      |\n" \
00217   << "|                                                                             |\n" \
00218   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00219   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00220   << "|                                                                             |\n" \
00221   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00222   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00223   << "|                                                                             |\n" \
00224   << "===============================================================================\n"\
00225   << "| TEST 17: solve 1D transport problem by approximating an infinite integral   |\n"\
00226   << "===============================================================================\n";
00227 
00228 
00229   // internal variables:
00230   int         errorFlag   = 0;
00231   long double TOL         = INTREPID_TOL;
00232   int dimension           = 8;
00233   std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_PATTERSON);
00234   std::vector<EIntrepidGrowth>   growth1D(dimension,GROWTH_FULLEXP);
00235   int                            maxLevel = 7;
00236   bool                           isNormalized = true;                  
00237   ASGdata<long double,StdVector<long double> > problem_data(dimension,
00238         rule1D,growth1D,maxLevel,isNormalized);
00239   Teuchos::RCP<std::vector<long double> > integralValue = 
00240     Teuchos::rcp(new std::vector<long double>(1,0.0));
00241   StdVector<long double> sol(integralValue); sol.Set(0.0);
00242   problem_data.init(sol);
00243 
00244   long double eta = adaptSG(sol,problem_data,TOL); 
00245   long double x = 0;
00246   long double gamma = 0.5;
00247   long double analyticInt = (1.0 - (1.0-gamma)*exp(gamma*(1.0-x)))/gamma;
00248   long double abstol = std::sqrt(INTREPID_TOL);
00249   long double absdiff = fabs(analyticInt-(*integralValue)[0]);
00250   try { 
00251     *outStream << "Adaptive Sparse Grid exited with global error " 
00252                << std::scientific << std::setprecision(16) << eta << "\n"
00253                << "Approx = " << std::scientific << std::setprecision(16) 
00254                << (*integralValue)[0] 
00255                << ",  Exact = " << std::scientific << std::setprecision(16) 
00256                << analyticInt << "\n"
00257                << "Error = " << std::scientific << std::setprecision(16) 
00258                << absdiff << "   " 
00259                << "<?" << "   " << abstol << "\n"; 
00260     if (absdiff > abstol) {
00261       errorFlag++;
00262       *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00263     }
00264   }
00265   catch (std::logic_error err) {    
00266     *outStream << err.what() << "\n";
00267     errorFlag = -1;
00268   };  
00269 
00270   if (errorFlag != 0)
00271     std::cout << "End Result: TEST FAILED\n";
00272   else
00273     std::cout << "End Result: TEST PASSED\n";
00274 
00275   // reset format state of std::cout
00276   std::cout.copyfmt(oldFormatState);
00277 
00278   return errorFlag;
00279 }