MAST
second_order_newmark_transient_solver.cpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2019 Manav Bhatia
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 // MAST includes
23 #include "base/elem_base.h"
25 #include "base/nonlinear_system.h"
26 
27 
28 // libMesh includes
29 #include "libmesh/numeric_vector.h"
30 
31 
34 beta(0.25),
35 gamma(0.5)
36 { }
37 
38 
40 { }
41 
42 
43 
44 
45 void
47 set_element_data(const std::vector<libMesh::dof_id_type>& dof_indices,
48  const std::vector<libMesh::NumericVector<Real>*>& sols){
49 
50  libmesh_assert_equal_to(sols.size(), 3);
51 
52  const unsigned int n_dofs = (unsigned int)dof_indices.size();
53 
54  // get the current state and velocity estimates
55  // also get the current discrete velocity replacement
57  sol = RealVectorX::Zero(n_dofs),
58  vel = RealVectorX::Zero(n_dofs),
59  accel = RealVectorX::Zero(n_dofs);
60 
61 
62  // get the references to current and previous sol and velocity
63  const libMesh::NumericVector<Real>
64  &sol_global = *sols[0],
65  &vel_global = *sols[1],
66  &acc_global = *sols[2];
67 
68  for (unsigned int i=0; i<n_dofs; i++) {
69 
70  sol(i) = sol_global(dof_indices[i]);
71  vel(i) = vel_global(dof_indices[i]);
72  accel(i) = acc_global(dof_indices[i]);
73  }
74 
75  _assembly_ops->set_elem_solution(sol);
76  _assembly_ops->set_elem_velocity(vel);
77  _assembly_ops->set_elem_acceleration(accel);
78 }
79 
80 
81 
82 void
84 extract_element_sensitivity_data(const std::vector<libMesh::dof_id_type>& dof_indices,
85  const std::vector<libMesh::NumericVector<Real>*>& sols,
86  std::vector<RealVectorX>& local_sols) {
87 
88  libmesh_assert_equal_to(sols.size(), 3);
89 
90  const unsigned int n_dofs = (unsigned int)dof_indices.size();
91 
92  local_sols.resize(3);
93 
95  &sol = local_sols[0],
96  &vel = local_sols[1],
97  &accel = local_sols[2];
98 
99  sol = RealVectorX::Zero(n_dofs);
100  vel = RealVectorX::Zero(n_dofs);
101  accel = RealVectorX::Zero(n_dofs);
102 
103 
104  // get the references to current and previous sol and velocity
105  const libMesh::NumericVector<Real>
106  &sol_global = *sols[0],
107  &vel_global = *sols[1],
108  &acc_global = *sols[2];
109 
110  for (unsigned int i=0; i<n_dofs; i++) {
111 
112  sol(i) = sol_global(dof_indices[i]);
113  vel(i) = vel_global(dof_indices[i]);
114  accel(i) = acc_global(dof_indices[i]);
115  }
116 
117 
118 }
119 
120 
121 
122 void
124 set_element_perturbed_data(const std::vector<libMesh::dof_id_type>& dof_indices,
125  const std::vector<libMesh::NumericVector<Real>*>& sols){
126 
127  libmesh_assert_equal_to(sols.size(), 3);
128 
129  const unsigned int n_dofs = (unsigned int)dof_indices.size();
130 
131  // get the current state and velocity estimates
132  // also get the current discrete velocity replacement
134  sol = RealVectorX::Zero(n_dofs),
135  vel = RealVectorX::Zero(n_dofs),
136  accel = RealVectorX::Zero(n_dofs);
137 
138 
139  // get the references to current and previous sol and velocity
140  const libMesh::NumericVector<Real>
141  &sol_global = *sols[0],
142  &vel_global = *sols[1],
143  &acc_global = *sols[2];
144 
145  for (unsigned int i=0; i<n_dofs; i++) {
146 
147  sol(i) = sol_global(dof_indices[i]);
148  vel(i) = vel_global(dof_indices[i]);
149  accel(i) = acc_global(dof_indices[i]);
150  }
151 
152  _assembly_ops->set_elem_perturbed_solution(sol);
153  _assembly_ops->set_elem_perturbed_velocity(vel);
154  _assembly_ops->set_elem_perturbed_acceleration(accel);
155 }
156 
157 
158 
159 
160 void
162 update_velocity(libMesh::NumericVector<Real>& vec,
163  const libMesh::NumericVector<Real>& sol) {
164 
165  const libMesh::NumericVector<Real>
166  &prev_sol = this->solution(1),
167  &prev_vel = this->velocity(1),
168  &prev_acc = this->acceleration(1);
169 
170  // first calculate the acceleration
171  vec.zero();
172  vec.add( 1., sol);
173  vec.add( -1., prev_sol);
174  vec.scale(gamma/beta/dt);
175  vec.add( 1.-gamma/beta, prev_vel);
176  vec.add((1.-gamma/2./beta)*dt, prev_acc);
177  vec.close();
178 }
179 
180 
181 
182 
183 void
185 update_acceleration(libMesh::NumericVector<Real>& vec,
186  const libMesh::NumericVector<Real>& sol) {
187 
188  const libMesh::NumericVector<Real>
189  &prev_sol = this->solution(1),
190  &prev_vel = this->velocity(1),
191  &prev_acc = this->acceleration(1);
192 
193  // first calculate the acceleration
194  vec.zero();
195  vec.add( 1, sol);
196  vec.add( -1., prev_sol);
197  vec.scale(1./beta/dt/dt);
198  vec.add( -1./beta/dt, prev_vel);
199  vec.add(-(.5-beta)/beta, prev_acc);
200  vec.close();
201 }
202 
203 
204 
205 void
207 update_sensitivity_velocity(libMesh::NumericVector<Real>& vec,
208  const libMesh::NumericVector<Real>& sol) {
209 
210  const libMesh::NumericVector<Real>
211  &prev_sol = this->solution_sensitivity(1),
212  &prev_vel = this->velocity_sensitivity(1),
213  &prev_acc = this->acceleration_sensitivity(1);
214 
215  // first calculate the acceleration
216  vec.zero();
217  vec.add( 1., sol);
218  vec.add( -1., prev_sol);
219  vec.scale(gamma/beta/dt);
220  vec.add( 1.-gamma/beta, prev_vel);
221  vec.add((1.-gamma/2./beta)*dt, prev_acc);
222  vec.close();
223 }
224 
225 
226 
227 void
229 update_sensitivity_acceleration(libMesh::NumericVector<Real>& vec,
230  const libMesh::NumericVector<Real>& sol) {
231 
232  const libMesh::NumericVector<Real>
233  &prev_sol = this->solution_sensitivity(1),
234  &prev_vel = this->velocity_sensitivity(1),
235  &prev_acc = this->acceleration_sensitivity(1);
236 
237  // first calculate the acceleration
238  vec.zero();
239  vec.add( 1, sol);
240  vec.add( -1., prev_sol);
241  vec.scale(1./beta/dt/dt);
242  vec.add( -1./beta/dt, prev_vel);
243  vec.add(-(.5-beta)/beta, prev_acc);
244  vec.close();
245 }
246 
247 
248 
249 
250 void
252 update_delta_velocity(libMesh::NumericVector<Real>& vec,
253  const libMesh::NumericVector<Real>& sol) {
254 
255  // first calculate the acceleration
256  vec.zero();
257  vec.add( gamma/beta/dt, sol);
258  vec.close();
259 }
260 
261 
262 
263 
264 void
266 update_delta_acceleration(libMesh::NumericVector<Real>& vec,
267  const libMesh::NumericVector<Real>& sol) {
268 
269  // first calculate the acceleration
270  vec.zero();
271  vec.add( 1./beta/dt/dt, sol);
272  vec.close();
273 }
274 
275 
276 
277 
278 
279 void
281 elem_calculations(bool if_jac,
282  RealVectorX& vec,
283  RealMatrixX& mat) {
284  // make sure that the assembly object is provided
285  libmesh_assert(_assembly);
286  unsigned int n_dofs = (unsigned int)vec.size();
287 
289  f_x = RealVectorX::Zero(n_dofs),
290  f_m = RealVectorX::Zero(n_dofs);
291 
293  f_m_jac_xddot = RealMatrixX::Zero(n_dofs, n_dofs),
294  f_m_jac_xdot = RealMatrixX::Zero(n_dofs, n_dofs),
295  f_m_jac = RealMatrixX::Zero(n_dofs, n_dofs),
296  f_x_jac_xdot = RealMatrixX::Zero(n_dofs, n_dofs),
297  f_x_jac = RealMatrixX::Zero(n_dofs, n_dofs);
298 
299  // perform the element assembly
300  _assembly_ops->elem_calculations(if_jac,
301  f_m, // mass vector
302  f_x, // forcing vector
303  f_m_jac_xddot, // Jac of mass wrt x_dotdot
304  f_m_jac_xdot, // Jac of mass wrt x_dot
305  f_m_jac, // Jac of mass wrt x
306  f_x_jac_xdot, // Jac of forcing vector wrt x_dot
307  f_x_jac); // Jac of forcing vector wrt x
308 
309 
310  if (_if_highest_derivative_solution) {
311 
312  //
313  // The residual here is modeled as
314  // r(x, xdot, xddot) = f_m(x, xdot, xddot) + f_x(x, xdot)= 0
315  //
316  // then, given x = x0, and xdot = xdot0, the residual can be used to
317  // evaluate xddot at the same time step as
318  //
319  // r(x0, xdot0, xddot) + dr/dxddot dxddot = 0
320  //
321 
322  // system residual
323  vec = (f_m + f_x);
324 
325  // system Jacobian
326  if (if_jac)
327  mat = f_m_jac_xddot;
328  }
329  else {
330 
331  // N-R solver: r(x) = 0, and provides estimates for x (current solution)
332  // in an iterative fashion.
333  // Newmark solver uses x and evaluates x_ddot and x_dot (acc and vel at
334  // current x).
335  //
336  // N-R requires r(x) and dr/dx (residual and Jacobian)
337  // N-R asks Newmark to provide these quantities
338  //
339  // Newmark asks elements (which provide physics) for contributions to
340  // f_m and f_x (using the current estimates of x_ddot, x_dot, x)
341  //
342  //
343  // The residual here is modeled as
344  // r = (f_m + f_x )= 0
345  // where, (for example)
346  // f_m = int_Omega phi u_dot [typical mass vector in conduction, for example]
347  // f_x = int_Omega phi_i u_i - int_Gamma phi q_n [typical conductance and heat flux combination, for example]
348  //
349  // This method assumes
350  // x = x0 + dt x0_dot + (1/2-beta) dt^2 x0_ddot + beta dt^2 x_ddot
351  // or, x_ddot = (x-x0)/beta/dt^2 - 1/beta/dt x0_dot - (1/2-beta)/beta x0_ddot
352  //
353  // x_dot = x0_dot + (1-gamma) dt x0_ddot + gamma dt x_ddot
354  // or, x_dot = x0_dot + (1-gamma) dt x0_ddot +
355  // gamma dt [(x-x0)/beta/dt^2 - 1/beta/dt x0_dot - (1/2-beta)/beta x0_ddot]
356  // = x0_dot + (1-gamma) dt x0_ddot +
357  // gamma/beta/dt (x-x0) - gamma/beta x0_dot - gamma (1/2-beta)/beta dt x0_ddot
358  // = gamma/beta/dt (x-x0) + (1 - gamma/beta) x0_dot + (1 - gamma/2/beta) dt x0_ddot
359  //
360  // Both f_m and f_x can be functions of x_dot and x. Then, the
361  // Jacobian is
362  // dr/dx = df_m/dx + df_x/dx +
363  // df_m/dx_dot dx_dot/dx + df_x/dx_dot dx_dot/dx
364  // = (df_m/dx + df_x/dx) +
365  // (df_m/dx_dot + df_x/dx_dot) (gamma/beta/dt) +
366  // df_m/dx_ddot (1/beta/dt/dt) +
367  //
368 
369 
370  // system residual
371  vec = (f_m + f_x);
372 
373  // system Jacobian
374  if (if_jac)
375  mat = ((1./beta/dt/dt) * f_m_jac_xddot +
376  (gamma/beta/dt)* (f_m_jac_xdot+f_x_jac_xdot) +
377  (f_m_jac + f_x_jac));
378  }
379 }
380 
381 
382 
383 void
386 
387  // make sure that the assembly object is provided
388  libmesh_assert(_assembly_ops);
389 
390  // perform the element assembly
391  _assembly_ops->linearized_jacobian_solution_product(vec);
392 }
393 
394 
395 
396 
397 void
400  RealVectorX& vec) {
401 
402  // make sure that the assembly object is provided
403  libmesh_assert(_assembly_ops);
404  unsigned int n_dofs = (unsigned int)vec.size();
405 
407  f_x = RealVectorX::Zero(n_dofs),
408  f_m = RealVectorX::Zero(n_dofs);
409 
410  // perform the element assembly
411  _assembly_ops->elem_sensitivity_calculations(f,
412  f_m, // mass vector
413  f_x); // forcing vector
414 
415  // system residual
416  vec = (f_m + f_x);
417 }
418 
419 
420 void
422 elem_sensitivity_contribution_previous_timestep(const std::vector<RealVectorX>& prev_sols,
423  RealVectorX& vec) {
424 
425  // make sure that the assembly object is provided
426  libmesh_assert(_assembly_ops);
427  libmesh_assert_equal_to(prev_sols.size(), 3);
428 
429  unsigned int n_dofs = (unsigned int)vec.size();
430 
431  const RealVectorX
432  &sol = prev_sols[0],
433  &vel = prev_sols[1],
434  &accel = prev_sols[2];
436  dummy_vec = RealVectorX::Zero(n_dofs);
437 
439  f_m_jac_xdot = RealMatrixX::Zero(n_dofs, n_dofs),
440  dummy_mat = RealMatrixX::Zero(n_dofs, n_dofs);
441 
442  // perform the element assembly
443  _assembly_ops->elem_calculations(true,
444  dummy_vec, // mass vector
445  dummy_vec, // forcing vector
446  f_m_jac_xdot, // Jac of mass wrt x_dot
447  dummy_mat, // Jac of mass wrt x
448  dummy_mat); // Jac of forcing vector wrt x
449 
450  libmesh_assert(false); // this expression is in error and needs to be set
451  vec -= f_m_jac_xdot * ( (1./beta/dt)*sol + (1.-beta)/beta * vel);
452 }
453 
454 
455 void
458  RealVectorX& vec) {
459 
460  libmesh_assert(false); // to be implemented
461 }
462 
463 
464 
465 void
469  RealVectorX& vec) {
470  libmesh_assert(false); // to be implemented
471 }
472 
473 
474 
virtual void set_element_data(const std::vector< libMesh::dof_id_type > &dof_indices, const std::vector< libMesh::NumericVector< Real > * > &sols)
provides the element with the transient data for calculations
virtual void elem_sensitivity_contribution_previous_timestep(const std::vector< RealVectorX > &prev_sols, RealVectorX &vec)
computes the contribution for this element from previous time step
virtual void update_delta_acceleration(libMesh::NumericVector< Real > &acc, const libMesh::NumericVector< Real > &sol)
update the perturbation in transient acceleration based on the current perturbed solution ...
virtual void update_acceleration(libMesh::NumericVector< Real > &vec, const libMesh::NumericVector< Real > &sol)
update the transient acceleration based on the current solution
virtual void update_sensitivity_velocity(libMesh::NumericVector< Real > &vel, const libMesh::NumericVector< Real > &sol)
update the transient sensitivity velocity based on the current sensitivity solution ...
virtual void elem_calculations(bool if_jac, RealVectorX &vec, RealMatrixX &mat)
performs the element calculations over elem, and returns the element vector and matrix quantities in ...
virtual void elem_sensitivity_calculations(const MAST::FunctionBase &f, RealVectorX &vec)=0
performs the element sensitivity calculations over elem, and returns the element residual sensitivity...
libMesh::NumericVector< Real > & solution_sensitivity(unsigned int prev_iter=0) const
libMesh::NumericVector< Real > & velocity(unsigned int prev_iter=0) const
virtual void update_delta_velocity(libMesh::NumericVector< Real > &vel, const libMesh::NumericVector< Real > &sol)
update the perturbation in transient velocity based on the current perturbed solution ...
Matrix< Real, Dynamic, Dynamic > RealMatrixX
virtual void elem_shape_sensitivity_calculations(const MAST::FunctionBase &f, RealVectorX &vec)=0
performs the element shape sensitivity calculations over elem, and returns the element residual sensi...
libMesh::NumericVector< Real > & velocity_sensitivity(unsigned int prev_iter=0) const
Matrix< Real, Dynamic, 1 > RealVectorX
virtual void elem_topology_sensitivity_calculations(const MAST::FunctionBase &f, const MAST::FieldFunction< RealVectorX > &vel, RealVectorX &vec)=0
performs the element topology sensitivity calculations over elem, and returns the element residual se...
virtual void update_sensitivity_acceleration(libMesh::NumericVector< Real > &acc, const libMesh::NumericVector< Real > &sol)
update the transient sensitivity acceleration based on the current sensitivity solution ...
libMesh::NumericVector< Real > & acceleration_sensitivity(unsigned int prev_iter=0) const
libMesh::NumericVector< Real > & solution(unsigned int prev_iter=0) const
virtual void update_velocity(libMesh::NumericVector< Real > &vec, const libMesh::NumericVector< Real > &sol)
update the transient velocity based on the current solution
virtual void set_element_perturbed_data(const std::vector< libMesh::dof_id_type > &dof_indices, const std::vector< libMesh::NumericVector< Real > * > &sols)
provides the element with the transient data for calculations
virtual void extract_element_sensitivity_data(const std::vector< libMesh::dof_id_type > &dof_indices, const std::vector< libMesh::NumericVector< Real > * > &sols, std::vector< RealVectorX > &local_sols)
provides the element with the sensitivity of transient data for calculations
virtual void elem_linearized_jacobian_solution_product(RealVectorX &vec)=0
performs the element calculations over elem, and returns the element vector quantity in vec...
libMesh::NumericVector< Real > & acceleration(unsigned int prev_iter=0) const