fluid.py 23.8 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 *
Jonathan Lambrechts's avatar
Jonathan Lambrechts 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
38
39
40
41
try :
    from .petsclsys import LinearSystem
except :
    print("PETSc4py not available, falling back to scipy linear system")
    from .scipylsys import LinearSystem
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
43
dir_path = os.path.dirname(os.path.realpath(__file__))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
44
45
lib2 = np.ctypeslib.load_library("libmbfluid2",dir_path)
lib3 = np.ctypeslib.load_library("libmbfluid3",dir_path)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
46
47
48

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
49
BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
50
51
52
53
54
55
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
56
57
        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
58
59
60
61
62
63
64
65
66
67
68
        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
69

70
71
72
73
74
75
def _np2c(a,dtype=float) :
    tmp = np.require(a,dtype,"C")
    r = c_void_p(tmp.ctypes.data)
    r.tmp = tmp
    return r

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
76
    
77
class FluidProblem :
Nathan Coppin's avatar
Nathan Coppin committed
78
    """Creates the numerical representation of the fluid."""
Matthieu Constant's avatar
Matthieu Constant committed
79

Matthieu Constant's avatar
Matthieu Constant committed
80
    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True):
Nathan Coppin's avatar
Nathan Coppin committed
81
        """Builds the fluid structure.
Matthieu Constant's avatar
Matthieu Constant committed
82
83

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
84
85
86
87
88
89
90
91
92
93
94
        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) 
        petsc_solver_type -- Optional argument to specify solver option (only when petsc is used)
        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)
        advection -- Enables advective terms (i.e. False = Stokes, True = Navier-Stokes)
Nathan Coppin's avatar
Nathan Coppin committed
95

96
97

        Raises:
Nathan Coppin's avatar
Nathan Coppin committed
98
99
        ValueError -- If the dimension differs from 2 or 3
        NameError -- If C builder cannot be found
Matthieu Constant's avatar
Matthieu Constant committed
100
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
101
        self.solver_options = petsc_solver_type
102
        self.strong_cb_ref = []
Matthieu Constant's avatar
Matthieu Constant committed
103
        self.weak_cb_ref = []
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
104
        self.sys = None
105
106
107
108
109
110
111
        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
112
        n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
113
        self._n_fluids = n_fluids
114
        self._dim = dim
Matthieu Constant's avatar
Matthieu Constant committed
115
        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)))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
116
        if self._fp == None :
117
            raise NameError("Cannot create fluid problem.")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
118
119

    def __del__(self):
Nathan Coppin's avatar
Nathan Coppin committed
120
        """Deletes the fluid structure."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
121
        if(self._fp is not None) :
122
            self._lib.fluid_problem_free(self._fp)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
123

124
    def load_msh(self, mesh_file_name) :
Nathan Coppin's avatar
Nathan Coppin committed
125
        """Sets the domain geometry for the fluid computation.
Matthieu Constant's avatar
Matthieu Constant committed
126
127

        Keyword argument:
Nathan Coppin's avatar
Nathan Coppin committed
128
        mesh_file_name -- Name of the mesh.msh file containing information about the domain
Matthieu Constant's avatar
Matthieu Constant committed
129
        """
130
        self._lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
131
        self.sys = None
132

133
    def set_wall_boundary(self, tag, pressure=None, velocity=None) :
Nathan Coppin's avatar
Nathan Coppin committed
134
        """Sets the weak boundaries (=normal fluxes) for the fluid problem.
Matthieu Constant's avatar
Matthieu Constant committed
135
136

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
137
138
139
        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
140
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
141
142
143
144
        if not _is_string(tag) :
            for t in tag :
                self.set_wall_boundary(t,pressure)
            return
145
        pid = -1
146
147
        vid = -1
        cb_or_value = None
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
148
        if pressure is not None:
149
            cb_or_value = [pressure]
150
            pid = 0
151
152
153
154
155
156
157
158
159
160
        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)
        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))
161

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

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
166
167
168
169
170
        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
171
172

        Raises:
Matthieu Constant's avatar
Matthieu Constant committed
173
174
        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
175
        NameError -- If the dimension of the velocity vector is not equal to the physical dimension of the problem
176
        """
Matthieu Constant's avatar
Matthieu Constant committed
177
178
        if porosity < 1e-3 :
            raise NameError("Inflow porosity too small!")
179
180
        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")
181
182
        if (velocity is not None and len(velocity) != self._dim)  :
            raise NameError("len(velocity) != dimension at open boundaries")
183
        if velocity is not None :
Matthieu Constant's avatar
Matthieu Constant committed
184
            cb_or_value = velocity+[porosity]
185
            vid = 0
Matthieu Constant's avatar
Matthieu Constant committed
186
            cid = self._dim
187
            pid = -1
Matthieu Constant's avatar
Matthieu Constant committed
188
            n = self._dim+1
189
        else :
Matthieu Constant's avatar
Matthieu Constant committed
190
            cb_or_value = [pressure]+[porosity]
191
            pid = 0
Matthieu Constant's avatar
Matthieu Constant committed
192
            cid = 1
193
            vid = -1
Matthieu Constant's avatar
Matthieu Constant committed
194
            n = 2
195
196
197
198
199
200
        if self._n_fluids == 2 :
            cb_or_value += [concentration]
            aid = n
        else :
            aid = -1

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
201
        bndcb = BNDCB(_Bnd(cb_or_value, self._dim).apply)
Matthieu Constant's avatar
Matthieu Constant committed
202
        self.weak_cb_ref.append(bndcb)
203
        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))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
204

205
    def set_strong_boundary(self, tag, field,callback_or_value, row=None) :
Nathan Coppin's avatar
Nathan Coppin committed
206
        """Sets the strong boundaries (=constrains) for the fluid problem.
Matthieu Constant's avatar
Matthieu Constant committed
207
208

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
209
210
211
212
        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
213
        """
214
215
        if row is None :
            row = field
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
216
        bndcb = BNDCB(_Bnd([callback_or_value], self._dim).apply)
217
        self.strong_cb_ref.append(bndcb)
218
        self._lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), c_int(field), c_int(row), bndcb)
219

Matthieu Constant's avatar
Matthieu Constant committed
220
    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
221
        """Adapts the mesh.
Matthieu Constant's avatar
Matthieu Constant committed
222
223

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
224
225
226
227
228
229
230
231
        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
232
        """
233
        self._lib.fluid_problem_adapt_mesh(self._fp, 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
234
        self.sys = None
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
235

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
236
237
238
239
240
241
242
243
244
245
246
247
    def _mesh_boundaries(self) :
        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

Matthieu Constant's avatar
Matthieu Constant committed
248
    def export_vtk(self, output_dir, t, it,stab=False) :
Nathan Coppin's avatar
Nathan Coppin committed
249
        """Writes output file for post-visualization.
Matthieu Constant's avatar
Matthieu Constant committed
250
251

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
252
253
254
255
        output_dir -- Output directory
        t -- Computation time
        it -- Computation iteration
        stab -- If True exports the stabilisation parametres in the output files
Matthieu Constant's avatar
Matthieu Constant committed
256
        """
257
        v = self.solution()[:,:self._dim]
Matthieu Constant's avatar
Matthieu Constant committed
258
        da = self.concentration_dg_grad()
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
259
        cell_data = []
260
        if self._dim == 2 :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
261
            v = np.insert(v,self._dim,0,axis=1)
Matthieu Constant's avatar
Matthieu Constant committed
262
            da = np.insert(da,self._dim,0,axis=1)
263
264
        data = [
            ("pressure",self.solution()[:,[self._dim]]),
265
            ("velocity",v),
266
            ("porosity",self.porosity()),
Matthieu Constant's avatar
Matthieu Constant committed
267
268
            ("old_porosity",self.old_porosity()),
            ("grad",da)
269
            ]
270
        cell_data.append(("curvature",self.curvature()))
271
        if self._n_fluids == 2 :
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
272
            cell_data.append(("concentration",self.concentration_dg()))
Matthieu Constant's avatar
Matthieu Constant committed
273
            cell_data.append(("stab",self._get_matrix("stab_param",self.n_elements(),1)))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
274
        field_data = [(b"Boundary %s"%(name), nodes) for name,nodes in self._mesh_boundaries().items()]
275
        connectivity = self.elements()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
276
        types = np.repeat([5 if self._dim == 2 else 10],connectivity.shape[0])
277
        offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
278
        VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
279

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
280
    def import_vtk(self, filename) :
Nathan Coppin's avatar
Nathan Coppin committed
281
        """Reads output file to reload computation.
Matthieu Constant's avatar
Matthieu Constant committed
282
283

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
284
        filename -- Name of the file to read
Matthieu Constant's avatar
Matthieu Constant committed
285
        """
286
        x,cells,data,cdata,fdata = VTK.read(filename)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
287
288
        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()))
289
        el = cells["connectivity"].reshape([-1,self._dim+1])
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
290
291
292
293
294
        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")
        self._lib.fluid_problem_set_elements(self._fp,
295
296
                c_int(x.shape[0]),_np2c(x,np.float64),
                c_int(el.shape[0]),_np2c(el,np.int32),
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
297
298
                c_int(bnds.shape[0]),c_void_p(bnds.ctypes.data),c_void_p(bnd_tags.ctypes.data),c_int(len(cbnd_names)),cbnd_names
                )
299
        sol = self.solution()
300
        sol[:,:self._dim] = data["velocity"][:,:self._dim]
301
        sol[:,[self._dim]] = data["pressure"]
302
        if self._n_fluids == 2 :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
303
            self.concentration_dg()[:] = cdata["concentration"]
304
305
306
        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
307
        self.sys = None
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
308

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

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
313
        dt -- Computation time step
Matthieu Constant's avatar
Matthieu Constant committed
314
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
315
        forces = np.ndarray([self.n_particles,self._dim],"d",order="C")
316
        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
317
318
        return forces

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
319
    def compute_node_torque(self, dt) :
Nathan Coppin's avatar
Nathan Coppin committed
320
321
        """Computes the angular drags to apply on each grain resulting from the fluid motion.
        Only in 2D
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
322
        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
323
        dt -- Computation time step
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
324
325
326
327
328
329
        """
        torques = np.ndarray([self.n_particles,1],"d",order="C")
        
        self._lib.fluid_problem_compute_node_particle_torque(self._fp, c_double(dt), c_void_p(torques.ctypes.data))
        return torques

330

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
331
    def implicit_euler(self, dt, check_residual_norm=-1, reduced_gravity=0, stab_param=0.) :
Nathan Coppin's avatar
Nathan Coppin committed
332
        """Solves the fluid equations.
Matthieu Constant's avatar
Matthieu Constant committed
333
334

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
335
336
337
338
        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
339
        """
340
341
        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
342
343
344
345
346
        self._lib.fluid_problem_apply_boundary_conditions(self._fp)
        solution_old = self.solution().copy()
        rhs = np.zeros_like(solution_old)
        self._lib.fluid_problem_reset_boundary_force(self._fp)
        self._lib.fluid_problem_compute_stab_parameters(self._fp,c_double(dt))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
347
348
349
350
        if self.sys is None :
            self.sys = LinearSystem(self.elements(),self.n_fields(),self.solver_options
                    )
        sys = self.sys
351
352
353
354
        localv = np.ndarray([sys.localsize*self.n_elements()])
        localm = np.ndarray([sys.localsize**2*self.n_elements()])
        self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
        sys.local_to_global(localv,localm,rhs)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
355
        if check_residual_norm > 0 :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
356
            norm0 = np.linalg.norm(rhs)
357
        self.solution()[:] -= sys.solve(rhs)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
358
        if check_residual_norm > 0 :
359
            sys.zero_matrix()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
360
            rhs[:] = 0
361
362
            self._lib.fluid_problem_assemble_system(self._fp,_np2c(rhs),_np2c(solution_old),c_double(dt),_np2c(localv),_np2c(localm))
            sys.local_to_global(localv,localm,rhs)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
363
364
            norm = np.linalg.norm(rhs)
            print("norm",norm)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
365
366
            if norm > check_residual_norm :
                raise ValueError("wrong derivative or linear system precision")
367
        self._lib.fluid_problem_node_force_volume(self._fp,_np2c(solution_old),c_double(dt),None,None)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
368

Matthieu Constant's avatar
Matthieu Constant committed
369
    def compute_cfl(self):
Nathan Coppin's avatar
Nathan Coppin committed
370
        """Computes the CFL number divided by the time step
Matthieu Constant's avatar
Matthieu Constant committed
371
372
373
374
375
376
377
378

        """
        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
379
    def advance_concentration(self,dt):
Nathan Coppin's avatar
Nathan Coppin committed
380
        """Solves the advection equation for the concentration using the current velocity field.
Matthieu Constant's avatar
Matthieu Constant committed
381
382

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
383
        dt -- Computation time step
Matthieu Constant's avatar
Matthieu Constant committed
384
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
385
        if self._n_fluids == 2 :
Matthieu Constant's avatar
Matthieu Constant committed
386
387
            cfl = self.compute_cfl()
            nsub = int(cfl*dt+1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
388
389
390
391
            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
392

393
    def set_particles(self, mass, volume, position, velocity, contact, init=False, reload = False):
Matthieu Constant's avatar
Matthieu Constant committed
394
395
396
        """Set location of the grains in the mesh and compute the porosity in each cell.

        Keyword arguments:
Nathan Coppin's avatar
Nathan Coppin committed
397
398
399
400
401
402
        mass -- List of particles mass
        volume -- List of particles volume
        position -- List of particles centre positions 
        velocity -- List of particles velocity
        init -- if True sets the old_porosity to zero (has to be used at the first call)
        reload -- Optional arguments specifying if the simulation is reloaded or if it is the first iteration
Matthieu Constant's avatar
Matthieu Constant committed
403
        """
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
404
        self.n_particles = mass.shape[0]
405
        self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),_np2c(mass),_np2c(volume),_np2c(position),_np2c(velocity),_np2c(contact),None,c_int(init),c_int(reload))
406

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
407

408
    def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :
409
        f = getattr(self._lib,"fluid_problem_"+f_name)
410
        f.restype = POINTER(typec)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
411
        return np.ctypeslib.as_array(f(self._fp),shape=(nrow,ncol))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
412
413

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
417
    def velocity(self) :
Nathan Coppin's avatar
Nathan Coppin committed
418
        """Gives access to the fluid velocity value at the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
419
420
421
        return self.solution()[:,:-1]

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

425
426
427
428
    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
429
    def coordinates(self) :
Nathan Coppin's avatar
Nathan Coppin committed
430
        """Gives access to the coordinate of the mesh nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
431
        return self._get_matrix("coordinates",self.n_nodes(), 3)
432

433
    def n_fluids(self) :
Nathan Coppin's avatar
Nathan Coppin committed
434
        """Returns the number of fluids."""
435
436
        return self._n_fluids

437
    def elements(self):
Nathan Coppin's avatar
Nathan Coppin committed
438
        """Gives read-only access to the elements of the mesh."""
439
440
441
        return self._get_matrix("elements", self.n_elements(), self._dim+1,c_int)

    def n_elements(self):
Nathan Coppin's avatar
Nathan Coppin committed
442
        """Returns the number of mesh nodes."""
443
444
        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
445
446
    
    def n_fields(self):
Nathan Coppin's avatar
Nathan Coppin committed
447
        """Returns the number of fluid fields."""
448
449
        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
450
451

    def n_nodes(self):
Nathan Coppin's avatar
Nathan Coppin committed
452
        """Returns the number of mesh nodes."""
453
454
        self._lib.fluid_problem_n_nodes.restype = c_int
        return self._lib.fluid_problem_n_nodes(self._fp);
455
456
    
    def porosity(self):
Nathan Coppin's avatar
Nathan Coppin committed
457
        """Returns the porosity at nodes"""
458
459
        return self._get_matrix("porosity", self.n_nodes(), 1)

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
460
    def set_concentration_cg(self,concentration):
Nathan Coppin's avatar
Nathan Coppin committed
461
        """Sets the concentration at nodes"""
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
462
        concentration = concentration.reshape((self.n_nodes(),1))
463
        self._lib.fluid_problem_set_concentration_cg(self._fp,_np2c(concentration))
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
464
465

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

Matthieu Constant's avatar
Matthieu Constant committed
469
    def concentration_dg_grad(self):
Nathan Coppin's avatar
Nathan Coppin committed
470
        """Returns the concentration at discontinuous nodes"""
Matthieu Constant's avatar
Matthieu Constant committed
471
472
        return self._get_matrix("concentration_dg_grad", self.n_nodes(), self._dim)

473
    def curvature(self):
Nathan Coppin's avatar
Nathan Coppin committed
474
        """Returns the porosity at previous time step"""
475
476
        return self._get_matrix("curvature", self.n_elements(), 1)

477
    def old_porosity(self):
Nathan Coppin's avatar
Nathan Coppin committed
478
        """Returns the porosity at previous time step"""
479
        return self._get_matrix("old_porosity", self.n_nodes(), 1)
480
481
    
    def volume_flux(self, bnd_tag):
Nathan Coppin's avatar
Nathan Coppin committed
482
        """Computes the integral of the (outgoing) normal velocity through boundary with tag bnd_tag"""
483
484
485
        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
486
    def particle_element_id(self) :
Nathan Coppin's avatar
Nathan Coppin committed
487
        """Returns the id of the mesh cell in which particles are located."""
488
        f = getattr(self._lib,"fluid_problem_particle_element_id")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
489
        f.restype = c_void_p
490
        fs = getattr(self._lib,"fluid_problem_n_particles")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
491
492
493
494
        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
495
        return np.ctypeslib.as_array(buf)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
496

497
    def particle_uvw(self) :
Nathan Coppin's avatar
Nathan Coppin committed
498
        """Returns the coordinates of the particles inside their element"""
499
500
        return self._get_matrix("particle_uvw",self.n_particles,self._dim)

Matthieu Constant's avatar
Matthieu Constant committed
501
    def particle_position(self) :
Nathan Coppin's avatar
Nathan Coppin committed
502
        """Gives access to the fluid field value at the mesh nodes."""
Matthieu Constant's avatar
Matthieu Constant committed
503
504
505
        return self._get_matrix("particle_position", self.n_particles, self._dim)

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

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
509
    def bulk_force(self) :
Nathan Coppin's avatar
Nathan Coppin committed
510
        """Gives access to the volume force at fluid nodes."""
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
511
512
        return self._get_matrix("bulk_force", self.n_nodes(), self._dim)

513
    def node_volume(self) :
Nathan Coppin's avatar
Nathan Coppin committed
514
        """Returns the volume associated with each node"""
515
516
        return self._get_matrix("node_volume",self.n_nodes(),1)

517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
    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))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
532
533
        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)
534
        return phys_name.value,phys_dim.value,phys_id.value,nodes,bnd
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
535
536
537

    def dimension(self) :
        return self._dim