Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
compute_coag.py
1import sys
2
3from .solver_DG import *
4
5
6def compute_coag_k0(eps, coeff_CFL, nbins, massgrid, gij, tensor_tabflux_coag, dthydro):
7 """
8 Function to compute coagulation solver for 1 hydro time-step
9
10 DG scheme k=0, piecewise constant approximation
11
12 Parameters
13 ----------
14 eps : scalar, type -> float
15 minimum value for mass distribution approximation gij
16 coeff_CFL : scalar, type -> float
17 timestep coefficient for stability of the SSPRK order 3 scheme
18 nbins : scalar, type -> integer
19 number of dust bins
20 massgrid : 1D array (dim = nbins+1), type -> float
21 grid of masses given borders value of mass bins
22 gij : 1D array (dim = nbins), type -> float
23 components of g on the polynomial basis
24 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
25 array to evaluate coagulation flux
26 dthydro : scalar, type -> float
27 hydro timestep, used as timestep to reach for coagulation process
28
29
30 Returns
31 -------
32 gij : 1D array (dim = nbins), type -> float
33 evolved components of g on the polynomial basis after 1 hydro timestep
34 nsub : scalar, type -> integer
35 number of subcycling coagulation timestep to reach dthydro
36 ndt : scalar, type -> integer
37 number of hydro timestep, when coagulation CFL > dthydro
38
39 """
40
41 nsub = 0
42 ndt = 0
43
44 # evaluate coagulation CFL
45 dtCFLsub = compute_CFL_coag_k0(eps, nbins, massgrid, gij, tensor_tabflux_coag)
46 dtCFLsub = coeff_CFL * dtCFLsub
47 # print("dtCFLsub=",dtCFLsub)
48
49 # compare hydro timestep and coagulation CFL
50 dt = min(dtCFLsub, dthydro)
51
52 # coagulation subcycling timesteps
53 if dt < dthydro:
54 dtsub = 0.0
55 while dtsub < dthydro and dthydro - dtsub > dtCFLsub:
56 dtsub += dtCFLsub
57
58 nsub += 1
59
60 gij = solver_coag_k0(eps, nbins, massgrid, gij, tensor_tabflux_coag, dtCFLsub)
61
62 dtCFLsub = compute_CFL_coag_k0(eps, nbins, massgrid, gij, tensor_tabflux_coag)
63 dtCFLsub = coeff_CFL * dtCFLsub
64
65 # print("gij=",gij)
66 # sys.exit(-1)
67
68 # last timestep to reach dthydro
69 dtlast = dthydro - dtsub
70 nsub += 1
71
72 gij = solver_coag_k0(eps, nbins, massgrid, gij, tensor_tabflux_coag, dtlast)
73
74 # when coagulation CFL > hydro timstep
75 else:
76 ndt += 1
77
78 gij = solver_coag_k0(eps, nbins, massgrid, gij, tensor_tabflux_coag, dt)
79
80 return gij, nsub, ndt
81
82
83def compute_coag_k0_kdv(eps, coeff_CFL, nbins, massgrid, gij, tensor_tabflux_coag, dv, dthydro):
84 """
85 Function to compute coagulation solver for 1 hydro time-step
86
87 DG scheme k=0, piecewise constant approximation
88
89 Function for ballistic kernel with differential velocities dv
90
91 Parameters
92 ----------
93 eps : scalar, type -> float
94 minimum value for mass distribution approximation gij
95 coeff_CFL : scalar, type -> float
96 timestep coefficient for stability of the SSPRK order 3 scheme
97 nbins : scalar, type -> integer
98 number of dust bins
99 massgrid : 1D array (dim = nbins+1), type -> float
100 grid of masses given borders value of mass bins
101 gij : 1D array (dim = nbins), type -> float
102 components of g on the polynomial basis
103 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
104 array to evaluate coagulation flux
105 dv : 2D array (dim = (nbins,nbins)), type -> float
106 array of the differential velocity between grains
107 dthydro : scalar, type -> float
108 hydro timestep, used as timestep to reach for coagulation process
109
110
111 Returns
112 -------
113 gij : 1D array (dim = nbins), type -> float
114 evolved components of g on the polynomial basis after 1 hydro timestep
115 nsub : scalar, type -> integer
116 number of subcycling coagulation timestep to reach dthydro
117 ndt : scalar, type -> integer
118 number of hydro timestep, when coagulation CFL > dthydro
119
120 """
121
122 nsub = 0
123 ndt = 0
124
125 # evaluate coagulation CFL
126 dtCFLsub = compute_CFL_coag_k0_kdv(eps, nbins, massgrid, gij, tensor_tabflux_coag, dv)
127 dtCFLsub = coeff_CFL * dtCFLsub
128 # print("dtCFLsub=",dtCFLsub)
129
130 # compare hydro timestep and coagulation CFL
131 dt = min(dtCFLsub, dthydro)
132
133 # coagulation subcycling timesteps
134 if dt < dthydro:
135 dtsub = 0.0
136 while dtsub < dthydro and dthydro - dtsub > dtCFLsub:
137 dtsub += dtCFLsub
138
139 nsub += 1
140
141 gij = solver_coag_k0_kdv(eps, nbins, massgrid, gij, tensor_tabflux_coag, dv, dtCFLsub)
142
143 dtCFLsub = compute_CFL_coag_k0_kdv(eps, nbins, massgrid, gij, tensor_tabflux_coag, dv)
144 dtCFLsub = coeff_CFL * dtCFLsub
145
146 # print("dtCFLsub=",dtCFLsub)
147
148 # last timestep to reach dthydro
149 dtlast = dthydro - dtsub
150 nsub += 1
151
152 gij = solver_coag_k0_kdv(eps, nbins, massgrid, gij, tensor_tabflux_coag, dv, dtlast)
153
154 # when coagulation CFL > hydro timstep
155 else:
156 ndt += 1
157
158 gij = solver_coag_k0_kdv(eps, nbins, massgrid, gij, tensor_tabflux_coag, dv, dt)
159
160 return gij, nsub, ndt
161
162
163def compute_coag(
164 eps,
165 coeff_CFL,
166 nbins,
167 kpol,
168 massgrid,
169 massbins,
170 gij,
171 tensor_tabflux_coag,
172 tensor_tabintflux_coag,
173 dthydro,
174):
175 """
176 Function to compute coagulation solver for 1 hydro time-step
177
178 DG scheme k>0, piecewise polynomial (order k) approximation
179
180 Parameters
181 ----------
182 eps : scalar, type -> float
183 minimum value for mass distribution approximation gij
184 coeff_CFL : scalar, type -> float
185 timestep coefficient for stability of the SSPRK order 3 scheme
186 nbins : scalar, type -> integer
187 number of dust bins
188 kpol : scalar, type -> integer
189 degree of polynomials for approximation
190 massgrid : 1D array (dim = nbins+1), type -> float
191 grid of masses given borders value of mass bins
192 massbins : 1D array (dim = nbins), type -> float
193 arithmetic mean of massgrid
194 gij : 2D array (dim = (nbins,kpol+1)), type -> float
195 components of g on the polynomial basis
196 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
197 array to evaluate coagulation flux
198 tensor_tabintflux_coag : 6D array (dim = (nbins,kpol+1,nbins,nbins,kpol+1,kpol+1)), type -> float
199 array to evaluate term including the integral of the coagulation flux
200 dthydro : scalar, type -> float
201 hydro timestep, used as timestep to reach for coagulation process
202
203
204 Returns
205 -------
206 gij : 2D array (dim = (nbins,kpol+1)), type -> float
207 evolve components of g on the polynomial basis after 1 hydro timestep
208 nsub : scalar, type -> integer
209 number of subcycling coagulation timestep to reach dthydro
210 ndt : scalar, type -> integer
211 number of hydro timestep, when coagulation CFL > dthydro
212
213 """
214
215 nsub = 0
216 ndt = 0
217
218 # evaluate coagulation CFL
219 dtCFLsub = compute_CFL_coag(eps, nbins, kpol, massgrid, gij, tensor_tabflux_coag)
220 dtCFLsub = coeff_CFL * dtCFLsub
221 # print("dtCFLsub = ",dtCFLsub)
222
223 # compare hydro timestep and coagulation CFL
224 dt = min(dtCFLsub, dthydro)
225
226 # coagulation subcycling timesteps
227 if dt < dthydro:
228 dtsub = 0.0
229 while dtsub < dthydro and dthydro - dtsub > dtCFLsub:
230 dtsub += dtCFLsub
231
232 nsub += 1
233
234 gij = solver_coag(
235 eps,
236 nbins,
237 kpol,
238 massgrid,
239 massbins,
240 gij,
241 tensor_tabflux_coag,
242 tensor_tabintflux_coag,
243 dtCFLsub,
244 )
245
246 dtCFLsub = compute_CFL_coag(eps, nbins, kpol, massgrid, gij, tensor_tabflux_coag)
247 dtCFLsub = coeff_CFL * dtCFLsub
248 # print("dtCFLsub = ",dtCFLsub)
249
250 # last timestep to reach dthydro
251 dtlast = dthydro - dtsub
252 nsub += 1
253
254 gij = solver_coag(
255 eps,
256 nbins,
257 kpol,
258 massgrid,
259 massbins,
260 gij,
261 tensor_tabflux_coag,
262 tensor_tabintflux_coag,
263 dtlast,
264 )
265
266 # when coagulation CFL > hydro timstep
267 else:
268 ndt += 1
269
270 gij = solver_coag(
271 eps,
272 nbins,
273 kpol,
274 massgrid,
275 massbins,
276 gij,
277 tensor_tabflux_coag,
278 tensor_tabintflux_coag,
279 dt,
280 )
281
282 return gij, nsub, ndt
283
284
285def compute_coag_kdv(
286 eps,
287 coeff_CFL,
288 nbins,
289 kpol,
290 massgrid,
291 massbins,
292 gij,
293 tensor_tabflux_coag,
294 tensor_tabintflux_coag,
295 dv,
296 dthydro,
297):
298 """
299 Function to compute coagulation solver for 1 hydro time-step
300
301 DG scheme k>0, piecewise polynomial (order k) approximation
302
303 Function for ballistic kernel with differential velocities dv
304
305 Parameters
306 ----------
307 eps : scalar, type -> float
308 minimum value for mass distribution approximation gij
309 coeff_CFL : scalar, type -> float
310 timestep coefficient for stability of the SSPRK order 3 scheme
311 nbins : scalar, type -> integer
312 number of dust bins
313 kpol : scalar, type -> integer
314 degree of polynomials for approximation
315 massgrid : 1D array (dim = nbins+1), type -> float
316 grid of masses given borders value of mass bins
317 massbins : 1D array (dim = nbins), type -> float
318 arithmetic mean of massgrid
319 gij : 2D array (dim = (nbins,kpol+1)), type -> float
320 components of g on the polynomial basis
321 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
322 array to evaluate coagulation flux
323 tensor_tabintflux_coag : 6D array (dim = (nbins,kpol+1,nbins,nbins,kpol+1,kpol+1)), type -> float
324 array to evaluate term including the integral of the coagulation flux
325 dv : 2D array (dim = (nbins,nbins)), type -> float
326 array of the differential velocity between grains
327 dthydro : scalar, type -> float
328 hydro timestep, used as timestep to reach for coagulation process
329
330
331 Returns
332 -------
333 gij : 2D array (dim = (nbins,kpol+1)), type -> float
334 evolve components of g on the polynomial basis after 1 hydro timestep
335 nsub : scalar, type -> integer
336 number of subcycling coagulation timestep to reach dthydro
337 ndt : scalar, type -> integer
338 number of hydro timestep, when coagulation CFL > dthydro
339
340 """
341
342 nsub = 0
343 ndt = 0
344
345 # evaluate coagulation CFL
346 dtCFLsub = compute_CFL_coag_kdv(eps, nbins, kpol, massgrid, gij, tensor_tabflux_coag, dv)
347 dtCFLsub = coeff_CFL * dtCFLsub
348 # print("dtCFLsub = ",dtCFLsub)
349
350 # compare hydro timestep and coagulation CFL
351 dt = min(dtCFLsub, dthydro)
352
353 # coagulation subcycling timesteps
354 if dt < dthydro:
355 dtsub = 0.0
356 while dtsub < dthydro and dthydro - dtsub > dtCFLsub:
357 dtsub += dtCFLsub
358
359 nsub += 1
360
361 gij = solver_coag_kdv(
362 eps,
363 nbins,
364 kpol,
365 massgrid,
366 massbins,
367 gij,
368 tensor_tabflux_coag,
369 tensor_tabintflux_coag,
370 dv,
371 dtCFLsub,
372 )
373
374 dtCFLsub = compute_CFL_coag_kdv(
375 eps, nbins, kpol, massgrid, gij, tensor_tabflux_coag, dv
376 )
377 dtCFLsub = coeff_CFL * dtCFLsub
378 # print("dtCFLsub = ",dtCFLsub)
379
380 # last timestep to reach dthydro
381 dtlast = dthydro - dtsub
382 nsub += 1
383
384 gij = solver_coag_kdv(
385 eps,
386 nbins,
387 kpol,
388 massgrid,
389 massbins,
390 gij,
391 tensor_tabflux_coag,
392 tensor_tabintflux_coag,
393 dv,
394 dtlast,
395 )
396
397 # when coagulation CFL > hydro timstep
398 else:
399 ndt += 1
400
401 gij = solver_coag_kdv(
402 eps,
403 nbins,
404 kpol,
405 massgrid,
406 massbins,
407 gij,
408 tensor_tabflux_coag,
409 tensor_tabintflux_coag,
410 dv,
411 dt,
412 )
413
414 return gij, nsub, ndt
Generated on Wed May 13 2026 13:21:48 for Shamrock by doxygen 1.9.8