fluidSolver.py 7.92 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
3
4
5
from .meshJacobian import MeshJacobian
from .particleFluid import ParticleFluid
from .Brinkman import Brinkman
from .dofManager import DofManager
from .Steady import Steady
Vincent Legat's avatar
Vincent Legat committed
6
from .ExplicitSolver import ExplicitSolver
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
7
8
9
10
11
from .mesh import Mesh
from .legendre import *
import time

class FluidSolver :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
12
    def __init__(self, mesh, dt, rhop, rho, mu, c, g, epsilon, TypeEMu, TypeEMp, dragFactor=0.3,timeScheme="Implicit", imex=True) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
13
14
15
16
17
18
        if isinstance(mesh, str):
            mesh = Mesh(mesh)
        self._mesh = mesh
        self._jac = MeshJacobian(mesh)
        self._pf = ParticleFluid(self._jac)
        self._imex = imex
Vincent Legat's avatar
Vincent Legat committed
19
        self._timeScheme = timeScheme
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
22
        self._dt = dt
        self._stabFactor = 1
        law = Brinkman(
Vincent Legat's avatar
Vincent Legat committed
23
            c =c,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
24
25
            mu = mu,
            pf = self._pf,
26
            ns = True,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
27
28
29
            g = np.array(g),
            rho = rho,
            rhop = rhop,
30
            epsilon = epsilon,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
31
            dragFactor=dragFactor,
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
32
33
34
35
            D = self._jac.dim,
            stab = dt)
        d = DofManager(mesh)
        if self._jac.dim == 2 :
36
37
38
            d.add_field(TypeEMu)
            d.add_field(TypeEMu)
            d.add_field(TypeEMp)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
39
        elif self._jac.dim == 3:
40
41
42
43
            d.add_field(TypeEMu)#TetrahedronP1Mini
            d.add_field(TypeEMu)
            d.add_field(TypeEMu)
            d.add_field(TypeEMp)#TetrahedronP1
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
        else :
            raise ValueError("invalid mesh dimension : %i" % self._jac.dim)
        self._dt = dt
        self._law = law
        law.dt = dt
        law.x = self._jac.x
        self._dof = d
        self._sol = None

    def add_boundary_condition(self, tag, field, *a, **b) :
        self._dof.add_boundary_condition(tag, field, *a, **b)

    def complete(self) :
        self._dof.complete()
        self._sol = self._dof.new_vector()
        self._sol[:] = 0
Vincent Legat's avatar
Vincent Legat committed
60
61
62
63
64
65
66
67
68
        if self._timeScheme == "Explicit" :
            self._solver = ExplicitSolver(self._law, self._dof, self._jac)
        elif self._timeScheme == "Implicit" :
            self._law.sol0 = self._dof.new_vector()
            self._law.sol0[:] = 0
            self._solver = Steady(self._law, self._dof, self._jac)
            self._solverEx = Steady(self._law, self._dof, self._jac)
        else :
            print("Error : unknown time scheme :-)")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
69
70
71

    def solve(self) :
        tic = time.clock()
Vincent Legat's avatar
Vincent Legat committed
72
        if self._timeScheme == "Explicit" :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
73
            self._law.stab = self._dt
Vincent Legat's avatar
Vincent Legat committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
            self._solver.solve(self._sol, self._dt)
        elif self._timeScheme == "Implicit" :
            if self._imex :
                v0 = self._pf.velocity.copy()
                self._law.stab = self._dt *0.5 * self._stabFactor
                self._sol[:] = self._solver.solve()
                self._pf.velocity += self._pf.forces * self._dt*0.5 / self._pf.mass
                self._law.stab = 0
                self._sol[:] = self._solverEx.solve()
                self._pf.velocity[:] = v0
            else :
                self._law.stab = self._dt * self._stabFactor
                self._sol[:] = self._solver.solve()
            self._law.sol0[:] = self._sol
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
88
        else :
Vincent Legat's avatar
Vincent Legat committed
89
            print("Error : unknown time scheme :-)")
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
90
91

        print("cpu : ", time.clock() - tic)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
92
93
94
        self.forces[:] = 0
        self.forces[self._pf.mask,:] = self._pf.forces
        return self.forces
Vincent Legat's avatar
Vincent Legat committed
95
       
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
    def set_stab_factor(self, s) :
        self._stabFactor = s

    def write_position(self, odir, oiter, time, filetype="msh"):
        if filetype == "msh" :
            basename = "%s/f_%05i" % (odir, oiter)
            if self._jac.dim == 2 :
                X = np.zeros(((self._jac.dim + 1) *3, self._jac.x.shape[0]))
                X[[0,1,3,4,6,7], :] = self._jac.x.reshape([-1, 6]).T
            else :
                X = self._jac.x.reshape(-1, 12).T
            self._jac.writeScalar(basename + "-xmesh.msh", "x", oiter, time, X)
        else :
            print("unknown export format : \"%s\"." % filetype)

    def write_solution(self, odir, oiter, time, filetype="msh") :
        dim = self._jac.dim
        if filetype == "msh" :
            basename = "%s/f_%05i" % (odir, oiter)
            self._jac.writeScalar(basename + "-porosity.msh", "porosity", oiter, time, self._jac.dofManager.getField(0, self._pf.porosity))
            self._jac.writeScalar(basename + "-pressure.msh", "pressure", oiter, time, self._dof.getField(dim, self._sol))
            X = np.zeros(((dim+1)*3, self._jac.x.shape[0]))
            for d in range (dim) :
Vincent Legat's avatar
Vincent Legat committed
119
                X[d::3, :] = self._dof.getField(d, self._sol)[:,:]
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
            self._jac.writeScalar(basename + "-velocity.msh", "velocity", oiter, time, X)
        elif filetype == "vtk" :
            output = open("%s/f_%05i.vtu" %(odir, oiter), "w")
            vertices = np.array(self._mesh.vertices)[:,:3]
            elements = np.array([[v[3] for v in e.vertices] for e in self._dof.getElements(0)], np.int32) -1
            v = np.zeros((vertices.shape[0], 3))
            for d in range (dim) :
                v[elements, d] = self._jac.x[:,:,d]
            output.write("<VTKFile type=\"UnstructuredGrid\" format=\"0.1\">\n")
            output.write("<UnstructuredGrid>\n")
            output.write("<Piece NumberOfPoints=\"%i\" NumberOfCells=\"%i\">\n" % (vertices.shape[0], elements.shape[0]))
            output.write("<Points>\n")
            output.write("<DataArray NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n")
            np.savetxt(output, v)
            output.write("</DataArray>\n")
            output.write("</Points>\n")
            output.write("<Cells>\n")
            output.write("<DataArray  type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n")
            np.savetxt(output, elements, fmt="%i")
            output.write("</DataArray>\n")
            output.write("<DataArray  type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n")
            np.savetxt(output, (np.array(range(elements.shape[0]), np.int32)+1)*(dim+1), fmt="%i", newline=" ")
            output.write("</DataArray>\n")
            output.write("<DataArray  type=\"Int32\" Name=\"types\" format=\"ascii\">\n")
            t = np.ndarray((elements.shape[0]), np.int32)
            t[:] = 5 if dim == 2 else 10
            np.savetxt(output, t, fmt="%i", newline=" ")
            output.write("</DataArray>\n")
            output.write("</Cells>\n")
            output.write("<PointData Scalars=\"Pressure Porosity\" Vectors=\"Velocity\">\n")
            output.write("<DataArray Name=\"Porosity\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
            np.savetxt(output, self._pf.porosity)
            output.write("</DataArray>\n")
            output.write("<DataArray Name=\"Pressure\" NumberOfComponents = \"1\" type=\"Float64\" format=\"ascii\">")
            p = np.ndarray((vertices.shape[0]))
            p[np.transpose(elements)] = self._dof.getField(dim, self._sol)
            np.savetxt(output, p)
            output.write("</DataArray>\n")
            output.write("<DataArray Name=\"Velocity\" NumberOfComponents = \"3\" type=\"Float64\" format=\"ascii\">")
            v = np.zeros((vertices.shape[0], 3))
            for d in range (self._jac.dim) :
                v[np.transpose(elements), d] = self._dof.getField(d, self._sol)[:self._jac.dim+1,:]
            np.savetxt(output, v)
            output.write("</DataArray>\n")
            output.write("</PointData>\n")
            output.write("</Piece>\n")
            output.write("</UnstructuredGrid>\n")
            output.write("</VTKFile>\n")
        else :
            print("unknown export format : \"%s\" (choose \"vtk\" or \"msh\")." % filetype)

    def mesh_position(self) :
        return self._jac.x
Vincent Legat's avatar
Vincent Legat committed
173
174
175
176
        
    def solution(self) :
        return self._sol

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
177
178
179
    def set_mesh_position(self, x) :
        self._jac.update(x)

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
180
    def set_particles(self,mass, volume, position, velocity) :
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
181
        self.forces = np.zeros_like(velocity)
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
182
183
        self._pf.set_particles(position, volume, velocity, mass)
        print("%i particles inside the domain" % (self._pf.eid.shape[0]))