Sacado Package Browser (Single Doxygen Collection) Version of the Day
fe_jac_fill_funcs.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                           Sacado Package
00005 //                 Copyright (2006) Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
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 David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00025 // (etphipp@sandia.gov).
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 
00030 #ifndef FE_JAC_FILL_FUNCS_HPP
00031 #define FE_JAC_FILL_FUNCS_HPP
00032 
00033 #include "Sacado.hpp"
00034 #include "Sacado_Fad_SimpleFad.hpp"
00035 #include "Sacado_CacheFad_DFad.hpp"
00036 #include "Sacado_CacheFad_SFad.hpp"
00037 #include "Sacado_CacheFad_SLFad.hpp"
00038 
00039 
00040 
00041 #include "Teuchos_Time.hpp"
00042 #include "Teuchos_CommandLineProcessor.hpp"
00043 
00044 // ADOL-C includes
00045 #ifdef HAVE_ADOLC
00046 #ifdef PACKAGE
00047 #undef PACKAGE
00048 #endif
00049 #ifdef PACKAGE_NAME
00050 #undef PACKAGE_NAME
00051 #endif
00052 #ifdef PACKAGE_BUGREPORT
00053 #undef PACKAGE_BUGREPORT
00054 #endif
00055 #ifdef PACKAGE_STRING
00056 #undef PACKAGE_STRING
00057 #endif
00058 #ifdef PACKAGE_TARNAME
00059 #undef PACKAGE_TARNAME
00060 #endif
00061 #ifdef PACKAGE_VERSION
00062 #undef PACKAGE_VERSION
00063 #endif
00064 #ifdef VERSION
00065 #undef VERSION
00066 #endif
00067 //#define ADOLC_TAPELESS
00068 #define NUMBER_DIRECTIONS 100
00069 #include "adolc/adouble.h"
00070 #include "adolc/drivers/drivers.h"
00071 #include "adolc/interfaces.h"
00072 #endif
00073 
00074 #ifdef HAVE_ADIC
00075 // We have included an ADIC differentiated version of the element fill 
00076 // routine to compare the speed of operator overloading to source
00077 // transformation.  To run this code, all that is necessary is to turn
00078 // on the ADIC TPL.  However to modify the code, it is necessary to
00079 // re-run the ADIC source transformation tool.  To do so, first update
00080 // the changes to adic_element_fill.c.  Then set the following environment
00081 // variables:
00082 //     export ADIC_ARCH=linux
00083 //     export ADIC=/home/etphipp/AD_libs/adic
00084 // Next run ADIC via in the tests/performance source directory:
00085 //     ${ADIC}/bin/linux/adiC -vd gradient -i adic_element_fill.init
00086 // Finally, copy the resulting differentiated function in adic_element_fill.ad.c
00087 // back into this file.  The function will need to be edited by changing
00088 // the allocation of s to a std::vector<DERIV_TYPE> (the compiler 
00089 // doesn't seem to like malloc), and commenting out the g_filenum lines.
00090 #define ad_GRAD_MAX 130
00091 #include "ad_deriv.h"
00092 #endif
00093 
00094 
00095 // A performance test that computes a finite-element-like Jacobian using
00096 // several Fad variants
00097 
00098 struct ElemData {
00099   static const unsigned int nqp = 2;
00100   static const unsigned int nnode = 2;
00101   double w[nqp], jac[nqp], phi[nqp][nnode], dphi[nqp][nnode];
00102   unsigned int gid[nnode];
00103 
00104   ElemData(double mesh_spacing) {
00105     // Quadrature points
00106     double xi[nqp];
00107     xi[0] = -1.0 / std::sqrt(3.0);
00108     xi[1] =  1.0 / std::sqrt(3.0);
00109     
00110     for (unsigned int i=0; i<nqp; i++) {
00111       // Weights
00112       w[i] = 1.0;
00113       
00114       // Element transformation Jacobian
00115       jac[i] = 0.5*mesh_spacing;
00116       
00117       // Shape functions & derivatives
00118       phi[i][0] = 0.5*(1.0 - xi[i]);
00119       phi[i][1] = 0.5*(1.0 + xi[i]);
00120       dphi[i][0] = -0.5;
00121       dphi[i][1] = 0.5;
00122     }
00123   }
00124 };
00125 
00126 template <class FadType>
00127 void fad_init_fill(const ElemData& e,
00128        unsigned int neqn,
00129        const std::vector<double>& x, 
00130        std::vector<FadType>& x_fad) {
00131   for (unsigned int node=0; node<e.nnode; node++)
00132     for (unsigned int eqn=0; eqn<neqn; eqn++)
00133       x_fad[node*neqn+eqn] = FadType(e.nnode*neqn, node*neqn+eqn, 
00134              x[e.gid[node]*neqn+eqn]);
00135 }
00136 
00137 template <class T>
00138 void template_element_fill(const ElemData& e, 
00139          unsigned int neqn,
00140          const std::vector<T>& x, 
00141          std::vector<T>& u, 
00142          std::vector<T>& du, 
00143          std::vector<T>& f) {
00144   // Construct element solution, derivative
00145   for (unsigned int qp=0; qp<e.nqp; qp++) {
00146     for (unsigned int eqn=0; eqn<neqn; eqn++) {
00147       u[qp*neqn+eqn] = 0.0;
00148       du[qp*neqn+eqn] = 0.0;
00149       for (unsigned int node=0; node<e.nnode; node++) {
00150   u[qp*neqn+eqn] += x[node*neqn+eqn] * e.phi[qp][node];
00151   du[qp*neqn+eqn] += x[node*neqn+eqn] * e.dphi[qp][node];
00152       }
00153     }
00154   }
00155 
00156   // Compute sum of equations for coupling
00157   std::vector<T> s(e.nqp);
00158   for (unsigned int qp=0; qp<e.nqp; qp++) {
00159     s[qp] = 0.0;
00160     for (unsigned int eqn=0; eqn<neqn; eqn++)
00161       s[qp] += u[qp*neqn+eqn]*u[qp*neqn+eqn];
00162   }
00163 
00164   // Evaluate element residual
00165   for (unsigned int node=0; node<e.nnode; node++) {
00166     for (unsigned int eqn=0; eqn<neqn; eqn++) {
00167       unsigned int row = node*neqn+eqn;
00168       f[row] = 0.0;
00169       for (unsigned int qp=0; qp<e.nqp; qp++) {
00170   double c1 = e.w[qp]*e.jac[qp];
00171   double c2 = -e.dphi[qp][node]/(e.jac[qp]*e.jac[qp]);
00172   f[row] += 
00173     c1*(c2*du[qp*neqn+eqn] + e.phi[qp][node]*s[qp]*exp(u[qp*neqn+eqn]));
00174       }
00175     }
00176   }
00177 }
00178 
00179 template <class FadType>
00180 void fad_process_fill(const ElemData& e,
00181           unsigned int neqn,
00182           const std::vector<FadType>& f_fad, 
00183           std::vector<double>& f,
00184           std::vector< std::vector<double> >& jac) {
00185   for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
00186     f[e.gid[0]*neqn+eqn_row] += f_fad[0*neqn+eqn_row].val();
00187     f[e.gid[1]*neqn+eqn_row] += f_fad[1*neqn+eqn_row].val();
00188     for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
00189       for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
00190   unsigned int col = node_col*neqn+eqn_col;
00191   unsigned int next_col = (node_col+1)*neqn+eqn_col;
00192   jac[e.gid[0]*neqn+eqn_row][next_col] += 
00193     f_fad[0*neqn+eqn_row].fastAccessDx(col);
00194   jac[e.gid[1]*neqn+eqn_row][col] += 
00195     f_fad[1*neqn+eqn_row].fastAccessDx(col);
00196       }
00197     }
00198   }
00199 }
00200 
00201 template <class FadType>
00202 double fad_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00203         double mesh_spacing) {
00204   ElemData e(mesh_spacing);
00205 
00206   // Solution vector, residual, jacobian
00207   std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
00208   std::vector< std::vector<double> > jac(num_nodes*num_eqns);
00209   for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
00210     for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
00211       unsigned int row = node_row*num_eqns + eqn_row;
00212       x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
00213       f[row] = 0.0;
00214       jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
00215       for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
00216   for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
00217     unsigned int col = node_col*num_eqns + eqn_col;
00218     jac[row][col] = 0.0;
00219   }
00220       }
00221     }
00222   }
00223 
00224   Teuchos::Time timer("FE Fad Jacobian Fill", false);
00225   timer.start(true);
00226   std::vector<FadType> x_fad(e.nnode*num_eqns), f_fad(e.nnode*num_eqns);
00227   std::vector<FadType> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
00228   for (unsigned int i=0; i<num_nodes-1; i++) {
00229     e.gid[0] = i;
00230     e.gid[1] = i+1;
00231     
00232     fad_init_fill(e, num_eqns, x, x_fad);
00233     template_element_fill(e, num_eqns, x_fad, u, du, f_fad);
00234     fad_process_fill(e, num_eqns, f_fad, f, jac);
00235   }
00236   timer.stop();
00237 
00238   // std::cout << "Fad Residual = " << std::endl;
00239   // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
00240   //   std::cout << "\t" << f[i] << std::endl;
00241 
00242   // std::cout.precision(8);
00243   // std::cout << "Fad Jacobian = " << std::endl;
00244   // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
00245   //   std::cout << "\t";
00246   //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
00247   //     std::cout << jac[i][j] << "\t";
00248   //   std::cout << std::endl;
00249   // }
00250 
00251   return timer.totalElapsedTime();
00252 }
00253 
00254 double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00255        double mesh_spacing);
00256 double residual_fill(unsigned int num_nodes, unsigned int num_eqns,
00257          double mesh_spacing);
00258 
00259 #ifdef HAVE_ADOLC
00260 #ifndef ADOLC_TAPELESS
00261 void adolc_init_fill(bool retape,
00262          int tag,
00263          const ElemData& e,
00264          unsigned int neqn,
00265          const std::vector<double>& x, 
00266          std::vector<double>& x_local,
00267          std::vector<adouble>& x_ad);
00268 
00269 void adolc_process_fill(bool retape,
00270       int tag, 
00271       const ElemData& e,
00272       unsigned int neqn,
00273       std::vector<double>& x_local,
00274       std::vector<adouble>& f_ad, 
00275       std::vector<double>& f,
00276       std::vector<double>& f_local,
00277       std::vector< std::vector<double> >& jac,
00278       double **seed,
00279       double **jac_local);
00280 
00281 double adolc_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00282           double mesh_spacing);
00283 
00284 double adolc_retape_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00285            double mesh_spacing);
00286 
00287 #else
00288 
00289 void adolc_tapeless_init_fill(const ElemData& e,
00290             unsigned int neqn,
00291             const std::vector<double>& x, 
00292             std::vector<adtl::adouble>& x_ad);
00293 
00294 void adolc_tapeless_process_fill(const ElemData& e,
00295          unsigned int neqn,
00296          const std::vector<adtl::adouble>& f_ad, 
00297          std::vector<double>& f,
00298          std::vector< std::vector<double> >& jac);
00299 
00300 double adolc_tapeless_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00301              double mesh_spacing);
00302 
00303 #endif
00304 #endif
00305 
00306 #ifdef HAVE_ADIC
00307 void adic_init_fill(const ElemData& e,
00308         unsigned int neqn,
00309         const std::vector<double>& x, 
00310         std::vector<DERIV_TYPE>& x_fad);
00311 void   adic_element_fill(ElemData  *e,unsigned int  neqn,DERIV_TYPE  *x,DERIV_TYPE  *u,DERIV_TYPE  *du,DERIV_TYPE  *f);
00312 void adic_process_fill(const ElemData& e,
00313            unsigned int neqn,
00314            const std::vector<DERIV_TYPE>& f_fad, 
00315            std::vector<double>& f,
00316            std::vector< std::vector<double> >& jac);
00317 double adic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00318          double mesh_spacing);
00319 inline void   AD_Init(int  arg0) {
00320   ad_AD_GradInit(arg0);  
00321 }
00322 inline void   AD_Final() {
00323   ad_AD_GradFinal();
00324 }
00325 #endif
00326 
00327 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines