+++ /dev/null
-# 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 <http://www.gnu.org/licenses/>.
-
-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))
