Fortran backend crash

[Omar-A-M]

Hey,

When running nonlinear site response, the Fortran backend crashes with an array bound error. This can happen for one horizontal component but not the other, even though both horizontal components have the same number of points. Is there a restriction to the number of points?

Please see the attached error below.

Thank you.

[jsh9]

Hi @Omar-A-M , could you kindly share a sample script/notebook and the input files for me to reproduce this error? Thank you!

[Omar-A-M]

Hey @jsh9,

Please find attached a sample script and the input files to reproduce the error.

soil_profile.txt
RSN283_primary.txt
analysis.py

Thanks

[jsh9]

Hi @Omar-A-M , thanks for sending those files.

I ran your script and was able to get results:

The input & output time histories look a bit strange though. But the bottom line is, the script runs through without issues.

I ran the script on my macOS computer, and I noticed you used Windows, and I think that's the reason.

Could you confirm for me these question:

• What Windows version are you using? (7, 10, or 11)
• Is your computer 32 bit or 64 bit?
• Is your computer x86_64 (such as Intel or AMD CPU), or ARM based?

[Omar-A-M]

Hey @jsh9,

Thank you so much for the quick response.

You are correct that the input file I provided seem strange as it was the truncated version (5% to 95% Arias intensity). The fact that the analysis runs on macOS but produces a flat-line output is a key finding, as that result is not realistic.

For your reference, I have attached the full, original ground motion file before truncation ("RSN283_primary_full"). Could you please check if the analysis with this full motion produces a more reasonable result on your machine? As I get the following results for the secondary horizontal direction. using the file ("RSN283_secondary"):

RSN283_primary_full.txt

RSN283_secondary.txt

To answer your questions, here are my system details:
Windows version: Windows 10 Education
SystemType: 64-bit operating system, x64-based processor
CPU: INTEL(R) XEON(R) PLATINUM 8562Y+ 2.80 GHz

Thank you again for looking into this. It seems there may be a numerical issue in the solver that is exacerbated by the Windows environment.

[jsh9]

I ran RSN283_primary_full.txt and got this result:

It looks much more "normal" to me now.

I think your guess is correct: the amplification is so small (or, the de-amplification is so big) that there's numerical stability issues on Windows.

[Omar-A-M]

Hey @jsh9,

I have managed to reproduce the error on macOS as well with the following ground motion.

RSN216_secondary.txt

[jsh9]

Interesting...
I was able to run through RSN216_secondary.txt too. I'm using M1 Macbook Air. Nothing fancy.

[Omar-A-M]

I am using M2 pro. The analysis runs on your end but the results doesn't make sense. The input acceleration time-history is a flat line as per your figures and the output is acceleration is 20000 m/s2.

[jsh9]

It appears flat because of scale.

When I plotted it by input_motion.plot(), this is what I saw:

This is the soil profile you provided:

There's definitely something wrong. Maybe this is the same numerical stability that you encountered previously.

Let me think about how to better diagnose this issue. I'll get back to you later.

[Omar-A-M]

Thank you for the update and for taking the time to look into this.

I will wait for your next update. Please let me know if you need any more information from my end.

[Omar-A-M]

Hey @jsh9,

I have some new information that I hope is helpful for your investigation.

I have used the SeismoSoil GUI for all my ground motions and I don't face these numerical stability issues. The difference in the results are within 5%.

However, some motions in PySeismoSoil give very weird outputs even when the Fortran solver doesn't crash. Here is an example where the GUI result makes more sense:
PySeismoSoil:

SeismoSoil:

Also, regarding the ground motion that was causing the numerical stability in your latest comment. The results from SeismoSoil looks more realistic:

It seems the numerical stability issue within the solver is what is causing this. I hope these direct comparisons are useful for your debugging.

Thanks again!

[jsh9]

Hi @Omar-A-M , thanks for sharing this information!

SeismoSoil and PySeismoSoil use identical binary files (NLHH.mac). The SHA value is: 1eb33ace0f43cf15d1a4b1e2bd3343334e58de2e (you can calculate it by shasum NLHH.mac).

So I think the root cause are the HH_x parameters. PySeismoSoil uses scipy's differential_evolution algorithm to find HH_x parameters, and MATLAB uses its proprietary genetic algorithm to find HH_x parameters.

So maybe you can use MATLAB to find HH_x parameters (Section 2.2.6 of SeismoSoil's manual), save them to the hard drive, and manually load them as an HH_Param object. You can instantiate an HH_Param using a Python dict (see Section 2.1 in this Jupyter notebook).

In theory, this should make PySeismoSoil produce identical results to SeismoSoil.

[jsh9]

And you can also save your HH_x parameters from PySeismoSoil and load them to SeismoSoil. In theory, SeismoSoil results should blow up in the same way.

(I'm not longer in academia so I don't have access to MATLAB any more.)

[Omar-A-M]

Hey @jsh9,

Sorry for the late reply, had some technical issues.

I wanted to confirm that your diagnosis was exactly right. I have tried using the HH_x parameters generated by the SeismoSoil GUI and loading them into my PySeismoSoil script.

This has solved the issue. The numerical stability issues are gone, and the results from PySeismoSoil now match the results from the GUI.