fluid.py 29.5 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Michel Henry's avatar
Michel Henry committed
4
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
5
# List of the contributors to the development of MigFlow: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
# Description and complete License: see LICENSE file.
Michel Henry's avatar
Michel Henry committed
7
8
9
#
# This program (MigFlow) is free software:
# you can redistribute it and/or modify it under the terms of the GNU Lesser General
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
Michel Henry's avatar
Michel Henry committed
12
#
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
13
14
15
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16
# GNU Lesser General Public License for more details.
Michel Henry's avatar
Michel Henry committed
17
#
18
# You should have received a copy of the GNU Lesser General Public License
Michel Henry's avatar
Michel Henry committed
19
# along with this program (see COPYING and COPYING.LESSER files).  If not,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

Matthieu Constant's avatar
Matthieu Constant committed
22
23
24
25
26
"""Model for Immersed Granular Flow -- Fluid user interface

    Contact: jonathan.lambrechts@uclouvain.be
    Webpage: www.migflow.be

Michel Henry's avatar
Michel Henry committed
27
    MigFlow computes immersed granular flows using an unresolved FEM-DEM model.
Matthieu Constant's avatar
Matthieu Constant committed
28
29
30
    The continuous phase (mixture) is computed by solving averaged Navier-Stokes equations on unstructured meshes with the continuous finite element method.
"""

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
from ctypes import *
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
32
import numpy as np
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
import os
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
34
import sys
35
from . import VTK
36
from ._tools import gmsh, get_linear_system_package, _np2c
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
37

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
38
dir_path = os.path.dirname(os.path.realpath(__file__))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
39
40
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
41

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42
BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
43
44
45
46
47
48
class _Bnd :
    def __init__(self, b, dim) :
        self._b = b
        self._dim = dim
    def apply(self, n, xp, vp) :
        nv = len(self._b)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
49
50
        x = np.frombuffer(cast(xp, POINTER(int(n)*self._dim*c_double)).contents).reshape([n,-1])
        v = np.frombuffer(cast(vp, POINTER(int(n)*nv*c_double)).contents).reshape([n,nv])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
51
52
53
54
55
56
57
58
59
60
61
        for i in range(nv):
            if callable(self._b[i]) :
                v[:,i] = self._b[i](x)
            else :
                v[:,i] = self._b[i]

def _is_string(s) :
    if sys.version_info >= (3, 0):
        return isinstance(s, str)
    else :
        return isinstance(s, basestring)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
62

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
63
def _load_msh(mesh_file_name, lib, dim) :
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
    """
    mesh_file_name -- Name of the mesh.msh file containing information about the domain
    """
    etype = 2 if dim == 2 else 4
    betype = 1 if dim == 2 else 2
    gmsh.model.add("tmp")
    gmsh.open(mesh_file_name)
    gmsh.model.mesh.renumberNodes()
    ntags, x, _ = gmsh.model.mesh.getNodes()
    x = x.reshape([-1,3])
    x[ntags-1] = x
    _, el = gmsh.model.mesh.getElementsByType(etype)
    el = (el-1).reshape([-1,dim+1])
    bnd = []
    btag = []
    bname = []
    periodic_nodes = []
    for edim, etag  in gmsh.model.getEntities() :
        ptag, cnodes, pnodes, _ = gmsh.model.mesh.getPeriodicNodes(edim, etag)
        if ptag == etag or len(cnodes) == 0: continue
        periodic_nodes.extend(zip(cnodes,pnodes))
    for ip,(pdim, ptag)  in enumerate(gmsh.model.getPhysicalGroups(dim-1)) :
        bname.append(gmsh.model.getPhysicalName(pdim,ptag))
        for etag in gmsh.model.getEntitiesForPhysicalGroup(pdim, ptag):
            for eltype, _, elnodes in zip(*gmsh.model.mesh.getElements(pdim,etag)):
                if eltype != betype: 
                    continue
                bnd.append(elnodes-1)
                btag.append(np.full((elnodes.shape[0]//dim), ip))
    periodic_nodes = np.array(periodic_nodes,dtype=np.int32).reshape([-1,2])-1
    bnd = np.hstack(bnd).reshape([-1,dim])
    btag = np.hstack(btag)
    cbname = (c_char_p*len(bname))(*(i.encode() for i in bname))

    # remove nodes and boundaries not connected to elements
    keepnodes = np.full(x.shape[0], False, bool)
    keepnodes[el] = True
    newidx = np.cumsum(keepnodes)-1
    el = newidx[el]
    x = x[keepnodes,:]
    periodic_nodes = newidx[periodic_nodes]
    keepbnd = np.full(bnd.shape[0], True, bool)
    for i in range(dim) :
        keepbnd = np.logical_and(keepbnd, keepnodes[bnd[:,i]])
    bnd = newidx[bnd][keepbnd,:]
    btag = btag[keepbnd]
    is_parent = np.full([x.shape[0]],True, bool)
    is_parent[periodic_nodes[:,0]] = False
    periodic_parent = np.cumsum(is_parent)-1
    periodic_parent[periodic_nodes[:,0]] = periodic_parent[periodic_nodes[:,1]]
    lib.mesh_new_from_elements.restype = c_void_p
    return c_void_p(lib.mesh_new_from_elements(
            c_int(x.shape[0]),_np2c(x,np.float64),
            c_int(el.shape[0]),_np2c(el,np.int32),
            c_int(bnd.shape[0]),_np2c(bnd,np.int32),
            _np2c(btag,np.int32),c_int(len(cbname)),cbname,
            _np2c(periodic_parent,np.int32)))
Michel Henry's avatar
Michel Henry committed
121

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
122

123
class FluidProblem :
Nathan Coppin's avatar
Nathan Coppin committed
124
    """Creates the numerical representation of the fluid."""
Matthieu Constant's avatar
Matthieu Constant committed
125

Jonathan Lambrechts's avatar
usolid    
Jonathan Lambrechts committed
126
    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., solver=None, solver_options="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True, usolid=False):
Nathan Coppin's avatar
Nathan Coppin committed
127
        """Builds the fluid structure.
Matthieu Constant's avatar
Matthieu Constant committed
128
129

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
130
131
132
133
134
135
        dim -- Dimension of the problem
        g -- Gravity vector
        viscosity -- List of fluid phases viscosities (this should be a vector whom dimension specify if there is one or two fluid)
        density -- List of fluid phases densities (this should be a vector whom dimension specify if there is one or two fluid)
        sigma -- Surface tension (only when there are two fluids)
        coeff_stab -- Optional argument used as coefficient for extra diffusivity added to stabilise the advection of concentration (only for two fluid problem) 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
136
137
        solver -- possible solver are "mumps", "petsc", "scipy"
        solver_options -- Optional argument to specify solver option (only when petsc is used)
Nathan Coppin's avatar
Nathan Coppin committed
138
139
140
        drag_in_stab -- States if the drag force is in the stabilisation term
        drag_coefficient_factor -- Factor multiplicating the drag coefficient
        temporal -- Enables d/dt (i.e. False = steady)
Jonathan Lambrechts's avatar
usolid    
Jonathan Lambrechts committed
141
142
143
        advection -- Enables advective terms (i.e. False = Stokes, True = Navier-Stokes
        usolid -- if False, div.u_solid is evaluated = (porosity-old_porosity)/dt
        )
Nathan Coppin's avatar
Nathan Coppin committed
144

145
146

        Raises:
Nathan Coppin's avatar
Nathan Coppin committed
147
148
        ValueError -- If the dimension differs from 2 or 3
        NameError -- If C builder cannot be found
Matthieu Constant's avatar
Matthieu Constant committed
149
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
150
151
        self.linear_system_package = get_linear_system_package(solver)
        self.solver_options = solver_options
152
        self.strong_cb_ref = []
Matthieu Constant's avatar
Matthieu Constant committed
153
        self.weak_cb_ref = []
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
154
        self.sys = None
155
        self.mean_pressure = None
156
157
158
159
160
161
162
        if dim == 2 :
            self._lib = lib2
        elif dim == 3 :
            self._lib = lib3
        else :
            raise ValueError("dim should be 2 or 3.")
        self._lib.fluid_problem_new.restype = c_void_p
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
163
        n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
164
        self._n_fluids = n_fluids
165
        self._dim = dim
Jonathan Lambrechts's avatar
usolid    
Jonathan Lambrechts committed
166
        self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_double(ip_factor),c_int(temporal), c_int(advection), c_int(usolid)))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
167
        if self._fp == None :
168
            raise NameError("Cannot create fluid problem.")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
169
170

    def __del__(self):
Nathan Coppin's avatar
Nathan Coppin committed
171
        """Deletes the fluid structure."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
172
        if(self._fp is not None) :
173
            self._lib.fluid_problem_free(self._fp)
174
175
176

    def set_coordinates(self, x):
        self._lib.fluid_problem_set_coordinates(self._fp, _np2c(x))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
177

178
    def load_msh(self, mesh_file_name) :
Nathan Coppin's avatar
Nathan Coppin committed
179
        """Sets the domain geometry for the fluid computation.
Matthieu Constant's avatar
Matthieu Constant committed
180
181

        Keyword argument:
Nathan Coppin's avatar
Nathan Coppin committed
182
        mesh_file_name -- Name of the mesh.msh file containing information about the domain
Matthieu Constant's avatar
Matthieu Constant committed
183
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
184
185
        mesh = _load_msh(mesh_file_name, self._lib, self.dimension())
        self._lib.fluid_problem_set_mesh(self._fp, mesh)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
186
        self.sys = None
187
        gmsh.model.remove()
188

189
190
191
192
    def set_mean_pressure(self, mean_pressure):
        self.mean_pressure = mean_pressure
        self.sys = None

Michel Henry's avatar
Michel Henry committed
193
    def set_wall_boundary(self, tag, pressure=None, velocity=None, compute_viscous_term=1) :
Nathan Coppin's avatar
Nathan Coppin committed
194
        """Sets the weak boundaries (=normal fluxes) for the fluid problem.
Matthieu Constant's avatar
Matthieu Constant committed
195
196

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
197
198
199
        tag -- Physical tag (or list of tags), set in the mesh.geo file, of the wall boundaries
        pressure -- The pressure value if imposed (callback or number)
        velocity -- The velocity value if imposed (callback or number)
Matthieu Constant's avatar
Matthieu Constant committed
200
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
201
202
203
204
        if not _is_string(tag) :
            for t in tag :
                self.set_wall_boundary(t,pressure)
            return
205
        pid = -1
206
207
        vid = -1
        cb_or_value = None
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
208
        if pressure is not None:
209
            cb_or_value = [pressure]
210
            pid = 0
211
212
213
214
215
216
217
218
219
        if velocity is not None:
            cb_or_value = velocity
            vid = 0
        if pressure is not None and velocity is not None:
            cb_or_value = velocity+[pressure]
            vid = 0
            pid = self._dim
        bndcb = BNDCB(_Bnd(cb_or_value,self._dim).apply)
        self.weak_cb_ref.append(bndcb)
Michel Henry's avatar
Michel Henry committed
220
        self._lib.fluid_problem_set_weak_boundary(self._fp,tag.encode(), c_int(0), bndcb, c_int(vid), c_int(pid), c_int(-1), c_int(-1), c_int(int(compute_viscous_term)))
221

Michel Henry's avatar
Michel Henry committed
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
    def set_symmetry_boundary(self, tag, pressure=None):
        """Sets the symmetry boundaries (=normal fluxes) for the fluid problem.

        Keyword arguments:
        tag -- Physical tag (or list of tags), set in the mesh.geo file, of the symmetry boundaries
        pressure -- The pressure value if imposed (callback or number)
        velocity -- The velocity value if imposed (callback or number)
        """
        if not _is_string(tag):
            for t in tag:
                self.set_symmetry_boundary(t,pressure)
            return
        pid = -1
        vid = -1
        cb_or_value = None
        if pressure is not None:
            cb_or_value = [pressure]
            pid = 0
        bndcb = BNDCB(_Bnd(cb_or_value, self._dim).apply)
        self.weak_cb_ref.append(bndcb)
242
        self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(), c_int(0), bndcb, c_int(vid), c_int(pid), c_int(-1), c_int(-1), c_int(0))
Michel Henry's avatar
Michel Henry committed
243
244

    def set_open_boundary(self, tag, velocity=None, pressure=None, porosity=1, concentration=1, compute_viscous_term=1):
Nathan Coppin's avatar
Nathan Coppin committed
245
        """Sets the weak boundaries (=normal fluxes) for the fluid problem.
246
247

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
248
249
250
251
252
        tag -- Physical tag (set in the mesh.geo file) of the boundary on which the flux type is added
        velocity -- The velocity value if imposed (callback or number)
        pressure -- The pressure value if imposed (callback or number)
        porosity -- Porosity value outside the boudary
        concentration -- Concentration outside the boundary
253
254

        Raises:
Matthieu Constant's avatar
Matthieu Constant committed
255
256
        NameError -- If the specified porosity outside the boundary is too small
        NameError -- If velocity and pressure are not specified or if they are both specified at the open boundary. It should be one or the other
Nathan Coppin's avatar
Nathan Coppin committed
257
        NameError -- If the dimension of the velocity vector is not equal to the physical dimension of the problem
258
        """
Matthieu Constant's avatar
Matthieu Constant committed
259
260
        if porosity < 1e-3 :
            raise NameError("Inflow porosity too small!")
261
262
        if (velocity is None and pressure is None) or (velocity is not None and pressure is not None)  :
            raise NameError("Pressure or Velocity (but not both) should be specified at open boundaries")
263
264
        if (velocity is not None and len(velocity) != self._dim)  :
            raise NameError("len(velocity) != dimension at open boundaries")
265
        if velocity is not None :
Matthieu Constant's avatar
Matthieu Constant committed
266
            cb_or_value = velocity+[porosity]
267
            vid = 0
Matthieu Constant's avatar
Matthieu Constant committed
268
            cid = self._dim
269
            pid = -1
Matthieu Constant's avatar
Matthieu Constant committed
270
            n = self._dim+1
271
        else :
Matthieu Constant's avatar
Matthieu Constant committed
272
            cb_or_value = [pressure]+[porosity]
273
            pid = 0
Matthieu Constant's avatar
Matthieu Constant committed
274
            cid = 1
275
            vid = -1
Matthieu Constant's avatar
Matthieu Constant committed
276
            n = 2
277
278
279
280
281
282
        if self._n_fluids == 2 :
            cb_or_value += [concentration]
            aid = n
        else :
            aid = -1

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
283
        bndcb = BNDCB(_Bnd(cb_or_value, self._dim).apply)
Matthieu Constant's avatar
Matthieu Constant committed
284
        self.weak_cb_ref.append(bndcb)
Michel Henry's avatar
Michel Henry committed
285
        self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),c_int(1),bndcb,c_int(vid),c_int(pid),c_int(cid),c_int(aid),c_int(int(compute_viscous_term)))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
286

287
    def set_strong_boundary(self, tag, field,callback_or_value, row=None) :
Nathan Coppin's avatar
Nathan Coppin committed
288
        """Sets the strong boundaries (=constrains) for the fluid problem.
Matthieu Constant's avatar
Matthieu Constant committed
289
290

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
291
292
293
294
        tag -- Physical tag (set in the mesh.geo file) of the boundary on which the constraint is added
        field -- Constrained field O: x-velocity; 1: y-velocity; 2: pressure
        callback_or_value -- Value assigned to the specified field on the specified boundary
        row -- Optional argument specifying the row on which the constrain has to be applied (if None the constrain is normally applied on the row field)
Matthieu Constant's avatar
Matthieu Constant committed
295
        """
296
297
        if row is None :
            row = field
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
298
        bndcb = BNDCB(_Bnd([callback_or_value], self._dim).apply)
299
        self.strong_cb_ref.append(bndcb)
300
        self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
301

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
302
    def adapt_mesh(self, lcmax, lcmin, n_el, old_n_particles, old_particle_position, old_particle_volume, cmax=1, cmin=1) :
Nathan Coppin's avatar
Nathan Coppin committed
303
        """Adapts the mesh.
Matthieu Constant's avatar
Matthieu Constant committed
304
305

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
306
307
308
309
310
311
312
313
        lcmax -- Maximum mesh radius
        lcmin -- Minimum mesh radius
        n_el -- Number of element wanted
        old_n_particles -- Number of particles at the previous time step
        old_particle_position -- Position of the particles at the previous time step
        old_particle_volume -- Volume of the particles at the previous time step
        cmax -- Optional argument to weigh maximum gradient used in the adaptation formula
        cmin -- Optional argument to weigh minimum gradient used in the adaptation formula
Matthieu Constant's avatar
Matthieu Constant committed
314
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
315
316
317
        self._lib.fluid_problem_adapt_gen_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin))
        mesh = _load_msh("adapt/mesh.msh", self._lib, self.dimension())
        self._lib.fluid_problem_adapt_mesh(self._fp, mesh, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin), c_int(old_n_particles), _np2c(old_particle_position), _np2c(old_particle_volume))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
318
        self.sys = None
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
319

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
320
    def mesh_boundaries(self) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
321
322
323
324
325
326
327
328
329
330
331
        n = self._lib.fluid_problem_mesh_n_boundaries(self._fp)
        bnds = {}
        for i in range(n) :
            bsize = c_int()
            bname = c_char_p()
            self._lib.fluid_problem_mesh_boundary_info(self._fp, c_int(i),byref(bname),byref(bsize))
            nodes = np.ndarray([bsize.value,self._dim],np.int32)
            self._lib.fluid_problem_mesh_boundary_interface_nodes(self._fp,c_int(i),c_void_p(nodes.ctypes.data))
            bnds[bname.value] = nodes
        return bnds

332
    def write_vtk(self, output_dir, it ,t, stab=False):
Nathan Coppin's avatar
Nathan Coppin committed
333
        """Writes output file for post-visualization.
Matthieu Constant's avatar
Matthieu Constant committed
334
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
335
336
        output_dir -- Output directory
        it -- Computation iteration
337
        t -- Computation time
Nathan Coppin's avatar
Nathan Coppin committed
338
        stab -- If True exports the stabilisation parametres in the output files
Matthieu Constant's avatar
Matthieu Constant committed
339
        """
340
        v = self.solution()[:,:self._dim]
Matthieu Constant's avatar
Matthieu Constant committed
341
        da = self.concentration_dg_grad()
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
342
        cell_data = []
343
        if self._dim == 2 :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
344
            v = np.insert(v,self._dim,0,axis=1)
Matthieu Constant's avatar
Matthieu Constant committed
345
            da = np.insert(da,self._dim,0,axis=1)
Jonathan Lambrechts's avatar
usolid    
Jonathan Lambrechts committed
346
347
        us = np.empty((self.n_nodes(),self.dimension()))
        self._lib.fluid_problem_u_solid(self._fp, _np2c(us))
348
349
        data = [
            ("pressure",self.solution()[:,[self._dim]]),
350
            ("velocity",v),
351
            ("porosity",self.porosity()),
Nathan Coppin's avatar
Nathan Coppin committed
352
            ("old_porosity",self.old_porosity()),
Jonathan Lambrechts's avatar
usolid    
Jonathan Lambrechts committed
353
            ("us", us),
Michel Henry's avatar
VTK    
Michel Henry committed
354
355
            ("grad",da),
            ("parent_node_id", self.parent_nodes())
356
            ]
357
        if self._n_fluids == 2 :
Michel Henry's avatar
VTK    
Michel Henry committed
358
            cell_data.append(("curvature",self.curvature()))
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
359
            cell_data.append(("concentration",self.concentration_dg()))
Matthieu Constant's avatar
Matthieu Constant committed
360
            cell_data.append(("stab",self._get_matrix("stab_param",self.n_elements(),1)))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
361
        field_data = [(b"Boundary %s"%(name), nodes) for name,nodes in self.mesh_boundaries().items()]
362
        connectivity = self.elements()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
363
        types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
364
        offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
365
        VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
366

367
    def read_vtk(self, dirname, i):
Nathan Coppin's avatar
Nathan Coppin committed
368
        """Reads output file to reload computation.
Matthieu Constant's avatar
Matthieu Constant committed
369
370

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
371
        filename -- Name of the file to read
Matthieu Constant's avatar
Matthieu Constant committed
372
        """
373
        filename = dirname + "/fluid_%05d.vtu"%i
374
        x,cells,data,cdata,fdata = VTK.read(filename)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
375
376
        mesh_boundaries = {name[9:]:nodes for name,nodes in fdata.items() if name.startswith("Boundary ")}
        cbnd_names = (c_char_p*len(mesh_boundaries))(*(i.encode() for i in mesh_boundaries.keys()))
377
        el = cells["connectivity"].reshape([-1,self._dim+1])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
378
379
380
381
        nbnd = len(mesh_boundaries)
        bnds = np.vstack(list(mesh_boundaries.values()))
        bnd_tags = np.repeat(list(range(nbnd)),list([v.shape[0] for v in mesh_boundaries.values()]))
        bnd_tags = np.require(bnd_tags,np.int32,"C")
Michel Henry's avatar
VTK    
Michel Henry committed
382
383
        self._lib.mesh_new_from_elements.restype = c_void_p
        _mesh = c_void_p(self._lib.mesh_new_from_elements(
384
385
                c_int(x.shape[0]),_np2c(x,np.float64),
                c_int(el.shape[0]),_np2c(el,np.int32),
Michel Henry's avatar
VTK    
Michel Henry committed
386
387
                c_int(bnds.shape[0]),_np2c(bnds,np.int32),
                _np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
388
                _np2c(data["parent_node_id"],np.int32) if "parent_node_id"  in data else None))
Michel Henry's avatar
VTK    
Michel Henry committed
389
        self._lib.fluid_problem_set_mesh(self._fp, _mesh)
390
        sol = self.solution()
391
        sol[:,:self._dim] = data["velocity"][:,:self._dim]
392
        sol[:,[self._dim]] = data["pressure"]
393
        if self._n_fluids == 2 :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
394
            self.concentration_dg()[:] = cdata["concentration"]
395
        self.porosity()[:] = data["porosity"]
Nathan Coppin's avatar
Nathan Coppin committed
396
        self.old_porosity()[:] = data["old_porosity"]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
397
        self.sys = None
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
398

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
399
    def compute_node_force(self, dt) :
Nathan Coppin's avatar
Nathan Coppin committed
400
        """Computes the forces to apply on each grain resulting from the fluid motion.
Matthieu Constant's avatar
Matthieu Constant committed
401
402

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
403
        dt -- Computation time step
Matthieu Constant's avatar
Matthieu Constant committed
404
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
405
        forces = np.ndarray([self.n_particles,self._dim],"d",order="C")
406
        self._lib.fluid_problem_compute_node_particle_force(self._fp, c_double(dt), c_void_p(forces.ctypes.data))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
407
408
        return forces

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
409
    def compute_node_torque(self, dt) :
Nathan Coppin's avatar
Nathan Coppin committed
410
411
        """Computes the angular drags to apply on each grain resulting from the fluid motion.
        Only in 2D
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
412
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
413
        dt -- Computation time step
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
414
415
        """
        torques = np.ndarray([self.n_particles,1],"d",order="C")
Michel Henry's avatar
Michel Henry committed
416

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
417
418
419
420
        self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data))
        return torques

    def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0.) :
Nathan Coppin's avatar
Nathan Coppin committed
421
        """Solves the fluid equations.
Matthieu Constant's avatar
Matthieu Constant committed
422
423

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
424
425
426
427
        dt -- Computation time step
        check_residual_norm -- If > 0, checks the residual norm after the system resolution
        reduced_graviy -- If True simulations are run with a reduced gravity (not to use in two fluids simulations)
        stab_param -- If zero pspg/supg/lsic stabilisation is computed otherwise the value is used as a coefficient for a pressure laplacian stabilisation term
Matthieu Constant's avatar
Matthieu Constant committed
428
        """
429
430
        self._lib.fluid_problem_set_reduced_gravity(self._fp,c_int(reduced_gravity))
        self._lib.fluid_problem_set_stab_param(self._fp,c_double(stab_param))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
431
432
433
434
        self._lib.fluid_problem_apply_boundary_conditions(self._fp)
        solution_old = self.solution().copy()
        self._lib.fluid_problem_reset_boundary_force(self._fp)
        self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
Michel Henry's avatar
Michel Henry committed
435

Michel Henry's avatar
Michel Henry committed
436
        periodic_map = self._get_matrix("periodic_mapping", self.n_nodes(), 1, c_int).flatten()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
437
        if self.sys is None :
438
439
440
441
442
443
444
            constraints = []
            if self.mean_pressure is not None:
                node_volume = self.node_volume()
                # todo periodic count !!!
                constraint = np.column_stack((np.arange(self.n_nodes()), np.repeat(self._dim, self.n_nodes())))
                constraints.append(constraint)
            self.sys = self.linear_system_package(periodic_map[self.elements()],self.n_fields(),self.solver_options, constraints)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
445
        sys = self.sys
446
447
        localv = np.ndarray([sys.localsize*self.n_elements()])
        localm = np.ndarray([sys.localsize**2*self.n_elements()])
448
449
450
451
452
453
        self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
        constraint_value = []
        if self.mean_pressure is not None:
            constraint_value.append((self.node_volume().reshape([-1]), -self.mean_pressure*np.sum(self.node_volume())))
        rhs = sys.local_to_global(localv,localm,self.solution().reshape([-1]),constraint_value)
        self.solution()[:] -= sys.solve(rhs).reshape([-1, self.n_fields()])[periodic_map]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
454
        if check_residual_norm > 0 :
455
456
            self._lib.fluid_problem_assemble_system(self._fp,_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
            rhs = sys.local_to_global(localv,localm,self.solution().reshape([-1]),constraint_value)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
457
458
            norm = np.linalg.norm(rhs)
            print("norm",norm)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
459
460
            if norm > check_residual_norm :
                raise ValueError("wrong derivative or linear system precision")
461
        self._lib.fluid_problem_node_force_volume(self._fp,_np2c(solution_old),c_double(dt),None,None)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
462

Matthieu Constant's avatar
Matthieu Constant committed
463
    def compute_cfl(self):
Nathan Coppin's avatar
Nathan Coppin committed
464
        """Computes the CFL number divided by the time step
Matthieu Constant's avatar
Matthieu Constant committed
465
466
467
468
469
470
471
472

        """
        nv = np.linalg.norm(self.solution()[:,:self._dim],axis=1)
        nvmax = np.max(nv[self.elements()],axis=1,keepdims=True)
        h = self._get_matrix("element_size",self.n_elements(),1)
        cfl = nvmax / (0.1*h)
        return np.max(cfl)

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
473
    def advance_concentration(self,dt):
Nathan Coppin's avatar
Nathan Coppin committed
474
        """Solves the advection equation for the concentration using the current velocity field.
Matthieu Constant's avatar
Matthieu Constant committed
475
476

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
477
        dt -- Computation time step
Matthieu Constant's avatar
Matthieu Constant committed
478
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
479
        if self._n_fluids == 2 :
Matthieu Constant's avatar
Matthieu Constant committed
480
481
            cfl = self.compute_cfl()
            nsub = int(cfl*dt+1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
482
483
484
485
            if (nsub != 1) :
                print("sub-iterating advection for cfl : %i sub-iterations"%nsub)
            for i in range(nsub) :
                self._lib.fluid_problem_advance_concentration(self._fp,c_double(dt/nsub))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
486

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
487
    def set_particles(self, mass, volume, position, velocity, contact):
Matthieu Constant's avatar
Matthieu Constant committed
488
489
490
        """Set location of the grains in the mesh and compute the porosity in each cell.

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
491
492
493
494
        mass -- List of particles mass
        volume -- List of particles volume
        position -- List of particles centre positions 
        velocity -- List of particles velocity
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
495
        contact -- List of particles contact resultant forces
Matthieu Constant's avatar
Matthieu Constant committed
496
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
497
        self.n_particles = mass.shape[0]
Jonathan Lambrechts's avatar
usolid    
Jonathan Lambrechts committed
498
        self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(position), _np2c(velocity),_np2c(contact))
Nathan Coppin's avatar
Nathan Coppin committed
499

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
500
501
    def move_particles(self, position, velocity, contact):
        """Set location of the grains in the mesh and compute the porosity in each cell.
502

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
503
504
505
506
507
508
        Keyword arguments:
        position -- List of particles centre positions 
        velocity -- List of particles velocity
        contact -- List of particles contact resultant forces
        """
        self._lib.fluid_problem_move_particles(self._fp,c_int(velocity.shape[0]),_np2c(position),_np2c(velocity),_np2c(contact))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
509

510
    def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
511
        f = getattr(self._lib,"fluid_problem_"+f_name)
512
        f.restype = POINTER(typec)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
513
        return np.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
514
515

    def solution(self) :
Nathan Coppin's avatar
Nathan Coppin committed
516
        """Gives access to the fluid field value at the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
517
518
        return self._get_matrix("solution", self.n_nodes(), self.n_fields())

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
519
    def velocity(self) :
Nathan Coppin's avatar
Nathan Coppin committed
520
        """Gives access to the fluid velocity value at the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
521
        return self.solution()[:,:-1]
Michel Henry's avatar
ALE    
Michel Henry committed
522
523
524
525
        
    def mesh_velocity(self) :
        """Give access to the mesh velocity value at the mesh nodes."""
        return self._get_matrix("mesh_velocity", self.n_nodes(), self._dim)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
526
527

    def pressure(self) :
Nathan Coppin's avatar
Nathan Coppin committed
528
        """Gives access to the fluid field value at the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
529
530
        return self.solution()[:,[-1]]

531
532
533
534
    def boundary_force(self) :
        """Give access to force exerted by the fluid on the boundaries."""
        return self._get_matrix("boundary_force", self._lib.fluid_problem_mesh_n_boundaries(self._fp), self._dim)

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
535
    def coordinates(self) :
Nathan Coppin's avatar
Nathan Coppin committed
536
        """Gives access to the coordinate of the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
537
        return self._get_matrix("coordinates",self.n_nodes(), 3)
538

Michel Henry's avatar
VTK    
Michel Henry committed
539
540
541
542
    def parent_nodes(self):
        """Gives access to the parent nodes of each node."""
        return self._get_matrix("periodic_mapping", self.n_nodes(),1, typec=c_int32)

543
    def n_fluids(self) :
Nathan Coppin's avatar
Nathan Coppin committed
544
        """Returns the number of fluids."""
545
546
        return self._n_fluids

547
    def elements(self):
Nathan Coppin's avatar
Nathan Coppin committed
548
        """Gives read-only access to the elements of the mesh."""
549
550
551
        return self._get_matrix("elements", self.n_elements(), self._dim+1,c_int)

    def n_elements(self):
Nathan Coppin's avatar
Nathan Coppin committed
552
        """Returns the number of mesh nodes."""
553
        self._lib.fluid_problem_n_elements.restype = c_int
Michel Henry's avatar
VTK    
Michel Henry committed
554
        return self._lib.fluid_problem_n_elements(self._fp)
Michel Henry's avatar
Michel Henry committed
555

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
556
    def n_fields(self):
Nathan Coppin's avatar
Nathan Coppin committed
557
        """Returns the number of fluid fields."""
558
        self._lib.fluid_problem_n_fields.restype = c_int
Michel Henry's avatar
VTK    
Michel Henry committed
559
        return self._lib.fluid_problem_n_fields(self._fp)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
560
561

    def n_nodes(self):
Nathan Coppin's avatar
Nathan Coppin committed
562
        """Returns the number of mesh nodes."""
563
        self._lib.fluid_problem_n_nodes.restype = c_int
Michel Henry's avatar
VTK    
Michel Henry committed
564
        return self._lib.fluid_problem_n_nodes(self._fp)
Michel Henry's avatar
Michel Henry committed
565

566
    def porosity(self):
Nathan Coppin's avatar
Nathan Coppin committed
567
        """Returns the porosity at nodes"""
568
569
        return self._get_matrix("porosity", self.n_nodes(), 1)

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
570
    def set_concentration_cg(self,concentration):
Nathan Coppin's avatar
Nathan Coppin committed
571
        """Sets the concentration at nodes"""
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
572
        concentration = concentration.reshape((self.n_nodes(),1))
573
        self._lib.fluid_problem_set_concentration_cg(self._fp,_np2c(concentration))
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
574
575

    def concentration_dg(self):
Nathan Coppin's avatar
Nathan Coppin committed
576
        """Returns the concentration at discontinuous nodes"""
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
577
578
        return self._get_matrix("concentration_dg", self.n_elements(), self._dim+1)

Matthieu Constant's avatar
Matthieu Constant committed
579
    def concentration_dg_grad(self):
Nathan Coppin's avatar
Nathan Coppin committed
580
        """Returns the concentration at discontinuous nodes"""
Matthieu Constant's avatar
Matthieu Constant committed
581
582
        return self._get_matrix("concentration_dg_grad", self.n_nodes(), self._dim)

583
    def curvature(self):
Nathan Coppin's avatar
Nathan Coppin committed
584
        """Returns the porosity at previous time step"""
585
586
        return self._get_matrix("curvature", self.n_elements(), 1)

Nathan Coppin's avatar
Nathan Coppin committed
587
    def old_porosity(self):
Nathan Coppin's avatar
Nathan Coppin committed
588
        """Returns the porosity at previous time step"""
Nathan Coppin's avatar
Nathan Coppin committed
589
        return self._get_matrix("old_porosity", self.n_nodes(), 1)
Michel Henry's avatar
Michel Henry committed
590

591
    def volume_flux(self, bnd_tag):
Nathan Coppin's avatar
Nathan Coppin committed
592
        """Computes the integral of the (outgoing) normal velocity through boundary with tag bnd_tag"""
593
594
595
        self._lib.fluid_problem_volume_flux.restype = c_double
        return self._lib.fluid_problem_volume_flux(self._fp,bnd_tag.encode())

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
596
    def particle_element_id(self) :
Nathan Coppin's avatar
Nathan Coppin committed
597
        """Returns the id of the mesh cell in which particles are located."""
598
        f = getattr(self._lib,"fluid_problem_particle_element_id")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
599
        f.restype = c_void_p
600
        fs = getattr(self._lib,"fluid_problem_n_particles")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
601
602
603
604
        fs.restype = c_int
        ptr = f(self._fp)
        size = fs(self._fp)
        buf = (size*c_int).from_address(ptr)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
605
        return np.ctypeslib.as_array(buf)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
606

607
    def particle_uvw(self) :
Nathan Coppin's avatar
Nathan Coppin committed
608
        """Returns the coordinates of the particles inside their element"""
609
610
        return self._get_matrix("particle_uvw",self.n_particles,self._dim)

Matthieu Constant's avatar
Matthieu Constant committed
611
    def particle_position(self) :
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
612
        """Gives access to stored paricles position."""
Matthieu Constant's avatar
Matthieu Constant committed
613
614
        return self._get_matrix("particle_position", self.n_particles, self._dim)

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
615
616
617
618
    def particle_velocity(self) :
        """Gives access to the stored particle velocity."""
        return self._get_matrix("particle_velocity", self.n_particles, self._dim)

Matthieu Constant's avatar
Matthieu Constant committed
619
    def particle_volume(self) :
Nathan Coppin's avatar
Nathan Coppin committed
620
        """Gives access to the fluid field value at the mesh nodes."""
Matthieu Constant's avatar
Matthieu Constant committed
621
622
        return self._get_matrix("particle_volume", self.n_particles, 1)

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
623
    def bulk_force(self) :
Nathan Coppin's avatar
Nathan Coppin committed
624
        """Gives access to the volume force at fluid nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
625
626
        return self._get_matrix("bulk_force", self.n_nodes(), self._dim)

627
    def node_volume(self) :
Nathan Coppin's avatar
Nathan Coppin committed
628
        """Returns the volume associated with each node"""
629
630
        return self._get_matrix("node_volume",self.n_nodes(),1)

631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
    def _n_physicals(self):
        f = self._lib.fluid_problem_private_n_physical
        f.restype = c_int
        return f(self._fp)

    def _physical(self,ibnd):
        f = self._lib.fluid_problem_private_physical
        phys_dim = c_int()
        phys_id = c_int()
        phys_n_nodes = c_int()
        phys_nodes = POINTER(c_int)()
        phys_name = c_char_p()
        phys_n_boundaries = c_int()
        phys_boundaries = POINTER(c_int)()
        f(self._fp,c_int(ibnd),byref(phys_name),byref(phys_dim),byref(phys_id),byref(phys_n_nodes),byref(phys_nodes),byref(phys_n_boundaries),byref(phys_boundaries))
Michel Henry's avatar
Michel Henry committed
646
        nodes = np.ctypeslib.as_array(phys_nodes,shape=(phys_n_nodes.value,)) if phys_nodes else np.ndarray([0],np.int32)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
647
        bnd = np.ctypeslib.as_array(phys_boundaries,shape=(phys_n_boundaries.value,2)) if phys_boundaries else np.ndarray([2,0],np.int32)
648
        return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
649
650
651

    def dimension(self) :
        return self._dim
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
652