Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76075631
controller.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Aug 6, 02:20
Size
22 KB
Mime Type
text/x-python
Expires
Thu, Aug 8, 02:20 (2 d)
Engine
blob
Format
Raw Data
Handle
19657474
Attached To
rCTRACKER ctracker3
controller.py
View Options
###########################################################################
# #
# Copyright 2017 Andrea Cimatoribus #
# EPFL ENAC IIE ECOL #
# GR A1 435 (Batiment GR) #
# Station 2 #
# CH-1015 Lausanne #
# Andrea.Cimatoribus@epfl.ch #
# #
# Alexandre Serex #
# alexandre.serex@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/>. #
# #
###########################################################################
import
sys
from
os.path
import
join
import
importlib
import
configparser
import
json
import
numpy
as
np
from
itertools
import
product
from
threading
import
Thread
from
queue
import
Queue
from
netCDF4
import
Dataset
from
xmitgcm
import
open_mdsdataset
as
mitgcmds
# from tracker.mods import transform, transform_cartesian
# from tracker.mods import loop_nogil, loopC_nogil
from
core
import
parameterContainer
,
parameter
from
mods.grid
import
Grid
from
mods.particleContainer
import
ParticleContainer
from
tools.containers
import
GridData
,
ParticleContainerData
class
Controller
:
def
__init__
(
self
,
config_file_name
):
self
.
config_file_name
=
config_file_name
self
.
config
=
{}
self
.
config_methods
=
{
"json"
:
self
.
config_json
,
"ini"
:
self
.
config_ini
}
try
:
file_extension
=
self
.
config_file_name
.
split
(
'.'
)[
-
1
]
self
.
config_methods
.
get
(
file_extension
,
self
.
config_invalid
)()
except
:
print
(
"
\n
Unsupported configuration file extension."
)
raise
print
(
"
\n
Configuration loaded successfully."
)
self
.
_grid
=
Grid
(
"grid"
,
GridData
(
0
))
self
.
_part_container
=
ParticleContainer
(
"part_container"
,
ParticleContainerData
(
0
))
try
:
self
.
init_mitgcm
()
self
.
init_grid
()
self
.
init_sim
()
self
.
init_parallelism
()
except
:
print
(
"
\n
Unsupported configuration file content."
)
raise
try
:
self
.
check_config
()
except
:
print
(
"
\n
Invalid configuration setup."
)
raise
print
(
"
\n
Initialization successful."
)
def
check_config
(
self
):
'''
Here lies all the ways the configuration might go wrong,
therefore all the exceptions that can be raised through initialization.
'''
if
self
.
gcm_geometry
not
in
(
"curvilinear"
,
"cartesian"
):
raise
ValueError
(
"Unrecognised MITgcm geometry"
)
if
self
.
ff
not
in
[
1.0
,
-
1.0
]:
raise
ValueError
(
"ff controls the time direction, "
"it can only be +1 or -1"
)
if
self
.
ff
==
-
1.0
:
raise
NotImplementedError
(
"Backward integration not yet supported."
)
if
self
.
outfreq
==
"custom"
:
raise
NotImplementedError
(
"Outputing data at custom timestamps not yet supported."
)
if
self
.
outfreq
not
in
(
"gcmstep"
,
"always"
,
"cross"
,
"custom"
):
raise
ValueError
(
"Outfreq not recognised"
)
if
np
.
any
(
self
.
dxyz
<
0.0
):
raise
ValueError
(
"Cells with negative volume."
)
if
np
.
any
((
self
.
dxyz
==
0.0
)
&
(
self
.
grid
.
hFacC
[::
-
1
,
:,
:]
!=
0.0
)):
raise
ValueError
(
"Zero cell volumes."
)
if
(
self
.
nend
!=
self
.
iendw
.
size
)
or
\
(
self
.
nend
!=
self
.
jendn
.
size
)
or
\
(
self
.
nend
!=
self
.
jends
.
size
):
raise
ValueError
(
"Wrong kill area definition"
)
if
(
self
.
out_dt
%
self
.
gcm_dt
)
!=
0
:
raise
ValueError
(
"The GCM output interval must be "
"a multiple of the model time step."
)
if
self
.
config
[
"Sim"
][
"inds_seed"
]:
if
np
.
any
(
self
.
x_seed
>
self
.
imax
)
or
\
np
.
any
(
self
.
x_seed
<
0
)
or
\
np
.
any
(
self
.
y_seed
>
self
.
jmax
)
or
\
np
.
any
(
self
.
y_seed
<
0
)
or
\
np
.
any
(
self
.
z_seed
>
self
.
kmax
)
or
\
np
.
any
(
self
.
z_seed
<
0
):
raise
ValueError
(
"Seeding with indexes beyond grid limits"
)
if
self
.
iters
.
size
<
2
:
raise
ValueError
(
"To start the computation, the simulation "
"must span in time at least two GCM outputs"
)
if
(
self
.
seed_interval
%
self
.
out_dt
)
!=
0
:
raise
ValueError
(
"The seeding interval must be "
"a multiple of the GCM step."
)
def
init_mitgcm
(
self
):
'''
Initialize the gcm input / output
'''
self
.
gcm_geometry
=
self
.
config
[
"MITgcm"
][
"gcm_geometry"
]
if
not
self
.
config
[
"MITgcm"
][
"u_root"
]
.
endswith
(
'.'
):
self
.
config
[
"MITgcm"
][
"u_root"
]
+=
'.'
if
not
self
.
config
[
"MITgcm"
][
"v_root"
]
.
endswith
(
'.'
):
self
.
config
[
"MITgcm"
][
"v_root"
]
+=
'.'
if
not
self
.
config
[
"MITgcm"
][
"v_root"
]
.
endswith
(
'.'
):
self
.
config
[
"MITgcm"
][
"v_root"
]
+=
'.'
self
.
u_root
=
join
(
self
.
gcm_geometry
,
self
.
config
[
"MITgcm"
][
"u_root"
])
self
.
v_root
=
join
(
self
.
gcm_geometry
,
self
.
config
[
"MITgcm"
][
"v_root"
])
self
.
w_root
=
join
(
self
.
gcm_geometry
,
self
.
config
[
"MITgcm"
][
"w_root"
])
self
.
out_prec
=
self
.
config
[
"MITgcm"
][
"out_prec"
]
self
.
gcm_dt
=
self
.
config
[
"MITgcm"
][
"gcm_dt"
]
self
.
out_dt
=
self
.
config
[
"MITgcm"
][
"out_dt"
]
self
.
gcm_start
=
self
.
config
[
"MITgcm"
][
"gcm_start"
]
def
init_grid
(
self
):
'''
Initializes the grid data
'''
# open data directory to load grid data
self
.
gcm_directory
=
self
.
config
[
"MITgcm"
][
"gcm_directory"
]
self
.
grid
=
mitgcmds
(
self
.
gcm_directory
,
read_grid
=
True
,
iters
=
[],
prefix
=
[
self
.
config
[
"MITgcm"
][
"u_root"
],
self
.
config
[
"MITgcm"
][
"v_root"
],
self
.
config
[
"MITgcm"
][
"w_root"
]],
swap_dims
=
False
,
geometry
=
self
.
gcm_geometry
,
ref_date
=
self
.
gcm_start
,
delta_t
=
self
.
gcm_dt
)
# load metrics to compute the conversion from fractional index
# to physical coordinate
self
.
xG
=
self
.
grid
.
XG
.
to_masked_array
()
.
filled
(
0
)
.
astype
(
"float32"
)
self
.
yG
=
self
.
grid
.
YG
.
to_masked_array
()
.
filled
(
0
)
.
astype
(
"float32"
)
self
.
dX
=
self
.
grid
.
dxG
.
to_masked_array
()
.
filled
(
0
)
.
astype
(
"float32"
)
self
.
dY
=
self
.
grid
.
dyG
.
to_masked_array
()
.
filled
(
0
)
.
astype
(
"float32"
)
self
.
dxdy
=
self
.
grid
.
rA
.
to_masked_array
()
.
filled
(
0
)
.
astype
(
"float32"
)
# we invert to have at the first position the bottom
self
.
dzt
=
np
.
ascontiguousarray
((
self
.
grid
.
drF
*
self
.
grid
.
hFacC
)
.
to_masked_array
()
.
filled
(
0
)[::
-
1
,
:,
:])
self
.
grid_shape
=
self
.
dzt
.
shape
self
.
kmax
,
self
.
jmax
,
self
.
imax
=
self
.
grid_shape
zG
=
np
.
zeros
((
self
.
kmax
+
1
,
self
.
jmax
,
self
.
imax
))
zG
[
1
:,
...
]
=
np
.
cumsum
(
self
.
dzt
[::
-
1
],
axis
=
0
)
# tracmass has opposite Z order
zG
=
zG
[::
-
1
,
...
]
self
.
Z
=
np
.
ascontiguousarray
(
zG
)
.
astype
(
"float32"
)
self
.
dxyz
=
self
.
dxdy
*
self
.
dzt
self
.
dsmax
=
self
.
dtmax
/
self
.
dxyz
self
.
dzu
=
np
.
ascontiguousarray
((
self
.
grid
.
drF
*
self
.
grid
.
hFacW
*
self
.
grid
.
dyG
)
.
to_masked_array
()
.
filled
(
0
)[::
-
1
,
:,
:])
self
.
dzv
=
np
.
ascontiguousarray
((
self
.
grid
.
drF
*
self
.
grid
.
hFacS
*
self
.
grid
.
dxG
)
.
to_masked_array
()
.
filled
(
0
)[::
-
1
,
:,
:])
self
.
kmtb
=
(
self
.
grid
.
hFacC
.
sum
(
dim
=
"k"
))
.
to_masked_array
()
.
filled
(
0
)
self
.
kmt
=
np
.
ceil
(
self
.
kmtb
)
.
astype
(
"int8"
)
self
.
CS
=
1
self
.
SN
=
0
if
self
.
gcm_geometry
==
"curvilinear"
:
self
.
CS
=
self
.
grid
.
CS
\
.
to_masked_array
()
.
filled
(
0
)
.
astype
(
"float32"
)
self
.
SN
=
self
.
grid
.
SN
\
.
to_masked_array
()
.
filled
(
0
)
.
astype
(
"float32"
)
xGp1
,
yGp1
=
self
.
_get_geometry
()
self
.
grid
.
attrs
[
"MITgcm_dir"
]
=
self
.
gcm_directory
# we also store the full mesh with edges coordinates
self
.
grid
.
coords
[
"i_p1"
]
=
(
"i_p1"
,
np
.
arange
(
xGp1
.
shape
[
1
]))
self
.
grid
.
coords
[
"j_p1"
]
=
(
"j_p1"
,
np
.
arange
(
xGp1
.
shape
[
0
]))
self
.
grid
.
coords
[
"XG_p1"
]
=
((
"j_p1"
,
"i_p1"
),
xGp1
)
self
.
grid
.
coords
[
"YG_p1"
]
=
((
"j_p1"
,
"i_p1"
),
yGp1
)
def
init_sim
(
self
):
'''
Initializes the seeding frequency and position and seeding steps for the simulation
'''
outfile_name
=
self
.
config
[
"Sim"
][
"outfile"
]
.
rstrip
(
".nc"
)
self
.
inoutfile
=
outfile_name
+
"_inout.nc"
self
.
runfile
=
outfile_name
+
"_run.nc"
self
.
outfreq
=
self
.
config
[
"Sim"
][
"outfreq"
]
self
.
out_gcmstep
=
self
.
config
[
"Sim"
][
"out_gcmstep"
]
self
.
subiters
=
int
(
self
.
config
[
"Sim"
][
"subiters"
])
self
.
dstep
=
1.0
/
self
.
subiters
self
.
dtmax
=
self
.
out_dt
*
self
.
dstep
self
.
ff
=
float
(
self
.
config
[
"Sim"
][
"ff"
])
# list of seed points coordinates
# we seed along a transect off the Rhone outflow
zs
=
np
.
arange
(
self
.
config
[
"Sim"
][
"z0"
],
self
.
config
[
"Sim"
][
"z1"
],
self
.
config
[
"Sim"
][
"dz"
])
xs
=
np
.
linspace
(
self
.
config
[
"Sim"
][
"p0x"
],
self
.
config
[
"Sim"
][
"p1x"
],
self
.
config
[
"Sim"
][
"npts"
])
ys
=
np
.
linspace
(
self
.
config
[
"Sim"
][
"p0y"
],
self
.
config
[
"Sim"
][
"p1y"
],
self
.
config
[
"Sim"
][
"npts"
])
X
,
Z
=
np
.
meshgrid
(
xs
,
zs
)
Y
,
Z
=
np
.
meshgrid
(
ys
,
zs
)
self
.
x_seed
=
X
.
ravel
()
self
.
y_seed
=
Y
.
ravel
()
self
.
z_seed
=
Z
.
ravel
()
self
.
gcm_start
=
np
.
datetime64
(
self
.
gcm_start
)
print
(
"
\n
Reference date of GCM:"
,
self
.
gcm_start
)
self
.
start
=
np
.
datetime64
(
self
.
config
[
"Sim"
][
"start"
])
print
(
"
\n
Begin particle tracking simulation at:"
,
self
.
start
)
self
.
end
=
np
.
datetime64
(
self
.
config
[
"Sim"
][
"end"
])
print
(
"
\n
End particle tracking simulation at:"
,
self
.
end
)
self
.
seed_start
=
np
.
datetime64
(
self
.
config
[
"Sim"
][
"seed_start"
])
print
(
"
\n
Start seeding particles at:"
,
self
.
seed_start
)
self
.
seed_end
=
np
.
datetime64
(
self
.
config
[
"Sim"
][
"seed_end"
])
print
(
"
\n
Stop seeding particles at:"
,
self
.
seed_end
)
self
.
is2D
=
self
.
config
[
"Sim"
][
"is2D"
];
print
(
"
\n
Simulation in 2D : "
,
self
.
is2D
)
#
# prepare the simulation
#
print
(
"
\n
Identify seeding points"
)
self
.
inds_seed
=
self
.
config
[
"Sim"
][
"inds_seed"
]
self
.
ijk_indseed
,
self
.
ijk_seed
=
\
self
.
_ijkseed
(
self
.
x_seed
,
self
.
y_seed
,
self
.
z_seed
,
self
.
inds_seed
)
self
.
xyz_seed
=
np
.
zeros
(
self
.
ijk_seed
.
shape
)
if
self
.
gcm_geometry
==
"curvilinear"
:
print
(
"
\n
Curvilinear simulation !"
)
'''transform(self.ijk_seed[:, 0],
self.ijk_seed[:, 1],
self.ijk_seed[:, 2],
self.xG, self.yG, self.dX, self.dY,
self.CS, self.SN, self.Z,
self.xyz_seed)'''
else
:
print
(
"
\n
Cartesian simulation !"
)
'''transform_cartesian(self.ijk_seed[:, 0],
self.ijk_seed[:, 1],
self.ijk_seed[:, 2],
self.xG, self.yG, self.dX, self.dY,
self.Z, self.xyz_seed)'''
self
.
n_seed
=
self
.
ijk_seed
.
shape
[
0
]
print
(
"
\n
Number of seeding points:
%d
"
%
self
.
n_seed
)
self
.
iende
=
np
.
asarray
(
self
.
config
[
"Sim"
][
"iende"
])
self
.
iendw
=
np
.
asarray
(
self
.
config
[
"Sim"
][
"iendw"
])
self
.
jendn
=
np
.
asarray
(
self
.
config
[
"Sim"
][
"jendn"
])
self
.
jends
=
np
.
asarray
(
self
.
config
[
"Sim"
][
"jends"
])
self
.
nend
=
self
.
iende
.
size
# first GCM iteration
iter_i
=
np
.
floor
(
(
self
.
start
-
self
.
gcm_start
)
.
astype
(
"timedelta64[s]"
)
.
astype
(
float
)
/
self
.
gcm_dt
)
.
astype
(
"int32"
)
# last GCM iteration
iter_e
=
np
.
ceil
(
(
self
.
end
-
self
.
gcm_start
)
.
astype
(
"timedelta64[s]"
)
.
astype
(
float
)
/
self
.
gcm_dt
)
.
astype
(
"int32"
)
self
.
iters
=
np
.
arange
(
iter_i
,
iter_e
+
1
,
int
(
self
.
out_dt
/
self
.
gcm_dt
))
self
.
seed_interval
=
self
.
config
[
"Sim"
][
"seed_interval"
]
print
(
"
\n
Seeding interval:
%d
s"
%
self
.
seed_interval
)
# seeding steps
seed_i
=
np
.
floor
(
(
self
.
seed_start
-
self
.
gcm_start
)
.
astype
(
"timedelta64[s]"
)
.
astype
(
float
)
/
self
.
gcm_dt
)
.
astype
(
"int32"
)
seed_e
=
np
.
ceil
(
(
self
.
seed_end
-
self
.
gcm_start
)
.
astype
(
"timedelta64[s]"
)
.
astype
(
float
)
/
self
.
gcm_dt
)
.
astype
(
"int32"
)
iter_seed
=
np
.
arange
(
seed_i
,
seed_e
+
1
,
int
(
self
.
seed_interval
/
self
.
gcm_dt
))
self
.
iter_seed
=
np
.
intersect1d
(
self
.
iters
,
iter_seed
)
print
(
"
\n
Number of seeding time steps:
%d
"
%
self
.
iter_seed
.
size
)
self
.
ntrackmax
=
self
.
iter_seed
.
size
*
self
.
n_seed
print
(
"
\n
Total number of particles to be released:
%d
"
%
self
.
ntrackmax
)
def
init_parallelism
(
self
):
# configure parallelism
self
.
n_threads
=
int
(
self
.
config
[
"Parallelism"
][
"n_threads"
])
if
self
.
n_threads
<=
1
:
self
.
n_threads
=
1
print
(
"
\n
Running numerics single threaded"
)
else
:
print
(
"
\n
Running numerics on {0} threads max"
.
format
(
self
.
n_threads
))
self
.
n_min_part_per_thread
=
self
.
config
[
"Parallelism"
][
"n_min_part_per_thread"
]
print
(
"
\n
Running with {0} particles per thread min"
.
format
(
self
.
n_min_part_per_thread
))
def
config_invalid
(
self
,
*
args
,
**
kwargs
):
raise
EnvironmentError
(
"Unsupported configuration file extension."
)
def
config_json
(
self
,
*
args
,
**
kwargs
):
with
open
(
self
.
config_file_name
)
as
f
:
self
.
config
=
json
.
load
(
f
)
if
isinstance
(
self
.
config
[
"Sim"
][
"out_custom"
],
list
):
self
.
timestamps_out
=
self
.
config
[
"Sim"
][
"out_custom"
]
else
:
name_out
=
self
.
config
[
"Sim"
][
"out_custom"
]
.
rstrip
(
".py"
)
self
.
timestamps_out
=
getattr
(
importlib
.
__import__
(
name_out
,
fromlist
=
[
name_out
]),
name_out
)
if
isinstance
(
self
.
config
[
"Sim"
][
"in_custom"
],
list
):
self
.
timestamps_in
=
self
.
config
[
"Sim"
][
"in_custom"
]
else
:
name_in
=
self
.
config
[
"Sim"
][
"in_custom"
]
.
rstrip
(
".py"
)
self
.
timestamps_in
=
getattr
(
importlib
.
__import__
(
name_in
,
fromlist
=
[
name_in
]),
name_in
)
def
config_ini
(
self
,
*
args
,
**
kwargs
):
config
=
configparser
.
ConfigParser
()
config
.
read
(
self
.
config_file_name
)
# Store in self.config the values of the config file.
# Example:
# [DEFAULT] hello = world
# will translate to :
# >>> self.config["DEFAULT"]["hello"]
# world
for
section_name
in
config
.
sections
():
section
=
config
[
section_name
]
self
.
config
[
section_name
]
=
{}
for
key
in
section
:
self
.
config
[
section_name
][
key
]
=
config
[
section_name
][
key
]
name_out
=
config
[
"Sim"
][
"out_custom"
]
.
rstrip
(
".py"
)
name_in
=
config
[
"Sim"
][
"in_custom"
]
.
rstrip
(
".py"
)
self
.
timestamps_out
=
getattr
(
importlib
.
__import__
(
name_out
,
fromlist
=
[
name_out
]),
name_out
)
self
.
timestamps_in
=
getattr
(
importlib
.
__import__
(
name_in
,
fromlist
=
[
name_in
]),
name_in
)
if
self
.
config
[
"Sim"
][
"outfreq"
]
==
"custom"
:
raise
NotImplementedError
(
"Loading of timestamps file is not yet available."
)
def
_ijkseed
(
self
,
ii
,
jj
,
kk
,
inds
):
'''
Define the ii, jj, kk indices from which particles will be
released.
ii, jj, kk: list, array, or scalar with indices (i, j, k) from which
particles will be released (must be integers).
inds: if True (default), ii, jj and kk are interpreted
as indices, otherwise they are interpreted as lists of
exact release positions'''
ii
=
np
.
atleast_1d
(
np
.
squeeze
(
ii
))
jj
=
np
.
atleast_1d
(
np
.
squeeze
(
jj
))
kk
=
np
.
atleast_1d
(
np
.
squeeze
(
kk
))
if
(
ii
.
size
!=
jj
.
size
)
or
(
ii
.
size
!=
kk
.
size
):
raise
ValueError
(
"ii, jj and kk must have all the same dimension."
)
if
inds
:
return
np
.
atleast_2d
(
np
.
asarray
(
zip
(
np
.
int16
(
ii
),
np
.
int16
(
jj
),
np
.
int16
(
kk
)))),
\
np
.
atleast_2d
(
np
.
asarray
(
zip
(
ii
,
jj
,
kk
)))
else
:
# if MITgcm coordinates are passed, we have to load the model grid
# and then translate into the "normalised" index coordinates of
# tracmass
def
xy2grid
(
x
,
y
,
dX
,
dY
,
C
,
S
):
nx
=
(
C
*
x
+
S
*
y
)
/
dX
ny
=
(
-
S
*
x
+
C
*
y
)
/
dY
return
nx
,
ny
xG
=
self
.
xG
yG
=
self
.
yG
dX
=
self
.
dX
dY
=
self
.
dY
if
self
.
gcm_geometry
==
"curvilinear"
:
cs
=
self
.
CS
sn
=
self
.
SN
dZ
=
self
.
dzt
zG
=
self
.
Z
iout
=
np
.
zeros
(
ii
.
size
)
*
np
.
nan
jout
=
np
.
zeros
(
ii
.
size
)
*
np
.
nan
kout
=
np
.
zeros
(
ii
.
size
)
*
np
.
nan
iindout
=
np
.
zeros
(
ii
.
size
,
dtype
=
"int16"
)
jindout
=
np
.
zeros
(
ii
.
size
,
dtype
=
"int16"
)
kindout
=
np
.
zeros
(
ii
.
size
,
dtype
=
"int16"
)
for
nn
,
(
xx
,
yy
,
zz
)
in
enumerate
(
zip
(
ii
,
jj
,
kk
)):
jseed
,
iseed
=
np
.
unravel_index
(
np
.
argmin
((
xx
-
xG
)
**
2
+
(
yy
-
yG
)
**
2
),
xG
.
shape
)
if
self
.
gcm_geometry
==
"curvilinear"
:
nx
,
ny
=
xy2grid
(
xx
-
xG
[
jseed
,
iseed
],
yy
-
yG
[
jseed
,
iseed
],
dX
[
jseed
,
iseed
],
dY
[
jseed
,
iseed
],
cs
[
jseed
,
iseed
],
sn
[
jseed
,
iseed
])
else
:
nx
=
(
xx
-
xG
[
jseed
,
iseed
])
/
dX
[
jseed
,
iseed
]
ny
=
(
yy
-
yG
[
jseed
,
iseed
])
/
dY
[
jseed
,
iseed
]
# check the computed nx and ny, and in case seek in
# nearby cell
if
(
nx
<
0
)
or
(
nx
>
1
)
or
(
ny
<
0
)
or
(
ny
>
1
):
for
ijtry
in
product
(
range
(
-
5
,
6
),
range
(
-
5
,
6
)):
inew
=
iseed
+
ijtry
[
0
]
jnew
=
jseed
+
ijtry
[
1
]
if
(
inew
<
0
)
or
(
jnew
<
0
)
\
or
(
inew
>=
self
.
imax
)
or
(
jnew
>=
self
.
jmax
):
continue
if
self
.
gcm_geometry
==
"curvilinear"
:
nx
,
ny
=
xy2grid
(
xx
-
xG
[
jnew
,
inew
],
yy
-
yG
[
jnew
,
inew
],
dX
[
jnew
,
inew
],
dY
[
jnew
,
inew
],
cs
[
jnew
,
inew
],
sn
[
jnew
,
inew
])
else
:
nx
=
(
xx
-
xG
[
jnew
,
inew
])
/
dX
[
jnew
,
inew
]
ny
=
(
yy
-
yG
[
jnew
,
inew
])
/
dY
[
jnew
,
inew
]
if
(
nx
>=
0
)
and
(
nx
<
1
)
and
(
ny
>=
0
)
and
(
ny
<
1
):
iseed
=
inew
jseed
=
jnew
break
else
:
raise
ValueError
(
"Could not find the point "
"(x=
%f
, y=
%f
)"
%
(
xx
,
yy
))
z_here
=
zG
[:,
jseed
,
iseed
]
if
(
zz
>
z_here
.
max
())
or
(
zz
<=
z_here
.
min
()):
print
(
"
\n
Point outside vertical bounds at "
"x,y,z=
%.2f
,
%.2f
,
%.2f
"
%
(
xx
,
yy
,
zz
))
continue
kk
=
np
.
where
(
z_here
>
zz
)[
0
][
-
1
]
nz
=
(
zz
-
z_here
[
kk
])
/
(
z_here
[
kk
+
1
]
-
z_here
[
kk
])
iindout
[
nn
]
=
iseed
jindout
[
nn
]
=
jseed
kindout
[
nn
]
=
kk
iout
[
nn
]
=
iseed
+
nx
jout
[
nn
]
=
jseed
+
ny
kout
[
nn
]
=
kk
+
nz
return
np
.
atleast_2d
(
np
.
asarray
(
zip
(
iindout
[
np
.
isfinite
(
iout
)],
jindout
[
np
.
isfinite
(
jout
)],
kindout
[
np
.
isfinite
(
kout
)]))),
\
np
.
atleast_2d
(
np
.
asarray
(
zip
(
iout
[
np
.
isfinite
(
iout
)],
jout
[
np
.
isfinite
(
jout
)],
kout
[
np
.
isfinite
(
kout
)])))
# end _ijkseed
def
run
(
self
):
pass
Event Timeline
Log In to Comment