tcfd_tests_darcy.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_darcy.py (7587B)
---
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.defineWorldBoundaries([1.0,1.0,1.0], dx=0.1)
18 orig.defineWorldBoundaries([0.4,0.3,0.4], dx=0.1)
19 orig.initFluid(cfd_solver = 1)
20 #orig.initFluid(mu = 8.9e-4)
21 orig.initTemporal(total = 0.2, file_dt = 0.01, dt = 1.0e-7)
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 py = sphere.sim(sid = orig.sid, fluid = True)
27 orig.run(verbose=False)
28 #orig.run(verbose=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/3)
35 it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
36 compare(it[:,1].sum(), 0.0, "Convergence rate (1/3):\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 print("# Pressure gradient")
51 orig.cleanup()
52 orig.p_f[:,:,-1] = 1.0
53 orig.initTemporal(total = 0.5, file_dt = 0.01, dt = 1.0e-6)
54 #orig.time_dt[0] *= 0.01
55 #orig.time_file_dt = orig.time_dt*0.99
56 #orig.time_total = orig.time_dt*1
57 #orig.run(device=2, verbose=False)
58 orig.run(verbose=False)
59 py.readlast(verbose = False)
60 ideal_grad_p_z = numpy.linspace(orig.p_f[0,0,0], orig.p_f[0,0,-1], orig.num[2])
61 #orig.writeVTKall()
62 compareNumpyArraysClose(numpy.zeros((1,orig.num[2])),\
63 ideal_grad_p_z - py.p_f[0,0,:],\
64 "Pressure gradient:\t", tolerance=1.0e-1)
65 #"Pressure gradient:\t", tolerance=1.0e-2)
66
67 # Fluid flow direction, opposite of gradient (i.e. towards -z)
68 if ((py.v_f[:,:,:,2] < 0.0).all() and (py.v_f[:,:,:,0:1] < 1.0e-7).all()):
69 print("Flow field:\t\t" + passed())
70 else:
71 print("Flow field:\t\t" + failed())
72 raise Exception("Failed")
73
74 # Convergence rate (2/3)
75 it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
76 if ((it[0:6,1] < 1000).all() and (it[6:,1] < 20).all()):
77 print("Convergence rate (2/3):\t" + passed())
78 else:
79 print("Convergence rate (2/3):\t" + failed())
80 #'''
81
82 # Long test
83 '''
84 #orig.p_f[:,:,-1] = 1.1
85 orig.time_total[0] = 0.1
86 orig.time_file_dt[0] = orig.time_total[0]/10.0
87 orig.run(verbose=False)
88 py.readlast(verbose = False)
89 ideal_grad_p_z = numpy.linspace(orig.p_f[0,0,0], orig.p_f[0,0,-1], orig.num[2])
90 #py.writeVTKall()
91 compareNumpyArraysClose(numpy.zeros((1,orig.num[2])),\
92 ideal_grad_p_z - py.p_f[0,0,:],\
93 "Pressure gradient (long test):", tolerance=1.0e-2)
94
95 # Fluid flow direction, opposite of gradient (i.e. towards -z)
96 if ((py.v_f[:,:,:,2] < 0.0).all() and (py.v_f[:,:,:,0:1] < 1.0e-7).all()):
97 print("Flow field:\t\t" + passed())
98 else:
99 print("Flow field:\t\t" + failed())
100
101 # Convergence rate (3/3)
102 # This test passes with BETA=0.0 and tolerance=1.0e-9
103 it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
104 if (it[0,1] < 700 and it[1,1] < 250 and (it[2:,1] < 20).all()):
105 print("Convergence rate (3/3):\t" + passed())
106 else:
107 print("Convergence rate (3/3):\t" + failed())
108 '''
109
110 '''
111 # Slow pressure modulation test
112 orig.cleanup()
113 orig.time_total[0] = 1.0e-1
114 orig.time_file_dt[0] = 0.101*orig.time_total[0]
115 orig.setFluidPressureModulation(A=1.0, f=1.0/orig.time_total[0])
116 #orig.plotPrescribedFluidPressures()
117 orig.run(verbose=False)
118 #py.readlast()
119 #py.plotConvergence()
120 #py.plotFluidDiffAdvPresZ()
121 #py.writeVTKall()
122 for it in range(1,py.status()): # gradient should be smooth in all output files
123 py.readstep(it, verbose=False)
124 ideal_grad_p_z =\
125 numpy.linspace(py.p_f[0,0,0], py.p_f[0,0,-1], py.num[2])
126 compareNumpyArraysClose(numpy.zeros((1,py.num[2])),\
127 ideal_grad_p_z - py.p_f[0,0,:],\
128 'Slow pressure modulation (' +
129 str(it+1) + '/' + str(py.status()) + '):', tolerance=1.0e-1)
130 '''
131
132 #'''
133 print("# Fast pressure modulation")
134 orig.time_total[0] = 1.0e-2
135 orig.time_file_dt[0] = 0.101*orig.time_total[0]
136 orig.setFluidPressureModulation(A=1.0, f=1.0/orig.time_total[0])
137 #orig.plotPrescribedFluidPressures()
138 orig.run(verbose=False)
139 #py.plotConvergence()
140 #py.plotFluidDiffAdvPresZ()
141 #py.writeVTKall()
142 for it in range(1,py.status()+1): # gradient should be smooth in all output files
143 py.readstep(it, verbose=False)
144 #py.plotFluidDiffAdvPresZ()
145 ideal_grad_p_z =\
146 numpy.linspace(py.p_f[0,0,0], py.p_f[0,0,-1], py.num[2])
147 compareNumpyArraysClose(numpy.zeros((1,py.num[2])),\
148 ideal_grad_p_z - py.p_f[0,0,:],\
149 'Fast pressure modulation (' +
150 str(it) + '/' + str(py.status()) + '):', tolerance=5.0e-1)
151 #'''
152
153 print("# Pressure perturbation")
154 orig = sphere.sim(np = 0, nd = 3, nw = 0, sid = "cfdtest", fluid = True)
155 cleanup(orig)
156 orig.defaultParams()
157 orig.defineWorldBoundaries([1.0,1.0,1.0], dx=0.1)
158 #orig.defineWorldBoundaries([0.4,0.3,0.4], dx=0.1)
159 orig.initFluid(cfd_solver = 1)
160 #orig.initFluid(mu = 8.9e-4)
161 orig.initTemporal(total = 0.2, file_dt = 0.01, dt = 1.0e-7)
162 #orig.g[2] = -10.0
163 orig.time_file_dt = orig.time_dt*0.99
164 orig.time_total = orig.time_dt*10
165 #orig.run(dry=True)
166 orig.p_f[4,2,5] = 2.0
167 #orig.run(verbose=False)
168 orig.run(verbose=True)
169 py = sphere.sim(sid = orig.sid, fluid = True)
170 #py.writeVTKall()
171
172
173 #ones = numpy.ones((orig.num))
174 #py.readlast(verbose = False)
175 #compareNumpyArrays(ones, py.p_f, "Conservation of pressure:")
176
177 # Convergence rate (1/3)
178 #it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
179 #compare(it[:,1].sum(), 0.0, "Convergence rate (1/3):\t")
180
181 # Fluid flow should be very small
182 #if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
183 # print("Flow field:\t\t" + passed())
184 #else:
185 # print("Flow field:\t\t" + failed())
186 # print(numpy.min(py.v_f))
187 # print(numpy.mean(py.v_f))
188 # print(numpy.max(py.v_f))
189 # raise Exception("Failed")
190
191
192
193 print('## Flux BC tests')
194 print('# Flux top BC test')
195 orig = sphere.sim(np = 0, nd = 3, nw = 0, sid = "cfdtest", fluid = True)
196 cleanup(orig)
197 orig.defaultParams()
198 orig.defineWorldBoundaries([1.0,1.0,1.0], dx=0.1)
199 #orig.defineWorldBoundaries([0.4,0.3,0.4], dx=0.1)
200 orig.initFluid(cfd_solver = 1)
201 #orig.initFluid(mu = 8.9e-4)
202 orig.initTemporal(total = 0.2, file_dt = 0.01, dt = 1.0e-7)
203 #orig.g[2] = -10.0
204 orig.time_file_dt = orig.time_dt*0.99
205 orig.time_total = orig.time_dt*10
206 #orig.run(dry=True)
207 orig.setFluidTopFixedFlux(1.0)
208 #orig.run(verbose=False)
209 orig.run(verbose=True)
210 py = sphere.sim(sid = orig.sid, fluid = True)
211 #py.writeVTKall()
212
213
214
215 # Add horizontal pressure gradient along X
216 print("# Pressure gradient along X")
217 orig.cleanup()
218 orig.p_f[:,:,:] = 0.0
219 orig.p_f[0,:,:] = 2.0
220 orig.setFluidXFixedPressure()
221 orig.setFluidYNoFlow()
222 #orig.setFluidTopFixedPressure()
223 orig.setFluidTopNoFlow()
224 orig.setFluidBottomNoFlow()
225 orig.initTemporal(total = 0.5, file_dt = 0.01, dt = 1.0e-6)
226 #orig.time_dt[0] *= 0.01
227 #orig.time_file_dt = orig.time_dt*0.99
228 #orig.time_total = orig.time_dt*1
229 #orig.run(device=2, verbose=False)
230 orig.run(verbose=False)
231 py.readlast(verbose = False)
232
233 # Fluid flow direction, opposite of gradient (i.e. towards +x)
234 if ((py.v_f[:,:,:,0] > 0.0).all()):
235 print("Flow field (X):\t\t" + passed())
236 else:
237 print("Flow field (X):\t\t" + failed())
238 raise Exception("Failed")
239
240
241
242 #cleanup(orig)