time_integration.py 9 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# MigFlow - Copyright (C) <2010-2018>
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	
# List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# 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
# GNU Lesser General Public License for more details.
# 
# 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, 
# see <http://www.gnu.org/licenses/>.
21
22
import numpy as np

Nathan Coppin's avatar
Nathan Coppin committed
23
def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_sub_iter=None,max_split=1) :
24
25
26
27
28
29
30
31
32
    """Contact solver for the grains

    Keyword arguments:
    particles -- particles structure
    f -- external forces on particles
    dt -- time step
    min_nsub -- minimal nsub for the particles iterations
    contact_tol -- tolerance for the contact solver
    iteration -- internal argument -- used to keep track of the number of time step splitting
Nathan Coppin's avatar
Nathan Coppin committed
33
34
    after_sub_iter -- callback to execute once a sub iteration has been made
    max_split -- maximum number of times the time step can be further split if conergence is not reached
35
    """
36
    # Compute free velocities of the grains
37
    v = particles.velocity()
38
    if f is None:
39
        f = np.zeros_like(v)
40
    # Estimation of the solid sub-time steps
Nathan Coppin's avatar
Nathan Coppin committed
41
    if particles.r() is not None and particles.r().size != 0:
42
      vn = v + f*dt / particles.mass()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
43
      vmax = np.max(np.linalg.norm(vn,axis=1))
44
45
46
      nsub = max(min_nsub, int(np.ceil((vmax * dt * 8)/min(particles.r()))))
    else:
      nsub = 1
Nathan Coppin's avatar
Nathan Coppin committed
47
48
      vmax = 0
    print("NSUB", nsub,"VMAX",vmax, "VMAX * dt", vmax * dt, "r", min(particles.r()) if (particles.r() is not None and particles.r().size !=0) else 0 )
49
    #If time step was split too many times, ignore the convergence and move the grains
Nathan Coppin's avatar
Nathan Coppin committed
50
    if iteration == max_split:
51
52
      for i in range(nsub) :
        particles.iterate(dt/nsub, f, tol=contact_tol,force_motion=1)
53
        if after_sub_iter :
54
            after_sub_iter(nsub)
55
56
      return
    # For each sub-time step solve contacts (do NLGS iterations) and move the grains 
57
    for i in range(nsub) :
58
59
        if not particles.iterate(dt/nsub, f, tol=contact_tol):
          print("Splitting time-step to level %d..."%(iteration+1))
60
          _advance_particles(particles,f,dt/nsub,2,contact_tol/2,iteration=iteration+1,after_sub_iter=after_sub_iter)
61
          print("Convergence reached for subinvervals of level %d"%(iteration+1))
62
63
        else :
            if after_sub_iter :
64
                after_sub_iter(nsub)
65

66
      
Nathan Coppin's avatar
Nathan Coppin committed
67
def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, alpha=0.5,after_sub_iter=None,max_split=1) :  
68
69
70
71
72
73
74
75
76
77
    """Predictor-corrector scheme to solve fluid and grains.

    Keyword arguments:
    fluid -- fluid structure
    particles -- particles structure
    dt -- time step
    min_nsub -- minimal nsub for the particles iterations
    contact_tol -- tolerance for the contact solver
    external_particles_forces -- external forces applied on the particles
    alpha -- parametre of the predictor-corrector scheme [alpha*f(n)+(1-alpha)*f(n+1)]
Nathan Coppin's avatar
Nathan Coppin committed
78
79
    after_sub_iter -- callback to execute once a sub iteration has been made
    max_split -- maximum number of times the time step can be further split if conergence is not reached
80
    """
81
    particles.save_state()
82
    if external_particles_forces is None:
83
        external_particles_forces = np.zeros_like(particles.velocity())
84
85
86
    # Copy the fluid solution of the previous time step before computing the prediction 
    sf0 = np.copy(fluid.solution())
    # Copy the solid solution of the previous time step before computing the prediction
87
    x0 = np.copy(particles.position())
88
    cp0 = np.copy(particles.contact_forces())
89
90
    seg0 = np.copy(particles.segments())
    disk0 = np.copy(particles.disks())
Nathan Coppin's avatar
Nathan Coppin committed
91
    if particles.dim() == 3 :
92
93
        tri0 = np.copy(particles.triangles())

94
95
96
97
98
99

    # Predictor
    #
    # Compute the fluid solution and keep a copy of this solution and the forces applied by the fluid on the grains
    fluid.implicit_euler(dt)
    sf1 = np.copy(fluid.solution())
100
    f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
Nathan Coppin's avatar
Nathan Coppin committed
101
    _advance_particles(particles,f0,dt,min_nsub,contact_tol,max_split=max_split)
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
102
    particles.contact_forces()[:] = cp0
103

104
    # Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
105
    fluid.set_particles(particles.mass(), particles.volume(), x0, particles.velocity(), particles.contact_forces(), reload=True)
106
107
108
109
110
111
112
113
114

    # Corrector
    #
    # Compute the fluid solution from the previous time-step one using the predicted solid velocities in the fluid-grain interaction force
    fluid.solution()[:,:] = sf0
    fluid.implicit_euler(dt)
    # Fluid solution is a weighted average of the predicted one and the corrected one.
    # The alpha parametre is the weight
    fluid.solution()[:,:] = (alpha*fluid.solution()+(1-alpha)*sf1)
115
    f1 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
116
117
    # For two fluids flows, advance the concentration field with the fluid velocity.
    # The number of sub-time steps for the advection is automatically computed.
118
    if fluid.n_fluids() == 2 :
119
120
        fluid.advance_concentration(dt)
    # Reset solid velocities
121
    particles.restore_state()
122
123
    particles.segments()[:] = seg0
    particles.disks()[:] = disk0
Nathan Coppin's avatar
Nathan Coppin committed
124
    if particles.dim() == 3 :
125
        particles.triangles()[:] = tri0
Nathan Coppin's avatar
Nathan Coppin committed
126
    _advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
127
    # Give to the fluid the solid information
128
    fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces()) 
129

Nathan Coppin's avatar
Nathan Coppin committed
130
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_split=1) :  
131
132
133
134
135
136
137
138
139
140
    """Migflow solver: the solver type depends on the fluid and particles given as arguments.

    Keyword arguments:
    fluid -- fluid structure
    particles -- particles structure
    dt -- time step
    min_nsub -- minimal nsub for the particles iterations
    contact_tol -- tolerance for the contact solver
    external_particles_forces -- vector of external forces applied on the particles
    fixed_grains -- boolean variable specifying if the grains are fixed in the fluid
Nathan Coppin's avatar
Nathan Coppin committed
141
142
    after_sub_iter -- callback to execute once a sub iteration has been made
    max_split -- maximum number of times the time step can be further split if conergence is not reached
143

144
145
146
147
148
149
    Raises:
        ValueError -- fluid and particles cannot be both None
        ValueError -- external_particles_forces must have shape (number of particles,dimension)
    """
    if particles is None and fluid is None:
        raise ValueError("fluid and particles are both None: not possible to iterate!")
Matthieu Constant's avatar
Matthieu Constant committed
150
    if particles is not None and external_particles_forces is not None:
Jonathan Lambrechts's avatar
fix    
Jonathan Lambrechts committed
151
        if external_particles_forces.shape!=(particles.n_particles(),particles.dim()):
152
            raise ValueError("external_particles_forces must have shape (number of particles,dimension!")
Nathan Coppin's avatar
Nathan Coppin committed
153
    if (particles is None or particles.r() is None or particles.r().size == 0) and fluid is not None:
154
        fluid.implicit_euler(dt)
155
156
157
158
        # For two fluids flows, advance the concentration field with the fluid velocity.
        # The number of sub-time steps for the advection is automatically computed.
        if fluid.n_fluids() == 2 :
            fluid.advance_concentration(dt)
159
    if particles is not None and fluid is None:
Nathan Coppin's avatar
Nathan Coppin committed
160
        _advance_particles(particles,external_particles_forces,dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
Nathan Coppin's avatar
Nathan Coppin committed
161
    if fluid is not None and particles is not None and particles.r() is not None and particles.r().size !=0:
162
163
        if fixed_grains:
            fluid.implicit_euler(dt)
164
            fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
165
        else:
Nathan Coppin's avatar
Nathan Coppin committed
166
            predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_split=max_split)
167
168

def particle_changed(fluid,p) :
169
170
    """ Modifies the porosity and fluid velocity fields when particles have been removed or added in the fluid domain
    """
171
172
    co = fluid.old_porosity().copy()
    c = fluid.porosity().copy()
173
    fluid.set_particles(p.mass(),p.volume(), p.position(), p.velocity(), p.contact_forces())
174
175
    c2 = fluid.porosity()
    fluid.old_porosity()[:] = co+c2-c
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
176
    fluid.velocity()[:] *= c2/c
177