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