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