On Fri, Aug 14, 2020 at 1:32 PM Stephen J. Turnbull <
turnbull.stephen...@u.tsukuba.ac.jp> wrote:
> Ricky Teachey writes:
>
> > I'll try to speak to [what would engineers use a matrix library
> > for] below.
>
> I don't get it. Among other things, the object uwt_per_ft isn't
> defined. I think it's a typo for the one-dimensional array
> γ_per_foot, but not sure. I don't see any advantage to numpy or a
> matrix library in this code, unless you're desperate for speed (see my
> last version where z is a Python list).
>
Sorry, yes, the code I wrote had a couple typos. I'll put the correct
version (tested) at the very bottom for clarity.
> This operation:
> >
> > q_s = γ_prime * z
> >
> > And these operations:
> >
> > P_n_per_foot = μ_s * q_s * C_shaft
> > P_n = np.sum(P_n_per_foot)
>
> Are you saying that the operations above would be clearer in a
> four-function calculator version of a matrix library? I don't see how
> the above is less clear than
>
> q_s = γ_prime * z
> P_n_per_foot = μ_s * q_s * C_shaft
> P_n = math.fsum(P_n_per_foot)
>
> or either is more clear than
>
> z = list(range(1, 21))
> q_s = [γ_prime * x for x in z]
> P_n_per_foot = [μ_s * x * C_shaft for x in q_s]
> P_n = math.fsum(P_n_per_foot)
>
> although they might be faster for large enough values of "20".
>
Well keep in mind that generally an engineer submits calculations-- not
code-- to be reviewed. This leaves me with a few options:
1. Write code the works and then submit a "cleaned up" version of the
code, with parts that don't look like math omitted
2. Write the code, and then write MORE code that prints out the clean,
looks-like-math version
3. Use a notebook (in practice this is not much different than #1)
I often end up doing #1 because it is actually easier and faster.
Hiding the portion of code that looks too much like programming (and not
enough like equations) takes some time, but if I do not do this it will
inevitably lead to a 40 minute phone call from an official from the NY
Department of Transportation, delays, and more hours billed to my client.
One easy way to do this code hiding is to simply copy and paste the *result* of
the calculation and delete the code that created it, kind of like this:
print(f"γ_per_ft = ({', '.join(f'{v:.0f}' for v in γ_per_ft)}) pcf")
Anyway typing the mathematical equations, and not a list comprehension, is
far clearer for the engineers who will be reading my calculations, and who
are not programmers.
Otherwise I end up typing everything twice, just like I end up having to do
using Excel: the code that actually does the calculation, and a "pretty
version" to display to the reviewer what the calculation is doing.
Something like the equation P_n = μ_s * q_s * C_shaft is directly out of
the AASHTO LRFD and NY DoT highway code:
# this looks like gobbledegook to, for example, New York State DoT official
reviewing my calculations
# so i will hide this part of the code, somehow
P_n_per_foot = [μ_s * x * C_shaft for x in q_s]
# this is what the official is expecting to see
P_n_per_foot = μ_s * q_s * C_shaft
Furthermore, the same NYDoT reviewer might come back and say:
*"Rejected: calculation does not take into account change in friction angle
with depth."*
The official would be saying that friction (μ_s) is not a scalar, it should
be another array variable, because the soils apply different friction along
the shaft depth.
So now, because the code expects μ_s to be a scalar, I have to go change
not just the value of μ_s but also the code to turn it into a
comprehension that will work for multiple lists, instead of one list:
μ_s = [ <some list of effective friction coefficients> ]
P_n_per_foot = [μ_s * x * C_shaft for f,x in zip(μ_s, q_s)]
If I could write the code like this, I could just change μ_s to an array,
and everything works fine as-is, AND I would not have to take steps to hide
the code, I can just print it out with no changes, which would be wonderful:
P_n_per_foot = μ_s * q_s * C_shaft
> > That being said, I do think it would be nice to have a "matrix"
> > type in the std lib that supported both element-wise operations
> > like this and basic matrix arithmetic.
>
> If the answer to the question above is "no", then nice for what? If
> "yes", do you think the stdlib 'matrix' module should have a
> 'linspace' class/function? Or that
>
> z = matrix.Matrix(columns=1, data=range(1, 21))
>
> is in general an adequate substitute?
>
The answer is yes. It's a good question on the linspace() function -- I
don't think so. Something like this would do fine:
## depths below grade to 20 ft in 1 foot increments
z = matrix(range(1, 21))
However I would certainly hide this part of the code, and include the copy
and pasted result from this:
print(f"γ_per_ft = ({', '.join(f'{v:.0f}' for v in z)}) pcf")
------------
Sidebar: it would be pretty fantastic if a native matrix library could
display values using f-strings like this:
>>> x = 1.1
>>> print(f"{x=:.0f} pcf")
1 pcf
>>> γ_per_ft = matrix( [120, 122, 118, 117] )
>>> print(f"{γ_per_ft=:.0f} pcf")
(120, 122, 118, 117) pcf
The idea here is that python knows the datatype of γ_per_ft is float, so
it can format all the members of the array *one at a time*.
End of sidebar.
------------
Here is the (now tested) code snippet. I'm sorry it was incorrect before:
import numpy as np
# Nominal skin friction capacity, Pn, of drilled shafts/piers
## depths below grade to 20 ft in 1 foot increments
z = np.linspace(1, 20, num=20)
## create array containing unit weight of soil at each depth
### unit weight log
### for looking up unit weights (pcf) based on depth (ft)
### of form {depth: unit_weight}
log = {3: 120, 9: 118, 12: 122, 15: 106, 18: 102, 21: 110}
### fill in gaps in unit weight log
#### assume first layer extends to ground surface
log[0] = log[next(iter(log))]
#### assume unit weight change occurs after mid point between each layer
log_depths = sorted(list(log.keys()))
log_filled = {}
for depth, next_depth in zip(log_depths[:-1], log_depths[1:]):
# layer boundary
log_filled[depth] = log[depth]
# layer change occurs after midpoint of boundary
log_filled[ (depth+next_depth)/2 ] = log[depth]
#### don't forget bottom of last layer
log_filled[next_depth] = log[next_depth]
#### final per foot unit weights, γ
γ_per_ft = np.array([])
for depth in z:
log_depth = next(d for d in log_filled if d >= depth)
γ_per_ft = np.append(γ_per_ft, log_filled[log_depth])
## calculate per foot stress
### unit weight of water (pcf)
γ_w = 62.4
### effective unit weight given high water table
γ_prime = γ_per_ft - γ_w
q_s = γ_prime * z
## Nominal capacity
### shaft diameter (ft)
D_shaft = 6
### shaft circumference
C_shaft = D_shaft * np.pi
# friction coefficient
friction_angle = 17 # degrees
μ_s = np.tan( friction_angle * np.pi/180 )
# nominal vertical capacity
P_n_per_foot = μ_s * q_s * C_shaft
P_n = np.sum(P_n_per_foot)
print(f"{P_n/2000:0.0f} tonf")
Output:
>>> 30 tonf
---
Ricky.
"I've never met a Kentucky man who wasn't either thinking about going home
or actually going home." - Happy Chandler
_______________________________________________
Python-ideas mailing list -- python-ideas@python.org
To unsubscribe send an email to python-ideas-le...@python.org
https://mail.python.org/mailman3/lists/python-ideas.python.org/
Message archived at
https://mail.python.org/archives/list/python-ideas@python.org/message/WISJ3FW5MQCT3C3SCD2FL4253675TEOY/
Code of Conduct: http://python.org/psf/codeofconduct/
