slimPre.py 30.8 KB
Newer Older
1
2
3
4
5
6
7
8
"""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.

    slimPre contains the tools writing the input files necessary for the simulations (for simulations, see slim.py file), such as partitionned 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

    get_data_from_netcdf:
        The function returns the data array contained in the netcdf file given as an argument

    get_distance:
48
        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. 
49
50

    fetch_ftp:
51
        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/
52
53
54
55
56
57
58
59
60

    set_global_data_directory:
        This functions gives the possibility to edit the global_data_directory

    make_directory:
        The function makes a directory if it does not exist
    exit:
        The function closes the program
"""
Philippe Delandmeter's avatar
Philippe Delandmeter committed
61
62
63
64
65
import dgpy
import slim_private
import numpy as np

class Mesh :
66
    """Create the slimPre Mesh"""
Philippe Delandmeter's avatar
Philippe Delandmeter committed
67
    def __init__(self, mesh_file_name, mesh_proj=None) :
68
69
70
71
72
        """keyword arguments:
        
        * mesh_file_name            
            path to the mesh file (.msh format). The partitioned mesh will automatically be loaded in multiprocessing
        * mesh_proj           
73
            mesh projection system (Default: None). 
74
            This argument must not be given if the mesh is writen in 3D (on a sphere), or if all the data are in the same projection system than the mesh.
75
        """
76
77
78
        if dgpy.Msg.GetCommSize() > 1 :
            dgpy.dgMeshPartition(mesh_file_name, dgpy.Msg.GetCommSize())
            mesh_file_name = mesh_file_name[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh'
Philippe Delandmeter's avatar
Philippe Delandmeter committed
79
80
        self._groups = dgpy.dgGroupCollection(mesh_file_name)
        self._mesh_proj = mesh_proj
81
        
Philippe Delandmeter's avatar
Philippe Delandmeter committed
82

Philippe Delandmeter's avatar
Philippe Delandmeter committed
83
class Region :
84
85
86
87
    """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
88
    def __init__(self, mesh, physical_tags='') :
89
90
91
        """keyword arguments:
        
        * mesh
92
            the slimPre Mesh object
93
        * physical_tags
94
95
96
            the list of all the tags on which put some data (optional)
            
            If physical_tags are not set, all the mesh nodes are gathered
97
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
98
        self._mesh = mesh
99
100
        if slim_private._is_string(physical_tags):
            physical_tags = [physical_tags]
Philippe Delandmeter's avatar
Philippe Delandmeter committed
101
102
103
104
105
        self._physical_tags = physical_tags
        self._nodes = []
        self._groups = []
        self._node_map = {}
        self._node_ids = []
106
        self._empty = False
Philippe Delandmeter's avatar
Philippe Delandmeter committed
107
108
        num_ = 0
        self._type = "interface"
109
        for iGroup in range(mesh._groups.getNbElementGroups()+mesh._groups.getNbGhostGroups()) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
110
            group = mesh._groups.getElementGroup(iGroup)
111
            if ((self._physical_tags != ['']) and (not group.getPhysicalTag() in self._physical_tags)) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
                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
127
            self.coordinates = self._evaluateFunctor(dgpy.function.getCoordinates(), 3)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
128
129
130
131
            return

        for iFaceGroup in range(mesh._groups.getNbFaceGroups()) :
            faceGroup = mesh._groups.getFaceGroup(iFaceGroup)
132
            if (not faceGroup.physicalTag() in self._physical_tags) :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
                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)
150
151
        if len(self._nodes) == 0:
            self._empty = True
Philippe Delandmeter's avatar
Philippe Delandmeter committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

    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 :
169
170
171
172
    """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)
173
        Coordinate_system.rotate3D ... 
174
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
175
    def __init__(self, region, data_proj) :
176
177
178
        """keyword arguments:
        
        * region
179
            the slimPre Region object
180
        * data_proj
181
182
            data projection system.

183
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
184
        self._region = region
185
186
        if region._empty:
            return
187
188
189
190
191
192
193
194
195
        print('Projection system for tags:')
        if region._physical_tags == '':
            print('    the entire mesh')
        else:
            if(slim_private._is_string(region._physical_tags)):
                print('    '+region._physical_tags)
            else:
                for tag in region._physical_tags:
                    print('    '+tag)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
196
        mesh_proj = region._mesh._mesh_proj
Philippe Delandmeter's avatar
Philippe Delandmeter committed
197
        if mesh_proj is not None :
198
            print('  Projection used:')
199
200
            print('    mesh proj: %s' % mesh_proj)
            print('    data proj: %s' % data_proj)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
201
202
203
            _mesh_proj = dgpy.pj_init_plus(mesh_proj)
            _data_proj = dgpy.pj_init_plus(data_proj)
        elif mesh_proj is None :
204
205
            print('  Projection used:')
            print('    mesh proj: no projection. Coordinates are 3D')
206
            print('    data proj: %s' % data_proj)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
207
208
209
210
            mesh_proj2 = "+proj=latlong +ellps=WGS84"
            _mesh_proj = dgpy.pj_init_plus(mesh_proj2)
            _data_proj = dgpy.pj_init_plus(data_proj)

211
        xyz = region.coordinates      
Philippe Delandmeter's avatar
Philippe Delandmeter committed
212
        if mesh_proj :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
213
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, xyz)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
214
215
216
217
218
219
220
221
222
223
224
225
            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
226
            self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, lonlat)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
227
228
229
230
            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)
231
            self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
232
233
234
235
236
237
238

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

    def rotate(self, u, v):
239
        """Return the rotated vector data array (u,v). Returned vector will be oriented in the mesh projection system
240
241
242
243
244
245
        
        keyword arguments:
        
        * u
            first component of the vector (in the data coordinate system)
        * v
246
            second component of the vector (in the data coordinate system)
247
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
248
249
250
        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)
251
        
252
    def rotate3D(self, F_lon, F_colat, F_r):
253
254
255
256
257
258
259
260
261
262
263
        """Return ...
        
        keyword arguments:
        
        * F_lon
            this var ...
        * F_colat
            this var ... 
        * F_r
            this var ...
        """
264
265
        lon = self.coordinates[0]
        colat = 0.5*np.pi - self.coordinates[1]
266
267
268
269
        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
270
271
272


class Time :
273
274
    """Create the slimPre Time"""

275
    def __init__(self, time_vector=None, initial_time='1970-01-01 00:00:00', final_time=None, time_step=None, periodic=False) :
276
277
278
        """keyword arguments:

        * time_vector
279
            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
280
        * initial_time
281
            initial_time given in the format 'YYYY-MM-DD hh:mm:ss'. If not set, initial_time is '1970-01-01 00:00:00'
282
        * final_time
283
            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.
284
        * time_step
285
286
287
288
289
            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.
290
        """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
291
292
        self._periodic = periodic
        fmt = '%Y-%m-%d %H:%M:%S'
293
294
295
        if not final_time and not time_step :
            self._time = np.array(time_vector)
            fmt = '%Y-%m-%d %H:%M:%S'
296
297
            date_ref  = slim_private.datetime.datetime.strptime('1970-01-01 00:00:00', fmt)
            date_user = slim_private.datetime.datetime.strptime(initial_time, fmt)
298
299
300
            date_shift = (date_user-date_ref)
            time_shift = date_shift.days * 86400 + date_shift.seconds
            self._time += time_shift
301
        elif not time_vector and initial_time and final_time and time_step:
302
303
304
305
            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])
306
            self._time = np.arange(initial_time, final_time+time_step*0.9999, time_step)
307
308
        else:
            dgpy.Msg.Fatal("Wrong combination of inputs to create a Time object")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
309
310


Philippe Delandmeter's avatar
Philippe Delandmeter committed
311
def write_file(file_name, region=None, time=None, data=[]):
312
    """Create a SLIM netcdf file, standard input for SLIM simulations.
313
314
315
316
317
318

    keyword arguments:

    * file_name
        name of the netcdf file to write
    * region
319
        slimPre Region object (optional).
320
    * time
321
        slimPre Time object (optional).
322
    * data
323
324
        data is a list of data tuples. Each tuple follows the form 
        
325
        ("data_name", data_array)
326
        
327
328
        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.
    """
329
330
    if dgpy.Msg.GetCommSize() > 1 :
        file_name = file_name+"_"+str(dgpy.Msg.GetCommRank())
Philippe Delandmeter's avatar
Philippe Delandmeter committed
331
    oF = dgpy.slimNetCDFWrap(file_name)
332
333
334
    space_dep = region is not None
    time_dep = time is not None
    if space_dep :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
335
        nbNodes = len(region._node_ids)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
336
337
        nodeDim = oF.create_dimension('node', nbNodes)
        nodeVar = oF.create_variable_int('node', [nodeDim])
338
    if time_dep :
Philippe Delandmeter's avatar
Philippe Delandmeter committed
339
340
341
342
343
        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:
344
            timeUnits = 'seconds since the beginning of the period'
345
346
        timeDim = oF.create_dimension('time', nTimes)
        timeVar = oF.create_variable('time', [timeDim], timeUnits, time._periodic)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
347
348
349
    nFields = len(data)
    valVar = [0]*nFields
    for iField in range(nFields) :
350
        if space_dep and time_dep:
351
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim, nodeDim], None, False)
352
353
354
        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:
355
            valVar[iField] = oF.create_variable(data[iField][0], [timeDim], None, False)
356
357
        else:
            valVar[iField] = oF.create_variable(data[iField][0], [], None, False)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
358
    oF.enddef()
359
360
361
362
    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
363
364
365
366
367
368
369
370
371
372
373
374
375
376
    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)


377
def netcdf_to_msh(mesh_file_name, nc_file_name, variable_name, export_file_name, time = 0):
378
    """Export a field in the netcdf format to msh format.
Philippe Delandmeter's avatar
Philippe Delandmeter committed
379
380

    keyword arguments:
381
382
383
384
385
386
387
388
389
390
391
    
    * 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
392
393
    """

394
395
396
397
    if dgpy.Msg.GetCommSize() > 1 :
        dgpy.dgMeshPartition(mesh_file_name, dgpy.Msg.GetCommSize())
        mesh_file_name = mesh_file_name[:-4] + '_' + str(dgpy.Msg.GetCommSize()) + '.msh'

Philippe Delandmeter's avatar
Philippe Delandmeter committed
398
    groups = dgpy.dgGroupCollection(mesh_file_name)
399
400
401
402
    if dgpy.Msg.GetCommSize() == 1 :
        f = dgpy.slimInputRead(nc_file_name, variable_name)
    else :
        f = dgpy.slimInputRead(nc_file_name+'_'+str(dgpy.Msg.GetCommRank()), variable_name)
403
    dgpy.function.getTime().set(time)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
404
405
    groups.exportFunctionMsh(f, export_file_name)

Philippe Delandmeter's avatar
Philippe Delandmeter committed
406
def tpxo_tide(region, time, data_file_name=None): 
407
408
409
410
    """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
411
412
    
    keyword arguments:
413
    
414
    * region
415
        slimPre Region object on which put the data
416
    * time
417
        Time object of the instant at which put the data
418
    * data_file_name
419
        name of the netcdf data file where tides are exported (optional)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
420
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
421
    lonlat = Coordinate_system(region, data_proj='+proj=latlong +ellps=WGS84')
422
    h = u = v = np.empty((1,1))
Philippe Delandmeter's avatar
Philippe Delandmeter committed
423
424
    h_file_name = 'h_tpxo7.2.nc'
    u_file_name = 'u_tpxo7.2.nc'
425
426
    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
427
428
429
    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')
430
431
432
433
    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
434

435
436
437
438
439
440
441
        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
442

443
        (u, v) = lonlat.rotate(u,v)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
444
445

    if data_file_name is not None:
Philippe Delandmeter's avatar
Philippe Delandmeter committed
446
        write_file(data_file_name, region=region, time=time, data=[('h', h), ('u', u), ('v', v)])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
447
448
449
450


    return (h, u, v)

Philippe Delandmeter's avatar
Philippe Delandmeter committed
451
def residual_flow(region, flow, bathymetry, data_file_name=None): 
452
453
454
    """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.
455
456
457
    The physical tags of the region object must be boundary physical tags.
    
    keyword arguments:
458
    
459
    * region
460
        slimPre Region object on which put the data
461
462
463
464
465
    * 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
466
        name of the netcdf data file where the residual flow is exported (optional)
467
    """
Philippe Delandmeter's avatar
Philippe Delandmeter committed
468

Philippe Delandmeter's avatar
Philippe Delandmeter committed
469
    groups = region._mesh._groups
Philippe Delandmeter's avatar
Philippe Delandmeter committed
470
471
472
    integNormals = dgpy.dgFunctionIntegratorInterface(groups, dgpy.function.getNormals())
    nx = 0
    ny = 0
Philippe Delandmeter's avatar
Philippe Delandmeter committed
473
    for tag in region._physical_tags:
474
475
476
477
478
        #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
479
480
481
482
483
484
485
486
487
488
489
    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
490
    sectionFM = dgpy.fullMatrixDouble(1,1,True)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
491
    for tag in region._physical_tags:
492
493
494
        #section += integInterface.compute(tag)
        integInterface.compute(tag, sectionFM)
        section += sectionFM.get(0,0)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
495
496
497
498
499

    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
500
    vals = region._evaluateFunctor(f, 2)
501
502
    if region._empty:
        vals = np.empty((1,2))
503
    if data_file_name:
504
        write_file(data_file_name, region, None, [("uflow", vals[:,0]), ("vflow", vals[:,1])])
Philippe Delandmeter's avatar
Philippe Delandmeter committed
505
506
    return vals

507
def smooth_bathymetry(mesh, bathymetry, output_file_name, coefficient=0.5, transform_p0=True):
508
    """Smooth the bathymetry in order to locally reduce the noise in the momentum
Philippe Delandmeter's avatar
Philippe Delandmeter committed
509
510
    
    keyword arguments:
511
512
513
    
    * mesh 
        slimPre mesh object
514
515
    * bathymetry
        bathymetry file
516
        
517
        Format:
518
519
520
521
            * 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).
522
523
524
    * 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
525
        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
526
527
528
    """
    groups = mesh._groups
    dof = dgpy.dgDofContainer(groups, 1)
Valentin Vallaeys's avatar
Valentin Vallaeys committed
529
    slim_private._load(dof, bathymetry)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
530
531
532
533
534
535
536
537
538
539
    #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) :
540
541
542
                if v >= n:
                    n *= 2
                    data.resize((n))
Philippe Delandmeter's avatar
Philippe Delandmeter committed
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
                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
571
572
        #TODO better print
        print(count,emax)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
573
574
        if emax < .2 :
            break
575
    #read data
Philippe Delandmeter's avatar
Philippe Delandmeter committed
576
577
    for ig in range(groups.getNbElementGroups()) :
        g = groups.getElementGroup(ig)
578
        gdata = dof.getGroupProxy(ig)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
579
580
581
        for ie in range(g.getNbElements()):
            vs = m.vertices(ig, ie)
            for i, v in enumerate(vs) :
582
                gdata.set(i,ie, data[v])
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
    
    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":
598
            dgpy.Msg.Info("The bathymetry is dg and saved in "+output_file_name[:-3]+"/"+slim_private.path.basename(output_file_name[:-3])+".idx")
599
600
601
602
            output_file_name = output_file_name[:-3]
        exporter = dgpy.dgIdxExporter(bath_P0, output_file_name)
        exporter.exportIdx()
        
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
def interpolate_on_structured_grid(x, y, ox, oy, dx, dy, data, fill_vallue=-9999.):
    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
    val = np.where(index, fill_vallue, \
          (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
        
def interpolate_from_structured_grid(nodes, structured_file, fill_value=-9999, transpose_matrix=False):
    """
    """
    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:
            data = np.fromfile(f, "d", nx * ny, sep=' ').reshape((ny, nx))        
        np.save("data.npy",data)
    f.close()
    x = nodes[:,0]
    y = nodes[:,1]
    data = np.transpose(data)
    val = interpolate_on_structured_grid(x, y, ox, oy, dx, dy, data, fill_vallue=fill_value)    
    return val
        
Philippe Delandmeter's avatar
Philippe Delandmeter committed
675
def get_data_from_netcdf(nc_file_name, variable_name):
676
    """Return the array containing the values of the variable in the nectdf file (SLIM format)
Philippe Delandmeter's avatar
Philippe Delandmeter committed
677
678

    keyword arguments:
679
    
680
681
682
683
    * 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
684
685
    """
    #TODO use a wrapper instead of netCDF4
686
    dgpy.Msg.Info("Loading netCDF4 for python, might not be installed on your computer...")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
687
    from netCDF4 import Dataset
688
689
690
691
    if (dgpy.Msg.GetCommSize() > 1):
        f = Dataset(nc_file_name+'_'+str(dgpy.Msg.GetCommRank()),"r")
    else:
        f = Dataset(nc_file_name,"r")
Philippe Delandmeter's avatar
Philippe Delandmeter committed
692
693
694
    val = f.variables[variable_name][:]
    return val

695
def get_distance(region, physical_tags):
696
    """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 
697
698
    
    keyword arguments:
699
    
700
    * region
701
        slimPre Region object of the nodes on which compute the distance
702
703
    * physical_tags
        vector containing the physical tags from which the minimum distance will be computed
704
705
    """
    xyz = region.coordinates
706
    dist = slim_private.distance.distance(region._mesh, physical_tags, evaluationPoints=xyz)
707
708
    return np.array(dist)

709
def fetch_ftp(file_name, data_dir=None):
710
711
    """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. 

712
713
    global_data_directory is by default ~/.cache/slim_data/
    """
714
    file_name =  slim_private.dgftp.get(file_name)
715
716
    if data_dir: 
        if dgpy.Msg.GetCommRank() == 0:
Valentin Vallaeys's avatar
Valentin Vallaeys committed
717
            slim_private.shutil.copyfile(file_name, data_dir+slim_private.path.basename(file_name))
718
            file_name = data_dir+slim_private.path.basename(file_name)
719
720
721
        dgpy.Msg.Barrier()
    return file_name

722
def set_global_data_directory(path):
723
    """Edit the global_data_directory"""
724
725
    slim_private.dgftp.global_data_local_directory = path

726
def make_directory(data_dir):
727
    """Create data_dir directory if it does not exist yet"""
728
    if dgpy.Msg.GetCommRank() == 0:
729
730
        if not slim_private.path.isdir(data_dir) :
            slim_private.makedirs(data_dir)
731
    dgpy.Msg.Barrier()
732
733

def exit(i):
734
    """Exit properly the program"""
735
    dgpy.Msg.Exit(i)