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