Difference between Complex and Real linear solvers

Numerical methods and mathematical models of Elmer
Post Reply
supreet
Posts: 19
Joined: 22 Nov 2022, 16:44
Antispam: Yes

Difference between Complex and Real linear solvers

Post by supreet »

I need to create complex valued matrices for an eigenvalue problem I am solving. I have modified WaveSolver.f90 to suit my requirements, mainly by declaring damping matrix and sound reaction term as complex. I've attached the modified solver module, the sif and mesh files in the zip file below.

To verify the validity of the solver, I am solving a damped eigenvalue problem with reaction term equal to 1.0. The in-built solver correctly calculates the eigenvalues, which I know from analytical solution. The reaction term is hard-coded as a complex value (1,0) in the custorm solver. The eigenvalues from custom solver are a bit off but close.

I'd be very thankful to anyone who points out the issue.
Attachments
files.zip
(661.02 KiB) Downloaded 26 times
kevinarden
Posts: 2237
Joined: 25 Jan 2019, 01:28
Antispam: Yes

Re: Difference between Complex and Real linear solvers

Post by kevinarden »

When I make this one change, two places, then your solver gets the same answer as the built in solver

REAL(KIND=dp) :: STIFF(nd,nd), FORCE(nd), MASS(nd,nd), DAMP(nd,nd)
! COMPLEX(KIND=dp) :: DAMP(nd,nd)
supreet
Posts: 19
Joined: 22 Nov 2022, 16:44
Antispam: Yes

Re: Difference between Complex and Real linear solvers

Post by supreet »

Hi Kevin,

I had noticed that too. But I need to have a complex damping matrix.

Is this a fault of the complex solver?
kevinarden
Posts: 2237
Joined: 25 Jan 2019, 01:28
Antispam: Yes

Re: Difference between Complex and Real linear solvers

Post by kevinarden »

If I make these changes in two places, the answer becomes much closer

Complex(KIND=dp) :: MASS(nd,nd), STIFF(nd,nd), FORCE(nd)

MASS(:,:) = (0._dp, 0._dp)
STIFF(:,:) = (0._dp, 0._dp)
FORCE(:) = (0._dp, 0._dp)
DAMP(:,:) = (0._dp, 0._dp)
supreet
Posts: 19
Joined: 22 Nov 2022, 16:44
Antispam: Yes

Re: Difference between Complex and Real linear solvers

Post by supreet »

Hi Kevin,

It's much better making all matrices complex.

Do you know what the problem is?
supreet
Posts: 19
Joined: 22 Nov 2022, 16:44
Antispam: Yes

Re: Difference between Complex and Real linear solvers

Post by supreet »

kevinarden wrote: 21 Dec 2022, 17:40 If I make these changes in two places, the answer becomes much closer

Complex(KIND=dp) :: MASS(nd,nd), STIFF(nd,nd), FORCE(nd)

MASS(:,:) = (0._dp, 0._dp)
STIFF(:,:) = (0._dp, 0._dp)
FORCE(:) = (0._dp, 0._dp)
DAMP(:,:) = (0._dp, 0._dp)
The real part is much closer. But the imaginary part is way, way off.
raback
Site Admin
Posts: 4812
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

Re: Difference between Complex and Real linear solvers

Post by raback »

Hi,

Most transient solvers of Elmer can be run in harmonic and eigenmode mode too. Assume any equation of the type

Mx'' + Dx' + Kx' = f

where M would be mass matrix, D damping matrix and K the stiffness matrix following naming convention from structucal mechanics. If you solve the system as transient, you need to replace time derivatives x' and x'' with suitable time discretization. If you solve harmonic or eigenmode you make an ansatz x=x(r)*exp(iwt), where i=sqrt(-1), resulting to time derivative being replaced as: d/dt -> iw.

In harmonic w is known and f is assumed to be in that frequency. In eigenmodes f is zero and we need to find the eigenvalues as part of the problem. As you can see in the form the resulting system will be complex if D is not zero.

Unfortunately the path you chose is not really required. The above described treatment of the matrices M, D and K is done fully on the library level. As an example you can see the WaveSolver having two test cases: WaveEqu and WaveEquEigen. The latter is an eigen system of the problem.

Now, depending on whether D is considered or not x will be real or complex. The truly complex eigensystem will be more costly to solve. Hence there is user control to decide how the system is treated:

Eigen System Complex = True / False

-Peter
Post Reply