Homec4science
Diffusion ctracker3 (master) Browse

ctracker3/_old/tracker/modules/loop.pyx
16061f733e12master

Last Change
Raw File

loop.pyx
View Options

###########################################################################
# #
# Copyright 2017 Andrea Cimatoribus #
# EPFL ENAC IIE ECOL #
# GR A1 435 (Batiment GR) #
# Station 2 #
# CH-1015 Lausanne #
# Andrea.Cimatoribus@epfl.ch #
# #
# This file is part of ctracker #
# #
# ctracker is free software: you can redistribute it and/or modify it #
# under the terms of the GNU General Public License as published by #
# the Free Software Foundation, either version 3 of the License, or #
# (at your option) any later version. #
# #
# ctracker is distributed in the hope that it will be useful, #
# but WITHOUT ANY WARRANTY; without even the implied warranty #
# of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. #
# See the GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
# along with ctracker. If not, see <http://www.gnu.org/licenses/>. #
# #
###########################################################################
from cython cimport wraparound, boundscheck, cdivision
from libc.math cimport log, exp, abs, fmin
from libc.stdio cimport printf
from libc.stdlib cimport exit, EXIT_FAILURE
@wraparound(False)
@boundscheck(False)
@cdivision(True)
cdef (double, double, double, double) \
cross(int ijk,
int ia,
int ja,
int ka,
double i_futr,
double i_past,
double r0,
double[:, :, :, ::1] flux) nogil:
"""
Compute crossing times
"""
# um: flux entering the cell
# uu: flux exiting the cell
# i_futr: interpolation coefficient for "future" GCM data
# i_past: interpolation coefficient for "past" GCM data
cdef double uu, um, ba, sp, sn, Du, ii
cdef double UNDEF = 1.0e20
um = (i_futr * flux[1, ka, ja, ia] +
i_past * flux[0, ka, ja, ia])
if (ijk == 1):
ii = ia
uu = (i_futr * flux[1, ka, ja, ia + 1] +
i_past * flux[0, ka, ja, ia + 1])
elif (ijk == 2):
ii = ja
uu = (i_futr * flux[1, ka, ja + 1, ia] +
i_past * flux[0, ka, ja + 1, ia])
elif (ijk == 3):
ii = ka
uu = (i_futr * flux[1, ka + 1, ja, ia] +
i_past * flux[0, ka + 1, ja, ia])
# east, north or upward crossing
if (uu > 0.0) and (r0 != (ii + 1)):
if um != uu:
Du = (uu - um)
ba = (r0 - ii) * Du + um
if ba > 0.0:
# eq. 7.9 & 7.6 Doos et al. 2013
sp = (log(uu) - log(ba)) / Du
# It can still be that there is a
# zero flux somewhere between r0 and
# the face
if (sp <= 0.0):
sp = UNDEF
else:
sp = UNDEF
else:
sp = (1.0 + ii - r0) / uu
else:
sp = UNDEF
# west, south or downward crossing
if (um < 0.0) and (r0 != ii):
if um != uu:
# negative velocity, so we go from exit (uu)
# to enter (um)
Du = (um - uu)
# (ii)-----(r0)---(ii+1)
# (um)-----(ba)---(uu)
ba = (1 + ii - r0) * Du + uu
if ba < 0.0:
sn = (log(-ba) - log(-um)) / Du
# It can still be that there is a
# zero flux somewhere between r0 and
# the face
if (sn <= 0.0):
sn = UNDEF
else:
sn = UNDEF
else:
sn = (ii - r0) / uu
else:
sn = UNDEF
return uu, um, sp, sn
@cdivision(True)
cdef double newpos(int ii,
double r0,
double ds,
double uu,
double um) nogil:
"""
Compute new position along ijk direction
"""
# um: flux entering the cell
# uu: flux exiting the cell
cdef double ba, Du, Dum
if um != uu:
Du = uu - um
mDu = um / Du
# eq. 7.8 Doos et al 2013
return (r0 - ii + mDu) * exp(Du * ds) + ii - mDu
else:
return r0 + uu * ds
@wraparound(False)
@boundscheck(False)
def transform(double[:] ii,
double[:] jj,
double[:] kk,
float[:, ::1] X,
float[:, ::1] Y,
float[:, ::1] dX,
float[:, ::1] dY,
float[:, ::1] CS,
float[:, ::1] SN,
float[:, :, ::1] Z,
double[:, ::1] stxyz):
cdef int m
cdef int N = ii.size
for m in range(N):
transform_one(ii[m], jj[m], kk[m],
X, Y, dX, dY, CS, SN, Z,
stxyz[m, :])
@wraparound(False)
@boundscheck(False)
cdef void transform_one(
double i,
double j,
double k,
float[:, ::1] X,
float[:, ::1] Y,
float[:, ::1] dX,
float[:, ::1] dY,
float[:, ::1] CS,
float[:, ::1] SN,
float[:, :, ::1] Z,
double[:] stxyz) nogil:
cdef int ii_int, jj_int, kk_int
cdef double nx, ny, px, py, pxpy, pxny, nxpy, nxny, nz
cdef double cosij, sinij, dxij, dyij
ii_int = int(i)
jj_int = int(j)
kk_int = int(k)
nx = i % 1
ny = j % 1
nz = k % 1
# inside the cell
if (nx >= 1e-9) & (ny >= 1e-9):
cosij = CS[jj_int, ii_int]
sinij = SN[jj_int, ii_int]
dxij = dX[jj_int, ii_int]
dyij = dY[jj_int, ii_int]
# x
stxyz[0] = X[jj_int, ii_int] + \
cosij * dxij * nx - sinij * dyij * ny
# y
stxyz[1] = Y[jj_int, ii_int] + \
sinij * dxij * nx + cosij * dyij * ny
# z
if nz >= 1e-9:
stxyz[2] = (1.0 - nz) * Z[kk_int, jj_int, ii_int] + \
nz * Z[kk_int+1, jj_int, ii_int]
else:
stxyz[2] = Z[kk_int, jj_int, ii_int]
# at the SW corner
elif (nx < 1e-9) & (ny < 1e-9):
# x
stxyz[0] = X[jj_int, ii_int]
# y
stxyz[1] = Y[jj_int, ii_int]
# z
if nz >= 1e-9:
stxyz[2] = (1.0 - nz) * Z[kk_int, jj_int, ii_int] + \
nz * Z[kk_int+1, jj_int, ii_int]
else:
stxyz[2] = Z[kk_int, jj_int, ii_int]
# on W face
elif nx < 1e-9:
cosij = CS[jj_int, ii_int]
sinij = SN[jj_int, ii_int]
dyij = dY[jj_int, ii_int]
# x
stxyz[0] = X[jj_int, ii_int] - sinij * dyij * ny
# y
stxyz[1] = Y[jj_int, ii_int] + cosij * dyij * ny
# z
if nz >= 1e-9:
stxyz[2] = (1.0 - nz) * Z[kk_int, jj_int, ii_int] + \
nz * Z[kk_int+1, jj_int, ii_int]
else:
stxyz[2] = Z[kk_int, jj_int, ii_int]
# on S face
else:
cosij = CS[jj_int, ii_int]
sinij = SN[jj_int, ii_int]
dxij = dX[jj_int, ii_int]
# x
stxyz[0] = X[jj_int, ii_int] + cosij * dxij * nx
# y
stxyz[1] = Y[jj_int, ii_int] + sinij * dxij * nx
# z
if nz >= 1e-9:
stxyz[2] = (1.0 - nz) * Z[kk_int, jj_int, ii_int] + \
nz * Z[kk_int+1, jj_int, ii_int]
else:
stxyz[2] = Z[kk_int, jj_int, ii_int]
@wraparound(False)
@boundscheck(False)
@cdivision(True)
cdef void loop_particle(
double[:, ::1] xyz,
short[:, ::1] ijk,
double[:, :, ::1] out_tijk,
double[:, :, ::1] out_xyz,
short[:] out_code,
int imax,
int jmax,
int kmax,
double out_dt,
float[:, :, ::1] Adsmax,
float[:, :, ::1] Adxyz,
double[:, :, :, ::1] uflux,
double[:, :, :, ::1] vflux,
double[:, :, :, ::1] wflux,
double dtmax,
double dstep,
int nend,
double[:] iendw,
double[:] iende,
double[:] jends,
double[:] jendn,
int outfreq,
int nvals,
int subiters,
float[:, ::1] xG,
float[:, ::1] yG,
float[:, ::1] dX,
float[:, ::1] dY,
float[:, ::1] CS,
float[:, ::1] SN,
float[:, :, ::1] Z) nogil:
cdef double x0, y0, z0, x1, y1, z1, ts, tt, dt
cdef double i_futr, i_past, in_futr, in_past
cdef double dsmax, dxyz
cdef double uu1, um1, vu1, vm1, wu1, wm1, uu
cdef double ds, dse, dsw, dsn, dss, dsu, dsd
cdef int ia, ja, ka, ib, jb, kb, rec_out, kk
cdef int scrivi, kill_area
cdef long ntrac, ntact = xyz.shape[0]
for ntrac in range(ntact):
# initialize stuff
ts = 0.0 # normalised time between two GCM time steps
tt = 0.0 # time in seconds between two GCM time steps
# output init
for kk in range(nvals):
out_tijk[ntrac, kk, 0] = -1e20
rec_out = 0
# take coord
# the coordinate system is as follows:
# cell i
# |---------------------|
# i i <= x < i+1 i+1
x1 = xyz[ntrac, 0]
y1 = xyz[ntrac, 1]
z1 = xyz[ntrac, 2]
ib = ijk[ntrac, 0]
jb = ijk[ntrac, 1]
kb = ijk[ntrac, 2]
scrivi = 0
while True:
# if we reached the new GCM output
if tt == out_dt:
xyz[ntrac, 0] = x1
xyz[ntrac, 1] = y1
xyz[ntrac, 2] = z1
ijk[ntrac, 0] = ib
ijk[ntrac, 1] = jb
ijk[ntrac, 2] = kb
break
# time interpolation constants
i_futr = ts % 1.0
i_past = 1.0 - i_futr
x0 = x1
y0 = y1
z0 = z1
ia = ib
ja = jb
ka = kb
dsmax = Adsmax[kb, jb, ib]
dxyz = Adxyz[kb, jb, ib]
# make sure we do not start from land
#TODO do we still need this?
if dxyz == 0.0:
printf("\nx0: %f\n", x0)
printf("y0: %f\n", y0)
printf("z0: %f\n", z0)
printf("ia: %d\n", ia)
printf("ja: %d\n", ja)
printf("ka: %d\n", ka)
printf("x1: %f\n", x1)
printf("y1: %f\n", y1)
printf("z1: %f\n", z1)
printf("ib: %d\n", ib)
printf("jb: %d\n", jb)
printf("kb: %d\n", kb)
printf("ds: %f\n", ds)
printf("dsmax: %f\n", dsmax)
printf("dxyz: %f\n", dxyz)
printf("dsn: %f\n", dsn)
printf("dss: %f\n", dss)
printf("dse: %f\n", dse)
printf("dsw: %f\n", dsw)
printf("dsu: %f\n", dsu)
printf("dsd: %f\n", dsd)
printf("Particle on land.")
exit(EXIT_FAILURE)
#
# calculate the 3 crossing times over the box
# choose the shortest time and calculate the
# new positions
#
# solving the differential equations
# note:
# space variables (x,...) are dimensionless
# time variables (ds,...) are in seconds/m^3
#
uu1, um1, dse, dsw = \
cross(1, ia, ja, ka, i_futr, i_past, x0, uflux)
vu1, vm1, dsn, dss = \
cross(2, ia, ja, ka, i_futr, i_past, y0, vflux)
wu1, wm1, dsu, dsd = \
cross(3, ia, ja, ka, i_futr, i_past, z0, wflux)
ds = fmin(fmin(fmin(fmin(fmin(fmin(dse, dsw), dsn), dss),
dsd), dsu), dsmax)
if (ds == 0.0) or (ds == 1e20):
printf("\nx0: %f\n", x0)
printf("y0: %f\n", y0)
printf("z0: %f\n", z0)
printf("ia: %d\n", ia)
printf("ja: %d\n", ja)
printf("ka: %d\n", ka)
printf("x1: %f\n", x1)
printf("y1: %f\n", y1)
printf("z1: %f\n", z1)
printf("ib: %d\n", ib)
printf("jb: %d\n", jb)
printf("kb: %d\n", kb)
printf("ds: %f\n", ds)
printf("dsmax: %f\n", dsmax)
printf("dxyz: %f\n", dxyz)
printf("dsn: %f\n", dsn)
printf("dss: %f\n", dss)
printf("dse: %f\n", dse)
printf("dsw: %f\n", dsw)
printf("dsu: %f\n", dsu)
printf("dsd: %f\n", dsd)
printf("Cannot integrate track for unknown reasons\n")
exit(EXIT_FAILURE)
# now update the time
dt = ds * dxyz
# if time step makes the integration time
# very close to the GCM output
if (tt + dt) >= out_dt:
dt = out_dt - tt
tt = out_dt
ts = 1.0
ds = dt / dxyz
end_loop = True
else:
end_loop = False
tt += dt
if dt == dtmax:
ts += dstep
else:
ts += dt / out_dt
# compute new time interpolation constant
in_futr = ts % 1.0
in_past = 1.0 - in_futr
# now we actually compute the new position
# eastward grid-cell exit
if ds == dse:
scrivi = 1
# if at the "new time" the flow at the cell
# face is positive the particle will exit
uu = (in_futr * uflux[1, ka, ja, ia + 1] +
in_past * uflux[0, ka, ja, ia + 1])
if uu > 0.0:
ib = ia + 1
x1 = float(ia + 1)
y1 = newpos(ja, y0, ds, vu1, vm1)
z1 = newpos(ka, z0, ds, wu1, wm1)
# deal with corner cases (including rounding error)
#TODO: corner cases are extremely rare, but
# may happen. Is there a better way to deal with them?
if (y1 - ja) >= 1.0: # y
uu = (in_futr * vflux[1, ka, ja + 1, ia] +
in_past * vflux[0, ka, ja + 1, ia])
if uu > 0.0:
jb = ja + 1
y1 = float(ja + 1)
elif (y1 - ja) <= 0.0:
uu = (in_futr * vflux[1, ka, ja, ia] +
in_past * vflux[0, ka, ja, ia])
if uu < 0.0:
jb = ja - 1
y1 = float(ja)
if (z1 - ka) >= 1.0: # z
uu = (in_futr * wflux[1, ka + 1, ja, ia] +
in_past * wflux[0, ka + 1, ja, ia])
if uu > 0.0:
kb = ka + 1
z1 = float(ka + 1)
elif (z1 - ka) <= 0.0:
uu = (in_futr * wflux[1, ka, ja, ia] +
in_past * wflux[0, ka, ja, ia])
if uu < 0.0:
kb = ka - 1
z1 = float(ka)
# westward grid-cell exit
elif ds == dsw:
scrivi = 1
uu = (in_futr * uflux[1, ka, ja, ia] +
in_past * uflux[0, ka, ja, ia])
if uu < 0.0:
ib = ia - 1
x1 = float(ia)
y1 = newpos(ja, y0, ds, vu1, vm1)
z1 = newpos(ka, z0, ds, wu1, wm1)
# deal with corner cases (including rounding error)
if (y1 - ja) >= 1.0: # y
uu = (in_futr * vflux[1, ka, ja + 1, ia] +
in_past * vflux[0, ka, ja + 1, ia])
if uu > 0.0:
jb = ja + 1
y1 = float(ja + 1)
elif (y1 - ja) <= 0.0:
uu = (in_futr * vflux[1, ka, ja, ia] +
in_past * vflux[0, ka, ja, ia])
if uu < 0.0:
jb = ja - 1
y1 = float(ja)
if (z1 - ka) >= 1.0: # z
uu = (in_futr * wflux[1, ka + 1, ja, ia] +
in_past * wflux[0, ka + 1, ja, ia])
if uu > 0.0:
kb = ka + 1
z1 = float(ka + 1)
elif (z1 - ka) <= 0.0:
uu = (in_futr * wflux[1, ka, ja, ia] +
in_past * wflux[0, ka, ja, ia])
if uu < 0.0:
kb = ka - 1
z1 = float(ka)
# northward grid-cell exit
elif ds == dsn:
scrivi = 1
uu = (in_futr * vflux[1, ka, ja + 1, ia] +
in_past * vflux[0, ka, ja + 1, ia])
if uu > 0.0:
jb = ja + 1
x1 = newpos(ia, x0, ds, uu1, um1)
y1 = float(ja + 1)
z1 = newpos(ka, z0, ds, wu1, wm1)
# deal with corner cases (including rounding error)
if (x1 - ia) >= 1.0: # x
uu = (in_futr * uflux[1, ka, ja, ia + 1] +
in_past * uflux[0, ka, ja, ia + 1])
if uu > 0.0:
ib = ia + 1
x1 = float(ia + 1)
elif (x1 - ia) <= 0.0:
uu = (in_futr * uflux[1, ka, ja, ia] +
in_past * uflux[0, ka, ja, ia])
if uu < 0.0:
ib = ia - 1
x1 = float(ia)
if (z1 - ka) >= 1.0: # z
uu = (in_futr * wflux[1, ka + 1, ja, ia] +
in_past * wflux[0, ka + 1, ja, ia])
if uu > 0.0:
kb = ka + 1
z1 = float(ka + 1)
elif (z1 - ka) <= 0.0:
uu = (in_futr * wflux[1, ka, ja, ia] +
in_past * wflux[0, ka, ja, ia])
if uu < 0.0:
kb = ka - 1
z1 = float(ka)
# southward grid-cell exit
elif ds == dss:
scrivi = 1
uu = (in_futr * vflux[1, ka, ja, ia] +
in_past * vflux[0, ka, ja, ia])
if uu < 0.0:
jb = ja - 1
x1 = newpos(ia, x0, ds, uu1, um1)
y1 = float(ja)
z1 = newpos(ka, z0, ds, wu1, wm1)
# deal with corner cases (including rounding error)
if (x1 - ia) >= 1.0: # x
uu = (in_futr * uflux[1, ka, ja, ia + 1] +
in_past * uflux[0, ka, ja, ia + 1])
if uu > 0.0:
ib = ia + 1
x1 = float(ia + 1)
elif (x1 - ia) <= 0.0:
uu = (in_futr * uflux[1, ka, ja, ia] +
in_past * uflux[0, ka, ja, ia])
if uu < 0.0:
ib = ia - 1
x1 = float(ia)
if (z1 - ka) >= 1.0: # z
uu = (in_futr * wflux[1, ka + 1, ja, ia] +
in_past * wflux[0, ka + 1, ja, ia])
if uu > 0.0:
kb = ka + 1
z1 = float(ka + 1)
elif (z1 - ka) <= 0.0:
uu = (in_futr * wflux[1, ka, ja, ia] +
in_past * wflux[0, ka, ja, ia])
if uu < 0.0:
kb = ka - 1
z1 = float(ka)
# upward grid-cell exit
elif ds == dsu:
scrivi = 1
uu = (in_futr * wflux[1, ka + 1, ja, ia] +
in_past * wflux[0, ka + 1, ja, ia])
if uu > 0.0:
kb = ka + 1
x1 = newpos(ia, x0, ds, uu1, um1)
y1 = newpos(ja, y0, ds, vu1, vm1)
z1 = float(ka + 1)
# deal with corner cases (including rounding error)
if (x1 - ia) >= 1.0: # x
uu = (in_futr * uflux[1, ka, ja, ia + 1] +
in_past * uflux[0, ka, ja, ia + 1])
if uu > 0.0:
ib = ia + 1
x1 = float(ia + 1)
elif (x1 - ia) <= 0.0:
uu = (in_futr * uflux[1, ka, ja, ia] +
in_past * uflux[0, ka, ja, ia])
if uu < 0.0:
ib = ia - 1
x1 = float(ia)
if (y1 - ja) >= 1.0: # y
uu = (in_futr * vflux[1, ka, ja + 1, ia] +
in_past * vflux[0, ka, ja + 1, ia])
if uu > 0.0:
jb = ja + 1
y1 = float(ja + 1)
elif (y1 - ja) <= 0.0:
uu = (in_futr * vflux[1, ka, ja, ia] +
in_past * vflux[0, ka, ja, ia])
if uu < 0.0:
jb = ja - 1
y1 = float(ja)
# downward grid-cell exit
elif ds == dsd:
scrivi = 1
uu = (in_futr * wflux[1, ka, ja, ia] +
in_past * wflux[0, ka, ja, ia])
if uu < 0.0:
kb = ka - 1
x1 = newpos(ia, x0, ds, uu1, um1)
y1 = newpos(ja, y0, ds, vu1, vm1)
z1 = float(ka)
# deal with corner cases (including rounding error)
if (x1 - ia) >= 1.0: # x
uu = (in_futr * uflux[1, ka, ja, ia + 1] +
in_past * uflux[0, ka, ja, ia + 1])
if uu > 0.0:
ib = ia + 1
x1 = float(ia + 1)
elif (x1 - ia) <= 0.0:
uu = (in_futr * uflux[1, ka, ja, ia] +
in_past * uflux[0, ka, ja, ia])
if uu < 0.0:
ib = ia - 1
x1 = float(ia)
if (y1 - ja) >= 1.0: # y
uu = (in_futr * vflux[1, ka, ja + 1, ia] +
in_past * vflux[0, ka, ja + 1, ia])
if uu > 0.0:
jb = ja + 1
y1 = float(ja + 1)
elif (y1 - ja) <= 0.0:
uu = (in_futr * vflux[1, ka, ja, ia] +
in_past * vflux[0, ka, ja, ia])
if uu < 0.0:
jb = ja - 1
y1 = float(ja)
elif end_loop or (ds == dsmax):
scrivi = 1
x1 = newpos(ia, x0, ds, uu1, um1)
y1 = newpos(ja, y0, ds, vu1, vm1)
z1 = newpos(ka, z0, ds, wu1, wm1)
#TODO for 2D advection we will need to include
# the check for convergence/divergence zones
else:
printf("\nUnrecognised ds\n")
printf("ds: %f\n", ds)
printf("dsmax: %f\n", dsmax)
printf("dxyz: %f\n", dxyz)
printf("dsn: %f\n", dsn)
printf("dss: %f\n", dss)
printf("dse: %f\n", dse)
printf("dsw: %f\n", dsw)
printf("dsu: %f\n", dsu)
printf("dsd: %f\n", dsd)
exit(EXIT_FAILURE)
# check if particle entered a kill area
kill_area = 0
for kk in range(nend):
if (iendw[kk] <= x1) and \
(x1 <= iende[kk]) and \
(jends[kk] <= y1) and \
(y1 <= jendn[kk]):
# we store the region where it was killed
kill_area = 10 + kk
break # break the for loop over kill areas
if kill_area > 0:
# write to buffer
if rec_out >= nvals:
printf("Kill region: buffer finished, make it larger!")
exit(EXIT_FAILURE)
out_tijk[ntrac, rec_out, 0] = tt
out_tijk[ntrac, rec_out, 1] = x1
out_tijk[ntrac, rec_out, 2] = y1
out_tijk[ntrac, rec_out, 3] = z1
out_code[ntrac] = kill_area
transform_one(x1, y1, z1,
xG, yG, dX, dY, CS, SN, Z,
out_xyz[ntrac, rec_out, :])
rec_out += 1
break # break particle loop
# check if the particle is still
# inside the domain
if (ib > imax) or (ib < 0) or \
(jb > jmax) or (jb < 0):
# write to buffer
if rec_out >= nvals:
printf("Exit domain: buffer finished, make it larger!")
exit(EXIT_FAILURE)
out_tijk[ntrac, rec_out, 0] = tt
out_tijk[ntrac, rec_out, 1] = x1
out_tijk[ntrac, rec_out, 2] = y1
out_tijk[ntrac, rec_out, 3] = z1
out_code[ntrac] = 1
transform_one(x1, y1, z1,
xG, yG, dX, dY, CS, SN, Z,
out_xyz[ntrac, rec_out, :])
rec_out += 1
break
# check if the particle exited through the top
if (kb >= kmax):
# write to buffer
#TODO there should be a better way to deal with the
# top boundary condition
if rec_out >= nvals:
printf("Airborne: buffer finished, make it larger!")
exit(EXIT_FAILURE)
out_tijk[ntrac, rec_out, 0] = tt
out_tijk[ntrac, rec_out, 1] = x1
out_tijk[ntrac, rec_out, 2] = y1
out_tijk[ntrac, rec_out, 3] = z1
out_code[ntrac] = 2
transform_one(x1, y1, z1,
xG, yG, dX, dY, CS, SN, Z,
out_xyz[ntrac, rec_out, :])
rec_out += 1
break
# At the end of the cycle we can write the
# position of the particle in buffer
# gcmstep => outfreq == 1
# always => outfreq == 2
# cross => outfreq == 3
if ((tt == out_dt) and (outfreq == 1)) \
or ((outfreq == 3) and (scrivi == 1)) \
or ((outfreq == 2) and (ts > 0.0)):
# write to buffer
if rec_out >= nvals:
printf("Normal write: buffer finished, make it larger!")
exit(EXIT_FAILURE)
out_tijk[ntrac, rec_out, 0] = tt
out_tijk[ntrac, rec_out, 1] = x1
out_tijk[ntrac, rec_out, 2] = y1
out_tijk[ntrac, rec_out, 3] = z1
transform_one(x1, y1, z1,
xG, yG, dX, dY, CS, SN, Z,
out_xyz[ntrac, rec_out, :])
rec_out += 1
def loop_nogil(
double[:, ::1] xyz,
short[:, ::1] ijk,
double[:, :, ::1] out_tijk,
double[:, :, ::1] out_xyz,
short[:] out_code,
int imax,
int jmax,
int kmax,
double out_dt,
float[:, :, ::1] Adsmax,
float[:, :, ::1] Adxyz,
double[:, :, :, ::1] uflux,
double[:, :, :, ::1] vflux,
double[:, :, :, ::1] wflux,
double dtmax,
double dstep,
int nend,
double[:] iendw,
double[:] iende,
double[:] jends,
double[:] jendn,
int outfreq,
int nvals,
int subiters,
float[:, ::1] xG,
float[:, ::1] yG,
float[:, ::1] dX,
float[:, ::1] dY,
float[:, ::1] CS,
float[:, ::1] SN,
float[:, :, ::1] Z):
"""
Kill the GIL
"""
with nogil:
loop_particle(
xyz,
ijk,
out_tijk,
out_xyz,
out_code,
imax,
jmax,
kmax,
out_dt,
Adsmax,
Adxyz,
uflux,
vflux,
wflux,
dtmax,
dstep,
nend,
iendw,
iende,
jends,
jendn,
outfreq,
nvals,
subiters,
xG,
yG,
dX,
dY,
CS,
SN,
Z)
c4science · Help