slimPre.py 33.7 KB
Newer Older
1
2
3
4
5
6
7
"""Second-generation Louvain-la-Neuve Ice-ocean model - User interface

    Contact: jonathan.lambrechts@uclouvain.be
    webpage: www.climate.be/slim

    SLIM solves the governing equations on adaptative unstructured meshes using the discontinuous Galerkin finite element method.

8
    slimPre contains the tools writing the input files necessary for the simulations (for simulations, see slim.py file), such as partitioned mesh, bathymetry, coriolis, wind, tides, ...
9
10
11
12
    
    To generate the input files, run the script my_prepro_script.py:
        - On a single core: rundgpy my_prepro_script.py
        - On multiple cores: mpirun -np N rundgpy my_prepro_script.py (replace N by the number of cores)
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

    Mesh:
        A Mesh object contains a mesh, loaded from a .msh file, and the mesh projection system

    Region:
        A Region object contains a Mesh object, and a list of physical tags. The coordinates of all the mesh nodes contained on the physical tags is loaded in variable coordinates

    Coordinate_system:
        A Coordinate_system object contains a Region object, and the projection system of some data needed by the user. The variable coordinates gives the node coordinates of the Region in the data reference system.
        Vector field data can be rotate using the "rotate" and "rotate3D" functions

    Time:
        A Time object contains an array of times at which the data will be given. Time can be built using either a vector of instants relative to some initial time, or an initial and a final time and a time step. Some periodic window can also be built.

    write_file:
28
        The write_file function writes a netcdf file, using a Region object (giving the nodes for which there is data), a time object (giving the times of the data), and a list of data tuples.
29
30
        The first variable of a data tuple is the data name.
        The secod variable of a data tuple is the data array, with as many lines as the number of nodes, and as many columns as the number of time steps
31
32
33
    
    netcdf_to_msh:
        The function transforms a netcdf file to a .msh file.
34
35
36
37
38
39
40
41

    tpxo_tide:
        The function returns the data array of sea surface elevation, and the two horizontal transports (rotated in the coordinate system of the mesh, in m^2 s^-1) of the tides, using TPXO 7.2

    residual_flow:
        The function returns the data array of the two horizontal transports (in m^2 s^-1) at the nodes of some physical tags, such that the velocity is constant at all the nodes and that the total transport (in m^3 s^-1) is equal to the one prescribed by the user.

    smooth_bathymetry:
42
        The function smooths the given bathymetry, following a smoothing criterion, and writes it to a file. It can also write a bathymetry constant per element (P0), for a better stability. If the bathymetry is P0, it will be written in a idx file. Otherwise, it will be written in a netcdf file.
43
44
45
46
47
48
49
    
    interpolate_from_structured_grid
        The function interpolates from structured data 
    
    interpolate_from_structured_file
        The function interpolates structured data from a GMSH-file
    
50
51
52
53
    get_data_from_netcdf:
        The function returns the data array contained in the netcdf file given as an argument

    get_distance:
54
        The function returns an array with the minimum distance from each node of a given Region object to the physical tags given by the user. 
55
56

    fetch_ftp:
57
        This function downloads some data from the SLIM servers if available. The data is downloaded in the user global_data_directory. Data is only downloaded if it does not exist locally or is older than the server data. global_data_directory is by default ~/.cache/slim_data/
58
59
60
61
62

    set_global_data_directory:
        This functions gives the possibility to edit the global_data_directory

    make_directory:
63
        This function makes a directory if it does not exist
64
        
65
66
67
    partition_mesh:
        This function partitions the mesh (if run in multiprocessing) 
    
68
69
70
71
72
73
    partition_id:
        This function returns the number of the current partition as a string
    
    partition_nb:
        This function returns the number of partitions used as a string
    
74
    exit:
75
        This function closes the program
76
"""
Philippe Delandmeter's avatar
Philippe Delandmeter committed
77
import dgpy
78
from dgpy.scripts import slim_private
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
79
from dgpy.scripts import distance
Philippe Delandmeter's avatar
Philippe Delandmeter committed
80
81
82
import numpy as np

class Mesh :
83
    """Create the slimPre Mesh"""
84
    def __init__(self, mesh_file_name, mesh_proj=None) :
85
86
87
        """keyword arguments:
        
        * mesh_file_name            
88
            path to the mesh file (.msh format)
89
        * mesh_proj           
90
            mesh projection system (Default: None). 
91
92
            
            This argument must not be given if the mesh is written in 3D (on a sphere), or if all the data are in the same projection system than the mesh.
93
        """
94
95
96
97
98
99
        isPart = slim_private._check_if_partionned(mesh_file_name)
        if (dgpy.Msg.GetCommSize() > 1) and not isPart:
            dgpy.Msg.Fatal("Calling slimPre.Mesh in multiprocessing requires providing a partitioned mesh (first use slimPre.partition_mesh).")
        elif (dgpy.Msg.GetCommSize() == 1) and isPart:
            dgpy.Msg.Fatal("Calling slimPre.Mesh on a single core requires providing a non-partitioned mesh.")
        self._groups = dgpy.dgGroupCollection(mesh_file_name)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
100
        self._mesh_proj = mesh_proj
101
        
Philippe Delandmeter's avatar
Philippe Delandmeter committed
102

Philippe Delandmeter's avatar
Philippe Delandmeter committed
103
class Region :
104
105
106
107
    """Create the slimPre Region

        region.coordinates is an array of shape (n,3), n being the number of nodes in the region, with the coordinates of all the region nodes
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
108
    def __init__(self, mesh, physical_tags='') :
109
110
111
        """keyword arguments:
        
        * mesh
112
            the slimPre Mesh object
113
        * physical_tags
114
115
116
            the list of all the tags on which put some data (optional)
            
            If physical_tags are not set, all the mesh nodes are gathered
117
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
118
        self._mesh = mesh
119
120
        if slim_private._is_string(physical_tags):
            physical_tags = [physical_tags]
Philippe Delandmeter's avatar
Philippe Delandmeter committed
121
122
123
124
125
        self._physical_tags = physical_tags
        self._nodes = []
        self._groups = []
        self._node_map = {}
        self._node_ids = []
126
        self._empty = False
Philippe Delandmeter's avatar
Philippe Delandmeter committed
127
128
        num_ = 0
        self._type = "interface"
129
        for iGroup in range(mesh._groups.getNbElementGroups()+mesh._groups.getNbGhostGroups()) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
130
            group = mesh._groups.getElementGroup(iGroup)
131
            if ((self._physical_tags != ['']) and (not group.getPhysicalTag() in self._physical_tags)) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
                continue
            self._type = "volume"
            self._groups.append(group)
            nbNodes = group.getNbNodes()
            ig = len(self._groups) - 1
            for iElement in range(group.getNbElements()) :
                element = group.getElement(iElement)
                for iNode in range(nbNodes) :
                    nodeId = mesh._groups.mesh().gmshVertexId(element.vertex(iNode))
                    if not nodeId in self._node_map :
                        self._nodes.append((ig, nbNodes*iElement+iNode))
                        num_ += 1
                        self._node_map[nodeId] = num_
                        self._node_ids.append(nodeId)
        if self._type == "volume" :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
147
            self.coordinates = self._evaluateFunctor(dgpy.function.getCoordinates(), 3)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
148
149
150
151
            return

        for iFaceGroup in range(mesh._groups.getNbFaceGroups()) :
            faceGroup = mesh._groups.getFaceGroup(iFaceGroup)
152
            if (not faceGroup.physicalTag() in self._physical_tags) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
                continue
            self._groups.append(faceGroup)
            nbNodes = faceGroup.getNbNodes()
            ig = len(self._groups) - 1
            for iFace in range(faceGroup.size()) :
                group = faceGroup.elementGroup(0)
                iElement = faceGroup.elementId(iFace, 0)
                element = group.getElement(iElement)
                cl = faceGroup.closure(iFace, 0)
                for iNode in range(nbNodes) :
                    nodeId = mesh._groups.mesh().gmshVertexId(element.vertex(cl[iNode]))
                    if not nodeId in self._node_map :
                        self._nodes.append((ig, nbNodes*iFace+iNode))
                        num_ += 1
                        self._node_map[nodeId] = num_
                        self._node_ids.append(nodeId)
        self.coordinates = self._evaluateFunctor(dgpy.function.getCoordinates(), 3)
170
171
        if len(self._nodes) == 0:
            self._empty = True
Philippe Delandmeter's avatar
Philippe Delandmeter committed
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188

    def _evaluateFunctor(self, f, size) :
        fc = dgpy.functorCache(dgpy.functorCache.NODE_GROUP_MODE, self._mesh._groups)
        val = np.empty((len(self._nodes), size))
        valg = []
        for g in self._groups :
            if self._type == "volume":
                fc.setGroup(g)
            else:
                fc.setInterfaceGroup(g)
            valg.append(np.copy(fc.get(f)))
        for i,(iG, iN) in enumerate(self._nodes) :
            val[i,:] = valg[iG][iN]
        return val


class Coordinate_system :
189
190
191
192
    """Create the slimPre Coordinate_system

        Coordinate_system.coordinates is an array of shape (n,3), n being the number of nodes in the region, with the coordinates of all the region nodes in the data projection system.
        Coordinate_system.rotate(u,v) returns the rotated vector array (u,v)
193
        Coordinate_system.rotate3D ... 
194
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
195
    def __init__(self, region, data_proj) :
196
197
198
        """keyword arguments:
        
        * region
199
            the slimPre Region object
200
        * data_proj
201
202
            data projection system.

203
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
204
        self._region = region
205
206
        if region._empty:
            return
207
        slim_private._print('Projection system for tags:')
208
        if region._physical_tags == '':
209
            slim_private._print('    the entire mesh')
210
211
        else:
            if(slim_private._is_string(region._physical_tags)):
212
                slim_private._print('    '+region._physical_tags)
213
214
            else:
                for tag in region._physical_tags:
215
                    slim_private._print('    '+tag)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
216
        mesh_proj = region._mesh._mesh_proj
Philippe Delandmeter's avatar
Philippe Delandmeter committed
217
        if mesh_proj is not None :
218
219
220
            slim_private._print('  Projection used:')
            slim_private._print('    mesh proj: %s' % mesh_proj)
            slim_private._print('    data proj: %s' % data_proj)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
221
222
223
            _mesh_proj = dgpy.pj_init_plus(mesh_proj)
            _data_proj = dgpy.pj_init_plus(data_proj)
        elif mesh_proj is None :
224
225
226
            slim_private._print('  Projection used:')
            slim_private._print('    mesh proj: no projection. Coordinates are 3D')
            slim_private._print('    data proj: %s' % data_proj)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
227
228
229
230
            mesh_proj2 = "+proj=latlong +ellps=WGS84"
            _mesh_proj = dgpy.pj_init_plus(mesh_proj2)
            _data_proj = dgpy.pj_init_plus(data_proj)

231
        xyz = region.coordinates      
Philippe Delandmeter's avatar
Philippe Delandmeter committed
232
        if mesh_proj :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
233
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, xyz)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
234
235
236
237
238
239
240
241
242
243
244
245
            epsVect = np.zeros(xyz.shape) 
            epsVect[:,0] = np.hypot(xyz[1,0]-xyz[0,0],xyz[1,1]-xyz[0,1]) / 1000.
            xyz0 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz-epsVect)
            xyz1 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz+epsVect)
            if dgpy.pj_is_latlong(_data_proj) :
                xyz1[:,0] = np.where( xyz1[:,0]<xyz0[:,0], xyz1[:,0]+2*np.pi, xyz1[:,0])
            self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])
        else :
            lonlat = np.zeros(xyz.shape)
            lonlat[:,0] = np.arctan2(xyz[:,1], xyz[:,0])
            lonlat[:,1] = np.arcsin(xyz[:,2] / np.sqrt(xyz[:,0]*xyz[:,0]+xyz[:,1]*xyz[:,1]+xyz[:,2]*xyz[:,2]) )
            lonlat[:,2] = np.sqrt(xyz[:,0]*xyz[:,0]+xyz[:,1]*xyz[:,1]+xyz[:,2]*xyz[:,2]) 
Philippe Delandmeter's avatar
Philippe Delandmeter committed
246
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, lonlat)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
247
248
249
250
            epsVect = np.zeros(xyz.shape) 
            epsVect[:,0] = 0.01
            xyz0 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz-epsVect)
            xyz1 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz+epsVect)
251
            self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
252
253
254
255
256
257
258

        if _mesh_proj:
            dgpy.pj_free(_mesh_proj)
        if _data_proj:
            dgpy.pj_free(_data_proj)

    def rotate(self, u, v):
259
        """Return the rotated vector data array (u,v). Returned vector will be oriented in the mesh projection system
260
261
262
263
264
265
        
        keyword arguments:
        
        * u
            first component of the vector (in the data coordinate system)
        * v
266
            second component of the vector (in the data coordinate system)
267
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
268
269
270
        u2 =  u * np.cos(self._angle) + v * np.sin(self._angle)
        v2 = -u * np.sin(self._angle) + v * np.cos(self._angle)
        return (u2, v2)
271
        
272
    def rotate3D(self, F_lon, F_colat, F_r):
273
274
275
276
277
278
279
280
281
282
283
        """Return ...
        
        keyword arguments:
        
        * F_lon
            this var ...
        * F_colat
            this var ... 
        * F_r
            this var ...
        """
284
285
        lon = self.coordinates[0]
        colat = 0.5*np.pi - self.coordinates[1]
286
287
288
289
        Fx = np.sin(colat)*np.cos(lon)*F_r + np.cos(colat)*np.cos(lon)*F_colat - np.sin(lon)*F_lon
        Fy = np.sin(colat)*np.sin(lon)*F_r + np.cos(colat)*np.sin(lon)*F_colat + np.cos(colat)*F_lon
        Fz = np.cos(colat)*F_r             - np.sin(colat)*F_colat
        return (Fx,Fy,Fz)    
Philippe Delandmeter's avatar
Philippe Delandmeter committed
290
291
292


class Time :
293
294
    """Create the slimPre Time"""

295
    def __init__(self, time_vector=None, initial_time='1970-01-01 00:00:00', final_time=None, time_step=None, periodic=False) :
296
297
298
        """keyword arguments:

        * time_vector
299
            vector with all the time steps. Time steps are given in seconds from initial_time. time_vector cannot be given if final_time and/or time_step are set
300
        * initial_time
301
            initial_time given in the format 'YYYY-MM-DD hh:mm:ss'. If not set, initial_time is '1970-01-01 00:00:00'
302
        * final_time
303
            final_time given in the format 'YYYY-MM-DD hh:mm:ss'. final_time cannot be given if the time_vector is set. If final_time is set, time_step must be given.
304
        * time_step
305
306
307
308
309
            time_step in seconds. time_step cannot be given if the time_vector is set. If time_step is set, final_time must be given.
        * periodic
            flag to define a periodic variable (Default: False)
            
            If periodic, the first time should be 0 and the last time should be the period (e.g. time_vector should be from 0 to T), and the first value is repeated.
310
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
311
312
        self._periodic = periodic
        fmt = '%Y-%m-%d %H:%M:%S'
313
314
315
        if not final_time and not time_step :
            self._time = np.array(time_vector)
            fmt = '%Y-%m-%d %H:%M:%S'
316
317
            date_ref  = slim_private.datetime.datetime.strptime('1970-01-01 00:00:00', fmt)
            date_user = slim_private.datetime.datetime.strptime(initial_time, fmt)
318
319
320
            date_shift = (date_user-date_ref)
            time_shift = date_shift.days * 86400 + date_shift.seconds
            self._time += time_shift
321
        elif not time_vector and initial_time and final_time and time_step:
322
323
324
325
            date  = slim_private.datetime.datetime.strptime(initial_time, fmt)
            initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
            date  = slim_private.datetime.datetime.strptime(final_time, fmt)
            final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
326
            self._time = np.arange(initial_time, final_time+time_step*0.9999, time_step)
327
328
        else:
            dgpy.Msg.Fatal("Wrong combination of inputs to create a Time object")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
329
330


Philippe Delandmeter's avatar
Philippe Delandmeter committed
331
def write_file(file_name, region=None, time=None, data=[]):
332
    """Create a SLIM netcdf file, standard input for SLIM simulations.
333
334
335
336
337
338

    keyword arguments:

    * file_name
        name of the netcdf file to write
    * region
339
        slimPre Region object (optional).
340
    * time
341
        slimPre Time object (optional).
342
    * data
343
344
        data is a list of data tuples. Each tuple follows the form 
        
345
        ("data_name", data_array)
346
        
347
348
        where data_array is an array of n lines and m columns, n being the number of nodes in the region and m being the number of time steps in the time.
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
349
    oF = dgpy.slimNetCDFWrap(file_name)
350
351
352
    space_dep = region is not None
    time_dep = time is not None
    if space_dep :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
353
        nbNodes = len(region._node_ids)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
354
355
        nodeDim = oF.create_dimension('node', nbNodes)
        nodeVar = oF.create_variable_int('node', [nodeDim])
356
    if time_dep :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
357
358
359
360
361
        nTimes = np.shape(time._time)[0]
        timeArray = np.array([time._time])
        if not time._periodic:
            timeUnits = 'seconds since 1970-01-01 00:00:00'
        else:
362
            timeUnits = 'seconds since the beginning of the period'
363
364
        timeDim = oF.create_dimension('time', nTimes)
        timeVar = oF.create_variable('time', [timeDim], timeUnits, time._periodic)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
365
366
367
    nFields = len(data)
    valVar = [0]*nFields
    for iField in range(nFields) :
368
        if space_dep and time_dep:
369
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim, nodeDim], None, False)
370
371
372
        elif space_dep and (not time_dep):
            valVar[iField] = oF.create_variable(data[iField][0], [nodeDim], None, False)
        elif (not space_dep) and time_dep:
373
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim], None, False)
374
375
        else:
            valVar[iField] = oF.create_variable(data[iField][0], [], None, False)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
376
    oF.enddef()
377
378
379
380
    if space_dep :
        oF.set_variable_int(nodeVar, region._node_ids)        
    if time_dep :
        oF.set_variable(timeVar, timeArray)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
381
382
383
384
385
386
387
388
389
390
391
392
393
394
    for iField in range(nFields) :
        if not isinstance(data[iField][1], np.ndarray) :
            try :
                values = np.array([float(data[iField][1])])
            except :
                values = np.array(data[iField][1])
        else :
            values = data[iField][1]
        if len(values.shape) == 1:
            values = values[None,:]

        oF.set_variable(valVar[iField], values)


395
def netcdf_to_msh(mesh_file_name, nc_file_name, variable_name, export_file_name, time = 0):
396
    """Export a field in the netcdf format to msh format.
Philippe Delandmeter's avatar
Philippe Delandmeter committed
397
398

    keyword arguments:
399
400
401
402
403
404
405
406
407
408
409
    
    * mesh_file_name
        path to the mesh file (.msh format)
    * nc_file_name    
        path to the netcdf file (.nc format)
    * variable_name 
        name of the variable in the netcdf file
    * export_file_name 
        path to the field name (.msh format)
    * time  
        time at which the data will be exported in msh format (default: 0)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
410
411
    """
    groups = dgpy.dgGroupCollection(mesh_file_name)
412
    f = slim_private._load_function((nc_file_name,variable_name),groups)
413
    dgpy.function.getTime().set(time)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
414
415
    groups.exportFunctionMsh(f, export_file_name)

Philippe Delandmeter's avatar
Philippe Delandmeter committed
416
def tpxo_tide(region, time, data_file_name=None): 
417
418
419
420
    """Give TPXO tides data
    
    This function will download the files tpxo and generate the data at the requested time steps.
    It will return the tpxo tides as variables (ssh, ssu, ssv)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
421
422
    
    keyword arguments:
423
    
424
    * region
425
        slimPre Region object on which put the data
426
    * time
427
        Time object of the instant at which put the data
428
    * data_file_name
429
        name of the netcdf data file where tides are exported (optional)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
430
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
431
    lonlat = Coordinate_system(region, data_proj='+proj=latlong +ellps=WGS84')
432
    h = u = v = np.empty((1,1))
Philippe Delandmeter's avatar
Philippe Delandmeter committed
433
434
    h_file_name = 'h_tpxo7.2.nc'
    u_file_name = 'u_tpxo7.2.nc'
435
436
    h_file_name = slim_private.dgftp.get( ("/slim_data/tides/%s" % h_file_name) )
    u_file_name = slim_private.dgftp.get( ("/slim_data/tides/%s" % u_file_name) )
Philippe Delandmeter's avatar
Philippe Delandmeter committed
437
438
439
    tide_ssh = dgpy.slimFunctionTpxo(h_file_name,'ha','hp','lon_z','lat_z')
    tide_u = dgpy.slimFunctionTpxo(u_file_name,'Ua','up','lon_u','lat_u')
    tide_v = dgpy.slimFunctionTpxo(u_file_name,'Va','vp','lon_v','lat_v')
440
441
442
443
    if not region._empty:
        lonlat_coordinates = lonlat.coordinates * 180 / np.pi
        n = np.size(time._time)
        m = len(lonlat_coordinates)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
444

445
446
447
448
449
450
451
        h = np.empty((n,m))
        u = np.empty((n,m))
        v = np.empty((n,m))
        for i in range(n):
            h[i,:] = tide_ssh.getAtPoints(lonlat_coordinates, time._time[i])[:,0]
            u[i,:] = tide_u.getAtPoints(lonlat_coordinates, time._time[i])[:,0]
            v[i,:] = tide_v.getAtPoints(lonlat_coordinates, time._time[i])[:,0]
Philippe Delandmeter's avatar
Philippe Delandmeter committed
452

453
        (u, v) = lonlat.rotate(u,v)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
454
455

    if data_file_name is not None:
Philippe Delandmeter's avatar
Philippe Delandmeter committed
456
        write_file(data_file_name, region=region, time=time, data=[('h', h), ('u', u), ('v', v)])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
457
458
459
460


    return (h, u, v)

Philippe Delandmeter's avatar
Philippe Delandmeter committed
461
def residual_flow(region, flow, bathymetry, data_file_name=None): 
462
463
464
    """Set a residual flow
    
    It will return the velocity (u, v) such that speeds at all the nodes are constant and the velocity is normal to the boundary. The total transport (in m^3 s^-1) is equal to the given flow.
465
466
467
    The physical tags of the region object must be boundary physical tags.
    
    keyword arguments:
468
    
469
    * region
470
        slimPre Region object on which put the data
471
472
473
474
475
    * flow
        value of the flow crossing the region (must be a boundary region), in m^3 s^-1
    * bathymetry
        bathymetry in netcdf format or msh format. It is necessary to compute the transport
    * data_file_name
476
        name of the netcdf data file where the residual flow is exported (optional)
477
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
478

Philippe Delandmeter's avatar
Philippe Delandmeter committed
479
    groups = region._mesh._groups
Philippe Delandmeter's avatar
Philippe Delandmeter committed
480
481
482
    integNormals = dgpy.dgFunctionIntegratorInterface(groups, dgpy.function.getNormals())
    nx = 0
    ny = 0
Philippe Delandmeter's avatar
Philippe Delandmeter committed
483
    for tag in region._physical_tags:
484
485
486
487
488
        #normals = integNormals.computeVector(tag)
        normals = dgpy.fullMatrixDouble(1,3,True)
        integNormals.compute(tag, normals)
        nx += normals.get(0,0)
        ny += normals.get(0,1)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
489
490
491
492
493
494
495
496
497
498
499
    norm = np.sqrt(nx*nx+ny*ny)
    nx = nx/norm
    ny = ny/norm

    bathFunc = slim_private._load_function(bathymetry, groups)
    def orientedSurfNumpy(cmap, val, h, n) : 
        val[:] = h[:] * (n[:,0]*nx + n[:,1]*ny)
    f = dgpy.functionNumpy(1, orientedSurfNumpy, [bathFunc, dgpy.function.getNormals()])

    integInterface = dgpy.dgFunctionIntegratorInterface(groups, f)
    section = 0
500
    sectionFM = dgpy.fullMatrixDouble(1,1,True)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
501
    for tag in region._physical_tags:
502
503
504
        #section += integInterface.compute(tag)
        integInterface.compute(tag, sectionFM)
        section += sectionFM.get(0,0)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
505
506
507
508
509

    def transportNumpy(cmap, val, h) : # eta is neglected here
        val[:,0] = -h[:] * flow / section * nx
        val[:,1] = -h[:] * flow / section * ny
    f = dgpy.functionNumpy(2, transportNumpy, [bathFunc])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
510
    vals = region._evaluateFunctor(f, 2)
511
512
    if region._empty:
        vals = np.empty((1,2))
513
    if data_file_name:
514
        write_file(data_file_name, region, None, [("uflow", vals[:,0]), ("vflow", vals[:,1])])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
515
516
    return vals

517
def smooth_bathymetry(mesh, bathymetry, output_file_name, coefficient=0.5, transform_p0=True):
518
    """Smooth the bathymetry in order to locally reduce the noise in the momentum
Philippe Delandmeter's avatar
Philippe Delandmeter committed
519
520
    
    keyword arguments:
521
522
523
    
    * mesh 
        slimPre mesh object
524
525
    * bathymetry
        bathymetry file
526
        
527
        Format:
528
529
530
531
            * tuple containing the raw bathymetry data file (.nc format) and the name of the variable bathymetry 
            * bathymetry mesh file (.msh or .idx) 
    * output_file_name
        name of the file name containing the smooth bathymetry. The format is netcdf (if transform_p0 = False) or idx (if transform_p0 = True).
532
533
534
    * coefficient
        coefficient for the smoothing algorithm (default: 0.5). This forbids the bathymetry from moving from more than 100*coefficient % on a single element
    * transform_p0
535
        flag to transform the bathymetry to a bathymetry constant per element. This gives a good stability and has negligible impact the results (default: True)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
536
537
538
    """
    groups = mesh._groups
    dof = dgpy.dgDofContainer(groups, 1)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
539
    slim_private._load(dof, bathymetry)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
540
541
542
543
544
545
546
547
548
549
    #load data
    m = groups.cgMesh()
    n = m.nVertices()
    data = np.ndarray((n),"d")
    for ig in range(groups.getNbElementGroups()) :
        g = groups.getElementGroup(ig)
        gdata = dof.getGroupProxy(ig)
        for ie in range(g.getNbElements()):
            vs = m.vertices(ig, ie)
            for i, v in enumerate(vs) :
550
551
552
                if v >= n:
                    n *= 2
                    data.resize((n))
Philippe Delandmeter's avatar
Philippe Delandmeter committed
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
                data[v] = gdata(i,ie)
    #smoothing
    count = 0
    factor = coefficient
    while True:
        count = count+1
        emax = 0
        for ig in range(groups.getNbElementGroups()):
            g = groups.getElementGroup(ig)
            for ie in range(g.getNbElements()):
                vs = m.vertices(ig, ie)
                vs = np.array(vs)
                bathe = data[vs]
                dmax = bathe.max() * factor
                for i in range(len(vs)):
                    for j in range(i+1, len(vs)):
                        d = abs(bathe[i] - bathe[j])
                        if (d > dmax):
                            mean = (bathe[i] + bathe[j])/2
                            if (d - dmax > emax):
                                emax = d - dmax
                            if (bathe[i] > bathe[j]):
                                 bathe[i] = mean+dmax/2
                                 bathe[j] = mean-dmax/2
                            else:
                                 bathe[i] = mean-dmax/2
                                 bathe[j] = mean+dmax/2
                data[vs,] = bathe
581
582
        #TODO better print
        print(count,emax)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
583
584
        if emax < .2 :
            break
585
    #read data
Philippe Delandmeter's avatar
Philippe Delandmeter committed
586
587
    for ig in range(groups.getNbElementGroups()) :
        g = groups.getElementGroup(ig)
588
        gdata = dof.getGroupProxy(ig)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
589
590
591
        for ie in range(g.getNbElements()):
            vs = m.vertices(ig, ie)
            for i, v in enumerate(vs) :
592
                gdata.set(i,ie, data[v])
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
    
    if (not transform_p0):
        region = Region(mesh)
        data = region._evaluateFunctor(dof.getFunction(), 1)
        if not slim_private._is_string(bathymetry):
            if bathymetry[0][-3:] == ".nc":
                write_file(output_file_name, region=region, time=None, data=[(bathymetry[1], data)])
        else:
            write_file(output_file_name, region=region, time=None, data=[('smoothed_bath', data)])
    else:
        nodalVolume = dgpy.dgNodalVolume(groups)
        converter = dgpy.dgPrismDofConverter(nodalVolume.get())
        bath_P0 = dgpy.dgDofContainer(groups,1)
        converter.convertToPZero2d(dof,bath_P0)
        if output_file_name[-3:] == ".nc":
608
            slim_private._print("The bathymetry is dg and saved in "+output_file_name[:-3]+"/"+slim_private.path.basename(output_file_name[:-3])+".idx")
609
610
611
612
            output_file_name = output_file_name[:-3]
        exporter = dgpy.dgIdxExporter(bath_P0, output_file_name)
        exporter.exportIdx()
        
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
def interpolate_from_structured_grid(x, y, ox, oy, dx, dy, data, fill_value=-9999.):
    """Interpolates the values from a structured grid

    keyword arguments:
    
    * x
        list of the x-coordinates of all the points to interpolate at
    * y
        list of the y-coordinates of all the points to interpolate at
    * ox
        x-origin of the structured data
    * oy
        y-origin of the structured data
    * dx
        x-step of the structured data
    * dy
        y-step of the structured data
    * data
        structured grid
    * fill_value
        default value given to a point outside the structured grid
    """    
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
    ny = data.shape[0]
    nx = data.shape[1]
    rx = x - ox
    ry = y - oy
    ix = (rx / dx).astype(int)
    iy = (ry / dy).astype(int)
    index=np.empty(x.shape)
    if (dx > 0):
        index = (x < ox) + (x > ox + (nx-1)*dx)
    else:
        index = (x > ox) + (x < ox + (nx-1)*dx)
    if (dy > 0):
        index+= (y < oy) + (y > oy + (ny-1)*dy)
    else:
        index+= (y > oy) + (y < oy + (ny-1)*dy)
    ix = np.minimum(np.maximum(ix,0),nx-1)
    iy = np.minimum(np.maximum(iy,0),ny-1)
    xsi = (rx - ix*dx) / dx
    eta = (ry - iy*dy) / dy
654
    val = np.where(index, fill_value, \
655
656
657
658
659
660
          (1-xsi)*(1-eta) * data[ iy      ,  ix      ] + \
          (  xsi)*(1-eta) * data[ iy      , (ix+1)%nx] + \
          (  xsi)*(  eta) * data[(iy+1)%ny, (ix+1)%nx] + \
          (1-xsi)*(  eta) * data[(iy+1)%ny,  ix      ])
    return val
        
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
def interpolate_from_structured_file(nodes, structured_file, fill_value=-9999, transpose_matrix=False):
    """Interpolates the values from a structured file of gmsh type:
    
    ox oy oz
    
    dx dy dz
    
    nx ny nz
    
    val11 val12 ...
    
    val21 val22 ...
    
    ...

676
677
    This module is only available in 2D !

678
679
680
681
682
683
684
685
686
687
688
    keyword arguments:
    
    * nodes
        coordinates of the nodes to interpolate at (from Region or Coordinate_system)
    * structured_file
        name of the file containing the structured data (binary files should end with '.bin' extension) 
    * fill_value
        default value given to a point outside the structured grid
    * transpose_matrix
        boolean assessing if the matrix should be transposed (Default: False) (True if nx is the number of lines and ny is the number of columns)
    """   
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
    if structured_file[-4:] == '.bin':
        import struct
        f = open(structured_file,'rb')
        O = struct.unpack("ddd", f.read(24))
        ox = O[0]
        oy = O[1]
        D = struct.unpack("ddd", f.read(24))
        dx = D[0]
        dy = D[1]
        N = struct.unpack("iii", f.read(12))
        nx = N[0]
        ny = N[1]
        if transpose_matrix:
            data = np.fromfile(f, "d", nx * ny).reshape((nx, ny))
        else:
            data = np.fromfile(f, "d", nx * ny).reshape((ny, nx))
    else:
        f = open(structured_file,'r')
        line = f.readline()
        origin = line.split(' ')
        ox = float(origin[0])
        oy = float(origin[1])
        line = f.readline()
        delta = line.split(' ')
        dx = float(delta[0])
        dy = float(delta[1])
        line = f.readline()
        number = line.split(' ')
        nx = int(number[0])
        ny = int(number[1])
        if transpose_matrix:
            data = np.fromfile(f, "d", nx * ny, sep=' ').reshape((nx, ny))
        else:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
722
            data = np.fromfile(f, "d", nx * ny, sep=' ').reshape((ny, nx))
723
724
725
726
    f.close()
    x = nodes[:,0]
    y = nodes[:,1]
    data = np.transpose(data)
727
    val = interpolate_from_structured_grid(x, y, ox, oy, dx, dy, data, fill_value=fill_value)    
728
729
    return val
        
Philippe Delandmeter's avatar
Philippe Delandmeter committed
730
def get_data_from_netcdf(nc_file_name, variable_name):
731
    """Return the array containing the values of the variable in the nectdf file (SLIM format)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
732
733

    keyword arguments:
734
    
735
736
737
738
    * nc_file_name
        path to the netcdf file (.nc format)
    * variable_name
        name of the variable in the netcdf file
Philippe Delandmeter's avatar
Philippe Delandmeter committed
739
740
    """
    #TODO use a wrapper instead of netCDF4
741
    slim_private._print("Loading netCDF4 for python, might not be installed on your computer...")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
742
    from netCDF4 import Dataset
743
    f = Dataset(nc_file_name,"r")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
744
745
746
    val = f.variables[variable_name][:]
    return val

747
def get_distance(region, physical_tags):
748
    """Return the minimum distance from each node on the region to the physical tags. The index of the 1-D array returned is the nodal index, same as Region.coordinates 
749
750
    
    keyword arguments:
751
    
752
    * region
753
        slimPre Region object of the nodes on which compute the distance
754
755
    * physical_tags
        vector containing the physical tags from which the minimum distance will be computed
756
757
    """
    xyz = region.coordinates
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
758
    dist = distance.distance(region._mesh, physical_tags, evaluationPoints=xyz)
759
760
    return np.array(dist)

761
def fetch_ftp(file_name, data_dir=None):
762
763
    """This function downloads some data from the SLIM servers if available. The data is downloaded in the user global_data_directory. Data is only downloaded if it does not exist locally or is older than the server data. 

764
765
    global_data_directory is by default ~/.cache/slim_data/
    """
766
    file_name =  slim_private.dgftp.get(file_name)
767
768
    if data_dir: 
        if dgpy.Msg.GetCommRank() == 0:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
769
            slim_private.shutil.copyfile(file_name, data_dir+slim_private.path.basename(file_name))
770
            file_name = data_dir+slim_private.path.basename(file_name)
771
772
773
        dgpy.Msg.Barrier()
    return file_name

774
def set_global_data_directory(path):
775
    """Edit the global_data_directory"""
776
777
    slim_private.dgftp.global_data_local_directory = path

778
def make_directory(data_dir):
779
    """Create data_dir directory if it does not exist yet"""
780
781
    if not slim_private.path.isdir(data_dir) :
        slim_private.makedirs(data_dir)
782
    dgpy.Msg.Barrier()
783
    
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
def partition_mesh(mesh_file_name, partitioned_mesh_file_name):
    """Partition the mesh 
    
    keyword arguments:
    
    * mesh_file_name            
            path to the mesh file (.msh format)
    * partitioned_mesh_file_name
            name where the partitioned mesh is stored.
    """
    isPart = slim_private._check_if_partionned(mesh_file_name)
    if isPart:
        dgpy.Msg.Warning("The mesh "+mesh_file_name+" is already partitioned !")
        slim_private.shutil.copyfile(mesh_file_name, partitioned_mesh_file_name)
    elif dgpy.Msg.GetCommSize() == 1:
        dgpy.Msg.Info("No partition done on a single processor, the partition will be done later, if needed")
        slim_private.shutil.copyfile(mesh_file_name, partitioned_mesh_file_name)
    else:
        dgpy.dgMeshPartition(mesh_file_name, dgpy.Msg.GetCommSize())
        mesh_name = mesh_file_name[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh'
        slim_private.shutil.copyfile(mesh_name, partitioned_mesh_file_name)          
        dgpy.Msg.Barrier()
        if dgpy.Msg.GetCommRank() == 0 :
            slim_private.remove(mesh_name)
    return partitioned_mesh_file_name

810
811
812
813
814
815
816
817
818
def partition_id():
    """Return the number of the current partition as a string
    """
    return str(dgpy.Msg.GetCommRank())

def partition_nb():
    """Return the number of partitions used as a string
    """
    return str(dgpy.Msg.GetCommSize())
819
820

def exit(i):
821
    """Exit properly the program"""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
822
    dgpy.Msg.Exit(i)