 Timestamp:
 Oct 8, 2018, 3:34:53 PM (20 months ago)
 Branches:
 master, dev, feature/makefiles, task/languageediting
 Children:
 0e576fc
 Parents:
 5bad6ec
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

source/python/mods/disaggregation.py
r25b14be rae88f7d 40 40 #  dapoly 41 41 #  darain 42 #  IA3 42 43 # 43 44 #******************************************************************************* … … 140 141 141 142 return nfield 143 144 def IA3(g): 145 """ 146 147 *************************************************************************** 148 * Copyright 2017 * 149 * Sabine Hittmeir, Anne Philipp, Petra Seibert * 150 * * 151 * This work is licensed under the Creative Commons Attribution 4.0 * 152 * International License. To view a copy of this license, visit * 153 * http://creativecommons.org/licenses/by/4.0/ or send a letter to * 154 * Creative Commons, PO Box 1866, Mountain View, CA 94042, USA. * 155 *************************************************************************** 156 157 @Description 158 The given data series will be interpolated with a nonnegative geometric 159 mean based algorithm. The original grid is reconstructed by adding two 160 sampling points in each data series interval. This subgrid is used to 161 keep all information during the interpolation within the associated 162 interval. Additionally, an advanced monotonicity filter is applied to 163 improve the monotonicity properties of the series. 164 165 For more information see article: 166 Hittmeir, S.; Philipp, A.; Seibert, P. (2017): A conservative 167 interpolation scheme for extensive quantities with application to the 168 Lagrangian particle dispersion model FLEXPART., 169 Geoscientific Model Development 170 171 @Input 172 g: list of float values 173 A list of float values which represents the complete data series that 174 will be interpolated having the dimension of the original raw series. 175 176 @Return 177 f: list of float values 178 The interpolated data series with additional subgrid points. 179 Its dimension is equal to the length of the input data series 180 times three. 181 """ 182 183 ####################### variable description ############################# 184 # # 185 # i  index variable for looping over the data series # 186 # g  input data series # 187 # f  interpolated and filtered data series with additional # 188 # grid points # 189 # fi  function value at position i, f_i # 190 # fi1  first subgrid function value f_i^1 # 191 # fi2  second subgrid function value f_i^2 # 192 # fip1  next function value at position i+1, f_(i+1) # 193 # dt  time step # 194 # fmon  monotonicity filter # 195 # # 196 ########################################################################### 197 198 199 import numpy as np 200 201 # time step 202 dt=1.0 203 204 ############### Nonnegative Geometric Mean Based Algorithm ############### 205 206 # for the left boundary the following boundary condition is valid: 207 # the value at t=0 of the interpolation algorithm coincides with the 208 # first data value according to the persistence hypothesis 209 f=[g[0]] 210 211 # compute two first subgrid intervals without monotonicity check 212 # go through the data series and extend each interval by two subgrid 213 # points and interpolate the corresponding data values 214 # except for the last interval due to boundary conditions 215 for i in range(0,2): 216 217 # as a requirement: 218 # if there is a zero data value such that g[i]=0, then the whole 219 # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0 220 # according to Eq. (6) 221 if g[i]==0.: 222 f.extend([0.,0.,0.]) 223 224 # otherwise the subgrid values are calculated and added to the list 225 else: 226 # temporal save of last value in interpolated list 227 # since it is the left boundary and hence the new (fi) value 228 fi = f[1] 229 230 # the value at the end of the interval (fip1) is prescribed by the 231 # geometric mean, restricted such that nonnegativity is guaranteed 232 # according to Eq. (25) 233 fip1=min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) ) 234 235 # the function value at the first subgrid point (fi1) is determined 236 # according to the equal area condition with Eq. (19) 237 fi1=3./2.*g[i]5./12.*fip11./12.*fi 238 239 # the function value at the second subgrid point (fi2) is determined 240 # according Eq. (18) 241 fi2=fi1+1./3.*(fip1fi) 242 243 # add next interval of interpolated (sub)grid values 244 f.append(fi1) 245 f.append(fi2) 246 f.append(fip1) 247 248 # compute rest of the data series intervals 249 # go through the data series and extend each interval by two subgrid 250 # points and interpolate the corresponding data values 251 # except for the last interval due to boundary conditions 252 for i in range(2,len(g)1): 253 254 # as a requirement: 255 # if there is a zero data value such that g[i]=0, then the whole 256 # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0 257 # according to Eq. (6) 258 if g[i]==0.: 259 # apply monotonicity filter for interval before 260 # check if there is "M" or "W" shape 261 if np.sign(f[5]f[6]) * np.sign(f[4]f[5])==1 \ 262 and np.sign(f[4]f[5]) * np.sign(f[3]f[4])==1 \ 263 and np.sign(f[3]f[4]) * np.sign(f[2]f[3])==1: 264 265 # the monotonicity filter corrects the value at (fim1) by 266 # substituting (fim1) with (fmon), see Eq. (27), (28) and (29) 267 fmon = min(3.*g[i2], \ 268 3.*g[i1], \ 269 np.sqrt(max(0,(18./13.*g[i2]  5./13.*f[7]) * 270 (18./13.*g[i1]  5./13.*f[1])))) 271 272 # recomputation of the subgrid interval values while the 273 # interval boundaries (fi) and (fip2) remains unchanged 274 # see Eq. (18) and (19) 275 f[4]=fmon 276 f[6]=3./2.*g[i2]5./12.*fmon1./12.*f[7] 277 f[5]=f[6]+(fmonf[7])/3. 278 f[3]=3./2.*g[i1]5./12.*f[1]1./12.*fmon 279 f[2]=f[3]+(f[1]fmon)/3. 280 281 f.extend([0.,0.,0.]) 282 283 # otherwise the subgrid values are calculated and added to the list 284 else: 285 # temporal save of last value in interpolated list 286 # since it is the left boundary and hence the new (fi) value 287 fi = f[1] 288 289 # the value at the end of the interval (fip1) is prescribed by the 290 # geometric mean, restricted such that nonnegativity is guaranteed 291 # according to Eq. (25) 292 fip1=min( 3.*g[i] , 3.*g[i+1] , np.sqrt(g[i+1]*g[i]) ) 293 294 # the function value at the first subgrid point (fi1) is determined 295 # according to the equal area condition with Eq. (19) 296 fi1=3./2.*g[i]5./12.*fip11./12.*fi 297 298 # the function value at the second subgrid point (fi2) is determined 299 # according Eq. (18) 300 fi2=fi1+1./3.*(fip1fi) 301 302 # apply monotonicity filter for interval before 303 # check if there is "M" or "W" shape 304 if np.sign(f[5]f[6]) * np.sign(f[4]f[5])==1 \ 305 and np.sign(f[4]f[5]) * np.sign(f[3]f[4])==1 \ 306 and np.sign(f[3]f[4]) * np.sign(f[2]f[3])==1: 307 308 # the monotonicity filter corrects the value at (fim1) by 309 # substituting (fim1) with fmon, see Eq. (27), (28) and (29) 310 fmon = min(3.*g[i2], \ 311 3.*g[i1], \ 312 np.sqrt(max(0,(18./13.*g[i2]  5./13.*f[7]) * 313 (18./13.*g[i1]  5./13.*f[1])))) 314 315 # recomputation of the subgrid interval values while the 316 # interval boundaries (fi) and (fip2) remains unchanged 317 # see Eq. (18) and (19) 318 f[4]=fmon 319 f[6]=3./2.*g[i2]5./12.*fmon1./12.*f[7] 320 f[5]=f[6]+(fmonf[7])/3. 321 f[3]=3./2.*g[i1]5./12.*f[1]1./12.*fmon 322 f[2]=f[3]+(f[1]fmon)/3. 323 324 # add next interval of interpolated (sub)grid values 325 f.append(fi1) 326 f.append(fi2) 327 f.append(fip1) 328 329 # separate treatment of the final interval 330 331 # as a requirement: 332 # if there is a zero data value such that g[i]=0, then the whole 333 # interval in f has to be zero to such that f[i+1]=f[i+2]=f[i+3]=0 334 # according to Eq. (6) 335 if g[1]==0.: 336 # apply monotonicity filter for interval before 337 # check if there is "M" or "W" shape 338 if np.sign(f[5]f[6]) * np.sign(f[4]f[5])==1 \ 339 and np.sign(f[4]f[5]) * np.sign(f[3]f[4])==1 \ 340 and np.sign(f[3]f[4]) * np.sign(f[2]f[3])==1: 341 342 # the monotonicity filter corrects the value at (fim1) by 343 # substituting (fim1) with (fmon), see Eq. (27), (28) and (29) 344 fmon = min(3.*g[3], \ 345 3.*g[2], \ 346 np.sqrt(max(0,(18./13.*g[3]  5./13.*f[7]) * 347 (18./13.*g[2]  5./13.*f[1])))) 348 349 # recomputation of the subgrid interval values while the 350 # interval boundaries (fi) and (fip2) remains unchanged 351 # see Eq. (18) and (19) 352 f[4]=fmon 353 f[6]=3./2.*g[3]5./12.*fmon1./12.*f[7] 354 f[5]=f[6]+(fmonf[7])/3. 355 f[3]=3./2.*g[2]5./12.*f[1]1./12.*fmon 356 f[2]=f[3]+(f[1]fmon)/3. 357 358 f.extend([0.,0.,0.]) 359 360 # otherwise the subgrid values are calculated and added to the list 361 # using the persistence hypothesis as boundary condition 362 else: 363 # temporal save of last value in interpolated list 364 # since it is the left boundary and hence the new (fi) value 365 fi = f[1] 366 # since last interval in series, last value is also fip1 367 fip1 = g[1] 368 # the function value at the first subgrid point (fi1) is determined 369 # according to the equal area condition with Eq. (19) 370 fi1 = 3./2.*g[1]5./12.*fip11./12.*fi 371 # the function value at the second subgrid point (fi2) is determined 372 # according Eq. (18) 373 fi2 = fi1+dt/3.*(fip1fi) 374 375 # apply monotonicity filter for interval before 376 # check if there is "M" or "W" shape 377 if np.sign(f[5]f[6]) * np.sign(f[4]f[5])==1 \ 378 and np.sign(f[4]f[5]) * np.sign(f[3]f[4])==1 \ 379 and np.sign(f[3]f[4]) * np.sign(f[2]f[3])==1: 380 381 # the monotonicity filter corrects the value at (fim1) by 382 # substituting (fim1) with (fmon), see Eq. (27), (28) and (29) 383 fmon = min(3.*g[3], \ 384 3.*g[2], \ 385 np.sqrt(max(0,(18./13.*g[3]  5./13.*f[7]) * 386 (18./13.*g[2]  5./13.*f[1])))) 387 388 # recomputation of the subgrid interval values while the 389 # interval boundaries (fi) and (fip2) remains unchanged 390 # see Eq. (18) and (19) 391 f[4]=fmon 392 f[6]=3./2.*g[3]5./12.*fmon1./12.*f[7] 393 f[5]=f[6]+(fmonf[7])/3. 394 f[3]=3./2.*g[2]5./12.*f[1]1./12.*fmon 395 f[2]=f[3]+(f[1]fmon)/3. 396 397 # add next interval of interpolated (sub)grid values 398 f.append(fi1) 399 f.append(fi2) 400 f.append(fip1) 401 402 return f
Note: See TracChangeset
for help on using the changeset viewer.