cnr-engineering / pyteltools Goto Github PK

In comp_ADCP_t2d the language is forced to "en". We should add a args.lang

Maillage.zip

Here is 3 files:

• "input.slf": File created by QGIS after mesh edition
• "pyteltools_out.slf": output file of below script
• "FUDAA_export.slf": export of pyteltools_out.slf with 1 timestep sith FUDAA

Python script

from pyteltools.slf import Serafin
with Serafin.Read("input.slf", "fr") as resin:
    resin.read_header()
    output_header = resin.header.copy()
    output_header.title = bytes("MAILLAGE", Serafin.SLF_EIT)
    output_header.empty_variables()
    output_header.add_variable_from_ID("B")
    output_header.add_variable_from_ID("W")
    output_values = np.ones(
        (output_header.nb_var, output_header.nb_nodes), dtype=output_header.np_float_type
    )

with Serafin.Write("pyteltools_out.slf", "fr" "", overwrite=True) as resout:
    resout.write_header(output_header)
    resout.write_entire_frame(output_header, 0.0, output_values)

Note np.ones

"Error" when reading pyteltools_out.slf with FUDAA

Frame 0 with 1 everywhere

Frame 1 (???) with 0

It seems harmless.

If "pyteltools_out.slf" is set as GEOMETRY FILE and computation is run in parallel, partel fails to split the file.

partel < partel_T2DGEO.par

 +-------------------------------------------------+
   PARTEL/PARRES: TELEMAC METISOLOGIC PARTITIONER

   REBEKKA KOPMANN & JACEK A. JANKOWSKI (BAW)
                  JEAN-MICHEL HERVOUET (LNHE)
                  CHRISTOPHE DENIS     (SINETICS)
                  YOANN AUDOUIN        (LNHE)
   PARTEL (C) COPYRIGHT 2000-2002
   BUNDESANSTALT FUER WASSERBAU, KARLSRUHE

   METIS 5.0.2 (C) COPYRIGHT 2012
   REGENTS OF THE UNIVERSITY OF MINNESOTA

   BIEF V8P5   (C) COPYRIGHT 2012 EDF
 +-------------------------------------------------+


   MAXIMUM NUMBER OF PARTITIONS:       100000

 +--------------------------------------------------+

 --INPUT FILE NAME <INPUT_NAME>:
 INPUT: T2DGEO
 --INPUT FILE FORMAT <INPFORMAT> [MED,SERAFIN,SERAFIND]:
  INPUT: SERAFIND
 --BOUNDARY CONDITIONS FILE NAME:
 INPUT: T2DCLI
--NUMBER OF PARTITIONS <NPARTS> [2 -100000]:
  INPUT:            2
  PARTITIONING METHOD <PMETHOD>                                                                                     [1 (METIS) OR 2 (SCOTCH)]:
 --INPUT:            1
 --CONTROL SECTIONS FILE NAME (OR RETURN) :
  NO SECTIONS
 --CONTROL ZONES FILE NAME (OR RETURN) :
  NO ZONES
 --WEIR FILE NAME (OR RETURN) :
  NO WEIRS
 --GEOMETRY FILE NAME <INPUT_NAME>:
 INPUT: T2DGEO
 --GEOMETRY FILE FORMAT <GEOFORMAT> [MED,SERAFIN,SERAFIND]:
  INPUT: SERAFIND
 --CONCATENATE FILES <YES-NO>:
 CONCATENATE: NO
 +---- PARTEL: BEGINNING -------------+
 FICHIER:T2DGEO                                                                                                                                                                        


 READ_MESH_INFO: TITLE= MAILLAGE
            NUMBER OF ELEMENTS:    26907
            NUMBER OF POINTS:    13628

            TYPE OF ELEMENT: TRIANGLE
            TYPE OF BND ELEMENT: POINT

            DOUBLE PRECISION FORMAT (R8)


 ONE-LEVEL MESH.
 NDP NODES PER ELEMENT:                    3
 ELEMENT TYPE :                           10
 NPOIN NUMBER OF MESH NODES:           13628
 NELEM NUMBER OF MESH ELEMENTS:        26907

 THE INPUT FILE ASSUMED TO BE 2D
 THERE ARE            1  TIME-DEPENDENT RECORDINGS

 THERE IS     3 ISLAND(S) IN THE DOMAIN

 THERE IS     2 LIQUID BOUNDARIES:

 BOUNDARY    1 :
  BEGINS AT BOUNDARY POINT:      125 , WITH GLOBAL NUMBER:         3
  AND COORDINATES:     350085.9           6247820.
  ENDS AT BOUNDARY POINT:      145 , WITH GLOBAL NUMBER:       319
  AND COORDINATES:     350253.9           6247411.

 BOUNDARY    2 :
  BEGINS AT BOUNDARY POINT:      237 , WITH GLOBAL NUMBER:      1017
  AND COORDINATES:     350444.2           6246046.
  ENDS AT BOUNDARY POINT:        1 , WITH GLOBAL NUMBER:      1218
  AND COORDINATES:     350472.9           6246000.

 THERE IS     5 SOLID BOUNDARIES:

 BOUNDARY    1 :
  BEGINS AT BOUNDARY POINT:        1 , WITH GLOBAL NUMBER:      1218
  AND COORDINATES:     350472.9           6246000.
  ENDS AT BOUNDARY POINT:      125 , WITH GLOBAL NUMBER:         3
  AND COORDINATES:     350085.9           6247820.

 BOUNDARY    2 :
  BEGINS AT BOUNDARY POINT:      145 , WITH GLOBAL NUMBER:       319
  AND COORDINATES:     350253.9           6247411.
  ENDS AT BOUNDARY POINT:      237 , WITH GLOBAL NUMBER:      1017
  AND COORDINATES:     350444.2           6246046.

 BOUNDARY    3 :
  BEGINS AT BOUNDARY POINT:      246 , WITH GLOBAL NUMBER:      1421
  AND COORDINATES:     350485.4           6246302.
  ENDS AT BOUNDARY POINT:      246 , WITH GLOBAL NUMBER:      1421
  AND COORDINATES:     350485.4           6246302.

 BOUNDARY    4 :
  BEGINS AT BOUNDARY POINT:      277 , WITH GLOBAL NUMBER:      3076
  AND COORDINATES:     350505.9           6246295.
  ENDS AT BOUNDARY POINT:      277 , WITH GLOBAL NUMBER:      3076
  AND COORDINATES:     350505.9           6246295.

 BOUNDARY    5 :
  BEGINS AT BOUNDARY POINT:      313 , WITH GLOBAL NUMBER:      6052
  AND COORDINATES:     350524.2           6246299.
  ENDS AT BOUNDARY POINT:      313 , WITH GLOBAL NUMBER:      6052
  AND COORDINATES:     350524.2           6246299.
  THE MESH PARTITIONING STEP STARTS
 BEGIN PARTITIONING WITH METIS
  RUNTIME OF METIS    0.00000000      SECONDS
  THE MESH PARTITIONING STEP HAS FINISHED
 ERROR IN BOUNDARIES J=         236
 WHILE NPTFR_P(I)=           0
 WHILE L=           0



 PLANTE: PROGRAM STOPPED AFTER AN ERROR
 RETURNING EXIT CODE:            2
STOP 1

With "FUDAA_export.slf", partel is OK

Hi Luc,

I'm using the duration over threshold in the GUI and what I'm finding is paradoxical.

I calculate the water depth's duration over 0.30 m in hours and then compute the water depth's duration over 0.60 m in hours and for the same node, the duration over the higher threshold (0.60 m) is longer than the duration spent over > 0.30 m. This seems physically impossible.

Am I misunderstanding the function of this tool perhaps?

Thanks,

When working with a big serafin file (10-15GB) the reading of the time using get_time() may take several seconds. Actually we need to read the whole file twice: ones to get time and then again to get 'real' data. So, it would be nice to add a possibility to compute the time instead of reading it from the file. The following code works nicely :)

    def get_time_fast(self):
        """!
        @brief Read the time from first two frames in the Serafin file and compute the time 
        using a time step between first and second frames
        """
        if self.header is None:
            raise SerafinRequestError('Cannot read time without any header (forgot read_header ?)')
        logger.debug('Reading the time series from the file')
        self.file.seek(self.header.header_size, 0)
        for n in range(self.header.nb_frames):
            if n < 2:
                self.file.read(4)
                self.time.append(self.header.unpack_float(self.file.read(self.header.float_size), 1)[0])
                self.file.read(4)
                self.file.seek(self.header.frame_size - 8 - self.header.float_size, 1)
            else:
                self.time.append(self.time[0] + n * (self.time[1] - self.time[0]))

Hello Luc, Using the workflow gui (and also via the source code) I see various references to geometry.Transformations. After having performed a simulation in geographic space, I'd like to perform some flux calculations and would like to first reproject the results into a UTM coordinate system.

The workflow GUI has a Transformation node but I'm not sure what the format (or how to create) the input file for it. Could you provide some more details there? Thank you.

When attempting to plot a variable (vertical cross section) whose min and max values are equal (e.g. both 0), PyTelTools crashes. Here is the traceback:

Traceback (most recent call last):
  File "D:\SOFTWARE\PyTelTools\src\pyteltools\PyTelTools\workflow\mono_gui.py", line 148, in mouseDoubleClickEvent
    node.configure()
  File "D:\SOFTWARE\PyTelTools\src\pyteltools\PyTelTools\workflow\nodes_vis.py", line 989, in configure
    if not self._prepare():
  File "D:\SOFTWARE\PyTelTools\src\pyteltools\PyTelTools\workflow\nodes_vis.py", line 1023, in _prepare
    self.plot_viewer.get_data(input_data, sections, line_interpolators_internal, section_indices)
  File "D:\SOFTWARE\PyTelTools\src\pyteltools\PyTelTools\workflow\util.py", line 1174, in get_data
    self.replot()
  File "D:\SOFTWARE\PyTelTools\src\pyteltools\PyTelTools\workflow\util.py", line 1142, in replot
    extend='both')
  File "D:\SOFTWARE\PyTelTools\lib\site-packages\matplotlib\axes\_axes.py", line 7677, in tricontourf
    return mtri.tricontourf(self, *args, **kwargs)
  File "D:\SOFTWARE\PyTelTools\lib\site-packages\matplotlib\tri\tricontour.py", line 283, in tricontourf
    return TriContourSet(ax, *args, **kwargs)
  File "D:\SOFTWARE\PyTelTools\lib\site-packages\matplotlib\tri\tricontour.py", line 40, in __init__
    ContourSet.__init__(self, ax, *args, **kwargs)
  File "D:\SOFTWARE\PyTelTools\lib\site-packages\matplotlib\contour.py", line 853, in __init__
    kwargs = self._process_args(*args, **kwargs)
  File "D:\SOFTWARE\PyTelTools\lib\site-packages\matplotlib\tri\tricontour.py", line 51, in _process_args
    tri, z = self._contour_args(args, kwargs)
  File "D:\SOFTWARE\PyTelTools\lib\site-packages\matplotlib\tri\tricontour.py", line 94, in _contour_args
    self._contour_level_args(z, args[1:])
  File "D:\SOFTWARE\PyTelTools\lib\site-packages\matplotlib\contour.py", line 1190, in _contour_level_args
    raise ValueError("Contour levels must be increasing")
ValueError: Contour levels must be increasing

I was able to work around this by adding some code to artificially increase the max value in replot(), but I'm sure there is a more elegant way of handling this.

It consists to move to routine in TatooineMesher here
You can assign it to me, it's on my todolist

Hey Luc,

I'm working on some 3D TELEMAC modeling applications and I would like to produce the file used to force the boundary. In particular the file this script produces.

http://docs.opentelemac.org/doxypydocs/v8p2r0/html/namespacepretel_1_1convert__to__bnd.html

I've already successfully interpolated 3D CMEMs temperature, salinity, and the general circulation to an extruded 2D unstructured grid using PyTelTools so essentially I just can grab the values at the liquid boundary and write those to this file format.

I'm in the process of trying to write a script that uses PyTelTools instead of the official Python TELEMAC API to write the file, however I'm not confident it will work because it seems to be mixing 2D/3D data structures in the same file. For instance, it seems to have a header with both ikle2 (2D triangulation connectivity) and ikle3 (3D triangulation connectivity). The IKLE2 has been modified to only represent boundary segments through some very confusing numpy logic.

Any advice? Do you think this could be possible to create this file using the existing code in PyTelTools or would it require additional data structures that would need to be added to PyTelTools?

Thank you,

Keith

Hi Luc,

I'm using the flux calculation node in PyTelTools. I'm curious what's the orientation convention for the positive and negative flux values?

For example is a positive flux going northwards and a negative flux passing southwards (assuming my feature is orientated east-west)?

(y-axis is m^3/s per m overflowing/overtopping rate).

Hi, I have configured the following workflow, the node Vertical Cross Section indicates the line is outside of the mesh but the node Locate Open Lines shows it inside. Any idea what may be incorrect or missing here? Thanks!

mesh = MeshInterpolator(slf.header, True)
gives
TypeError: () got an unexpected keyword argument 'unit'

mesh2d.py line 15
def init(self, input_header, construct_index=False, iter_pbar=lambda x: x):
probably should look like
def init(self, input_header, construct_index=False, iter_pbar=lambda x, unit: x):

Thanks

Hey Luc,

Would it be possible to save a workflow created in a GUI and then run it from the terminal?

This would be an awesome feature, if it doesn't exist yet.

Improve VTK export:

• add binary (in addition to ascii)
• support recent version of VTK file format

Paraview 5.5 sometimes indeed report this kind of error in parsing a VTK file generated by PyTelTools:

ERROR: In C:\bbd\7cc78367\build\superbuild\paraview\src\VTK\IO\Legacy\vtkDataReader.cxx, line 939
vtkUnstructuredGridReader (00000000128D60A0): Unsupported point attribute type: 1.31162e-41 for file: XXX.vtk

Hello,

Is there a API support to multiply two variables in one Serafin file and produce a new input file with the result?

Thanks

Hello Luc,

I would like to know if you have any scripts to interpolate a time-varying 3D result stored in a Serafin file onto a fixed vertical cross sectional plane?

Thanks for any help you can provide!

Keith

Add a layer selector to visualize scalars and vectors from a Serafin 3D input file.
Default layer selection should be top layer.

Add an efficient method for SerafinHeader to build IPOBO array from IKLE connectivity table.
This method is necessary to write a Serafin file from scratch.

• I have a mesh with a large amount of islands ~35k and I export it to a Serafin file for use with TELEMAC it takes hours
• The performance bottleneck is this loop when building ipobo.

  PyTelTools/pyteltools/slf/Serafin.py

  Line 623 in cd6cd48

  def build_ipobo(self):

Serafin.Write(output,'en',overwrite=True) as resout:
        output_header = Serafin.SerafinHeader(title='Converted from ADCIRC mesh format') 
        ikle = mesh.elements.triangles.to_numpy()
        xy = mesh.get_xy().to_numpy()
        output_header.from_triangulation(xy, ikle)

• Perhaps could be a option to pass your own IPOBO table to the Serafin header class constructor?

First of all, thanks for making this great project available. I had some difficulties installing it under Windows, but now it's working.
At first however, I was unable to load any Serafin files and the program would crash with a python error message regarding var_unit.decode('utf-8')
I found out that the problem was caused by the variable TEMPERATURE, whose unit is °C.
On a whim, I replaced all occurrences in the code of unit.decode('utf-8') with unit.decode('iso-8859-1'), and I can now load Serafin result files just fine.
Is this perhaps a Windows-specific problem?

Hi Luc,

I use PyTelTools in Qgis plugins and to avoid unintentionnal breaking behavior caused by new developpement, it would be nice if release could be created each time there is improvment or bug fix.

Thanks for the codes and hard work! Is there a script or a command that I can call from the Python interface to produce spatial maps of a given variable? I.e., visualize the free surface evolution in the entire domain?I looked through the codes and GUI but couldn't locate it.

I have the problem that my river is much flatter than it is wide. The ratio width to height is about 1:20. The square display of the section is therefore unfavorable for the evaluation. Can the height and width of the image be manually adjusted by the editor, e.g. 5 cm high and 15 cm wide?
Is the elevation and distance of the section itself changeable? The elevation and distance of my cross sections change as the flow progresses, so I would like to manually adjust them as well, e.g. the elevation of 96-108m to 96-102m.
before:

after:

Hello Luc,

I'm sorry to keep bothering you but I promise it's a good sign!

I'm using the code to compute volume in an enclosed polygonal area with the GUI.

If the grid and results are in geographic space (lat-lon), what are the units of volume?

Should I reproject the result into a UTM space first then do this calculation?

Thanks,

Keith

The temporal mean is not working (output Serafin file is empty).
It should be fixed in the workflow interface (Mean operator) and in the classic interface (Max/Min/Mean/Arrival/Duration).

Using a Serafin result file produced using Telemac v7p2: the Vertical Cross Section always displays only a single value (seemingly the initial value), see the screenshot. Also, when I select a different section, the x and y axis (elevation) are not adjusted to the new min-max values of that section.
Other visualizations (Vertical Temporal Profile 3D) and processes (Write vtk) work fine and show that there are indeed values in the Serafin file.
The problem does not occur using Serafin result files produced by Telemac v6p1, so I'm not sure what the issue is.

Hi Luc, I'm trying to execute the SynchMax node but it's force closing with this error message in the terminal. Any idea what I may be doing wrong? Thanks again for all your help, I really think this GUI can be useful there's just a learning curve.

The max operator works correctly. Another note: the node SynchMax takes almost no time to run.

Is there existing functionality in PyTelTools to depth-average a 3D Serafin file producing a 2D file? Thank you