X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDRO%2FTELEMAC2D.py;fp=src%2FHYDRO%2FTELEMAC2D.py;h=85f1f01700d416dc68f244c43fc59a984a62d87e;hb=52fa018273f662e7b5ec8e82b047a87b28e39189;hp=0000000000000000000000000000000000000000;hpb=8facaaa341432a098c7f2b73388ae31282b3ad5a;p=modules%2Fhydrosolver.git diff --git a/src/HYDRO/TELEMAC2D.py b/src/HYDRO/TELEMAC2D.py new file mode 100644 index 0000000..85f1f01 --- /dev/null +++ b/src/HYDRO/TELEMAC2D.py @@ -0,0 +1,579 @@ +# Copyright (C) 2012-2013 EDF +# +# This file is part of SALOME HYDRO module. +# +# SALOME HYDRO module 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. +# +# SALOME HYDRO module 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 SALOME HYDRO module. If not, see . + +import sys +import os +import logging +import threading +import inspect +import traceback +import cPickle +import math +import types + +import salome +import HYDROSOLVER_ORB__POA +import SALOME_ComponentPy +import SALOME +import calcium +import dsccalcium + +from salome.kernel.logger import Logger +from salome.kernel import termcolor +logger = Logger("TELEMAC2D", color = termcolor.BLUE) +logger.setLevel(logging.DEBUG) + +from salome.kernel.parametric.compo_utils import \ + create_input_dict, create_normal_parametric_output, create_error_parametric_output + +from TelApy.api.t2d import Telemac2D +from TelApy.tools.polygon import get_list_points_in_polygon +from salome.hydro.study import jdc_to_dict, get_jdc_dict_var_as_tuple, HydroStudyEditor + +################################################ + +class TELEMAC2D(HYDROSOLVER_ORB__POA.TELEMAC2D, dsccalcium.PyDSCComponent): + + lock = threading.Lock() + + """ + Pour etre un composant SALOME cette classe Python + doit avoir le nom du composant et heriter de la + classe HYDRO_ORB__POA.MASCARET issue de la compilation de l'idl + par omniidl et de la classe SALOME_ComponentPy_i + qui porte les services generaux d'un composant SALOME + """ + def __init__ ( self, orb, poa, contID, containerName, instanceName, + interfaceName ): + logger.info("__init__: " + containerName + ' ; ' + instanceName) + dsccalcium.PyDSCComponent.__init__(self, orb, poa,contID, + containerName,instanceName,interfaceName) + self.t2d = None + self.last_start_time = None + self.saved_state = None + + def init_service(self,service): + return service in ("Compute",) + + def _raiseSalomeError(self): + message = "Error in component %s running in container %s." % (self._instanceName, self._containerName) + logger.exception(message) + message += " " + traceback.format_exc() + exc = SALOME.ExceptionStruct(SALOME.INTERNAL_ERROR, message, + inspect.stack()[1][1], inspect.stack()[1][2]) + raise SALOME.SALOME_Exception(exc) + + def Init(self, studyId, detCaseEntry): + try: + self.beginService("TELEMAC2D.Init") + self.jdc_dict = self.get_jdc_dict_from_case(studyId, detCaseEntry) + + # Ugly hack to check if we are in OpenTURNS context or Mascaret coupling context + if "VARIABLE_ENTREE" in self.jdc_dict: + # Create variable mappings for OpenTURNS + self.input_var_map = {} + for input_var in get_jdc_dict_var_as_tuple(self.jdc_dict, "VARIABLE_ENTREE"): + name = input_var["NOM"] + if not self.input_var_map.has_key(name): + self.input_var_map[name] = [] + self.input_var_map[name] += [Telemac2DInputVariable(input_var["VARIABLE_MODELE_T2D"])] + + self.output_var_map = {} + for output_var in get_jdc_dict_var_as_tuple(self.jdc_dict, "VARIABLE_SORTIE"): + t2d_output_var = Telemac2DOutputVariable(output_var["VARIABLE_T2D"]) + self.output_var_map[output_var["NOM"]] = t2d_output_var + + else: + # Initialize before launching the coupling with Mascaret + self.init_t2d_from_jdc_dict(self.jdc_dict) + + self.endService("TELEMAC2D.Init") + except: + self._raiseSalomeError() + + def get_jdc_dict_from_case(self, study_id, case_entry): + TELEMAC2D.lock.acquire() + salome.salome_init() + TELEMAC2D.lock.release() + + # Get filenames from the case in Salome study + ed = HydroStudyEditor(study_id) + sobj = ed.editor.study.FindObjectID(case_entry) + jdcpath = sobj.GetComment() + with open(jdcpath) as jdcfile: + jdc = jdcfile.read() + return jdc_to_dict(jdc, ["TELEMAC2D", "_F"]) + + def init_t2d_from_jdc_dict(self, jdc_dict): + steering_filename = jdc_dict["FICHIER_CAS"] + user_fortran = jdc_dict.get("FICHIER_FORTRAN_UTILISATEUR") + dico_filename = jdc_dict.get("FICHIER_DICO") + geo_filename = jdc_dict.get("FICHIER_GEOMETRIE") + bc_filename = jdc_dict.get("FICHIER_CONDITIONS_LIMITES") + result_filename = jdc_dict.get("RESULTAT") + + logger.debug("Initializing Telemac2D API") + self.t2d = Telemac2D(user_fortran) + logger.debug("Telemac2D API initialization OK") + os.chdir(os.path.dirname(steering_filename)) + self.t2d.run_read_case_t2d(steering_filename, dico_filename) + if bc_filename is not None: + self.t2d.set_string('MODEL.BCFILE', bc_filename) + if geo_filename is not None: + self.t2d.set_string('MODEL.GEOMETRYFILE', geo_filename) + if result_filename is not None: + self.t2d.set_string('MODEL.RESULTFILE', result_filename) + ierr = self.t2d.api.run_allocation_t2d(self.t2d.id) + if ierr != 0: + raise Exception('Error in run_allocation_t2d:', ierr) + ierr = self.t2d.api.run_init_t2d(self.t2d.id) + if ierr != 0: + raise Exception('Error in run_init_t2d:', ierr) + + def Exec(self, paramInput): + try: + self.beginService("TELEMAC2D.Exec") + + # Save / restore state + # Not used yet because time can not be reset + #if self.saved_state is None: + # self.saved_state = self.t2d.save_state() + #else: + # self.t2d.restore_state(self.saved_state) + + # Get id and execution mode + id = "" + exec_mode = "" + for parameter in paramInput.specificParameters: + if parameter.name == "id": + id = parameter.value + if parameter.name == "executionMode": + exec_mode = parameter.value + logger.debug("ID: %s" % id) + logger.debug("Execution mode: %s" % exec_mode) + + # Append id to result file name + result_filename = self.jdc_dict.get("RESULTAT") or "r2d.slf" + (result_filename_wo_ext, ext) = os.path.splitext(result_filename) + current_jdc_dict = self.jdc_dict.copy() + current_jdc_dict["RESULTAT"] = "%s_%s%s" % (result_filename_wo_ext, id, ext) + + # Reload model + self.init_t2d_from_jdc_dict(current_jdc_dict) + + # Set the input values + logger.debug("inputVarList: %s" % paramInput.inputVarList) + logger.debug("outputVarList: %s" % paramInput.outputVarList) + logger.debug("inputValues: %s" % paramInput.inputValues) + + input_dict = create_input_dict({}, paramInput) + logger.debug("input_dict = %s" % input_dict) + + for (key, value) in input_dict.iteritems(): + for t2d_var in self.input_var_map[key]: + t2d_var.set_value(self.t2d, value) + + # Compute + self.t2d.run_all_timesteps() + + # Get the output values + output_dict = {} + for output_var in paramInput.outputVarList: + output_dict[output_var] = self.output_var_map[output_var].get_value(self.t2d) + + param_output = create_normal_parametric_output(output_dict, paramInput) + logger.debug("outputValues: %s" % param_output.outputValues) + + # Finalize T2D + self.t2d.finalize() + self.t2d = None + + self.endService("TELEMAC2D.Exec") + return param_output + except: + self._raiseSalomeError() + + def InitBorders(self, bordersPickled): + """ + This method initializes the border(s) for coupling + """ + try: + self.beginService("TELEMAC2D.InitBorders") + borders = cPickle.loads(bordersPickled) + self.borders = {} + for border in borders: + if border["model1"] is not None and border["model1"]["code"] == "TELEMAC2D": + self.borders[border["name"]] = border["model1"] + elif border["model2"] is not None and border["model2"]["code"] == "TELEMAC2D": + self.borders[border["name"]] = border["model2"] + self.endService("TELEMAC2D.InitBorders") + except: + self._raiseSalomeError() + + def ExecStep(self, timeDataPickled, inputDataPickled): + """ + This method computes a single time step + """ + try: + self.beginService("TELEMAC2D.ExecStep") + timeData = cPickle.loads(timeDataPickled) + inputData = cPickle.loads(inputDataPickled) + + if timeData["start_time"] == self.last_start_time: + # Convergence iteration, restore saved state + self.t2d.restore_state(self.saved_state) + else: + # First iteration of a time step, save state + self.saved_state = self.t2d.save_state(self.saved_state) + + self.last_start_time = timeData["start_time"] + + if inputData is not None: + for (name, desc) in self.borders.iteritems(): + if name in inputData: + z = None + v = None + q = None + if desc["position"] == "upstream": + if inputData[name]["Froude"] < 1.0: # Fluvial + z = inputData[name]["Z"] + else: + v = inputData[name]["V"] + q = inputData[name]["Q"] + if inputData[name]["Froude"] > 1.0: # Torrential + z = inputData[name]["Z"] + self.set_z_on_border(desc, z) + #self.set_v_on_border(desc, v) + self.set_q_on_border(desc, q) + + # Compute + # TODO: Use time values from coupling + ierr = self.t2d.api.run_timestep_t2d(self.t2d.id) + if ierr != 0: + raise Exception('Error in run_timestep_t2d:', ierr) + + # Get the output values + outputData = {} + for (name, desc) in self.borders.iteritems(): + outputData[name] = {} + (outputData[name]["Z"], idx) = self.get_z_on_border(desc) + if idx is None: + outputData[name]["V"] = 0.0 + else: + (outputData[name]["V"], froude) = self.get_v_at_point(idx) + outputData[name]["Q"] = self.get_q_on_border(desc) + outputDataPickled = cPickle.dumps(outputData, -1) + + self.endService("TELEMAC2D.ExecStep") + return outputDataPickled + except: + self._raiseSalomeError() + + def get_z_on_border(self, border): + """ + Find the water elevation (i.e. water height + bottom elevation) on a border. + Return a tuple (water_elevation, index_of_model_point_with_min_water_elevation). + If we use method 2 (compute mean water elevation), index is set to None. + """ + # Method 1: Get point with minimum elevation + """ + idx = -1 + bor_idx = -1 + z_min = sys.float_info.max + (first, after_last) = border["points_id"] + k = first + while k != after_last: + nbor_k = self.t2d.get_integer('MODEL.NBOR', k) + h_k = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k) + z_bottom_k = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_k) + z_k = h_k + z_bottom_k + logger.debug("extract z for border point %d, h:%f, zf:%f, z:%f" % (k, h_k, z_bottom_k, z_k)) + if (z_k < z_min): + z_min = z_k + bor_idx = k + idx = nbor_k + k = self.t2d.get_integer('MODEL.KP1BOR', k) + logger.debug("min z for border %d with method 1 is %f at border point %d (model point %d)" % + (border["border_id"], z_min, bor_idx, idx)) + z = z_min + """ + # Method 2: Compute mean elevation, ignoring dry zones (method from PALM coupling) + z_mean = sys.float_info.max + for i in range(10): + h_tot = 0.0 + l_tot = 0.0 + zf_min = sys.float_info.max + (first, after_last) = border["points_id"] + k = first + while k != after_last: + knext = self.t2d.get_integer('MODEL.KP1BOR', k) + nbor_k = self.t2d.get_integer('MODEL.NBOR', k) + nbor_knext = self.t2d.get_integer('MODEL.NBOR', knext) + zf1 = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_k) + zf2 = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_knext) + zf_min = min(zf_min, zf1) + h1 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k) + h2 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_knext) + z1 = h1 + zf1 + z2 = h2 + zf2 + x1 = self.t2d.get_double('MODEL.X', nbor_k) + x2 = self.t2d.get_double('MODEL.X', nbor_knext) + y1 = self.t2d.get_double('MODEL.Y', nbor_k) + y2 = self.t2d.get_double('MODEL.Y', nbor_knext) + l = math.sqrt((x1-x2)**2 + (y1-y2)**2) + if zf1 <= z_mean and zf2 <= z_mean: + h_tot += (z1+z2) * 0.5 * l + l_tot += l + k = knext + + new_z_mean = h_tot / l_tot + delta = abs(z_mean - new_z_mean) + z_mean = new_z_mean + if delta < 1.0e-6: + break + + logger.debug("mean z for border %d with method 2 is %f" % (border["border_id"], z_mean)) + z = z_mean + idx = None + + return (z, idx) + + def get_v_at_point(self, idx): + """ + Find the water velocity at the model point idx. + Return a tuple (velocity, froude number) + """ + u_i = self.t2d.get_double('MODEL.VELOCITYU', idx) + v_i = self.t2d.get_double('MODEL.VELOCITYV', idx) + h_i = self.t2d.get_double('MODEL.WATERDEPTH', idx) + v = math.sqrt(u_i**2 + v_i**2) + ##### TODO: Remove that and find the real error + if h_i <= 0.0: + print 'Error: water depth is negative or zero on point', idx + #raise Exception('Null water depth value') + h_i = 0.01 + froude = v / (h_i * 9.81) + return (v, froude) + + def get_q_on_border(self, border): + """ + Find the flow on a border. + The returned value is positive if the flow goes in the natural direction + according to the borders description (from upstream to downstream). + """ + # Method 1: Integrate the flow on the whole border + """ + q = 0.0 + (first, after_last) = border["points_id"] + k = first + while k != after_last: + knext = self.t2d.get_integer('MODEL.KP1BOR', k) + nbor_k = self.t2d.get_integer('MODEL.NBOR', k) + nbor_knext = self.t2d.get_integer('MODEL.NBOR', knext) + h_i = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k) + h_i2 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_knext) + h2d1 = 0.5 * (h_i + h_i2) + if (h2d1 > 0.0): + u_i = self.t2d.get_double('MODEL.VELOCITYU', nbor_k) + u_i2 = self.t2d.get_double('MODEL.VELOCITYU', nbor_knext) + v_i = self.t2d.get_double('MODEL.VELOCITYV', nbor_k) + v_i2 = self.t2d.get_double('MODEL.VELOCITYV', nbor_knext) + v2d1 = 0.5 * (math.sqrt(u_i**2+v_i**2) + math.sqrt(u_i2**2+v_i2**2)) + x2d1 = self.t2d.get_double('MODEL.X', nbor_k) + x2d2 = self.t2d.get_double('MODEL.X', nbor_knext) + y2d1 = self.t2d.get_double('MODEL.Y', nbor_k) + y2d2 = self.t2d.get_double('MODEL.Y', nbor_knext) + q += v2d1 * h2d1 * math.sqrt((x2d1-x2d2)**2 + (y2d1-y2d2)**2) + k = knext + logger.debug("flow for border %d with method 1 is %f" % (border["border_id"], q)) + + # Method 2: Get the flow from Telemac model + q = self.t2d.get_double('MODEL.FLUX_BOUNDARIES', border["border_id"]) + if border["position"] == "downstream": + q = -q + logger.debug("flow for border %d with method 2 is %f" % (border["border_id"], q)) + """ + + # Method 3: Method from PALM coupling + q = 0.0 + (first, after_last) = border["points_id"] + k = first + while k != after_last: + knext = self.t2d.get_integer('MODEL.KP1BOR', k) + nbor_k = self.t2d.get_integer('MODEL.NBOR', k) + nbor_knext = self.t2d.get_integer('MODEL.NBOR', knext) + h_i = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k) + h_i2 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_knext) + u_i = self.t2d.get_double('MODEL.VELOCITYU', nbor_k) + u_i2 = self.t2d.get_double('MODEL.VELOCITYU', nbor_knext) + v_i = self.t2d.get_double('MODEL.VELOCITYV', nbor_k) + v_i2 = self.t2d.get_double('MODEL.VELOCITYV', nbor_knext) + x2d1 = self.t2d.get_double('MODEL.X', nbor_k) + x2d2 = self.t2d.get_double('MODEL.X', nbor_knext) + y2d1 = self.t2d.get_double('MODEL.Y', nbor_k) + y2d2 = self.t2d.get_double('MODEL.Y', nbor_knext) + NX=y2d1-y2d2 + NY=x2d2-x2d1 + UN1=u_i*NX+v_i*NY + UN2=u_i2*NX+v_i2*NY + q += ((h_i+h_i2)*(UN1+UN2)+h_i2*UN2+h_i*UN1)/6.0 + k = knext + if border["position"] == "upstream": + q = -q + logger.debug("flow for border %d with method 3 is %f" % (border["border_id"], q)) + + return q + + def set_z_on_border(self, border, z): + """ + Set the water elevation (i.e. water height + bottom elevation) on a border. + If z is None, then the water elevation for this border is free (not imposed). + """ + # Method 1: Impose water elevation point by point + """ + # Deactivate water level imposition in bord.f (requires a patch in bord.f) + self.t2d.set_double('MODEL.COTE', -1.0, border["border_id"]) + + (first, after_last) = border["points_id"] + k = first + while k != after_last: + if z is None: + self.t2d.set_integer('MODEL.LIHBOR', 4, k) # Free elevation on that border + else: + self.t2d.set_integer('MODEL.LIHBOR', 5, k) # Imposed elevation on that border + nbor_k = self.t2d.get_integer('MODEL.NBOR', k) + bottom_elevation = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_k) + hbor_k = max(0.0, z - bottom_elevation) + self.t2d.set_double('MODEL.HBOR', hbor_k, k) + logger.debug("set z for border point %d, h:%f, zf:%f, z:%f" % (k, hbor_k, bottom_elevation, z)) + k = self.t2d.get_integer('MODEL.KP1BOR', k) + """ + + # Method 2: Just tell Telemac what is the imposed elevation for the border + if z is not None: + self.t2d.set_double('MODEL.COTE', z, border["border_id"]) + (first, after_last) = border["points_id"] + k = first + while k != after_last: + if z is None: + self.t2d.set_integer('MODEL.LIHBOR', 4, k) # Free elevation on that border + else: + self.t2d.set_integer('MODEL.LIHBOR', 5, k) # Imposed elevation on that border + k = self.t2d.get_integer('MODEL.KP1BOR', k) + + def set_v_on_border(self, border, v): + """ + Set the water velocity on a border. + If v is None, then the water velocity for this border is free (not imposed). + """ + # This method works only if the water level imposition is deactivated in bord.f + # (requires a patch in bord.f) + (first, after_last) = border["points_id"] + k = first + while k != after_last: + if v is None: + self.t2d.set_integer('MODEL.LIUBOR', 4, k) # Free velocity on that border + self.t2d.set_integer('MODEL.LIVBOR', 4, k) + else: + self.t2d.set_integer('MODEL.LIUBOR', 6, k) # Imposed velocity on that border + self.t2d.set_integer('MODEL.LIVBOR', 6, k) + xnebor_k = self.t2d.get_double('MODEL.XNEBOR', k) + ubor_k = -xnebor_k * v + ynebor_k = self.t2d.get_double('MODEL.YNEBOR', k) + vbor_k = -ynebor_k * v + self.t2d.set_double('MODEL.UBOR', ubor_k, k) + self.t2d.set_double('MODEL.VBOR', vbor_k, k) + k = self.t2d.get_integer('MODEL.KP1BOR', k) + + def set_q_on_border(self, border, q): + """ + Set the water flow on a border. + If q is None, then the water flow for this border is free (not imposed). + """ + if q is not None: + self.t2d.set_double('MODEL.DEBIT', q, border["border_id"]) + logger.debug("flow for border %d is set to %f" % (border["border_id"], q)) + (first, after_last) = border["points_id"] + k = first + while k != after_last: + if q is None: + self.t2d.set_integer('MODEL.LIUBOR', 4, k) # Free velocity on that border + self.t2d.set_integer('MODEL.LIVBOR', 4, k) + else: + self.t2d.set_integer('MODEL.LIUBOR', 5, k) # Imposed discharge on that border + self.t2d.set_integer('MODEL.LIVBOR', 5, k) + k = self.t2d.get_integer('MODEL.KP1BOR', k) + + def Finalize(self): + try: + self.beginService("TELEMAC2D.Finalize") + if self.t2d: + self.t2d.finalize() + self.t2d = None + self.endService("TELEMAC2D.Finalize") + except: + self._raiseSalomeError() + + def Compute(self, studyId, caseEntry): + try: + self.beginService("TELEMAC2D.Compute") + + jdc_dict = self.get_jdc_dict_from_case(studyId, caseEntry) + self.init_t2d_from_jdc_dict(jdc_dict) + + # Compute + self.t2d.run_all_timesteps() + + # Cleanup + self.t2d.finalize() + + self.endService("TELEMAC2D.Compute") + except: + self._raiseSalomeError() + + +class Telemac2DOutputVariable: + + def __init__(self, varstr): + par1 = varstr.find("(") + par2 = varstr.find(")") + self.varname = varstr[:par1] + self.idx = int(varstr[par1+1:par2]) + + def get_value(self, inst): + return inst.get_double(self.varname, self.idx) + + +class Telemac2DInputVariable: + + def __init__(self, vardict): + self.varname = vardict["NOM_VARIABLE_MODELE"] + def_area_dict = vardict["ZONE_DE_DEFINITION"] + self.idx = def_area_dict.get("INDICE") + self.polygon = def_area_dict.get("POLYGONE") + self.points_in_polygon_list = None + + def set_value(self, inst, value): + if self.idx is not None: + inst.set_double(self.varname, value, self.idx) + elif self.polygon is not None: + if self.points_in_polygon_list is None: + self.points_in_polygon_list = get_list_points_in_polygon(inst, self.polygon) + for idx in self.points_in_polygon_list: + inst.set_double(self.varname, value, idx) + else: + raise Exception("Invalid area definition for input variable {0}".format(self.varname))