Skip to content

Commit

Permalink
Fix Simpson backward compatbility
Browse files Browse the repository at this point in the history
  • Loading branch information
@AnkitBarik
AnkitBarik committed Jul 30, 2024
1 parent d06d699 commit 2a16ba9
Show file tree
Hide file tree
Showing 4 changed files with 38 additions and 38 deletions.
4 changes: 2 additions & 2 deletions python/magic/bLayers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -478,7 +478,7 @@ def __init__(self, iplot=False, quiet=False):
y = RolC[self.radius >= self.ro-self.bcTopSlope]
x = self.radius[self.radius >= self.ro-self.bcTopSlope]
try:
self.rolTop = simps(3.*y*x**2, x)/(self.ro**3-(self.ro-self.bcTopSlope)**3)
self.rolTop = simps(3.*y*x**2, x=x)/(self.ro**3-(self.ro-self.bcTopSlope)**3)
except IndexError:
self.rolTop = 0.

Expand All @@ -500,7 +500,7 @@ def __init__(self, iplot=False, quiet=False):

y = RolC[self.radius <= self.ri+self.bcBotSlope]
x = self.radius[self.radius <= self.ri+self.bcBotSlope]
self.rolBot = simps(3.*y*x**2, x)/((self.ri+self.bcBotSlope)**3-self.ri**3)
self.rolBot = simps(3.*y*x**2, x=x)/((self.ri+self.bcBotSlope)**3-self.ri**3)
print('reynols bc, reynolds bulk', self.rebl, self.rebulk)
print('reh bc, reh bulk', self.rehbl, self.rehbulk)
print('rolbc, rolbulk, roltop, rolbot', self.rolbl, self.rolbulk,
Expand Down
2 changes: 1 addition & 1 deletion python/magic/butterfly.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,7 +284,7 @@ def fourier2D(self, renorm=False):
self.omega = w[1:nt//2+1]
self.amp1D = np.zeros_like(self.omega)
for i in range(len(self.omega)):
self.amp1D[i] = simps(self.amp[:, i], self.theta)
self.amp1D[i] = simps(self.amp[:, i], x=self.theta)

fig = plt.figure()
ax = fig.add_subplot(211)
Expand Down
4 changes: 2 additions & 2 deletions python/magic/libmagic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1201,12 +1201,12 @@ def horizontal_mean(field, colat, std=False):
"""

field_m = field.mean(axis=0) # Azimuthal average
field_mean = 0.5 * simps(field_m.T*np.sin(colat), colat, axis=-1)
field_mean = 0.5 * simps(field_m.T*np.sin(colat), x=colat, axis=-1)

if std:
dat = (field-field_mean)**2
dat_m = dat.mean(axis=0)
dat_mean = 0.5 * simps(dat_m.T*np.sin(colat), colat, axis=-1)
dat_mean = 0.5 * simps(dat_m.T*np.sin(colat), x=colat, axis=-1)
field_std = np.sqrt(dat_mean)
return field_mean, field_std
else:
Expand Down
66 changes: 33 additions & 33 deletions python/magic/thHeat.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,8 +149,8 @@ def __init__(self, iplot=False, angle=10, pickleName='thHeat.pickle',
for i in range(self.ntheta):
th2D[i, :] = self.colat[i]

self.temprmmean = 0.5*simps(self.tempmean*np.sin(th2D), th2D, axis=0)
self.temprmstd = 0.5*simps(self.tempstd*np.sin(th2D), th2D, axis=0)
self.temprmmean = 0.5*simps(self.tempmean*np.sin(th2D), x=th2D, axis=0)
self.temprmstd = 0.5*simps(self.tempstd*np.sin(th2D), x=th2D, axis=0)
sinTh = np.sin(self.colat)

# Conducting temperature profile (Boussinesq only!)
Expand All @@ -177,46 +177,46 @@ def __init__(self, iplot=False, angle=10, pickleName='thHeat.pickle',
# Close to the equator
mask2D = (th2D>=np.pi/2.-angle/2.)*(th2D<=np.pi/2+angle/2.)
mask = (self.colat>=np.pi/2.-angle/2.)*(self.colat<=np.pi/2+angle/2.)
fac = 1./simps(sinTh[mask], self.colat[mask])
self.nussBotEq = fac*simps(self.nussbotmean[mask]*sinTh[mask], self.colat[mask])
self.nussTopEq = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
fac = 1./simps(sinTh[mask], x=self.colat[mask])
self.nussBotEq = fac*simps(self.nussbotmean[mask]*sinTh[mask], x=self.colat[mask])
self.nussTopEq = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])
sinC = sinTh.copy()
sinC[~mask] = 0.
fac = 1./simps(sinC, self.colat)
fac = 1./simps(sinC, x=self.colat)
tempC = self.tempmean.copy()
tempC[~mask2D] = 0.
self.tempEqmean = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
self.tempEqmean = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)
tempC = self.tempstd.copy()
tempC[~mask2D] = 0.
self.tempEqstd = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
self.tempEqstd = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)

dtempEq = rderavg(self.tempEqmean, self.radius)
self.betaEq = dtempEq[len(dtempEq)//2]

# 45\deg inclination
mask2D = (th2D>=np.pi/4.-angle/2.)*(th2D<=np.pi/4+angle/2.)
mask = (self.colat>=np.pi/4.-angle/2.)*(self.colat<=np.pi/4+angle/2.)
fac = 1./simps(np.sin(self.colat[mask]), self.colat[mask])
nussBot45NH = fac*simps(self.nussbotmean[mask]*sinTh[mask], self.colat[mask])
nussTop45NH = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
fac = 1./simps(np.sin(self.colat[mask]), x=self.colat[mask])
nussBot45NH = fac*simps(self.nussbotmean[mask]*sinTh[mask], x=self.colat[mask])
nussTop45NH = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])
sinC = sinTh.copy()
sinC[~mask] = 0.
fac = 1./simps(sinC, self.colat)
fac = 1./simps(sinC, x=self.colat)
tempC = self.tempmean.copy()
tempC[~mask2D] = 0.
temp45NH = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
temp45NH = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)

mask2D = (th2D>=3.*np.pi/4.-angle/2.)*(th2D<=3.*np.pi/4+angle/2.)
mask = (self.colat>=3.*np.pi/4.-angle/2.)*(self.colat<=3.*np.pi/4+angle/2.)
fac = 1./simps(np.sin(self.colat[mask]), self.colat[mask])
nussBot45SH = fac*simps(self.nussbotmean[mask]*sinTh[mask], self.colat[mask])
nussTop45SH = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
fac = 1./simps(np.sin(self.colat[mask]), x=self.colat[mask])
nussBot45SH = fac*simps(self.nussbotmean[mask]*sinTh[mask], x=self.colat[mask])
nussTop45SH = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])
sinC = sinTh.copy()
sinC[~mask] = 0.
fac = 1./simps(sinC, self.colat)
fac = 1./simps(sinC, x=self.colat)
tempC = self.tempmean.copy()
tempC[~mask2D] = 0.
temp45SH = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
temp45SH = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)

self.nussTop45 = 0.5*(nussTop45NH+nussTop45SH)
self.nussBot45 = 0.5*(nussBot45NH+nussBot45SH)
Expand All @@ -228,33 +228,33 @@ def __init__(self, iplot=False, angle=10, pickleName='thHeat.pickle',
# Polar regions
mask2D = (th2D<=angle/2.)
mask = (self.colat<=angle/2.)
fac = 1./simps(np.sin(self.colat[mask]), self.colat[mask])
nussBotPoNH = fac*simps(self.nussbotmean[mask]*sinTh[mask], self.colat[mask])
nussTopPoNH = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
fac = 1./simps(np.sin(self.colat[mask]), x=self.colat[mask])
nussBotPoNH = fac*simps(self.nussbotmean[mask]*sinTh[mask], x=self.colat[mask])
nussTopPoNH = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])
sinC = sinTh.copy()
sinC[~mask] = 0.
fac = 1./simps(sinC, self.colat)
fac = 1./simps(sinC, x=self.colat)
tempC = self.tempmean.copy()
tempC[~mask2D] = 0.
tempPolNHmean = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
tempPolNHmean = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)
tempC = self.tempstd.copy()
tempC[~mask2D] = 0.
tempPolNHstd = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
tempPolNHstd = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)

mask2D = (th2D>=np.pi-angle/2.)
mask = (self.colat>=np.pi-angle/2.)
fac = 1./simps(np.sin(self.colat[mask]), self.colat[mask])
nussBotPoSH = fac*simps(self.nussbotmean[mask]*sinTh[mask], self.colat[mask])
nussTopPoSH = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
nussBotPoSH = fac*simps(self.nussbotmean[mask]*sinTh[mask], x=self.colat[mask])
nussTopPoSH = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])
sinC = sinTh.copy()
sinC[~mask] = 0.
fac = 1./simps(sinC, self.colat)
fac = 1./simps(sinC, x=self.colat)
tempC = self.tempmean.copy()
tempC[~mask2D] = 0.
tempPolSHmean = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
tempPolSHmean = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)
tempC = self.tempstd.copy()
tempC[~mask2D] = 0.
tempPolSHstd = fac*simps(tempC*np.sin(th2D), th2D, axis=0)
tempPolSHstd = fac*simps(tempC*np.sin(th2D), x=th2D, axis=0)

self.nussBotPo = 0.5*(nussBotPoNH+nussBotPoSH)
self.nussTopPo = 0.5*(nussTopPoNH+nussTopPoSH)
Expand All @@ -268,20 +268,20 @@ def __init__(self, iplot=False, angle=10, pickleName='thHeat.pickle',
angleTC = np.arcsin(self.ri/self.ro)
mask2D = (th2D<=angleTC)
mask = (self.colat<=angleTC)
fac = 1./simps(np.sin(self.colat[mask]), self.colat[mask])
nussITC_NH = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
fac = 1./simps(np.sin(self.colat[mask]), x=self.colat[mask])
nussITC_NH = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])

mask2D = (th2D>=np.pi-angleTC)
mask = (self.colat>=np.pi-angleTC)
fac = 1./simps(np.sin(self.colat[mask]), self.colat[mask])
nussITC_SH = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
nussITC_SH = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])

self.nussITC = 0.5*(nussITC_NH+nussITC_SH)

mask2D = (th2D>=angleTC)*(th2D<=np.pi-angleTC)
mask = (self.colat>=angleTC)*(self.colat<=np.pi-angleTC)
fac = 1./simps(sinTh[mask], self.colat[mask])
self.nussOTC = fac*simps(self.nusstopmean[mask]*sinTh[mask], self.colat[mask])
self.nussOTC = fac*simps(self.nusstopmean[mask]*sinTh[mask], x=self.colat[mask])

if iplot:
self.plot()
Expand Down

0 comments on commit 2a16ba9

Please sign in to comment.