1
+ from LBM_python import *
2
+ import matplotlib .pyplot as plt
3
+ #parameters
4
+ nx = 200
5
+ ny = 100
6
+ center_x = 10
7
+ center_y = 5
8
+ radius = 2
9
+ omega1 = 1
10
+ omega2 = 1
11
+ #define geometry
12
+ circle = Indicator .circle (center_x ,center_y ,radius )
13
+ botBox = Indicator .cuboid (80 ,1 ,120 ,40 )
14
+ bottomPlate = Indicator .cuboid (0 ,0 ,nx ,0 )
15
+
16
+ cGeometry = CuboidGeometry2D (0 ,0 ,nx ,ny )
17
+ cGeometry .setPeriodicity ()
18
+ superG = SuperGeometry (cGeometry )
19
+ superG .rename (0 ,1 )
20
+ superG .rename (1 ,3 ,bottomPlate )
21
+ superG .rename (1 ,2 ,botBox )
22
+
23
+ print ('================================================' )
24
+ #lattice
25
+ rho = 1
26
+ rhoZero = 0
27
+ u = [0 ,0. ]
28
+ sLattice1 = SuperLattice2D (superG )
29
+ sLattice2 = SuperLattice2D (superG )
30
+
31
+ sLattice1 .defineRhoU (superG ,1 ,rho ,u )
32
+ sLattice1 .defineRhoU (superG ,2 ,rhoZero ,u )
33
+ #sLattice1.defineRhoU(superG,3,0.5,u)
34
+
35
+
36
+ sLattice2 .defineRhoU (superG ,1 ,rhoZero ,u )
37
+ sLattice2 .defineRhoU (superG ,2 ,rho ,u )
38
+ #sLattice2.defineRhoU(superG,3,0.5,u)
39
+
40
+ bulk1 = ShanChenBGKdynamics (omega1 )
41
+ bulk2 = ShanChenBGKdynamics (omega2 )
42
+
43
+
44
+ bounceback1 = BBwall (0.6 )
45
+ bounceback2 = BBwall (0.4 )
46
+ sLattice1 .defineDynamics (superG ,1 ,bulk1 )# SuperGeometry, materialNum, dynamics
47
+ sLattice1 .defineDynamics (superG ,2 ,bulk1 )
48
+ sLattice1 .defineDynamics (superG ,3 ,bounceback1 )
49
+
50
+ sLattice2 .defineDynamics (superG ,1 ,bulk2 )
51
+ sLattice2 .defineDynamics (superG ,2 ,bulk2 )
52
+ sLattice2 .defineDynamics (superG ,3 ,bounceback2 )
53
+
54
+
55
+ outputDirectory = 'data'
56
+ if not os .path .exists (outputDirectory ):
57
+ os .makedirs (outputDirectory )
58
+
59
+
60
+ # maxVelocity = numpy.array([0.01,0])
61
+ # poV = Poiseuille2D(superG,3,maxVelocity,0.5).getVelocityField() #SuperGeometry,materialNum,maxVelocity,distance2Wall
62
+ # poV = numpy.array(maxVelocity)
63
+
64
+ G = 3
65
+
66
+ fig = plt .figure ()
67
+ # MAIN LOOP
68
+
69
+
70
+ sLattice1 .addLatticeCoupling (superG ,G ,sLattice2 )
71
+
72
+ print ('initial average rho1: {0:.5f}' .format (sLattice1 .getAverageRho ())) #initialize
73
+ print ('initial average rho2: {0:.5f}' .format (sLattice2 .getAverageRho ()))
74
+
75
+ #ok now
76
+
77
+ for iT in numpy .arange ( 2000 ):
78
+
79
+ # sLattice.openBC(superG,3)
80
+ # sLattice1.updateRhoU() #13 #needed for defineU_BC / defineRho_BC
81
+ # sLattice2.updateRhoU() #13
82
+
83
+ # sLattice.defineU_BC(superG,4,poV)
84
+ # sLattice.defineRho_BC(superG,3,1)
85
+
86
+ if iT % 50 == 0 :
87
+ im = plt .imshow (sLattice1 .rhoMap , animated = True )
88
+ plt .draw ()
89
+ if iT == 0 :
90
+ plt .pause (5 )
91
+ else :
92
+ plt .pause (0.00001 )
93
+ # numpy.savetxt('{}/rho1_{}'.format(outputDirectory,iT),sLattice1.getRhoMap())
94
+ print ('{}/2000' .format (iT ))
95
+
96
+
97
+
98
+ sLattice1 .prepareCoupling ()
99
+ sLattice2 .prepareCoupling ()
100
+
101
+ sLattice1 .executeCoupling (sLattice2 )
102
+
103
+ sLattice1 .collide ()
104
+ sLattice2 .collide ()
105
+ sLattice1 .stream ()
106
+ sLattice2 .stream ()
107
+
108
+ # print('===============================final Ux:\n{}\n'.format(sLattice1.getUxMap()))
109
+
110
+ print ('final average rho1: {0:.5f}' .format (sLattice1 .getAverageRho ()))
111
+ print ('final average rho2: {0:.5f}' .format (sLattice2 .getAverageRho ()))
112
+
113
+ # for i in numpy.arange(sLattice1.rhoMap.shape[0]):
114
+ # for j in numpy.arange(sLattice1.rhoMap.shape[1]):
115
+ # print('%10.4f' %sLattice1.rhoMap[i,j],end='')
116
+ # print('')
117
+
118
+ print ('=========================================================================================' )
119
+ # for i in numpy.arange(sLattice1.rhoMap.shape[0]):
120
+ # for j in numpy.arange(sLattice1.rhoMap.shape[1]):
121
+ # print('%10.4f' %sLattice2.rhoMap[i,j],end='')
122
+ # print('')
0 commit comments