source: codes/icosagcm/devel/Python/dynamico/advect.py @ 615

Last change on this file since 615 was 615, checked in by dubos, 6 years ago

devel/unstructured : import Python layer from HEAT
File size: 7.2 KB
Line 
1import numpy as np
2from scipy import sparse as sparse
3
4#-------------------------------- Strang-split advection --------------------------------
5
6class SplitTransport:
7    def __init__(self, metric, Horiz, Vert):
8        # Horiz and Vert are classes for horizontal and vertical transport, respectively
9        self.horiz = Horiz(metric)
10        self.vert = Vert(metric)
11        self.di, self.dk = metric.di, metric.dk
12    def next(self, Fjk, Fil,  # time-integrated mass fluxes
13             mik,             # mass field before advection
14             scalars):        # tuple of 2D mass-weighted scalars to be advected
15        scalars =  [Qik for Qik in self.vert.advect(.5*Fil,mik,scalars)]
16        mik = mik - .5*self.dk(Fil)
17        scalars =  [Qik for Qik in self.horiz.advect(Fjk,mik,scalars)]
18        mik = mik - self.di(Fjk)
19        scalars =  [Qik for Qik in self.vert.advect(.5*Fil,mik,scalars)]
20        mik = mik - .5*self.dk(Fil)
21        return mik, scalars
22   
23#---------------------------- Horizontal advection schemes ----------------------------
24
25class HorizGodunov:
26    def __init__(self,metric):
27        self.di, self.dj = metric.di, metric.dj
28        self.mi, self.mj = metric.mi, metric.mj
29    # member function reconstruct() (to be overloaded by derived classes)
30    # yields for each point j the mean and difference of right and left reconstructions
31    def reconstruct(self,Fjk,maik,scalars):
32        for Qik in scalars:
33            qik = Qik/maik
34            yield Qik, self.mj(qik), self.dj(qik) # piecewise-constant reconstruction
35    def advect(self,Fjk,maik,scalars):
36        for Qik, mean_qr_ql, qr_minus_ql in self.reconstruct(Fjk,maik,scalars):
37            qjk = mean_qr_ql - 0.5*np.sign(Fjk)*qr_minus_ql  # Upwind sub-grid reconstruction
38            yield Qik - self.di(qjk*Fjk)
39           
40class HorizVanLeer(HorizGodunov):
41    def reconstruct(self,Fjk,mik,scalars):
42        fjk = .5*Fjk/self.mj(mik) # displacement, same for all scalars
43        for Qik in scalars:
44            qik = Qik/mik
45            delqik = self.di(self.mj(qik))
46            mean_qr_ql =   self.mj(qik) - .25*self.dj(delqik) - fjk*self.mj(delqik)
47            qr_minus_ql =  self.dj(qik) - self.mj(delqik)     - fjk*self.dj(delqik)
48            yield Qik, mean_qr_ql, qr_minus_ql
49           
50#----------------------------- Vertical advection schemes ---------------------------------
51
52class VertGodunov:
53    def __init__(self,metric):
54        self.dk, self.dl, self.ml = metric.dk, metric.dl, metric.ml_int
55    # member function reconstruct() (to be overloaded by derived classes)
56    # yields for each point j the mean and difference of right and left reconstructions
57    def advect(self,Fil,mik,scalars):
58        for Qik, mean_qr_ql, qr_minus_ql in self.reconstruct(Fil,mik,scalars):
59            qil = mean_qr_ql - 0.5*np.sign(Fil)*qr_minus_ql  # Upwind sub-grid reconstruction
60            yield Qik - self.dk(qil*Fil)
61    def reconstruct(self,Fil,mik,scalars):
62        for Qik in scalars:
63            qik = Qik/mik
64            yield Qik, self.ml(qik), self.dl(0.,qik,0.)
65
66class VertVanLeer(VertGodunov):
67    def reconstruct(self,Fil,mik,scalars):
68        fil = .5*np.array(Fil/self.ml(mik)) # displacement, same for all scalars
69        for Qik in scalars:
70            qik = Qik/mik
71            qil = self.ml(qik)
72            delqik = self.dk(qil) # subgrid slope, FIXME at top and bottom
73            mean_qr_ql  = self.ml(qik)     - .25*self.dl(0,delqik,0) - fil*self.ml(delqik)
74            qr_minus_ql = self.dl(0,qik,0) - self.ml(delqik)         - fil*self.dl(0,delqik,0)
75            yield Qik, mean_qr_ql, qr_minus_ql
76   
77#--------------------- Sarvesh original code -------------------------
78
79class adv:
80    def __init__(self,metric,maik,qik,Uajk,Fail):
81        self.metric = metric
82        self.maik=maik
83        self.qik=qik
84        self.Uajk=Uajk
85        self.Fail=Fail
86        self.di, self.dj = metric.di, metric.dj
87        self.dk, self.dl = metric.dk, metric.dl
88        self.mi, self.mj = metric.mi, metric.mj
89        Nx, Ny = np.shape(maik)
90        ones=np.ones(Nx);
91        left = sparse.spdiags([ones],[-1],Nx,Nx, format='csr')
92        left[0,-1]=1
93        right = sparse.spdiags([ones],[1],Nx,Nx, format='csr')
94        right[-1,0]=1
95        self.left =lambda f:left.dot(f)
96        self.right =lambda f:right.dot(f)
97        self.mk_ext, self.ml_ext = metric.mk_ext, metric.ml_ext
98        self.mk_int, self.ml_int = metric.mk_int, metric.ml_int
99
100    def adv_horizon(self,maik,qik,Uajk):
101        maik_old = maik  # Mass at the start of the advection (before N-adv dynamical steps) 
102        qik_r = self.right(qik)
103        maik_r = self.right(maik)
104        qik_l = self.left(qik)
105        delqik = 0.5*(qik_r-qik_l)
106        qr = qik + 0.5*(1 - Uajk/maik)*delqik
107        ql = qik_r - 0.5*(1 + Uajk/maik_r)*self.right(delqik)
108        qUPjk= 0.5*((qr+ql) + np.sign(Uajk)*(qr-ql))  ##Upwind sub-grid reconstruction
109        qxflx = qUPjk*Uajk
110        maik = maik - self.di(Uajk) ## Updated Mass
111        qik = (qik*maik_old - self.di(qxflx))/maik
112        return maik,qik
113   
114    def adv_horizon1(self,maik,qik,Uajk):
115        maik_old = maik  # Mass at the start of the advection (before N-adv dynamical steps) 
116        qjk = self.mj(qik)
117        delqik = self.di(qjk)
118        qik_r = qik + 0.5*(1 - Uajk/maik)*delqik
119        qik_l = np.zeros(np.shape(qik))
120        qik_l[:-1,:] = qik[1:,:]-0.5*(1.0 - Uajk[:-1,:]/maik[1:,:])*delqik[1:,:]
121        qik_l[-1,:] = qik[0,:] - 0.5*(1.0 - Uajk[-1,:]/maik[0,:])*delqik[0,:]
122        qUPjk= 0.5*((qik_r+qik_l) + np.sign(Uajk[:,:])*(qik_r-qik_l))
123        qxflx = qUPjk*Uajk
124        maik = maik - self.di(Uajk) ## Updated Mass
125        qik = (qik*maik_old - self.di(qxflx))/maik
126        return maik,qik
127   
128    def adv_godnov(self,maik,qik,Uajk):
129        maik_old = maik  # Mass at the start of the advection (before N-adv dynamical steps) 
130        qjk = self.mj(qik)
131        qik_r = qik
132        qik_l = np.zeros(np.shape(qik))
133        qik_l[:-1,:] = qik[1:,:]
134        qik_l[-1,:] = qik[0,:] 
135        qUPjk= 0.5*((qik_r+qik_l) + np.sign(Uajk[:,:])*(qik_r-qik_l))
136        qxflx = qUPjk*Uajk
137        maik = maik - self.di(Uajk) ## Updated Mass
138        qik = (qik*maik_old - self.di(qxflx))/maik
139        return maik,qik
140   
141    def adv_vert(self,maik,qik,Fil):
142        maik_old = maik
143        qil = self.ml_int(qik)
144        delqik = self.dk(qil) # subgrid slope
145        qrik = qik[:,:-1] + 0.5*(1 - Fil[:,1:-1]/maik[:,:-1])*delqik[:,:-1]
146        qlik = qik[:,1:] - 0.5*(1 - Fil[:,1:-1]/maik[:,1:])*delqik[:,1:]
147     
148        qUPil = np.zeros(np.shape(Fil))
149        qUPil[:,1:-1]= 0.5*((qrik+qlik) + np.sign(Fil[:,1:-1])*(qrik-qlik))
150        qUPil[:,0]= 0.0 #qik[:,0] - 0.5*(1 - Fil[:,0]/maik[:,0])*delqik[:,0]
151        qUPil[:,-1]= 0.0 #qik[:,-1] + 0.5*(1 - Fil[:,-1]/maik[:,-1])*delqik[:,-1]
152        qzflx = qUPil*Fil
153        maik = maik - self.dk(Fil)
154        qik = (qik*maik_old - self.dk(qzflx))/maik
155        return maik,qik
156   
157    def adv_vgod(self,maik,qik,Fil):
158        maik_old = maik
159        qrik = qik[:,:-1]
160        qlik = qik[:,1:]       
161        qUPil = np.zeros(np.shape(Fil))
162        qUPil[:,1:-1]= 0.5*((qrik+qlik) + np.sign(Fil[:,1:-1])*(qrik-qlik))
163        qUPil[:,0]= 0.0
164        qUPil[:,-1]= 0.0
165        qzflx = qUPil*Fil
166        maik = maik - self.dk(Fil)
167        qik = (qik - self.dk(qzflx)) #/maik
168        return maik,qik
Note: See TracBrowser for help on using the repository browser.