1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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
19
20 #include "libmesh/libmesh_config.h"
21 #ifdef LIBMESH_HAVE_SLEPC
22
23 #include "libmesh/diff_system.h"
24 #include "libmesh/eigen_time_solver.h"
25 #include "libmesh/eigen_solver.h"
26 #include "libmesh/sparse_matrix.h"
27 #include "libmesh/enum_eigen_solver_type.h"
28
29 namespace libMesh
30 {
31
32 // Constructor
EigenTimeSolver(sys_type & s)33 EigenTimeSolver::EigenTimeSolver (sys_type & s)
34 : Parent(s),
35 eigen_solver (EigenSolver<Number>::build(s.comm())),
36 tol(pow(TOLERANCE, 5./3.)),
37 maxits(1000),
38 n_eigenpairs_to_compute(5),
39 n_basis_vectors_to_use(3*n_eigenpairs_to_compute),
40 n_converged_eigenpairs(0),
41 n_iterations_reqd(0)
42 {
43 libmesh_experimental();
44 eigen_solver->set_eigenproblem_type(GHEP);//or GNHEP
45 eigen_solver->set_position_of_spectrum(LARGEST_MAGNITUDE);
46 }
47
~EigenTimeSolver()48 EigenTimeSolver::~EigenTimeSolver ()
49 {
50 }
51
reinit()52 void EigenTimeSolver::reinit ()
53 {
54 // empty...
55 }
56
init()57 void EigenTimeSolver::init ()
58 {
59 // Add matrix "B" to _system if not already there.
60 // The user may have already added a matrix "B" before
61 // calling the System initialization. This would be
62 // necessary if e.g. the System originally started life
63 // with a different type of TimeSolver and only later
64 // had its TimeSolver changed to an EigenTimeSolver.
65 if (!_system.have_matrix("B"))
66 _system.add_matrix("B");
67 }
68
solve()69 void EigenTimeSolver::solve ()
70 {
71 // The standard implementation is basically to call:
72 // _diff_solver->solve();
73 // which internally assembles (when necessary) matrices and vectors
74 // and calls linear solver software while also doing Newton steps (see newton_solver.C)
75 //
76 // The element_residual and side_residual functions below control
77 // what happens in the interior of the element assembly loops.
78 // We have a system reference, so it's possible to call _system.assembly()
79 // ourselves if we want to...
80 //
81 // Interestingly, for the EigenSolver we don't need residuals...just Jacobians.
82 // The Jacobian should therefore always be requested, and always return
83 // jacobian_computed as being true.
84
85 // The basic plan of attack is:
86 // .) Construct the Jacobian using _system.assembly(true,true) as we
87 // would for a steady system. Use a flag in this class to
88 // control behavior in element_residual and side_residual
89 // .) Swap _system.matrix to matrix "B" (be sure to add this extra matrix during init)
90 // .) Call _system.assembly(true,true) again, using the flag in element_residual
91 // and side_residual to only get the mass matrix terms.
92 // .) Send A and B to Steffen's EigenSolver interface.
93
94 // Assemble the spatial part (matrix A) of the operator
95 if (!this->quiet)
96 libMesh::out << "Assembling matrix A." << std::endl;
97 _system.matrix = &( _system.get_matrix ("System Matrix") );
98 this->now_assembling = Matrix_A;
99 _system.assembly(true, true);
100 //_system.matrix->print_matlab("matrix_A.m");
101
102 // Point the system's matrix at B, call assembly again.
103 if (!this->quiet)
104 libMesh::out << "Assembling matrix B." << std::endl;
105 _system.matrix = &( _system.get_matrix ("B") );
106 this->now_assembling = Matrix_B;
107 _system.assembly(true, true);
108 //_system.matrix->print_matlab("matrix_B.m");
109
110 // Send matrices A, B to Steffen's SlepcEigenSolver interface
111 //libmesh_here();
112 if (!this->quiet)
113 libMesh::out << "Calling the EigenSolver." << std::endl;
114 std::pair<unsigned int, unsigned int> solve_data =
115 eigen_solver->solve_generalized (_system.get_matrix ("System Matrix"),
116 _system.get_matrix ("B"),
117 n_eigenpairs_to_compute,
118 n_basis_vectors_to_use,
119 tol,
120 maxits);
121
122 this->n_converged_eigenpairs = solve_data.first;
123 this->n_iterations_reqd = solve_data.second;
124 }
125
126
127
element_residual(bool request_jacobian,DiffContext & context)128 bool EigenTimeSolver::element_residual(bool request_jacobian,
129 DiffContext & context)
130 {
131 // The EigenTimeSolver always computes jacobians!
132 libmesh_assert (request_jacobian);
133
134 context.elem_solution_rate_derivative = 1;
135 context.elem_solution_derivative = 1;
136
137 // Assemble the operator for the spatial part.
138 if (now_assembling == Matrix_A)
139 {
140 bool jacobian_computed =
141 _system.get_physics()->element_time_derivative(request_jacobian, context);
142
143 // The user shouldn't compute a jacobian unless requested
144 libmesh_assert(request_jacobian || !jacobian_computed);
145
146 bool jacobian_computed2 =
147 _system.get_physics()->element_constraint(jacobian_computed, context);
148
149 // The user shouldn't compute a jacobian unless requested
150 libmesh_assert (jacobian_computed || !jacobian_computed2);
151
152 return jacobian_computed && jacobian_computed2;
153
154 }
155
156 // Assemble the mass matrix operator
157 else if (now_assembling == Matrix_B)
158 {
159 bool mass_jacobian_computed =
160 _system.get_physics()->mass_residual(request_jacobian, context);
161
162 // Scale Jacobian by -1 to get positive matrix from new negative
163 // mass_residual convention
164 context.get_elem_jacobian() *= -1.0;
165
166 return mass_jacobian_computed;
167 }
168
169 else
170 libmesh_error_msg("Unrecognized value now_assembling = " << now_assembling);
171 }
172
173
174
side_residual(bool request_jacobian,DiffContext & context)175 bool EigenTimeSolver::side_residual(bool request_jacobian,
176 DiffContext & context)
177 {
178 // The EigenTimeSolver always requests jacobians?
179 //libmesh_assert (request_jacobian);
180
181 context.elem_solution_rate_derivative = 1;
182 context.elem_solution_derivative = 1;
183
184 // Assemble the operator for the spatial part.
185 if (now_assembling == Matrix_A)
186 {
187 bool jacobian_computed =
188 _system.get_physics()->side_time_derivative(request_jacobian, context);
189
190 // The user shouldn't compute a jacobian unless requested
191 libmesh_assert (request_jacobian || !jacobian_computed);
192
193 bool jacobian_computed2 =
194 _system.get_physics()->side_constraint(jacobian_computed, context);
195
196 // The user shouldn't compute a jacobian unless requested
197 libmesh_assert (jacobian_computed || !jacobian_computed2);
198
199 return jacobian_computed && jacobian_computed2;
200
201 }
202
203 // There is now a "side" equivalent for the mass matrix
204 else if (now_assembling == Matrix_B)
205 {
206 bool mass_jacobian_computed =
207 _system.get_physics()->side_mass_residual(request_jacobian, context);
208
209 // Scale Jacobian by -1 to get positive matrix from new negative
210 // mass_residual convention
211 context.get_elem_jacobian() *= -1.0;
212
213 return mass_jacobian_computed;
214 }
215
216 else
217 libmesh_error_msg("Unrecognized value now_assembling = " << now_assembling);
218 }
219
220
221
nonlocal_residual(bool request_jacobian,DiffContext & context)222 bool EigenTimeSolver::nonlocal_residual(bool request_jacobian,
223 DiffContext & context)
224 {
225 // The EigenTimeSolver always requests jacobians?
226 //libmesh_assert (request_jacobian);
227
228 // Assemble the operator for the spatial part.
229 if (now_assembling == Matrix_A)
230 {
231 bool jacobian_computed =
232 _system.get_physics()->nonlocal_time_derivative(request_jacobian, context);
233
234 // The user shouldn't compute a jacobian unless requested
235 libmesh_assert (request_jacobian || !jacobian_computed);
236
237 bool jacobian_computed2 =
238 _system.get_physics()->nonlocal_constraint(jacobian_computed, context);
239
240 // The user shouldn't compute a jacobian unless requested
241 libmesh_assert (jacobian_computed || !jacobian_computed2);
242
243 return jacobian_computed && jacobian_computed2;
244
245 }
246
247 // There is now a "side" equivalent for the mass matrix
248 else if (now_assembling == Matrix_B)
249 {
250 bool mass_jacobian_computed =
251 _system.get_physics()->nonlocal_mass_residual(request_jacobian, context);
252
253 // Scale Jacobian by -1 to get positive matrix from new negative
254 // mass_residual convention
255 context.get_elem_jacobian() *= -1.0;
256
257 return mass_jacobian_computed;
258 }
259
260 else
261 libmesh_error_msg("Unrecognized value now_assembling = " << now_assembling);
262 }
263
264 } // namespace libMesh
265
266 #endif // LIBMESH_HAVE_SLEPC
267