Eigenproblem and creating a preconditioner out of a linear operator, Eigensolver Selection

[Andreas Hegendörfer]

 Hello,

I want to carry out a modal analysis of a full coupled system (piezoelectric system with displacements and potential as unknowns), which renders the eigenvalues and eigenvectors of the system. 

The following blockmatrices describe the system. The first line describes the mechanical displacements while the second line is for the electric potential: 

  _                                                                            _

 I  I¯  M 0  ¯I    *eigenvalue     +  I¯    Kmm Kme ¯I I  *eigenvector =0

 I  I_   0 0   _I                               I_     Kme Kee   _I I

  ¯                                                                            ¯

When I express the potential with the displacements I get for the eigenvalueproblem:

potential=-(Kee)^-1*Kem *displacements 

-> [M*eigenvalue+(Kmm-Kme*(Kee)^-1*Kem)]*eigenvector=0;

I implemented this with linear operators in the following way:

const auto &Kmm=stiffness.block(0,0);    const auto op_Kmm = linear_operator(Kmm);
const auto &Kme=stiffness.block(0,1);    const auto op_Kme = linear_operator(Kme);
const auto &Kem=stiffness.block(1,0);    const auto op_Kem = linear_operator(Kem);
const auto &Kee=stiffness.block(1,1);    const auto op_Kee = linear_operator(Kee);

const auto &M = mass.block(0,0); const auto op_M = linear_operator(M);

 //Kee^1
ReductionControl     reduction_control_Kee(200000, 1.0e-18, 1.0e-10);
SolverCG<>           solver_Kee(reduction_control_Kee);

PreconditionJacobi<> preconditioner_Kee;
preconditioner_Kee.initialize(Kee);
const auto op_Kee_inv = inverse_operator(op_Kee, solver_Kee, preconditioner_Kee);

 //K_Piezo=Kmm-Kme*(Kee)^-1*Kem

const auto op_K_Piezo=op_Kmm-op_Kme*op_Kee_inv*op_Kem;

//(K_Piezo)^-1 for Arpack solver

SolverControl     solver_control_K_Piezo(200000, 1.0e-10);
SolverCG<>           solver_K_Piezo(solver_control_K_Piezo);
const auto op_K_Piezo_inv=inverse_operator(op_K_Piezo,  solver_K_Piezo, PreconditionIdentity());        Here I need a preconditioner

//solve Eigenproblem

SolverControl       solver_control(dof_handler.n_dofs(), 1e-10);
ArpackSolver eigensolver(solver_control);//, additional_data);
eigensolver.solve(op_K_Piezo, op_M,  op_K_Piezo_inv, eigenvalues_Arpack, eigenfunctions_Arpack, eigenvalues_Arpack.size());

This works but the code is very slow because I dont know how to build a preconditioner out of a linear operator. The following Questions arises: 

1. Is it possible to build a preconditioner out of a linear operator in my application? 

2. Are there better ways to solve this Eigenproblem, maybe with another solver?

Thank you very much in advance, every help is welcomed!

Andreas

[Daniel Garcia-Sanchez]

Hi Andreas,

On Thursday, July 4, 2019 at 7:41:31 PM UTC+2, Andreas Hegendörfer wrote:

 Hello,

2. Are there better ways to solve this Eigenproblem, maybe with another solver?

Have you tried an spectral transformation? I'm not familiar with Arpack, I use SLEPc. I think that you could use the function ArpackSolver::set_shift(). sigma is your guess for the eigenvalue.

Note that if you use an spectral transformation, you can not use an iterative solver. You have to use a direct solver, but it considerably accelerates the eigenvalue calculation.

If your eigenvalue is from the interior of the spectrum and you don't do an spectral transformation, it can take very long.

Best,

Daniel

[Andreas Hegendörfer]

Hi Daniel,

thank you very much for your answer! 

I also tried a spectral transformation with the Arpack solver with the same result as without spectral transformation. I am interested in the smallest real eigenvalues. I know from previous calculations with the Krylov Schur solver form SLEPc that using or not using a spectral transformation makes a very big difference here. 

I am also solving the mechanical eigenproblem only: [M*eigenvalue+Kmm]*eigenvector=0

Here I am able to build a preconditioner out of the BlockSparseMatrix Kmm. The calculations with a preconditioner are really fast. 

My situation now is, that I use the same preconditioner for the mechanical-only system as for the coupled piezosystem (both build only out of Kmm). I think this should work as long as Kme*(Kee)^-1*Kem is relatively small with respect to Kmm but I will test this in further computations.

I have also tried another approach using matrices instead of linear operators. This works better but I think using linear operators is the more stylish way (this is my first contact with linear operators and I want to learn to handle them).

I do not want to use SLEPc here because I think handling the PETSc Matrices and vectors is too uncomfortable for my application. Am I right at this point? What do you think about using SLEPc here?

Thank you very much in advance, I appreciate every help!

Andreas

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a157c16e-4e0e-4b13-bd8d-ca3b2258fd30%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Daniel Garcia-Sanchez]

Hi Andreas,

On Sunday, July 7, 2019 at 4:26:39 PM UTC+2, Andreas Hegendörfer wrote:

I also tried a spectral transformation with the Arpack solver with the same result as without spectral transformation. I am interested in the smallest real eigenvalues. I know from previous calculations with the Krylov Schur solver form SLEPc that using or not using a spectral transformation makes a very big difference here. 

Yes, using a spectral transformation makes a huge difference if your eigenvalues are from the interior of the spectrum, although if you are interested in the smallest real eigenvalue, the difference might not be that important.

The drawback of an spectral transformation is that you have to use a direct solver.

You can use an iterative solver with an spectral transformation for very large problems. However, using an iterative linear solver with an spectral transformation makes the overall solution process less robust.

Although the direct solver approach may seem too costly, the factorization is only carried out at the beginning of the eigenvalue calculation and this cost is amortized in each subsequent application of the operator.

I think that the drawback of an iterative solver is that it requires lot of memory and does not scale very well. This figure can give you an idea of the scalability of MUMPS (direct parallel solver).

https://www.researchgate.net/figure/Strong-scalability-of-an-iterative-solver-with-different-preconditioners-versus-MUMPS_fig2_282172435

 

I do not want to use SLEPc here because I think handling the PETSc Matrices and vectors is too uncomfortable for my application. Am I right at this point? What do you think about using SLEPc here?

I'm writing a tutorial about how to calculate the eigenmodes of an electromagnetic cavity. So far, I've written the code, I'm writing now the documentation. You can take a look. You can find the code in this pull request:

https://github.com/dealii/dealii/pull/8345

The tutorial uses MPI, and SLEPc with std::complex.

The eigenvalues in this tutorial are from the interior of the spectrum, therefore I have to use an spectral transformation with a direct solver.

Note that PETSc has to be compiled with MUMPS if you want to run that tutorial.

Best,

Daniel

[Andreas Hegendörfer]

Hi Daniel,

thank you very much for your great help! I think your tutorial is very helpful. I think one important part of your tutorial is, that you use  SLEPcWrappers::TransformationShiftInvert shift.

When I was new to the topic of eigenvalues calculation in deal.ii a few weaks ago I only built my code out of the tutorial step 36. The Problem here was, that my code was very slow, because here is no spectral transformation applied. 

It costs me some hours to figure out what to do to solve this issue...

My code runs very well with the Arpack solver and a preconditioner built out of Kmm only. I am doing my calculations now and later on I will try using SLEPc. If problems appear I can not solve by myself I am going to ask here. 

Thank you very much for the informations and your help, Daniel!

Best,

Andreas

 

--

The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/5399942a-80d6-43b3-a2ea-7236e19e4842%40googlegroups.com.

[Anton Ermakov]

Dear Daniel,

I am working on a similar problem (planetary acoustic oscillations). I am interested in looking at your tutorial for an electromagnetic cavity, but it seems that the pull request was deleted. I wonder if it can be revived.

Thank you,

Anton.