Brinkman.py 9.15 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
21
22
23
24
25
26
27
28
import numpy as np
from . import quadrature

_u_ = 0
_v_ = 1
_w_ = 2

def derivate(f, x, nd = 1, eps = 1e-8) :
    f0 = f(x)
    xd = x.reshape(x.shape[:x.ndim-nd]+(-1,))
    d = np.ndarray((f0.shape +(xd.shape[-1],)))
    for i in range(xd.shape[-1]) :
        xd[..., i] += eps
        d[..., i] = (f(xd.reshape(x.shape)) - f0) / eps
        xd[..., i] -= eps
    return d.reshape(f0.shape + x.shape[x.ndim - nd:])


class Brinkman :

    def dot(self, a, b) :
        return a[..., 0] * b[..., 0] + a[..., 1] * b[..., 1] + (0 if self.D == 2 else a[..., 2] * b[..., 2])

    def psiTerm0(self, sol, solgrad):
        _p_ = self.D
        _U_ = range(self.D)
        val = np.zeros_like(sol)
        #self.g[None, None, :] * self.ca * self.rhof
Vincent Legat's avatar
Vincent Legat committed
29
30
        # c=None == schema implicite classique
        # sinon compression artificielle
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
        nu = self.mu/self.rhof
Vincent Legat's avatar
Vincent Legat committed
32
        if self.c is not None :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
33
            val[..., _U_] = - solgrad[..., _p_, :]/self.rhof \
matthieu's avatar
matthieu committed
34
                - nu * (self.dot(solgrad[...,_U_,:], self.dca)-sol[...,_U_]/self.ca**2*self.dot(self.dca, self.dca))
Vincent Legat's avatar
Vincent Legat committed
35
            val[..., _p_] = - self.c * (solgrad[..., _u_, 0] + solgrad[..., _v_, 1] + (0 if self.D == 2 else solgrad[..., _w_, 2]))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
36
            val[..., _p_] += - self.c * (self.dvs[..., 0, 0] + self.dvs[..., 1, 1] + (0 if self.D == 2 else self.dvs[..., 2, 2]))
Vincent Legat's avatar
Vincent Legat committed
37
        else :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
38
39
            val[..., _U_] =   - (sol[...,_U_] - self.sol0qp[...,_U_]) * self.rhof/self.dt - solgrad[..., _p_, :] \
                + self.mu * (self.dot(solgrad[...,_U_,:], self.dca)-sol[...,_U_]/self.ca**2*self.dot(self.dca, self.dca))
Vincent Legat's avatar
Vincent Legat committed
40
41
            val[..., _p_] = - solgrad[..., _u_, 0] - solgrad[..., _v_, 1] - (0 if self.D == 2 else solgrad[..., _w_, 2])
            val[..., _p_] += - (self.dvs[..., 0, 0] + self.dvs[..., 1, 1] + (0 if self.D == 2 else self.dvs[..., 2, 2]))
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
42
43
44
45
46
47
48
        return val

    def psiTerm1(self, sol, solgrad):
        _p_ = self.D
        _U_ = range(self.D)

        val = np.ndarray(sol.shape + (self.D,))
49
        nu = self.mu/self.rhof
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
50
51
        val[..., _U_, :] = -solgrad[..., _U_, :] * nu
        val[..., _p_, :] =  0 #self.c * (self.vs + sol[..., _U_])
52
        val[..., _p_, :] = -self.epsilon * solgrad[..., _p_, :]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
53
        if self.ns :
54
55
            val[..., _u_, :] += sol[..., [_u_]] * sol[..., _U_] / self.ca
            val[..., _v_, :] += sol[..., [_v_]] * sol[..., _U_] / self.ca
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
56
            if self.D == 3:
57
                val[..., _w_, :] += sol[..., [_w_]] * sol[..., _U_] / self.ca
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
58
59
60
61
62
63
64
65
        return val

    # forces by the particle ON the fluid
    def particleForces(self, sol, solgrad, fonpart = None) :
        _p_ = self.D
        _U_ = range(self.D)
        val = np.zeros_like(sol)
        if self.D == 2 :
66
            #diametre
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
67
            d = 2 * (self.pf.volume[0, 0] / np.pi)**0.5
68
69
            #aire projetee sur le plan perpendiculaire a l'ecoulement
            Ap=d
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
70
71
        else :
            d = 2 * (3./(4 * np.pi) * self.pf.volume[0, 0])**(1./3)
72
73
74
            #aire projetee sur le plan perpendiculaire a l'ecoulement
            Ap=np.pi*(d/2)**2
        #difference des vitesses solide et fluide u_p-u_f
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
75
76
        U = self.vs - sol[..., _U_]/self.ca
        normU = np.sqrt(self.dot(U, U))[...,np.newaxis]
77
        #Reynolds/|u_p-u_f|
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
78
        Re_O_u = d * self.ca/self.mu * self.rhof
79
80
        #coefficient du drag
        Cd_u = (0.63 * normU**0.5 + 4.8/(Re_O_u**0.5))**2
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
81
82
        f = self.ca**(-1.8)
        F = solgrad[..., _p_, :]
83
84
        #drag
        GoU = f * Cd_u * self.rhof / 2 * Ap
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
85
        reaction = (F * self.pf.volume + (GoU*U) / (1 + self.stab/self.pf.mass * GoU))*self.dragFactor
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
86
        if fonpart is not None:
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
87
            fonpart[...] = -reaction - self.g * self.pf.volume * self.rhof
matthieu's avatar
matthieu committed
88
        val[..., _U_] =  reaction / self.rhof 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
89
90
91
        val[..., _p_] = 0
        return val

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
92
    def __init__(self, mu, pf, g, rho, rhop, epsilon, ns = False, stab = 0, D = 2, c=1, dragFactor=0.3) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
93
        self.D = D
Vincent Legat's avatar
Vincent Legat committed
94
        self.c = c
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
95
        self.ns = ns
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
96
        self.dragFactor = dragFactor
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
97
98
99
100
101
102
        self.rhof = rho
        self.rhop = rhop
        self.stab = stab
        self.mu = mu
        self.pf = pf
        self.g = np.array(g)
103
        self.epsilon = epsilon
matthieu's avatar
matthieu committed
104
        self.Ecf = 0
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
105
106
107
108
109
110
111
112
113
114
115
116

    def particleInteractionTerm(self, meshJ, solution, res, jac) :
        _p_ = self.D
        _U_ = range(self.D)
        dofm = solution.dof_manager()
        self.ca = self.pf.porosity.at_points(self.pf.eid, self.pf.uvw)
        self.vs = self.pf.velocity
        sol = solution.at_points(self.pf.eid, self.pf.uvw)
        gradsol = solution.gradient_at_points(self.pf.eid, self.pf.uvw, meshJ)
        pTerm = self.particleForces(sol, gradsol, self.pf.forces)
        for i in range(dofm.n_fields()) :
            np.add.at(res, dofm._fields[i][0][:, self.pf.eid], (dofm.get_field_basis(i).f(self.pf.uvw)*pTerm[:, [i]]).T)
Vincent Legat's avatar
Vincent Legat committed
117
118
        if  jac is None :
            return
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
        d00 = derivate(lambda s : self.particleForces(s, gradsol), sol, 1)
        d01 = derivate(lambda s : self.particleForces(sol, s), gradsol, 2)
        for f in range(dofm.n_fields()) :
            ff = dofm.get_field_basis(f).f(self.pf.uvw)
            for g in range(dofm.n_fields()) :
                pp = d00[:, f, g, None] * dofm.get_field_basis(g).f(self.pf.uvw)
                d01xi = meshJ.multDXiDXLeft(d01[:, f, g, :], 1, self.pf.eid)
                pp += (d01xi[:, None, :] * dofm.get_field_basis(g).df(self.pf.uvw)).sum(2)
                jac.add_to_field(f, g, pp[:, None, :] * ff[:, :, None], self.pf.eid)

    def fluidTerm(self, meshJ, solution, res, jac) :
        dofm = solution.dof_manager()
        qp, qw = quadrature.get_triangle(3) if self.D == 2 else quadrature.get_tetrahedron(3)
        self.vs = self.pf.momentum.at_qp(qp)
        self.dvs = self.pf.momentum.gradient_at_qp(qp, meshJ)
        self.ca = self.pf.porosity.at_qp(qp)
        self.dca = self.pf.porosity.gradient_at_qp(qp, meshJ)
        sol = solution.at_qp(qp)
        gradsol = solution.gradient_at_qp(qp, meshJ)
Vincent Legat's avatar
Vincent Legat committed
138
139
140
       
        if self.c is None :
            self.sol0qp = self.sol0.at_qp(qp)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
141
142
143
144
145
146
147
148
149
150
151
152
        #f = meshJ.dofManager.get_field_basis(0).f(qp)
        #self.xqp = np.tensordot(self.x,f,[1,1]).swapaxes(1,2)

        x0 = meshJ.multJ(self.psiTerm0(sol, gradsol))
        x1 = meshJ.multJ(meshJ.multDXiDXLeft(self.psiTerm1(sol, gradsol), 3))
        for i in range(dofm.n_fields()) :
            basis = dofm.get_field_basis(i)
            psiW = basis.f(qp) * qw[:, None]
            r0 = np.tensordot(psiW, x0[:, :, i], [0, 1])
            dPsiW = basis.df(qp) * qw[:, None, None]
            r1 = np.tensordot(dPsiW, x1[:, :, i, :], [[0, 2], [1, 2]])
            np.add.at(res, dofm._fields[i][0], r0 + r1)
matthieu's avatar
matthieu committed
153
154
155
156
157
158
159
160
        
        #test pour l'energie
        _U_ = range(self.D)
        solu = solution.at_qp(qp)
        vL = (solu[...,_u_]**2 +solu[...,_v_]**2)* 0.5
        Vg = meshJ.multJ(vL)
        Ecf1 = np.tensordot(qw,Vg[:,:],[0,1]) * self.rhof
        self.Ecf = sum(Ecf1)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
161

Vincent Legat's avatar
Vincent Legat committed
162
163
        if jac is None: 
            return
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
        d00 = derivate(lambda s : self.psiTerm0(s, gradsol), sol, 1)
        d01 = derivate(lambda s : self.psiTerm0(sol, s), gradsol, 2)
        d10 = derivate(lambda s : self.psiTerm1(s, gradsol) , sol, 1)
        d11 = derivate(lambda s : self.psiTerm1(sol, s), gradsol, 2)
        for f in range(dofm.n_fields()):
            basisf = dofm.get_field_basis(f)
            for g in range(dofm.n_fields()):
                basisg = dofm.get_field_basis(g)
                if (np.max(np.abs(d00[...,f,g])) > 1e-8) :
                    psiPsiW = np.einsum("qn, qm, q-> qnm", basisf.f(qp), basisg.f(qp), qw)
                    jac.add_to_field(f, g, np.tensordot(meshJ.multJ(d00[:,:,f,g]), psiPsiW, 1))
                if (np.max(np.abs(d01[:,:,f,g,:])) > 1e-8) :
                    dPsiPsiW = np.einsum("qmd, qn, q -> qdnm", basisg.df(qp), basisf.f(qp), qw)
                    jac.add_to_field(f, g, np.tensordot(meshJ.multJ(meshJ.multDXiDXLeft(d01[:,:,f,g,:], 2)), dPsiPsiW, 2))
                if (np.max(np.abs(d10[:,:,f,:,g])) > 1e-8) :
                    psiDPsiW = np.einsum("qmd, qn, q -> qdmn", basisf.df(qp), basisg.f(qp), qw)
                    jac.add_to_field(f, g, np.tensordot(meshJ.multJ(meshJ.multDXiDXLeft(d10[:,:,f,:,g], 2)), psiDPsiW, 2))
                if (np.max(np.abs(d11[:,:,f,:,g,:])) > 1e-8) :
                    dPsiDPsiW = np.einsum("qmx, qny, q -> qxymn", basisf.df(qp), basisg.df(qp), qw)
                    jac.add_to_field(f, g, np.tensordot(meshJ.multDXiDXLeft(meshJ.multDXiDXLeft(meshJ.multJ(d11[:,:,f,:,g,:]), 3), 2), dPsiDPsiW, 3))
matthieu's avatar
matthieu committed
184
185
186
187
188
189
190
        
        #test pour l'energie
        sol = solution.at_qp(qp)
        vL = sol[...,_U_]**2 * 0.5
        Vg = meshJ.multJ(vL)
        self.Ecf = np.tensordot(qw,Vg[:,:],[0,1]) * self.rhof
        print("ECFFFF",self.Ecf)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
191
192
193
194

    def massTerm(self, alpha, solution, res, jac) :
        pass

Vincent Legat's avatar
Vincent Legat committed
195
    def volumeTerm(self, meshJ, solution, jacNeeded=True) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
196
        dofm = solution.dof_manager()
Vincent Legat's avatar
Vincent Legat committed
197
        jac = dofm.new_matrix() if jacNeeded else None
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
198
        res = dofm.new_vector()
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
199
200
        if self.pf.eid.shape[0] > 0 :
            self.particleInteractionTerm(meshJ, solution, res, jac)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
201
202
        self.fluidTerm(meshJ, solution, res, jac)
        #self.massTerm(meshJ, alphamass, solution, res, jac)
Vincent Legat's avatar
Vincent Legat committed
203
204
        if jac is None:
            return res
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
205
        return res, jac.to_csr_matrix()