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