Skip to content

Commit 80f21f2

Browse files
SierdSierd
committed
Added 4th generation sweep solver for testing
1 parent c45aeb2 commit 80f21f2

File tree

1 file changed

+246
-0
lines changed

1 file changed

+246
-0
lines changed

‎aeolis/utils.py‎

Lines changed: 246 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1034,3 +1034,249 @@ def sweep3(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w):
10341034
# pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]<0,0].sum()
10351035

10361036
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+
1282+
return Ct, pickup

0 commit comments

Comments
 (0)