From bcb7c8c26711da063221c7b0b38753b04e8d51b3 Mon Sep 17 00:00:00 2001 From: Sierd Date: Fri, 7 Feb 2025 17:20:40 +0100 Subject: [PATCH 1/3] updated circular boundaries for sweep solver --- aeolis/model.py | 15 +++++++++----- aeolis/netcdf.py | 4 ++-- aeolis/run_console.py | 5 ++++- aeolis/utils.py | 48 +++++++++++++++++++++++++++++++++++++++++-- 4 files changed, 62 insertions(+), 10 deletions(-) diff --git a/aeolis/model.py b/aeolis/model.py index cafb4427..c723c365 100644 --- a/aeolis/model.py +++ b/aeolis/model.py @@ -207,16 +207,21 @@ def initialize(self)-> None: if 1: # this is the new quasi 2D stuff # this is where we make the 1D grid into a 2D grid to ensure process compatibility - self.p['xgrid_file'] = np.transpose(np.stack((self.p['xgrid_file'], self.p['xgrid_file'], self.p['xgrid_file']),axis=1)) + + # define size of 2D grid 3 is minumum, variable defined for debugging purposes + qnr = 3 + + self.p['xgrid_file'] = np.transpose(np.stack([self.p['xgrid_file'] for i in range (qnr)],axis=1)) + + # redefine shape self.p['ny'], self.p['nx'] = self.p['xgrid_file'].shape + # repeat the above for ygrid_file assume dy is equal to dx dy = self.p['xgrid_file'][1,2]-self.p['xgrid_file'][1,1] - - # repeat the above for ygrid_file - 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)) + self.p['ygrid_file'] = np.transpose(np.stack([self.p['ygrid_file'] + i * dy for i in range(qnr)], axis=1)) # repeat the above for bed_file - self.p['bed_file'] = np.transpose(np.stack((self.p['bed_file'], self.p['bed_file'], self.p['bed_file']),axis=1)) + self.p['bed_file'] = np.transpose(np.stack([self.p['bed_file'] for i in range (qnr)],axis=1)) # change from number of points to number of cells self.p['nx'] -= 1 diff --git a/aeolis/netcdf.py b/aeolis/netcdf.py index 44b8d968..558e39af 100644 --- a/aeolis/netcdf.py +++ b/aeolis/netcdf.py @@ -297,8 +297,8 @@ def initialize(outputfile, outputvars, s, p, dimensions): # this is where a quasi 2D model is stored in 1D if p['ny']+1 == 3: nc.variables['n'][:] = np.arange(1) - nc.variables['x'][:,:] = s['x'][0,:] - nc.variables['y'][:,:] = s['y'][0,:] + nc.variables['x'][:,:] = s['x'][1,:] + nc.variables['y'][:,:] = s['y'][1,:] else: nc.variables['n'][:] = np.arange(p['ny']+1) nc.variables['x'][:,:] = s['x'] diff --git a/aeolis/run_console.py b/aeolis/run_console.py index 690e49ef..195b49f5 100644 --- a/aeolis/run_console.py +++ b/aeolis/run_console.py @@ -6,8 +6,11 @@ def main()-> None: '''Runs AeoLiS model in debugging mode.''' # configfile = r'c:\Users\weste_bt\aeolis\Tests\RotatingWind\Barchan_Grid270\aeolis.txt' - configfile = r'C:\Users\svries\Documents\GitHub\OE_aeolis-python\aeolis\examples\grainsizevariations\aeolis_horizontalgradient1.txt' + configfile = r'C:\Users\svries\Documents\GitHub\Bart_mass\aeolis_duran.txt' + # configfile = r'C:\Users\svries\Documents\GitHub\Bart_mass\aeolis_windspeed.txt' + aeolis_debug(configfile) + if __name__ == '__main__': main() diff --git a/aeolis/utils.py b/aeolis/utils.py index ece00753..3b9030e9 100644 --- a/aeolis/utils.py +++ b/aeolis/utils.py @@ -823,13 +823,47 @@ def sweep3(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w): ufn[1:-1,:, :] = 0.5*un[:-1,:, :] + 0.5*un[1:,:, :] # print(ufs[5,:,0]) - - #boundary values + + # boundary values ufs[:,0, :] = us[:,0, :] ufs[:,-1, :] = us[:,-1, :] ufn[0,:, :] = un[0,:, :] ufn[-1,:, :] = un[-1,:, :] + # first lets take the average of the top and bottom and left/right boundary cells + # apply the average to the boundary cells + ufs[:,0,:] = (ufs[:,0,:]+ufs[:,-1,:])/2 + ufs[:,-1,:] = ufs[:,0,:] + ufs[0,:,:] = (ufs[0,:,:]+ufs[-1,:,:])/2 + ufs[-1,:,:] = ufs[0,:,:] + + ufn[-1,:,:] = ufn[0,:,:] + + # now make sure that there is no gradients at the bondares + ufs[:,1,:] = ufs[:,0,:] + ufs[:,-2,:] = ufs[:,-1,:] + ufs[1,:,:] = ufs[0,:,:] + ufs[-2,:,:] = ufs[-1,:,:] + + ufn[:,1,:] = ufn[:,0,:] + ufn[:,-2,:] = ufn[:,-1,:] + ufn[1,:,:] = ufn[0,:,:] + ufn[-2,:,:] = ufn[-1,:,:] + + # ufn[:,:,:] = ufn[-2,:,:] + + # also correct for the potential gradients at the boundary cells in the equilibrium concentrations + Cu[:,0,:] = Cu[:,1,:] + Cu[:,-1,:] = Cu[:,-2,:] + Cu[0,:,:] = Cu[1,:,:] + Cu[-1,:,:] = Cu[-2,:,:] + + # #boundary values + # ufs[:,0, :] = us[:,0, :] + # ufs[:,-1, :] = us[:,-1, :] + + # ufn[0,:, :] = un[0,:, :] + # ufn[-1,:, :] = un[-1,:, :] Ct_last = Ct.copy() while k==0 or np.any(np.abs(Ct[:,:,i]-Ct_last[:,:,i])>1e-10): @@ -982,8 +1016,18 @@ def sweep3(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w): k+=1 + + # print(k) # print("q1 = " + str(np.sum(q==1)) + " q2 = " + str(np.sum(q==2)) \ # + " q3 = " + str(np.sum(q==3)) + " q4 = " + str(np.sum(q==4)) \ # + " q5 = " + str(np.sum(q==5))) + # print("pickup deviation percentage = " + str(pickup.sum()/pickup[pickup>0].sum()*100) + " %") + # print("pickup deviation percentage = " + str(pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]>0,0].sum()*100) + " %") + # print("pickup maximum = " + str(pickup.max()) + " mass max = " + str(mass.max())) + # print("pickup minimum = " + str(pickup.min())) + # print("pickup average = " + str(pickup.mean())) + # print("number of cells for pickup maximum = " + str((pickup == mass.max()).sum())) + # pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]<0,0].sum() + return Ct, pickup From c45aeb2b7f177d3609fe3d4e7f67bfad77436dab Mon Sep 17 00:00:00 2001 From: Sierd Date: Fri, 7 Feb 2025 17:21:41 +0100 Subject: [PATCH 2/3] Update utils.py --- aeolis/utils.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/aeolis/utils.py b/aeolis/utils.py index 3b9030e9..1b700c48 100644 --- a/aeolis/utils.py +++ b/aeolis/utils.py @@ -837,6 +837,9 @@ def sweep3(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w): ufs[0,:,:] = (ufs[0,:,:]+ufs[-1,:,:])/2 ufs[-1,:,:] = ufs[0,:,:] + ufn[:,0,:] = (ufn[:,0,:]+ufn[:,-1,:])/2 + ufn[:,-1,:] = ufn[:,0,:] + ufn[0,:,:] = (ufn[0,:,:]+ufn[-1,:,:])/2 ufn[-1,:,:] = ufn[0,:,:] # now make sure that there is no gradients at the bondares From 80f21f23f885fbfa4cfbf149456132d92adfd46e Mon Sep 17 00:00:00 2001 From: Sierd Date: Fri, 7 Feb 2025 17:23:02 +0100 Subject: [PATCH 3/3] Added 4th generation sweep solver for testing --- aeolis/utils.py | 246 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 246 insertions(+) diff --git a/aeolis/utils.py b/aeolis/utils.py index 1b700c48..d801a19f 100644 --- a/aeolis/utils.py +++ b/aeolis/utils.py @@ -1034,3 +1034,249 @@ def sweep3(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w): # pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]<0,0].sum() return Ct, pickup + + +def sweep4(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w): + # this is where is make full circular boundaries + + pickup = np.zeros(Cu.shape) + i=0 + k=0 + + nf = np.shape(Ct)[2] + + # Are the lateral boundary conditions circular? + circ_lateral = False + if Ct[0,1,0]==-1: + circ_lateral = True + Ct[0,:,0] = 0 + Ct[-1,:,0] = 0 + + circ_offshore = False + if Ct[1,0,0]==-1: + circ_offshore = True + Ct[:,0,0] = 0 + Ct[:,-1,0] = 0 + + recirc_offshore = False + if Ct[1,0,0]==-2: + recirc_offshore = True + Ct[:,0,0] = 0 + Ct[:,-1,0] = 0 + + + ufs = np.zeros((np.shape(us)[0], np.shape(us)[1]+1, np.shape(us)[2])) + ufn = np.zeros((np.shape(un)[0]+1, np.shape(un)[1], np.shape(un)[2])) + + # define fluxes + ufs[:,1:-1, :] = 0.5*us[:,:-1, :] + 0.5*us[:,1:, :] + ufn[1:-1,:, :] = 0.5*un[:-1,:, :] + 0.5*un[1:,:, :] + + # print(ufs[5,:,0]) + + ufs[:,0, :] = us[:,0, :] + ufs[:,-1, :] = us[:,-1, :] + + ufn[0,:, :] = un[0,:, :] + ufn[-1,:, :] = un[-1,:, :] + + # boundary values circular speed, taking the average of the top and bottom and left/right boundary cells + ufs[:,0,:] = (ufs[:,0,:]+ufs[:,-1,:])/2 + ufs[:,-1,:] = ufs[:,0,:] + ufs[0,:,:] = (ufs[0,:,:]+ufs[-1,:,:])/2 + ufs[-1,:,:] = ufs[0,:,:] + + ufn[:,0,:] = (ufn[:,0,:]+ufn[:,-1,:])/2 + ufn[:,-1,:] = ufn[:,0,:] + ufn[0,:,:] = (un[0,:,:]+ufn[-1,:,:])/2 + ufn[-1,:,:] = ufn[0,:,:] + + + + # ufs[:,0, :] = (us[:,0, :]+us[:,-1, :])/2 + # ufs[:,-1, :] = ufs[:,0, :] + + # ufs[:,0, :] = (us[:,0, :]+us[:,-1, :])/2 + # ufs[:,-1, :] = ufs[:,0, :] + + # ufs[:,0, :] = us[:,0, :] + # ufs[:,-1, :] = us[:,-1, :] + + # ufn[0,:, :] = un[0,:, :] + # ufn[-1,:, :] = un[-1,:, :] + + + Ct_last = Ct.copy() + while k==0 or np.any(np.abs(Ct[:,:,i]-Ct_last[:,:,i])>1e-10): + # while k==0 or np.any(np.abs(Ct[:,:,i]-Ct_last[:,:,i])!=0): + Ct_last = Ct.copy() + + # lateral boundaries circular + if circ_lateral: + Ct[0,:,0],Ct[-1,:,0] = Ct[-1,:,0].copy(),Ct[0,:,0].copy() + # ufn[0,:,0],ufn[-1,:,0] = ufn[-1,:,0].copy(),ufn[0,:,0].copy() + if circ_offshore: + Ct[:,0,0],Ct[:,-1,0] = Ct[:,-1,0].copy(),Ct[:,0,0].copy() + # ufs[0,:,0],ufs[-1,:,0] = ufs[-1,:,0].copy(),ufs[0,:,0].copy() + + if recirc_offshore: + # print(Ct[:,1,0]) + # print(Ct[:,-2,0]) + Ct[:,0,0],Ct[:,-1,0] = np.average(Ct[:,-2,0]),np.average(Ct[:,1,0]) + # print(Ct[:,0,0]) + # print(Ct[:,-1,0]) + + # make an array with a bolean operator. This keeps track of considerd cells. We start with all False (not considered) + q = np.zeros(Cu.shape[:2]) + + ######################################################################################## + # in this sweeping algorithm we sweep over the 4 quadrants + # assuming that most cells have no converging/divering charactersitics. + # In the last quadrant we take converging and diverging cells into account. + + # The First quadrant + + nn = Ct.shape[0] + ns = Ct.shape[1] + + for n in range(0,nn): + for s in range(0,ns): + 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): + Ct[n,s,:] = (+ (Ct[n-1,s,:] * ufn[n,s,:] * ds[n,s]) \ + + (Ct[n,s-1,:] * ufs[n,s,:] * dn[n,s]) \ + + w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \ + / ( + (ufn[n+1,s,:] * ds[n,s]) \ + + (ufs[n,s+1,:] * dn[n,s]) \ + + (ds[n,s] * dn [n,s] / Ts) ) + #calculate pickup + pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:] - Ct[n,s,:]) * dt/Ts + #check for supply limitations and re-iterate concentration to account for supply limitations + for nfs in range(0,nf): + if pickup[n,s,nfs]>mass[n,s,0,nfs]: + pickup[n,s,nfs] = mass[n,s,0,nfs] + Ct[n,s,nfs] = (+ (Ct[n-1,s,nfs] * ufn[n,s,nfs] * ds[n,s]) \ + + (Ct[n,s-1,nfs] * ufs[n,s,nfs] * dn[n,s]) \ + + pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \ + / (+(ufn[n+1,s,nfs] * ds[n,s]) \ + + (ufs[n,s+1,nfs] * dn[n,s])) + q[n,s]=1 + + # The second quadrant + for n in range(0,nn): + for s in range(ns-1,-1,-1): + 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): + Ct[n,s,:] = (+ (Ct[n-1,s,:] * ufn[n,s,:] * ds[n,s]) \ + + ( -Ct[n,(s+1) % ns,:] * ufs[n,s+1,:] * dn[n,s]) \ + + w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \ + / ( + (ufn[n+1,s,:] * ds[n,s]) \ + + (-ufs[n,s,:] * dn[n,s]) \ + + (ds[n,s] * dn [n,s] / Ts) ) + #calculate pickup + pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts + #check for supply limitations and re-iterate concentration to account for supply limitations + for nfs in range(0,nf): + if pickup[n,s,nfs]>mass[n,s,0,nfs]: + pickup[n,s,nfs] = mass[n,s,0,nfs] + Ct[n,s,nfs] = (+ (Ct[n-1,s,nfs] * ufn[n,s,nfs] * ds[n,s]) \ + + ( -Ct[n,(s+1) % ns,nfs] * ufs[n,s+1,nfs] * dn[n,s]) \ + + pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \ + / ( + (ufn[n+1,s,nfs] * ds[n,s]) \ + + (-ufs[n,s,nfs] * dn[n,s])) + q[n,s]=2 + # The third quadrant + for n in range(nn-1,-1,-1): + for s in range(ns-1,-1,-1): + 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): + Ct[n,s,:] = (+ ( -Ct[(n+1) % nn,s,:] * ufn[n+1,s,:] * dn[n,s]) \ + + ( -Ct[n,(s+1) % ns,:] * ufs[n,s+1,:] * dn[n,s]) \ + + w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \ + / ( + (-ufn[n,s,:] * dn[n,s]) \ + + (-ufs[n,s,:] * dn[n,s]) \ + + (ds[n,s] * dn [n,s] / Ts) ) + #calculate pickup + pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts + #check for supply limitations and re-iterate concentration to account for supply limitations + for nfs in range(0,nf): + if pickup[n,s,nfs]>mass[n,s,0,nfs]: + pickup[n,s,nfs] = mass[n,s,0,nfs] + Ct[n,s,nfs] = (+ ( -Ct[(n+1) % nn,s,nfs] * ufn[n+1,s,nfs] * dn[n,s]) \ + + ( -Ct[n,(s+1) % ns,nfs] * ufs[n,s+1,nfs] * dn[n,s]) \ + + pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \ + / ( + (-ufn[n,s,nfs] * dn[n,s]) \ + + (-ufs[n,s,nfs] * dn[n,s])) + q[n,s]=3 + # The fourth guadrant including all remainnig unadressed cells + for n in range(nn-1,-1,-1): + for s in range(0,ns): + if (not q[n,s]): + 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): + # this is the fourth quadrant + Ct[n,s,:] = (+ (Ct[n,s-1,:] * ufs[n,s,:] * dn[n,s]) \ + + ( -Ct[(n+1) % nn,s,:] * ufn[n+1,s,:] * dn[n,s]) \ + + w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \ + / ( + (ufs[n,s+1,:] * dn[n,s]) \ + + (-ufn[n,s,:] * dn[n,s]) \ + + (ds[n,s] * dn [n,s] / Ts) ) + #calculate pickup + pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts + #check for supply limitations and re-iterate concentration to account for supply limitations + for nfs in range(0,nf): + if pickup[n,s,nfs]>mass[n,s,0,nfs]: + pickup[n,s,nfs] = mass[n,s,0,nfs] + Ct[n,s,nfs] = (+ (Ct[n,s-1,nfs] * ufs[n,s,nfs] * dn[n,s]) \ + + ( -Ct[(n+1) % nn,s,nfs] * ufn[n+1,s,nfs] * dn[n,s]) \ + + pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \ + / ( + (ufs[n,s+1,nfs] * dn[n,s]) \ + + (-ufn[n,s,nfs] * dn[n,s])) + q[n,s]=4 + else: + if True : + # if (not n==0) and (not s==Ct.shape[1]-1): + # This is where we apply a generic stencil where all posible directions on the cell boundaries are solved for. + # all remaining cells will be calculated for and q=5 is assigned. + # this stencil is nested in the q4 loop which is the final quadrant. + # grid boundaries are filtered in both if statements. + Ct[n,s,:] = (+ (ufn[n,s,0]>0) * (Ct[n-1,s,:] * ufn[n,s,:] * ds[n,s]) \ + + (ufs[n,s,0]>0) * (Ct[n,s-1,:] * ufs[n,s,:] * dn[n,s]) \ + + (ufn[n+1,s,0]<0) * ( -Ct[(n+1) % nn,s,:] * ufn[n+1,s,:] * dn[n,s]) \ + + (ufs[n,s+1,0]<0) * ( -Ct[n,(s+1) % ns,:] * ufs[n,s+1,:] * dn[n,s]) \ + + w[n,s,:] * Cu[n,s,:] * ds[n,s] * dn [n,s] / Ts ) \ + / ( + (ufn[n+1,s,0]>0) * (ufn[n+1,s,:] * ds[n,s]) \ + + (ufs[n,s+1,0]>0) * (ufs[n,s+1,:] * dn[n,s]) \ + + (ufn[n,s,0]<0) * (-ufn[n,s,:] * dn[n,s]) \ + + (ufs[n,s,0]<0) * (-ufs[n,s,:] * dn[n,s]) \ + + (ds[n,s] * dn [n,s] / Ts) ) + #calculate pickup + pickup[n,s,:] = (w[n,s,:] * Cu[n,s,:]-Ct[n,s,:]) * dt/Ts + #check for supply limitations and re-iterate concentration to account for supply limitations + for nfs in range(0,nf): + if pickup[n,s,nfs]>mass[n,s,0,nfs]: + pickup[n,s,nfs] = mass[n,s,0,nfs] + Ct[n,s,nfs] = (+ (ufn[n,s,0]>0) * (Ct[n-1,s,nfs] * ufn[n,s,nfs] * ds[n,s]) \ + + (ufs[n,s,0]>0) * (Ct[n,s-1,nfs] * ufs[n,s,nfs] * dn[n,s]) \ + + (ufn[n+1,s,0]<0) * ( -Ct[(n+1) % nn,s,nfs] * ufn[n+1,s,nfs] * dn[n,s]) \ + + (ufs[n,s+1,0]<0) * ( -Ct[n,(s+1) % ns,nfs] * ufs[n,s+1,nfs] * dn[n,s]) \ + + pickup[n,s,nfs] * ds[n,s] * dn [n,s] / dt ) \ + / ( + (ufn[n+1,s,0]>0) * (ufn[n+1,s,nfs] * ds[n,s]) \ + + (ufs[n,s+1,0]>0) * (ufs[n,s+1,nfs] * dn[n,s]) \ + + (ufn[n,s,0]<0) * (-ufn[n,s,nfs] * dn[n,s]) \ + + (ufs[n,s,0]<0) * (-ufs[n,s,nfs] * dn[n,s])) + q[n,s]=5 + + + k+=1 + + print(k) + + print("q1 = " + str(np.sum(q==1)) + " q2 = " + str(np.sum(q==2)) \ + + " q3 = " + str(np.sum(q==3)) + " q4 = " + str(np.sum(q==4)) \ + + " q5 = " + str(np.sum(q==5))) + print("pickup deviation percentage = " + str(pickup.sum()/pickup[pickup>0].sum()*100) + " %") + print("pickup deviation percentage = " + str(pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]>0,0].sum()*100) + " %") + print("pickup maximum = " + str(pickup.max()) + " mass max = " + str(mass.max())) + print("pickup minimum = " + str(pickup.min())) + print("pickup average = " + str(pickup.mean())) + print("number of cells for pickup maximum = " + str((pickup == mass.max()).sum())) + # pickup[1,:,0].sum()/pickup[1,pickup[1,:,0]<0,0].sum() + + return Ct, pickup