bubble.py 5.31 KB
Newer Older
dazerki's avatar
dazerki committed
1
2
3
4
5
6
7
8
9
10
# 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/>.

#!/usr/bin/env python
dazerki's avatar
dazerki committed
11
12
13
14
# import sys
# sys.path.append('/home/Documents/migflow/migflow/python')

# import python
dazerki's avatar
dazerki committed
15
16
17
18
19
20
21
22
23
24
from migflow import fluid as mbfluid
from migflow import time_integration

import numpy as np
import os
import subprocess
import time
import shutil
import random
import unittest
dazerki's avatar
dazerki committed
25
import xmesh2d
dazerki's avatar
dazerki committed
26
import gmsh
dazerki's avatar
dazerki committed
27
28
29
30
31
32
33
34
35

class Poiseuille(unittest.TestCase) :
    def runTest(self) :
        dir_path = os.path.dirname(os.path.realpath(__file__))
        os.chdir(dir_path)
        outputdir = "output"
        if not os.path.isdir(outputdir) :
            os.makedirs(outputdir)

tleyssens6's avatar
tleyssens6 committed
36
        subprocess.call(["gmsh", "-2", "mesh2.geo","-clscale", "2"])
dazerki's avatar
dazerki committed
37
38
39
40
        t = 0
        ii = 0

        #physical parameters
tleyssens6's avatar
BCs    
tleyssens6 committed
41
        g =  np.array([0,-9.81])                                      # gravity
dazerki's avatar
dazerki committed
42
        rho0 = 1                                     # fluid density
dazerki's avatar
dazerki committed
43
        rho1 = 1000
dazerki's avatar
dazerki committed
44
45
46
        nu0 = 1e-3                                   # kinematic viscosity
        nu1 = 1e-3                                   # kinematic viscosity
        tEnd = 100000                                     # final time
dazerki's avatar
dazerki committed
47
48
49
50
51
52
        Uin = 1
        def boundaryV(x):
            if t%1<0.4:
                return Uin
            else:
                return 0
tleyssens6's avatar
tleyssens6 committed
53
54
55
56
57
58
59
60
61
62
63
64
65
        
        def setInitialConcentration(r):
            concentration = np.full(fluid.n_nodes(),0)
            #tags, x, _ = gmsh.model.mesh.get_nodes(includeBoundary=False)
            x = fluid.coordinates()
            print(x.shape)
            #x = x.reshape([-1,3])
            # circle centered at (0,0.15), radius r
            for i in range(fluid.n_nodes()):
                if (x[i,0]-0)**2 + (x[i,1]-0.15)**2 < r**2 :
                    concentration[i] = 1
            fluid.set_concentration_cg(concentration)
        
dazerki's avatar
dazerki committed
66
67
68

        #numerical parameters
        lcmin = .1                                  # mesh size
tleyssens6's avatar
tleyssens6 committed
69
        dt = 0.001                                    # time step
dazerki's avatar
dazerki committed
70
71
        alpha = 1e-4                                    # stabilization coefficient

tleyssens6's avatar
tleyssens6 committed
72
        shutil.copy("mesh2.msh", outputdir +"/mesh2.msh")
dazerki's avatar
dazerki committed
73
74
75
76
77

        outf = 1                                       # number of iterations between output files

        #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
        fluid = mbfluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1])
tleyssens6's avatar
tleyssens6 committed
78
        fluid.load_msh("mesh2.msh")
dazerki's avatar
dazerki committed
79

tleyssens6's avatar
tleyssens6 committed
80
81
82
83
84
        #fluid.set_open_boundary("BottomMiddle", velocity=[0, boundaryV], concentration=1)
        #fluid.set_wall_boundary("BottomLeft",velocity=[0,0])
        #fluid.set_wall_boundary("BottomRight",velocity=[0,0])
        
        fluid.set_wall_boundary("BottomMiddle", velocity=[0,0])
dazerki's avatar
dazerki committed
85
86
        fluid.set_wall_boundary("BottomLeft",velocity=[0,0])
        fluid.set_wall_boundary("BottomRight",velocity=[0,0])
tleyssens6's avatar
tleyssens6 committed
87
88
89
90
91
        
        #fluid.set_open_boundary("BottomMiddle", velocity=[0, boundaryV], concentration=1)
        #fluid.set_open_boundary("BottomLeft",velocity=[0, boundaryV], concentration=1)
        #fluid.set_open_boundary("BottomRight",velocity=[0, boundaryV], concentration=1)
        
dazerki's avatar
dazerki committed
92
        # fluid.set_open_boundary("Insert", velocity=[0, boundaryV], concentration=1)
dazerki's avatar
dazerki committed
93
        fluid.set_open_boundary("Top",pressure=0, concentration=0)
tleyssens6's avatar
BCs    
tleyssens6 committed
94
95
        fluid.set_wall_boundary("Right",velocity=[0,0])
        fluid.set_wall_boundary("Left",velocity=[0,0])
tleyssens6's avatar
tleyssens6 committed
96
97
98
99
        
        setInitialConcentration(0.04)
        
        
dazerki's avatar
dazerki committed
100
101
102
103
104
105
106
        ii = 0
        t = 0

        #set initial_condition

        fluid.write_vtk(outputdir,0,0)

tleyssens6's avatar
tleyssens6 committed
107
        mesh,x = xmesh2d.Mesh.import_from_file("mesh2.msh")
dazerki's avatar
dazerki committed
108
109
        tracker = xmesh2d.FrontTracker(mesh,x)

dazerki's avatar
dazerki committed
110
111
        ii = 0
        tic = time.time()
tleyssens6's avatar
tleyssens6 committed
112
        while ii <  200:
dazerki's avatar
dazerki committed
113
114
115
            #Fluid solver
            time_integration.iterate(fluid,None,dt)
            t += dt
dazerki's avatar
dazerki committed
116

dazerki's avatar
dazerki committed
117
            concentration = fluid.concentration_cg().reshape((fluid.n_nodes(),))
tleyssens6's avatar
tleyssens6 committed
118
            print(concentration)
dazerki's avatar
dazerki committed
119
            level = concentration - 0.5;
tleyssens6's avatar
tleyssens6 committed
120
            #level = np.where(concentration < 0.5, -1, 1)
dazerki's avatar
dazerki committed
121
122
123

            newx = np.copy(fluid.coordinates()[:])
            newx = tracker.move_front(newx, level[:,])
tleyssens6's avatar
tleyssens6 committed
124
            newx = tracker.relax(newx, 1)
dazerki's avatar
dazerki committed
125
126
127
128
            print(x[:,:2].shape)
            print(np.sum(newx-fluid.coordinates()))
            fluid.mesh_velocity()[:] = (newx[:,:2] - fluid.coordinates()[:,:2])/dt
            fluid.coordinates()[:] = newx[:]
dazerki's avatar
dazerki committed
129

dazerki's avatar
dazerki committed
130
131
132
133
134
135
136
137
138
            #Output files writting
            if ii %outf == 0 :
                ioutput = int(ii/outf) + 1
                fluid.write_vtk(outputdir, ioutput, t)
            ii += 1
            print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))

        s = fluid.solution()
        x = fluid.coordinates()
dazerki's avatar
dazerki committed
139
140
        # vel = (s[:,0]-1/(20*nu0*rho0)*x[:,1]*(1-x[:, 1]))**2
        # vS = np.sum((1/(20*nu0*rho0)*x[:,1]*(1-x[:, 1]))**2)
dazerki's avatar
dazerki committed
141

dazerki's avatar
dazerki committed
142
        # print('Error', (vel.sum())**.5)
dazerki's avatar
dazerki committed
143

dazerki's avatar
dazerki committed
144
        # self.assertLess(vel.sum()**.5,(vS**0.5)/50, "error is too large in Poiseuille")