Skip to content

Commit

Permalink
Rearchitect and jit CS fatigue (#2984)
Browse files Browse the repository at this point in the history 
* Rearchitect and jit CS fatigue

Speeds up my convergence by x10 🤯

* black format

* Update cs_fatigue.py
  • Loading branch information
@je-cook
je-cook authored Nov 3, 2023
1 parent ae12fc9 commit 7e5e043
Show file tree
Hide file tree
Showing 2 changed files with 113 additions and 108 deletions.
198 changes: 107 additions & 91 deletions process/cs_fatigue.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
from numba import njit
from process.fortran import constants
from process.fortran import cs_fatigue_variables as csfv
import numpy
Expand Down Expand Up @@ -47,45 +48,42 @@ def ncycle(
delta = 1.0e-4

# Initialise number of cycles
n_cycle = 0.0
N_pulse = 0.0
Kmax = 0.0

# factor 2 taken as saftey factors in the crack sizes
# Default CS steel undergoes fast fracture when SIF > 200 MPa, under a saftey factor 1.5 we use 133MPa
pi_2_arr = numpy.array([numpy.pi / 2.0e0, 0])
while (
(a <= t_structural_vertical / csfv.sf_vertical_crack)
and (c <= t_structural_radial / csfv.sf_radial_crack)
and (Kmax <= csfv.fracture_toughness / csfv.sf_fast_fracture)
):
# find SIF max from SIF_a and SIF_c
Ka = self.surface_stress_intensity_factor(
Ka, Kc = self.surface_stress_intensity_factor(
hoop_stress_MPa,
t_structural_vertical,
t_structural_radial,
a,
c,
numpy.pi / 2.0e0,
)
Kc = self.surface_stress_intensity_factor(
hoop_stress_MPa, t_structural_vertical, t_structural_radial, a, c, 0.0e0
pi_2_arr,
)
Kmax = max(Ka, Kc)

# run euler_method and find number of cycles needed to give crack increase
deltaN = delta / (CR * (Kmax**csfv.paris_power_law))

# update a and c, N
# update a and c, N (+= doesnt work for fortran (?) reasons)
a = a + delta * (Ka / Kmax) ** csfv.paris_power_law
c = c + delta * (Kc / Kmax) ** csfv.paris_power_law
N_pulse = N_pulse + deltaN

# two pulses - ramp to Vsmax and ramp down per cycle
n_cycle = N_pulse / 2.0e0

return n_cycle, t_crack_radial
return N_pulse / 2.0e0, t_crack_radial

def embedded_stress_intensity_factor(self, hoop_stress, t, w, a, c, phi):
@staticmethod
@njit(cache=True)
def embedded_stress_intensity_factor(hoop_stress, t, w, a, c, phi):
# ! Assumes an embedded elliptical efect geometry
# ! geometric quantities
# ! hoop_stress - change in hoop stress over cycle
Expand All @@ -99,47 +97,53 @@ def embedded_stress_intensity_factor(self, hoop_stress, t, w, a, c, phi):
# ! Ref: C. Jong, Magnet Structural Design
# ! Criteria Part 1: Main Structural Components and Welds 2012

# reuse of calc
a_c = a / c
a_t = a / t
cos_phi = numpy.cos(phi)
cos_phi_2 = cos_phi**2.0e0
sin_phi_2 = numpy.sin(phi) ** 2.0e0

if a <= c:
Q = 1.0e0 + 1.464e0 * (a / c) ** 1.65e0
Q = 1.0e0 + 1.464e0 * a_c**1.65e0
m1 = 1.0e0
m2 = 0.05e0 / (0.11e0 + (a / c) ** (3.0e0 / 2.0e0))
m3 = 0.29e0 / (0.23e0 + (a / c) ** (3.0e0 / 2.0e0))
g = 1.0e0 - (
(a / t) ** 4.0e0
* numpy.sqrt(2.6e0 - (2.0e0 * a / t))
/ (1.0e0 + 4.0e0 * (a / c))
) * abs(numpy.cos(phi))
f_phi = (
(a / c) ** 2.0e0 * (numpy.cos(phi)) ** 2.0e0 + (numpy.sin(phi) ** 2.0)
) ** (1.0e0 / 4.0e0)
f_w = numpy.sqrt(
1.0e0 / numpy.cos(numpy.sqrt(a / t) * numpy.pi * c / (2.0e0 * w))
)
elif a > c:
Q = 1.0e0 + 1.464e0 * (c / a) ** 1.65e0
m1 = numpy.sqrt(c / a)
m2 = 0.05e0 / (0.11e0 + (a / c) ** (3.0e0 / 2.0e0))
m3 = 0.29e0 / (0.23e0 + (a / c) ** (3.0e0 / 2.0e0))
g = 1.0e0 - (
(a / t) ** 4.0e0
* numpy.sqrt(2.6e0 - (2.0e0 * a / t))
/ (1.0e0 + 4.0e0 * (a / c))
) * abs(numpy.cos(phi))
f_phi = (
(c / a) ** 2.0e0 * (numpy.sin(phi)) ** 2.0e0 + (numpy.cos(phi) ** 2.0e0)
) ** (1.0e0 / 4.0e0)
f_w = numpy.sqrt(
1.0e0 / numpy.cos(numpy.sqrt(a / t) * numpy.pi * c / (2.0e0 * w))
)
f_phi = (a_c**2.0e0 * cos_phi_2 + sin_phi_2) ** 0.25e0
else: # elif a > c:
c_a = c / a
Q = 1.0e0 + 1.464e0 * c_a**1.65e0
m1 = numpy.sqrt(c_a)
f_phi = (c_a**2.0e0 * sin_phi_2 + cos_phi_2) ** 0.25e0

# compute the unitless geometric correction
f = (m1 + m2 * (a / t) ** 2.0e0 + m3 * (a / t) ** 4.0e0) * g * f_phi * f_w

# compute the stress intensity factor
k = hoop_stress * f * numpy.sqrt(numpy.pi * a / Q)
return k
return (
hoop_stress
* ( # f
(
m1
+ (0.05e0 / (0.11e0 + a_c**1.5e0)) * a_t**2.0e0 # m2 *a_t^2
+ (0.29e0 / (0.23e0 + a_c**1.5e0)) * a_t**4.0e0 # m3 *a_t^4
)
* ( # g
1.0e0
- (
a_t**4.0e0
* numpy.sqrt(2.6e0 - (2.0e0 * a_t))
/ (1.0e0 + 4.0e0 * a_c)
)
* abs(cos_phi)
)
* f_phi
* numpy.sqrt( # f_w
1.0e0 / numpy.cos(numpy.sqrt(a_t) * numpy.pi * c / (2.0e0 * w))
)
)
* numpy.sqrt(numpy.pi * a / Q)
)

def surface_stress_intensity_factor(self, hoop_stress, t, w, a, c, phi):
@staticmethod
@njit(cache=True)
def surface_stress_intensity_factor(hoop_stress, t, w, a, c, phi):
# ! Assumes an surface semi elliptical defect geometry
# ! geometric quantities
# ! hoop_stress - change in hoop stress over cycle
Expand All @@ -155,56 +159,68 @@ def surface_stress_intensity_factor(self, hoop_stress, t, w, a, c, phi):

bending_stress = 0.0e0 # * 3.0 * M / (w*d**2.0)

# reuse of calc
a_t = a / t
a_t_2 = a_t**2.0e0
sin_phi = numpy.sin(phi)
cos_phi_2 = numpy.cos(phi) ** 2.0e0

if a <= c:
Q = 1.0e0 + 1.464e0 * (a / c) ** 1.65e0
m1 = 1.13e0 - 0.09e0 * a / c
m2 = -0.54e0 + 0.89e0 / (0.2e0 + a / c)
m3 = 0.5e0 - 1.0e0 / (0.65e0 + a / c) + 14.0e0 * (1 - a / c) ** 24.0e0
g = (
1.0e0
+ (0.1e0 + 0.35e0 * (a / c) ** 2.0e0)
* (1.0e0 - numpy.sin(phi)) ** 2.0e0
# reuse of calc
a_c = a / c

Q = 1.0e0 + 1.464e0 * a_c**1.65e0
m1 = 1.13e0 - 0.09e0 * a_c
m2 = -0.54e0 + 0.89e0 / (0.2e0 + a_c)
m3 = 0.5e0 - 1.0e0 / (0.65e0 + a_c) + 14.0e0 * (1 - a_c) ** 24.0e0
g = 1.0e0 + (0.1e0 + 0.35e0 * a_c**2.0e0) * (1.0e0 - sin_phi) ** 2.0e0
f_phi = (a_c**2.0e0 * cos_phi_2 + sin_phi**2.0e0) ** 0.25e0
p = 0.2e0 + a_c + 0.6e0 * a_t
H1 = 1.0e0 - 0.34e0 * a_t - 0.11e0 * a * a / (c * t)
H2 = (
1.0
+ (-1.22e0 - 0.12e0 * a_c) * a_t # G21 * a / t
+ (0.55e0 - 1.05e0 * a_c**0.75e0 + 0.47e0 * a_c**1.5e0) # G22
* a_t_2
)
f_phi = (
(a / c) ** 2.0e0 * (numpy.cos(phi)) ** 2.0e0 + (numpy.sin(phi)) ** 2.0e0
) ** (1.0e0 / 4.0e0)
f_w = numpy.sqrt(
1.0e0 / numpy.cos(numpy.sqrt(a / t) * numpy.pi * c / (2.0e0 * w))
)
p = 0.2e0 + a / c + 0.6e0 * a / t
G21 = -1.22e0 - 0.12e0 * a / c
G22 = 0.55e0 - 1.05e0 * (a / c) ** 0.75e0 + 0.47e0 * (a / c) ** 1.5e0
H1 = 1.0e0 - 0.34e0 * a / t - 0.11e0 * a * a / (c * t)
H2 = 1.0e0 + G21 * a / t + G22 * (a / t) ** 2.0e0
elif a > c:
Q = 1.0e0 + 1.464e0 * (c / a) ** 1.65e0
m1 = numpy.sqrt(c / a) * (1.0e0 + 0.04e0 * c / a)
m2 = 0.2e0 * (c / a) ** 4.0e0
m3 = -0.11e0 * (c / a) ** 4.0e0
g = (
else: # elif a > c:
c_a = c / a
c_a_4 = c_a**4.0e0

Q = 1.0e0 + 1.464e0 * c_a**1.65e0
m1 = numpy.sqrt(c_a) * (1.0e0 + 0.04e0 * c_a)

m2 = 0.2e0 * c_a_4
m3 = -0.11e0 * c_a_4

g = 1.0e0 + (0.1e0 + 0.35e0 * c_a * a_t_2) * (1.0e0 - sin_phi) ** 2.0e0
f_phi = (c_a**2.0e0 * sin_phi**2.0e0 + cos_phi_2) ** 0.25e0
p = 0.2e0 + c_a + 0.6e0 * a_t
H1 = (
1.0e0
+ (0.1e0 + 0.35e0 * (c / a) * (a / t) ** 2.0e0)
* (1.0e0 - numpy.sin(phi)) ** 2.0e0
+ (-0.04e0 - 0.41e0 * c_a) * a_t # G11 * a / t
+ (0.55e0 - 1.93e0 * c_a**0.75e0 + 1.38e0 * c_a**1.5e0) # G12
* a_t_2
)
f_phi = (
(c / a) ** 2.0e0 * (numpy.sin(phi)) ** 2.0e0 + (numpy.cos(phi)) ** 2.0e0
) ** (1.0e0 / 4.0e0)
f_w = numpy.sqrt(
1.0e0 / numpy.cos(numpy.sqrt(a / t) * numpy.pi * c / (2.0e0 * w))
H2 = (
1.0e0
+ (-2.11e0 + 0.77e0 * c_a) * a_t # G21 * a / t
+ (0.55e0 - 0.72e0 * c_a * 0.75e0 + 0.14e0 * c_a * 1.5e0) * a_t_2 # G22
)
p = 0.2e0 + c / a + 0.6e0 * a / t
G11 = -0.04e0 - 0.41e0 * c / a
G12 = 0.55e0 - 1.93e0 * (c / a) ** 0.75e0 + 1.38e0 * (c / a) ** 1.5e0
G21 = -2.11e0 + 0.77e0 * c / a
G22 = 0.55e0 - 0.72e0 * (c / a) * 0.75e0 + 0.14e0 * (c / a) * 1.5e0
H1 = 1.0e0 + G11 * a / t + G12 * (a / t) ** 2.0e0
H2 = 1.0e0 + G21 * a / t + G22 * (a / t) ** 2.0e0

# compute the unitless geometric correction
Hs = H1 + (H2 - H1) * (numpy.sin(phi)) ** p
f = (m1 + m2 * (a / t) ** 2.0e0 + m3 * (a / t) ** 4.0e0) * g * f_phi * f_w

# compute the stress intensity factor
k = (hoop_stress + Hs * bending_stress) * f * numpy.sqrt(numpy.pi * a / Q)

return k
return (
( # hoop_stress + Hs * bending_stress
hoop_stress + (H1 + (H2 - H1) * sin_phi**p) * bending_stress
)
* ( # f
(m1 + m2 * a_t_2 + m3 * a_t**4.0e0)
* g
* f_phi
* numpy.sqrt( # f_w
1.0e0 / numpy.cos(numpy.sqrt(a_t) * numpy.pi * c / (2.0e0 * w))
)
)
* numpy.sqrt(numpy.pi * a / Q)
)
23 changes: 6 additions & 17 deletions process/pfcoil.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,24 +72,13 @@ def pfcoil(self):
pcls0 = np.zeros(pfv.ngrpmx, dtype=int)
ncls0 = np.zeros(pfv.ngrpmx + 2, dtype=int)

pf.rcls0 = np.zeros((pfv.ngrpmx, pfv.nclsmx), order="F")
pf.zcls0 = np.zeros((pfv.ngrpmx, pfv.nclsmx), order="F")
pf.rcls0, pf.zcls0 = np.zeros((2, pfv.ngrpmx, pfv.nclsmx), order="F")
pf.ccls0 = np.zeros(int(pfv.ngrpmx / 2))
sigma = np.zeros(pfv.ngrpmx)
work2 = np.zeros(pfv.ngrpmx)
rc = np.zeros(pfv.nclsmx)
zc = np.zeros(pfv.nclsmx)
cc = np.zeros(pfv.nclsmx)
xc = np.zeros(pfv.nclsmx)
brin = np.zeros(pfv.nptsmx)
bzin = np.zeros(pfv.nptsmx)
rpts = np.zeros(pfv.nptsmx)
zpts = np.zeros(pfv.nptsmx)
bfix = np.zeros(lrow1)
bvec = np.zeros(lrow1)
gmat = np.zeros((lrow1, lcol1), order="F")
umat = np.zeros((lrow1, lcol1), order="F")
vmat = np.zeros((lrow1, lcol1), order="F")
sigma, work2 = np.zeros((2, pfv.ngrpmx))
rc, zc, cc, xc = np.zeros((4, pfv.nclsmx))
brin, bzin, rpts, zpts = np.zeros((4, pfv.nptsmx))
bfix, bvec = np.zeros((2, lrow1))
gmat, umat, vmat = np.zeros((3, lrow1, lcol1), order="F")
signn = np.zeros(2)
aturn = np.zeros(pfv.ngc2)

Expand Down

0 comments on commit 7e5e043

Please sign in to comment.