Sacado Package Browser (Single Doxygen Collection)  Version of the Day
fe_jac_fill_funcs.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Sacado Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25 // (etphipp@sandia.gov).
26 //
27 // ***********************************************************************
28 // @HEADER
29 
30 #ifndef FE_JAC_FILL_FUNCS_HPP
31 #define FE_JAC_FILL_FUNCS_HPP
32 
33 #include "Sacado_No_Kokkos.hpp"
34 #include "Sacado_Fad_SimpleFad.hpp"
35 #include "Sacado_CacheFad_DFad.hpp"
36 #include "Sacado_CacheFad_SFad.hpp"
38 
39 
40 
41 #include "Teuchos_Time.hpp"
43 
44 // ADOL-C includes
45 #ifdef HAVE_ADOLC
46 #ifdef PACKAGE
47 #undef PACKAGE
48 #endif
49 #ifdef PACKAGE_NAME
50 #undef PACKAGE_NAME
51 #endif
52 #ifdef PACKAGE_BUGREPORT
53 #undef PACKAGE_BUGREPORT
54 #endif
55 #ifdef PACKAGE_STRING
56 #undef PACKAGE_STRING
57 #endif
58 #ifdef PACKAGE_TARNAME
59 #undef PACKAGE_TARNAME
60 #endif
61 #ifdef PACKAGE_VERSION
62 #undef PACKAGE_VERSION
63 #endif
64 #ifdef VERSION
65 #undef VERSION
66 #endif
67 //#define ADOLC_TAPELESS
68 #define NUMBER_DIRECTIONS 100
69 #include "adolc/adouble.h"
70 #include "adolc/drivers/drivers.h"
71 #include "adolc/interfaces.h"
72 #endif
73 
74 #ifdef HAVE_ADIC
75 // We have included an ADIC differentiated version of the element fill
76 // routine to compare the speed of operator overloading to source
77 // transformation. To run this code, all that is necessary is to turn
78 // on the ADIC TPL. However to modify the code, it is necessary to
79 // re-run the ADIC source transformation tool. To do so, first update
80 // the changes to adic_element_fill.c. Then set the following environment
81 // variables:
82 // export ADIC_ARCH=linux
83 // export ADIC=/home/etphipp/AD_libs/adic
84 // Next run ADIC via in the tests/performance source directory:
85 // ${ADIC}/bin/linux/adiC -vd gradient -i adic_element_fill.init
86 // Finally, copy the resulting differentiated function in adic_element_fill.ad.c
87 // back into this file. The function will need to be edited by changing
88 // the allocation of s to a std::vector<DERIV_TYPE> (the compiler
89 // doesn't seem to like malloc), and commenting out the g_filenum lines.
90 #define ad_GRAD_MAX 130
91 #include "ad_deriv.h"
92 #endif
93 
94 
95 // A performance test that computes a finite-element-like Jacobian using
96 // several Fad variants
97 
98 struct ElemData {
99  static const unsigned int nqp = 2;
100  static const unsigned int nnode = 2;
101  double w[nqp], jac[nqp], phi[nqp][nnode], dphi[nqp][nnode];
102  unsigned int gid[nnode];
103 
104  ElemData(double mesh_spacing) {
105  // Quadrature points
106  double xi[nqp];
107  xi[0] = -1.0 / std::sqrt(3.0);
108  xi[1] = 1.0 / std::sqrt(3.0);
109 
110  for (unsigned int i=0; i<nqp; i++) {
111  // Weights
112  w[i] = 1.0;
113 
114  // Element transformation Jacobian
115  jac[i] = 0.5*mesh_spacing;
116 
117  // Shape functions & derivatives
118  phi[i][0] = 0.5*(1.0 - xi[i]);
119  phi[i][1] = 0.5*(1.0 + xi[i]);
120  dphi[i][0] = -0.5;
121  dphi[i][1] = 0.5;
122  }
123  }
124 };
125 
126 template <class FadType>
127 void fad_init_fill(const ElemData& e,
128  unsigned int neqn,
129  const std::vector<double>& x,
130  std::vector<FadType>& x_fad) {
131  for (unsigned int node=0; node<e.nnode; node++)
132  for (unsigned int eqn=0; eqn<neqn; eqn++)
133  x_fad[node*neqn+eqn] = FadType(e.nnode*neqn, node*neqn+eqn,
134  x[e.gid[node]*neqn+eqn]);
135 }
136 
137 template <class T>
139  unsigned int neqn,
140  const std::vector<T>& x,
141  std::vector<T>& u,
142  std::vector<T>& du,
143  std::vector<T>& f) {
144  // Construct element solution, derivative
145  for (unsigned int qp=0; qp<e.nqp; qp++) {
146  for (unsigned int eqn=0; eqn<neqn; eqn++) {
147  u[qp*neqn+eqn] = 0.0;
148  du[qp*neqn+eqn] = 0.0;
149  for (unsigned int node=0; node<e.nnode; node++) {
150  u[qp*neqn+eqn] += x[node*neqn+eqn] * e.phi[qp][node];
151  du[qp*neqn+eqn] += x[node*neqn+eqn] * e.dphi[qp][node];
152  }
153  }
154  }
155 
156  // Compute sum of equations for coupling
157  std::vector<T> s(e.nqp);
158  for (unsigned int qp=0; qp<e.nqp; qp++) {
159  s[qp] = 0.0;
160  for (unsigned int eqn=0; eqn<neqn; eqn++)
161  s[qp] += u[qp*neqn+eqn]*u[qp*neqn+eqn];
162  }
163 
164  // Evaluate element residual
165  for (unsigned int node=0; node<e.nnode; node++) {
166  for (unsigned int eqn=0; eqn<neqn; eqn++) {
167  unsigned int row = node*neqn+eqn;
168  f[row] = 0.0;
169  for (unsigned int qp=0; qp<e.nqp; qp++) {
170  double c1 = e.w[qp]*e.jac[qp];
171  double c2 = -e.dphi[qp][node]/(e.jac[qp]*e.jac[qp]);
172  f[row] +=
173  c1*(c2*du[qp*neqn+eqn] + e.phi[qp][node]*s[qp]*exp(u[qp*neqn+eqn]));
174  }
175  }
176  }
177 }
178 
179 template <class FadType>
181  unsigned int neqn,
182  const std::vector<FadType>& f_fad,
183  std::vector<double>& f,
184  std::vector< std::vector<double> >& jac) {
185  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
186  f[e.gid[0]*neqn+eqn_row] += f_fad[0*neqn+eqn_row].val();
187  f[e.gid[1]*neqn+eqn_row] += f_fad[1*neqn+eqn_row].val();
188  for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
189  for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
190  unsigned int col = node_col*neqn+eqn_col;
191  unsigned int next_col = (node_col+1)*neqn+eqn_col;
192  jac[e.gid[0]*neqn+eqn_row][next_col] +=
193  f_fad[0*neqn+eqn_row].fastAccessDx(col);
194  jac[e.gid[1]*neqn+eqn_row][col] +=
195  f_fad[1*neqn+eqn_row].fastAccessDx(col);
196  }
197  }
198  }
199 }
200 
201 template <class FadType>
202 double fad_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
203  double mesh_spacing) {
204  ElemData e(mesh_spacing);
205 
206  // Solution vector, residual, jacobian
207  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
208  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
209  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
210  for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
211  unsigned int row = node_row*num_eqns + eqn_row;
212  x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
213  f[row] = 0.0;
214  jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
215  for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
216  for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
217  unsigned int col = node_col*num_eqns + eqn_col;
218  jac[row][col] = 0.0;
219  }
220  }
221  }
222  }
223 
224  Teuchos::Time timer("FE Fad Jacobian Fill", false);
225  timer.start(true);
226  std::vector<FadType> x_fad(e.nnode*num_eqns), f_fad(e.nnode*num_eqns);
227  std::vector<FadType> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
228  for (unsigned int i=0; i<num_nodes-1; i++) {
229  e.gid[0] = i;
230  e.gid[1] = i+1;
231 
232  fad_init_fill(e, num_eqns, x, x_fad);
233  template_element_fill(e, num_eqns, x_fad, u, du, f_fad);
234  fad_process_fill(e, num_eqns, f_fad, f, jac);
235  }
236  timer.stop();
237 
238  // std::cout << "Fad Residual = " << std::endl;
239  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
240  // std::cout << "\t" << f[i] << std::endl;
241 
242  // std::cout.precision(8);
243  // std::cout << "Fad Jacobian = " << std::endl;
244  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
245  // std::cout << "\t";
246  // for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
247  // std::cout << jac[i][j] << "\t";
248  // std::cout << std::endl;
249  // }
250 
251  return timer.totalElapsedTime();
252 }
253 
254 double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
255  double mesh_spacing);
256 double residual_fill(unsigned int num_nodes, unsigned int num_eqns,
257  double mesh_spacing);
258 
259 #ifdef HAVE_ADOLC
260 #ifndef ADOLC_TAPELESS
261 void adolc_init_fill(bool retape,
262  int tag,
263  const ElemData& e,
264  unsigned int neqn,
265  const std::vector<double>& x,
266  std::vector<double>& x_local,
267  std::vector<adouble>& x_ad);
268 
269 void adolc_process_fill(bool retape,
270  int tag,
271  const ElemData& e,
272  unsigned int neqn,
273  std::vector<double>& x_local,
274  std::vector<adouble>& f_ad,
275  std::vector<double>& f,
276  std::vector<double>& f_local,
277  std::vector< std::vector<double> >& jac,
278  double **seed,
279  double **jac_local);
280 
281 double adolc_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
282  double mesh_spacing);
283 
284 double adolc_retape_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
285  double mesh_spacing);
286 
287 #else
288 
289 void adolc_tapeless_init_fill(const ElemData& e,
290  unsigned int neqn,
291  const std::vector<double>& x,
292  std::vector<adtl::adouble>& x_ad);
293 
294 void adolc_tapeless_process_fill(const ElemData& e,
295  unsigned int neqn,
296  const std::vector<adtl::adouble>& f_ad,
297  std::vector<double>& f,
298  std::vector< std::vector<double> >& jac);
299 
300 double adolc_tapeless_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
301  double mesh_spacing);
302 
303 #endif
304 #endif
305 
306 #ifdef HAVE_ADIC
307 void adic_init_fill(const ElemData& e,
308  unsigned int neqn,
309  const std::vector<double>& x,
310  std::vector<DERIV_TYPE>& x_fad);
311 void adic_element_fill(ElemData *e,unsigned int neqn,DERIV_TYPE *x,DERIV_TYPE *u,DERIV_TYPE *du,DERIV_TYPE *f);
312 void adic_process_fill(const ElemData& e,
313  unsigned int neqn,
314  const std::vector<DERIV_TYPE>& f_fad,
315  std::vector<double>& f,
316  std::vector< std::vector<double> >& jac);
317 double adic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
318  double mesh_spacing);
319 inline void AD_Init(int arg0) {
320  ad_AD_GradInit(arg0);
321 }
322 inline void AD_Final() {
323  ad_AD_GradFinal();
324 }
325 #endif
326 
327 #endif
void f()
double residual_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
void AD_Init(int)
void template_element_fill(const ElemData &e, unsigned int neqn, const std::vector< T > &x, std::vector< T > &u, std::vector< T > &du, std::vector< T > &f)
InactiveDouble * jac
InactiveDouble * w
void AD_Final()
void adic_element_fill(ElemData *e, unsigned int neqn, const DERIV_TYPE *x, DERIV_TYPE *u, DERIV_TYPE *du, DERIV_TYPE *f)
Sacado::Fad::DFad< double > FadType
SimpleFad< ValueT > sqrt(const SimpleFad< ValueT > &a)
ElemData(double mesh_spacing)
void fad_init_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &x, std::vector< FadType > &x_fad)
SimpleFad< ValueT > exp(const SimpleFad< ValueT > &a)
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
void start(bool reset=false)
double stop()
double totalElapsedTime(bool readCurrentTime=false) const
InactiveDouble ** phi
InactiveDouble ** dphi
void fad_process_fill(const ElemData &e, unsigned int neqn, const std::vector< FadType > &f_fad, std::vector< double > &f, std::vector< std::vector< double > > &jac)
double fad_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)