Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
solver_DG.py
1import sys
2import time
3
4import numpy as np
5from scipy.special import eval_legendre
6
7from .generate_flux_intflux import *
8from .limiter import *
9from .utils_polynomials import *
10
11
12# for analytical collision kernel
13def compute_CFL_coag_k0(eps, nbins, massgrid, gij, tensor_tabflux_coag):
14 """
15 Function to compute coagulation CFL for DG scheme k=0 piecewise constant approximation
16
17 CFL formulation from Filbet & Laurencot 2004, dt <= mean_g * dm/dF
18
19 Parameters
20 ----------
21 eps : scalar, type -> float
22 minimum value for mass distribution approximation gij
23 nbins : scalar, type -> integer
24 number of dust bins
25 massgrid : 1D array (dim = nbins+1), type -> float
26 grid of masses given borders value of mass bins
27 gij : 1D array (dim = nbins), type -> float
28 components of g on the polynomial basis
29 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
30 array to evaluate coagulation flux
31
32
33 Returns
34 -------
35 CFL_k0 : scalar, type -> float
36 CFL restriction for coagulation
37
38 """
39
40 tabdflux = np.zeros(nbins)
41 tabdtCFL = np.zeros(nbins)
42
43 flux = compute_flux_coag_k0(gij, tensor_tabflux_coag)
44
45 if np.array_equal(flux, np.zeros(nbins)):
46 # calculs continue with no coag
47 CFL_k0 = 1e3
48
49 else:
50 tabdflux[0] = flux[0]
51 tabdtCFL[0] = np.abs(gij[0] * (massgrid[1] - massgrid[0]) / (tabdflux[0]))
52
53 for j in range(1, nbins):
54 hj = massgrid[j + 1] - massgrid[j]
55 tabdflux[j] = flux[j] - flux[j - 1]
56
57 if gij[j] > eps:
58 tabdtCFL[j] = np.abs(gij[j] * hj / tabdflux[j])
59 else:
60 tabdtCFL[j] = 1e3
61
62 CFL_k0 = np.min(tabdtCFL)
63
64 if CFL_k0 == 0.0:
65 print("gij=", gij)
66 print("hj = ", massgrid[1:] - massgrid[:nbins])
67 print("tabdtCFL = ", tabdtCFL)
68 print("CFL_k0 = ", CFL_k0)
69 print("issue in CFL coagulation")
70 sys.exit()
71
72 return CFL_k0
73
74
75def L_func_coag_k0(nbins, massgrid, gij, tensor_tabflux_coag):
76 """
77 Function to compute the DG operator L for piecewise constant approximation (see Lombart et al., 2021)
78 It is used for the time solver
79
80 Parameters
81 ----------
82 nbins : scalar, type -> integer
83 number of dust bins
84 massgrid : 1D array (dim = nbins+1), type -> float
85 grid of masses given borders value of mass bins
86 gij : 1D array (dim = nbins), type -> float
87 components of g on the polynomial basis
88 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
89 array to evaluate coagulation flux
90
91
92 Returns
93 -------
94 Lk0 : 1D array (dim = nbins), type -> float
95 DG operator for piecewise constant approximation in each bin
96
97 """
98
99 flux = compute_flux_coag_k0(gij, tensor_tabflux_coag)
100
101 Lk0 = np.zeros(nbins)
102 Lk0[0] = -flux[0] / (massgrid[1] - massgrid[0])
103 Lk0[1:] = -(flux[1:] - flux[0 : nbins - 1]) / (massgrid[2:] - massgrid[1:nbins])
104
105 return Lk0
106
107
108def solver_coag_k0(eps, nbins, massgrid, gij, tensor_tabflux_coag, dt):
109 """
110 Function to compute SSPRK order 3 time solver with piecewise constant approximation
111 See Zhang & Shu 2010 and Lombart et al., 2021
112
113
114 Parameters
115 ----------
116 eps : scalar, type -> float
117 minimum value for mass distribution approximation gij
118 nbins : scalar, type -> integer
119 number of dust bins
120 massgrid : 1D array (dim = nbins+1), type -> float
121 grid of masses given borders value of mass bins
122 gij : 1D array (dim = nbins), type -> float
123 components of g on the polynomial basis
124 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
125 array to evaluate coagulation flux
126 dt : scalar, type -> float
127 timestep
128
129
130 Returns
131 -------
132 gijnew : 1D array (dim = nbins), type -> float
133 evolved components of g on the polynomial basis after 1 timestep
134
135 """
136
137 # SSPRK3 algo (Zhang & Shu 2010)
138 # step 1
139 Lk0 = L_func_coag_k0(nbins, massgrid, gij, tensor_tabflux_coag)
140
141 gij1 = gij + dt * Lk0
142
143 # limit gij1 to eps value
144 if len(gij1[gij1 < 0.0]) > 0:
145 print("gij", gij)
146 print("dt*Lk0 = ", dt * Lk0)
147 print("gij1 = ", gij1)
148 sys.exit()
149 else:
150 gij1[gij1 < eps] = eps
151
152 Lk0_1 = L_func_coag_k0(nbins, massgrid, gij1, tensor_tabflux_coag)
153
154 gij2 = 3.0 * gij / 4.0 + (gij1 + dt * Lk0_1) / 4.0
155
156 # limit gij2 to eps value
157 if len(gij2[gij2 < 0.0]) > 0:
158 print("gi1", gi1)
159 print("dt*Lk0_1 = ", dt * Lk0_1)
160 print("gij2 = ", gij2)
161 sys.exit()
162 else:
163 gij2[gij2 < eps] = eps
164
165 # step 3
166 Lk0_2 = L_func_coag_k0(nbins, massgrid, gij2, tensor_tabflux_coag)
167
168 gijnew = gij / 3.0 + 2.0 * (gij2 + dt * Lk0_2) / 3.0
169
170 # limit gijnew to eps value
171 if len(gijnew[gijnew < 0.0]) > 0:
172 print("gij2", gij2)
173 print("dt*Lk0_2 = ", dt * Lk0_2)
174 print("gijnew = ", gijnew)
175 sys.exit()
176 else:
177 gijnew[gijnew < eps] = eps
178
179 return gijnew
180
181
182# for ballistic collision kernel with non analytical dv
183def compute_CFL_coag_k0_kdv(eps, nbins, massgrid, gij, tensor_tabflux_coag, dv):
184 """
185 Function to compute coagulation CFL for DG scheme k=0 piecewise constant approximation.
186
187 Function for ballistic kernel with differential velocities dv
188
189 CFL formulation from Filbet & Laurencot 2004, dt <= mean_g * dm/dF
190
191 Parameters
192 ----------
193 eps : scalar, type -> float
194 minimum value for mass distribution approximation gij
195 nbins : scalar, type -> integer
196 number of dust bins
197 massgrid : 1D array (dim = nbins+1), type -> float
198 grid of masses given borders value of mass bins
199 gij : 1D array (dim = nbins), type -> float
200 components of g on the polynomial basis
201 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
202 array to evaluate coagulation flux
203 dv : 2D array (dim = (nbins,nbins)), type -> float
204 array of the differential velocity between grains
205
206
207 Returns
208 -------
209 CFL_k0 : scalar, type -> float
210 CFL restriction for coagulation
211
212 """
213
214 tabdflux = np.zeros(nbins)
215 tabdtCFL = np.zeros(nbins)
216
217 flux = compute_flux_coag_k0_kdv(gij, tensor_tabflux_coag, dv)
218
219 if np.array_equal(flux, np.zeros(nbins)):
220 # calculs continue with no coag
221 CFL_k0 = 1e3
222
223 else:
224 tabdflux[0] = flux[0]
225 tabdtCFL[0] = np.abs(gij[0] * (massgrid[1] - massgrid[0]) / (tabdflux[0]))
226
227 for j in range(1, nbins):
228 hj = massgrid[j + 1] - massgrid[j]
229 tabdflux[j] = flux[j] - flux[j - 1]
230
231 if gij[j] > eps:
232 tabdtCFL[j] = np.abs(gij[j] * hj / tabdflux[j])
233 else:
234 tabdtCFL[j] = 1e3
235
236 # print("gij=",gij)
237 # print("tabdtCFL =",tabdtCFL)
238
239 CFL_k0 = np.min(tabdtCFL)
240
241 if CFL_k0 == 0.0:
242 print("gij=", gij)
243 print("hj = ", massgrid[1:] - massgrid[:nbins])
244 print("tabdtCFL = ", tabdtCFL)
245 print("CFL_k0 = ", CFL_k0)
246 print("issue in CFL coagulation")
247 sys.exit()
248
249 return CFL_k0
250
251
252# for ballistic collision kernel with non analytical dv
253def L_func_coag_k0_kdv(nbins, massgrid, gij, tensor_tabflux_coag, dv):
254 """
255 Function to compute the DG operator L for piecewise constant approximation (see Lombart et al., 2021)
256 Function for ballistic kernel with differential velocities dv
257 It is used for the time solver
258
259 Parameters
260 ----------
261 nbins : scalar, type -> integer
262 number of dust bins
263 massgrid : 1D array (dim = nbins+1), type -> float
264 grid of masses given borders value of mass bins
265 gij : 1D array (dim = nbins), type -> float
266 components of g on the polynomial basis
267 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
268 array to evaluate coagulation flux
269 dv : 2D array (dim = (nbins,nbins)), type -> float
270 array of the differential velocity between grains
271
272
273 Returns
274 -------
275 Lk0 : 1D array (dim = nbins), type -> float
276 DG operator for piecewise constant approximation in each bin
277
278 """
279
280 flux = compute_flux_coag_k0_kdv(gij, tensor_tabflux_coag, dv)
281
282 Lk0 = np.zeros(nbins)
283 Lk0[0] = -flux[0] / (massgrid[1] - massgrid[0])
284 Lk0[1:] = -(flux[1:] - flux[0 : nbins - 1]) / (massgrid[2:] - massgrid[1:nbins])
285
286 return Lk0
287
288
289# for ballistic collision kernel with non analytical dv
290def solver_coag_k0_kdv(eps, nbins, massgrid, gij, tensor_tabflux_coag, dv, dt):
291 """
292 Function to compute SSPRK order 3 time solver with piecewise constant approximation
293 Function for ballistic kernel with differential velocities dv
294 See Zhang & Shu 2010 and Lombart et al., 2021
295
296
297 Parameters
298 ----------
299 eps : scalar, type -> float
300 minimum value for mass distribution approximation gij
301 nbins : scalar, type -> integer
302 number of dust bins
303 massgrid : 1D array (dim = nbins+1), type -> float
304 grid of masses given borders value of mass bins
305 gij : 1D array (dim = nbins), type -> float
306 components of g on the polynomial basis
307 tensor_tabflux_coag : 3D array (dim = (nbins,nbins,nbins)), type -> float
308 array to evaluate coagulation flux
309 dv : 2D array (dim = (nbins,nbins)), type -> float
310 array of the differential velocity between grains
311 dt : scalar, type -> float
312 timestep
313
314
315 Returns
316 -------
317 gijnew : 1D array (dim = nbins), type -> float
318 evolved components of g on the polynomial basis after 1 timestep
319
320 """
321
322 # SSPRK3 algo (Zhang & Shu 2010)
323 # step 1
324 Lk0 = L_func_coag_k0_kdv(nbins, massgrid, gij, tensor_tabflux_coag, dv)
325
326 gij1 = gij + dt * Lk0
327
328 # limit gij1 to eps value
329 if len(gij1[gij1 < 0.0]) > 0:
330 print("gij", gij)
331 print("dt*Lk0 = ", dt * Lk0)
332 print("gij1 = ", gij1)
333 sys.exit()
334 else:
335 gij1[gij1 < eps] = eps
336
337 Lk0_1 = L_func_coag_k0_kdv(nbins, massgrid, gij1, tensor_tabflux_coag, dv)
338
339 gij2 = 3.0 * gij / 4.0 + (gij1 + dt * Lk0_1) / 4.0
340
341 # limit gij2 to eps value
342 if len(gij2[gij2 < 0.0]) > 0:
343 print("gi1", gi1)
344 print("dt*Lk0_1 = ", dt * Lk0_1)
345 print("gij2 = ", gij2)
346 sys.exit()
347 else:
348 gij2[gij2 < eps] = eps
349
350 # step 3
351 Lk0_2 = L_func_coag_k0_kdv(nbins, massgrid, gij2, tensor_tabflux_coag, dv)
352
353 gijnew = gij / 3.0 + 2.0 * (gij2 + dt * Lk0_2) / 3.0
354
355 # limit gijnew to eps value
356 if len(gijnew[gijnew < 0.0]) > 0:
357 print("gij2", gij2)
358 print("dt*Lk0_2 = ", dt * Lk0_2)
359 print("gijnew = ", gijnew)
360 sys.exit()
361 else:
362 gijnew[gijnew < eps] = eps
363
364 return gijnew
365
366
367def compute_CFL_coag(eps, nbins, kpol, massgrid, gij, tensor_tabflux_coag):
368 """
369 Function to compute coagulation CFL for DG scheme k>0 piecewise polynomial approximation
370
371 CFL formulation from Filbet & Laurencot 2004, dt <= mean_g * dm/dF
372
373 Parameters
374 ----------
375 eps : scalar, type -> float
376 minimum value for mass distribution approximation gij
377 nbins : scalar, type -> integer
378 number of dust bins
379 kpol : scalar, type -> integer
380 degree of polynomials for approximation
381 massgrid : 1D array (dim = nbins+1), type -> float
382 grid of masses given borders value of mass bins
383 gij : 2D array (dim = (nbins,kpol+1)), type -> float
384 components of g on the polynomial basis
385 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
386 array to evaluate coagulation flux
387
388
389 Returns
390 -------
391 CFL : scalar, type -> float
392 CFL restriction for coagulation
393
394 """
395
396 tabdflux = np.zeros(nbins)
397 tabdtCFL = np.zeros(nbins)
398
399 flux = compute_flux_coag(gij, tensor_tabflux_coag)
400
401 if np.array_equal(flux, np.zeros(nbins)):
402 # calculs continue with no coag
403 CFL_k0 = 1e3
404
405 else:
406 tabdflux[0] = flux[0]
407 tabdtCFL[0] = np.abs(gij[0, 0] * (massgrid[1] - massgrid[0]) / (tabdflux[0]))
408
409 for j in range(1, nbins):
410 hj = massgrid[j + 1] - massgrid[j]
411 tabdflux[j] = flux[j] - flux[j - 1]
412
413 if gij[j, 0] > eps:
414 tabdtCFL[j] = np.abs(gij[j, 0] * hj / tabdflux[j])
415 else:
416 tabdtCFL[j] = 1e3
417
418 CFL = np.min(tabdtCFL)
419
420 if CFL <= 0.0:
421 print("gij[:,0]=", gij[:, 0])
422 print("hj = ", massgrid[1:] - massgrid[: nbins + 1])
423 print("tabdtCFL = ", tabdtCFL)
424 print("CFL = ", CFL)
425 print("issue in CFL")
426 sys.exit()
427
428 return CFL
429
430
431def L_func_coag(nbins, kpol, massgrid, gij, tensor_tabflux_coag, tensor_tabintflux_coag):
432 """
433 Function to compute the DG operator L for piecewise polynomial approximation (see Lombart et al., 2021)
434 It is used for the time solver
435
436 Parameters
437 ----------
438 nbins : scalar, type -> integer
439 number of dust bins
440 kpol : scalar, type -> integer
441 degree of polynomials for approximation
442 massgrid : 1D array (dim = nbins+1), type -> float
443 grid of masses given borders value of mass bins
444 gij : 2D array (dim = (nbins,kpol+1)), type -> float
445 components of g on the polynomial basis
446 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
447 array to evaluate coagulation flux
448 tensor_tabintflux_coag : 6D array (dim = (nbins,kpol+1,nbins,nbins,kpol+1,kpol+1)), type -> float
449 array to evaluate the term including the integral of coagulation flux
450
451
452 Returns
453 -------
454 Lk : 2D array (dim = (nbins,kpol+1)), type -> float
455 DG operator for piecewise polynomials approximation in each bin
456
457 """
458
459 flux = compute_flux_coag(gij, tensor_tabflux_coag)
460 intflux = compute_intflux_coag(gij, tensor_tabintflux_coag)
461
462 mat_intflux = intflux.reshape(nbins, kpol + 1)
463
464 LegPright = eval_legendre(range(kpol + 1), 1.0)
465 flux_right = np.outer(flux, LegPright)
466
467 flux_left = np.zeros((nbins, kpol + 1))
468 LegPleft = eval_legendre(range(kpol + 1), -1.0)
469 flux_left[1:, :] = np.outer(flux[0 : nbins - 1], LegPleft)
470
471 mat_coeff_norm = np.outer(
472 2.0 / (massgrid[1:] - massgrid[0:nbins]), 1.0 / coeff_norm_leg(np.arange(kpol + 1))
473 )
474
475 Lk = mat_coeff_norm * (mat_intflux - (flux_right - flux_left))
476
477 return Lk
478
479
480def solver_coag(
481 eps, nbins, kpol, massgrid, massbins, gij, tensor_tabflux_coag, tensor_tabintflux_coag, dt
482):
483 """
484 Function to compute SSPRK order 3 time solver with piecewise polynomial approximation
485 See Zhang & Shu 2010 and Lombart et al., 2021
486
487
488 Parameters
489 ----------
490 eps : scalar, type -> float
491 minimum value for mass distribution approximation gij
492 nbins : scalar, type -> integer
493 number of dust bins
494 kpol : scalar, type -> integer
495 degree of polynomials for approximation
496 massgrid : 1D array (dim = nbins+1), type -> float
497 grid of masses given borders value of mass bins
498 massbins : 1D array (dim = nbins), type -> float
499 arithmetic mean value of massgrid for each mass bins
500 gij : 2D array (dim = (nbins,kpol+1)), type -> float
501 components of g on the polynomial basis
502 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
503 array to evaluate coagulation flux
504 tensor_tabintflux_coag : 6D array (dim = (nbins,kpol+1,nbins,nbins,kpol+1,kpol+1)), type -> float
505 array to evaluate the term including the integral of coagulation flux
506 dt : scalar, type -> float
507 timestep
508
509
510 Returns
511 -------
512 gijnew : 2D array (dim = (nbins,kpol+1)), type -> float
513 evolved components of g on the polynomial basis after 1 timestep
514
515 """
516
517 # SSPRK3 algo (Zhang & Shu 2010)
518 # step 1
519 Lk = L_func_coag(nbins, kpol, massgrid, gij, tensor_tabflux_coag, tensor_tabintflux_coag)
520
521 gij1 = gij + dt * Lk
522
523 # apply limiter coefficient
524 tabgamma = gammafunction(eps, nbins, kpol, massgrid, massbins, gij1)
525 gij1[:, 1:] = tabgamma.reshape(nbins, 1) * gij1[:, 1:]
526
527 # limit to eps value
528 if gij1[:, 0][gij1[:, 0] < 0.0].any():
529 print("gij", gij)
530 print("dt*Lk = ", dt * Lk)
531 print("gij1 = ", gij1)
532 sys.exit()
533 else:
534 gij1[:, 1:][gij1[:, 0] < eps] = 0.0
535 gij1[:, 0][gij1[:, 0] < eps] = eps
536
537 Lk_1 = L_func_coag(nbins, kpol, massgrid, gij1, tensor_tabflux_coag, tensor_tabintflux_coag)
538
539 gij2 = 3.0 * gij / 4.0 + (gij1 + dt * Lk_1) / 4.0
540
541 # apply limiter coefficient
542 tabgamma = gammafunction(eps, nbins, kpol, massgrid, massbins, gij2)
543 gij2[:, 1:] = tabgamma.reshape(nbins, 1) * gij2[:, 1:]
544
545 # limit to eps value
546
547 if gij2[:, 0][gij2[:, 0] < 0.0].any():
548 print("gij", gij)
549 print("dt*Lk_1 = ", dt * Lk_1)
550 print("gij2 = ", gij2)
551 sys.exit()
552 else:
553 gij2[:, 1:][gij2[:, 0] < eps] = 0.0
554 gij2[:, 0][gij2[:, 0] < eps] = eps
555
556 # step 3
557 Lk_2 = L_func_coag(nbins, kpol, massgrid, gij2, tensor_tabflux_coag, tensor_tabintflux_coag)
558
559 gijnew = gij / 3.0 + 2.0 * (gij2 + dt * Lk_2) / 3.0
560
561 # apply limiter coefficient
562 tabgamma = gammafunction(eps, nbins, kpol, massgrid, massbins, gijnew)
563 gijnew[:, 1:] = tabgamma.reshape(nbins, 1) * gijnew[:, 1:]
564
565 # limit to eps value
566
567 if gijnew[:, 0][gijnew[:, 0] < 0.0].any():
568 print("gij", gij)
569 print("dt*Lk_1 = ", dt * Lk_1)
570 print("gijnew = ", gijnew)
571 sys.exit()
572 else:
573 gijnew[:, 1:][gijnew[:, 0] < eps] = 0.0
574 gijnew[:, 0][gijnew[:, 0] < eps] = eps
575
576 return gijnew
577
578
579# for ballistic collision kernel with non analytical dv
580def compute_CFL_coag_kdv(eps, nbins, kpol, massgrid, gij, tensor_tabflux_coag, dv):
581 """
582 Function to compute coagulation CFL for DG scheme k>0 piecewise polynomial approximation
583 Function for ballistic kernel with differential velocities dv
584 CFL formulation from Filbet & Laurencot 2004, dt <= mean_g * dm/dF
585
586 Parameters
587 ----------
588 eps : scalar, type -> float
589 minimum value for mass distribution approximation gij
590 nbins : scalar, type -> integer
591 number of dust bins
592 kpol : scalar, type -> integer
593 degree of polynomials for approximation
594 massgrid : 1D array (dim = nbins+1), type -> float
595 grid of masses given borders value of mass bins
596 gij : 2D array (dim = (nbins,kpol+1)), type -> float
597 components of g on the polynomial basis
598 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
599 array to evaluate coagulation flux
600 dv : 2D array (dim = (nbins,nbins)), type -> float
601 array of the differential velocity between grains
602
603
604 Returns
605 -------
606 CFL : scalar, type -> float
607 CFL restriction for coagulation
608
609 """
610 tabdflux = np.zeros(nbins)
611 tabdtCFL = np.zeros(nbins)
612
613 flux = compute_flux_coag_kdv(gij, tensor_tabflux_coag, dv)
614
615 if np.array_equal(flux, np.zeros(nbins)):
616 # calculs continue with no coag
617 CFL = 1e3
618
619 else:
620 tabdflux[0] = flux[0]
621 tabdtCFL[0] = np.abs(gij[0, 0] * (massgrid[1] - massgrid[0]) / (tabdflux[0]))
622
623 for j in range(1, nbins):
624 hj = massgrid[j + 1] - massgrid[j]
625 tabdflux[j] = flux[j] - flux[j - 1]
626
627 if gij[j, 0] > eps:
628 tabdtCFL[j] = np.abs(gij[j, 0] * hj / tabdflux[j])
629 else:
630 tabdtCFL[j] = 1e3
631
632 CFL = np.min(tabdtCFL)
633
634 if CFL <= 0.0:
635 print("gij[:,0]=", gij[:, 0])
636 print("hj = ", massgrid[1:] - massgrid[: nbins + 1])
637 print("tabdtCFL = ", tabdtCFL)
638 print("CFL = ", CFL)
639 print("issue in CFL")
640 sys.exit()
641
642 return CFL
643
644
645# for ballistic collision kernel with non analytical dv
646def L_func_coag_kdv(nbins, kpol, massgrid, gij, tensor_tabflux_coag, tensor_tabintflux_coag, dv):
647 """
648 Function to compute the DG operator L for piecewise polynomial approximation (see Lombart et al., 2021)
649 Function for ballistic kernel with differential velocities dv
650 It is used for the time solver
651
652 Parameters
653 ----------
654 nbins : scalar, type -> integer
655 number of dust bins
656 kpol : scalar, type -> integer
657 degree of polynomials for approximation
658 massgrid : 1D array (dim = nbins+1), type -> float
659 grid of masses given borders value of mass bins
660 gij : 2D array (dim = (nbins,kpol+1)), type -> float
661 components of g on the polynomial basis
662 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
663 array to evaluate coagulation flux
664 tensor_tabintflux_coag : 6D array (dim = (nbins,kpol+1,nbins,nbins,kpol+1,kpol+1)), type -> float
665 array to evaluate the term including the integral of coagulation flux
666 dv : 2D array (dim = (nbins,nbins)), type -> float
667 array of the differential velocity between grains
668
669 Returns
670 -------
671 Lk : 2D array (dim = (nbins,kpol+1)), type -> float
672 DG operator for piecewise polynomials approximation in each bin
673
674 """
675
676 flux = compute_flux_coag_kdv(gij, tensor_tabflux_coag, dv)
677 intflux = compute_intflux_coag_kdv(gij, tensor_tabintflux_coag, dv)
678
679 mat_intflux = intflux.reshape(nbins, kpol + 1)
680
681 LegPright = eval_legendre(range(kpol + 1), 1.0)
682 flux_right = np.outer(flux, LegPright)
683
684 flux_left = np.zeros((nbins, kpol + 1))
685 LegPleft = eval_legendre(range(kpol + 1), -1.0)
686 flux_left[1:, :] = np.outer(flux[0 : nbins - 1], LegPleft)
687
688 mat_coeff_norm = np.outer(
689 2.0 / (massgrid[1:] - massgrid[0:nbins]), 1.0 / coeff_norm_leg(np.arange(kpol + 1))
690 )
691
692 Lk = mat_coeff_norm * (mat_intflux - (flux_right - flux_left))
693
694 return Lk
695
696
697# for ballistic collision kernel with non analytical dv
698def solver_coag_kdv(
699 eps, nbins, kpol, massgrid, massbins, gij, tensor_tabflux_coag, tensor_tabintflux_coag, dv, dt
700):
701 """
702 Function to compute SSPRK order 3 time solver with piecewise constant approximation
703 Function for ballistic kernel with differential velocities dv
704 See Zhang & Shu 2010 and Lombart et al., 2021
705
706
707 Parameters
708 ----------
709 eps : scalar, type -> float
710 minimum value for mass distribution approximation gij
711 nbins : scalar, type -> integer
712 number of dust bins
713 kpol : scalar, type -> integer
714 degree of polynomials for approximation
715 massgrid : 1D array (dim = nbins+1), type -> float
716 grid of masses given borders value of mass bins
717 massbins : 1D array (dim = nbins), type -> float
718 arithmetic mean value of massgrid for each mass bins
719 gij : 2D array (dim = (nbins,kpol+1)), type -> float
720 components of g on the polynomial basis
721 tensor_tabflux_coag : 5D array (dim = (nbins,nbins,nbins,kpol+1,kpol+1)), type -> float
722 array to evaluate coagulation flux
723 tensor_tabintflux_coag : 6D array (dim = (nbins,kpol+1,nbins,nbins,kpol+1,kpol+1)), type -> float
724 array to evaluate the term including the integral of coagulation flux
725 dv : 2D array (dim = (nbins,nbins)), type -> float
726 array of the differential velocity between grains
727 dt : scalar, type -> float
728 timestep
729
730
731 Returns
732 -------
733 gijnew : 2D array (dim = (nbins,kpol+1)), type -> float
734 evolved components of g on the polynomial basis after 1 timestep
735
736 """
737
738 # SSPRK3 algo (Zhang & Shu 2010)
739 # step 1
740 Lk = L_func_coag_kdv(
741 nbins, kpol, massgrid, gij, tensor_tabflux_coag, tensor_tabintflux_coag, dv
742 )
743
744 gij1 = gij + dt * Lk
745
746 # apply limiter coefficient
747 tabgamma = gammafunction(eps, nbins, kpol, massgrid, massbins, gij1)
748 gij1[:, 1:] = tabgamma.reshape(nbins, 1) * gij1[:, 1:]
749
750 # limit to eps value
751 if gij1[:, 0][gij1[:, 0] < 0.0].any():
752 print("gij", gij)
753 print("dt*Lk = ", dt * Lk)
754 print("gij1 = ", gij1)
755 print("Negative value ")
756 sys.exit()
757 else:
758 gij1[:, 1:][gij1[:, 0] < eps] = 0.0
759 gij1[:, 0][gij1[:, 0] < eps] = eps
760
761 Lk_1 = L_func_coag_kdv(
762 nbins, kpol, massgrid, gij1, tensor_tabflux_coag, tensor_tabintflux_coag, dv
763 )
764
765 gij2 = 3.0 * gij / 4.0 + (gij1 + dt * Lk_1) / 4.0
766
767 # apply limiter coefficient
768 tabgamma = gammafunction(eps, nbins, kpol, massgrid, massbins, gij2)
769 gij2[:, 1:] = tabgamma.reshape(nbins, 1) * gij2[:, 1:]
770
771 # limit to eps value
772
773 if gij2[:, 0][gij2[:, 0] < 0.0].any():
774 print("gij", gij)
775 print("dt*Lk_1 = ", dt * Lk_1)
776 print("gij2 = ", gij2)
777 print("Negative value ")
778 sys.exit()
779 else:
780 gij2[:, 1:][gij2[:, 0] < eps] = 0.0
781 gij2[:, 0][gij2[:, 0] < eps] = eps
782
783 # step 3
784 Lk_2 = L_func_coag_kdv(
785 nbins, kpol, massgrid, gij2, tensor_tabflux_coag, tensor_tabintflux_coag, dv
786 )
787
788 gijnew = gij / 3.0 + 2.0 * (gij2 + dt * Lk_2) / 3.0
789
790 # apply limiter coefficient
791 tabgamma = gammafunction(eps, nbins, kpol, massgrid, massbins, gijnew)
792 gijnew[:, 1:] = tabgamma.reshape(nbins, 1) * gijnew[:, 1:]
793
794 # limit to eps value
795
796 if gijnew[:, 0][gijnew[:, 0] < 0.0].any():
797 print("gij", gij)
798 print("dt*Lk_1 = ", dt * Lk_1)
799 print("gijnew = ", gijnew)
800 print("Negative value ")
801 sys.exit()
802 else:
803 gijnew[:, 1:][gijnew[:, 0] < eps] = 0.0
804 gijnew[:, 0][gijnew[:, 0] < eps] = eps
805
806 return gijnew
Generated on Wed May 13 2026 13:21:48 for Shamrock by doxygen 1.9.8