geosharma/ex_g_02_siteresponse_totalstress_tcl2py.py

## ex_g_02_siteresponse_totalstress_tcl2py.py
#!/usr/bin/env python
# -*-coding:utf-8 -*-

import numpy as np
import openseespy.opensees as ops
from pathlib import Path

# Example problem:
# Description: 1D Total Stress Site Response Analysis, Single soil layer
# Original authors of this example: Christopher McGann and Pedro Arduino, University of Washington
# ref: https://opensees.berkeley.edu/wiki/index.php/Site_Response_Analysis_of_a_Layered_Soil_Column_(Total_Stress_Analysis)
# Units: kN, m, sec


def roundup(num):
    """Check if a float is a whole number if not round up"""
    return int(num) if num.is_integer() else int(np.ceil(num))

def data_recorders(fpath, fname_prefix, number_of_nodes, number_of_elements, motiondt):
    # this is a deviation from the original solution, since only ground
    # accelerations, velocity and displacements are considered, for now only ground
    # node values will be recorded

    # fmt: off
    # recorders for nodal displacements, velocity and acceleration at each time step
    restypes = ["disp", "vel", "accel"]
    for restype in restypes:
        fname = fname_prefix + restype + ".out"
        print(f"Filename: {fname}")
        fout = str(fpath.joinpath(fname))
        ops.recorder(
        "Node",
        "-file", fout,
        "-time",
        "-dT", motiondt,
        "-nodeRange", 1, number_of_nodes,
        "-dof", *[1, 2],
        restype,
    )

    # # ref: http://opensees.berkeley.edu/wiki/index.php/PressureIndependMultiYield_Material
    # # For 2D problems, the stress output follows this order: σxx, σyy, σzz, σxy, ηr, where ηr
    # # is the ratio between the shear (deviatoric) stress and peak shear strength at the current
    # # confinement (0<=ηr<=1.0). The strain output follows this order: εxx, εyy, γxy.
    # stresses and strains at each Gauss point in the soil elements
    gausspoints = ["1", "2", "3", "4"]
    for gpt in gausspoints:
        fname = fname_prefix + "_stress" + gpt + ".out"
        fout = str(fpath.joinpath(fname))
        ops.recorder(
            "Element",
            "-file", fout,
            "-time",
            "-dT", motiondt,
            "-eleRange", 1, number_of_elements,
            "material", gpt,
            "stress",
        )
        fname = fname_prefix + "_strain" + gpt + ".out"
        fout = str(fpath.joinpath(fname))
        ops.recorder(
            "Element",
            "-file", fout,
            "-time",
            "-dT", motiondt,
            "-eleRange", 1, number_of_elements,
            "material", gpt,
            "strain",
        )
    # fmt: on

def main():
    # file path and filenames
    path_root = Path(r"./")
    path_data = path_root / "data"

    # general filename for recorder outputs
    fname_prefix = "ex_g_02_siteresponse_totalstress_"

    # velocity time history file
    velfile = "velocityhistory.out"
    fvel = str(path_data.joinpath(velfile))
    print(f"Ground motion file path: {fvel}")
    fvel = "data/velocityHistory.out"

    # 1. Define analysis parameters

    # soil profile, thickness of the soil layer [m]
    depth_soillayer = 40.0

    # acceleration due to gravity [m/s2]
    accg = 9.81

    # material properties
    # # soil saturated density [Mg/m3], or soil mass density
    soil_satdensity = 1.7
    # soil shear wave velocity [m/s]
    soil_vs = 250.0
    # Poisson's ratio of the soil
    nu = 0.0
    # the undrained shear strength, c of Mohr-Coulomb failure envelop [kPa]
    cohesion = 95.0
    # soil friction angle, phi of Mohr-Coulomb failure envelop [deg]
    friction_angle = 0.0
    # peak shear strain
    peak_shear_strain = 0.05
    # reference pressure [kPa]
    reference_pressure = 80.0
    # pressure dependency coefficient
    pressure_dependence = 0.0
    # bedrock shear wave velocity [m/s]
    rock_vs = 760.0
    # bedrock mass density [Mg/m^3]
    rock_density = 2.4
    # soil shear modulus [kPa]
    modG = soil_satdensity * soil_vs ** 2
    # soil elastic modulus [kPa]
    modE = 2 * modG * (1 + nu)
    # soil bulk modulus [kPa]
    modK = modE / (3 * (1 - 2 * nu))
    # number of yield surfaces
    nsurf = 25

    print(f"Soil elastic properites:")
    print(f"G: {modG/1000:6.1f} MPa, E: {modE/1000:6.1f} MPa, K: {modK/1000:6.1f} MPa")

    # ground motion paramters
    # time step in ground motion record [s]
    gm_dt = 0.005
    # number of steps in ground motion record
    # this depends on the number of points in the ground motion file
    gm_nsteps = 7990

    # wavelength parameters
    # highest frequency desired to be well resolved [Hz]
    fmax = 100.0
    # wavelength of the highest resolved frequency
    lamdamin = soil_vs / fmax
    # number of elements per one wavelength
    nele_wave = 10

    # rayleigh damping parameters (read on this)
    # damping ratio
    damping_ratio = 0.02
    # lower frequency [Hz]
    f_lower = 0.2
    omega_lower = 2 * np.pi * f_lower
    # upper frequency [Hz]
    f_upper = 20
    omega_upper = 2 * np.pi * f_upper
    # rayleigh damping coefficients
    alpham = 2.0 * damping_ratio * omega_lower * omega_upper / (omega_lower + omega_upper)
    betak = 2.0 * damping_ratio / (omega_lower + omega_upper)
    betakinit = 0.0
    betakcomm = 0.0
    print(f"Damping coefficients: alpham = {alpham:6.3f}, betak = {betak:6.5f}")

    # Newmark parameters
    newmark_gamma = 0.5
    newmark_beta = 0.25

    # acceleration due to gravity in the x and y direction
    accg_x = 0.0
    accg_y = -accg
    wgt_x = 0.0
    wgt_y = accg_y * soil_satdensity

    # 2. define mesh geometry
    # number of element in x and y direction
    nx = 1
    # trial for vertical element size
    # initial element size
    h_trial = lamdamin / nele_wave
    # if the number of elements is not an integer round up
    ny = roundup(depth_soillayer / h_trial)
    print(f"Number of elements: nx = {nx}, ny = {ny}")

    # size of elements in x and y direction of the nodes
    ly = float(depth_soillayer / ny)
    lx = ly
    print(f"Element size: lx = {lx:6.3f} m, ly = {ly:6.3f} m")

    # soil element thickness
    thickness = 1.0
    # analysis type, "PlaneStrain" or "PlaneStress"
    analysis_type = "PlaneStrain"
    # surface pressure on elements
    surface_pressure = 0.0

    # start node number
    start_node_number = 1

    # number of nodes in the vertical direction
    num_node_y = 2 * (ny + 1)

    # wipe previous analyses
    ops.wipe()

    # corner nodes are 2D with 2 DOF (displacements in x1 and x2)
    ops.model("basic", "-ndm", 2, "-ndf", 2)

    ycoord = 0.0
    count = 0

    # loop over nodes
    a = range(1, num_node_y + 1, 2)
    for j in range(1, num_node_y, 2):
        ycrd = ycoord + count * ly
        ops.node(j, 0.0, ycrd)
        ops.node(j+1, lx, ycrd)
        print(f"Node: {j}, x: 0.0,  y: {ycrd}")
        print(f"Node: {j + 1}, x: {lx}, y: {ycrd}")
        count += 1

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

    # 4. Define dashpot nodes
    ops.node(2000, 0.0, 0.0)
    ops.node(2001, 0.0, 0.0)
    print("Finished creating dashpot nodes....")

    # 5. Define boundary conditions and equalDOF
    # definie fixity of the base nodes
    ops.fix(1, 0, 1)
    ops.fix(2, 0, 1)

    # define fixity of the dashpot nodes
    ops.fix(2000, 1, 1)
    ops.fix(2001, 0, 1)

    # define equalDOF for simple shear deformation of soil elements
    for k in range(3, num_node_y, 2):
        ops.equalDOF(k, k + 1, 1, 2)
        print(f"equalDOF: rnode: {k}, cnode: {k + 1}, dofs: 1 2")

    # define equalDOF for dashpot and base soil nodes
    ops.equalDOF(1, 2, 1)
    ops.equalDOF(1, 2001, 1)

    print(f"Created all boundary conditions and equalDOF:")

    # 6. Define soil materials
    soil_mattag = 1
    # nd = 2 for plane strain and 3 for 3-D
    nd = 2
    # number of dimensions (nd) 2 for plain-strain 3 for 3D
    ops.nDMaterial(
        "PressureIndependMultiYield",
        soil_mattag,
        nd,
        soil_satdensity,
        modG,
        modK,
        cohesion,
        peak_shear_strain,
        friction_angle,
        reference_pressure,
        pressure_dependence,
        nsurf,
    )
    print(f"Soil material created")

    element_mass_density = 0.0
    print(
        "{0:5s}, {1:6s}, {2:6s}, {3:6s}, {4:6s}".format(
            "Ele", "Node i", "Node j", "Node k", "Node l"
        ))

    for j in range(1, ny + 1):
        nodei = 2 * j - 1
        nodej = nodei + 1
        nodek = nodei + 3
        nodel = nodei + 2
        ops.element("quad", j, nodei, nodej, nodek, nodel, 1.0, "PlaneStrain", 1.0, 0.0, 0.0, wgt_x, wgt_y)
        print("{0:5d}, {1:6d}, {2:6d}, {3:6d}, {4:6d}".format(j, nodei, nodej, nodek, nodel))

    print(f"All soil elements created")

    # material and elements for viscous damping
    # dashpot coefficients, again figure out how this coefficient comes about
    dashpot_c = lx * rock_density * rock_vs
    # for linear damping, power factor = 1
    dashpot_alpha = 1.0
    # dashpot material
    dashpot_mattag = 4000
    ops.uniaxialMaterial("Viscous", 4000, dashpot_c, 1.0)
    # element, connecting two dashpot nodes
    dashpot_elenum = 5000

    ops.element(
        "zeroLength",
        5000,
        *[2000, 2001],
        "-mat",
        4000,
        "-dir",
        *[1],
    )

    print("Finished creating dashpot material and element...")
    print(f"Gravity analysis:")
    # update material stage to 0 for elastic behavior
    # ops.updateMaterialStage("-material", soil_mattag, "-stage", 0)
    print(f"Material updated for elastic behavior during gravity loading")

    # gravity loading
    ops.constraints("Transformation")
    ops.test("NormDispIncr", 1e-5, 30, 1)
    ops.algorithm("Newton")
    ops.numberer("RCM")
    ops.system("ProfileSPD")
    ops.integrator("Newmark", newmark_gamma, newmark_beta)
    ops.analysis("Transient")
    ops.analyze(10, 500.0)
    print(f"Elastic gravity analysis over")

    # update material stage to include plastic analysis
    ops.updateMaterialStage("-material", soil_mattag, "-stage", 1)
    ops.analyze(40, 500.0)
    print(f"Plastic gravity analysis over")
    ops.setTime(0.0)
    ops.wipeAnalysis()


    print(f"Create data recorders")
    # specify the node and element at which the acceleration, velocity, displacement,
    # stress and strain time histories are required, for surface node
    last_node_number = 322
    last_element_number = 160
    data_recorders(
        path_data, fname_prefix, last_node_number, last_element_number, gm_dt
    )
    print(f"Finished creating data recorders")

    # destory previous analysis object and start consolidation
    print(
        f"Wiping gravity analysis, resetting time to 0.0, and starting dynamic analysis:"
    )


    # define constant factor for applied velocity, find out how this comes about?
    # this is what I found, the link might not work the next time
    # http://www.bgu.ac.il/geol/hatzor/articles/Bao_Hatzor_Huang_2012.pdf
    cfactor = dashpot_c

    # load pattern and timeseries for applied force history
    print(f"Ground motion file: {fvel}, dt: {gm_dt}, cfactor: {cfactor}")

    # load timeseries tag
    tstag_eq = 20
    # load pattern tag
    pattag_eq = 10
    # node at which the excitation is applied
    nodetag_eq = 1
    loadvals = [1.0, 0.0]
    ops.timeSeries(
        "Path", tstag_eq, "-dt", gm_dt, "-filePath", fvel, "-factor", cfactor
    )

    ops.pattern("Plain", pattag_eq, tstag_eq)
    ops.load(nodetag_eq, *loadvals)

    # determine the analysis time  step using CFL condition
    # CFL stands for Courant-Friedrichs-Levy
    # ref: http://web.mit.edu/16.90/BackUp/www/pdfs/Chapter14.pdf
    # duration of groundmotion [s]
    duration = gm_dt * gm_nsteps

    # trial analysis step, I do not know how this came about
    ktrial = lx / np.sqrt(rock_vs)
    if gm_dt < ktrial:
        nsteps = gm_nsteps
        dt = gm_dt
    else:
        nsteps = int(floor(duration / ktrial) + 1)
        dt = duration / nsteps

    print(f"Number of steps: {nsteps}, time step, dt: {dt}")
    print(f"Newmark: gamma: {newmark_gamma}, beta: {newmark_beta}")
    print(f"Rayleigh: alpham: {alpham}, betak: {betak}, betakinit: {betakinit}, betakcomm: {betakcomm}")
    ops.constraints("Transformation")
    ops.test("NormDispIncr", 1.0e-3, 15, 5)
    ops.algorithm("Newton")
    ops.numberer("RCM")
    ops.system("ProfileSPD")
    ops.integrator("Newmark", newmark_gamma, newmark_beta)
    ops.rayleigh(alpham, betak, betakinit, betakcomm)
    ops.analysis("Transient")
    ops.analyze(nsteps, dt)

    print(f"Dynamic analysis over")

    ops.wipe()


if __name__ == "__main__":
    main()
	#!/usr/bin/env python
	# --coding:utf-8 --

	import numpy as np
	import openseespy.opensees as ops
	from pathlib import Path

	# Example problem:
	# Description: 1D Total Stress Site Response Analysis, Single soil layer
	# Original authors of this example: Christopher McGann and Pedro Arduino, University of Washington
	# ref: https://opensees.berkeley.edu/wiki/index.php/Site_Response_Analysis_of_a_Layered_Soil_Column_(Total_Stress_Analysis)
	# Units: kN, m, sec


	def roundup(num):
	"""Check if a float is a whole number if not round up"""
	return int(num) if num.is_integer() else int(np.ceil(num))

	def data_recorders(fpath, fname_prefix, number_of_nodes, number_of_elements, motiondt):
	# this is a deviation from the original solution, since only ground
	# accelerations, velocity and displacements are considered, for now only ground
	# node values will be recorded

	# fmt: off
	# recorders for nodal displacements, velocity and acceleration at each time step
	restypes = ["disp", "vel", "accel"]
	for restype in restypes:
	fname = fname_prefix + restype + ".out"
	print(f"Filename: {fname}")
	fout = str(fpath.joinpath(fname))
	ops.recorder(
	"Node",
	"-file", fout,
	"-time",
	"-dT", motiondt,
	"-nodeRange", 1, number_of_nodes,
	"-dof", *[1, 2],
	restype,
	)

	# # ref: http://opensees.berkeley.edu/wiki/index.php/PressureIndependMultiYield_Material
	# # For 2D problems, the stress output follows this order: σxx, σyy, σzz, σxy, ηr, where ηr
	# # is the ratio between the shear (deviatoric) stress and peak shear strength at the current
	# # confinement (0<=ηr<=1.0). The strain output follows this order: εxx, εyy, γxy.
	# stresses and strains at each Gauss point in the soil elements
	gausspoints = ["1", "2", "3", "4"]
	for gpt in gausspoints:
	fname = fname_prefix + "_stress" + gpt + ".out"
	fout = str(fpath.joinpath(fname))
	ops.recorder(
	"Element",
	"-file", fout,
	"-time",
	"-dT", motiondt,
	"-eleRange", 1, number_of_elements,
	"material", gpt,
	"stress",
	)
	fname = fname_prefix + "_strain" + gpt + ".out"
	fout = str(fpath.joinpath(fname))
	ops.recorder(
	"Element",
	"-file", fout,
	"-time",
	"-dT", motiondt,
	"-eleRange", 1, number_of_elements,
	"material", gpt,
	"strain",
	)
	# fmt: on

	def main():
	# file path and filenames
	path_root = Path(r"./")
	path_data = path_root / "data"

	# general filename for recorder outputs
	fname_prefix = "ex_g_02_siteresponse_totalstress_"

	# velocity time history file
	velfile = "velocityhistory.out"
	fvel = str(path_data.joinpath(velfile))
	print(f"Ground motion file path: {fvel}")
	fvel = "data/velocityHistory.out"

	# 1. Define analysis parameters

	# soil profile, thickness of the soil layer [m]
	depth_soillayer = 40.0

	# acceleration due to gravity [m/s2]
	accg = 9.81

	# material properties
	# # soil saturated density [Mg/m3], or soil mass density
	soil_satdensity = 1.7
	# soil shear wave velocity [m/s]
	soil_vs = 250.0
	# Poisson's ratio of the soil
	nu = 0.0
	# the undrained shear strength, c of Mohr-Coulomb failure envelop [kPa]
	cohesion = 95.0
	# soil friction angle, phi of Mohr-Coulomb failure envelop [deg]
	friction_angle = 0.0
	# peak shear strain
	peak_shear_strain = 0.05
	# reference pressure [kPa]
	reference_pressure = 80.0
	# pressure dependency coefficient
	pressure_dependence = 0.0
	# bedrock shear wave velocity [m/s]
	rock_vs = 760.0
	# bedrock mass density [Mg/m^3]
	rock_density = 2.4
	# soil shear modulus [kPa]
	modG = soil_satdensity * soil_vs ** 2
	# soil elastic modulus [kPa]
	modE = 2 * modG * (1 + nu)
	# soil bulk modulus [kPa]
	modK = modE / (3 * (1 - 2 * nu))
	# number of yield surfaces
	nsurf = 25

	print(f"Soil elastic properites:")
	print(f"G: {modG/1000:6.1f} MPa, E: {modE/1000:6.1f} MPa, K: {modK/1000:6.1f} MPa")

	# ground motion paramters
	# time step in ground motion record [s]
	gm_dt = 0.005
	# number of steps in ground motion record
	# this depends on the number of points in the ground motion file
	gm_nsteps = 7990

	# wavelength parameters
	# highest frequency desired to be well resolved [Hz]
	fmax = 100.0
	# wavelength of the highest resolved frequency
	lamdamin = soil_vs / fmax
	# number of elements per one wavelength
	nele_wave = 10

	# rayleigh damping parameters (read on this)
	# damping ratio
	damping_ratio = 0.02
	# lower frequency [Hz]
	f_lower = 0.2
	omega_lower = 2 * np.pi * f_lower
	# upper frequency [Hz]
	f_upper = 20
	omega_upper = 2 * np.pi * f_upper
	# rayleigh damping coefficients
	alpham = 2.0 * damping_ratio * omega_lower * omega_upper / (omega_lower + omega_upper)
	betak = 2.0 * damping_ratio / (omega_lower + omega_upper)
	betakinit = 0.0
	betakcomm = 0.0
	print(f"Damping coefficients: alpham = {alpham:6.3f}, betak = {betak:6.5f}")

	# Newmark parameters
	newmark_gamma = 0.5
	newmark_beta = 0.25

	# acceleration due to gravity in the x and y direction
	accg_x = 0.0
	accg_y = -accg
	wgt_x = 0.0
	wgt_y = accg_y * soil_satdensity

	# 2. define mesh geometry
	# number of element in x and y direction
	nx = 1
	# trial for vertical element size
	# initial element size
	h_trial = lamdamin / nele_wave
	# if the number of elements is not an integer round up
	ny = roundup(depth_soillayer / h_trial)
	print(f"Number of elements: nx = {nx}, ny = {ny}")

	# size of elements in x and y direction of the nodes
	ly = float(depth_soillayer / ny)
	lx = ly
	print(f"Element size: lx = {lx:6.3f} m, ly = {ly:6.3f} m")

	# soil element thickness
	thickness = 1.0
	# analysis type, "PlaneStrain" or "PlaneStress"
	analysis_type = "PlaneStrain"
	# surface pressure on elements
	surface_pressure = 0.0

	# start node number
	start_node_number = 1

	# number of nodes in the vertical direction
	num_node_y = 2 * (ny + 1)

	# wipe previous analyses
	ops.wipe()

	# corner nodes are 2D with 2 DOF (displacements in x1 and x2)
	ops.model("basic", "-ndm", 2, "-ndf", 2)

	ycoord = 0.0
	count = 0

	# loop over nodes
	a = range(1, num_node_y + 1, 2)
	for j in range(1, num_node_y, 2):
	ycrd = ycoord + count * ly
	ops.node(j, 0.0, ycrd)
	ops.node(j+1, lx, ycrd)
	print(f"Node: {j}, x: 0.0, y: {ycrd}")
	print(f"Node: {j + 1}, x: {lx}, y: {ycrd}")
	count += 1

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

	# 4. Define dashpot nodes
	ops.node(2000, 0.0, 0.0)
	ops.node(2001, 0.0, 0.0)
	print("Finished creating dashpot nodes....")

	# 5. Define boundary conditions and equalDOF
	# definie fixity of the base nodes
	ops.fix(1, 0, 1)
	ops.fix(2, 0, 1)

	# define fixity of the dashpot nodes
	ops.fix(2000, 1, 1)
	ops.fix(2001, 0, 1)

	# define equalDOF for simple shear deformation of soil elements
	for k in range(3, num_node_y, 2):
	ops.equalDOF(k, k + 1, 1, 2)
	print(f"equalDOF: rnode: {k}, cnode: {k + 1}, dofs: 1 2")

	# define equalDOF for dashpot and base soil nodes
	ops.equalDOF(1, 2, 1)
	ops.equalDOF(1, 2001, 1)

	print(f"Created all boundary conditions and equalDOF:")

	# 6. Define soil materials
	soil_mattag = 1
	# nd = 2 for plane strain and 3 for 3-D
	nd = 2
	# number of dimensions (nd) 2 for plain-strain 3 for 3D
	ops.nDMaterial(
	"PressureIndependMultiYield",
	soil_mattag,
	nd,
	soil_satdensity,
	modG,
	modK,
	cohesion,
	peak_shear_strain,
	friction_angle,
	reference_pressure,
	pressure_dependence,
	nsurf,
	)
	print(f"Soil material created")

	element_mass_density = 0.0
	print(
	"{0:5s}, {1:6s}, {2:6s}, {3:6s}, {4:6s}".format(
	"Ele", "Node i", "Node j", "Node k", "Node l"
	))

	for j in range(1, ny + 1):
	nodei = 2 * j - 1
	nodej = nodei + 1
	nodek = nodei + 3
	nodel = nodei + 2
	ops.element("quad", j, nodei, nodej, nodek, nodel, 1.0, "PlaneStrain", 1.0, 0.0, 0.0, wgt_x, wgt_y)
	print("{0:5d}, {1:6d}, {2:6d}, {3:6d}, {4:6d}".format(j, nodei, nodej, nodek, nodel))

	print(f"All soil elements created")

	# material and elements for viscous damping
	# dashpot coefficients, again figure out how this coefficient comes about
	dashpot_c = lx * rock_density * rock_vs
	# for linear damping, power factor = 1
	dashpot_alpha = 1.0
	# dashpot material
	dashpot_mattag = 4000
	ops.uniaxialMaterial("Viscous", 4000, dashpot_c, 1.0)
	# element, connecting two dashpot nodes
	dashpot_elenum = 5000

	ops.element(
	"zeroLength",
	5000,
	*[2000, 2001],
	"-mat",
	4000,
	"-dir",
	*[1],
	)

	print("Finished creating dashpot material and element...")
	print(f"Gravity analysis:")
	# update material stage to 0 for elastic behavior
	# ops.updateMaterialStage("-material", soil_mattag, "-stage", 0)
	print(f"Material updated for elastic behavior during gravity loading")

	# gravity loading
	ops.constraints("Transformation")
	ops.test("NormDispIncr", 1e-5, 30, 1)
	ops.algorithm("Newton")
	ops.numberer("RCM")
	ops.system("ProfileSPD")
	ops.integrator("Newmark", newmark_gamma, newmark_beta)
	ops.analysis("Transient")
	ops.analyze(10, 500.0)
	print(f"Elastic gravity analysis over")

	# update material stage to include plastic analysis
	ops.updateMaterialStage("-material", soil_mattag, "-stage", 1)
	ops.analyze(40, 500.0)
	print(f"Plastic gravity analysis over")
	ops.setTime(0.0)
	ops.wipeAnalysis()


	print(f"Create data recorders")
	# specify the node and element at which the acceleration, velocity, displacement,
	# stress and strain time histories are required, for surface node
	last_node_number = 322
	last_element_number = 160
	data_recorders(
	path_data, fname_prefix, last_node_number, last_element_number, gm_dt
	)
	print(f"Finished creating data recorders")

	# destory previous analysis object and start consolidation
	print(
	f"Wiping gravity analysis, resetting time to 0.0, and starting dynamic analysis:"
	)


	# define constant factor for applied velocity, find out how this comes about?
	# this is what I found, the link might not work the next time
	# http://www.bgu.ac.il/geol/hatzor/articles/Bao_Hatzor_Huang_2012.pdf
	cfactor = dashpot_c

	# load pattern and timeseries for applied force history
	print(f"Ground motion file: {fvel}, dt: {gm_dt}, cfactor: {cfactor}")

	# load timeseries tag
	tstag_eq = 20
	# load pattern tag
	pattag_eq = 10
	# node at which the excitation is applied
	nodetag_eq = 1
	loadvals = [1.0, 0.0]
	ops.timeSeries(
	"Path", tstag_eq, "-dt", gm_dt, "-filePath", fvel, "-factor", cfactor
	)

	ops.pattern("Plain", pattag_eq, tstag_eq)
	ops.load(nodetag_eq, *loadvals)

	# determine the analysis time step using CFL condition
	# CFL stands for Courant-Friedrichs-Levy
	# ref: http://web.mit.edu/16.90/BackUp/www/pdfs/Chapter14.pdf
	# duration of groundmotion [s]
	duration = gm_dt * gm_nsteps

	# trial analysis step, I do not know how this came about
	ktrial = lx / np.sqrt(rock_vs)
	if gm_dt < ktrial:
	nsteps = gm_nsteps
	dt = gm_dt
	else:
	nsteps = int(floor(duration / ktrial) + 1)
	dt = duration / nsteps

	print(f"Number of steps: {nsteps}, time step, dt: {dt}")
	print(f"Newmark: gamma: {newmark_gamma}, beta: {newmark_beta}")
	print(f"Rayleigh: alpham: {alpham}, betak: {betak}, betakinit: {betakinit}, betakcomm: {betakcomm}")
	ops.constraints("Transformation")
	ops.test("NormDispIncr", 1.0e-3, 15, 5)
	ops.algorithm("Newton")
	ops.numberer("RCM")
	ops.system("ProfileSPD")
	ops.integrator("Newmark", newmark_gamma, newmark_beta)
	ops.rayleigh(alpham, betak, betakinit, betakcomm)
	ops.analysis("Transient")
	ops.analyze(nsteps, dt)

	print(f"Dynamic analysis over")

	ops.wipe()


	if __name__ == "__main__":
	main()