Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "MatrixEigenvalueCheck.h"
      11             : 
      12             : #include <petscmat.h>
      13             : #include <slepceps.h>
      14             : 
      15             : registerMooseObject("NavierStokesApp", MatrixEigenvalueCheck);
      16             : 
      17             : InputParameters
      18         288 : MatrixEigenvalueCheck::validParams()
      19             : {
      20         288 :   InputParameters params = GeneralPostprocessor::validParams();
      21         288 :   params.addClassDescription("Report the number of zero eigenvalues of a matrix.");
      22         576 :   params.addRequiredParam<std::string>("mat", "The petsc binary mat file containing the matrix");
      23         576 :   params.addParam<Real>("zero_tol", 1e-8, "The tolerance for determining zero values");
      24         576 :   params.addParam<bool>("print", false, "Whether to print the eigensolves to the console");
      25         288 :   return params;
      26           0 : }
      27             : 
      28         144 : MatrixEigenvalueCheck::MatrixEigenvalueCheck(const InputParameters & parameters)
      29             :   : GeneralPostprocessor(parameters),
      30         144 :     _zero_tol(getParam<Real>("zero_tol")),
      31         288 :     _mat_name(getParam<std::string>("mat")),
      32         432 :     _print(getParam<bool>("print"))
      33             : {
      34         144 : }
      35             : 
      36             : void
      37         108 : MatrixEigenvalueCheck::initialize()
      38             : {
      39             :   PetscViewer matviewer;
      40             : 
      41         108 :   LibmeshPetscCallA(_communicator.get(), MatCreate(_communicator.get(), &_petsc_mat));
      42         108 :   LibmeshPetscCallA(
      43             :       _communicator.get(),
      44             :       PetscViewerBinaryOpen(_communicator.get(), _mat_name.c_str(), FILE_MODE_READ, &matviewer));
      45         108 :   LibmeshPetscCallA(_communicator.get(), MatLoad(_petsc_mat, matviewer));
      46         108 :   LibmeshPetscCallA(_communicator.get(), PetscViewerDestroy(&matviewer));
      47             : 
      48         108 :   _num_zero_eigenvalues = 0;
      49         108 : }
      50             : 
      51             : void
      52         108 : MatrixEigenvalueCheck::execute()
      53             : {
      54         108 :   if (_print)
      55           0 :     _console << std::endl << "Conducting eigen-solve for " << _mat_name << std::endl;
      56             : 
      57             :   PetscScalar kr, ki;
      58             :   PetscReal error, re, im;
      59             :   PetscInt nconv, nev, i;
      60             :   EPS eps;
      61             : 
      62         108 :   LibmeshPetscCallA(_communicator.get(), EPSCreate(_communicator.get(), &eps));
      63         108 :   LibmeshPetscCallA(_communicator.get(), EPSSetOperators(eps, _petsc_mat, nullptr));
      64         108 :   LibmeshPetscCallA(_communicator.get(), EPSSetType(eps, EPSLAPACK));
      65         108 :   LibmeshPetscCallA(_communicator.get(), EPSSolve(eps));
      66         108 :   LibmeshPetscCallA(_communicator.get(), EPSGetDimensions(eps, &nev, nullptr, nullptr));
      67         108 :   if (_print)
      68             :   {
      69           0 :     LibmeshPetscCallA(_communicator.get(),
      70             :                       PetscPrintf(_communicator.get(),
      71             :                                   " Number of requested eigenvalues: %" PetscInt_FMT "\n",
      72             :                                   nev));
      73             :   }
      74             : 
      75         108 :   LibmeshPetscCallA(_communicator.get(), EPSGetConverged(eps, &nconv));
      76         108 :   if (_print)
      77             :   {
      78           0 :     LibmeshPetscCallA(_communicator.get(),
      79             :                       PetscPrintf(_communicator.get(),
      80             :                                   " Number of converged eigenpairs: %" PetscInt_FMT "\n\n",
      81             :                                   nconv));
      82             :   }
      83             : 
      84         108 :   if (nconv > 0)
      85             :   {
      86         108 :     if (_print)
      87             :     {
      88             :       /*
      89             :          Display eigenvalues and relative errors
      90             :       */
      91           0 :       LibmeshPetscCallA(_communicator.get(),
      92             :                         PetscPrintf(_communicator.get(),
      93             :                                     "           k          ||Ax-kx||/||kx||\n"
      94             :                                     "   ----------------- ------------------\n"));
      95             :     }
      96             : 
      97       35658 :     for (i = 0; i < nconv; i++)
      98             :     {
      99             :       /*
     100             :         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
     101             :         ki (imaginary part)
     102             :       */
     103       35550 :       LibmeshPetscCallA(_communicator.get(), EPSGetEigenpair(eps, i, &kr, &ki, nullptr, nullptr));
     104             :       /*
     105             :          Compute the relative error associated to each eigenpair
     106             :       */
     107       35550 :       LibmeshPetscCallA(_communicator.get(), EPSComputeError(eps, i, EPS_ERROR_RELATIVE, &error));
     108             : 
     109             : #if defined(PETSC_USE_COMPLEX)
     110             :       re = PetscRealPart(kr);
     111             :       im = PetscImaginaryPart(kr);
     112             : #else
     113       35550 :       re = kr;
     114       35550 :       im = ki;
     115             : #endif
     116       35550 :       const bool zero_im = std::abs(im) < _zero_tol;
     117       35550 :       const bool zero_re = std::abs(re) < _zero_tol;
     118       35550 :       if (_print)
     119             :       {
     120           0 :         if (!zero_im)
     121           0 :           LibmeshPetscCallA(
     122             :               _communicator.get(),
     123             :               PetscPrintf(
     124             :                   PETSC_COMM_WORLD, " %9f+%9fi %12g\n", (double)re, (double)im, (double)error));
     125             :         else
     126           0 :           LibmeshPetscCallA(
     127             :               _communicator.get(),
     128             :               PetscPrintf(_communicator.get(), "   %12f       %12g\n", (double)re, (double)error));
     129             :       }
     130       35550 :       if (zero_im && zero_re)
     131       13086 :         ++_num_zero_eigenvalues;
     132             :     }
     133             : 
     134         108 :     if (_print)
     135             :     {
     136           0 :       LibmeshPetscCallA(_communicator.get(), PetscPrintf(_communicator.get(), "\n"));
     137             :     }
     138             :   }
     139         108 : }
     140             : 
     141             : void
     142         108 : MatrixEigenvalueCheck::finalize()
     143             : {
     144         108 : }
     145             : 
     146             : Real
     147         108 : MatrixEigenvalueCheck::getValue() const
     148             : {
     149         108 :   return _num_zero_eigenvalues;
     150             : }