Skip to content

calc_GGmax_from_stress_strain() and calc_damping_from_stress_strain() do not check that strain and stress have the same length, and silently broadcast a length-1 array #108

Open

Description

@jsh9

Found while reviewing #74. Pre-existing on main (b764057); #74 keeps the same behavior.

Problem

calc_GGmax_from_stress_strain() checks that strain_in_unit_1 and stress are 1D arrays, but not that they have the same length. It then computes stress / strain_in_unit_1. So:

  • a length-1 stress with a longer strain, or a length-1 strain with a longer stress, is broadcast by numpy, and a G/Gmax curve is returned without an error;
  • any other mismatch fails inside numpy with "operands could not be broadcast together", which does not name the inputs.

hlp.assert_1D_numpy_array(strain_in_unit_1)
hlp.assert_1D_numpy_array(stress)
if strain_in_unit_1[0] == 0:
raise ValueError(
'`strain_in_unit_1` should start with a non-zero value.'
)
if Gmax is None:
Gmax = stress[0] / strain_in_unit_1[0]
G = stress / strain_in_unit_1 # secant modulus
return G / Gmax

calc_damping_from_stress_strain() has no checks of its own. It takes n = len(strain) and passes both arrays to calc_GGmax_from_stress_strain(), so it has the same gaps. With a length-1 stress it returns one damping value per strain. With a length-1 strain it returns [0.]. calc_damping_from_param() computes the stress from the strain itself, so it always passes matching lengths.

strain = strain_in_unit_1
n = len(strain)
G_Gmax = calc_GGmax_from_stress_strain(strain, stress, Gmax=Gmax)
area = np.zeros(n)
damping = np.zeros(n)
area[0] = 0.5 * (strain[0] * G_Gmax[0]) * strain[0]
damping[0] = 2.0 / np.pi * (2.0 * area[0] / G_Gmax[0] / strain[0] ** 2 - 1)
for i in range(1, n):
area[i] = area[i - 1] + 0.5 * (
strain[i - 1] * G_Gmax[i - 1] + strain[i] * G_Gmax[i]
) * (strain[i] - strain[i - 1])
damping[i] = (
2.0 / np.pi * (2 * area[i] / G_Gmax[i] / strain[i] ** 2 - 1)
)
return np.maximum(damping, 0.0) # make sure all damping values are >= 0

import numpy as np
import PySeismoSoil.helper_site_response as sr

strain = np.array([0.1, 0.2, 0.3])
stress = np.array([4.0, 6.0, 7.0])

# Matching lengths (the hand-calculated case in the existing test)
print(sr.calc_GGmax_from_stress_strain(strain, stress))
# [1.         0.75       0.58333333]

# 1. A length-1 stress is broadcast against the 3 strains
print(sr.calc_GGmax_from_stress_strain(strain, stress[:1]))
# [1.         0.5        0.33333333]   (expected: ValueError)

# 2. A length-1 strain is broadcast against the 3 stresses
print(sr.calc_GGmax_from_stress_strain(strain[:1], stress))
# [1.   1.5  1.75]   (expected: ValueError)

# 3. Other mismatches fail inside numpy
try:
    sr.calc_GGmax_from_stress_strain(strain, np.array([4.0, 6.0]))
except Exception as e:
    print(type(e).__name__, e)
# ValueError operands could not be broadcast together with shapes (2,) (3,)
# (expected: a ValueError that names the two inputs)

# Same gaps in calc_damping_from_stress_strain()
print(sr.calc_damping_from_stress_strain(strain, stress[:1], 40.0))
# [0.         0.31830989 0.42441318]   (expected: ValueError)
print(sr.calc_damping_from_stress_strain(strain[:1], stress, 40.0))
# [0.]   (expected: ValueError)
try:
    sr.calc_damping_from_stress_strain(strain, np.array([4.0, 6.0]), 40.0)
except Exception as e:
    print(type(e).__name__, e)
# ValueError operands could not be broadcast together with shapes (2,) (3,)
# (expected: a ValueError that names the two inputs)

The existing tests only use matching lengths: test_calc_GGmax_from_stress_strain_curve__linear_case_GGmax_is_1, test_calc_GGmax_from_stress_strain_curve__a_hand_calculated_case, test_calc_damping_from_stress_strain__case_1 and test_calc_damping_from_stress_strain__case_2 on main (in #74: test_calc_GGmax_from_stress_strain_curve[linear_case_GGmax_is_1], [a_hand_calculated_case], test_calc_damping_from_stress_strain[linear] and [elasto_perfectly_plastic]).

Fix

  1. In calc_GGmax_from_stress_strain(), after the two 1D checks, add:

         hlp.assert_1D_numpy_array(strain_in_unit_1)
         hlp.assert_1D_numpy_array(stress)
    
    +    if len(stress) != len(strain_in_unit_1):
    +        raise ValueError(
    +            '`strain_in_unit_1` and `stress` should have the same length.'
    +        )
    +
         if strain_in_unit_1[0] == 0:

    calc_damping_from_stress_strain() calls calc_GGmax_from_stress_strain() before it uses the arrays, so this one check covers both functions. With only this change on main, all six mismatched calls in the repro raise this ValueError, the matching case still gives [1. 0.75 0.58333333], and all 240 tests pass.

  2. Add a test that both functions raise ValueError for stress of length 1, strain of length 1, and lengths 3 and 2.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions