Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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"
35
36#include "Teuchos_Time.hpp"
37#include "Teuchos_CommandLineProcessor.hpp"
38
39// ADOL-C includes
40#ifdef HAVE_ADOLC
41#ifdef PACKAGE
42#undef PACKAGE
43#endif
44#ifdef PACKAGE_NAME
45#undef PACKAGE_NAME
46#endif
47#ifdef PACKAGE_BUGREPORT
48#undef PACKAGE_BUGREPORT
49#endif
50#ifdef PACKAGE_STRING
51#undef PACKAGE_STRING
52#endif
53#ifdef PACKAGE_TARNAME
54#undef PACKAGE_TARNAME
55#endif
56#ifdef PACKAGE_VERSION
57#undef PACKAGE_VERSION
58#endif
59#ifdef VERSION
60#undef VERSION
61#endif
62//#define ADOLC_TAPELESS
63#define NUMBER_DIRECTIONS 100
64#include "adolc/adouble.h"
65#include "adolc/drivers/drivers.h"
66#include "adolc/interfaces.h"
67#include "adolc/taping.h"
68#endif
69
70#ifdef HAVE_ADIC
71// We have included an ADIC differentiated version of the element fill
72// routine to compare the speed of operator overloading to source
73// transformation. To run this code, all that is necessary is to turn
74// on the ADIC TPL. However to modify the code, it is necessary to
75// re-run the ADIC source transformation tool. To do so, first update
76// the changes to adic_element_fill.c. Then set the following environment
77// variables:
78// export ADIC_ARCH=linux
79// export ADIC=/home/etphipp/AD_libs/adic
80// Next run ADIC via in the tests/performance source directory:
81// ${ADIC}/bin/linux/adiC -vd gradient -i adic_element_fill.init
82// Finally, copy the resulting differentiated function in adic_element_fill.ad.c
83// back into this file. The function will need to be edited by changing
84// the allocation of s to a std::vector<DERIV_TYPE> (the compiler
85// doesn't seem to like malloc), and commenting out the g_filenum lines.
86#define ad_GRAD_MAX 130
87#include "ad_deriv.h"
88#endif
89
90
91// A performance test that computes a finite-element-like Jacobian using
92// several Fad variants
93
94struct ElemData {
95 static const unsigned int nqp = 2;
96 static const unsigned int nnode = 2;
97 double w[nqp], jac[nqp], phi[nqp][nnode], dphi[nqp][nnode];
98 unsigned int gid[nnode];
99
100 ElemData(double mesh_spacing) {
101 // Quadrature points
102 double xi[nqp];
103 xi[0] = -1.0 / std::sqrt(3.0);
104 xi[1] = 1.0 / std::sqrt(3.0);
105
106 for (unsigned int i=0; i<nqp; i++) {
107 // Weights
108 w[i] = 1.0;
109
110 // Element transformation Jacobian
111 jac[i] = 0.5*mesh_spacing;
112
113 // Shape functions & derivatives
114 phi[i][0] = 0.5*(1.0 - xi[i]);
115 phi[i][1] = 0.5*(1.0 + xi[i]);
116 dphi[i][0] = -0.5;
117 dphi[i][1] = 0.5;
118 }
119 }
120};
121
122template <class FadType>
124 unsigned int neqn,
125 const std::vector<double>& x,
126 std::vector<FadType>& x_fad) {
127 for (unsigned int node=0; node<e.nnode; node++)
128 for (unsigned int eqn=0; eqn<neqn; eqn++)
129 x_fad[node*neqn+eqn] = FadType(e.nnode*neqn, node*neqn+eqn,
130 x[e.gid[node]*neqn+eqn]);
131}
132
133template <class T>
135 unsigned int neqn,
136 const std::vector<T>& x,
137 std::vector<T>& u,
138 std::vector<T>& du,
139 std::vector<T>& f) {
140 // Construct element solution, derivative
141 for (unsigned int qp=0; qp<e.nqp; qp++) {
142 for (unsigned int eqn=0; eqn<neqn; eqn++) {
143 u[qp*neqn+eqn] = 0.0;
144 du[qp*neqn+eqn] = 0.0;
145 for (unsigned int node=0; node<e.nnode; node++) {
146 u[qp*neqn+eqn] += x[node*neqn+eqn] * e.phi[qp][node];
147 du[qp*neqn+eqn] += x[node*neqn+eqn] * e.dphi[qp][node];
148 }
149 }
150 }
151
152 // Compute sum of equations for coupling
153 std::vector<T> s(e.nqp);
154 for (unsigned int qp=0; qp<e.nqp; qp++) {
155 s[qp] = 0.0;
156 for (unsigned int eqn=0; eqn<neqn; eqn++)
157 s[qp] += u[qp*neqn+eqn]*u[qp*neqn+eqn];
158 }
159
160 // Evaluate element residual
161 for (unsigned int node=0; node<e.nnode; node++) {
162 for (unsigned int eqn=0; eqn<neqn; eqn++) {
163 unsigned int row = node*neqn+eqn;
164 f[row] = 0.0;
165 for (unsigned int qp=0; qp<e.nqp; qp++) {
166 double c1 = e.w[qp]*e.jac[qp];
167 double c2 = -e.dphi[qp][node]/(e.jac[qp]*e.jac[qp]);
168 f[row] +=
169 c1*(c2*du[qp*neqn+eqn] + e.phi[qp][node]*s[qp]*exp(u[qp*neqn+eqn]));
170 }
171 }
172 }
173}
174
175template <class FadType>
177 unsigned int neqn,
178 const std::vector<FadType>& f_fad,
179 std::vector<double>& f,
180 std::vector< std::vector<double> >& jac) {
181 for (unsigned int eqn_row=0; eqn_row<neqn; eqn_row++) {
182 f[e.gid[0]*neqn+eqn_row] += f_fad[0*neqn+eqn_row].val();
183 f[e.gid[1]*neqn+eqn_row] += f_fad[1*neqn+eqn_row].val();
184 for (unsigned int node_col=0; node_col<e.nnode; node_col++) {
185 for (unsigned int eqn_col=0; eqn_col<neqn; eqn_col++) {
186 unsigned int col = node_col*neqn+eqn_col;
187 unsigned int next_col = (node_col+1)*neqn+eqn_col;
188 jac[e.gid[0]*neqn+eqn_row][next_col] +=
189 f_fad[0*neqn+eqn_row].fastAccessDx(col);
190 jac[e.gid[1]*neqn+eqn_row][col] +=
191 f_fad[1*neqn+eqn_row].fastAccessDx(col);
192 }
193 }
194 }
195}
196
197template <class FadType>
198double fad_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
199 double mesh_spacing) {
200 ElemData e(mesh_spacing);
201
202 // Solution vector, residual, jacobian
203 std::vector<double> x(num_nodes*num_eqns), f(num_nodes*num_eqns);
204 std::vector< std::vector<double> > jac(num_nodes*num_eqns);
205 for (unsigned int node_row=0; node_row<num_nodes; node_row++) {
206 for (unsigned int eqn_row=0; eqn_row<num_eqns; eqn_row++) {
207 unsigned int row = node_row*num_eqns + eqn_row;
208 x[row] = (mesh_spacing*node_row - 0.5)*(mesh_spacing*node_row - 0.5);
209 f[row] = 0.0;
210 jac[row] = std::vector<double>((e.nnode+1)*num_eqns);
211 for (unsigned int node_col=0; node_col<e.nnode+1; node_col++) {
212 for (unsigned int eqn_col=0; eqn_col<num_eqns; eqn_col++) {
213 unsigned int col = node_col*num_eqns + eqn_col;
214 jac[row][col] = 0.0;
215 }
216 }
217 }
218 }
219
220 Teuchos::Time timer("FE Fad Jacobian Fill", false);
221 timer.start(true);
222 std::vector<FadType> x_fad(e.nnode*num_eqns), f_fad(e.nnode*num_eqns);
223 std::vector<FadType> u(e.nqp*num_eqns), du(e.nqp*num_eqns);
224 for (unsigned int i=0; i<num_nodes-1; i++) {
225 e.gid[0] = i;
226 e.gid[1] = i+1;
227
228 fad_init_fill(e, num_eqns, x, x_fad);
229 template_element_fill(e, num_eqns, x_fad, u, du, f_fad);
230 fad_process_fill(e, num_eqns, f_fad, f, jac);
231 }
232 timer.stop();
233
234 // std::cout << "Fad Residual = " << std::endl;
235 // for (unsigned int i=0; i<num_nodes*num_eqns; i++)
236 // std::cout << "\t" << f[i] << std::endl;
237
238 // std::cout.precision(8);
239 // std::cout << "Fad Jacobian = " << std::endl;
240 // for (unsigned int i=0; i<num_nodes*num_eqns; i++) {
241 // std::cout << "\t";
242 // for (unsigned int j=0; j<(e.nnode+1)*num_eqns; j++)
243 // std::cout << jac[i][j] << "\t";
244 // std::cout << std::endl;
245 // }
246
247 return timer.totalElapsedTime();
248}
249
250double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
251 double mesh_spacing);
252double residual_fill(unsigned int num_nodes, unsigned int num_eqns,
253 double mesh_spacing);
254
255#ifdef HAVE_ADOLC
256#ifndef ADOLC_TAPELESS
257void adolc_init_fill(bool retape,
258 int tag,
259 const ElemData& e,
260 unsigned int neqn,
261 const std::vector<double>& x,
262 std::vector<double>& x_local,
263 std::vector<adouble>& x_ad);
264
265void adolc_process_fill(bool retape,
266 int tag,
267 const ElemData& e,
268 unsigned int neqn,
269 std::vector<double>& x_local,
270 std::vector<adouble>& f_ad,
271 std::vector<double>& f,
272 std::vector<double>& f_local,
273 std::vector< std::vector<double> >& jac,
274 double **seed,
275 double **jac_local);
276
277double adolc_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
278 double mesh_spacing);
279
280double adolc_retape_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
281 double mesh_spacing);
282
283#else
284
285void adolc_tapeless_init_fill(const ElemData& e,
286 unsigned int neqn,
287 const std::vector<double>& x,
288 std::vector<adtl::adouble>& x_ad);
289
290void adolc_tapeless_process_fill(const ElemData& e,
291 unsigned int neqn,
292 const std::vector<adtl::adouble>& f_ad,
293 std::vector<double>& f,
294 std::vector< std::vector<double> >& jac);
295
296double adolc_tapeless_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
297 double mesh_spacing);
298
299#endif
300#endif
301
302#ifdef HAVE_ADIC
303void adic_init_fill(const ElemData& e,
304 unsigned int neqn,
305 const std::vector<double>& x,
306 std::vector<DERIV_TYPE>& x_fad);
307void adic_element_fill(ElemData *e,unsigned int neqn,DERIV_TYPE *x,DERIV_TYPE *u,DERIV_TYPE *du,DERIV_TYPE *f);
308void adic_process_fill(const ElemData& e,
309 unsigned int neqn,
310 const std::vector<DERIV_TYPE>& f_fad,
311 std::vector<double>& f,
312 std::vector< std::vector<double> >& jac);
313double adic_jac_fill(unsigned int num_nodes, unsigned int num_eqns,
314 double mesh_spacing);
315inline void AD_Init(int arg0) {
316 ad_AD_GradInit(arg0);
317}
318inline void AD_Final() {
319 ad_AD_GradFinal();
320}
321#endif
322
323#endif
exp(expr.val())
void AD_Final()
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)
Sacado::Fad::DFad< double > FadType
void fad_init_fill(const ElemData &e, unsigned int neqn, const std::vector< double > &x, std::vector< FadType > &x_fad)
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 residual_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
double fad_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
double analytic_jac_fill(unsigned int num_nodes, unsigned int num_eqns, double mesh_spacing)
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
ElemData(double mesh_spacing)
InactiveDouble ** phi
InactiveDouble ** dphi