1 # Copyright (C) 2012-2013 EDF
3 # This file is part of SALOME HYDRO module.
5 # SALOME HYDRO module is free software: you can redistribute it and/or modify
6 # it under the terms of the GNU General Public License as published by
7 # the Free Software Foundation, either version 3 of the License, or
8 # (at your option) any later version.
10 # SALOME HYDRO module is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 # GNU General Public License for more details.
15 # You should have received a copy of the GNU General Public License
16 # along with SALOME HYDRO module. If not, see <http://www.gnu.org/licenses/>.
29 import HYDROSOLVER_ORB__POA
30 import SALOME_ComponentPy
35 from salome.kernel.logger import Logger
36 from salome.kernel import termcolor
37 logger = Logger("TELEMAC2D", color = termcolor.BLUE)
38 logger.setLevel(logging.DEBUG)
40 from salome.kernel.parametric.compo_utils import \
41 create_input_dict, create_normal_parametric_output, create_error_parametric_output
43 from salome.hydro.telemac2d import Telemac2D
44 from salome.hydro.telemac2d.polygon import get_list_points_in_polygon
45 from salome.hydro.study import jdc_to_dict, get_jdc_dict_var_as_tuple, HydroStudyEditor
47 ################################################
49 class TELEMAC2D(HYDROSOLVER_ORB__POA.TELEMAC2D, dsccalcium.PyDSCComponent):
51 lock = threading.Lock()
54 Pour etre un composant SALOME cette classe Python
55 doit avoir le nom du composant et heriter de la
56 classe HYDRO_ORB__POA.MASCARET issue de la compilation de l'idl
57 par omniidl et de la classe SALOME_ComponentPy_i
58 qui porte les services generaux d'un composant SALOME
60 def __init__ ( self, orb, poa, contID, containerName, instanceName,
62 logger.info("__init__: " + containerName + ' ; ' + instanceName)
63 dsccalcium.PyDSCComponent.__init__(self, orb, poa,contID,
64 containerName,instanceName,interfaceName)
66 self.last_start_time = None
67 self.saved_state = None
69 def init_service(self,service):
70 return service in ("Compute",)
72 def _raiseSalomeError(self):
73 message = "Error in component %s running in container %s." % (self._instanceName, self._containerName)
74 logger.exception(message)
75 message += " " + traceback.format_exc()
76 exc = SALOME.ExceptionStruct(SALOME.INTERNAL_ERROR, message,
77 inspect.stack()[1][1], inspect.stack()[1][2])
78 raise SALOME.SALOME_Exception(exc)
80 def Init(self, studyId, detCaseEntry):
82 self.beginService("TELEMAC2D.Init")
83 self.jdc_dict = self.get_jdc_dict_from_case(studyId, detCaseEntry)
85 # Ugly hack to check if we are in OpenTURNS context or Mascaret coupling context
86 if "VARIABLE_ENTREE" in self.jdc_dict:
87 # Create variable mappings for OpenTURNS
88 self.input_var_map = {}
89 for input_var in get_jdc_dict_var_as_tuple(self.jdc_dict, "VARIABLE_ENTREE"):
90 name = input_var["NOM"]
91 if not self.input_var_map.has_key(name):
92 self.input_var_map[name] = []
93 self.input_var_map[name] += [Telemac2DInputVariable(input_var["VARIABLE_MODELE_T2D"])]
95 self.output_var_map = {}
96 for output_var in get_jdc_dict_var_as_tuple(self.jdc_dict, "VARIABLE_SORTIE"):
97 t2d_output_var = Telemac2DOutputVariable(output_var["VARIABLE_T2D"])
98 self.output_var_map[output_var["NOM"]] = t2d_output_var
101 # Initialize before launching the coupling with Mascaret
102 self.init_t2d_from_jdc_dict(self.jdc_dict)
104 self.endService("TELEMAC2D.Init")
106 self._raiseSalomeError()
108 def get_jdc_dict_from_case(self, study_id, case_entry):
109 TELEMAC2D.lock.acquire()
111 TELEMAC2D.lock.release()
113 # Get filenames from the case in Salome study
114 ed = HydroStudyEditor(study_id)
115 sobj = ed.editor.study.FindObjectID(case_entry)
116 jdcpath = sobj.GetComment()
117 with open(jdcpath) as jdcfile:
119 return jdc_to_dict(jdc, ["TELEMAC2D", "_F"])
121 def init_t2d_from_jdc_dict(self, jdc_dict):
122 steering_filename = jdc_dict["FICHIER_CAS"]
123 user_fortran = jdc_dict.get("FICHIER_FORTRAN_UTILISATEUR")
124 dico_filename = jdc_dict.get("FICHIER_DICO")
125 geo_filename = jdc_dict.get("FICHIER_GEOMETRIE")
126 bc_filename = jdc_dict.get("FICHIER_CONDITIONS_LIMITES")
127 result_filename = jdc_dict.get("RESULTAT")
129 logger.debug("Initializing Telemac2D API")
130 self.t2d = Telemac2D(user_fortran)
131 logger.debug("Telemac2D API initialization OK")
132 os.chdir(os.path.dirname(steering_filename))
133 self.t2d.run_read_case_t2d(steering_filename, dico_filename)
134 if bc_filename is not None:
135 self.t2d.set_string('MODEL.BCFILE', bc_filename)
136 if geo_filename is not None:
137 self.t2d.set_string('MODEL.GEOMETRYFILE', geo_filename)
138 if result_filename is not None:
139 self.t2d.set_string('MODEL.RESULTFILE', result_filename)
140 ierr = self.t2d.api.run_allocation_t2d(self.t2d.id)
142 raise Exception('Error in run_allocation_t2d:', ierr)
143 ierr = self.t2d.api.run_init_t2d(self.t2d.id)
145 raise Exception('Error in run_init_t2d:', ierr)
147 def Exec(self, paramInput):
149 self.beginService("TELEMAC2D.Exec")
151 # Save / restore state
152 # Not used yet because time can not be reset
153 #if self.saved_state is None:
154 # self.saved_state = self.t2d.save_state()
156 # self.t2d.restore_state(self.saved_state)
158 # Get id and execution mode
161 for parameter in paramInput.specificParameters:
162 if parameter.name == "id":
164 if parameter.name == "executionMode":
165 exec_mode = parameter.value
166 logger.debug("ID: %s" % id)
167 logger.debug("Execution mode: %s" % exec_mode)
169 # Append id to result file name
170 result_filename = self.jdc_dict.get("RESULTAT") or "r2d.slf"
171 (result_filename_wo_ext, ext) = os.path.splitext(result_filename)
172 current_jdc_dict = self.jdc_dict.copy()
173 current_jdc_dict["RESULTAT"] = "%s_%s%s" % (result_filename_wo_ext, id, ext)
176 self.init_t2d_from_jdc_dict(current_jdc_dict)
178 # Set the input values
179 logger.debug("inputVarList: %s" % paramInput.inputVarList)
180 logger.debug("outputVarList: %s" % paramInput.outputVarList)
181 logger.debug("inputValues: %s" % paramInput.inputValues)
183 input_dict = create_input_dict({}, paramInput)
184 logger.debug("input_dict = %s" % input_dict)
186 for (key, value) in input_dict.iteritems():
187 for t2d_var in self.input_var_map[key]:
188 t2d_var.set_value(self.t2d, value)
191 self.t2d.run_all_timesteps()
193 # Get the output values
195 for output_var in paramInput.outputVarList:
196 output_dict[output_var] = self.output_var_map[output_var].get_value(self.t2d)
198 param_output = create_normal_parametric_output(output_dict, paramInput)
199 logger.debug("outputValues: %s" % param_output.outputValues)
205 self.endService("TELEMAC2D.Exec")
208 self._raiseSalomeError()
210 def InitBorders(self, bordersPickled):
212 This method initializes the border(s) for coupling
215 self.beginService("TELEMAC2D.InitBorders")
216 borders = cPickle.loads(bordersPickled)
218 for border in borders:
219 if border["model1"] is not None and border["model1"]["code"] == "TELEMAC2D":
220 self.borders[border["name"]] = border["model1"]
221 elif border["model2"] is not None and border["model2"]["code"] == "TELEMAC2D":
222 self.borders[border["name"]] = border["model2"]
223 self.endService("TELEMAC2D.InitBorders")
225 self._raiseSalomeError()
227 def ExecStep(self, timeDataPickled, inputDataPickled):
229 This method computes a single time step
232 self.beginService("TELEMAC2D.ExecStep")
233 timeData = cPickle.loads(timeDataPickled)
234 inputData = cPickle.loads(inputDataPickled)
236 if timeData["start_time"] == self.last_start_time:
237 # Convergence iteration, restore saved state
238 self.t2d.restore_state(self.saved_state)
240 # First iteration of a time step, save state
241 self.saved_state = self.t2d.save_state(self.saved_state)
243 self.last_start_time = timeData["start_time"]
245 if inputData is not None:
246 for (name, desc) in self.borders.iteritems():
247 if name in inputData:
251 if desc["position"] == "upstream":
252 if inputData[name]["Froude"] < 1.0: # Fluvial
253 z = inputData[name]["Z"]
255 v = inputData[name]["V"]
256 q = inputData[name]["Q"]
257 if inputData[name]["Froude"] > 1.0: # Torrential
258 z = inputData[name]["Z"]
259 self.set_z_on_border(desc, z)
260 #self.set_v_on_border(desc, v)
261 self.set_q_on_border(desc, q)
264 # TODO: Use time values from coupling
265 ierr = self.t2d.api.run_timestep_t2d(self.t2d.id)
267 raise Exception('Error in run_timestep_t2d:', ierr)
269 # Get the output values
271 for (name, desc) in self.borders.iteritems():
272 outputData[name] = {}
273 (outputData[name]["Z"], idx) = self.get_z_on_border(desc)
275 outputData[name]["V"] = 0.0
277 (outputData[name]["V"], froude) = self.get_v_at_point(idx)
278 outputData[name]["Q"] = self.get_q_on_border(desc)
279 outputDataPickled = cPickle.dumps(outputData, -1)
281 self.endService("TELEMAC2D.ExecStep")
282 return outputDataPickled
284 self._raiseSalomeError()
286 def get_z_on_border(self, border):
288 Find the water elevation (i.e. water height + bottom elevation) on a border.
289 Return a tuple (water_elevation, index_of_model_point_with_min_water_elevation).
290 If we use method 2 (compute mean water elevation), index is set to None.
292 # Method 1: Get point with minimum elevation
296 z_min = sys.float_info.max
297 (first, after_last) = border["points_id"]
299 while k != after_last:
300 nbor_k = self.t2d.get_integer('MODEL.NBOR', k)
301 h_k = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k)
302 z_bottom_k = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_k)
303 z_k = h_k + z_bottom_k
304 logger.debug("extract z for border point %d, h:%f, zf:%f, z:%f" % (k, h_k, z_bottom_k, z_k))
309 k = self.t2d.get_integer('MODEL.KP1BOR', k)
310 logger.debug("min z for border %d with method 1 is %f at border point %d (model point %d)" %
311 (border["border_id"], z_min, bor_idx, idx))
314 # Method 2: Compute mean elevation, ignoring dry zones (method from PALM coupling)
315 z_mean = sys.float_info.max
319 zf_min = sys.float_info.max
320 (first, after_last) = border["points_id"]
322 while k != after_last:
323 knext = self.t2d.get_integer('MODEL.KP1BOR', k)
324 nbor_k = self.t2d.get_integer('MODEL.NBOR', k)
325 nbor_knext = self.t2d.get_integer('MODEL.NBOR', knext)
326 zf1 = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_k)
327 zf2 = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_knext)
328 zf_min = min(zf_min, zf1)
329 h1 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k)
330 h2 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_knext)
333 x1 = self.t2d.get_double('MODEL.X', nbor_k)
334 x2 = self.t2d.get_double('MODEL.X', nbor_knext)
335 y1 = self.t2d.get_double('MODEL.Y', nbor_k)
336 y2 = self.t2d.get_double('MODEL.Y', nbor_knext)
337 l = math.sqrt((x1-x2)**2 + (y1-y2)**2)
338 if zf1 <= z_mean and zf2 <= z_mean:
339 h_tot += (z1+z2) * 0.5 * l
343 new_z_mean = h_tot / l_tot
344 delta = abs(z_mean - new_z_mean)
349 logger.debug("mean z for border %d with method 2 is %f" % (border["border_id"], z_mean))
355 def get_v_at_point(self, idx):
357 Find the water velocity at the model point idx.
358 Return a tuple (velocity, froude number)
360 u_i = self.t2d.get_double('MODEL.VELOCITYU', idx)
361 v_i = self.t2d.get_double('MODEL.VELOCITYV', idx)
362 h_i = self.t2d.get_double('MODEL.WATERDEPTH', idx)
363 v = math.sqrt(u_i**2 + v_i**2)
364 ##### TODO: Remove that and find the real error
366 print 'Error: water depth is negative or zero on point', idx
367 #raise Exception('Null water depth value')
369 froude = v / (h_i * 9.81)
372 def get_q_on_border(self, border):
374 Find the flow on a border.
375 The returned value is positive if the flow goes in the natural direction
376 according to the borders description (from upstream to downstream).
378 # Method 1: Integrate the flow on the whole border
381 (first, after_last) = border["points_id"]
383 while k != after_last:
384 knext = self.t2d.get_integer('MODEL.KP1BOR', k)
385 nbor_k = self.t2d.get_integer('MODEL.NBOR', k)
386 nbor_knext = self.t2d.get_integer('MODEL.NBOR', knext)
387 h_i = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k)
388 h_i2 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_knext)
389 h2d1 = 0.5 * (h_i + h_i2)
391 u_i = self.t2d.get_double('MODEL.VELOCITYU', nbor_k)
392 u_i2 = self.t2d.get_double('MODEL.VELOCITYU', nbor_knext)
393 v_i = self.t2d.get_double('MODEL.VELOCITYV', nbor_k)
394 v_i2 = self.t2d.get_double('MODEL.VELOCITYV', nbor_knext)
395 v2d1 = 0.5 * (math.sqrt(u_i**2+v_i**2) + math.sqrt(u_i2**2+v_i2**2))
396 x2d1 = self.t2d.get_double('MODEL.X', nbor_k)
397 x2d2 = self.t2d.get_double('MODEL.X', nbor_knext)
398 y2d1 = self.t2d.get_double('MODEL.Y', nbor_k)
399 y2d2 = self.t2d.get_double('MODEL.Y', nbor_knext)
400 q += v2d1 * h2d1 * math.sqrt((x2d1-x2d2)**2 + (y2d1-y2d2)**2)
402 logger.debug("flow for border %d with method 1 is %f" % (border["border_id"], q))
404 # Method 2: Get the flow from Telemac model
405 q = self.t2d.get_double('MODEL.FLUX_BOUNDARIES', border["border_id"])
406 if border["position"] == "downstream":
408 logger.debug("flow for border %d with method 2 is %f" % (border["border_id"], q))
411 # Method 3: Method from PALM coupling
413 (first, after_last) = border["points_id"]
415 while k != after_last:
416 knext = self.t2d.get_integer('MODEL.KP1BOR', k)
417 nbor_k = self.t2d.get_integer('MODEL.NBOR', k)
418 nbor_knext = self.t2d.get_integer('MODEL.NBOR', knext)
419 h_i = self.t2d.get_double('MODEL.WATERDEPTH', nbor_k)
420 h_i2 = self.t2d.get_double('MODEL.WATERDEPTH', nbor_knext)
421 u_i = self.t2d.get_double('MODEL.VELOCITYU', nbor_k)
422 u_i2 = self.t2d.get_double('MODEL.VELOCITYU', nbor_knext)
423 v_i = self.t2d.get_double('MODEL.VELOCITYV', nbor_k)
424 v_i2 = self.t2d.get_double('MODEL.VELOCITYV', nbor_knext)
425 x2d1 = self.t2d.get_double('MODEL.X', nbor_k)
426 x2d2 = self.t2d.get_double('MODEL.X', nbor_knext)
427 y2d1 = self.t2d.get_double('MODEL.Y', nbor_k)
428 y2d2 = self.t2d.get_double('MODEL.Y', nbor_knext)
433 q += ((h_i+h_i2)*(UN1+UN2)+h_i2*UN2+h_i*UN1)/6.0
435 if border["position"] == "upstream":
437 logger.debug("flow for border %d with method 3 is %f" % (border["border_id"], q))
441 def set_z_on_border(self, border, z):
443 Set the water elevation (i.e. water height + bottom elevation) on a border.
444 If z is None, then the water elevation for this border is free (not imposed).
446 # Method 1: Impose water elevation point by point
448 # Deactivate water level imposition in bord.f (requires a patch in bord.f)
449 self.t2d.set_double('MODEL.COTE', -1.0, border["border_id"])
451 (first, after_last) = border["points_id"]
453 while k != after_last:
455 self.t2d.set_integer('MODEL.LIHBOR', 4, k) # Free elevation on that border
457 self.t2d.set_integer('MODEL.LIHBOR', 5, k) # Imposed elevation on that border
458 nbor_k = self.t2d.get_integer('MODEL.NBOR', k)
459 bottom_elevation = self.t2d.get_double('MODEL.BOTTOMELEVATION', nbor_k)
460 hbor_k = max(0.0, z - bottom_elevation)
461 self.t2d.set_double('MODEL.HBOR', hbor_k, k)
462 logger.debug("set z for border point %d, h:%f, zf:%f, z:%f" % (k, hbor_k, bottom_elevation, z))
463 k = self.t2d.get_integer('MODEL.KP1BOR', k)
466 # Method 2: Just tell Telemac what is the imposed elevation for the border
468 self.t2d.set_double('MODEL.COTE', z, border["border_id"])
469 (first, after_last) = border["points_id"]
471 while k != after_last:
473 self.t2d.set_integer('MODEL.LIHBOR', 4, k) # Free elevation on that border
475 self.t2d.set_integer('MODEL.LIHBOR', 5, k) # Imposed elevation on that border
476 k = self.t2d.get_integer('MODEL.KP1BOR', k)
478 def set_v_on_border(self, border, v):
480 Set the water velocity on a border.
481 If v is None, then the water velocity for this border is free (not imposed).
483 # This method works only if the water level imposition is deactivated in bord.f
484 # (requires a patch in bord.f)
485 (first, after_last) = border["points_id"]
487 while k != after_last:
489 self.t2d.set_integer('MODEL.LIUBOR', 4, k) # Free velocity on that border
490 self.t2d.set_integer('MODEL.LIVBOR', 4, k)
492 self.t2d.set_integer('MODEL.LIUBOR', 6, k) # Imposed velocity on that border
493 self.t2d.set_integer('MODEL.LIVBOR', 6, k)
494 xnebor_k = self.t2d.get_double('MODEL.XNEBOR', k)
495 ubor_k = -xnebor_k * v
496 ynebor_k = self.t2d.get_double('MODEL.YNEBOR', k)
497 vbor_k = -ynebor_k * v
498 self.t2d.set_double('MODEL.UBOR', ubor_k, k)
499 self.t2d.set_double('MODEL.VBOR', vbor_k, k)
500 k = self.t2d.get_integer('MODEL.KP1BOR', k)
502 def set_q_on_border(self, border, q):
504 Set the water flow on a border.
505 If q is None, then the water flow for this border is free (not imposed).
508 self.t2d.set_double('MODEL.DEBIT', q, border["border_id"])
509 logger.debug("flow for border %d is set to %f" % (border["border_id"], q))
510 (first, after_last) = border["points_id"]
512 while k != after_last:
514 self.t2d.set_integer('MODEL.LIUBOR', 4, k) # Free velocity on that border
515 self.t2d.set_integer('MODEL.LIVBOR', 4, k)
517 self.t2d.set_integer('MODEL.LIUBOR', 5, k) # Imposed discharge on that border
518 self.t2d.set_integer('MODEL.LIVBOR', 5, k)
519 k = self.t2d.get_integer('MODEL.KP1BOR', k)
523 self.beginService("TELEMAC2D.Finalize")
527 self.endService("TELEMAC2D.Finalize")
529 self._raiseSalomeError()
531 def Compute(self, studyId, caseEntry):
533 self.beginService("TELEMAC2D.Compute")
535 jdc_dict = self.get_jdc_dict_from_case(studyId, caseEntry)
536 self.init_t2d_from_jdc_dict(jdc_dict)
539 self.t2d.run_all_timesteps()
544 self.endService("TELEMAC2D.Compute")
546 self._raiseSalomeError()
549 class Telemac2DOutputVariable:
551 def __init__(self, varstr):
552 par1 = varstr.find("(")
553 par2 = varstr.find(")")
554 self.varname = varstr[:par1]
555 self.idx = int(varstr[par1+1:par2])
557 def get_value(self, inst):
558 return inst.get_double(self.varname, self.idx)
561 class Telemac2DInputVariable:
563 def __init__(self, vardict):
564 self.varname = vardict["NOM_VARIABLE_MODELE"]
565 def_area_dict = vardict["ZONE_DE_DEFINITION"]
566 self.idx = def_area_dict.get("INDICE")
567 self.polygon = def_area_dict.get("POLYGONE")
568 self.points_in_polygon_list = None
570 def set_value(self, inst, value):
571 if self.idx is not None:
572 inst.set_double(self.varname, value, self.idx)
573 elif self.polygon is not None:
574 if self.points_in_polygon_list is None:
575 self.points_in_polygon_list = get_list_points_in_polygon(inst, self.polygon)
576 for idx in self.points_in_polygon_list:
577 inst.set_double(self.varname, value, idx)
579 raise Exception("Invalid area definition for input variable {0}".format(self.varname))