Segmentation fault running dynamic analysis

Forum for asking and answering questions related to use of the OpenSeesPy module

Moderators: silvia, selimgunay, Moderators

Post Reply
AnnaKowal
Posts: 14
Joined: Thu Oct 25, 2018 12:54 pm
Location: University of Otago

Segmentation fault running dynamic analysis

Post by AnnaKowal »

Hi,

I got a segmentation fault (core dumped) on anlyze step when running site response analysis of a soil deposit with a parabolically varying shear wave velocity profile on an elastic half-space in OpenSeesPy. The finite rigidity of the underlying medium is considered through the use of a viscous damper at the base of the soil column.

The analysis went smoothly through gravity analysis and return Segmentation fault during the dynamic analysis on command:
analyze (motionSteps, motionDT),
where motionSteps and motionDT are defined previously, 14000 and 0.005 respectively.

Can anyone help me to address that issue,
Many thanks,
Anna
Top
mhscott
Posts: 880
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:
Contact mhscott

Re: Segmentation fault running dynamic analysis

Post by mhscott »

You'll have to provide your code. Could be so many things. If you post it here, you should use the code block (the </> button above) so that Python indentation is preserved.
Top
AnnaKowal
Posts: 14
Joined: Thu Oct 25, 2018 12:54 pm
Location: University of Otago

Re: Segmentation fault running dynamic analysis

Post by AnnaKowal »

Hi,

Thanks for the response. The script is quite massive - see below. Too long to put here as it is originally - remove most of the nodes and elements definition. It fails on analyze (motionSteps, motionDT) - line 1109

Many thanks,
Anna

Code: Select all

###########################################################
#                                                         #
# Site response analysis of a soil deposit with a         # 
# parabolically varying shear wave velocity profile on an #
# elastic half-space.  The finite rigidity of the         #
# underlying medium is considered through the use of a    #
# viscous damper at the base of the soil column.          #
#                                                         #
#   Created by:  Chris McGann                             #
#                HyungSuk Shin                            #
#                Pedro Arduino                            #
#                Peter Mackenzie-Helnwein                 #
#              --University of Washington--               #
#                                                         #
# ---> Basic units are kN and m   (unless specified)      #
#                                                         #
###########################################################

from openseespy.opensees import *
from math import pi, floor
import datetime

def iter(count):
    return range(1, count + 1)

def mkLayerArr(init=0):
    global numLayers
    return [init] * (numLayers + 1)

wipe()

#-----------------------------------------------------------------------------------------------------------
#  1. CREATE SOIL NODES AND FIXITIES
#-----------------------------------------------------------------------------------------------------------
# soil nodes are created in 2 dimensions, with 3 dof (2 translational, 1 porePressure)
model('basic', '-ndm', 2, '-ndf', 2)

#define soil nodes
node(1,0.00000,64.00000)
#.......
node(17302,440.59610,6.78066)
node(17303,386.10930,11.34698)
print (f"Finished defyning all nodes...")

# define fixities for soil nodes
# define fixities for model base
fix(2,1,1)
fix(18,1,1)
fix(98,1,1)
fix(99,1,1)
fix(100,1,1)
fix(101,1,1)
fix(102,1,1)
fix(103,1,1)
fix(104,1,1)
fix(105,1,1)
fix(106,1,1)
fix(107,1,1)
fix(108,1,1)
fix(109,1,1)
fix(110,1,1)
fix(111,1,1)
fix(112,1,1)
fix(113,1,1)
fix(114,1,1)
fix(115,1,1)
fix(116,1,1)
fix(117,1,1)
fix(118,1,1)
fix(119,1,1)
fix(120,1,1)
fix(121,1,1)
fix(122,1,1)
fix(123,1,1)
fix(124,1,1)
fix(125,1,1)
fix(126,1,1)
fix(127,1,1)
fix(128,1,1)
fix(129,1,1)
fix(130,1,1)
fix(131,1,1)
fix(132,1,1)
fix(133,1,1)
fix(134,1,1)
fix(135,1,1)
fix(136,1,1)
fix(137,1,1)
fix(138,1,1)
fix(139,1,1)
fix(140,1,1)
fix(141,1,1)
fix(142,1,1)
fix(143,1,1)
fix(144,1,1)
fix(145,1,1)
fix(146,1,1)
fix(147,1,1)
fix(148,1,1)
fix(149,1,1)
fix(150,1,1)
fix(151,1,1)
fix(152,1,1)
fix(153,1,1)
fix(154,1,1)
fix(155,1,1)
fix(156,1,1)
fix(157,1,1)
fix(158,1,1)
fix(159,1,1)
fix(160,1,1)
fix(161,1,1)
fix(162,1,1)
fix(163,1,1)
fix(164,1,1)
fix(165,1,1)
fix(166,1,1)
fix(167,1,1)
fix(168,1,1)
fix(169,1,1)
fix(170,1,1)
fix(171,1,1)
fix(172,1,1)
fix(173,1,1)
fix(174,1,1)
fix(175,1,1)
fix(176,1,1)
fix(177,1,1)
fix(178,1,1)
fix(179,1,1)
fix(180,1,1)
fix(181,1,1)
fix(182,1,1)
fix(183,1,1)
fix(184,1,1)
fix(185,1,1)
fix(186,1,1)
fix(187,1,1)
fix(188,1,1)
fix(189,1,1)
fix(190,1,1)
fix(191,1,1)
fix(192,1,1)
fix(193,1,1)
fix(194,1,1)
fix(195,1,1)
fix(196,1,1)
fix(197,1,1)
fix(198,1,1)
fix(199,1,1)
fix(200,1,1)

print (f"Finished defining fixities for model base...")
print (f"Finished creating all -ndf 2 boundary conditions...")

# define equal degrees of freedom for free-field columns
equalDOF(19,1,1,2)
equalDOF(18,2,1,2)
equalDOF(203,94,1,2)
equalDOF(202,90,1,2)
equalDOF(201,91,1,2)
equalDOF(205,95,1,2)
equalDOF(200,92,1,2)
equalDOF(206,96,1,2)
equalDOF(204,93,1,2)
print (f"Finished creating equalDOF constraints for soil columns...")

#-----------------------------------------------------------------------------------------------------------
#  2. DEFINE ANALYSIS PARAMETERS
#-----------------------------------------------------------------------------------------------------------

#---SOIL GEOMETRY
# number of soil layers
numLayers       = 6
# layer thicknesses (m)
layerThick = mkLayerArr()
#layerThick[6] = 5.0
#layerThick[5] = 10.0
#layerThick[4] = 10.0
#layerThick[3] = 25.0
#layerThick[2] = 50.0
#layerThick[1] = 100.0

layerThick[6] = 1.0
layerThick[5] = 1.0
layerThick[4] = 1.0
layerThick[3] = 1.0
layerThick[2] = 1.0
layerThick[1] = 1.0

#---MATERIAL PROPERTIES
# soil mass density (Mg/m^3)
rho = mkLayerArr()
rho[6]   =       1.6
rho[5]   =       2.0
rho[4]   =       2.2
rho[3]   =       2.3
rho[2]   =       2.575
rho[1]   =       2.6

# soil shear wave velocity for each layer(m/s)
Vs = mkLayerArr()
Vs[6]           = 135.0
Vs[5]           = 180.0
Vs[4]           = 280.0
Vs[3]           = 900.0
Vs[2]           = 2300.0
Vs[1]           = 2800.0

# soil shear modulus for each layer (kPa)
G = mkLayerArr()
G[6]           = rho[6] * Vs[6] * Vs[6]
G[5]           = rho[5] * Vs[5] * Vs[5]
G[4]           = rho[4] * Vs[4] * Vs[4]
G[3]           = rho[3] * Vs[3] * Vs[3]
G[2]           = rho[2] * Vs[2] * Vs[2]
G[1]           = rho[1] * Vs[1] * Vs[1]
# poisson's ratio of soil
nu              = 0.25
# soil elastic modulus for each layer (kPa)
E = [2 * Gi * (1 + nu) for Gi in G]
# soil bulk modulus for each layer (kPa)
bulk = [Ei / (3 * (1 - 2 * nu)) for Ei in E]
# soil friction angle
phi           = 35.0
# peak shear strain
gammaPeak       = 0.1
# reference pressure
refPress        = 80.0
# pressure dependency coefficient
pressCoeff = 0.5 
# phase transformation angle
phaseAng        = 27.0
# contraction
contract        = 0.03
# dilation coefficients
dilate1         = 0.8
dilate2         =   5
# liquefaction coefficients
liq1            = 0.0 
liq2            = 0.0 
liq3            = 0.0

print (f"Finished defining parameters for all soil materials...")

#-----------------------------------------------------------------------------------------------------------
#  3. DEFINE SOIL MATERIALS
#-----------------------------------------------------------------------------------------------------------                                      
nDMaterial ('PressureDependMultiYield02', 6, 2.0, rho[6], G[6], bulk[6], phi, gammaPeak,
                                      refPress, pressCoeff, phaseAng, 0.067, 0.23, 0.06,
                                      0.27, 20.0, 5.0, 3.0, 1.0, 
                                      0.0, 0.73, 0.9, 0.02, 0.7, 101.0, 15.0)                                 

nDMaterial('PressureDependMultiYield02', 5, 2.0, rho[5], G[5], bulk[5], phi, gammaPeak,
                                     refPress, pressCoeff, phaseAng, 0.067, 0.23, 0.06,
                                    0.27, 20.0, 5.0, 3.0, 1.0, 
                                      0.0, 0.73, 0.9, 0.02, 0.7, 101.0, 15.0)
                                      
nDMaterial('PressureDependMultiYield02', 4, 2.0, rho[4], G[4], bulk[4], phi, gammaPeak,
                                     refPress, pressCoeff, phaseAng, 0.067, 0.23, 0.06,
                                    0.27, 20.0, 5.0, 3.0, 1.0, 
                                      0.0, 0.73, 0.9, 0.02, 0.7, 101.0, 15.0)
                                      
nDMaterial('PressureDependMultiYield02', 3, 2.0, rho[3], G[3], bulk[3], phi, gammaPeak,
                                     refPress, pressCoeff, phaseAng, 0.067, 0.23, 0.06,
                                    0.27, 20.0, 5.0, 3.0, 1.0, 
                                      0.0, 0.73, 0.9, 0.02, 0.7, 101.0, 15.0)  
                                      
nDMaterial('PressureDependMultiYield02', 2, 2.0, rho[2], G[2], bulk[2], phi, gammaPeak,
                                     refPress, pressCoeff, phaseAng, 0.067, 0.23, 0.06,
                                    0.27, 20.0, 5.0, 3.0, 1.0, 
                                      0.0, 0.73, 0.9, 0.02, 0.7, 101.0, 15.0)                                                                          
                                      	
nDMaterial('ElasticIsotropic',1, E[1], 0.25, rho[1])	

print("Finished creating all soil materials...")

#-----------------------------------------------------------------------------------------------------------
#  4. DEFINE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------------------------
wgtX = mkLayerArr()
wgtX1            = 0.0 
wgtX2            = 0.0 
wgtX3            = 0.0 
wgtX4            = 0.0 
wgtX5            = 0.0 
wgtX6            = 0.0 
#for k in iter(numLayers):
#		wgtY  = -9.81*rho[k]

wgtY1            = -9.81*rho[1]
wgtY2            = -9.81*rho[2]
wgtY3            = -9.81*rho[3]
wgtY4            = -9.81*rho[4]
wgtY5            = -9.81*rho[5] 
wgtY6            = -9.81*rho[6]
#count = 0
element('quad',1,115,149,3908,3907,layerThick[1],'PlaneStrain',1, 0.0, 0.0, wgtX1, wgtY1)
#.....................
element('quad',16834,3648,16949,17303,3679,layerThick[6],'PlaneStrain',6, 0.0, 0.0, wgtX6, wgtY6)
print(f"Finished creating all soil elements..")


#-----------------------------------------------------------------------------------------------
#  5. LYSMER DASHPOT
#-----------------------------------------------------------------------------------------------

# Bottom boundary
# define dashpot master nodes
node(20000,4105.00000,-200.00000)
node(20001,4065.52900,-200.00000)
node(20002,4026.05800,-200.00000)
node(20003,3986.58700,-200.00000)
node(20004,3947.11500,-200.00000)
node(20005,3907.64400,-200.00000)
node(20006,3868.17300,-200.00000)
node(20007,3828.70200,-200.00000)
node(20008,3789.23100,-200.00000)
node(20009,3749.76000,-200.00000)
node(20010,3710.28800,-200.00000)
node(20011,3670.81700,-200.00000)
node(20012,3631.34600,-200.00000)
node(20013,3591.87500,-200.00000)
node(20014,3552.40400,-200.00000)
node(20015,3512.93300,-200.00000)
node(20016,3473.46200,-200.00000)
node(20017,3433.99000,-200.00000)
node(20018,3394.51900,-200.00000)
node(20019,3355.04800,-200.00000)
node(20020,3315.57700,-200.00000)
node(20021,3276.10600,-200.00000)
node(20022,3236.63500,-200.00000)
node(20023,3197.16300,-200.00000)
node(20024,3157.69200,-200.00000)
node(20025,3118.22100,-200.00000)
node(20026,3078.75000,-200.00000)
node(20027,3039.27900,-200.00000)
node(20028,2999.80800,-200.00000)
node(20029,2960.33700,-200.00000)
node(20030,2920.86500,-200.00000)
node(20031,2881.39400,-200.00000)
node(20032,2841.92300,-200.00000)
node(20033,2802.45200,-200.00000)
node(20034,2762.98100,-200.00000)
node(20035,2723.51000,-200.00000)
node(20036,2684.03800,-200.00000)
node(20037,2644.56700,-200.00000)
node(20038,2605.09600,-200.00000)
node(20039,2565.62500,-200.00000)
node(20040,2526.15400,-200.00000)
node(20041,2486.68300,-200.00000)
node(20042,2447.21200,-200.00000)
node(20043,2407.74000,-200.00000)
node(20044,2368.26900,-200.00000)
node(20045,2328.79800,-200.00000)
node(20046,2289.32700,-200.00000)
node(20047,2249.85600,-200.00000)
node(20048,2210.38500,-200.00000)
node(20049,2170.91300,-200.00000)
node(20050,2131.44200,-200.00000)
node(20051,2091.97100,-200.00000)
node(20052,2052.50000,-200.00000)
node(20053,2013.02900,-200.00000)
node(20054,1973.55800,-200.00000)
node(20055,1934.08700,-200.00000)
node(20056,1894.61500,-200.00000)
node(20057,1855.14400,-200.00000)
node(20058,1815.67300,-200.00000)
node(20059,1776.20200,-200.00000)
node(20060,1736.73100,-200.00000)
node(20061,1697.26000,-200.00000)
node(20062,1657.78800,-200.00000)
node(20063,1618.31700,-200.00000)
node(20064,1578.84600,-200.00000)
node(20065,1539.37500,-200.00000)
node(20066,1499.90400,-200.00000)
node(20067,1460.43300,-200.00000)
node(20068,1420.96200,-200.00000)
node(20069,1381.49000,-200.00000)
node(20070,1342.01900,-200.00000)
node(20071,1302.54800,-200.00000)
node(20072,1263.07700,-200.00000)
node(20073,1223.60600,-200.00000)
node(20074,1184.13500,-200.00000)
node(20075,1144.66300,-200.00000)
node(20076,1105.19200,-200.00000)
node(20077,1065.72100,-200.00000)
node(20078,1026.25000,-200.00000)
node(20079,986.77880,-200.00000)
node(20080,947.30770,-200.00000)
node(20081,907.83650,-200.00000)
node(20082,868.36540,-200.00000)
node(20083,828.89420,-200.00000)
node(20084,789.42310,-200.00000)
node(20085,749.95190,-200.00000)
node(20086,710.48080,-200.00000)
node(20087,671.00960,-200.00000)
node(20088,631.53850,-200.00000)
node(20089,592.06730,-200.00000)
node(20090,552.59620,-200.00000)
node(20091,513.12500,-200.00000)
node(20092,473.65380,-200.00000)
node(20093,434.18270,-200.00000)
node(20094,394.71150,-200.00000)
node(20095,355.24040,-200.00000)
node(20096,315.76920,-200.00000)
node(20097,276.29810,-200.00000)
node(20098,236.82690,-200.00000)
node(20099,197.35580,-200.00000)
node(20100,157.88460,-200.00000)
node(20101,118.41350,-200.00000)
node(20102,78.94231,-200.00000)
node(20103,39.47115,-200.00000)
node(20104,0.00000,-200.00000)

print("Finished creating dashpot nodes...")

# define fixity of dashpot nodes
fix(20000,1,1)
fix(20001,1,1)
fix(20002,1,1)
fix(20003,1,1)
fix(20004,1,1)
fix(20005,1,1)
fix(20006,1,1)
fix(20007,1,1)
fix(20008,1,1)
fix(20009,1,1)
fix(20010,1,1)
fix(20011,1,1)
fix(20012,1,1)
fix(20013,1,1)
fix(20014,1,1)
fix(20015,1,1)
fix(20016,1,1)
fix(20017,1,1)
fix(20018,1,1)
fix(20019,1,1)
fix(20020,1,1)
fix(20021,1,1)
fix(20022,1,1)
fix(20023,1,1)
fix(20024,1,1)
fix(20025,1,1)
fix(20026,1,1)
fix(20027,1,1)
fix(20028,1,1)
fix(20029,1,1)
fix(20030,1,1)
fix(20031,1,1)
fix(20032,1,1)
fix(20033,1,1)
fix(20034,1,1)
fix(20035,1,1)
fix(20036,1,1)
fix(20037,1,1)
fix(20038,1,1)
fix(20039,1,1)
fix(20040,1,1)
fix(20041,1,1)
fix(20042,1,1)
fix(20043,1,1)
fix(20044,1,1)
fix(20045,1,1)
fix(20046,1,1)
fix(20047,1,1)
fix(20048,1,1)
fix(20049,1,1)
fix(20050,1,1)
fix(20051,1,1)
fix(20052,1,1)
fix(20053,1,1)
fix(20054,1,1)
fix(20055,1,1)
fix(20056,1,1)
fix(20057,1,1)
fix(20058,1,1)
fix(20059,1,1)
fix(20060,1,1)
fix(20061,1,1)
fix(20062,1,1)
fix(20063,1,1)
fix(20064,1,1)
fix(20065,1,1)
fix(20066,1,1)
fix(20067,1,1)
fix(20068,1,1)
fix(20069,1,1)
fix(20070,1,1)
fix(20071,1,1)
fix(20072,1,1)
fix(20073,1,1)
fix(20074,1,1)
fix(20075,1,1)
fix(20076,1,1)
fix(20077,1,1)
fix(20078,1,1)
fix(20079,1,1)
fix(20080,1,1)
fix(20081,1,1)
fix(20082,1,1)
fix(20083,1,1)
fix(20084,1,1)
fix(20085,1,1)
fix(20086,1,1)
fix(20087,1,1)
fix(20088,1,1)
fix(20089,1,1)
fix(20090,1,1)
fix(20091,1,1)
fix(20092,1,1)
fix(20093,1,1)
fix(20094,1,1)
fix(20095,1,1)
fix(20096,1,1)
fix(20097,1,1)
fix(20098,1,1)
fix(20099,1,1)
fix(20100,1,1)
fix(20101,1,1)
fix(20102,1,1)
fix(20103,1,1)
fix(20104,1,1)
print("Finished creating dashpot nodes and boundary conditions...")

#dashpot area
baseArea = 39.4715

#dashpot coefficient
dashpotCoeff_s = 7280.00
dashpotCoeff_p = 12740.00

# dashpot coefficient 
mCs    = baseArea*dashpotCoeff_s
mCp    = baseArea*dashpotCoeff_p

# dashpot material
# material
uniaxialMaterial('Viscous', 3, 0.5*mCs, 1.0)
uniaxialMaterial('Viscous', 4, mCs,  1.0)
uniaxialMaterial('Viscous', 5, mCp,  1.0)
uniaxialMaterial('Viscous', 6, 0.5*mCp,  1.0)

# define dashpot elements
element('zeroLength',17000,20000,18,'-mat',3,5,'-dir',1,2)
element('zeroLength',17001,20001,172,'-mat',4,6,'-dir',1,2)
element('zeroLength',17002,20002,147,'-mat',4,6,'-dir',1,2)
element('zeroLength',17003,20003,197,'-mat',4,6,'-dir',1,2)
element('zeroLength',17004,20004,146,'-mat',4,6,'-dir',1,2)
element('zeroLength',17005,20005,152,'-mat',4,6,'-dir',1,2)
element('zeroLength',17006,20006,145,'-mat',4,6,'-dir',1,2)
element('zeroLength',17007,20007,151,'-mat',4,6,'-dir',1,2)
element('zeroLength',17008,20008,144,'-mat',4,6,'-dir',1,2)
element('zeroLength',17009,20009,153,'-mat',4,6,'-dir',1,2)
element('zeroLength',17010,20010,143,'-mat',4,6,'-dir',1,2)
element('zeroLength',17011,20011,171,'-mat',4,6,'-dir',1,2)
element('zeroLength',17012,20012,142,'-mat',4,6,'-dir',1,2)
element('zeroLength',17013,20013,148,'-mat',4,6,'-dir',1,2)
element('zeroLength',17014,20014,141,'-mat',4,6,'-dir',1,2)
element('zeroLength',17015,20015,162,'-mat',4,6,'-dir',1,2)
element('zeroLength',17016,20016,140,'-mat',4,6,'-dir',1,2)
element('zeroLength',17017,20017,195,'-mat',4,6,'-dir',1,2)
element('zeroLength',17018,20018,139,'-mat',4,6,'-dir',1,2)
element('zeroLength',17019,20019,178,'-mat',4,6,'-dir',1,2)
element('zeroLength',17020,20020,138,'-mat',4,6,'-dir',1,2)
element('zeroLength',17021,20021,194,'-mat',4,6,'-dir',1,2)
element('zeroLength',17022,20022,137,'-mat',4,6,'-dir',1,2)
element('zeroLength',17023,20023,159,'-mat',4,6,'-dir',1,2)
element('zeroLength',17024,20024,136,'-mat',4,6,'-dir',1,2)
element('zeroLength',17025,20025,198,'-mat',4,6,'-dir',1,2)
element('zeroLength',17026,20026,135,'-mat',4,6,'-dir',1,2)
element('zeroLength',17027,20027,150,'-mat',4,6,'-dir',1,2)
element('zeroLength',17028,20028,134,'-mat',4,6,'-dir',1,2)
element('zeroLength',17029,20029,157,'-mat',4,6,'-dir',1,2)
element('zeroLength',17030,20030,133,'-mat',4,6,'-dir',1,2)
element('zeroLength',17031,20031,185,'-mat',4,6,'-dir',1,2)
element('zeroLength',17032,20032,132,'-mat',4,6,'-dir',1,2)
element('zeroLength',17033,20033,181,'-mat',4,6,'-dir',1,2)
element('zeroLength',17034,20034,131,'-mat',4,6,'-dir',1,2)
element('zeroLength',17035,20035,167,'-mat',4,6,'-dir',1,2)
element('zeroLength',17036,20036,130,'-mat',4,6,'-dir',1,2)
element('zeroLength',17037,20037,175,'-mat',4,6,'-dir',1,2)
element('zeroLength',17038,20038,129,'-mat',4,6,'-dir',1,2)
element('zeroLength',17039,20039,161,'-mat',4,6,'-dir',1,2)
element('zeroLength',17040,20040,128,'-mat',4,6,'-dir',1,2)
element('zeroLength',17041,20041,192,'-mat',4,6,'-dir',1,2)
element('zeroLength',17042,20042,127,'-mat',4,6,'-dir',1,2)
element('zeroLength',17043,20043,188,'-mat',4,6,'-dir',1,2)
element('zeroLength',17044,20044,126,'-mat',4,6,'-dir',1,2)
element('zeroLength',17045,20045,163,'-mat',4,6,'-dir',1,2)
element('zeroLength',17046,20046,125,'-mat',4,6,'-dir',1,2)
element('zeroLength',17047,20047,174,'-mat',4,6,'-dir',1,2)
element('zeroLength',17048,20048,124,'-mat',4,6,'-dir',1,2)
element('zeroLength',17049,20049,166,'-mat',4,6,'-dir',1,2)
element('zeroLength',17050,20050,123,'-mat',4,6,'-dir',1,2)
element('zeroLength',17051,20051,183,'-mat',4,6,'-dir',1,2)
element('zeroLength',17052,20052,122,'-mat',4,6,'-dir',1,2)
element('zeroLength',17053,20053,189,'-mat',4,6,'-dir',1,2)
element('zeroLength',17054,20054,121,'-mat',4,6,'-dir',1,2)
element('zeroLength',17055,20055,158,'-mat',4,6,'-dir',1,2)
element('zeroLength',17056,20056,120,'-mat',4,6,'-dir',1,2)
element('zeroLength',17057,20057,177,'-mat',4,6,'-dir',1,2)
element('zeroLength',17058,20058,119,'-mat',4,6,'-dir',1,2)
element('zeroLength',17059,20059,169,'-mat',4,6,'-dir',1,2)
element('zeroLength',17060,20060,118,'-mat',4,6,'-dir',1,2)
element('zeroLength',17061,20061,179,'-mat',4,6,'-dir',1,2)
element('zeroLength',17062,20062,117,'-mat',4,6,'-dir',1,2)
element('zeroLength',17063,20063,182,'-mat',4,6,'-dir',1,2)
element('zeroLength',17064,20064,116,'-mat',4,6,'-dir',1,2)
element('zeroLength',17065,20065,155,'-mat',4,6,'-dir',1,2)
element('zeroLength',17066,20066,115,'-mat',4,6,'-dir',1,2)
element('zeroLength',17067,20067,156,'-mat',4,6,'-dir',1,2)
element('zeroLength',17068,20068,114,'-mat',4,6,'-dir',1,2)
element('zeroLength',17069,20069,184,'-mat',4,6,'-dir',1,2)
element('zeroLength',17070,20070,113,'-mat',4,6,'-dir',1,2)
element('zeroLength',17071,20071,180,'-mat',4,6,'-dir',1,2)
element('zeroLength',17072,20072,112,'-mat',4,6,'-dir',1,2)
element('zeroLength',17073,20073,168,'-mat',4,6,'-dir',1,2)
element('zeroLength',17074,20074,111,'-mat',4,6,'-dir',1,2)
element('zeroLength',17075,20075,176,'-mat',4,6,'-dir',1,2)
element('zeroLength',17076,20076,110,'-mat',4,6,'-dir',1,2)
element('zeroLength',17077,20077,160,'-mat',4,6,'-dir',1,2)
element('zeroLength',17078,20078,109,'-mat',4,6,'-dir',1,2)
element('zeroLength',17079,20079,191,'-mat',4,6,'-dir',1,2)
element('zeroLength',17080,20080,108,'-mat',4,6,'-dir',1,2)
element('zeroLength',17081,20081,187,'-mat',4,6,'-dir',1,2)
element('zeroLength',17082,20082,107,'-mat',4,6,'-dir',1,2)
element('zeroLength',17083,20083,164,'-mat',4,6,'-dir',1,2)
element('zeroLength',17084,20084,106,'-mat',4,6,'-dir',1,2)
element('zeroLength',17085,20085,173,'-mat',4,6,'-dir',1,2)
element('zeroLength',17086,20086,105,'-mat',4,6,'-dir',1,2)
element('zeroLength',17087,20087,165,'-mat',4,6,'-dir',1,2)
element('zeroLength',17088,20088,104,'-mat',4,6,'-dir',1,2)
element('zeroLength',17089,20089,186,'-mat',4,6,'-dir',1,2)
element('zeroLength',17090,20090,103,'-mat',4,6,'-dir',1,2)
element('zeroLength',17091,20091,190,'-mat',4,6,'-dir',1,2)
element('zeroLength',17092,20092,102,'-mat',4,6,'-dir',1,2)
element('zeroLength',17093,20093,196,'-mat',4,6,'-dir',1,2)
element('zeroLength',17094,20094,101,'-mat',4,6,'-dir',1,2)
element('zeroLength',17095,20095,154,'-mat',4,6,'-dir',1,2)
element('zeroLength',17096,20096,100,'-mat',4,6,'-dir',1,2)
element('zeroLength',17097,20097,199,'-mat',4,6,'-dir',1,2)
element('zeroLength',17098,20098,99,'-mat',4,6,'-dir',1,2)
element('zeroLength',17099,20099,149,'-mat',4,6,'-dir',1,2)
element('zeroLength',17100,20100,98,'-mat',4,6,'-dir',1,2)
element('zeroLength',17101,20101,170,'-mat',4,6,'-dir',1,2)
element('zeroLength',17102,20102,97,'-mat',4,6,'-dir',1,2)
element('zeroLength',17103,20103,193,'-mat',4,6,'-dir',1,2)
element('zeroLength',17104,20104,2,'-mat',3,5,'-dir',1,2)
print("Finished creating dashpot material and elements...")

#-----------------------------------------------------------------------------------------
#  6. CREATE GRAVITY RECORDERS
#-----------------------------------------------------------------------------------------

# record nodal displacment, velocities and acceleration at each time step
#recorder('Node', '-file', 'Gdisplacement_Dunedin_stBeach_Hyde.out', '-time',  '-nodeRange', 1, 17303, '-dof', 1, 2,  'disp')
#recorder('Node', '-file', 'Gvelocity_Dunedin_stBeach_Hyde.out'    , '-time',  '-nodeRange', 1, 17303, '-dof', 1, 2,  'vel')
#recorder('Node', '-file', 'Gacceleration_Dunedin_stBeach_Hyde.out', '-time',  '-nodeRange', 1, 17303, '-dof', 1, 2,  'accel')
#record nodal reaction
#recorder('Node', '-file', 'Greaction_Dunedin_stBeach_Hyde.out', '-time',  '-nodeRange', 1,17303, '-dof', 1, 2,  'reaction')

# record stress and strain at each gauss point in the soil elements
#recorder('Element', '-file' ,'Gstress1_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 1, 'stress')
#recorder('Element', '-file' ,'Gstress2_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 2, 'stress')
#recorder('Element', '-file' ,'Gstress3_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 3, 'stress')
#recorder('Element', '-file' ,'Gstress4_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 4, 'stress')
#recorder('Element', '-file' ,'Gstress5_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 5, 'stress')
#recorder('Element', '-file' ,'Gstress6_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 6, 'stress')

#recorder('Element', '-file' ,'Gstrain1_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 1, 'strain')
#recorder('Element', '-file' ,'Gstrain2_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 2, 'strain')
#recorder('Element', '-file' ,'Gstrain3_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 3, 'strain')
#recorder('Element', '-file' ,'Gstrain4_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 4, 'strain')
#recorder('Element', '-file' ,'Gstrain5_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834,   'material', 5, 'strain')
#recorder('Element', '-file' ,'Gstrain6_Dunedin_stBeach_Hyde.out',   '-time',  '-eleRange',  1,   16834       ,   'material', 6, 'strain')
print("Finished creating gravity recorders...")

#-----------------------------------------------------------------------------------------
#  7. DEFINE ANALYSIS PARAMETERS
#-----------------------------------------------------------------------------------------

#---GROUND MOTION PARAMETERS
# time step in ground motion record
motionDT        = 0.005
# number of steps in ground motion record
motionSteps     = 18000

#---WAVELENGTH PARAMETERS
# highest frequency desired to be well resolved (Hz)
## fMax     = 100.0
# shear wave velocity desired to be well resolved
## vsMin    = Vs[6]
# wavelength of highest resolved frequency
## wave     = vsMin/fMax
# number of elements per one wavelength
## nEle     = 8

#---RAYLEIGH DAMPING PARAMETERS
# damping ratio
damp    = 0.01
# lower frequency
omega1  = 2 * pi * 0.0
# upper frequency
omega2  = 2 * pi * 15
# damping coefficients
a0      = 2 * damp * omega1 * omega2 / (omega1 + omega2)
a1      = 2 * damp / (omega1 + omega2)
print(f"damping coefficients: a_0 = {a0};  a_1 = {a1}")

dt = motionDT

#---ANALYSIS PARAMETERS
# Newmark parameters
gamma           = 0.5
beta            = 0.25

print("Finished defining analysis parameters...")

#-----------------------------------------------------------------------------------------------------------
#  8.GRAVITY ANALYSIS
#-----------------------------------------------------------------------------------------------------------

# update materials to consider plastic behavior
#for j in iter(numLayers):
#    updateMaterialStage('-material', j, '-stage', 0)

updateMaterialStage('-material', 1, '-stage', 0)
updateMaterialStage('-material', 2, '-stage', 0)
updateMaterialStage('-material', 3, '-stage', 0)
updateMaterialStage('-material', 4, '-stage', 0)
updateMaterialStage('-material', 5, '-stage', 0)
updateMaterialStage('-material', 6, '-stage', 0)

#constraints ('Penalty', alphaS = 1.e18, alphaM = 1.e18 )
constraints ('Penalty',  1.e18,  1.e18 )
test        ('NormDispIncr', 1e-6, 35, 1)
algorithm   ('KrylovNewton')
numberer    ('RCM')
system      ('ProfileSPD')
integrator  ('Newmark', gamma, beta)
rayleigh    (0.87, 0.022, 0.0, 0.0)
analysis    ('Transient')

startT = datetime.datetime.now()
analyze     (10, 5.0e2)
analyze     (10, 5.0e3)

print("Finished with elastic gravity analysis...")

# update materials to consider plastic behavior
#for j in iter(numLayers):
#    updateMaterialStage('-material', j, '-stage', 1)
updateMaterialStage('-material', 1, '-stage', 1)
updateMaterialStage('-material', 2, '-stage', 1)
updateMaterialStage('-material', 3, '-stage', 1)
updateMaterialStage('-material', 4, '-stage', 1)
updateMaterialStage('-material', 5, '-stage', 1)
updateMaterialStage('-material', 6, '-stage', 1)

# plastic gravity loading
#analyze     (40, 5.0e2)
analyze     (10, 5.0e-3)
analyze     (10, 1.0e-2)
analyze     (10, 5.0e-2)
analyze     (10, 1.0e-1)
analyze     (10, 5.0e-1)
analyze     (10, 5.0e0)
analyze     (10, 5.0e1)
analyze     (10, 5.0e10)

#Import the boundary nodes from text file
nodesBdry = []
print("Reading nodes from file")

with open("nodes_boundary_St_Beach_revised.txt") as file:
	p = list(file.read().splitlines())
	p= list(map(int,p))
	
NnodesBdry = len(nodesBdry)
#Define list for nodal reaction:
React1 = []
React2 = []
	
#Loop over boundary nodes; obtain current nodal reaction;
#Release the nodal fixity; apply the nodal reaction
for i in range(0,NnodesBdry):
	tmp_React = nodeReaction(nodesBdry[i])
	#Release the nodal fixity:
	remove('sp',nodesBdry [i],1)
	remove('sp',nodesBdry [i],2)
	#store the nodal reaction:
	React1.append(tmp_React[0])
	React2.append(tmp_React[1])
	
timeSeries('Constant',1001,'-factor', 1.0)
pattern ('Plain',1001, 1001, '-fact', 1.0)
for i in range(0,NnodesBdry):
	load(nodesBdry[i],React1[i],React2[i])
	
print("Base reaction created...")

analyze(10,5.0e-1)
analyze(10,5.0e0)
analyze(10,5.0e1)

print("Finished with plastic gravity analysis...")

#-----------------------------------------------------------------------------------------------------------
#  10. CREATE POST-GRAVITY RECORDERS
#-----------------------------------------------------------------------------------------------------------

# reset time and analysis
#setTime (0.0)
wipeAnalysis()

# record nodal displacments, velocities, and accelerations at each time step
#recorder('Node', '-file', 'displacement_Dunedin_stBeach_Hyde.out', '-time', '-dT', motionDT, '-nodeRange', 1, 17303, '-dof', 1, 2,  'disp')
#recorder('Node', '-file', 'velocity_Dunedin_stBeach_Hyde.out',     '-time', '-dT', motionDT, '-nodeRange', 1, 17303, '-dof', 1, 2,  'vel')

#recorder('Node', '-file', '/home/kowal/Desktop/2D_test/acceleration_Dunedin_StBeach_Hyde.out', '-time', '-dT', motionDT, '-nodeRange', 1, 17303, '-dof', 1, 2,  'accel')
recorder ('Node', '-file', 'acceleration_Dunedin_StBeach_90s_Hyde_GMSstation.out', '-time', '-dt', motionDT,'-node', 559, '-dof', 1, 2,  'accel')
recorder ('Node', '-file', 'acceleration_Dunedin_StBeach_90s_Hyde_node1.out', '-time', '-dt', motionDT, '-node', 1, '-dof', 1, 2,  'accel')
recorder ('Node', '-file', 'acceleration_Dunedin_StBeach_90s_Hyde_node89.out' '-time', 'dt', motionDT, '-node', 89, '-dof', 1, 2,  'accel')
recorder ('Node', '-file', 'acceleration_Dunedin_StBeach_90s_Hyde_node9.out', '-time', '-dt', motionDT, '-node', 9, '-dof', 1, 2,  'accel')
recorder ('Node', '-file', 'acceleration_Dunedin_StBeach_90s_Hyde_node15.out', '-time', '-dt', motionDT, '-node', 15, '-dof', 1, 2,  'accel')
recorder ('Node', '-file', 'acceleration_Dunedin_StBeach_90s_Hyde_node19.out', '-time', '-dt', motionDT, '-node', 19, '-dof', 1, 2,  'accel')

# record stress and strain at each gauss point in the soil elements
#recorder('Element', '-file' ,'stress1_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 1, 'stress')
#recorder('Element', '-file' ,'stress2_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 2, 'stress')
#recorder('Element', '-file' ,'stress3_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 3, 'stress')
#recorder('Element', '-file' ,'stress4_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 4, 'stress')
#recorder('Element', '-file' ,'stress5_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 5, 'stress')
recorder('Element', '-file' ,'stress6_Dunedin_stBeach_90s_Hyde.out',   '-time', '-dt', motionDT, '-eleRange',  13184,   16834,   'material', 6, 'stress')

#recorder('Element', '-file' ,'strain1_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 1, 'strain')
#recorder('Element', '-file' ,'strain2_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 2, 'strain')
#recorder('Element', '-file' ,'strain3_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 3, 'strain')
#recorder('Element', '-file' ,'strain4_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 4, 'strain')
#recorder('Element', '-file' ,'strain5_Dunedin_stBeach_Hyde.out',   '-time', '-dT', motionDT, '-eleRange',  1,   16834,   'material', 5, 'strain')
recorder('Element', '-file' ,'strain6_Dunedin_stBeach_90s_Hyde.out',   '-time', '-dt', motionDT, '-eleRange',  13184,   16834,   'material', 6, 'strain')

print("Finished creating all recorders...")

#-----------------------------------------------------------------------------------------------------------
#  11. DYNAMIC ANALYSIS  
#-----------------------------------------------------------------------------------------------------------

# define constant factor for applied velocity
cFactor = 1.00*baseArea*dashpotCoeff_s
#cFactorV = 0.00*baseArea*dashpotCoeff_p

# define velocity time history file
velocityFile = '/home/kowal/2D_site_analysis_20202501/deconvolution/Dunedin_St_Beach/deconv_03_a_Dunedin_St_Beach_Hyde_NS_shorter.out'
#velocityFileV = '/home/kowal/2D_site_analysis_20202501/deconvolution/Dunedin_St_Beach/deconv_01_a_Dunedin_St_Beach_Akatore_NS_70s.out'

# timeseries object for force history
#timeSeries('Path', 1, '-filePath', velocityFile, '-dt', motionDT, '-factor', cFactor)
# timeseries object for force history
timeSeries('Path', 10, '-filePath', velocityFile, '-dt', motionDT, '-factor', cFactor)
#timeSeries('Path', 1, '-filePath', velocityFile, '-dt', motionDT, '-factor', cFactor)
#timeSeries('Path', 10, '-filePath', velocityFile, '-dT', motionDT, '-factor', cFactorH)
#timeSeriesH('Path', 10, '-filePath', velocityFile, '-dt', motionDT, '-factor', cFactor)
#timeSeries('Path', 11, '-filePath', velocityFile, '-dT', motionDT, '-factor', cFactorV)
cFactor = 0.00*baseArea*dashpotCoeff_p
timeSeries('Path', 11, '-filePath', velocityFile, '-dt', motionDT, '-factor', cFactor)
#timeSeriesV('Path', 11, '-filePath', velocityFileV, '-dt', motionDT, '-factor', cFactorV)

# loading object
pattern('Plain', 10, 10, '-fact', 1.0)
load(2,0.50,0.00)
load(18,0.50,0.00)
load(98,1.00,0.00)
load(99,1.00,0.00)
load(100,1.00,0.00)
load(101,1.00,0.00)
load(102,1.00,0.00)
load(103,1.00,0.00)
load(104,1.00,0.00)
load(105,1.00,0.00)
load(106,1.00,0.00)
load(107,1.00,0.00)
load(108,1.00,0.00)
load(109,1.00,0.00)
load(110,1.00,0.00)
load(111,1.00,0.00)
load(112,1.00,0.00)
load(113,1.00,0.00)
load(114,1.00,0.00)
load(115,1.00,0.00)
load(116,1.00,0.00)
load(117,1.00,0.00)
load(118,1.00,0.00)
load(119,1.00,0.00)
load(120,1.00,0.00)
load(121,1.00,0.00)
load(122,1.00,0.00)
load(123,1.00,0.00)
load(124,1.00,0.00)
load(125,1.00,0.00)
load(126,1.00,0.00)
load(127,1.00,0.00)
load(128,1.00,0.00)
load(129,1.00,0.00)
load(130,1.00,0.00)
load(131,1.00,0.00)
load(132,1.00,0.00)
load(133,1.00,0.00)
load(134,1.00,0.00)
load(135,1.00,0.00)
load(136,1.00,0.00)
load(137,1.00,0.00)
load(138,1.00,0.00)
load(139,1.00,0.00)
load(140,1.00,0.00)
load(141,1.00,0.00)
load(142,1.00,0.00)
load(143,1.00,0.00)
load(144,1.00,0.00)
load(145,1.00,0.00)
load(146,1.00,0.00)
load(147,1.00,0.00)
load(148,1.00,0.00)
load(149,1.00,0.00)
load(150,1.00,0.00)
load(151,1.00,0.00)
load(152,1.00,0.00)
load(153,1.00,0.00)
load(154,1.00,0.00)
load(155,1.00,0.00)
load(156,1.00,0.00)
load(157,1.00,0.00)
load(158,1.00,0.00)
load(159,1.00,0.00)
load(160,1.00,0.00)
load(161,1.00,0.00)
load(162,1.00,0.00)
load(163,1.00,0.00)
load(164,1.00,0.00)
load(165,1.00,0.00)
load(166,1.00,0.00)
load(167,1.00,0.00)
load(168,1.00,0.00)
load(169,1.00,0.00)
load(170,1.00,0.00)
load(171,1.00,0.00)
load(172,1.00,0.00)
load(173,1.00,0.00)
load(174,1.00,0.00)
load(175,1.00,0.00)
load(176,1.00,0.00)
load(177,1.00,0.00)
load(178,1.00,0.00)
load(179,1.00,0.00)
load(180,1.00,0.00)
load(181,1.00,0.00)
load(182,1.00,0.00)
load(183,1.00,0.00)
load(184,1.00,0.00)
load(185,1.00,0.00)
load(186,1.00,0.00)
load(187,1.00,0.00)
load(188,1.00,0.00)
load(189,1.00,0.00)
load(190,1.00,0.00)
load(191,1.00,0.00)
load(192,1.00,0.00)
load(193,1.00,0.00)
load(194,1.00,0.00)
load(195,1.00,0.00)
load(196,1.00,0.00)
load(197,1.00,0.00)
load(198,1.00,0.00)
load(199,1.00,0.00)
load(200,1.00,0.00)

pattern('Plain', 11, 11, '-fact', 1.0)
load(2,0.00,0.50)
load(18,0.00,0.50)
load(98,0.00,1.00)
load(99,0.00,1.00)
load(100,0.00,1.00)
load(101,0.00,1.00)
load(102,0.00,1.00)
load(103,0.00,1.00)
load(104,0.00,1.00)
load(105,0.00,1.00)
load(106,0.00,1.00)
load(107,0.00,1.00)
load(108,0.00,1.00)
load(109,0.00,1.00)
load(110,0.00,1.00)
load(111,0.00,1.00)
load(112,0.00,1.00)
load(113,0.00,1.00)
load(114,0.00,1.00)
load(115,0.00,1.00)
load(116,0.00,1.00)
load(117,0.00,1.00)
load(118,0.00,1.00)
load(119,0.00,1.00)
load(120,0.00,1.00)
load(121,0.00,1.00)
load(122,0.00,1.00)
load(123,0.00,1.00)
load(124,0.00,1.00)
load(125,0.00,1.00)
load(126,0.00,1.00)
load(127,0.00,1.00)
load(128,0.00,1.00)
load(129,0.00,1.00)
load(130,0.00,1.00)
load(131,0.00,1.00)
load(132,0.00,1.00)
load(133,0.00,1.00)
load(134,0.00,1.00)
load(135,0.00,1.00)
load(136,0.00,1.00)
load(137,0.00,1.00)
load(138,0.00,1.00)
load(139,0.00,1.00)
load(140,0.00,1.00)
load(141,0.00,1.00)
load(142,0.00,1.00)
load(143,0.00,1.00)
load(144,0.00,1.00)
load(145,0.00,1.00)
load(146,0.00,1.00)
load(147,0.00,1.00)
load(148,0.00,1.00)
load(149,0.00,1.00)
load(150,0.00,1.00)
load(151,0.00,1.00)
load(152,0.00,1.00)
load(153,0.00,1.00)
load(154,0.00,1.00)
load(155,0.00,1.00)
load(156,0.00,1.00)
load(157,0.00,1.00)
load(158,0.00,1.00)
load(159,0.00,1.00)
load(160,0.00,1.00)
load(161,0.00,1.00)
load(162,0.00,1.00)
load(163,0.00,1.00)
load(164,0.00,1.00)
load(165,0.00,1.00)
load(166,0.00,1.00)
load(167,0.00,1.00)
load(168,0.00,1.00)
load(169,0.00,1.00)
load(170,0.00,1.00)
load(171,0.00,1.00)
load(172,0.00,1.00)
load(173,0.00,1.00)
load(174,0.00,1.00)
load(175,0.00,1.00)
load(176,0.00,1.00)
load(177,0.00,1.00)
load(178,0.00,1.00)
load(179,0.00,1.00)
load(180,0.00,1.00)
load(181,0.00,1.00)
load(182,0.00,1.00)
load(183,0.00,1.00)
load(184,0.00,1.00)
load(185,0.00,1.00)
load(186,0.00,1.00)
load(187,0.00,1.00)
load(188,0.00,1.00)
load(189,0.00,1.00)
load(190,0.00,1.00)
load(191,0.00,1.00)
load(192,0.00,1.00)
load(193,0.00,1.00)
load(194,0.00,1.00)
load(195,0.00,1.00)
load(196,0.00,1.00)
load(197,0.00,1.00)
load(198,0.00,1.00)
load(199,0.00,1.00)
load(200,0.00,1.00)
print("Dynamic loading created...")

print(f"number of steps in analysis: {motionSteps}")
print(f"analysis time step: {motionDT }")

# analysis objects
constraints ('Penalty',  1.e18,  1.e18 )
print(f"analysis constrains defined")
test        ('NormDispIncr', 1e-3, 35, 1)
print(f"analysis test defined")
algorithm   ('KrylovNewton')
print(f"analysis algorithm defined")
numberer    ('RCM')
print(f"analysis numberer defined")
system      ('ProfileSPD')
print(f"analysis system defined")
integrator  ('Newmark', gamma, beta)
print(f"integrator defined")
rayleigh    (a0, a1, 0.0, 0.0)
print(f"analysis rayleigh damping defined")
analysis    ('Transient')
print(f"analysis constrains defined")

analyze     (motionSteps,  motionDT)

print(f"analyze performed")

#ok = analyze (motionSteps,  motionDT)
ok = analyze (motionSteps,  motionDT)
if ok != 0:
	print("did not converge, reducing time step")
	curTime = getTime()
	print("curTime:{0}".format(curTime))
	curStep=curTime/motionDT 
	print("curStep:{0}".format(curStep))
	remStep=int((remStep-curStep)*2.0)
	print("remStep:{0}".format(remStep))
	dt=dt/2.0
	print("dT:{0}".format(motionDT ))
	analyze(remStep,motionDT )
	
	
endT = datetime.datetime.now()
print("Finished with dynamic analysis...")
print("Analysis execution time:{0}".format(str(endt-startT)))

wipe()

wipe()
Top
mhscott
Posts: 880
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:
Contact mhscott

Re: Segmentation fault running dynamic analysis

Post by mhscott »

Please email the file to me (michael.scott@oregonstate.edu).

Note: please change the extension from .py to .txt before attaching.
Top
mhscott
Posts: 880
Joined: Tue Jul 06, 2004 3:38 pm
Location: Corvallis, Oregon USA
Contact:
Contact mhscott

Re: Segmentation fault running dynamic analysis

Post by mhscott »

A few things I noticed

The ProfileSPD solver is VERY slow. Use Mumps instead, it's like 20x faster than ProfileSPD on your model.

Code: Select all

system('Mumps')
Your penalty numbers are VERY large

Code: Select all

constraints ('Penalty',  1.e18,  1.e18 )
You have many zero length elements with non-zero length

Code: Select all

WARNING ZeroLength::setDomain(): Element 17100 has L= 78.9423, which is greater than the tolerance
WARNING ZeroLength::setDomain(): Element 17101 has L= 2131.44, which is greater than the tolerance
WARNING ZeroLength::setDomain(): Element 17102 has L= 244.117, which is greater than the tolerance
WARNING ZeroLength::setDomain(): Element 17103 has L= 3078.75, which is greater than the tolerance
There is a warning with material stage parameter

Code: Select all

WARNING: MaterialStageParameter::setDomain() - no effect with material tag 1
The segmentation fault is from the time series with the load pattern attempting to read past the values in the time series file because the time in the domain is very large. You should uncomment the setTime command.

Code: Select all

# reset time and analysis
setTime (0.0)
wipeAnalysis()
Top
Post Reply

Return to “OpenSeesPy”