Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #include "libmesh/newmark_solver.h"
19 :
20 : #include "libmesh/diff_solver.h"
21 : #include "libmesh/diff_system.h"
22 : #include "libmesh/dof_map.h"
23 : #include "libmesh/numeric_vector.h"
24 :
25 : namespace libMesh
26 : {
27 426 : NewmarkSolver::NewmarkSolver (sys_type & s)
28 : : SecondOrderUnsteadySolver(s),
29 402 : _beta(0.25),
30 402 : _gamma(0.5),
31 402 : _is_accel_solve(false),
32 426 : _initial_accel_set(false)
33 426 : {}
34 :
35 804 : NewmarkSolver::~NewmarkSolver () = default;
36 :
37 0 : Real NewmarkSolver::error_order() const
38 : {
39 0 : if (_gamma == 0.5)
40 0 : return 2.;
41 0 : return 1.;
42 : }
43 :
44 3976 : void NewmarkSolver::advance_timestep ()
45 : {
46 : // We need to update velocity and acceleration before
47 : // we update the nonlinear solution (displacement) and
48 : // delta_t
49 :
50 : NumericVector<Number> & old_solution_rate =
51 3976 : _system.get_vector("_old_solution_rate");
52 :
53 : NumericVector<Number> & old_solution_accel =
54 3976 : _system.get_vector("_old_solution_accel");
55 :
56 3976 : if (!first_solve)
57 : {
58 : NumericVector<Number> & old_nonlinear_soln =
59 3550 : _system.get_vector("_old_nonlinear_solution");
60 :
61 : NumericVector<Number> & nonlinear_solution =
62 3550 : *(_system.solution);
63 :
64 : // We need to cache the new solution_rate before updating the old_solution_rate
65 : // so we can update acceleration with the proper old_solution_rate
66 : // v_{n+1} = gamma/(beta*Delta t)*(x_{n+1}-x_n)
67 : // - ((gamma/beta)-1)*v_n
68 : // - (gamma/(2*beta)-1)*(Delta t)*a_n
69 3650 : std::unique_ptr<NumericVector<Number>> new_solution_rate = nonlinear_solution.clone();
70 3550 : (*new_solution_rate) -= old_nonlinear_soln;
71 3550 : (*new_solution_rate) *= (_gamma/(_beta*_system.deltat));
72 3550 : new_solution_rate->add( (1.0-_gamma/_beta), old_solution_rate );
73 3550 : new_solution_rate->add( (1.0-_gamma/(2.0*_beta))*_system.deltat, old_solution_accel );
74 :
75 : // a_{n+1} = (1/(beta*(Delta t)^2))*(x_{n+1}-x_n)
76 : // - 1/(beta*Delta t)*v_n
77 : // - (1-1/(2*beta))*a_n
78 3650 : std::unique_ptr<NumericVector<Number>> new_solution_accel = old_solution_accel.clone();
79 3550 : (*new_solution_accel) *= -(1.0/(2.0*_beta)-1.0);
80 3550 : new_solution_accel->add( -1.0/(_beta*_system.deltat), old_solution_rate );
81 3550 : new_solution_accel->add( 1.0/(_beta*_system.deltat*_system.deltat), nonlinear_solution );
82 3550 : new_solution_accel->add( -1.0/(_beta*_system.deltat*_system.deltat), old_nonlinear_soln );
83 :
84 : // Now update old_solution_rate
85 3650 : old_solution_rate = (*new_solution_rate);
86 3650 : old_solution_accel = (*new_solution_accel);
87 3350 : }
88 :
89 : // Localize updated vectors
90 : old_solution_rate.localize
91 3976 : (*_old_local_solution_rate,
92 3976 : _system.get_dof_map().get_send_list());
93 :
94 : old_solution_accel.localize
95 3976 : (*_old_local_solution_accel,
96 3976 : _system.get_dof_map().get_send_list());
97 :
98 : // Now we can finish advancing the timestep
99 3976 : UnsteadySolver::advance_timestep();
100 3976 : }
101 :
102 0 : void NewmarkSolver::adjoint_advance_timestep ()
103 : {
104 0 : libmesh_not_implemented();
105 : }
106 :
107 426 : void NewmarkSolver::compute_initial_accel()
108 : {
109 : // We need to compute the initial acceleration based off of
110 : // the initial position and velocity and, thus, acceleration
111 : // is the unknown in diff_solver and not the displacement. So,
112 : // We swap solution and acceleration. NewmarkSolver::_general_residual
113 : // will check _is_accel_solve and provide the correct
114 : // values to the FEMContext assuming this swap was made.
115 426 : this->_is_accel_solve = true;
116 :
117 : //solution->accel, accel->solution
118 438 : _system.solution->swap(_system.get_vector("_old_solution_accel"));
119 426 : _system.update();
120 :
121 426 : this->_diff_solver->solve();
122 :
123 : // solution->solution, accel->accel
124 438 : _system.solution->swap(_system.get_vector("_old_solution_accel"));
125 426 : _system.update();
126 :
127 : // We're done, so no longer doing an acceleration solve
128 426 : this->_is_accel_solve = false;
129 :
130 426 : this->set_initial_accel_avail(true);
131 426 : }
132 :
133 0 : void NewmarkSolver::project_initial_accel(FunctionBase<Number> * f,
134 : FunctionBase<Gradient> * g)
135 : {
136 : NumericVector<Number> & old_solution_accel =
137 0 : _system.get_vector("_old_solution_accel");
138 :
139 0 : _system.project_vector( old_solution_accel, f, g );
140 :
141 0 : this->set_initial_accel_avail(true);
142 0 : }
143 :
144 426 : void NewmarkSolver::set_initial_accel_avail( bool initial_accel_set )
145 : {
146 426 : _initial_accel_set = initial_accel_set;
147 426 : }
148 :
149 3550 : void NewmarkSolver::solve ()
150 : {
151 : // First, check that the initial accel was set one way or another
152 3550 : libmesh_error_msg_if(!_initial_accel_set,
153 : "ERROR: Must first set initial acceleration using one of:\n"
154 : "NewmarkSolver::compute_initial_accel()\n"
155 : "NewmarkSolver::project_initial_accel()\n");
156 :
157 : // That satisfied, now we can solve
158 3550 : UnsteadySolver::solve();
159 3550 : }
160 :
161 973808 : bool NewmarkSolver::element_residual (bool request_jacobian,
162 : DiffContext & context)
163 : {
164 1055152 : return this->_general_residual(request_jacobian,
165 : context,
166 : &DifferentiablePhysics::mass_residual,
167 : &DifferentiablePhysics::damping_residual,
168 : &DifferentiablePhysics::_eulerian_time_deriv,
169 : &DifferentiablePhysics::element_constraint,
170 973808 : &DiffContext::elem_reinit);
171 : }
172 :
173 :
174 :
175 638720 : bool NewmarkSolver::side_residual (bool request_jacobian,
176 : DiffContext & context)
177 : {
178 691968 : return this->_general_residual(request_jacobian,
179 : context,
180 : &DifferentiablePhysics::side_mass_residual,
181 : &DifferentiablePhysics::side_damping_residual,
182 : &DifferentiablePhysics::side_time_derivative,
183 : &DifferentiablePhysics::side_constraint,
184 638720 : &DiffContext::elem_side_reinit);
185 : }
186 :
187 :
188 :
189 0 : bool NewmarkSolver::nonlocal_residual (bool request_jacobian,
190 : DiffContext & context)
191 : {
192 0 : return this->_general_residual(request_jacobian,
193 : context,
194 : &DifferentiablePhysics::nonlocal_mass_residual,
195 : &DifferentiablePhysics::nonlocal_damping_residual,
196 : &DifferentiablePhysics::nonlocal_time_derivative,
197 : &DifferentiablePhysics::nonlocal_constraint,
198 0 : &DiffContext::nonlocal_reinit);
199 : }
200 :
201 :
202 :
203 1612528 : bool NewmarkSolver::_general_residual (bool request_jacobian,
204 : DiffContext & context,
205 : ResFuncType mass,
206 : ResFuncType damping,
207 : ResFuncType time_deriv,
208 : ResFuncType constraint,
209 : ReinitFuncType reinit_func)
210 : {
211 134592 : unsigned int n_dofs = context.get_elem_solution().size();
212 :
213 : // We might need to save the old jacobian in case one of our physics
214 : // terms later is unable to update it analytically.
215 1881712 : DenseMatrix<Number> old_elem_jacobian(n_dofs, n_dofs);
216 :
217 : // Local velocity at old time step
218 1477936 : DenseVector<Number> old_elem_solution_rate(n_dofs);
219 37409680 : for (unsigned int i=0; i != n_dofs; ++i)
220 35797152 : old_elem_solution_rate(i) =
221 38782384 : old_solution_rate(context.get_dof_indices()[i]);
222 :
223 : // The user is computing the initial acceleration
224 : // So upstream we've swapped _system.solution and _old_local_solution_accel
225 : // So we need to give the context the correct entries since we're solving for
226 : // acceleration here.
227 1612528 : if (_is_accel_solve)
228 : {
229 : // System._solution is actually the acceleration right now so we need
230 : // to reset the elem_solution to the right thing, which in this case
231 : // is the initial guess for displacement, which got swapped into
232 : // _old_solution_accel vector
233 1176296 : DenseVector<Number> old_elem_solution(n_dofs);
234 30683232 : for (unsigned int i=0; i != n_dofs; ++i)
235 29400000 : old_elem_solution(i) =
236 31850000 : old_solution_accel(context.get_dof_indices()[i]);
237 :
238 1283232 : context.elem_solution_derivative = 0.0;
239 1283232 : context.elem_solution_rate_derivative = 0.0;
240 1283232 : context.elem_solution_accel_derivative = 1.0;
241 :
242 : // Acceleration is currently the unknown so it's already sitting
243 : // in elem_solution() thanks to FEMContext::pre_fe_reinit
244 106936 : context.get_elem_solution_accel() = context.get_elem_solution();
245 :
246 : // Now reset elem_solution() to what the user is expecting
247 106936 : context.get_elem_solution() = old_elem_solution;
248 :
249 106936 : context.get_elem_solution_rate() = old_elem_solution_rate;
250 :
251 : // The user's Jacobians will be targeting derivatives w.r.t. u_{n+1}.
252 : // Although the vast majority of cases will have the correct analytic
253 : // Jacobians in this iteration, since we reset elem_solution_derivative*,
254 : // if there are coupled/overlapping problems, there could be
255 : // mismatches in the Jacobian. So we force finite differencing for
256 : // the first iteration.
257 106936 : request_jacobian = false;
258 : }
259 : // Otherwise, the unknowns are the displacements and everything is straight
260 : // forward and is what you think it is
261 : else
262 : {
263 329296 : if (request_jacobian)
264 150820 : old_elem_jacobian.swap(context.get_elem_jacobian());
265 :
266 : // Local displacement at old timestep
267 329296 : DenseVector<Number> old_elem_solution(n_dofs);
268 6726448 : for (unsigned int i=0; i != n_dofs; ++i)
269 6397152 : old_elem_solution(i) =
270 6932384 : old_nonlinear_solution(context.get_dof_indices()[i]);
271 :
272 : // Local acceleration at old time step
273 329296 : DenseVector<Number> old_elem_solution_accel(n_dofs);
274 6726448 : for (unsigned int i=0; i != n_dofs; ++i)
275 6397152 : old_elem_solution_accel(i) =
276 6932384 : old_solution_accel(context.get_dof_indices()[i]);
277 :
278 : // Convenience
279 329296 : libMesh::Real dt = _system.deltat;
280 :
281 329296 : context.elem_solution_derivative = 1.0;
282 :
283 : // Local velocity at current time step
284 : // v_{n+1} = gamma/(beta*Delta t)*(x_{n+1}-x_n)
285 : // + (1-(gamma/beta))*v_n
286 : // + (1-gamma/(2*beta))*(Delta t)*a_n
287 329296 : context.elem_solution_rate_derivative = (_gamma/(_beta*dt));
288 :
289 27656 : context.get_elem_solution_rate() = context.get_elem_solution();
290 301640 : context.get_elem_solution_rate() -= old_elem_solution;
291 329296 : context.get_elem_solution_rate() *= context.elem_solution_rate_derivative;
292 329296 : context.get_elem_solution_rate().add( (1.0-_gamma/_beta), old_elem_solution_rate);
293 329296 : context.get_elem_solution_rate().add( (1.0-_gamma/(2.0*_beta))*dt, old_elem_solution_accel);
294 :
295 :
296 :
297 : // Local acceleration at current time step
298 : // a_{n+1} = (1/(beta*(Delta t)^2))*(x_{n+1}-x_n)
299 : // - 1/(beta*Delta t)*v_n
300 : // - (1/(2*beta)-1)*a_n
301 329296 : context.elem_solution_accel_derivative = 1.0/(_beta*dt*dt);
302 :
303 27656 : context.get_elem_solution_accel() = context.get_elem_solution();
304 301640 : context.get_elem_solution_accel() -= old_elem_solution;
305 329296 : context.get_elem_solution_accel() *= context.elem_solution_accel_derivative;
306 329296 : context.get_elem_solution_accel().add(-1.0/(_beta*dt), old_elem_solution_rate);
307 329296 : context.get_elem_solution_accel().add(-(1.0/(2.0*_beta)-1.0), old_elem_solution_accel);
308 :
309 : // Move the mesh into place first if necessary, set t = t_{n+1}
310 329296 : (context.*reinit_func)(1.);
311 : }
312 :
313 : // If a fixed solution is requested, we'll use x_{n+1}
314 1612528 : if (_system.use_fixed_solution)
315 0 : context.get_elem_fixed_solution() = context.get_elem_solution();
316 :
317 : // Get the time derivative at t_{n+1}, F(u_{n+1})
318 3090464 : bool jacobian_computed = (_system.get_physics()->*time_deriv)(request_jacobian, context);
319 :
320 : // Damping at t_{n+1}, C(u_{n+1})
321 3090464 : jacobian_computed = (_system.get_physics()->*damping)(jacobian_computed, context) &&
322 : jacobian_computed;
323 :
324 : // Mass at t_{n+1}, M(u_{n+1})
325 3090464 : jacobian_computed = (_system.get_physics()->*mass)(jacobian_computed, context) &&
326 : jacobian_computed;
327 :
328 : // Add the constraint term
329 3090464 : jacobian_computed = (_system.get_physics()->*constraint)(jacobian_computed, context) &&
330 : jacobian_computed;
331 :
332 : // Add back (or restore) the old jacobian
333 1612528 : if (request_jacobian)
334 : {
335 164648 : if (jacobian_computed)
336 150820 : context.get_elem_jacobian() += old_elem_jacobian;
337 : else
338 0 : context.get_elem_jacobian().swap(old_elem_jacobian);
339 : }
340 :
341 1747120 : return jacobian_computed;
342 1343344 : }
343 :
344 : }
|
