MAST
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
check_sensitivity.h
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 #ifndef __mast_check_sensitivity_h__
21 #define __mast_check_sensitivity_h__
22 
23 // MAST includes
24 #include "tests/base/test_comparisons.h"
26 #include "elasticity/structural_discipline.h"
30 #include "base/parameter.h"
33 
34 
35 // libMesh includes
36 #include "libmesh/numeric_vector.h"
37 
38 
39 namespace MAST {
40 
41  template <typename ValType>
42  void check_sensitivity(ValType& v) {
43 
44 
45  const Real
46  delta = 1.e-5,
47  tol = 1.e-3;
48 
49  // verify the sensitivity solution of this system
51  sol,
52  dsol,
53  dsol_fd;
54 
55  const libMesh::NumericVector<Real>& sol_vec = v.solve();
56 
57  // make sure that each stress object has a single stored value
58  for (unsigned int i=0; i<v._outputs.size(); i++) {
59  BOOST_CHECK((v._outputs[i]->n_elem_in_storage() == 1));
60  }
61 
62  const unsigned int
63  n_dofs = sol_vec.size(),
64  n_elems = v._mesh->n_elem();
65 
66  const Real
67  p_val = 2.;
68 
69  // store the stress values for sensitivity calculations
70  // make sure that the current setup has specified one stress output
71  // per element
72  libmesh_assert(v._outputs.size() == n_elems);
74  stress0 = RealVectorX::Zero(n_elems),
75  dstressdp = RealVectorX::Zero(n_elems),
76  dstressdp_fd = RealVectorX::Zero(n_elems);
77 
78  for (unsigned int i=0; i<n_elems; i++) {
79  // the call to all elements should actually include a single element only
80  // the p-norm used is for p=2.
81  stress0(i) = v._outputs[i]->von_Mises_p_norm_functional_for_all_elems(p_val);
82  }
83 
84 
85  sol = RealVectorX::Zero(n_dofs);
86  dsol = RealVectorX::Zero(n_dofs);
87 
88  // copy the solution for later use
89  for (unsigned int i=0; i<n_dofs; i++) {
90  sol(i) = sol_vec(i);
91  }
92 
93  // now clear the stress data structures
94  v.clear_stresss();
95 
96  // now iterate over all the parameters and calculate the analytical
97  // sensitivity and compare with the numerical sensitivity
98 
99  Real
100  p0 = 0.,
101  dp = 0.;
102 
104  // now evaluate the direct sensitivity
106 
107  for (unsigned int i=0; i<v._params_for_sensitivity.size(); i++ ) {
108 
109  MAST::Parameter& f = *v._params_for_sensitivity[i];
110 
111  dsol = RealVectorX::Zero(n_dofs);
112  dsol_fd = RealVectorX::Zero(n_dofs);
113  dstressdp = RealVectorX::Zero(n_elems);
114  dstressdp_fd = RealVectorX::Zero(n_elems);
115 
116  // calculate the analytical sensitivity
117  // analysis is required at the baseline before sensitivity solution
118  // and the solution has changed after the previous perturbed solution
119  v.solve();
120  const libMesh::NumericVector<Real>& dsol_vec = v.sensitivity_solve(f);
121 
122  // make sure that each stress object has a single stored value
123  for (unsigned int i=0; i<v._outputs.size(); i++) {
124  BOOST_CHECK((v._outputs[i]->n_elem_in_storage() == 1));
125  }
126 
127  // copy the sensitivity solution
128  for (unsigned int j=0; j<n_dofs; j++) {
129  dsol(j) = dsol_vec(j);
130  }
131 
132  // copy the analytical sensitivity of stress values
133  for (unsigned int j=0; j<n_elems; j++) {
134  dstressdp(j) =
135  v._outputs[j]->von_Mises_p_norm_functional_sensitivity_for_all_elems
136  (p_val, &f);
137  }
138 
139  // now clear the stress data structures
140  v.clear_stresss();
141 
142  // now calculate the finite difference sensitivity
143 
144  // identify the perturbation in the parameter
145  p0 = f();
146  (fabs(p0) > 0)? dp=delta*p0 : dp=delta;
147  f() += dp;
148 
149  // solve at the perturbed parameter value
150  const libMesh::NumericVector<Real>& sol_vec1 = v.solve();
151 
152  // make sure that each stress object has a single stored value
153  for (unsigned int i=0; i<v._outputs.size(); i++) {
154  BOOST_CHECK((v._outputs[i]->n_elem_in_storage() == 1));
155  }
156 
157  // copy the perturbed solution
158  for (unsigned int i=0; i<n_dofs; i++) {
159  dsol_fd(i) = sol_vec1(i);
160  }
161 
162  // calculate the finite difference sensitivity for solution
163  dsol_fd -= sol;
164  dsol_fd /= dp;
165 
166  // copy the perturbed stress values
167  for (unsigned int j=0; j<n_elems; j++) {
168  dstressdp_fd(j) =
169  v._outputs[j]->von_Mises_p_norm_functional_for_all_elems(p_val);
170  }
171 
172  // calculate the finite difference sensitivity for stress
173  dstressdp_fd -= stress0;
174  dstressdp_fd /= dp;
175 
176  // now clear the stress data structures
177  v.clear_stresss();
178 
179  // reset the parameter value
180  f() = p0;
181 
182  // now compare the solution sensitivity
183  BOOST_TEST_MESSAGE(" ** dX/dp (total) wrt : " << f.name() << " **");
184  BOOST_CHECK(MAST::compare_vector( dsol_fd, dsol, tol));
185 
186  // now compare the stress sensitivity
187  BOOST_TEST_MESSAGE(" ** dvm-stress/dp (total) wrt : " << f.name() << " **");
188  BOOST_CHECK(MAST::compare_vector( dstressdp_fd, dstressdp, tol));
189  }
190 
191  v.clear_stresss();
192 
194  // now evaluate the adjoint sensitivity
196 
197 
198  }
199 
200 
201 }
202 
203 
204 #endif // __mast_check_sensitivity_h__
const std::string & name() const
returns the name of this function
Definition: function_base.h:60
bool compare_vector(const RealVectorX &v0, const RealVectorX &v, const Real tol)
This is a scalar function whose value can be changed and one that can be used as a design variable in...
Definition: parameter.h:35
void check_sensitivity(ValType &v)
libMesh::Real Real
Matrix< Real, Dynamic, 1 > RealVectorX