bubble.py

# 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
# import sys
# sys.path.append('/home/Documents/migflow/migflow/python')

# import python
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
import xmesh2d
import gmsh

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)

        subprocess.call(["gmsh", "-2", "mesh2.geo","-clscale", "2"])
        t = 0
        ii = 0

        #physical parameters
        g =  np.array([0,-9.81])                                      # gravity
        rho0 = 1                                     # fluid density
        rho1 = 1000
        nu0 = 1e-3                                   # kinematic viscosity
        nu1 = 1e-3                                   # kinematic viscosity
        tEnd = 100000                                     # final time
        Uin = 1
        def boundaryV(x):
            if t%1<0.4:
                return Uin
            else:
                return 0
        
        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)
        

        #numerical parameters
        lcmin = .1                                  # mesh size
        dt = 0.001                                    # time step
        alpha = 1e-4                                    # stabilization coefficient

        shutil.copy("mesh2.msh", outputdir +"/mesh2.msh")

        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])
        fluid.load_msh("mesh2.msh")

        #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])
        fluid.set_wall_boundary("BottomLeft",velocity=[0,0])
        fluid.set_wall_boundary("BottomRight",velocity=[0,0])
        
        #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)
        
        # fluid.set_open_boundary("Insert", velocity=[0, boundaryV], concentration=1)
        fluid.set_open_boundary("Top",pressure=0, concentration=0)
        fluid.set_wall_boundary("Right",velocity=[0,0])
        fluid.set_wall_boundary("Left",velocity=[0,0])
        
        setInitialConcentration(0.04)
        
        
        ii = 0
        t = 0

        #set initial_condition

        fluid.write_vtk(outputdir,0,0)

        mesh,x = xmesh2d.Mesh.import_from_file("mesh2.msh")
        tracker = xmesh2d.FrontTracker(mesh,x)

        ii = 0
        tic = time.time()
        while ii <  200:
            #Fluid solver
            time_integration.iterate(fluid,None,dt)
            t += dt

            concentration = fluid.concentration_cg().reshape((fluid.n_nodes(),))
            print(concentration)
            level = concentration - 0.5;
            #level = np.where(concentration < 0.5, -1, 1)

            newx = np.copy(fluid.coordinates()[:])
            newx = tracker.move_front(newx, level[:,])
            newx = tracker.relax(newx, 1)
            print(x[:,:2].shape)
            print(np.sum(newx-fluid.coordinates()))
            fluid.mesh_velocity()[:] = (newx[:,:2] - fluid.coordinates()[:,:2])/dt
            fluid.coordinates()[:] = newx[:]

            #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()
        # 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)

        # print('Error', (vel.sum())**.5)

        # self.assertLess(vel.sum()**.5,(vS**0.5)/50, "error is too large in Poiseuille")