MAST
generalized_alpha_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 rho_inf(0.),
35 alpha_m(0.),
36 alpha_f(0.) {
37 
38  this->update_coefficient(0.2);
39 }
40 
41 
43 { }
44 
45 
46 
47 void
49 
50  this->rho_inf = rho_infinity;
51 
52  // values of alpha_m and alpha_f based on rho_inf
53  this->alpha_m = (2.*rho_inf - 1.)/(rho_inf + 1.);
54  this->alpha_f = rho_inf/(rho_inf + 1.);
55 
56  // for second order accuracy
57  this->gamma = 0.5 - alpha_m + alpha_f;
58 
59  // for maximizing high frequency dissipation
60  this->beta = .25 * std::pow(1. - alpha_m + alpha_f, 2);
61 }
62 
63 
64 
65 void
67 set_element_data(const std::vector<libMesh::dof_id_type>& dof_indices,
68  const std::vector<libMesh::NumericVector<Real>*>& sols){
69 
70  libmesh_assert_equal_to(sols.size(), 3);
71 
72  const unsigned int n_dofs = (unsigned int)dof_indices.size();
73 
74  // get the current state and velocity estimates
75  // also get the current discrete velocity replacement
77  sol = RealVectorX::Zero(n_dofs),
78  vel = RealVectorX::Zero(n_dofs),
79  accel = RealVectorX::Zero(n_dofs);
80 
81 
82  // get the references to current and previous sol and velocity
83  const libMesh::NumericVector<Real>
84  &sol_global = *sols[0],
85  &vel_global = *sols[1],
86  &acc_global = *sols[2],
87  &prev_sol = this->solution(1),
88  &prev_vel = this->velocity(1),
89  &prev_acc = this->acceleration(1);
90 
91 
92  for (unsigned int i=0; i<n_dofs; i++) {
93 
94  sol(i) = ((1.-alpha_f)*sol_global(dof_indices[i]) +
95  alpha_f*prev_sol(dof_indices[i]));
96  vel(i) = ((1.-alpha_f)*vel_global(dof_indices[i]) +
97  alpha_f*prev_vel(dof_indices[i]));
98  accel(i) = ((1.-alpha_m)*acc_global(dof_indices[i]) +
99  alpha_m*prev_acc(dof_indices[i]));
100  }
101 
102  _assembly_ops->set_elem_solution(sol);
103  _assembly_ops->set_elem_velocity(vel);
104  _assembly_ops->set_elem_acceleration(accel);
105 }
106 
107 
108 
109 void
111 extract_element_sensitivity_data(const std::vector<libMesh::dof_id_type>& dof_indices,
112  const std::vector<libMesh::NumericVector<Real>*>& sols,
113  std::vector<RealVectorX>& local_sols) {
114 
115  libmesh_assert_equal_to(sols.size(), 3);
116 
117  const unsigned int n_dofs = (unsigned int)dof_indices.size();
118 
119  local_sols.resize(3);
120 
122  &sol = local_sols[0],
123  &vel = local_sols[1],
124  &accel = local_sols[2];
125 
126  sol = RealVectorX::Zero(n_dofs);
127  vel = RealVectorX::Zero(n_dofs);
128  accel = RealVectorX::Zero(n_dofs);
129 
130 
131  // get the references to current and previous sol and velocity
132  const libMesh::NumericVector<Real>
133  &sol_global = *sols[0],
134  &vel_global = *sols[1],
135  &acc_global = *sols[2],
136  &prev_sol = this->solution_sensitivity(1),
137  &prev_vel = this->velocity_sensitivity(1),
138  &prev_acc = this->acceleration_sensitivity(1);
139 
140  for (unsigned int i=0; i<n_dofs; i++) {
141 
142  sol(i) = ((1.-alpha_f)*sol_global(dof_indices[i]) +
143  alpha_f*prev_sol(dof_indices[i]));
144  vel(i) = ((1.-alpha_f)*vel_global(dof_indices[i]) +
145  alpha_f*prev_vel(dof_indices[i]));
146  accel(i) = ((1.-alpha_m)*acc_global(dof_indices[i]) +
147  alpha_m*prev_acc(dof_indices[i]));
148  }
149 
150 
151 }
152 
153 
154 
155 void
157 set_element_perturbed_data(const std::vector<libMesh::dof_id_type>& dof_indices,
158  const std::vector<libMesh::NumericVector<Real>*>& sols){
159 
160  libmesh_assert_equal_to(sols.size(), 3);
161 
162  const unsigned int n_dofs = (unsigned int)dof_indices.size();
163 
164  // get the current state and velocity estimates
165  // also get the current discrete velocity replacement
167  sol = RealVectorX::Zero(n_dofs),
168  vel = RealVectorX::Zero(n_dofs),
169  accel = RealVectorX::Zero(n_dofs);
170 
171 
172  // get the references to current and previous sol and velocity
173  const libMesh::NumericVector<Real>
174  &sol_global = *sols[0],
175  &vel_global = *sols[1],
176  &acc_global = *sols[2];
177 
178  for (unsigned int i=0; i<n_dofs; i++) {
179 
180  sol(i) = sol_global(dof_indices[i]);
181  vel(i) = vel_global(dof_indices[i]);
182  accel(i) = acc_global(dof_indices[i]);
183  }
184 
185  _assembly_ops->set_elem_perturbed_solution(sol);
186  _assembly_ops->set_elem_perturbed_velocity(vel);
187  _assembly_ops->set_elem_perturbed_acceleration(accel);
188 }
189 
190 
191 
192 
193 void
195 elem_calculations(bool if_jac,
196  RealVectorX& vec,
197  RealMatrixX& mat) {
198  // make sure that the assembly object is provided
199  libmesh_assert(_assembly);
200  unsigned int n_dofs = (unsigned int)vec.size();
201 
203  f_x = RealVectorX::Zero(n_dofs),
204  f_m = RealVectorX::Zero(n_dofs);
205 
207  f_m_jac_xddot = RealMatrixX::Zero(n_dofs, n_dofs),
208  f_m_jac_xdot = RealMatrixX::Zero(n_dofs, n_dofs),
209  f_m_jac = RealMatrixX::Zero(n_dofs, n_dofs),
210  f_x_jac_xdot = RealMatrixX::Zero(n_dofs, n_dofs),
211  f_x_jac = RealMatrixX::Zero(n_dofs, n_dofs);
212 
213  // perform the element assembly
214  _assembly_ops->elem_calculations(if_jac,
215  f_m, // mass vector
216  f_x, // forcing vector
217  f_m_jac_xddot, // Jac of mass wrt x_dotdot
218  f_m_jac_xdot, // Jac of mass wrt x_dot
219  f_m_jac, // Jac of mass wrt x
220  f_x_jac_xdot, // Jac of forcing vector wrt x_dot
221  f_x_jac); // Jac of forcing vector wrt x
222 
223 
224  if (_if_highest_derivative_solution) {
225 
226  //
227  // The residual here is modeled as
228  // r(x, xdot, xddot) = f_m(\tilde{x}, \tilde{xdot}, \tilde{xddot}) +
229  // f_x(\tilde{x}, \tilde{xdot})= 0
230  // \tilde{x} = (1-a_f) x1 + a_f x0
231  // \tilde{xdot} = (1-a_f) xdot1 + a_f xdot0
232  // \tilde{xddot} = (1-a_m) xddot1 + a_m xddot0
233  //
234  // then, given x = x0, and xdot = xdot0, the residual can be used to
235  // evaluate xddot at the same time step as
236  //
237  // r(x0, xdot0, xddot) + dr/dxddot dxddot = 0
238  //
239 
240  // system residual
241  vec = (f_m + f_x);
242 
243  // system Jacobian
244  if (if_jac)
245  mat = f_m_jac_xddot;
246  }
247  else {
248 
249  // N-R solver: r(x) = 0, and provides estimates for x (current solution)
250  // in an iterative fashion.
251  // Newmark solver uses x and evaluates x_ddot and x_dot (acc and vel at
252  // current x).
253  //
254  // N-R requires r(x) and dr/dx (residual and Jacobian)
255  // N-R asks Newmark to provide these quantities
256  //
257  // Newmark asks elements (which provide physics) for contributions to
258  // f_m and f_x (using the current estimates of x_ddot, x_dot, x)
259  //
260  //
261  // The residual here is modeled as
262  // r = (f_m + f_x )= 0
263  // where, (for example)
264  // f_m = int_Omega phi u_dot [typical mass vector in conduction, for example]
265  // f_x = int_Omega phi_i u_i - int_Gamma phi q_n [typical conductance and heat flux combination, for example]
266  //
267  // This method assumes
268  // x = x0 + dt x0_dot + (1/2-beta) dt^2 x0_ddot + beta dt^2 x_ddot
269  // or, x_ddot = (x-x0)/beta/dt^2 - 1/beta/dt x0_dot - (1/2-beta)/beta x0_ddot
270  //
271  // x_dot = x0_dot + (1-gamma) dt x0_ddot + gamma dt x_ddot
272  // or, x_dot = x0_dot + (1-gamma) dt x0_ddot +
273  // gamma dt [(x-x0)/beta/dt^2 - 1/beta/dt x0_dot - (1/2-beta)/beta x0_ddot]
274  // = x0_dot + (1-gamma) dt x0_ddot +
275  // gamma/beta/dt (x-x0) - gamma/beta x0_dot - gamma (1/2-beta)/beta dt x0_ddot
276  // = gamma/beta/dt (x-x0) + (1 - gamma/beta) x0_dot + (1 - gamma/2/beta) dt x0_ddot
277  //
278  // Both f_m and f_x can be functions of x_dot and x. Then, the
279  // Jacobian is
280  // dr/dx = (df_m/dx + df_x/dx) (1-a_f) +
281  // df_m/dx_dot dx_dot/dx (1-a_f) +
282  // df_x/dx_dot dx_dot/dx (1-a_m)
283  // = (df_m/dx + df_x/dx) (1-a_f) +
284  // (df_m/dx_dot + df_x/dx_dot) (gamma/beta/dt) (1-a_f) +
285  // df_m/dx_ddot (1/beta/dt/dt) (1-a_m) +
286  //
287 
288 
289  // system residual
290  vec = (f_m + f_x);
291 
292  // system Jacobian
293  if (if_jac)
294  mat = ((1./beta/dt/dt) * (1.-alpha_m) * f_m_jac_xddot +
295  (gamma/beta/dt) * (1.-alpha_f) * (f_m_jac_xdot+f_x_jac_xdot) +
296  (1.-alpha_f) * (f_m_jac + f_x_jac));
297  }
298 }
299 
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 ...
This class implements the Newmark solver for solution of a second-order ODE.
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
void update_coefficient(Real rho_infinity)
Computes the value of coefficients basde on value of .
libMesh::Real Real
libMesh::NumericVector< Real > & solution_sensitivity(unsigned int prev_iter=0) const
libMesh::NumericVector< Real > & velocity(unsigned int prev_iter=0) const
Matrix< Real, Dynamic, Dynamic > RealMatrixX
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
libMesh::NumericVector< Real > & velocity_sensitivity(unsigned int prev_iter=0) const
Matrix< Real, Dynamic, 1 > RealVectorX
libMesh::NumericVector< Real > & acceleration_sensitivity(unsigned int prev_iter=0) const
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
libMesh::NumericVector< Real > & solution(unsigned int prev_iter=0) const
libMesh::NumericVector< Real > & acceleration(unsigned int prev_iter=0) const