fluid.py 15.9 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 *
Nathan Coppin's avatar
Nathan Coppin committed
32
import numpy as np
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__))
Nathan Coppin's avatar
Nathan Coppin 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
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)
Nathan Coppin's avatar
Nathan Coppin committed
51
52
        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
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) :
Nathan Coppin's avatar
Nathan Coppin committed
86
            tmp = np.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
Nathan Coppin's avatar
Nathan Coppin committed
97
        n_fluids = np.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
Nathan Coppin's avatar
Nathan Coppin committed
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
        """
Nathan Coppin's avatar
Nathan Coppin committed
195
196
197
        v = self.solution()[:,:self._dim]
        if self._dim == 2 :
            v = np.insert(v,self._dim,0,axis=1)
198
199
        data = [
            ("pressure",self.solution()[:,[self._dim]]),
Nathan Coppin's avatar
Nathan Coppin committed
200
            ("velocity",v),
201
202
203
            ("porosity",self.porosity()),
            ("old_porosity",self.old_porosity())
            ]
204
        if self._n_fluids == 2 :
205
206
207
208
209
210
            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))
Nathan Coppin's avatar
Nathan Coppin committed
211
212
213
214
        connectivity = self.elements()
        types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
        offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0]))
        VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data)
215

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
216
    def import_vtk(self, filename) :
217
        """Read output file to reload computation.
Matthieu Constant's avatar
Matthieu Constant committed
218
219
220
221

        Keyword arguments:
        filename -- name of the file to read
        """
222
        def np2c(a,dtype) :
Nathan Coppin's avatar
Nathan Coppin committed
223
            tmp = np.require(a,dtype,"C")
224
225
226
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
Nathan Coppin's avatar
Nathan Coppin committed
227
228
        x,cells,data,cdata,fdata = VTK.read(filename)
        el = cells["connectivity"].reshape([-1,self._dim+1])
229
        self._lib.fluid_problem_private_set_elements(self._fp,
Nathan Coppin's avatar
Nathan Coppin committed
230
231
                c_int(x.shape[0]),np2c(x,np.float64),
                c_int(el.shape[0]),np2c(el,np.int32))
232
233
        phys = {}
        bnds = {}
Nathan Coppin's avatar
Nathan Coppin committed
234
        for n,d in fdata.items() :
235
236
237
238
239
240
241
242
243
            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")),
Nathan Coppin's avatar
Nathan Coppin committed
244
245
                    c_int(nodes.shape[0]), np2c(nodes,np.int32),
                    c_int(bnd.shape[0]), np2c(bnd,np.int32))
246
        sol = self.solution()
Nathan Coppin's avatar
Nathan Coppin committed
247
        sol[:,:self._dim] = data["velocity"][:,:self._dim]
248
        sol[:,[self._dim]] = data["pressure"]
249
        if self._n_fluids == 2 :
250
251
252
253
            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
254

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
255
    def compute_node_force(self, dt) :
Matthieu Constant's avatar
Matthieu Constant committed
256
257
258
259
260
        """Compute the forces to apply on each grain resulting from the fluid motion.

        Keyword arguments:
        dt -- computation time step
        """
Nathan Coppin's avatar
Nathan Coppin committed
261
        forces = np.ndarray([self.n_particles,self._dim],"d",order="C")
262
        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
263
264
        return forces

Matthieu Constant's avatar
Matthieu Constant committed
265
    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
266
267
268
269
        """Solve the fluid equations.

        Keyword arguments:
        dt -- computation time step
270
271
272
        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
273
        """
Matthieu Constant's avatar
Matthieu Constant committed
274
        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
275

Matthieu Constant's avatar
Matthieu Constant committed
276
    def set_particles(self, mass, volume, position, velocity, reload = False):
Matthieu Constant's avatar
Matthieu Constant committed
277
278
279
280
281
282
283
284
285
        """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
        """
286
        def np2c(a) :
Nathan Coppin's avatar
Nathan Coppin committed
287
            tmp = np.require(a,"float","C")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
288
289
290
            r = c_void_p(tmp.ctypes.data)
            r.tmp = tmp
            return r
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
291
        self.n_particles = mass.shape[0]
292
        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
293

294
    def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
295
        f = getattr(self._lib,"fluid_problem_"+f_name)
296
        f.restype = POINTER(typec)
Nathan Coppin's avatar
Nathan Coppin committed
297
        return np.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
298
299

    def solution(self) :
Matthieu Constant's avatar
Matthieu Constant committed
300
        """Give access to the fluid field value at the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
301
302
303
        return self._get_matrix("solution", self.n_nodes(), self.n_fields())

    def coordinates(self) :
Matthieu Constant's avatar
Matthieu Constant committed
304
        """Give access to the coordinate of the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
305
        return self._get_matrix("coordinates",self.n_nodes(), 3)
306
307
308
309
310
311
312
313
314

    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
315
316
    
    def n_fields(self):
Matthieu Constant's avatar
Matthieu Constant committed
317
        """Return the number of fluid fields."""
318
319
        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
320
321

    def n_nodes(self):
Matthieu Constant's avatar
Matthieu Constant committed
322
        """Return the number of mesh nodes."""
323
324
        self._lib.fluid_problem_n_nodes.restype = c_int
        return self._lib.fluid_problem_n_nodes(self._fp);
325
326
327
328
329
330
331
332
    
    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
333

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
334
    def particle_element_id(self) :
Matthieu Constant's avatar
Matthieu Constant committed
335
        """Return the id of the mesh cell in which particles are located."""
336
        f = getattr(self._lib,"fluid_problem_particle_element_id")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
337
        f.restype = c_void_p
338
        fs = getattr(self._lib,"fluid_problem_n_particles")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
339
340
341
342
        fs.restype = c_int
        ptr = f(self._fp)
        size = fs(self._fp)
        buf = (size*c_int).from_address(ptr)
Nathan Coppin's avatar
Nathan Coppin committed
343
        return np.ctypeslib.as_array(buf)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
344

345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
    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))
Nathan Coppin's avatar
Nathan Coppin committed
360
361
        nodes = np.ctypeslib.as_array(phys_nodes,shape=(phys_n_nodes.value,)) if phys_nodes else np.ndarray([0],np.int32) 
        bnd = np.ctypeslib.as_array(phys_boundaries,shape=(phys_n_boundaries.value,2)) if phys_boundaries else np.ndarray([2,0],np.int32)
362
363
        return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd