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.
Difference between Complex and Real linear solvers
Difference between Complex and Real linear solvers
- Attachments
-
- files.zip
- (661.02 KiB) Downloaded 26 times
-
- Posts: 2237
- Joined: 25 Jan 2019, 01:28
- Antispam: Yes
Re: Difference between Complex and Real linear solvers
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)
REAL(KIND=dp) :: STIFF(nd,nd), FORCE(nd), MASS(nd,nd), DAMP(nd,nd)
! COMPLEX(KIND=dp) :: DAMP(nd,nd)
Re: Difference between Complex and Real linear solvers
Hi Kevin,
I had noticed that too. But I need to have a complex damping matrix.
Is this a fault of the complex solver?
I had noticed that too. But I need to have a complex damping matrix.
Is this a fault of the complex solver?
-
- Posts: 2237
- Joined: 25 Jan 2019, 01:28
- Antispam: Yes
Re: Difference between Complex and Real linear solvers
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)
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)
Re: Difference between Complex and Real linear solvers
Hi Kevin,
It's much better making all matrices complex.
Do you know what the problem is?
It's much better making all matrices complex.
Do you know what the problem is?
Re: Difference between Complex and Real linear solvers
The real part is much closer. But the imaginary part is way, way off.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)
-
- Site Admin
- Posts: 4812
- Joined: 22 Aug 2009, 11:57
- Antispam: Yes
- Location: Espoo, Finland
- Contact:
Re: Difference between Complex and Real linear solvers
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
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