cfd_tests.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
cfd_tests.py (7678B)
---
1 #!/usr/bin/env python
2 from pytestutils import *
3
4 import sphere
5 import sys
6 import numpy
7
8 print("### CFD tests - Dirichlet BCs ###")
9
10 # Iteration and conservation of mass test
11 # No gravity, no pressure gradients => no flow
12 print('# No forcing')
13 orig = sphere.sim(np = 0, nd = 3, nw = 0, sid = "cfdtest", fluid = True)
14 cleanup(orig)
15 orig.defaultParams()
16 orig.addParticle([0.5,0.5,0.5], 0.05)
17 orig.defineWorldBoundaries([1.0,1.0,1.0])
18 orig.initFluid(mu = 0.0)
19 #orig.initFluid(mu = 8.9e-4)
20 orig.initTemporal(total = 0.2, file_dt = 0.01)
21 #orig.g[2] = -10.0
22 orig.time_file_dt = orig.time_dt*0.99
23 orig.time_total = orig.time_dt*10
24 #orig.run(dry=True)
25 orig.run(verbose=False)
26 #orig.run(verbose=True)
27 py = sphere.sim(sid = orig.sid, fluid = True)
28
29 zeros = numpy.zeros((orig.num))
30 py.readlast(verbose = False)
31 compareNumpyArrays(zeros, py.p_f, "Conservation of pressure:")
32
33 # Convergence rate (1/2)
34 it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
35 compare(it[:,1].sum(), 0.0, "Convergence rate (1/2):\t")
36
37 # Fluid flow should be very small
38 if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
39 print("Flow field:\t\t" + passed())
40 else:
41 print("Flow field:\t\t" + failed())
42 print(numpy.min(py.v_f))
43 print(numpy.mean(py.v_f))
44 print(numpy.max(py.v_f))
45 raise Exception("Failed")
46
47
48 # Add pressure gradient
49 # This test passes with BETA=0.0 and tolerance=1.0e-9
50 print('# Pressure gradient')
51 orig.p_f[:,:,-1] = 1.1
52 orig.run(verbose=False)
53 #orig.run(verbose=True)
54 py.readlast(verbose = False)
55 ideal_grad_p_z = numpy.linspace(orig.p_f[0,0,0], orig.p_f[0,0,-1], orig.num[2])
56 #orig.writeVTKall()
57 compareNumpyArraysClose(numpy.zeros((1,orig.num[2])),\
58 ideal_grad_p_z - py.p_f[0,0,:],\
59 "Pressure gradient:\t", tolerance=1.0e-1)
60 #"Pressure gradient:\t", tolerance=1.0e-2)
61
62 # Fluid flow direction, opposite of gradient (i.e. towards -z)
63 if ((py.v_f[:,:,:,2] < 0.0).all() and (py.v_f[:,:,:,0:1] < 1.0e-7).all()):
64 print("Flow field:\t\t" + passed())
65 else:
66 print("Flow field:\t\t" + failed())
67 raise Exception("Failed")
68
69 # Convergence rate (2/2)
70 # This test passes with BETA=0.0 and tolerance=1.0e-9
71 it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
72 if ((it[0:6,1] < 1000).all() and (it[6:,1] < 20).all()):
73 print("Convergence rate (2/2):\t" + passed())
74 else:
75 print("Convergence rate (2/2):\t" + failed())
76
77 '''
78 # Long test
79 # This test passes with BETA=0.0 and tolerance=1.0e-9
80 orig.p_f[:,:,-1] = 1.1
81 orig.time_total[0] = 5.0
82 orig.time_file_dt[0] = orig.time_total[0]/10.0
83 orig.run(verbose=True)
84 py.readlast(verbose = False)
85 ideal_grad_p_z = numpy.linspace(orig.p_f[0,0,0], orig.p_f[0,0,-1], orig.num[2])
86 py.writeVTKall()
87 compareNumpyArraysClose(numpy.zeros((1,orig.num[2])),\
88 ideal_grad_p_z - py.p_f[0,0,:],\
89 "Pressure gradient (long test):", tolerance=1.0e-2)
90
91 # Fluid flow direction, opposite of gradient (i.e. towards -z)
92 if ((py.v_f[:,:,:,2] < 0.0).all() and (py.v_f[:,:,:,0:1] < 1.0e-7).all()):
93 print("Flow field:\t\t" + passed())
94 else:
95 print("Flow field:\t\t" + failed())
96
97 # Convergence rate (2/2)
98 # This test passes with BETA=0.0 and tolerance=1.0e-9
99 it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
100 if (it[0,1] < 700 and it[1,1] < 250 and (it[2:,1] < 20).all()):
101 print("Convergence rate (2/2):\t" + passed())
102 else:
103 print("Convergence rate (2/2):\t" + failed())
104 '''
105 # Add viscosity which will limit the fluid flow. Used to test the stress tensor
106 # in the fluid velocity prediction
107 #print(numpy.mean(py.v_f[:,:,:,2]))
108 print('# Viscid flow')
109 orig.time_file_dt[0] = 1.0e-4
110 orig.time_total[0] = 1.0e-3
111 orig.initFluid(mu = 8.9-4) # water at 25 deg C
112 orig.p_f[:,:,-1] = 2.0
113 orig.run(verbose=False)
114 #orig.writeVTKall()
115
116 #py.plotConvergence()
117
118 py.readsecond(verbose=False)
119 #py.plotFluidDiffAdvPresZ()
120
121 # The v_z values are read from sb.v_f[0,0,:,2]
122 dz = py.L[2]/py.num[2]
123 rho = 1000.0 # fluid density
124
125 # Central difference gradients
126 dvz_dz = (py.v_f[0,0,1:,2] - py.v_f[0,0,:-1,2])/(2.0*dz)
127 dvzvz_dz = (py.v_f[0,0,1:,2]**2 - py.v_f[0,0,:-1,2]**2)/(2.0*dz)
128
129 # Diffusive contribution to velocity change
130 dvz_diff = 2.0*py.mu/rho*dvz_dz*py.time_dt
131
132 # Advective contribution to velocity change
133 dvz_adv = dvzvz_dz*py.time_dt
134
135 # Diffusive and advective terms should have opposite terms
136 if ((numpy.sign(dvz_diff) == numpy.sign(-dvz_adv)).all()):
137 print("Diffusion-advection (1/2):" + passed())
138 else:
139 print("Diffusion-advection (1/2):" + failed())
140 raise Exception("Failed")
141
142
143 py.readlast(verbose=False)
144 #py.plotFluidDiffAdvPresZ()
145
146 # The v_z values are read from sb.v_f[0,0,:,2]
147 dz = py.L[2]/py.num[2]
148 rho = 1000.0 # fluid density
149
150 # Central difference gradients
151 dvz_dz = (py.v_f[0,0,1:,2] - py.v_f[0,0,:-1,2])/(2.0*dz)
152 dvzvz_dz = (py.v_f[0,0,1:,2]**2 - py.v_f[0,0,:-1,2]**2)/(2.0*dz)
153
154 # Diffusive contribution to velocity change
155 dvz_diff = 2.0*py.mu/rho*dvz_dz*py.time_dt
156
157 # Advective contribution to velocity change
158 dvz_adv = dvzvz_dz*py.time_dt
159
160 # Diffusive and advective terms should have opposite terms
161 if ((numpy.sign(dvz_diff) == numpy.sign(-dvz_adv)).all()):
162 print("Diffusion-advection (2/2):" + passed())
163 else:
164 print("Diffusion-advection (2/2):" + failed())
165
166
167 # Slow pressure modulation test
168 '''
169 orig.time_total[0] = 1.0e-1
170 orig.time_file_dt[0] = 0.101*orig.time_total[0]
171 orig.mu[0] = 0.0 # dont let diffusion add transient effects
172 orig.setFluidPressureModulation(A=1.0, f=1.0/orig.time_total[0])
173 #orig.plotPrescribedFluidPressures()
174 orig.run(verbose=False)
175 #py.readlast()
176 #py.plotConvergence()
177 #py.plotFluidDiffAdvPresZ()
178 #py.writeVTKall()
179 for it in range(1,py.status()): # gradient should be smooth in all output files
180 py.readstep(it)
181 ideal_grad_p_z =\
182 numpy.linspace(py.p_f[0,0,0], py.p_f[0,0,-1], py.num[2])
183 compareNumpyArraysClose(numpy.zeros((1,py.num[2])),\
184 ideal_grad_p_z - py.p_f[0,0,:],\
185 'Slow pressure modulation (' +
186 str(it+1) + '/' + str(py.status()) + '):', tolerance=1.0e-1)
187 '''
188
189 # Fast pressure modulation test
190 print('# Fast pressure modulation')
191 orig.time_total[0] = 1.0e-2
192 orig.time_file_dt[0] = 0.101*orig.time_total[0]
193 orig.mu[0] = 0.0 # dont let diffusion add transient effects
194 orig.setFluidPressureModulation(A=1.0, f=1.0/orig.time_total[0])
195 #orig.plotPrescribedFluidPressures()
196 orig.run(verbose=False)
197 #py.plotConvergence()
198 #py.plotFluidDiffAdvPresZ()
199 #py.writeVTKall()
200 for it in range(1,py.status()+1): # gradient should be smooth in all output files
201 py.readstep(it, verbose=False)
202 #py.plotFluidDiffAdvPresZ()
203 ideal_grad_p_z =\
204 numpy.linspace(py.p_f[0,0,0], py.p_f[0,0,-1], py.num[2])
205 compareNumpyArraysClose(numpy.zeros((1,py.num[2])),\
206 ideal_grad_p_z - py.p_f[0,0,:],\
207 'Fast pressure modulation (' +
208 str(it) + '/' + str(py.status()) + '):', tolerance=5.0e-1)
209
210 '''
211 # Top: Dirichlet, bot: Neumann
212 orig.disableFluidPressureModulation()
213 orig.time_total[0] = 1.0e-2
214 orig.time_file_dt = orig.time_total/20
215 orig.p_f[:,:,-1] = 1.0
216 orig.g[2] = -1.0
217 orig.mu[0] = 8.9e-4 # water
218 orig.bc_bot[0] = 1 # No-flow BC at bottom
219 #orig.run(dry=True)
220 orig.run(verbose=False)
221 orig.writeVTKall()
222 py.readlast(verbose = False)
223 #ideal_grad_p_z = numpy.linspace(orig.p_f[0,0,0], orig.p_f[0,0,-1], orig.num[2])
224 #compareNumpyArraysClose(numpy.zeros((1,orig.num[2])),\
225 #ideal_grad_p_z - py.p_f[0,0,:],\
226 #"Pressure gradient:\t", tolerance=1.0e-2)
227
228 # Fluid flow direction, opposite of gradient (i.e. towards -z)
229 #if ((py.v_f[:,:,:,2] < 0.0).all() and (py.v_f[:,:,:,0:1] < 1.0e-7).all()):
230 #print("Flow field:\t\t" + passed())
231 #else:
232 #print("Flow field:\t\t" + failed())
233 #raise Exception("Failed")
234
235 '''
236 cleanup(orig)