Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fe_jac_fill_funcs.cpp
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 #include "fe_jac_fill_funcs.hpp"
30 
31 void analytic_init_fill(const ElemData& e,
32  unsigned int neqn,
33  const std::vector<double>& x,
34  std::vector<double>& x_local) {
35  for (unsigned int node=0; node<e.nnode; node++)
36  for (unsigned int eqn=0; eqn<neqn; eqn++)
37  x_local[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
38 }
39 
40 void analytic_element_fill(const ElemData& e,
41  unsigned int neqn,
42  const std::vector<double>& x,
43  std::vector<double>& u,
44  std::vector<double>& du,
45  std::vector<double>& f,
46  std::vector< std::vector<double> >& jac) {
47  // Construct element solution, derivative
48  for (unsigned int qp=0; qp<e.nqp; qp++) {
49  for (unsigned int eqn=0; eqn<neqn; eqn++) {
50  u[qp*neqn+eqn] = 0.0;
51  du[qp*neqn+eqn] = 0.0;
52  for (unsigned int node=0; node<e.nnode; node++) {
53  u[qp*neqn+eqn] += x[node*neqn+eqn] * e.phi[qp][node];
54  du[qp*neqn+eqn] += x[node*neqn+eqn] * e.dphi[qp][node];
55  }
56  }
57  }
58 
59  // Compute sum of equations for coupling
60  std::vector<double> s(e.nqp);
61  for (unsigned int qp=0; qp<e.nqp; qp++) {
62  s[qp] = 0.0;
63  for (unsigned int eqn=0; eqn<neqn; eqn++)
64  s[qp] += u[qp*neqn+eqn]*u[qp*neqn+eqn];
65  }
66 
67  // Evaluate element residual
68  for (unsigned int node_row=0; node_row<e.nnode; node_row++) {
69  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
70  unsigned int row = node_row*neqn+eqn_row;
71  f[row] = 0.0;
72  for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
73  for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
74  unsigned int col = node_col*neqn+eqn_col;
75  jac[row][col] = 0.0;
76  }
77  }
78  for (unsigned int qp=0; qp<e.nqp; qp++) {
79  double c1 = e.w[qp]*e.jac[qp];
80  double c2 = -e.dphi[qp][node_row]/(e.jac[qp]*e.jac[qp]);
81  double c3 = exp(u[qp*neqn+eqn_row]);
82  double c4 = e.phi[qp][node_row]*s[qp]*c3;
83  f[row] += c1*(c2*du[qp*neqn+eqn_row] + c4);
84  for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
85  jac[row][node_col*neqn+eqn_row] +=
86  c1*(c2*e.dphi[qp][node_col] + c4*e.phi[qp][node_col]);
87  for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++)
88  jac[row][node_col*neqn+eqn_col] +=
89  2.0*c1*e.phi[qp][node_row]*u[qp*neqn+eqn_col]*e.phi[qp][node_col]*c3;
90  }
91  }
92  }
93  }
94 }
95 
96 void analytic_process_fill(const ElemData& e,
97  unsigned int neqn,
98  const std::vector<double>& f_local,
99  const std::vector< std::vector<double> >& jac_local,
100  std::vector<double>& f,
101  std::vector< std::vector<double> >& jac) {
102  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
103  f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
104  f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
105  for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
106  for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
107  unsigned int col = node_col*neqn+eqn_col;
108  unsigned int next_col = (node_col+1)*neqn+eqn_col;
109  jac[e.gid[0]*neqn+eqn_row][next_col] += jac_local[0*neqn+eqn_row][col];
110  jac[e.gid[1]*neqn+eqn_row][col] += jac_local[1*neqn+eqn_row][col];
111  }
112  }
113  }
114 }
115 
116 double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
117  double mesh_spacing) {
118  ElemData e(mesh_spacing);
119 
120  // Solution vector, residual, jacobian
121  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
122  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
123  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
124  for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
125  unsigned int row = node_row*num_eqns + eqn_row;
126  x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
127  f[row] = 0.0;
128  jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
129  for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
130  for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
131  unsigned int col = node_col*num_eqns + eqn_col;
132  jac[row][col] = 0.0;
133  }
134  }
135  }
136  }
137 
138  Teuchos::Time timer("FE Analytic Jacobian Fill", false);
139  timer.start(true);
140  std::vector<double> x_local(e.nnode*num_eqns), f_local(e.nnode*num_eqns);
141  std::vector< std::vector<double> > jac_local(e.nnode*num_eqns);
142  std::vector<double> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
143  for (unsigned int i=0; i<e.nnode*num_eqns; i++)
144  jac_local[i] = std::vector<double>(e.nnode*num_eqns);
145  for (unsigned int i=0; i<num_nodes-1; i++) {
146  e.gid[0] = i;
147  e.gid[1] = i+1;
148 
149  analytic_init_fill(e, num_eqns, x, x_local);
150  analytic_element_fill(e, num_eqns, x_local, u, du, f_local, jac_local);
151  analytic_process_fill(e, num_eqns, f_local, jac_local, f, jac);
152  }
153  timer.stop();
154 
155  // std::cout.precision(8);
156  // std::cout << "Analytic Residual = " << std::endl;
157  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
158  // std::cout << "\t" << f[i] << std::endl;
159 
160  // std::cout.precision(8);
161  // std::cout.setf(std::ios::scientific);
162  // std::cout << "Analytic Jacobian = " << std::endl;
163  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
164  // std::cout << "\t";
165  // for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
166  // std::cout << jac[i][j] << "\t";
167  // std::cout << std::endl;
168  // }
169 
170  return timer.totalElapsedTime();
171 }
172 
173 void residual_process_fill(const ElemData& e,
174  unsigned int neqn,
175  const std::vector<double>& f_local,
176  std::vector<double>& f) {
177  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
178  f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
179  f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
180  }
181 }
182 
183 double residual_fill(unsigned int num_nodes, unsigned int num_eqns,
184  double mesh_spacing) {
185  ElemData e(mesh_spacing);
186 
187  // Solution vector, residual, jacobian
188  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
189  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
190  for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
191  unsigned int row = node_row*num_eqns + eqn_row;
192  x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
193  f[row] = 0.0;
194  }
195  }
196 
197  Teuchos::Time timer("FE Residual Fill", false);
198  timer.start(true);
199  std::vector<double> x_local(e.nnode*num_eqns), f_local(e.nnode*num_eqns);
200  std::vector<double> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
201  for (unsigned int i=0; i<num_nodes-1; i++) {
202  e.gid[0] = i;
203  e.gid[1] = i+1;
204 
205  analytic_init_fill(e, num_eqns, x, x_local);
206  template_element_fill(e, num_eqns, x_local, u, du, f_local);
207  residual_process_fill(e, num_eqns, f_local, f);
208  }
209  timer.stop();
210 
211  // std::cout.precision(8);
212  // std::cout << "Analytic Residual = " << std::endl;
213  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
214  // std::cout << "\t" << f[i] << std::endl;
215 
216  return timer.totalElapsedTime();
217 }
218 
219 #ifdef HAVE_ADOLC
220 #ifndef ADOLC_TAPELESS
221 void adolc_init_fill(bool retape,
222  int tag,
223  const ElemData& e,
224  unsigned int neqn,
225  const std::vector<double>& x,
226  std::vector<double>& x_local,
227  std::vector<adouble>& x_ad) {
228  if (retape)
229  trace_on(tag);
230  for (unsigned int node=0; node<e.nnode; node++)
231  for (unsigned int eqn=0; eqn<neqn; eqn++) {
232  x_local[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
233  if (retape)
234  x_ad[node*neqn+eqn] <<= x_local[node*neqn+eqn];
235  }
236 }
237 
238 void adolc_process_fill(bool retape,
239  int tag,
240  const ElemData& e,
241  unsigned int neqn,
242  std::vector<double>& x_local,
243  std::vector<adouble>& f_ad,
244  std::vector<double>& f,
245  std::vector<double>& f_local,
246  std::vector< std::vector<double> >& jac,
247  double **seed,
248  double **jac_local) {
249  if (retape) {
250  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++)
251  f_ad[0*neqn+eqn_row] >>= f_local[0*neqn+eqn_row];
252  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++)
253  f_ad[1*neqn+eqn_row] >>= f_local[1*neqn+eqn_row];
254  trace_off();
255  }
256 
257  //jacobian(tag, neqn*e.nnode, neqn*e.nnode, &x_local[0], jac_local);
258  forward(tag, neqn*e.nnode, neqn*e.nnode, neqn*e.nnode, &x_local[0],
259  seed, &f_local[0], jac_local);
260 
261  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
262  f[e.gid[0]*neqn+eqn_row] += f_local[0*neqn+eqn_row];
263  f[e.gid[1]*neqn+eqn_row] += f_local[1*neqn+eqn_row];
264  for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
265  for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
266  unsigned int col = node_col*neqn+eqn_col;
267  unsigned int next_col = (node_col+1)*neqn+eqn_col;
268  jac[e.gid[0]*neqn+eqn_row][next_col] += jac_local[0*neqn+eqn_row][col];
269  jac[e.gid[1]*neqn+eqn_row][col] += jac_local[1*neqn+eqn_row][col];
270  }
271  }
272  }
273 }
274 
275 double adolc_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
276  double mesh_spacing) {
277  ElemData e(mesh_spacing);
278 
279  // Solution vector, residual, jacobian
280  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
281  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
282  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
283  for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
284  unsigned int row = node_row*num_eqns + eqn_row;
285  x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
286  f[row] = 0.0;
287  jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
288  for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
289  for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
290  unsigned int col = node_col*num_eqns + eqn_col;
291  jac[row][col] = 0.0;
292  }
293  }
294  }
295  }
296 
297  Teuchos::Time timer("FE ADOL-C Jacobian Fill", false);
298  timer.start(true);
299  std::vector<adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
300  std::vector<adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
301  std::vector<double> x_local(e.nnode*num_eqns);
302  std::vector<double> f_local(e.nnode*num_eqns);
303  double **jac_local = new double*[e.nnode*num_eqns];
304  double **seed = new double*[e.nnode*num_eqns];
305  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
306  jac_local[i] = new double[e.nnode*num_eqns];
307  seed[i] = new double[e.nnode*num_eqns];
308  for (unsigned int j=0; j<e.nnode*num_eqns; j++)
309  seed[i][j] = 0.0;
310  seed[i][i] = 1.0;
311  }
312 
313  // Tape first element
314  e.gid[0] = 0;
315  e.gid[1] = 1;
316  adolc_init_fill(true, 0, e, num_eqns, x, x_local, x_ad);
317  template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
318  adolc_process_fill(true, 0, e, num_eqns, x_local, f_ad, f, f_local,
319  jac, seed, jac_local);
320 
321  // Now do remaining fills reusing tape
322  for (unsigned int i=1; i<num_nodes-1; i++) {
323  e.gid[0] = i;
324  e.gid[1] = i+1;
325 
326  adolc_init_fill(false, 0, e, num_eqns, x, x_local, x_ad);
327  adolc_process_fill(false, 0, e, num_eqns, x_local, f_ad, f, f_local,
328  jac, seed, jac_local);
329  }
330  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
331  delete [] jac_local[i];
332  delete [] seed[i];
333  }
334  delete [] jac_local;
335  delete [] seed;
336  timer.stop();
337 
338  // std::cout << "ADOL-C Residual = " << std::endl;
339  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
340  // std::cout << "\t" << f[i] << std::endl;
341 
342  // std::cout.precision(8);
343  // std::cout << "ADOL-C Jacobian = " << std::endl;
344  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
345  // std::cout << "\t";
346  // for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
347  // std::cout << jac[i][j] << "\t";
348  // std::cout << std::endl;
349  // }
350 
351  return timer.totalElapsedTime();
352 }
353 
354 double adolc_retape_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
355  double mesh_spacing) {
356  ElemData e(mesh_spacing);
357 
358  // Solution vector, residual, jacobian
359  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
360  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
361  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
362  for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
363  unsigned int row = node_row*num_eqns + eqn_row;
364  x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
365  f[row] = 0.0;
366  jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
367  for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
368  for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
369  unsigned int col = node_col*num_eqns + eqn_col;
370  jac[row][col] = 0.0;
371  }
372  }
373  }
374  }
375 
376  Teuchos::Time timer("FE ADOL-C Retape Jacobian Fill", false);
377  timer.start(true);
378  std::vector<adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
379  std::vector<adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
380  std::vector<double> x_local(e.nnode*num_eqns);
381  std::vector<double> f_local(e.nnode*num_eqns);
382  double **jac_local = new double*[e.nnode*num_eqns];
383  double **seed = new double*[e.nnode*num_eqns];
384  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
385  jac_local[i] = new double[e.nnode*num_eqns];
386  seed[i] = new double[e.nnode*num_eqns];
387  for (unsigned int j=0; j<e.nnode*num_eqns; j++)
388  seed[i][j] = 0.0;
389  seed[i][i] = 1.0;
390  }
391  for (unsigned int i=0; i<num_nodes-1; i++) {
392  e.gid[0] = i;
393  e.gid[1] = i+1;
394 
395  adolc_init_fill(true, 1, e, num_eqns, x, x_local, x_ad);
396  template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
397  adolc_process_fill(true, 1, e, num_eqns, x_local, f_ad, f, f_local,
398  jac, seed, jac_local);
399  }
400  for (unsigned int i=0; i<e.nnode*num_eqns; i++) {
401  delete [] jac_local[i];
402  delete [] seed[i];
403  }
404  delete [] jac_local;
405  delete [] seed;
406  timer.stop();
407 
408  // std::cout << "ADOL-C Residual (retaped) = " << std::endl;
409  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
410  // std::cout << "\t" << f[i] << std::endl;
411 
412  // std::cout.precision(8);
413  // std::cout << "ADOL-C Jacobian (retaped) = " << std::endl;
414  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
415  // std::cout << "\t";
416  // for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
417  // std::cout << jac[i][j] << "\t";
418  // std::cout << std::endl;
419  // }
420 
421  return timer.totalElapsedTime();
422 }
423 
424 #else
425 
426 void adolc_tapeless_init_fill(const ElemData& e,
427  unsigned int neqn,
428  const std::vector<double>& x,
429  std::vector<adtl::adouble>& x_ad) {
430  for (unsigned int node=0; node<e.nnode; node++)
431  for (unsigned int eqn=0; eqn<neqn; eqn++) {
432  x_ad[node*neqn+eqn] = x[e.gid[node]*neqn+eqn];
433  for (unsigned int i=0; i<neqn*e.nnode; i++)
434  x_ad[node*neqn+eqn].setADValue(i, 0.0);
435  x_ad[node*neqn+eqn].setADValue(node*neqn+eqn, 1.0);
436  }
437 
438 }
439 
440 void adolc_tapeless_process_fill(const ElemData& e,
441  unsigned int neqn,
442  const std::vector<adtl::adouble>& f_ad,
443  std::vector<double>& f,
444  std::vector< std::vector<double> >& jac) {
445  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
446  f[e.gid[0]*neqn+eqn_row] += f_ad[0*neqn+eqn_row].getValue();
447  f[e.gid[1]*neqn+eqn_row] += f_ad[1*neqn+eqn_row].getValue();
448  for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
449  for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
450  unsigned int col = node_col*neqn+eqn_col;
451  unsigned int next_col = (node_col+1)*neqn+eqn_col;
452  jac[e.gid[0]*neqn+eqn_row][next_col] +=
453  f_ad[0*neqn+eqn_row].getADValue(col);
454  jac[e.gid[1]*neqn+eqn_row][col] +=
455  f_ad[1*neqn+eqn_row].getADValue(col);
456  }
457  }
458  }
459 }
460 
461 double adolc_tapeless_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
462  double mesh_spacing) {
463  ElemData e(mesh_spacing);
464 
465  // Solution vector, residual, jacobian
466  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
467  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
468  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
469  for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
470  unsigned int row = node_row*num_eqns + eqn_row;
471  x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
472  f[row] = 0.0;
473  jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
474  for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
475  for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
476  unsigned int col = node_col*num_eqns + eqn_col;
477  jac[row][col] = 0.0;
478  }
479  }
480  }
481  }
482 
483  Teuchos::Time timer("FE Tapeless ADOL-C Jacobian Fill", false);
484  timer.start(true);
485  std::vector<adtl::adouble> x_ad(e.nnode*num_eqns), f_ad(e.nnode*num_eqns);
486  std::vector<adtl::adouble> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
487  for (unsigned int i=0; i<num_nodes-1; i++) {
488  e.gid[0] = i;
489  e.gid[1] = i+1;
490 
491  adolc_tapeless_init_fill(e, num_eqns, x, x_ad);
492  template_element_fill(e, num_eqns, x_ad, u, du, f_ad);
493  adolc_tapeless_process_fill(e, num_eqns, f_ad, f, jac);
494  }
495  timer.stop();
496 
497  // std::cout << "Tapeless ADOL-C Residual = " << std::endl;
498  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
499  // std::cout << "\t" << f[i] << std::endl;
500 
501  // std::cout.precision(8);
502  // std::cout << "Tapeless ADOL-C Jacobian = " << std::endl;
503  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
504  // std::cout << "\t";
505  // for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
506  // std::cout << jac[i][j] << "\t";
507  // std::cout << std::endl;
508  // }
509 
510  return timer.totalElapsedTime();
511 }
512 #endif
513 #endif
514 
515 #ifdef HAVE_ADIC
516 void adic_init_fill(const ElemData& e,
517  unsigned int neqn,
518  const std::vector<double>& x,
519  std::vector<DERIV_TYPE>& x_fad) {
520  static bool first = true;
521  for (unsigned int node=0; node<e.nnode; node++)
522  for (unsigned int eqn=0; eqn<neqn; eqn++) {
523  x_fad[node*neqn+eqn].value = x[e.gid[node]*neqn+eqn];
524  if (first)
525  ad_AD_SetIndep(x_fad[node*neqn+eqn]);
526  }
527  if (first) {
528  ad_AD_SetIndepDone();
529  first = false;
530  }
531 }
532 
533 /************************** DISCLAIMER ********************************/
534 /* */
535 /* This file was generated on 04/13/12 11:10:49 by the version of */
536 /* ADIC 1.2.3 compiled on 04/14/09 12:39:01 */
537 /* */
538 /* ADIC was prepared as an account of work sponsored by an */
539 /* agency of the United States Government and the University of */
540 /* Chicago. NEITHER THE AUTHOR(S), THE UNITED STATES GOVERNMENT */
541 /* NOR ANY AGENCY THEREOF, NOR THE UNIVERSITY OF CHICAGO, INCLUDING */
542 /* ANY OF THEIR EMPLOYEES OR OFFICERS, MAKES ANY WARRANTY, EXPRESS */
543 /* OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR */
544 /* THE ACCURACY, COMPLETENESS, OR USEFULNESS OF ANY INFORMATION OR */
545 /* PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT INFRINGE */
546 /* PRIVATELY OWNED RIGHTS. */
547 /* */
548 /**********************************************************************/
549 void adic_element_fill(ElemData *e,unsigned int neqn, DERIV_TYPE *x,DERIV_TYPE *u,DERIV_TYPE *du,DERIV_TYPE *f) {
550 unsigned int var_0, var_1, var_2, var_3, var_4, var_5, var_6, var_7;
551 DERIV_TYPE var_8;
552 double adji_0;
553  double loc_0;
554  double loc_1;
555  double loc_2;
556  double loc_3;
557  double loc_4;
558  double loc_5;
559  double loc_6;
560  double loc_7;
561  double loc_8;
562  double loc_9;
563  double adj_0;
564  double adj_1;
565  double adj_2;
566  double adj_3;
567 
568  // static int g_filenum = 0;
569  // if (g_filenum == 0) {
570  // adintr_ehsfid(&g_filenum, __FILE__, "adic_element_fill");
571  // }
572  for (unsigned int qp = 0; qp < e->nqp; ) {
573  for (unsigned int eqn = 0; eqn < neqn; ) {
574  {
575  ad_grad_axpy_0(&(u[qp * neqn + eqn]));
576  DERIV_val(u[qp * neqn + eqn]) = 0.0;
577  }
578  {
579  ad_grad_axpy_0(&(du[qp * neqn + eqn]));
580  DERIV_val(du[qp * neqn + eqn]) = 0.0;
581  }
582  for (unsigned int node = 0; node < e->nnode; ) {
583  {
584  loc_0 = DERIV_val(x[node * neqn + eqn]) * e->phi[qp][node];
585  loc_1 = DERIV_val(u[qp * neqn + eqn]) + loc_0;
586  ad_grad_axpy_2(&(u[qp * neqn + eqn]), 1.000000000000000e+00, &(u[qp * neqn + eqn]), e->phi[qp][node], &(x[node * neqn + eqn]));
587  DERIV_val(u[qp * neqn + eqn]) = loc_1;
588  }
589  {
590  loc_0 = DERIV_val(x[node * neqn + eqn]) * e->dphi[qp][node];
591  loc_1 = DERIV_val(du[qp * neqn + eqn]) + loc_0;
592  ad_grad_axpy_2(&(du[qp * neqn + eqn]), 1.000000000000000e+00, &(du[qp * neqn + eqn]), e->dphi[qp][node], &(x[node * neqn + eqn]));
593  DERIV_val(du[qp * neqn + eqn]) = loc_1;
594  }
595  var_2 = node++;
596  }
597  var_1 = eqn++;
598  }
599  var_0 = qp++;
600  }
601 //DERIV_TYPE *s = malloc(e->nqp * sizeof (DERIV_TYPE ));
602  std::vector<DERIV_TYPE> s(e->nqp);
603  for (unsigned int qp = 0; qp < e->nqp; ) {
604  {
605  ad_grad_axpy_0(&(s[qp]));
606  DERIV_val(s[qp]) = 0.0;
607  }
608  for (unsigned int eqn = 0; eqn < neqn; ) {
609  {
610  loc_0 = DERIV_val(u[qp * neqn + eqn]) * DERIV_val(u[qp * neqn + eqn]);
611  loc_1 = DERIV_val(s[qp]) + loc_0;
612  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]));
613  DERIV_val(s[qp]) = loc_1;
614  }
615  var_4 = eqn++;
616  }
617  var_3 = qp++;
618  }
619  for (unsigned int node = 0; node < e->nnode; ) {
620  for (unsigned int eqn = 0; eqn < neqn; ) {
621 unsigned int row = node * neqn + eqn;
622  {
623  ad_grad_axpy_0(&(f[row]));
624  DERIV_val(f[row]) = 0.0;
625  }
626  for (unsigned int qp = 0; qp < e->nqp; ) {
627  DERIV_val(var_8) = exp(( DERIV_val(u[qp * neqn + eqn])));
628  adji_0 = DERIV_val(var_8);
629  {
630  ad_grad_axpy_1(&(var_8), adji_0, &(u[qp * neqn + eqn]));
631  }
632  {
633  loc_0 = e->w[qp] * e->jac[qp];
634  loc_1 = -e->dphi[qp][node];
635  loc_2 = e->jac[qp] * e->jac[qp];
636  loc_3 = loc_1 / loc_2;
637  loc_4 = loc_3 * DERIV_val(du[qp * neqn + eqn]);
638  loc_5 = e->phi[qp][node] * DERIV_val(s[qp]);
639  loc_6 = loc_5 * DERIV_val(var_8);
640  loc_7 = loc_4 + loc_6;
641  loc_8 = loc_0 * loc_7;
642  loc_9 = DERIV_val(f[row]) + loc_8;
643  adj_0 = loc_5 * loc_0;
644  adj_1 = DERIV_val(var_8) * loc_0;
645  adj_2 = e->phi[qp][node] * adj_1;
646  adj_3 = loc_3 * loc_0;
647  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));
648  DERIV_val(f[row]) = loc_9;
649  }
650  var_7 = qp++;
651  }
652  var_6 = eqn++;
653  }
654  var_5 = node++;
655  }
656  //free(s);
657 }
658 
659 void adic_process_fill(const ElemData& e,
660  unsigned int neqn,
661  const std::vector<DERIV_TYPE>& f_fad,
662  std::vector<double>& f,
663  std::vector< std::vector<double> >& jac) {
664  for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
665  f[e.gid[0]*neqn+eqn_row] += f_fad[0*neqn+eqn_row].value;
666  f[e.gid[1]*neqn+eqn_row] += f_fad[1*neqn+eqn_row].value;
667  for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
668  for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
669  unsigned int col = node_col*neqn+eqn_col;
670  unsigned int next_col = (node_col+1)*neqn+eqn_col;
671  jac[e.gid[0]*neqn+eqn_row][next_col] +=
672  f_fad[0*neqn+eqn_row].grad[col];
673  jac[e.gid[1]*neqn+eqn_row][col] +=
674  f_fad[1*neqn+eqn_row].grad[col];
675  }
676  }
677  }
678 }
679 
680 double adic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
681  double mesh_spacing) {
682  AD_Init(0);
683  ElemData e(mesh_spacing);
684 
685  // Solution vector, residual, jacobian
686  std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
687  std::vector< std::vector<double> > jac(num_nodes*num_eqns);
688  for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
689  for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
690  unsigned int row = node_row*num_eqns + eqn_row;
691  x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
692  f[row] = 0.0;
693  jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
694  for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
695  for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
696  unsigned int col = node_col*num_eqns + eqn_col;
697  jac[row][col] = 0.0;
698  }
699  }
700  }
701  }
702 
703  Teuchos::Time timer("FE ADIC Jacobian Fill", false);
704  timer.start(true);
705  std::vector<DERIV_TYPE> x_fad(e.nnode*num_eqns), f_fad(e.nnode*num_eqns);
706  std::vector<DERIV_TYPE> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
707  for (unsigned int i=0; i<num_nodes-1; i++) {
708  e.gid[0] = i;
709  e.gid[1] = i+1;
710 
711  adic_init_fill(e, num_eqns, x, x_fad);
712  adic_element_fill(&e, num_eqns, &x_fad[0], &u[0], &du[0], &f_fad[0]);
713  adic_process_fill(e, num_eqns, f_fad, f, jac);
714  }
715  timer.stop();
716  AD_Final();
717 
718  // std::cout.precision(8);
719  // std::cout << "ADIC Residual = " << std::endl;
720  // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
721  // std::cout << "\t" << f[i] << std::endl;
722 
723  // std::cout.precision(8);
724  // std::cout << "ADIC Jacobian = " << std::endl;
725  // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
726  // std::cout << "\t";
727  // for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
728  // std::cout << jac[i][j] << "\t";
729  // std::cout << std::endl;
730  // }
731 
732  return timer.totalElapsedTime();
733 }
734 #endif
735 
f
void f()
AD_Init
void AD_Init(int)
Definition: adic_element_fill.ad.c:136
residual_fill
double residual_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
Definition: fe_jac_fill_funcs.cpp:183
template_element_fill
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)
Definition: fe_jac_fill_funcs.hpp:133
ElemData::jac
InactiveDouble * jac
Definition: adic_element_fill.ad.c:24
ElemData::w
InactiveDouble * w
Definition: adic_element_fill.ad.c:24
DERIV_val
#define DERIV_val(a)
Definition: ad_deriv.h:40
AD_Final
void AD_Final()
Definition: adic_element_fill.ad.c:140
adic_element_fill
void adic_element_fill(ElemData *e, unsigned int neqn, const DERIV_TYPE *x, DERIV_TYPE *u, DERIV_TYPE *du, DERIV_TYPE *f)
Definition: adic_element_fill.ad.c:28
analytic_jac_fill
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
Definition: fe_jac_fill_funcs.cpp:116
analytic_element_fill
void analytic_element_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &x, std::vector< double > &u, std::vector< double > &du, std::vector< double > &f, std::vector< std::vector< double > > &jac)
Definition: fe_jac_fill_funcs.cpp:40
analytic_process_fill
void analytic_process_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &f_local, const std::vector< std::vector< double > > &jac_local, std::vector< double > &f, std::vector< std::vector< double > > &jac)
Definition: fe_jac_fill_funcs.cpp:96
DERIV_TYPE
Definition: ad_deriv.h:35
Teuchos::Time::start
void start(bool reset=false)
analytic_init_fill
void analytic_init_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &x, std::vector< double > &x_local)
Definition: fe_jac_fill_funcs.cpp:31
Teuchos::Time::stop
double stop()
fe_jac_fill_funcs.hpp
residual_process_fill
void residual_process_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &f_local, std::vector< double > &f)
Definition: fe_jac_fill_funcs.cpp:173
ElemData::phi
InactiveDouble ** phi
Definition: adic_element_fill.ad.c:24
ElemData::dphi
InactiveDouble ** dphi
Definition: adic_element_fill.ad.c:24
Teuchos::Time::totalElapsedTime
double totalElapsedTime(bool readCurrentTime=false) const
exp
exp(expr.val())
Generated by   doxygen 1.8.5