Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
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
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
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
116double 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
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
183double 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
221void 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
238void 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
275double 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
354double 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
426void 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
440void 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
461double 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
516void 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/**********************************************************************/
549void adic_element_fill(ElemData *e,unsigned int neqn, DERIV_TYPE *x,DERIV_TYPE *u,DERIV_TYPE *du,DERIV_TYPE *f) {
550unsigned int var_0, var_1, var_2, var_3, var_4, var_5, var_6, var_7;
551DERIV_TYPE var_8;
552double 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; ) {
621unsigned 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
659void 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
680double 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
exp(expr.val())
void AD_Final()
#define DERIV_val(a)
Definition ad_deriv.h:40
void AD_Init(int)
void adic_element_fill(ElemData *e, unsigned int neqn, const DERIV_TYPE *x, DERIV_TYPE *u, DERIV_TYPE *du, DERIV_TYPE *f)
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)
double residual_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
void analytic_init_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &x, std::vector< double > &x_local)
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
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)
void residual_process_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &f_local, std::vector< double > &f)
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
InactiveDouble ** phi
InactiveDouble ** dphi