fluid.py 15.5 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2018>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts 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
7
# Description and complete License: see LICENSE file.
# 	
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
8
# This program (MigFlow) is free software: 
9
# 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
12
13
14
15
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
# 
# 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.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
# 
18
19
# You should have received a copy of the GNU Lesser General Public License
# 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
27
28
29
30
"""Model for Immersed Granular Flow -- Fluid user interface

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

    MigFlow computes immersed granular flows using an unresolved FEM-DEM model. 
    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 *
32
import numpy
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
34
import signal
import os
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
35
import sys
36
from . import VTK
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__))
39
40
lib2 = CDLL(os.path.join(dir_path,"libmbfluid2.so"))
lib3 = CDLL(os.path.join(dir_path,"libmbfluid3.so"))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
41
42
43

signal.signal(signal.SIGINT, signal.SIG_DFL)

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
44
BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
45
46
47
48
49
50
class _Bnd :
    def __init__(self, b, dim) :
        self._b = b
        self._dim = dim
    def apply(self, n, xp, vp) :
        nv = len(self._b)
51
52
        x = numpy.frombuffer(cast(xp, POINTER(int(n)*self._dim*c_double)).contents).reshape([n,-1])
        v = numpy.frombuffer(cast(vp, POINTER(int(n)*nv*c_double)).contents).reshape([n,nv])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
53
54
55
56
57
58
59
60
61
62
63
        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
64

65
class FluidProblem :
66
    """Create numerical representation of the fluid."""
Matthieu Constant's avatar
Matthieu Constant committed
67

68
    def __init__(self, dim, g, mu, rho, coeff_stab=0.01, petsc_solver_type="-pc-type lu"):
Matthieu Constant's avatar
Matthieu Constant committed
69
70
71
72
73
74
75
        """Build the fluid structure.

        Keyword arguments:
        dim -- dimension of the problem
        g -- gravity (scalar)
        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)
Frédéric Dubois's avatar
Frédéric Dubois committed
76
        coeff_stab -- optional argument used as coefficient for extra diffusivity added to stabilise the advection of concentration (only for two fluid problem) 
77
        petsc_solver_type -- optional argument to specify solver option (only when petsc is used)
78
79
80
81

        Raises:
        ValueError -- if the dimension differs from 2 or 3
        NameError -- if C builder cannot be found
Matthieu Constant's avatar
Matthieu Constant committed
82
        """
83
        self.strong_cb_ref = []
Matthieu Constant's avatar
Matthieu Constant committed
84
        self.weak_cb_ref = []
Matthieu Constant's avatar
Matthieu Constant committed
85
        def np2c(a) :
86
            tmp = numpy.require(a,"float","C")
Matthieu Constant's avatar
Matthieu Constant committed
87
88
89
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
90
91
92
93
94
95
96
        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
97
        n_fluids = numpy.require(rho,"float","C").reshape([-1]).shape[0]
98
        self._n_fluids = n_fluids
99
        self._dim = dim
100
        self._fp = c_void_p(self._lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho), c_double(coeff_stab), petsc_solver_type.encode()))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
        if self._fp == None :
102
            raise NameError("Cannot create fluid problem.")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
103
104

    def __del__(self):
Matthieu Constant's avatar
Matthieu Constant committed
105
        """Delete the fluid structure."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
106
        if(self._fp is not None) :
107
            self._lib.fluid_problem_free(self._fp)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
108

109
    def load_msh(self, mesh_file_name) :
Matthieu Constant's avatar
Matthieu Constant committed
110
111
112
113
114
        """Set the domain geometry for the fluid computation.

        Keyword argument:
        mesh_file_name -- name of the mesh.msh file containing information about the domain
        """
115
        self._lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())
116

117
    def set_wall_boundary(self, tag, pressure=None) :
Matthieu Constant's avatar
Matthieu Constant committed
118
119
120
        """Set the weak boundaries (=normal fluxes) for the fluid problem.

        Keyword arguments:
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121
122
        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)
Matthieu Constant's avatar
Matthieu Constant committed
123
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
124
125
126
127
128
        if not _is_string(tag) :
            for t in tag :
                self.set_wall_boundary(t,pressure)
            return
        pressurecb = None
129
        if pressure :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
130
131
132
            pressurecb = BNDCB(_Bnd(pressure,self._dim).apply)
            self.weak_cb_ref.append(pressurecb)
        self._lib.fluid_problem_set_wall_boundary(self._fp,tag.encode(), pressurecb)
133

134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
    def set_open_boundary(self, tag, velocity=None, pressure=None, porosity=1, concentration=1):
        """Set the weak boundaries (=normal fluxes) for the fluid problem.

        Keyword arguments:
        tag -- physical tag (set in the mesh.geo file) of the boundary on which the flux type is added
        bnd_type -- type of the boundary flux: wall, pressure, velocity

        Raises:
        ValueError -- if the number of values set in the list callback_or_value do not match what is expected for the boundary type : wall do not take any values; pressure takes a value for the pressure and concentration if any; velocity takes values for each velocity dimension and the concentration if any.
        """
        if porosity != 1 :
            raise NameError("Only porosity=1 is implemented")
        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")
        if velocity is not None :
            n_value = self._dim+self._n_fluids-1    # u, v, w, a
150
            cb_or_value = [*velocity,concentration] if self._n_fluids == 2 else velocity
151
152
153
154
155
            bnd_type = "velocity"
        else :
            n_value = self._n_fluids                # p, a
            cb_or_value = [pressure,concentration] if self._n_fluids == 2 else [pressure]
            bnd_type = "pressure"
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
156
        bndcb = BNDCB(_Bnd(cb_or_value, self._dim).apply)
Matthieu Constant's avatar
Matthieu Constant committed
157
        self.weak_cb_ref.append(bndcb)
Matthieu Constant's avatar
wip fix    
Matthieu Constant committed
158
        self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bnd_type.encode(),bndcb)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
159

160
    def set_strong_boundary(self, tag, field,callback_or_value, row=None) :
Matthieu Constant's avatar
Matthieu Constant committed
161
162
163
164
165
166
167
168
        """Set the strong boundaries (=constrains) for the fluid problem.

        Keyword arguments:
        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)
        """
169
170
        if row is None :
            row = field
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
171
        bndcb = BNDCB(_Bnd([callback_or_value], self._dim).apply)
172
        self.strong_cb_ref.append(bndcb)
173
        self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
174

175
    def adapt_mesh(self, lcmax, lcmin, n_el, cmax=1, cmin=1) :
Matthieu Constant's avatar
Matthieu Constant committed
176
177
178
179
180
181
182
183
184
        """Adapt the mesh.

        Keyword arguments:
        lcmax -- maximum mesh radius
        lcmin -- minimum mesh radius
        n_el -- number of element wanted
        cmax -- optional argument to weigh maximum gradient used in the adaptation formula
        cmin -- optional argument to weigh minimum gradient used in the adaptation formula
        """
185
        self._lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el), c_double(cmax), c_double(cmin))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
186

187
    def export_vtk(self, output_dir, t, it) :
Matthieu Constant's avatar
Matthieu Constant committed
188
189
190
191
192
193
194
        """Write output file for post-visualization.

        Keyword arguments:
        output_dir -- output directory
        t -- computation time
        it -- computation iteration
        """
195
196
        data = [
            ("pressure",self.solution()[:,[self._dim]]),
197
            ("velocity",self.solution()[:,:self._dim]),
198
199
200
            ("porosity",self.porosity()),
            ("old_porosity",self.old_porosity())
            ]
201
        if self._n_fluids == 2 :
202
203
204
205
206
207
            data.append(("concentration",self.solution()[:,[self._dim+1]]))
        field_data = []
        for i in range(self._n_physicals()) :
            name,dim,pid,nodes,bnd = self._physical(i)
            field_data.append((b"Physical %i %i %s" % (dim,pid,name),nodes[:,None]))
            field_data.append((b"Boundary %i %i %s" % (dim,pid,name),bnd))
208
        VTK.write(output_dir+"/fluid",it,t,self.elements(),self.coordinates(),data,field_data)
209

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
210
    def import_vtk(self, filename) :
211
        """Read output file to reload computation.
Matthieu Constant's avatar
Matthieu Constant committed
212
213
214
215

        Keyword arguments:
        filename -- name of the file to read
        """
216
        def np2c(a,dtype) :
217
            tmp = numpy.require(a,dtype,"C")
218
219
220
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
221
        x,el,data,fdata = VTK.read(filename)
222
        self._lib.fluid_problem_private_set_elements(self._fp,
223
224
                c_int(x.shape[0]),np2c(x,numpy.float64),
                c_int(el.shape[0]),np2c(el,numpy.int32))
225
226
        phys = {}
        bnds = {}
227
        for n,d in fdata :
228
229
230
231
232
233
234
235
236
            k,dim,pid,name = n.split(None,4)
            if k == "Physical" :
                phys[(dim,pid)] = (name,d)
            elif k == "Boundary" :
                bnds[(dim,pid)] = d
        for (dim,pid),(name,nodes) in phys.items() :
            bnd = bnds[(dim,pid)]
            self._lib.fluid_problem_private_add_physical(self._fp,
                    c_int(int(pid)), c_int(int(dim)), c_char_p(name.encode("utf-8")),
237
238
                    c_int(nodes.shape[0]), np2c(nodes,numpy.int32),
                    c_int(bnd.shape[0]), np2c(bnd,numpy.int32))
239
        sol = self.solution()
240
241
        data = dict(data)
        sol[:,:self._dim] = data["velocity"]
242
        sol[:,[self._dim]] = data["pressure"]
243
        if self._n_fluids == 2 :
244
245
246
247
            sol[:,[self._dim+1]] = data["concentration"]
        self.porosity()[:] = data["porosity"]
        self.old_porosity()[:] = data["old_porosity"]
        self._lib.fluid_problem_after_import(self._fp)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
248

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
249
    def compute_node_force(self, dt) :
Matthieu Constant's avatar
Matthieu Constant committed
250
251
252
253
254
        """Compute the forces to apply on each grain resulting from the fluid motion.

        Keyword arguments:
        dt -- computation time step
        """
255
        forces = numpy.ndarray([self.n_particles,self._dim],"d",order="C")
256
        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
257
258
        return forces

Matthieu Constant's avatar
Matthieu Constant committed
259
    def implicit_euler(self, dt, newton_atol=5e-6, newton_rtol=1e-5, newton_max_it=10, reduced_gravity=0, stab_param=0.) :
Matthieu Constant's avatar
Matthieu Constant committed
260
261
262
263
        """Solve the fluid equations.

        Keyword arguments:
        dt -- computation time step
264
265
266
        newton_atol -- optional argument to specify absolute tolerance for nonlinear solver
        newton_rtol -- optional argument to specify relative tolerance for nonlinear solver
        newton_max_it -- optional argument to specify maximum number of iterations for nonlinear solver
Matthieu Constant's avatar
Matthieu Constant committed
267
        """
Matthieu Constant's avatar
Matthieu Constant committed
268
        return self._lib.fluid_problem_implicit_euler(self._fp, c_double(dt), c_double(newton_atol), c_double(newton_rtol), c_int(newton_max_it), c_int(reduced_gravity), c_double(stab_param))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
269

Matthieu Constant's avatar
Matthieu Constant committed
270
    def set_particles(self, mass, volume, position, velocity, reload = False):
Matthieu Constant's avatar
Matthieu Constant committed
271
272
273
274
275
276
277
278
279
        """Set location of the grains in the mesh and compute the porosity in each cell.

        Keyword arguments:
        mass -- list of particles mass
        volume -- list of particles volume
        position -- list of particles centre positions 
        velocity -- list of particles velocity
        reload -- optional arguments specifying if the simulation is reloaded or if it is the first iteration
        """
280
        def np2c(a) :
281
            tmp = numpy.require(a,"float","C")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
282
283
284
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
285
        self.n_particles = mass.shape[0]
286
        self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),None,c_int(reload))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
287

288
    def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
289
        f = getattr(self._lib,"fluid_problem_"+f_name)
290
        f.restype = POINTER(typec)
291
        return numpy.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
292
293

    def solution(self) :
Matthieu Constant's avatar
Matthieu Constant committed
294
        """Give access to the fluid field value at the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
295
296
297
        return self._get_matrix("solution", self.n_nodes(), self.n_fields())

    def coordinates(self) :
Matthieu Constant's avatar
Matthieu Constant committed
298
        """Give access to the coordinate of the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
299
        return self._get_matrix("coordinates",self.n_nodes(), 3)
300
301
302
303
304
305
306
307
308

    def elements(self):
        """Read-only access to the elements of the mesh."""
        return self._get_matrix("elements", self.n_elements(), self._dim+1,c_int)

    def n_elements(self):
        """Return the number of mesh nodes."""
        self._lib.fluid_problem_n_elements.restype = c_int
        return self._lib.fluid_problem_n_elements(self._fp);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
309
310
    
    def n_fields(self):
Matthieu Constant's avatar
Matthieu Constant committed
311
        """Return the number of fluid fields."""
312
313
        self._lib.fluid_problem_n_fields.restype = c_int
        return self._lib.fluid_problem_n_fields(self._fp);
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
314
315

    def n_nodes(self):
Matthieu Constant's avatar
Matthieu Constant committed
316
        """Return the number of mesh nodes."""
317
318
        self._lib.fluid_problem_n_nodes.restype = c_int
        return self._lib.fluid_problem_n_nodes(self._fp);
319
320
321
322
323
324
325
326
    
    def porosity(self):
        """Return the porosity at nodes"""
        return self._get_matrix("porosity", self.n_nodes(), 1)

    def old_porosity(self):
        """Return the porosity at previous time step"""
        return self._get_matrix("old_porosity", self.n_nodes(), 1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
327

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
328
    def particle_element_id(self) :
Matthieu Constant's avatar
Matthieu Constant committed
329
        """Return the id of the mesh cell in which particles are located."""
330
        f = getattr(self._lib,"fluid_problem_particle_element_id")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
331
        f.restype = c_void_p
332
        fs = getattr(self._lib,"fluid_problem_n_particles")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
333
334
335
336
        fs.restype = c_int
        ptr = f(self._fp)
        size = fs(self._fp)
        buf = (size*c_int).from_address(ptr)
337
        return numpy.ctypeslib.as_array(buf)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
338

339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
    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))
354
355
        nodes = numpy.ctypeslib.as_array(phys_nodes,shape=(phys_n_nodes.value,)) if phys_nodes else numpy.ndarray([0],numpy.int32) 
        bnd = numpy.ctypeslib.as_array(phys_boundaries,shape=(phys_n_boundaries.value,2)) if phys_boundaries else numpy.ndarray([2,0],numpy.int32)
356
357
        return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd