|
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 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 }
1.7.6.1