Sacado Package Browser (Single Doxygen Collection) Version of the Day
fe_jac_fill_funcs.cpp
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 #include "fe_jac_fill_funcs.hpp"
00030 
00031 void analytic_init_fill(const ElemData& e,
00032       unsigned int neqn,
00033       const std::vector<double>& x, 
00034       std::vector<double>& x_local) {
00035   for (unsigned int node=0; node<e.nnode; node++) 
00036     for (unsigned int eqn=0; eqn<neqn; eqn++)
00037       x_local[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
00038 }
00039 
00040 void analytic_element_fill(const ElemData& e, 
00041          unsigned int neqn,
00042          const std::vector<double>& x, 
00043          std::vector<double>& u, 
00044          std::vector<double>& du, 
00045          std::vector<double>& f,
00046          std::vector< std::vector<double> >& jac) {
00047   // Construct element solution, derivative
00048   for (unsigned int qp=0; qp<e.nqp; qp++) {
00049     for (unsigned int eqn=0; eqn<neqn; eqn++) {
00050       u[qp*neqn+eqn] = 0.0;
00051       du[qp*neqn+eqn] = 0.0;
00052       for (unsigned int node=0; node<e.nnode; node++) {
00053   u[qp*neqn+eqn] += x[node*neqn+eqn] * e.phi[qp][node];
00054   du[qp*neqn+eqn] += x[node*neqn+eqn] * e.dphi[qp][node];
00055       }
00056     }
00057   }
00058 
00059   // Compute sum of equations for coupling
00060   std::vector<double> s(e.nqp);
00061   for (unsigned int qp=0; qp<e.nqp; qp++) {
00062     s[qp] = 0.0;
00063     for (unsigned int eqn=0; eqn<neqn; eqn++)
00064       s[qp] += u[qp*neqn+eqn]*u[qp*neqn+eqn];
00065   }
00066 
00067   // Evaluate element residual
00068   for (unsigned int node_row=0; node_row<e.nnode; node_row++) {
00069     for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
00070       unsigned int row = node_row*neqn+eqn_row;
00071       f[row] = 0.0;
00072       for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
00073   for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
00074     unsigned int col = node_col*neqn+eqn_col;
00075     jac[row][col] = 0.0;
00076   }
00077       }
00078       for (unsigned int qp=0; qp<e.nqp; qp++) {
00079   double c1 = e.w[qp]*e.jac[qp];
00080   double c2 = -e.dphi[qp][node_row]/(e.jac[qp]*e.jac[qp]);
00081   double c3 = exp(u[qp*neqn+eqn_row]);
00082   double c4 = e.phi[qp][node_row]*s[qp]*c3;
00083   f[row] += c1*(c2*du[qp*neqn+eqn_row] + c4);
00084   for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
00085     jac[row][node_col*neqn+eqn_row] += 
00086       c1*(c2*e.dphi[qp][node_col] + c4*e.phi[qp][node_col]);
00087     for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++)
00088         jac[row][node_col*neqn+eqn_col] += 
00089       2.0*c1*e.phi[qp][node_row]*u[qp*neqn+eqn_col]*e.phi[qp][node_col]*c3;
00090   }
00091       }
00092     }
00093   }
00094 }
00095 
00096 void analytic_process_fill(const ElemData& e,
00097          unsigned int neqn,
00098          const std::vector<double>& f_local, 
00099          const std::vector< std::vector<double> >& jac_local, 
00100          std::vector<double>& f,
00101          std::vector< std::vector<double> >& jac) {
00102   for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
00103     f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
00104     f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
00105     for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
00106       for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
00107   unsigned int col = node_col*neqn+eqn_col;
00108   unsigned int next_col = (node_col+1)*neqn+eqn_col;
00109   jac[e.gid[0]*neqn+eqn_row][next_col] += jac_local[0*neqn+eqn_row][col];
00110   jac[e.gid[1]*neqn+eqn_row][col] += jac_local[1*neqn+eqn_row][col];
00111       }
00112     }
00113   }
00114 }
00115 
00116 double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00117        double mesh_spacing) {
00118   ElemData e(mesh_spacing);
00119 
00120   // Solution vector, residual, jacobian
00121   std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
00122   std::vector< std::vector<double> > jac(num_nodes*num_eqns);
00123   for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
00124     for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
00125       unsigned int row = node_row*num_eqns + eqn_row;
00126       x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
00127       f[row] = 0.0;
00128       jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
00129       for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
00130   for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
00131     unsigned int col = node_col*num_eqns + eqn_col;
00132     jac[row][col] = 0.0;
00133   }
00134       }
00135     }
00136   }
00137 
00138   Teuchos::Time timer("FE Analytic Jacobian Fill", false);
00139   timer.start(true);
00140   std::vector<double> x_local(e.nnode*num_eqns), f_local(e.nnode*num_eqns);
00141   std::vector< std::vector<double> > jac_local(e.nnode*num_eqns);
00142   std::vector<double> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
00143   for (unsigned int i=0; i<e.nnode*num_eqns; i++)
00144     jac_local[i] = std::vector<double>(e.nnode*num_eqns);
00145   for (unsigned int i=0; i<num_nodes-1; i++) {
00146     e.gid[0] = i;
00147     e.gid[1] = i+1;
00148     
00149     analytic_init_fill(e, num_eqns, x, x_local);
00150     analytic_element_fill(e, num_eqns, x_local, u, du, f_local, jac_local);
00151     analytic_process_fill(e, num_eqns, f_local, jac_local, f, jac);
00152   }
00153   timer.stop();
00154 
00155   // std::cout.precision(8);
00156   // std::cout << "Analytic Residual = " << std::endl;
00157   // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
00158   //   std::cout << "\t" << f[i] << std::endl;
00159 
00160   // std::cout.precision(8);
00161   // std::cout.setf(std::ios::scientific);
00162   // std::cout << "Analytic Jacobian = " << std::endl;
00163   // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
00164   //   std::cout << "\t";
00165   //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
00166   //     std::cout << jac[i][j] << "\t";
00167   //   std::cout << std::endl;
00168   // }
00169 
00170   return timer.totalElapsedTime();
00171 }
00172 
00173 void residual_process_fill(const ElemData& e,
00174          unsigned int neqn,
00175          const std::vector<double>& f_local, 
00176          std::vector<double>& f) {
00177   for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
00178     f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
00179     f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
00180   }
00181 }
00182 
00183 double residual_fill(unsigned int num_nodes, unsigned int num_eqns,
00184          double mesh_spacing) {
00185   ElemData e(mesh_spacing);
00186 
00187   // Solution vector, residual, jacobian
00188   std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
00189   for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
00190     for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
00191       unsigned int row = node_row*num_eqns + eqn_row;
00192       x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
00193       f[row] = 0.0;
00194     }
00195   }
00196 
00197   Teuchos::Time timer("FE Residual Fill", false);
00198   timer.start(true);
00199   std::vector<double> x_local(e.nnode*num_eqns), f_local(e.nnode*num_eqns);
00200   std::vector<double> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
00201   for (unsigned int i=0; i<num_nodes-1; i++) {
00202     e.gid[0] = i;
00203     e.gid[1] = i+1;
00204     
00205     analytic_init_fill(e, num_eqns, x, x_local);
00206     template_element_fill(e, num_eqns, x_local, u, du, f_local);
00207     residual_process_fill(e, num_eqns, f_local, f);
00208   }
00209   timer.stop();
00210 
00211   // std::cout.precision(8);
00212   // std::cout << "Analytic Residual = " << std::endl;
00213   // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
00214   //   std::cout << "\t" << f[i] << std::endl;
00215 
00216   return timer.totalElapsedTime();
00217 }
00218 
00219 #ifdef HAVE_ADOLC
00220 #ifndef ADOLC_TAPELESS
00221 void adolc_init_fill(bool retape,
00222          int tag,
00223          const ElemData& e,
00224          unsigned int neqn,
00225          const std::vector<double>& x, 
00226          std::vector<double>& x_local,
00227          std::vector<adouble>& x_ad) {
00228   if (retape)
00229     trace_on(tag);
00230   for (unsigned int node=0; node<e.nnode; node++)
00231     for (unsigned int eqn=0; eqn<neqn; eqn++) {
00232       x_local[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
00233       if (retape)
00234   x_ad[node*neqn+eqn] <<= x_local[node*neqn+eqn];
00235     }
00236 }
00237 
00238 void adolc_process_fill(bool retape,
00239       int tag, 
00240       const ElemData& e,
00241       unsigned int neqn,
00242       std::vector<double>& x_local,
00243       std::vector<adouble>& f_ad, 
00244       std::vector<double>& f,
00245       std::vector<double>& f_local,
00246       std::vector< std::vector<double> >& jac,
00247       double **seed,
00248       double **jac_local) {
00249   if (retape) {
00250     for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++)
00251       f_ad[0*neqn+eqn_row] >>= f_local[0*neqn+eqn_row];
00252     for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++)
00253       f_ad[1*neqn+eqn_row] >>= f_local[1*neqn+eqn_row];
00254     trace_off();
00255   }
00256 
00257   //jacobian(tag, neqn*e.nnode, neqn*e.nnode, &x_local[0], jac_local);
00258   forward(tag, neqn*e.nnode, neqn*e.nnode, neqn*e.nnode, &x_local[0],
00259     seed, &f_local[0], jac_local);
00260   
00261   for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
00262     f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
00263     f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
00264     for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
00265       for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
00266   unsigned int col = node_col*neqn+eqn_col;
00267   unsigned int next_col = (node_col+1)*neqn+eqn_col;
00268   jac[e.gid[0]*neqn+eqn_row][next_col] += jac_local[0*neqn+eqn_row][col];
00269   jac[e.gid[1]*neqn+eqn_row][col] += jac_local[1*neqn+eqn_row][col];
00270       }
00271     }
00272   }
00273 }
00274 
00275 double adolc_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00276           double mesh_spacing) {
00277   ElemData e(mesh_spacing);
00278 
00279   // Solution vector, residual, jacobian
00280   std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
00281   std::vector< std::vector<double> > jac(num_nodes*num_eqns);
00282   for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
00283     for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
00284       unsigned int row = node_row*num_eqns + eqn_row;
00285       x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
00286       f[row] = 0.0;
00287       jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
00288       for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
00289   for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
00290     unsigned int col = node_col*num_eqns + eqn_col;
00291     jac[row][col] = 0.0;
00292   }
00293       }
00294     }
00295   }
00296 
00297   Teuchos::Time timer("FE ADOL-C Jacobian Fill", false);
00298   timer.start(true);
00299   std::vector<adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
00300   std::vector<adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
00301   std::vector<double> x_local(e.nnode*num_eqns);
00302   std::vector<double> f_local(e.nnode*num_eqns);
00303   double **jac_local = new double*[e.nnode*num_eqns];
00304   double **seed = new double*[e.nnode*num_eqns];
00305   for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
00306     jac_local[i] = new double[e.nnode*num_eqns];
00307     seed[i] = new double[e.nnode*num_eqns];
00308     for (unsigned int j=0; j<e.nnode*num_eqns; j++)
00309       seed[i][j] = 0.0;
00310     seed[i][i] = 1.0;
00311   }
00312   
00313   // Tape first element
00314   e.gid[0] = 0;
00315   e.gid[1] = 1;
00316   adolc_init_fill(true, 0, e, num_eqns, x, x_local, x_ad);
00317   template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
00318   adolc_process_fill(true, 0, e, num_eqns, x_local, f_ad, f, f_local,
00319          jac, seed, jac_local);
00320 
00321   // Now do remaining fills reusing tape
00322   for (unsigned int i=1; i<num_nodes-1; i++) {
00323     e.gid[0] = i;
00324     e.gid[1] = i+1;
00325     
00326     adolc_init_fill(false, 0, e, num_eqns, x, x_local, x_ad);
00327     adolc_process_fill(false, 0, e, num_eqns, x_local, f_ad, f, f_local,
00328            jac, seed, jac_local);
00329   }
00330   for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
00331     delete [] jac_local[i];
00332     delete [] seed[i];
00333   }
00334   delete [] jac_local;
00335   delete [] seed;
00336   timer.stop();
00337 
00338   // std::cout << "ADOL-C Residual = " << std::endl;
00339   // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
00340   //   std::cout << "\t" << f[i] << std::endl;
00341 
00342   // std::cout.precision(8);
00343   // std::cout << "ADOL-C Jacobian = " << std::endl;
00344   // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
00345   //   std::cout << "\t";
00346   //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
00347   //     std::cout << jac[i][j] << "\t";
00348   //   std::cout << std::endl;
00349   // }
00350 
00351   return timer.totalElapsedTime();
00352 }
00353 
00354 double adolc_retape_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00355            double mesh_spacing) {
00356   ElemData e(mesh_spacing);
00357 
00358   // Solution vector, residual, jacobian
00359   std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
00360   std::vector< std::vector<double> > jac(num_nodes*num_eqns);
00361   for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
00362     for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
00363       unsigned int row = node_row*num_eqns + eqn_row;
00364       x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
00365       f[row] = 0.0;
00366       jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
00367       for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
00368   for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
00369     unsigned int col = node_col*num_eqns + eqn_col;
00370     jac[row][col] = 0.0;
00371   }
00372       }
00373     }
00374   }
00375 
00376   Teuchos::Time timer("FE ADOL-C Retape Jacobian Fill", false);
00377   timer.start(true);
00378   std::vector<adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
00379   std::vector<adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
00380   std::vector<double> x_local(e.nnode*num_eqns);
00381   std::vector<double> f_local(e.nnode*num_eqns);
00382   double **jac_local = new double*[e.nnode*num_eqns];
00383   double **seed = new double*[e.nnode*num_eqns];
00384   for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
00385     jac_local[i] = new double[e.nnode*num_eqns];
00386     seed[i] = new double[e.nnode*num_eqns];
00387     for (unsigned int j=0; j<e.nnode*num_eqns; j++)
00388       seed[i][j] = 0.0;
00389     seed[i][i] = 1.0;
00390   }
00391   for (unsigned int i=0; i<num_nodes-1; i++) {
00392     e.gid[0] = i;
00393     e.gid[1] = i+1;
00394     
00395     adolc_init_fill(true, 1, e, num_eqns, x, x_local, x_ad);
00396     template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
00397     adolc_process_fill(true, 1, e, num_eqns, x_local, f_ad, f, f_local,
00398            jac, seed, jac_local);
00399   }
00400   for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
00401     delete [] jac_local[i];
00402     delete [] seed[i];
00403   }
00404   delete [] jac_local;
00405   delete [] seed;
00406   timer.stop();
00407 
00408   // std::cout << "ADOL-C Residual (retaped) = " << std::endl;
00409   // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
00410   //   std::cout << "\t" << f[i] << std::endl;
00411 
00412   // std::cout.precision(8);
00413   // std::cout << "ADOL-C Jacobian (retaped) = " << std::endl;
00414   // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
00415   //   std::cout << "\t";
00416   //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
00417   //     std::cout << jac[i][j] << "\t";
00418   //   std::cout << std::endl;
00419   // }
00420 
00421   return timer.totalElapsedTime();
00422 }
00423 
00424 #else
00425 
00426 void adolc_tapeless_init_fill(const ElemData& e,
00427             unsigned int neqn,
00428             const std::vector<double>& x, 
00429             std::vector<adtl::adouble>& x_ad) {
00430   for (unsigned int node=0; node<e.nnode; node++)
00431     for (unsigned int eqn=0; eqn<neqn; eqn++) {
00432       x_ad[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
00433       for (unsigned int i=0; i<neqn*e.nnode; i++)
00434   x_ad[node*neqn+eqn].setADValue(i, 0.0);
00435       x_ad[node*neqn+eqn].setADValue(node*neqn+eqn, 1.0);
00436     }
00437   
00438 }
00439 
00440 void adolc_tapeless_process_fill(const ElemData& e,
00441          unsigned int neqn,
00442          const std::vector<adtl::adouble>& f_ad, 
00443          std::vector<double>& f,
00444          std::vector< std::vector<double> >& jac) {
00445   for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
00446     f[e.gid[0]*neqn+eqn_row] += f_ad[0*neqn+eqn_row].getValue();
00447     f[e.gid[1]*neqn+eqn_row] += f_ad[1*neqn+eqn_row].getValue();
00448     for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
00449       for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
00450   unsigned int col = node_col*neqn+eqn_col;
00451   unsigned int next_col = (node_col+1)*neqn+eqn_col;
00452   jac[e.gid[0]*neqn+eqn_row][next_col] += 
00453     f_ad[0*neqn+eqn_row].getADValue(col);
00454   jac[e.gid[1]*neqn+eqn_row][col] += 
00455     f_ad[1*neqn+eqn_row].getADValue(col);
00456       }
00457     }
00458   }
00459 }
00460 
00461 double adolc_tapeless_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00462              double mesh_spacing) {
00463   ElemData e(mesh_spacing);
00464 
00465   // Solution vector, residual, jacobian
00466   std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
00467   std::vector< std::vector<double> > jac(num_nodes*num_eqns);
00468   for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
00469     for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
00470       unsigned int row = node_row*num_eqns + eqn_row;
00471       x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
00472       f[row] = 0.0;
00473       jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
00474       for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
00475   for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
00476     unsigned int col = node_col*num_eqns + eqn_col;
00477     jac[row][col] = 0.0;
00478   }
00479       }
00480     }
00481   }
00482 
00483   Teuchos::Time timer("FE Tapeless ADOL-C Jacobian Fill", false);
00484   timer.start(true);
00485   std::vector<adtl::adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
00486   std::vector<adtl::adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
00487   for (unsigned int i=0; i<num_nodes-1; i++) {
00488     e.gid[0] = i;
00489     e.gid[1] = i+1;
00490     
00491     adolc_tapeless_init_fill(e, num_eqns, x, x_ad);
00492     template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
00493     adolc_tapeless_process_fill(e, num_eqns, f_ad, f, jac);
00494   }
00495   timer.stop();
00496 
00497   // std::cout << "Tapeless ADOL-C Residual = " << std::endl;
00498   // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
00499   //   std::cout << "\t" << f[i] << std::endl;
00500 
00501   // std::cout.precision(8);
00502   // std::cout << "Tapeless ADOL-C Jacobian = " << std::endl;
00503   // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
00504   //   std::cout << "\t";
00505   //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
00506   //     std::cout << jac[i][j] << "\t";
00507   //   std::cout << std::endl;
00508   // }
00509 
00510   return timer.totalElapsedTime();
00511 }
00512 #endif
00513 #endif
00514 
00515 #ifdef HAVE_ADIC
00516 void adic_init_fill(const ElemData& e,
00517         unsigned int neqn,
00518         const std::vector<double>& x, 
00519         std::vector<DERIV_TYPE>& x_fad) {
00520   static bool first = true;
00521   for (unsigned int node=0; node<e.nnode; node++)
00522     for (unsigned int eqn=0; eqn<neqn; eqn++) {
00523       x_fad[node*neqn+eqn].value = x[e.gid[node]*neqn+eqn];
00524       if (first)
00525   ad_AD_SetIndep(x_fad[node*neqn+eqn]);
00526     }
00527   if (first) {
00528     ad_AD_SetIndepDone();
00529     first = false;
00530   }
00531 }
00532 
00533 /************************** DISCLAIMER ********************************/
00534 /*                                                                    */
00535 /*   This file was generated on 04/13/12 11:10:49 by the version of   */
00536 /*   ADIC 1.2.3 compiled on  04/14/09 12:39:01                        */
00537 /*                                                                    */
00538 /*   ADIC was prepared as an account of work sponsored by an          */
00539 /*   agency of the United States Government and the University of     */
00540 /*   Chicago.  NEITHER THE AUTHOR(S), THE UNITED STATES GOVERNMENT    */
00541 /*   NOR ANY AGENCY THEREOF, NOR THE UNIVERSITY OF CHICAGO, INCLUDING */
00542 /*   ANY OF THEIR EMPLOYEES OR OFFICERS, MAKES ANY WARRANTY, EXPRESS  */
00543 /*   OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR */
00544 /*   THE ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY INFORMATION OR  */
00545 /*   PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE */
00546 /*   PRIVATELY OWNED RIGHTS.                                          */
00547 /*                                                                    */
00548 /**********************************************************************/
00549 void   adic_element_fill(ElemData  *e,unsigned int  neqn, DERIV_TYPE  *x,DERIV_TYPE  *u,DERIV_TYPE  *du,DERIV_TYPE  *f) {
00550 unsigned int  var_0, var_1, var_2, var_3, var_4, var_5, var_6, var_7;
00551 DERIV_TYPE  var_8;
00552 double  adji_0;
00553     double  loc_0;
00554     double  loc_1;
00555     double  loc_2;
00556     double  loc_3;
00557     double  loc_4;
00558     double  loc_5;
00559     double  loc_6;
00560     double  loc_7;
00561     double  loc_8;
00562     double  loc_9;
00563     double  adj_0;
00564     double  adj_1;
00565     double  adj_2;
00566     double  adj_3;
00567 
00568         // static int g_filenum = 0;
00569         // if (g_filenum == 0) {
00570         //     adintr_ehsfid(&g_filenum, __FILE__, "adic_element_fill");
00571         // }
00572             for (unsigned int  qp = 0;     qp < e->nqp;     )    {
00573         for (unsigned int  eqn = 0;         eqn < neqn;         )        {
00574             {
00575                 ad_grad_axpy_0(&(u[qp * neqn + eqn]));
00576                 DERIV_val(u[qp * neqn + eqn]) = 0.0;
00577             }
00578             {
00579                 ad_grad_axpy_0(&(du[qp * neqn + eqn]));
00580                 DERIV_val(du[qp * neqn + eqn]) = 0.0;
00581             }
00582             for (unsigned int  node = 0;             node < e->nnode;             )            {
00583                 {
00584                     loc_0 = DERIV_val(x[node * neqn + eqn]) * e->phi[qp][node];
00585                     loc_1 = DERIV_val(u[qp * neqn + eqn]) + loc_0;
00586                     ad_grad_axpy_2(&(u[qp * neqn + eqn]), 1.000000000000000e+00, &(u[qp * neqn + eqn]), e->phi[qp][node], &(x[node * neqn + eqn]));
00587                     DERIV_val(u[qp * neqn + eqn]) = loc_1;
00588                 }
00589                 {
00590                     loc_0 = DERIV_val(x[node * neqn + eqn]) * e->dphi[qp][node];
00591                     loc_1 = DERIV_val(du[qp * neqn + eqn]) + loc_0;
00592                     ad_grad_axpy_2(&(du[qp * neqn + eqn]), 1.000000000000000e+00, &(du[qp * neqn + eqn]), e->dphi[qp][node], &(x[node * neqn + eqn]));
00593                     DERIV_val(du[qp * neqn + eqn]) = loc_1;
00594                 }
00595                 var_2 = node++;
00596             }
00597             var_1 = eqn++;
00598         }
00599         var_0 = qp++;
00600     }
00601 //DERIV_TYPE  *s = malloc(e->nqp *  sizeof (DERIV_TYPE ));
00602       std::vector<DERIV_TYPE> s(e->nqp);
00603     for (unsigned int  qp = 0;     qp < e->nqp;     )    {
00604         {
00605             ad_grad_axpy_0(&(s[qp]));
00606             DERIV_val(s[qp]) = 0.0;
00607         }
00608         for (unsigned int  eqn = 0;         eqn < neqn;         )        {
00609             {
00610                 loc_0 = DERIV_val(u[qp * neqn + eqn]) * DERIV_val(u[qp * neqn + eqn]);
00611                 loc_1 = DERIV_val(s[qp]) + loc_0;
00612                 ad_grad_axpy_3(&(s[qp]), 1.000000000000000e+00, &(s[qp]), DERIV_val(u[qp * neqn + eqn]), &(u[qp * neqn + eqn]), DERIV_val(u[qp * neqn + eqn]), &(u[qp * neqn + eqn]));
00613                 DERIV_val(s[qp]) = loc_1;
00614             }
00615             var_4 = eqn++;
00616         }
00617         var_3 = qp++;
00618     }
00619     for (unsigned int  node = 0;     node < e->nnode;     )    {
00620         for (unsigned int  eqn = 0;         eqn < neqn;         )        {
00621 unsigned int  row = node * neqn + eqn;
00622             {
00623                 ad_grad_axpy_0(&(f[row]));
00624                 DERIV_val(f[row]) = 0.0;
00625             }
00626             for (unsigned int  qp = 0;             qp < e->nqp;             )            {
00627      DERIV_val(var_8) = exp(( DERIV_val(u[qp * neqn + eqn])));
00628       adji_0 = DERIV_val(var_8);
00629                 {
00630                     ad_grad_axpy_1(&(var_8), adji_0, &(u[qp * neqn + eqn]));
00631                 }
00632                 {
00633                     loc_0 = e->w[qp] * e->jac[qp];
00634                     loc_1 =  -e->dphi[qp][node];
00635                     loc_2 = e->jac[qp] * e->jac[qp];
00636                     loc_3 = loc_1 / loc_2;
00637                     loc_4 = loc_3 * DERIV_val(du[qp * neqn + eqn]);
00638                     loc_5 = e->phi[qp][node] * DERIV_val(s[qp]);
00639                     loc_6 = loc_5 * DERIV_val(var_8);
00640                     loc_7 = loc_4 + loc_6;
00641                     loc_8 = loc_0 * loc_7;
00642                     loc_9 = DERIV_val(f[row]) + loc_8;
00643                     adj_0 = loc_5 * loc_0;
00644                     adj_1 = DERIV_val(var_8) * loc_0;
00645                     adj_2 = e->phi[qp][node] * adj_1;
00646                     adj_3 = loc_3 * loc_0;
00647                     ad_grad_axpy_4(&(f[row]), 1.000000000000000e+00, &(f[row]), adj_3, &(du[qp * neqn + eqn]), adj_2, &(s[qp]), adj_0, &(var_8));
00648                     DERIV_val(f[row]) = loc_9;
00649                 }
00650                 var_7 = qp++;
00651             }
00652             var_6 = eqn++;
00653         }
00654         var_5 = node++;
00655     }
00656     //free(s);
00657 }
00658 
00659 void adic_process_fill(const ElemData& e,
00660            unsigned int neqn,
00661            const std::vector<DERIV_TYPE>& f_fad, 
00662            std::vector<double>& f,
00663            std::vector< std::vector<double> >& jac) {
00664   for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
00665     f[e.gid[0]*neqn+eqn_row] += f_fad[0*neqn+eqn_row].value;
00666     f[e.gid[1]*neqn+eqn_row] += f_fad[1*neqn+eqn_row].value;
00667     for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
00668       for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
00669   unsigned int col = node_col*neqn+eqn_col;
00670   unsigned int next_col = (node_col+1)*neqn+eqn_col;
00671   jac[e.gid[0]*neqn+eqn_row][next_col] += 
00672     f_fad[0*neqn+eqn_row].grad[col];
00673   jac[e.gid[1]*neqn+eqn_row][col] += 
00674     f_fad[1*neqn+eqn_row].grad[col];
00675       }
00676     }
00677   }
00678 }
00679 
00680 double adic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
00681        double mesh_spacing) {
00682   AD_Init(0);
00683   ElemData e(mesh_spacing);
00684 
00685   // Solution vector, residual, jacobian
00686   std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
00687   std::vector< std::vector<double> > jac(num_nodes*num_eqns);
00688   for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
00689     for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) { 
00690       unsigned int row = node_row*num_eqns + eqn_row;
00691       x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
00692       f[row] = 0.0;
00693       jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
00694       for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
00695   for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) { 
00696     unsigned int col = node_col*num_eqns + eqn_col;
00697     jac[row][col] = 0.0;
00698   }
00699       }
00700     }
00701   }
00702 
00703   Teuchos::Time timer("FE ADIC Jacobian Fill", false);
00704   timer.start(true);
00705   std::vector<DERIV_TYPE> x_fad(e.nnode*num_eqns), f_fad(e.nnode*num_eqns);
00706   std::vector<DERIV_TYPE> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
00707   for (unsigned int i=0; i<num_nodes-1; i++) {
00708     e.gid[0] = i;
00709     e.gid[1] = i+1;
00710     
00711     adic_init_fill(e, num_eqns, x, x_fad);
00712     adic_element_fill(&e, num_eqns, &x_fad[0], &u[0], &du[0], &f_fad[0]);
00713     adic_process_fill(e, num_eqns, f_fad, f, jac);
00714   }
00715   timer.stop();
00716   AD_Final();
00717 
00718   // std::cout.precision(8);
00719   // std::cout << "ADIC Residual = " << std::endl;
00720   // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
00721   //   std::cout << "\t" << f[i] << std::endl;
00722 
00723   // std::cout.precision(8);
00724   // std::cout << "ADIC Jacobian = " << std::endl;
00725   // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
00726   //   std::cout << "\t";
00727   //   for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
00728   //     std::cout << jac[i][j] << "\t";
00729   //   std::cout << std::endl;
00730   // }
00731 
00732   return timer.totalElapsedTime();
00733 }
00734 #endif
00735 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines