1 | % ================================================================ |
---|
2 | % Appendix E : UBS scheme |
---|
3 | % ================================================================ |
---|
4 | \chapter{UBS scheme} |
---|
5 | \label{AN_E} |
---|
6 | \minitoc |
---|
7 | |
---|
8 | |
---|
9 | % ------------------------------------------------------------------------------------------------------------- |
---|
10 | % UBS scheme |
---|
11 | % ------------------------------------------------------------------------------------------------------------- |
---|
12 | \section{Upstream Biased Scheme (UBS) (\np{ln\_traadv\_ubs}=T)} |
---|
13 | \label{TRA_adv_ubs} |
---|
14 | |
---|
15 | The UBS advection scheme is an upstream biased third order scheme based on |
---|
16 | an upstream-biased parabolic interpolation. It is also known as Cell Averaged |
---|
17 | QUICK scheme (Quadratic Upstream Interpolation for Convective |
---|
18 | Kinematics). For example, in the $i$-direction : |
---|
19 | \begin{equation} \label{Eq_tra_adv_ubs2} |
---|
20 | \tau _u^{ubs} = \left\{ \begin{aligned} |
---|
21 | & \tau _u^{cen4} + \frac{1}{12} \,\tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ |
---|
22 | & \tau _u^{cen4} - \frac{1}{12} \,\tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 |
---|
23 | \end{aligned} \right. |
---|
24 | \end{equation} |
---|
25 | or equivalently, the advective flux is |
---|
26 | \begin{equation} \label{Eq_tra_adv_ubs2} |
---|
27 | U_{i+1/2} \ \tau _u^{ubs} |
---|
28 | =U_{i+1/2} \ \overline{ T_i - \frac{1}{6}\,\tau"_i }^{\,i+1/2} |
---|
29 | - \frac{1}{2}\, |U|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] |
---|
30 | \end{equation} |
---|
31 | where $U_{i+1/2} = e_{1u}\,e_{3u}\,u_{i+1/2}$ and |
---|
32 | $\tau "_i =\delta _i \left[ {\delta _{i+1/2} \left[ \tau \right]} \right]$. |
---|
33 | By choosing this expression for $\tau "$ we consider a fourth order approximation |
---|
34 | of $\partial_i^2$ with a constant i-grid spacing ($\Delta i=1$). |
---|
35 | |
---|
36 | Alternative choice: introduce the scale factors: |
---|
37 | $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} }\delta _{i+1/2}[\tau] \right]$. |
---|
38 | |
---|
39 | |
---|
40 | This results in a dissipatively dominant (i.e. hyper-diffusive) truncation |
---|
41 | error \citep{Sacha2005}. The overall performance of the |
---|
42 | advection scheme is similar to that reported in \cite{Farrow1995}. |
---|
43 | It is a relatively good compromise between accuracy and smoothness. It is |
---|
44 | not a \emph{positive} scheme meaning false extrema are permitted but the |
---|
45 | amplitude of such are significantly reduced over the centred second order |
---|
46 | method. Nevertheless it is not recommended to apply it to a passive tracer |
---|
47 | that requires positivity. |
---|
48 | |
---|
49 | The intrinsic diffusion of UBS makes its use risky in the vertical direction |
---|
50 | where the control of artificial diapycnal fluxes is of paramount importance. |
---|
51 | It has therefore been preferred to evaluate the vertical flux using the TVD |
---|
52 | scheme when \np{ln\_traadv\_ubs}=T. |
---|
53 | |
---|
54 | For stability reasons, in \eqref{Eq_tra_adv_ubs}, the first term which corresponds |
---|
55 | to a second order centred scheme is evaluated using the \textit{now} velocity |
---|
56 | (centred in time) while the second term which is the diffusive part of the scheme, |
---|
57 | is evaluated using the \textit{before} velocity (forward in time. This is discussed |
---|
58 | by \citet{Webb1998} in the context of the Quick advection scheme. UBS and QUICK |
---|
59 | schemes only differ by one coefficient. Substituting 1/6 with 1/8 in |
---|
60 | (\ref{Eq_tra_adv_ubs}) leads to the QUICK advection scheme \citep{Webb1998}. |
---|
61 | This option is not available through a namelist parameter, since the 1/6 |
---|
62 | coefficient is hard coded. Nevertheless it is quite easy to make the |
---|
63 | substitution in \mdl{traadv\_ubs} module and obtain a QUICK scheme |
---|
64 | |
---|
65 | NB 1: When a high vertical resolution $O(1m)$ is used, the model stability can |
---|
66 | be controlled by vertical advection (not vertical diffusion which is usually |
---|
67 | solved using an implicit scheme). Computer time can be saved by using a |
---|
68 | time-splitting technique on vertical advection. This possibility have been |
---|
69 | implemented and validated in ORCA05-L301. It is not currently offered in the |
---|
70 | current reference version. |
---|
71 | |
---|
72 | NB 2 : In a forthcoming release four options will be proposed for the |
---|
73 | vertical component used in the UBS scheme. $\tau _w^{ubs}$ will be |
---|
74 | evaluated using either \textit{(a)} a centred $2^{nd}$ order scheme , |
---|
75 | or \textit{(b)} a TVD scheme, or \textit{(c)} an interpolation based on conservative |
---|
76 | parabolic splines following \citet{Sacha2005} implementation of UBS in ROMS, |
---|
77 | or \textit{(d)} an UBS. The $3^{rd}$ case has dispersion properties similar to an |
---|
78 | eight-order accurate conventional scheme. |
---|
79 | |
---|
80 | NB 3 : It is straight forward to rewrite \eqref{Eq_tra_adv_ubs} as follows: |
---|
81 | \begin{equation} \label{Eq_tra_adv_ubs2} |
---|
82 | \tau _u^{ubs} = \left\{ \begin{aligned} |
---|
83 | & \tau _u^{cen4} + \frac{1}{12} \tau"_i & \quad \text{if }\ u_{i+1/2} \geqslant 0 \\ |
---|
84 | & \tau _u^{cen4} - \frac{1}{12} \tau"_{i+1} & \quad \text{if }\ u_{i+1/2} < 0 |
---|
85 | \end{aligned} \right. |
---|
86 | \end{equation} |
---|
87 | or equivalently |
---|
88 | \begin{equation} \label{Eq_tra_adv_ubs2} |
---|
89 | \begin{split} |
---|
90 | e_{2u} e_{3u}\,u_{i+1/2} \ \tau _u^{ubs} |
---|
91 | &= e_{2u} e_{3u}\,u_{i+1/2} \ \overline{ T - \frac{1}{6}\,\tau"_i }^{\,i+1/2} \\ |
---|
92 | & - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] |
---|
93 | \end{split} |
---|
94 | \end{equation} |
---|
95 | \eqref{Eq_tra_adv_ubs2} has several advantages. First it clearly evidence that |
---|
96 | the UBS scheme is based on the fourth order scheme to which is added an |
---|
97 | upstream biased diffusive term. Second, this emphasises that the $4^{th}$ order |
---|
98 | part have to be evaluated at \emph{now} time step, not only the $2^{th}$ order |
---|
99 | part as stated above using \eqref{Eq_tra_adv_ubs}. Third, the diffusive term is |
---|
100 | in fact a biharmonic operator with a eddy coefficient with is simply proportional |
---|
101 | to the velocity. |
---|
102 | |
---|
103 | laplacian diffusion: |
---|
104 | \begin{equation} \label{Eq_tra_ldf_lap} |
---|
105 | \begin{split} |
---|
106 | D_T^{lT} =\frac{1}{e_{1T} \; e_{2T}\; e_{3T} } &\left[ {\quad \delta _i |
---|
107 | \left[ {A_u^{lT} \frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2} |
---|
108 | \left[ T \right]} \right]} \right. |
---|
109 | \\ |
---|
110 | &\ \left. {+\; \delta _j \left[ |
---|
111 | {A_v^{lT} \left( {\frac{e_{1v} e_{3v} }{e_{2v} }\;\delta _{j+1/2} \left[ T |
---|
112 | \right]} \right)} \right]\quad } \right] |
---|
113 | \end{split} |
---|
114 | \end{equation} |
---|
115 | |
---|
116 | bilaplacian: |
---|
117 | \begin{equation} \label{Eq_tra_ldf_lap} |
---|
118 | \begin{split} |
---|
119 | D_T^{lT} =&-\frac{1}{e_{1T} \; e_{2T}\; e_{3T}} \\ |
---|
120 | & \delta _i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta _{i+1/2} |
---|
121 | \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} |
---|
122 | \delta _i \left[ \sqrt{A_u^{lT}}\ \frac{e_{2u}\,e_{3u}}{e_{1u}}\;\delta _{i+1/2} |
---|
123 | [T] \right] \right] \right] |
---|
124 | \end{split} |
---|
125 | \end{equation} |
---|
126 | with ${A_u^{lT}}^2 = \frac{1}{12} {e_{1u}}^3\ |u|$, |
---|
127 | $i.e.$ $A_u^{lT} = \frac{1}{\sqrt{12}} \,e_{1u}\ \sqrt{ e_{1u}\,|u|\,}$ |
---|
128 | it comes : |
---|
129 | \begin{equation} \label{Eq_tra_ldf_lap} |
---|
130 | \begin{split} |
---|
131 | D_T^{lT} =&-\frac{1}{12}\,\frac{1}{e_{1T} \; e_{2T}\; e_{3T}} \\ |
---|
132 | & \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta _{i+1/2} |
---|
133 | \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} |
---|
134 | \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}\,|u|\,}\;\delta _{i+1/2} |
---|
135 | [T] \right] \right] \right] |
---|
136 | \end{split} |
---|
137 | \end{equation} |
---|
138 | if the velocity is uniform ($i.e.$ $|u|=cst$) then the diffusive flux is |
---|
139 | \begin{equation} \label{Eq_tra_ldf_lap} |
---|
140 | \begin{split} |
---|
141 | F_u^{lT} = - \frac{1}{12} |
---|
142 | e_{2u}\,e_{3u}\,|u| \;\sqrt{ e_{1u}}\,\delta _{i+1/2} |
---|
143 | \left[ \frac{1}{e_{1T}\,e_{2T}\, e_{3T}} |
---|
144 | \delta _i \left[ e_{2u}\,e_{3u}\,\sqrt{ e_{1u}}\:\delta _{i+1/2} |
---|
145 | [T] \right] \right] |
---|
146 | \end{split} |
---|
147 | \end{equation} |
---|
148 | beurk.... reverte the logic: starting from the diffusive part of the advective flux it comes: |
---|
149 | |
---|
150 | \begin{equation} \label{Eq_tra_adv_ubs2} |
---|
151 | \begin{split} |
---|
152 | F_u^{lT} |
---|
153 | &= - \frac{1}{2} e_{2u} e_{3u}\,|u|_{i+1/2} \;\frac{1}{6} \;\delta_{i+1/2}[\tau"_i] |
---|
154 | \end{split} |
---|
155 | \end{equation} |
---|
156 | if the velocity is uniform ($i.e.$ $|u|=cst$) and choosing $\tau "_i =\frac{e_{1T}}{e_{2T}\,e_{3T}}\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right]$ |
---|
157 | |
---|
158 | sol 1 coefficient at T-point ( add $e_{1u}$ and $e_{1T}$ on both side of first $\delta$): |
---|
159 | \begin{equation} \label{Eq_tra_adv_ubs2} |
---|
160 | \begin{split} |
---|
161 | F_u^{lT} |
---|
162 | &= - \frac{1}{12} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{e_{1T}^3\,|u|}{e_{1T}e_{2T}\,e_{3T}}\,\delta _i \left[ \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right] \right] |
---|
163 | \end{split} |
---|
164 | \end{equation} |
---|
165 | which leads to ${A_T^{lT}}^2 = \frac{1}{12} {e_{1T}}^3\ \overline{|u|}^{\,i+1/2}$ |
---|
166 | |
---|
167 | sol 2 coefficient at u-point: split $|u|$ into $\sqrt{|u|}$ and $e_{1T}$ into $\sqrt{e_{1u}}$ |
---|
168 | \begin{equation} \label{Eq_tra_adv_ubs2} |
---|
169 | \begin{split} |
---|
170 | F_u^{lT} |
---|
171 | &= - \frac{1}{12} {e_{1u}}^1 \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{2T}\,e_{3T}}\,\delta _i \left[ \sqrt{e_{1u}|u|} \frac{e_{2u} e_{3u} }{e_{1u} } \delta _{i+1/2}[\tau] \right] \right] \\ |
---|
172 | &= - \frac{1}{12} e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u}}{e_{1u}}\;\delta_{i+1/2}\left[ \frac{1}{e_{1T}\,e_{2T}\,e_{3T}}\,\delta _i \left[ e_{1u} \sqrt{e_{1u}|u|\,} \frac{e_{2u} e_{3u} }{e_{1u}} \delta _{i+1/2}[\tau] \right] \right] |
---|
173 | \end{split} |
---|
174 | \end{equation} |
---|
175 | which leads to ${A_u^{lT}} = \frac{1}{12} {e_{1u}}^3\ |u|$ |
---|
176 | |
---|
177 | |
---|
178 | % ------------------------------------------------------------------------------------------------------------- |
---|
179 | % Leap-Frog energetic |
---|
180 | % ------------------------------------------------------------------------------------------------------------- |
---|
181 | \section{Leap-Frog energetic } |
---|
182 | \label{LF} |
---|
183 | |
---|
184 | We adopt the following semi-discrete notation for time derivative. Given the values of a variable $q$ at successive time step, the time derivation and averaging operators at the mid time step are: |
---|
185 | \begin{subequations} \label{dt_mt} |
---|
186 | \begin{align} |
---|
187 | \delta _{t+\Delta t/2} [q] &= \ \ \, q^{t+\Delta t} - q^{t} \\ |
---|
188 | \overline q^{\,t+\Delta t/2} &= \left\{ q^{t+\Delta t} + q^{t} \right\} \; / \; 2 |
---|
189 | \end{align} |
---|
190 | \end{subequations} |
---|
191 | As for space operator, the adjoint of the derivation and averaging time operators are |
---|
192 | $\delta_t^*=\delta_{t+\Delta t/2}$ and $\overline{\cdot}^{\,t\,*}= \overline{\cdot}^{\,t+\Delta/2}$ |
---|
193 | , respectively. |
---|
194 | |
---|
195 | The Leap-frog time stepping given by \eqref{Eq_DOM_nxt} can be defined as: |
---|
196 | \begin{equation} \label{LF} |
---|
197 | \frac{\partial q}{\partial t} |
---|
198 | \equiv \frac{1}{\Delta t} \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} |
---|
199 | = \frac{q^{t+\Delta t}-q^{t-\Delta t}}{2\Delta t} |
---|
200 | \end{equation} |
---|
201 | Note that \eqref{LF} shows that the leapfrog time step is $\Delta t$, not $2\Delta t$ |
---|
202 | as it can be found sometime in literature. |
---|
203 | The leap-Frog time stepping is a second order centered scheme. As such it respects |
---|
204 | the quadratic invariant in integral forms, $i.e.$ the following continuous property, |
---|
205 | \begin{equation} \label{Energy} |
---|
206 | \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} |
---|
207 | =\int_{t_0}^{t_1} {\frac{1}{2}\, \frac{\partial q^2}{\partial t} \;dt} |
---|
208 | = \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) , |
---|
209 | \end{equation} |
---|
210 | is satisfied in discrete form. Indeed, |
---|
211 | \begin{equation} \begin{split} |
---|
212 | \int_{t_0}^{t_1} {q\, \frac{\partial q}{\partial t} \;dt} |
---|
213 | &\equiv \sum\limits_{0}^{N} |
---|
214 | {\frac{1}{\Delta t} q^t \ \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} \ \Delta t} |
---|
215 | \equiv \sum\limits_{0}^{N} { q^t \ \overline{ \delta _{t+\Delta t/2}[q]}^{\,t} } \\ |
---|
216 | &\equiv \sum\limits_{0}^{N} { \overline{q}^{\,t+\Delta/2}{ \delta _{t+\Delta t/2}[q]}} |
---|
217 | \equiv \sum\limits_{0}^{N} { \frac{1}{2} \delta _{t+\Delta t/2}[q^2] }\\ |
---|
218 | &\equiv \sum\limits_{0}^{N} { \frac{1}{2} \delta _{t+\Delta t/2}[q^2] } |
---|
219 | \equiv \frac{1}{2} \left( {q_{t_1}}^2 - {q_{t_0}}^2 \right) |
---|
220 | \end{split} \end{equation} |
---|
221 | NB here pb of boundary condition when applying the adjoin! In space, setting to 0 |
---|
222 | the quantity in land area is sufficient to get rid of the boundary condition |
---|
223 | (equivalently of the boundary value of the integration by part). In time this boundary |
---|
224 | condition is not physical and \textbf{add something here!!!} |
---|
225 | |
---|
226 | |
---|
227 | |
---|