forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

Porting SLATEC Library -- Undocumented routines in SALFLIBC
Goto page Previous  1, 2, 3, 4  Next
 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support
View previous topic :: View next topic  
Author Message
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Mon Mar 13, 2023 3:42 pm    Post subject: Reply with quote

Thank you, Paul.
Back to top
View user's profile Send private message
Kenneth_Smith



Joined: 18 May 2012
Posts: 697
Location: Hamilton, Lanarkshire, Scotland.

PostPosted: Tue Mar 14, 2023 7:42 pm    Post subject: Reply with quote

mecej4,

Below is a link to another example (a well known problem in the electrical world) which uses the some of the SLATEC eigenvalue and eigenvector routines.

https://www.dropbox.com/s/l8mgxhbeiyjdix5/ieee1.f95?dl=0

I resisted the temptation to plot the mode shapes with %pl. The results all look to be correct to me.

Here is a better, tidied up version:

https://www.dropbox.com/s/58djudpp1fvkguf/ieee2.f95?dl=0
Back to top
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Wed Mar 15, 2023 4:50 am    Post subject: Reply with quote

Kenneth,

I think that a reference/link to the source of the "magic" constants in arrays HC and TS should be added. How about https://link.springer.com/content/pdf/bbm:978-1-4615-5633-6/1.pdf ? Table A3 there lists the constants.

The comments in the code refer to routines SPPDI and SDEEV, but there are no calls to those in that code.

The code would be a bit easier to read and understand if array operations were used. For example:
Code:
!     Do this to avoid very small negative numbers which cause problems with SQRT later
      where (abs(hm) < eps) hm = 0.0
      where (abs(km) < eps) km = 0.0
Back to top
View user's profile Send private message
Kenneth_Smith



Joined: 18 May 2012
Posts: 697
Location: Hamilton, Lanarkshire, Scotland.

PostPosted: Thu Mar 16, 2023 11:02 am    Post subject: Reply with quote

mecej4,

A revised version can be found here:

https://www.dropbox.com/s/9s6jenwi8qvc2xm/ieee3.f95?dl=0

At the end of the revised code, a clearwin+ window is opened which outputs all the results in graphical format, which is appropriate for a dll intended for use with FTN95

I have added a description of the general method for solving this problem, and more comments in the code, which is now separated into a number of distinct steps.

In this particular problem, the parameters for the rotating mass inertia and shaft stiffness (the “magic” numbers) are not in standard units, but expressed in the per unit system used by electrical engineers for dynamic studies. The definition of these units is explained, as is the necessary departures from the “general” method. As a consequence of this, the intermediate values of the eigenvalues reported by the program do not correspond to omega**2. The system base frequency is correctly factored in at the end when reporting the actual torsional frequencies in Hz, i.e. this takes us back to the results that would be obtained in we started with rotating mass inertia J and shaft torsional stiffness K in their SI units.

The parameters in the benchmark model correspond to the large turbine generators at Mohave, Modes 2, 4, and 5 can be easily excited if their complementary electrical frequencies are present in the connected power system, due to the high participation of the generator mass in these modes as show by the mode shapes (eigenvalues). Two shaft failures occurred at Mohave in the early 1970s, which led to the development of the benchmark model for analysis purposes. If the electrical system produces frequencies complementary to one or more the torsional modes you can very rapidly fatigue the mechanical string.

Upon reflection, this was not perhaps a good example, as somebody blindly reusing the code for a torsional problem with J in kg,^2 and K in Nm/rad will run into problems – system electrical frequency is not relevant with this approach. Nevertheless it does demonstrate the basic methodology, the use of some of the SLATEC routines, and was an interesting problem for the “wee small hours” of the last few days. A refresher – I’ve not looked at this problem in earnest for at least a decade. I’m a bit fatigued after this one.

PS I can see a way of fixing all of these issues. But it will be some time before I can implement this, as there a a few other urgent things on the to-do list at the moment.
Back to top
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Thu Mar 16, 2023 10:37 pm    Post subject: Reply with quote

Kenneth, thanks for writing up a detailed explanation. I found a set of slides for a presentation on the Mohave turbine failure:
https://www.geenergyconsulting.com/content/dam/Energy_Consulting/global/en_US/pdfs/GE_-_SSO_Risk_Analysis_Protection_and_Mitigation_Techniques.pdf

One of the figures in the presentation shows mode shapes, and it struck me that it was very similar to the graph that your program produces. Neat.

A more detailed discussion of SSR (Sub-Synchronous Resonance):

https://energiforskmedia.blob.core.windows.net/media/25107/torsional-vibrations-in-steam-turbine-shaft-trains-energiforskrapport-2018-522.pdf

Since the installation where the failure occurred is in the US, recasting the description into SI units could be considered sacrilegious.
Back to top
View user's profile Send private message
Kenneth_Smith



Joined: 18 May 2012
Posts: 697
Location: Hamilton, Lanarkshire, Scotland.

PostPosted: Fri Mar 17, 2023 3:06 pm    Post subject: Reply with quote

Mecej4,

Thanks for sharing both of these links.

The Energiforsk document appears to be mostly drawn from an earlier EPRI report which is on my desk, updated to reflect the more recent advances in SSR protection methods.

The GE document which you shared confirms something I have long expected. The use of series filters to block troublesome frequencies takes up a lot of physical space. The photograph of the filters at Navajo (plant now decommissioned) confirms this. A very good example of a purely analogue (R/L/C) solution to a complex problem.

Some pictures of the more conventional shunt filters on the IFA1 HVDC link on the UK side at Sellindge which also have a large footprint. The rack mounted capacitors and the cylindrical air cored reactors can be clearly seen.

https://www.geograph.org.uk/photo/6964452

My first attempts to understand the torsional problem were considered sacrilegious by the mechanical engineering team; I converted the mechanical system to its electrical equivalent circuit (series L to represent the shafts and C to ground to represent the inertias). The torsional frequencies appeared as local maxima/minima in the impedance vs. frequency scans at each of the nodes in the equivalent circuit, and from the off diagonal elements of the full impedance matrix at these specific frequencies the mode shapes were derived. Even although we agreed that the system differential equations were essentially the same, they remained sceptical for a long time!
Back to top
View user's profile Send private message Visit poster's website
DanRRight



Joined: 10 Mar 2008
Posts: 2813
Location: South Pole, Antarctica

PostPosted: Fri Mar 17, 2023 6:10 pm    Post subject: Reply with quote

There are a lot of fast Fourier transform subroutines in this library. Can anyone share experience ***and demo examples*** with these routines or any other libraries?

Need the ones which will deal with
1) higher order harmonics 10-100w
2) lower order around ~2w
3) and also those like 1/2, 3/2w etc
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Fri Mar 17, 2023 8:42 pm    Post subject: Reply with quote

Dan, it is nice to see you expressing interest in the SLATEC DLLs for FTN95.

Of the 54 tests that come with the package, Test51 is devoted to the FFT routines. The comments in the code indicate that the FFT routines originated in NCAR (U.S. Natl. Center for Atmospheric Research), and are dated 1985 or earlier. From the comment lines in the pertinent source file, FFTQX, we are able to note:
Code:
!             This Program Tests The Package Of Fast Fourier
!     Transforms For Both Complex And Real Periodic Sequences And
!     Certain Other Symmetric Sequences That Are Listed Below.
!
!     1.   Rffti     Initialize  Rfftf And Rfftb
!     2.   Rfftf     Forward Transform Of A Real Periodic Sequence
!     3.   Rfftb     Backward Transform Of A Real Coefficient Array
!
!     4.   Ezffti    Initialize Ezfftf And Ezfftb
!     5.   Ezfftf    A Simplified Real Periodic Forward Transform
!     6.   Ezfftb    A Simplified Real Periodic Backward Transform
!
!     7.   Sinti     Initialize Sint
!     8.   Sint      Sine Transform Of A Real Odd Sequence
!
!     9.   Costi     Initialize Cost
!     10.  Cost      Cosine Transform Of A Real Even Sequence
!
!     11.  Sinqi     Initialize Sinqf And Sinqb
!     12.  Sinqf     Forward Sine Transform With Odd Wave Numbers
!     13.  Sinqb     Unnormalized Inverse Of Sinqf
!
!     14.  Cosqi     Initialize Cosqf And Cosqb
!     15.  Cosqf     Forward Cosine Transform With Odd Wave Numbers
!     16.  Cosqb     Unnormalized Inverse Of Cosqf
!
!     17.  Cffti     Initialize Cfftf And Cfftb
!     18.  Cfftf     Forward Transform Of A Complex Periodic Sequence
!     19.  Cfftb     Unnormalized Inverse Of Cfftf


I do not think that SLATEC should be regarded as a replacement for specialized or more recent packages for solving problems in any of the numerous areas that it covers. On the other hand, it is free, high quality but old, and the availability of DLLs may make it easier for someone such as yourself to try it out. The download and documentation links that I provided above should get you started.

As for demos, other than the 54 quick checks included in the package, we have a handful, including those that Kenneth has contributed. If you can prepare one for FFTs, go for it! As Kenneth has demonstrated, less than 10 lines of Fortran code will do to produce very good graphs using FTN95. Such demos may persuade prospective users to continue with Fortran instead of switching to Python, Julia, etc., for the graphics.
Back to top
View user's profile Send private message
Kenneth_Smith



Joined: 18 May 2012
Posts: 697
Location: Hamilton, Lanarkshire, Scotland.

PostPosted: Sun Mar 19, 2023 11:52 am    Post subject: Reply with quote

mecej4,

Final version of the eigenvalue example can be found here.

https://www.dropbox.com/s/lggp8hvaufa7hm5/ieee4.f95?dl=0

Incorrect labels (mode numbers) in the graphics output corrected and the description now contains a physical interpretation of the mode shapes.
Back to top
View user's profile Send private message Visit poster's website
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Mon Mar 20, 2023 5:17 am    Post subject: Reply with quote

Thanks, Kenneth.

I am waiting for Paul's fix for the complex arithmetic bug to be released. When that occurs, I can rebuild the SLATEC DLLs and make them available, along with the modified sources.

In the meantime, I am replacing a large number of code segments in the SLATEC sources similar to
Code:
         Ws(iw+1) = zero
         call scopy(n,Ws(iw+1),0,Ws(iw+1),1)

by the modern Fortran equivalent

Code:
Ws(iw+1:iw+n) = Zero
Back to top
View user's profile Send private message
Kenneth_Smith



Joined: 18 May 2012
Posts: 697
Location: Hamilton, Lanarkshire, Scotland.

PostPosted: Mon Mar 20, 2023 11:50 am    Post subject: Reply with quote

mecej4,

I produced this using RFFTF1.

https://www.dropbox.com/s/zlc7dyrtbd86ldt/fft.png?dl=0

Basically extracting the harmonic content of a signal which is hidden within random noise.

Not yet completed, as the applied scaling factors associated with the FFT output are by inspection rather than calculated values. The description of the output from RFFTF1 is difficult to understand upon first reading.

Also, I discovered that the Silverfrost library routine random@ is not PURE and cannot be used in a FORALL construct.
Back to top
View user's profile Send private message Visit poster's website
DanRRight



Joined: 10 Mar 2008
Posts: 2813
Location: South Pole, Antarctica

PostPosted: Tue Mar 21, 2023 4:42 pm    Post subject: Reply with quote

Yes, examples of use of library is the key for its popularity. Such examples of course are hard to produce but they quickly show usability of the method. Without nice examples, specifically in some presentable visual forms like plots and tables any library will be just slowly forgotten. FFT is specifically good at showing spectra, harmonics etc
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Tue Mar 21, 2023 6:14 pm    Post subject: Reply with quote

New versions of the SLATEC DLLs may be downloaded from

https://drive.google.com/file/d/10spK3G9RHFcNevrJgDtbBONqyZshJSKs/view?usp=sharing

Please read README.txt before using the DLLs.

With the new DLLs, TEST04 of SLATEC, which failed with FTN95 /64 earlier, now works correctly.

Thanks to Kenneth Smith for contributing more examples of programs that use SLATEC.
Back to top
View user's profile Send private message
Kenneth_Smith



Joined: 18 May 2012
Posts: 697
Location: Hamilton, Lanarkshire, Scotland.

PostPosted: Wed Mar 22, 2023 9:31 am    Post subject: Reply with quote

Well done.

Proof that the statement by Ellis in his text “Fortran 77 Programming” published 33 years ago, “Programs written in Fortran 77 will have life well into the 21st century”, was not in error.
Back to top
View user's profile Send private message Visit poster's website
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Wed Mar 22, 2023 12:42 pm    Post subject: Reply with quote

That would be the 2nd edition. I wonder what he said in the 1st?
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Support All times are GMT + 1 Hour
Goto page Previous  1, 2, 3, 4  Next
Page 3 of 4

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group