Homec4science
Diffusion ctracker3 (master) Browse

ctracker3/mods/particleContainer.py
2cfb5a04471amaster

Last Change
Raw File

particleContainer.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 numpy as np
from queue import Queue
from core.parameterContainer import ParameterContainer
from parameters.particle import Particle
from tools.CFunctions import ext_loop
class ParticleContainer(ParameterContainer):
def __init__(self, name, data, grid):
super(ParticleContainer, self).__init__(name, data)
'''
Initializes the seeding frequency and position and seeding steps for the simulation
'''
self.outfreq = data.outfreq
self.ff = float(data.ff)
self.gcm_geometry = data.gcm_geometry
self.grid = grid
self.gcm_dt = data.gcm_dt
self.out_dt = data.out_dt
self.gcm_start = data.gcm_start
self.iende = np.asarray(data.iende)
self.iendw = np.asarray(data.iendw)
self.jendn = np.asarray(data.jendn)
self.jends = np.asarray(data.jends)
self.nend = self.iende.size
self.max_id = 0
# Total number of active particles
self.ntact = 0
# Number of particles which exited the domain
self.ntout = 0
# this is a copy in memory of the NETcdf time
# it is useful to find the index where to write
# data without reading from disk
self.nc_time = np.array([], dtype="float64")
self.queue_in = Queue(maxsize=3)
# self.queue_out = Queue(maxsize=50)
self.ntact = 0
self.active = np.array([], dtype=Particle)
self.exited = np.array([], dtype=Particle)
self._children = np.array([], dtype=Particle)
self.nvals = 1
self.is_running = False
self.inputs["default"] = [self.process_input]
self.outputs = ['default']
def process_input(self, *args, **kwargs):
self.queue_in.put(*args, **kwargs)
def process_output(self, *args, **kwargs):
pass
"""data = [1, 2, 3]
print("ParticleContainer process output with {0}".format(data))
self.output('default', data, *args, **kwargs)"""
def is_config_valid(self):
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.nend != self.iendw.size) or \
(self.nend != self.jendn.size) or \
(self.nend != self.jends.size):
raise ValueError("Wrong kill area definition")
return True
def run(self, seed_ctx):
self.is_running = True
#
# Start time loop over GCM output steps
#
# array with particle ids
active_ids = np.zeros(self.ntact, dtype="short")
out_code = np.zeros(self.ntact, dtype="short")
# input/output position of particles
xyz = np.zeros((0, 3), dtype="float64")
ijk = np.zeros((0, 3), dtype="int16")
self.nvals = seed_ctx.nvals
out_tijk = np.zeros((self.ntact, self.nvals, 4),
dtype="f8")
out_xyz = np.zeros((self.ntact, self.nvals, 3),
dtype="f8")
iters = seed_ctx.iters
timestamps_in = seed_ctx.timestamps_in
timestamps_out = seed_ctx.timestamps_out
all_stops = np.sort(np.unique(np.concatenate((iters, timestamps_in, timestamps_out))))
# initialize stuff
ts = 0.0 # normalised time between two GCM time steps
tt = 0.0 # time in seconds between two GCM time steps
for ngi in range(1, all_stops.size):
# we are starting from the start!
# even if we will read the following
# gcm output for interpolating
t = all_stops[ngi - 1]
ngiQ = uflux = vflux = wflux = 0
t0 = t * self.gcm_dt
t1 = all_stops[ngi] * self.gcm_dt
dt = t1 - t0
# check if we reached a timestamp at which we should extract flux data
if t in iters:
# get data from the queue
ngiQ, uflux, vflux, wflux = self.queue_in.get()
ts = tt = 0.0
if ngiQ != ngi:
raise ValueError("Queue provides wrong timestep")
# remove particles that went out of the domain or other
del_index = out_code == 0
n_remove = (~del_index).sum()
self.ntout += n_remove
self.ntact = active_ids.size
if n_remove:
active_ids = active_ids[del_index]
xyz = xyz[del_index, ...]
ijk = ijk[del_index, ...]
self.ntout += n_remove
self.ntact = active_ids.size
# the active particles number changed, so
# we update the buffer size
out_code = np.zeros(self.ntact, dtype="short")
out_tijk = np.zeros((self.ntact, seed_ctx.nvals, 4),
dtype="f8")
out_xyz = np.zeros((self.ntact, seed_ctx.nvals, 3),
dtype="f8")
if t in timestamps_in:
new_ids = np.arange(self.max_id + 1,
self.max_id + 1 + seed_ctx.n_seed,
dtype="int64")
active_ids = np.append(active_ids, new_ids)
xyz = np.vstack([xyz,
seed_ctx.ijk_seed])
ijk = np.vstack([ijk,
np.int16(seed_ctx.ijk_indseed)])
self.max_id = self.max_id + seed_ctx.n_seed
# we store the seeding position
self.output('default', (t0, None, None, new_ids.copy(), None, True))
self.ntact = active_ids.size
# the active particles number changed, so
# we update the buffer size
out_code = np.zeros(self.ntact, dtype="short")
out_tijk = np.zeros((self.ntact, seed_ctx.nvals, 4),
dtype="f8")
out_xyz = np.zeros((self.ntact, seed_ctx.nvals, 3),
dtype="f8")
print("\nStarting numerical loop at ",
np.datetime64(self.gcm_start) +
np.timedelta64(int(t0 * 1000), "[ms]"),
"\n|.... active: {:d} ....|.... out: {:d} ....|"
.format(self.ntact, self.ntout))
if self.outfreq == "gcmstep":
outfreq = 1
elif self.outfreq == "always":
outfreq = 2
elif self.outfreq == "cross":
outfreq = 3
elif self.outfreq == "custom":
outfreq = 4
else:
raise ValueError("outfreq {} not recognised"
.format(self.outfreq))
ext_loop(self.nvals, self.ntact,
xyz, ijk, out_xyz, out_tijk,
timestamps_in, timestamps_out, t, [ts, tt],
dt, outfreq,
self.nend, self.iendw, self.iende, self.jendn, self.jends,
self.grid.imax, self.grid.jmax, self.grid.kmax,
self.grid.dsmax, self.grid.dxyz, self.grid.dtmax,
self.grid.dstep, self.grid.CS, self.grid.SN, self.grid.dX,
self.grid.dY, self.grid.xG, self.grid.yG, self.grid.Z,
uflux, vflux, wflux, out_code)
# instead of doing the netCDF writing here, we rather
# put the output in a queue, a separate thread will
# deal with it
# Note that we must store copies, otherwise we would be
# changing the output before it's written
self.output('default',
(t0,
out_tijk.copy(),
out_xyz.copy(),
active_ids.copy(),
out_code.copy(),
False)
)
if t in timestamps_out:
timestamps_out.pop(0)
if t in timestamps_in:
timestamps_in.pop(0)
# end for cycle over iterations
# Tell the output the processing is finished
self.output('default', None)
# finish the output queue before closing the netcdf output file
# self.queue_out.join()
c4science · Help