Rythmos_RKButcherTableauHelpers.hpp

00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) 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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 
00030 #ifndef RYTHMOS_RK_BUTCHER_TABLEAU_HELPERS_HPP
00031 #define RYTHMOS_RK_BUTCHER_TABLEAU_HELPERS_HPP
00032 
00033 #include "Rythmos_Types.hpp"
00034 
00035 #include "Rythmos_RKButcherTableauBase.hpp"
00036 #include "Teuchos_Assert.hpp"
00037 #include "Thyra_ProductVectorBase.hpp"
00038 
00039 namespace Rythmos {
00040 
00041 /* \brief . */
00042 template<class Scalar>
00043 void assembleIRKState(
00044   const int stageIndex,
00045   const Teuchos::SerialDenseMatrix<int,Scalar> &A_in,
00046   const Scalar dt,
00047   const Thyra::VectorBase<Scalar> &x_base,
00048   const Thyra::ProductVectorBase<Scalar> &x_stage_bar,
00049   Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
00050   )
00051 {
00052 
00053   typedef ScalarTraits<Scalar> ST;
00054 
00055   const int numStages_in = A_in.numRows();
00056   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE( stageIndex, 0, numStages_in );
00057   TEUCHOS_ASSERT_EQUALITY( A_in.numRows(), numStages_in );
00058   TEUCHOS_ASSERT_EQUALITY( A_in.numCols(), numStages_in );
00059   TEUCHOS_ASSERT_EQUALITY( x_stage_bar.productSpace()->numBlocks(), numStages_in );
00060   Thyra::VectorBase<Scalar>& x_out = *x_out_ptr;
00061 
00062   V_V( outArg(x_out), x_base );
00063   for ( int j = 0; j < numStages_in; ++j ) {
00064     Vp_StV( outArg(x_out), dt * A_in(stageIndex,j), *x_stage_bar.getVectorBlock(j) );
00065   }
00066 
00067 }
00068 
00069 
00070 /* \brief . */
00071 template<class Scalar>
00072 void assembleIRKSolution(
00073   const Teuchos::SerialDenseVector<int,Scalar> &b_in,
00074   const Scalar dt,
00075   const Thyra::VectorBase<Scalar> &x_base,
00076   const Thyra::ProductVectorBase<Scalar> &x_stage_bar,
00077   Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
00078   )
00079 {
00080 
00081   typedef ScalarTraits<Scalar> ST;
00082 
00083   const int numStages_in = b_in.length();
00084   TEUCHOS_ASSERT_EQUALITY( b_in.length(), numStages_in );
00085   TEUCHOS_ASSERT_EQUALITY( x_stage_bar.productSpace()->numBlocks(), numStages_in );
00086   Thyra::VectorBase<Scalar>& x_out = *x_out_ptr;
00087 
00088   V_V( outArg(x_out), x_base );
00089   for ( int j = 0; j < numStages_in; ++j ) {
00090     Vp_StV( outArg(x_out), dt * b_in(j), *x_stage_bar.getVectorBlock(j) );
00091   }
00092 
00093 }
00094 
00095 /* \brief . */
00096 template<class Scalar>
00097 void assembleERKState(
00098   const int stageIndex,
00099   const Teuchos::SerialDenseMatrix<int,Scalar> &A_in,
00100   const Scalar dt,
00101   const Thyra::VectorBase<Scalar> &x_base,
00102   const Thyra::VectorBase<Scalar> &x_stage_bar,
00103   Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
00104   )
00105 {
00106   TEST_FOR_EXCEPT(true);
00107 }
00108 
00109 /* \brief . */
00110 template<class Scalar>
00111 void assembleERKSolution(
00112   const Teuchos::SerialDenseVector<int,Scalar> &b_in,
00113   const Scalar dt,
00114   const Thyra::VectorBase<Scalar> &x_base,
00115   const Thyra::VectorBase<Scalar> &x_stage_bar,
00116   Teuchos::Ptr<Thyra::VectorBase<Scalar> > x_out_ptr
00117   )
00118 {
00119   TEST_FOR_EXCEPT(true);
00120 }
00121 
00122 template<class Scalar>
00123 bool isEmptyRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt ) {
00124   typedef ScalarTraits<Scalar> ST;
00125   
00126   // Check that numStages > 0
00127   if (rkbt.numStages() == 0) {
00128     return true;
00129   }
00130 
00131   // Check that the b vector has _some_ non-zero entry
00132   int numNonZero = 0;
00133   int numStages_local = rkbt.numStages();
00134   const Teuchos::SerialDenseVector<int,Scalar> b_local = rkbt.b();
00135   for (int i=0 ; i<numStages_local ; ++i) {
00136     if (b_local(i) != ST::zero()) {
00137       numNonZero++;
00138     }
00139   }
00140   if (numNonZero == 0) {
00141     return true;
00142   }
00143   // There is no reason to check A and c because they can be zero and you're
00144   // producing an explicit method as long as b has something in it.
00145   return false;
00146 }
00147 
00148 
00149 template<class Scalar>
00150 void assertNonEmptyRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
00151 {
00152   TEST_FOR_EXCEPTION( isEmptyRKButcherTableau(rkbt), std::logic_error,
00153       "Error, this RKButcherTableau is either empty or the b vector is all zeros!\n"
00154       );
00155 }
00156 
00157 template<class Scalar>
00158 bool isDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
00159 {
00160   if (isEmptyRKButcherTableau(rkbt)) {
00161     return false;
00162   }
00163   typedef ScalarTraits<Scalar> ST;
00164   int numStages_local = rkbt.numStages();
00165   const Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
00166   for (int i=0 ; i<numStages_local ; ++i) {
00167     for (int j=0 ; j<numStages_local ; ++j) {
00168       if ((j>i) && (A_local(i,j) != ST::zero())) {
00169         return false;
00170       }
00171     }
00172   }
00173   return true;
00174 }
00175 
00176 template<class Scalar>
00177 bool isIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt ) 
00178 {
00179   if (isEmptyRKButcherTableau(rkbt)) {
00180     return false;
00181   }
00182   return true;
00183 }
00184 
00185 template<class Scalar>
00186 void validateIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
00187 {
00188   TEST_FOR_EXCEPTION( !isIRKButcherTableau(rkbt), std::logic_error,
00189     "Error!  This implicit RK Butcher Tableau is empty!\n"
00190     );
00191 }
00192 
00193 template<class Scalar>
00194 void validateDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
00195 {
00196   TEST_FOR_EXCEPTION( !isDIRKButcherTableau(rkbt), std::logic_error,
00197       "Error!  This Diagonal Implicit RK Butcher Tableau has non-zeros in the upper triangular part!\n" 
00198       );
00199 }
00200 
00201 template<class Scalar>
00202 bool isSDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt ) 
00203 {
00204   if (isEmptyRKButcherTableau(rkbt)) {
00205     return false;
00206   }
00207   if (!isDIRKButcherTableau(rkbt)) {
00208     return false;
00209   }
00210   // Verify the diagonal entries are all equal.
00211   typedef ScalarTraits<Scalar> ST;
00212   int numStages_local = rkbt.numStages();
00213   const Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
00214   Scalar val = A_local(0,0);
00215   for (int i=0 ; i<numStages_local ; ++i) {
00216     if (A_local(i,i) != val) {
00217       return false;
00218     }
00219   }
00220   return true;
00221 }
00222 
00223 template<class Scalar>
00224 void validateSDIRKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
00225 {
00226   TEST_FOR_EXCEPTION( !isSDIRKButcherTableau(rkbt), std::logic_error,
00227       "Error!  This Singly Diagonal Implicit RK Butcher Tableau does not have equal diagonal entries!\n"
00228       );
00229 }
00230 
00231 template<class Scalar>
00232 bool isERKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt) 
00233 {
00234   if (isEmptyRKButcherTableau(rkbt)) {
00235     return false;
00236   }
00237   // Verify the diagonal is zero and the upper triangular part is zero
00238   typedef ScalarTraits<Scalar> ST;
00239   int numStages_local = rkbt.numStages();
00240   const Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
00241   for (int i=0 ; i<numStages_local ; ++i) {
00242     for (int j=0 ; j<numStages_local ; ++j) {
00243       if ((j>=i) && ((A_local(i,j) != ST::zero()))) {
00244         return false;
00245       }
00246     }
00247   }
00248   const Teuchos::SerialDenseVector<int,Scalar> c_local = rkbt.c();
00249   if( c_local(0) != ST::zero() ) {
00250     return false;
00251   }
00252   // 08/13/08 tscoffe:  I'm not sure what else I can check for b & c...
00253   return true;
00254 }
00255 
00256 
00257 
00258 template<class Scalar>
00259 void validateERKButcherTableau( const RKButcherTableauBase<Scalar>& rkbt )
00260 {
00261   TEST_FOR_EXCEPTION( !isERKButcherTableau(rkbt), std::logic_error,
00262       "Error!  This ERK Butcher Tableau is not lower triangular or c(0) is not zero!\n" 
00263       );
00264 }
00265 
00266 /*
00267 template<class Scalar>
00268 void validateERKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
00269 {
00270   typedef ScalarTraits<Scalar> ST;
00271   Teuchos::SerialDenseMatrix<int,Scalar> A_local = rkbt.A();
00272   Teuchos::SerialDenseVector<int,Scalar> b_local = rkbt.b();
00273   Teuchos::SerialDenseVector<int,Scalar> c_local = rkbt.c();
00274   int N = rkbt.numStages();
00275   TEST_FOR_EXCEPT(N == 0);
00276 
00277   if (order_in == 3) {
00278     Scalar sum1 = ST::zero();
00279     Scalar sum2 = ST::zero();
00280     Scalar sum3 = ST::zero();
00281     Scalar sum4 = ST::zero();
00282     for (int j=0 ; j<N ; ++j) {
00283       sum1 += b_local(j);
00284       for (int k=0 ; k<N ; ++k) {
00285         sum2 += 2*b_local(j)*A_local(j,k);
00286         for (int l=0 ; l<N ; ++l) {
00287           sum3 += 3*b_local(j)*A_local(j,k)*A_local(j,l);
00288           sum4 += 6*b_local(j)*A_local(j,k)*A_local(k,l);
00289         }
00290       }
00291     }
00292     TEST_FOR_EXCEPTION(
00293         (
00294          ( sum1 != ST::one() ) || 
00295          ( sum2 != ST::one() ) ||
00296          ( sum3 != ST::one() ) ||
00297          ( sum4 != ST::one() )
00298         ), 
00299         std::logic_error,
00300         "Error!, this RK Butcher Tableau does not meet the order conditions for 3rd order\n"
00301         );
00302   } else {
00303     TEST_FOR_EXCEPTION( true, std::logic_error,
00304         "Error!  this function is only defined for order 3\n"
00305         );
00306   }
00307 }
00308 
00309 template<class Scalar>
00310 void validateIRKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
00311 {
00312   TEST_FOR_EXCEPT(true);
00313 }
00314 
00315 template<class Scalar>
00316 void validateDIRKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
00317 {
00318   TEST_FOR_EXCEPT(true);
00319 }
00320 
00321 template<class Scalar>
00322 void validateSDIRKOrder( RKButcherTableauBase<Scalar> rkbt, int order_in )
00323 {
00324   TEST_FOR_EXCEPT(true);
00325 }
00326 */
00327  
00328 enum E_RKButcherTableauTypes {
00329   RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_INVALID,
00330   RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_ERK,
00331   RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_IRK,
00332   RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK,
00333   RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK
00334 };
00335 
00336 template<class Scalar>
00337 E_RKButcherTableauTypes determineRKBTType(const RKButcherTableauBase<Scalar>& rkbt) {
00338   if (isEmptyRKButcherTableau(rkbt)) {
00339     return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_INVALID;
00340   }
00341   if (isERKButcherTableau(rkbt)) {
00342     return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_ERK;
00343   }
00344   if (rkbt.numStages() == 1) { 
00345     // In this case, its just an IRK method.
00346     return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_IRK;
00347   }
00348   if (isSDIRKButcherTableau(rkbt)) {
00349     return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_SDIRK;
00350   }
00351   if (isDIRKButcherTableau(rkbt)) {
00352     return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_DIRK;
00353   }
00354   if (isIRKButcherTableau(rkbt)) {
00355     return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_IRK;
00356   }
00357   return RYTHMOS_RK_BUTCHER_TABLEAU_TYPE_INVALID;
00358 }
00359 
00360 
00361 
00362 } // namespace Rythmos
00363 
00364 
00365 #endif // RYTHMOS_RK_BUTCHER_TABLEAU_HELPERS_HPP

Generated on Wed May 12 21:25:43 2010 for Rythmos - Transient Integration for Differential Equations by  doxygen 1.4.7