Skip to content

Commit 39681da

Browse files
SierdSierd
authored
Bart fix (#234)
* updated circular boundaries for sweep solver * Update utils.py * Added 4th generation sweep solver for testing
1 parent 6afe8ee commit 39681da

File tree

4 files changed

+311
-10
lines changed

4 files changed

+311
-10
lines changed

‎aeolis/model.py

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -207,16 +207,21 @@ def initialize(self)-> None:
207207
if 1:
208208
# this is the new quasi 2D stuff
209209
# this is where we make the 1D grid into a 2D grid to ensure process compatibility
210-
self.p['xgrid_file'] = np.transpose(np.stack((self.p['xgrid_file'], self.p['xgrid_file'], self.p['xgrid_file']),axis=1))
210+
211+
# define size of 2D grid 3 is minumum, variable defined for debugging purposes
212+
qnr = 3
213+
214+
self.p['xgrid_file'] = np.transpose(np.stack([self.p['xgrid_file'] for i in range (qnr)],axis=1))
215+
216+
# redefine shape
211217
self.p['ny'], self.p['nx'] = self.p['xgrid_file'].shape
212218

219+
# repeat the above for ygrid_file assume dy is equal to dx
213220
dy = self.p['xgrid_file'][1,2]-self.p['xgrid_file'][1,1]
214-
215-
# repeat the above for ygrid_file
216-
self.p['ygrid_file'] = np.transpose(np.stack((self.p['ygrid_file'], self.p['ygrid_file']+dy, self.p['ygrid_file']+2*dy),axis=1))
221+
self.p['ygrid_file'] = np.transpose(np.stack([self.p['ygrid_file'] + i * dy for i in range(qnr)], axis=1))
217222

218223
# repeat the above for bed_file
219-
self.p['bed_file'] = np.transpose(np.stack((self.p['bed_file'], self.p['bed_file'], self.p['bed_file']),axis=1))
224+
self.p['bed_file'] = np.transpose(np.stack([self.p['bed_file'] for i in range (qnr)],axis=1))
220225

221226
# change from number of points to number of cells
222227
self.p['nx'] -= 1

‎aeolis/netcdf.py

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -297,8 +297,8 @@ def initialize(outputfile, outputvars, s, p, dimensions):
297297
# this is where a quasi 2D model is stored in 1D
298298
if p['ny']+1 == 3:
299299
nc.variables['n'][:] = np.arange(1)
300-
nc.variables['x'][:,:] = s['x'][0,:]
301-
nc.variables['y'][:,:] = s['y'][0,:]
300+
nc.variables['x'][:,:] = s['x'][1,:]
301+
nc.variables['y'][:,:] = s['y'][1,:]
302302
else:
303303
nc.variables['n'][:] = np.arange(p['ny']+1)
304304
nc.variables['x'][:,:] = s['x']

‎aeolis/run_console.py

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,8 +6,11 @@ def main()-> None:
66
'''Runs AeoLiS model in debugging mode.'''
77

88
# configfile = r'c:\Users\weste_bt\aeolis\Tests\RotatingWind\Barchan_Grid270\aeolis.txt'
9-
configfile = r'C:\Users\svries\Documents\GitHub\OE_aeolis-python\aeolis\examples\grainsizevariations\aeolis_horizontalgradient1.txt'
9+
configfile = r'C:\Users\svries\Documents\GitHub\Bart_mass\aeolis_duran.txt'
10+
# configfile = r'C:\Users\svries\Documents\GitHub\Bart_mass\aeolis_windspeed.txt'
11+
1012
aeolis_debug(configfile)
1113

14+
1215
if __name__ == '__main__':
1316
main()

‎aeolis/utils.py

Lines changed: 295 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -823,13 +823,50 @@ def sweep3(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w):
823823
ufn[1:-1,:, :] = 0.5*un[:-1,:, :] + 0.5*un[1:,:, :]
824824

825825
# print(ufs[5,:,0])
826-
827-
#boundary values
826+
827+
# boundary values
828828
ufs[:,0, :] = us[:,0, :]
829829
ufs[:,-1, :] = us[:,-1, :]
830830

831831
ufn[0,:, :] = un[0,:, :]
832832
ufn[-1,:, :] = un[-1,:, :]
833+
# first lets take the average of the top and bottom and left/right boundary cells
834+
# apply the average to the boundary cells
835+
ufs[:,0,:] = (ufs[:,0,:]+ufs[:,-1,:])/2
836+
ufs[:,-1,:] = ufs[:,0,:]
837+
ufs[0,:,:] = (ufs[0,:,:]+ufs[-1,:,:])/2
838+
ufs[-1,:,:] = ufs[0,:,:]
839+
840+
ufn[:,0,:] = (ufn[:,0,:]+ufn[:,-1,:])/2
841+
ufn[:,-1,:] = ufn[:,0,:]
842+
ufn[0,:,:] = (ufn[0,:,:]+ufn[-1,:,:])/2
843+
ufn[-1,:,:] = ufn[0,:,:]
844+
845+
# now make sure that there is no gradients at the bondares
846+
ufs[:,1,:] = ufs[:,0,:]
847+
ufs[:,-2,:] = ufs[:,-1,:]
848+
ufs[1,:,:] = ufs[0,:,:]
849+
ufs[-2,:,:] = ufs[-1,:,:]
850+
851+
ufn[:,1,:] = ufn[:,0,:]
852+
ufn[:,-2,:] = ufn[:,-1,:]
853+
ufn[1,:,:] = ufn[0,:,:]
854+
ufn[-2,:,:] = ufn[-1,:,:]
855+
856+
# ufn[:,:,:] = ufn[-2,:,:]
857+
858+
# also correct for the potential gradients at the boundary cells in the equilibrium concentrations
859+
Cu[:,0,:] = Cu[:,1,:]
860+
Cu[:,-1,:] = Cu[:,-2,:]
861+
Cu[0,:,:] = Cu[1,:,:]
862+
Cu[-1,:,:] = Cu[-2,:,:]
863+
864+
# #boundary values
865+
# ufs[:,0, :] = us[:,0, :]
866+
# ufs[:,-1, :] = us[:,-1, :]
867+
868+
# ufn[0,:, :] = un[0,:, :]
869+
# ufn[-1,:, :] = un[-1,:, :]
833870

834871
Ct_last = Ct.copy()
835872
while k==0 or np.any(np.abs(Ct[:,:,i]-Ct_last[:,:,i])>1e-10):
@@ -982,8 +1019,264 @@ def sweep3(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w):
9821019

9831020

9841021
k+=1
1022+
1023+
# print(k)
9851024

9861025
# print("q1 = " + str(np.sum(q==1)) + " q2 = " + str(np.sum(q==2)) \
9871026
# + " q3 = " + str(np.sum(q==3)) + " q4 = " + str(np.sum(q==4)) \
9881027
# + " q5 = " + str(np.sum(q==5)))
1028+
# print("pickup deviation percentage = " + str(pickup.sum()/pickup[pickup>0].sum()*100) + " %")
1029+
# print("pickup deviation percentage = " + str(pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]>0,0].sum()*100) + " %")
1030+
# print("pickup maximum = " + str(pickup.max()) + " mass max = " + str(mass.max()))
1031+
# print("pickup minimum = " + str(pickup.min()))
1032+
# print("pickup average = " + str(pickup.mean()))
1033+
# print("number of cells for pickup maximum = " + str((pickup == mass.max()).sum()))
1034+
# pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]<0,0].sum()
1035+
1036+
return Ct, pickup
1037+
1038+
1039+
def sweep4(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w):
1040+
# this is where is make full circular boundaries
1041+
1042+
pickup = np.zeros(Cu.shape)
1043+
i=0
1044+
k=0
1045+
1046+
nf = np.shape(Ct)[2]
1047+
1048+
# Are the lateral boundary conditions circular?
1049+
circ_lateral = False
1050+
if Ct[0,1,0]==-1:
1051+
circ_lateral = True
1052+
Ct[0,:,0] = 0
1053+
Ct[-1,:,0] = 0
1054+
1055+
circ_offshore = False
1056+
if Ct[1,0,0]==-1:
1057+
circ_offshore = True
1058+
Ct[:,0,0] = 0
1059+
Ct[:,-1,0] = 0
1060+
1061+
recirc_offshore = False
1062+
if Ct[1,0,0]==-2:
1063+
recirc_offshore = True
1064+
Ct[:,0,0] = 0
1065+
Ct[:,-1,0] = 0
1066+
1067+
1068+
ufs = np.zeros((np.shape(us)[0], np.shape(us)[1]+1, np.shape(us)[2]))
1069+
ufn = np.zeros((np.shape(un)[0]+1, np.shape(un)[1], np.shape(un)[2]))
1070+
1071+
# define fluxes
1072+
ufs[:,1:-1, :] = 0.5*us[:,:-1, :] + 0.5*us[:,1:, :]
1073+
ufn[1:-1,:, :] = 0.5*un[:-1,:, :] + 0.5*un[1:,:, :]
1074+
1075+
# print(ufs[5,:,0])
1076+
1077+
ufs[:,0, :] = us[:,0, :]
1078+
ufs[:,-1, :] = us[:,-1, :]
1079+
1080+
ufn[0,:, :] = un[0,:, :]
1081+
ufn[-1,:, :] = un[-1,:, :]
1082+
1083+
# boundary values circular speed, taking the average of the top and bottom and left/right boundary cells
1084+
ufs[:,0,:] = (ufs[:,0,:]+ufs[:,-1,:])/2
1085+
ufs[:,-1,:] = ufs[:,0,:]
1086+
ufs[0,:,:] = (ufs[0,:,:]+ufs[-1,:,:])/2
1087+
ufs[-1,:,:] = ufs[0,:,:]
1088+
1089+
ufn[:,0,:] = (ufn[:,0,:]+ufn[:,-1,:])/2
1090+
ufn[:,-1,:] = ufn[:,0,:]
1091+
ufn[0,:,:] = (un[0,:,:]+ufn[-1,:,:])/2
1092+
ufn[-1,:,:] = ufn[0,:,:]
1093+
1094+
1095+
1096+
# ufs[:,0, :] = (us[:,0, :]+us[:,-1, :])/2
1097+
# ufs[:,-1, :] = ufs[:,0, :]
1098+
1099+
# ufs[:,0, :] = (us[:,0, :]+us[:,-1, :])/2
1100+
# ufs[:,-1, :] = ufs[:,0, :]
1101+
1102+
# ufs[:,0, :] = us[:,0, :]
1103+
# ufs[:,-1, :] = us[:,-1, :]
1104+
1105+
# ufn[0,:, :] = un[0,:, :]
1106+
# ufn[-1,:, :] = un[-1,:, :]
1107+
1108+
1109+
Ct_last = Ct.copy()
1110+
while k==0 or np.any(np.abs(Ct[:,:,i]-Ct_last[:,:,i])>1e-10):
1111+
# while k==0 or np.any(np.abs(Ct[:,:,i]-Ct_last[:,:,i])!=0):
1112+
Ct_last = Ct.copy()
1113+
1114+
# lateral boundaries circular
1115+
if circ_lateral:
1116+
Ct[0,:,0],Ct[-1,:,0] = Ct[-1,:,0].copy(),Ct[0,:,0].copy()
1117+
# ufn[0,:,0],ufn[-1,:,0] = ufn[-1,:,0].copy(),ufn[0,:,0].copy()
1118+
if circ_offshore:
1119+
Ct[:,0,0],Ct[:,-1,0] = Ct[:,-1,0].copy(),Ct[:,0,0].copy()
1120+
# ufs[0,:,0],ufs[-1,:,0] = ufs[-1,:,0].copy(),ufs[0,:,0].copy()
1121+
1122+
if recirc_offshore:
1123+
# print(Ct[:,1,0])
1124+
# print(Ct[:,-2,0])
1125+
Ct[:,0,0],Ct[:,-1,0] = np.average(Ct[:,-2,0]),np.average(Ct[:,1,0])
1126+
# print(Ct[:,0,0])
1127+
# print(Ct[:,-1,0])
1128+
1129+
# make an array with a bolean operator. This keeps track of considerd cells. We start with all False (not considered)
1130+
q = np.zeros(Cu.shape[:2])
1131+
1132+
########################################################################################
1133+
# in this sweeping algorithm we sweep over the 4 quadrants
1134+
# assuming that most cells have no converging/divering charactersitics.
1135+
# In the last quadrant we take converging and diverging cells into account.
1136+
1137+
# The First quadrant
1138+
1139+
nn = Ct.shape[0]
1140+
ns = Ct.shape[1]
1141+
1142+
for n in range(0,nn):
1143+
for s in range(0,ns):
1144+
if (not q[n,s]) and (ufn[n,s,0]>=0) and (ufs[n,s,0]>=0) and (ufn[n+1,s,0]>=0) and (ufs[n,s+1,0]>=0):
1145+
Ct[n,s,:] = (+ (Ct[n-1,s,:] * ufn[n,s,:] * ds[n,s]) \
1146+
+ (Ct[n,s-1,:] * ufs[n,s,:] * dn[n,s]) \
1147+
+ w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \
1148+
/ ( + (ufn[n+1,s,:] * ds[n,s]) \
1149+
+ (ufs[n,s+1,:] * dn[n,s]) \
1150+
+ (ds[n,s] * dn [n,s] / Ts) )
1151+
#calculate pickup
1152+
pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:] - Ct[n,s,:]) * dt/Ts
1153+
#check for supply limitations and re-iterate concentration to account for supply limitations
1154+
for nfs in range(0,nf):
1155+
if pickup[n,s,nfs]>mass[n,s,0,nfs]:
1156+
pickup[n,s,nfs] = mass[n,s,0,nfs]
1157+
Ct[n,s,nfs] = (+ (Ct[n-1,s,nfs] * ufn[n,s,nfs] * ds[n,s]) \
1158+
+ (Ct[n,s-1,nfs] * ufs[n,s,nfs] * dn[n,s]) \
1159+
+ pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \
1160+
/ (+(ufn[n+1,s,nfs] * ds[n,s]) \
1161+
+ (ufs[n,s+1,nfs] * dn[n,s]))
1162+
q[n,s]=1
1163+
1164+
# The second quadrant
1165+
for n in range(0,nn):
1166+
for s in range(ns-1,-1,-1):
1167+
if (not q[n,s]) and (ufn[n,s,0]>=0) and (ufs[n,s,0]<=0) and (ufn[n+1,s,0]>=0) and (ufs[n,s+1,0]<=0):
1168+
Ct[n,s,:] = (+ (Ct[n-1,s,:] * ufn[n,s,:] * ds[n,s]) \
1169+
+ ( -Ct[n,(s+1) % ns,:] * ufs[n,s+1,:] * dn[n,s]) \
1170+
+ w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \
1171+
/ ( + (ufn[n+1,s,:] * ds[n,s]) \
1172+
+ (-ufs[n,s,:] * dn[n,s]) \
1173+
+ (ds[n,s] * dn [n,s] / Ts) )
1174+
#calculate pickup
1175+
pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts
1176+
#check for supply limitations and re-iterate concentration to account for supply limitations
1177+
for nfs in range(0,nf):
1178+
if pickup[n,s,nfs]>mass[n,s,0,nfs]:
1179+
pickup[n,s,nfs] = mass[n,s,0,nfs]
1180+
Ct[n,s,nfs] = (+ (Ct[n-1,s,nfs] * ufn[n,s,nfs] * ds[n,s]) \
1181+
+ ( -Ct[n,(s+1) % ns,nfs] * ufs[n,s+1,nfs] * dn[n,s]) \
1182+
+ pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \
1183+
/ ( + (ufn[n+1,s,nfs] * ds[n,s]) \
1184+
+ (-ufs[n,s,nfs] * dn[n,s]))
1185+
q[n,s]=2
1186+
# The third quadrant
1187+
for n in range(nn-1,-1,-1):
1188+
for s in range(ns-1,-1,-1):
1189+
if (not q[n,s]) and (ufn[n,s,0]<=0) and (ufs[n,s,0]<=0) and (ufn[n+1,s,0]<=0) and (ufs[n,s+1,0]<=0):
1190+
Ct[n,s,:] = (+ ( -Ct[(n+1) % nn,s,:] * ufn[n+1,s,:] * dn[n,s]) \
1191+
+ ( -Ct[n,(s+1) % ns,:] * ufs[n,s+1,:] * dn[n,s]) \
1192+
+ w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \
1193+
/ ( + (-ufn[n,s,:] * dn[n,s]) \
1194+
+ (-ufs[n,s,:] * dn[n,s]) \
1195+
+ (ds[n,s] * dn [n,s] / Ts) )
1196+
#calculate pickup
1197+
pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts
1198+
#check for supply limitations and re-iterate concentration to account for supply limitations
1199+
for nfs in range(0,nf):
1200+
if pickup[n,s,nfs]>mass[n,s,0,nfs]:
1201+
pickup[n,s,nfs] = mass[n,s,0,nfs]
1202+
Ct[n,s,nfs] = (+ ( -Ct[(n+1) % nn,s,nfs] * ufn[n+1,s,nfs] * dn[n,s]) \
1203+
+ ( -Ct[n,(s+1) % ns,nfs] * ufs[n,s+1,nfs] * dn[n,s]) \
1204+
+ pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \
1205+
/ ( + (-ufn[n,s,nfs] * dn[n,s]) \
1206+
+ (-ufs[n,s,nfs] * dn[n,s]))
1207+
q[n,s]=3
1208+
# The fourth guadrant including all remainnig unadressed cells
1209+
for n in range(nn-1,-1,-1):
1210+
for s in range(0,ns):
1211+
if (not q[n,s]):
1212+
if (ufn[n,s,0]<=0) and (ufs[n,s,0]>=0) and (ufn[n+1,s,0]<=0) and (ufs[n,s+1,0]>=0):
1213+
# this is the fourth quadrant
1214+
Ct[n,s,:] = (+ (Ct[n,s-1,:] * ufs[n,s,:] * dn[n,s]) \
1215+
+ ( -Ct[(n+1) % nn,s,:] * ufn[n+1,s,:] * dn[n,s]) \
1216+
+ w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \
1217+
/ ( + (ufs[n,s+1,:] * dn[n,s]) \
1218+
+ (-ufn[n,s,:] * dn[n,s]) \
1219+
+ (ds[n,s] * dn [n,s] / Ts) )
1220+
#calculate pickup
1221+
pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts
1222+
#check for supply limitations and re-iterate concentration to account for supply limitations
1223+
for nfs in range(0,nf):
1224+
if pickup[n,s,nfs]>mass[n,s,0,nfs]:
1225+
pickup[n,s,nfs] = mass[n,s,0,nfs]
1226+
Ct[n,s,nfs] = (+ (Ct[n,s-1,nfs] * ufs[n,s,nfs] * dn[n,s]) \
1227+
+ ( -Ct[(n+1) % nn,s,nfs] * ufn[n+1,s,nfs] * dn[n,s]) \
1228+
+ pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \
1229+
/ ( + (ufs[n,s+1,nfs] * dn[n,s]) \
1230+
+ (-ufn[n,s,nfs] * dn[n,s]))
1231+
q[n,s]=4
1232+
else:
1233+
if True :
1234+
# if (not n==0) and (not s==Ct.shape[1]-1):
1235+
# This is where we apply a generic stencil where all posible directions on the cell boundaries are solved for.
1236+
# all remaining cells will be calculated for and q=5 is assigned.
1237+
# this stencil is nested in the q4 loop which is the final quadrant.
1238+
# grid boundaries are filtered in both if statements.
1239+
Ct[n,s,:] = (+ (ufn[n,s,0]>0) * (Ct[n-1,s,:] * ufn[n,s,:] * ds[n,s]) \
1240+
+ (ufs[n,s,0]>0) * (Ct[n,s-1,:] * ufs[n,s,:] * dn[n,s]) \
1241+
+ (ufn[n+1,s,0]<0) * ( -Ct[(n+1) % nn,s,:] * ufn[n+1,s,:] * dn[n,s]) \
1242+
+ (ufs[n,s+1,0]<0) * ( -Ct[n,(s+1) % ns,:] * ufs[n,s+1,:] * dn[n,s]) \
1243+
+ w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \
1244+
/ ( + (ufn[n+1,s,0]>0) * (ufn[n+1,s,:] * ds[n,s]) \
1245+
+ (ufs[n,s+1,0]>0) * (ufs[n,s+1,:] * dn[n,s]) \
1246+
+ (ufn[n,s,0]<0) * (-ufn[n,s,:] * dn[n,s]) \
1247+
+ (ufs[n,s,0]<0) * (-ufs[n,s,:] * dn[n,s]) \
1248+
+ (ds[n,s] * dn [n,s] / Ts) )
1249+
#calculate pickup
1250+
pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts
1251+
#check for supply limitations and re-iterate concentration to account for supply limitations
1252+
for nfs in range(0,nf):
1253+
if pickup[n,s,nfs]>mass[n,s,0,nfs]:
1254+
pickup[n,s,nfs] = mass[n,s,0,nfs]
1255+
Ct[n,s,nfs] = (+ (ufn[n,s,0]>0) * (Ct[n-1,s,nfs] * ufn[n,s,nfs] * ds[n,s]) \
1256+
+ (ufs[n,s,0]>0) * (Ct[n,s-1,nfs] * ufs[n,s,nfs] * dn[n,s]) \
1257+
+ (ufn[n+1,s,0]<0) * ( -Ct[(n+1) % nn,s,nfs] * ufn[n+1,s,nfs] * dn[n,s]) \
1258+
+ (ufs[n,s+1,0]<0) * ( -Ct[n,(s+1) % ns,nfs] * ufs[n,s+1,nfs] * dn[n,s]) \
1259+
+ pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \
1260+
/ ( + (ufn[n+1,s,0]>0) * (ufn[n+1,s,nfs] * ds[n,s]) \
1261+
+ (ufs[n,s+1,0]>0) * (ufs[n,s+1,nfs] * dn[n,s]) \
1262+
+ (ufn[n,s,0]<0) * (-ufn[n,s,nfs] * dn[n,s]) \
1263+
+ (ufs[n,s,0]<0) * (-ufs[n,s,nfs] * dn[n,s]))
1264+
q[n,s]=5
1265+
1266+
1267+
k+=1
1268+
1269+
print(k)
1270+
1271+
print("q1 = " + str(np.sum(q==1)) + " q2 = " + str(np.sum(q==2)) \
1272+
+ " q3 = " + str(np.sum(q==3)) + " q4 = " + str(np.sum(q==4)) \
1273+
+ " q5 = " + str(np.sum(q==5)))
1274+
print("pickup deviation percentage = " + str(pickup.sum()/pickup[pickup>0].sum()*100) + " %")
1275+
print("pickup deviation percentage = " + str(pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]>0,0].sum()*100) + " %")
1276+
print("pickup maximum = " + str(pickup.max()) + " mass max = " + str(mass.max()))
1277+
print("pickup minimum = " + str(pickup.min()))
1278+
print("pickup average = " + str(pickup.mean()))
1279+
print("number of cells for pickup maximum = " + str((pickup == mass.max()).sum()))
1280+
# pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]<0,0].sum()
1281+
9891282
return Ct, pickup

0 commit comments

Comments
 (0)