controller.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Tue, Aug 6, 02:20

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("\nUnsupported configuration file extension.")
	raise

	print("\nConfiguration 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("\nUnsupported configuration file content.")
	raise

	try:
	self.check_config()
	except:
	print("\nInvalid configuration setup.")
	raise

	print("\nInitialization 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("\nReference date of GCM:", self.gcm_start)
	self.start = np.datetime64(self.config["Sim"]["start"])
	print("\nBegin particle tracking simulation at:", self.start)
	self.end = np.datetime64(self.config["Sim"]["end"])
	print("\nEnd particle tracking simulation at:", self.end)
	self.seed_start = np.datetime64(self.config["Sim"]["seed_start"])
	print("\nStart seeding particles at:", self.seed_start)
	self.seed_end = np.datetime64(self.config["Sim"]["seed_end"])
	print("\nStop seeding particles at:", self.seed_end)

	self.is2D = self.config["Sim"]["is2D"];
	print("\nSimulation in 2D : ", self.is2D)


	#
	# prepare the simulation
	#
	print("\nIdentify 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("\nCurvilinear 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("\nCartesian 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("\nNumber 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("\nSeeding interval: %ds" % 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("\nNumber of seeding time steps: %d" % self.iter_seed.size)

	self.ntrackmax = self.iter_seed.size * self.n_seed
	print("\nTotal 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("\nRunning numerics single threaded")
	else:
	print("\nRunning numerics on {0} threads max".format(self.n_threads))
	self.n_min_part_per_thread = self.config["Parallelism"]["n_min_part_per_thread"]
	print("\nRunning 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("\nPoint 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

controller.pyNo OneTemporaryActions

File Metadata

controller.pyView Options

Event Timeline

controller.py
No OneTemporary
Actions

controller.py
View Options