Dust Coagulation v2 #92
Conversation
…he 4D par_for_outer.
Coagulation Refactor
| const Real &gdens = vmesh(b, gas::prim::density(0), k, j, i); | ||
| const Real &gsie = vmesh(b, gas::prim::sie(0), k, j, i); | ||
| const Real &gbulk = eos_d.BulkModulusFromDensityInternalEnergy(gdens, gsie); | ||
| const Real cs1 = std::sqrt(gbulk / gdens) * vel0; |
There was a problem hiding this comment.
Sound speed vs thermal speed? --- From #91 @adamdempsey90
There was a problem hiding this comment.
@shengtai can you advise here?
There was a problem hiding this comment.
The coagulation algorithm fromTil Birnstiel uses the sound speed to calculate the temperature in their code. The sound speed should be used here. I am not sure what is the difference between thermal speed and sound speed
There was a problem hiding this comment.
If the sound speed is used to compute temperature, we should just pass the temperature. If the sound speed is used for calculating a scale height, we should reformulate that to not use the scale height. If the sound speed is used to calculate a viscosity, we should pass the viscosity structs instead.
We should only be using this sound speed of something actually needs it. If anything uses PPD quantities (surface density, scale height, alpha viscosity, etc.), we need to replace them with their more generic counterparts (e.g.,
There was a problem hiding this comment.
If the sound speed is used to compute temperature, we should just pass the temperature. If the sound speed is used for calculating a scale height, we should reformulate that to not use the scale height. If the sound speed is used to calculate a viscosity, we should pass the viscosity structs instead.
We should only be using this sound speed of something actually needs it. If anything uses PPD quantities (surface density, scale height, alpha viscosity, etc.), we need to replace them with their more generic counterparts (e.g.,
There was a problem hiding this comment.
I checked their coagulation kernel and found c_s (sound speed) is used to calculate the scale height, temperature, relative velocities, etc. We have two options: pass all informations to coagulation kernel: sound speed, scale height, and temperature using our calculated value, or pass only sound speed and let coagulation kernel to calculate other info needed. I prefer to pass only the sound speed.
To my knowledge, Til's algorithm is designed for proto-planetary disk applications. That is why it needs alpha, omega, and scale height etc. Its temperature calculation may also differ our temperature calculation in artemis. If we need to rewrite coagulation for general purpose (not only for disk applications), we can change their algorithm for other applications.
| } | ||
|
|
||
| if (i - pgrid - 1 >= 0) { | ||
| Kokkos::atomic_add(&source(i - 1), sum0); |
There was a problem hiding this comment.
ScatterView --- From #91 @adamdempsey90
There was a problem hiding this comment.
The i thread is updating the i-1 thread with the sum it created. Can this i loop be refactored so that thread i updates source(i) instead? I think the answer is yes and you just replace i with i+1 in the sum calculation above?
| for (int j = 0; j <= i; j++) { | ||
| coag3d(iphifrag, j, i) = std::pow(mass_grid(j), frag_slope); | ||
| sum_pF += coag3d(iphifrag, j, i); | ||
| } | ||
| // normalization | ||
| for (int j = 0; j <= i; j++) { | ||
| coag3d(iphifrag, j, i) /= sum_pF; // switch (i,j) from fortran | ||
| } |
There was a problem hiding this comment.
I'm seeing a mix of (j,i) and (i,j) usage throughout this file. Should there be one common indexing pattern, or do you really want to reverse the order for some?
| parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, mbr, 0, mimax, [&](const int i) { | ||
| for (int j = 0; j <= i; j++) { | ||
| // Calculate the Rate | ||
| const Real fett_t = CoagulationRate(i, j, kernel4, vel, stoppingTime, coag, 0); | ||
| const Real Rc1 = distri(i) * distri(j) * fett_t; | ||
| for (int nz = 0; nz < 4; nz++) { | ||
| const int k = coag.cpod_notzero(i, j, nz); | ||
| if (k < 0) continue; | ||
| const Real src_coag0 = coag.cpod_short(i, j, nz) * Rc1; | ||
| Kokkos::atomic_add(&source(k), src_coag0); | ||
| } | ||
| } | ||
| }); |
There was a problem hiding this comment.
The atomic can be replaced with an inner reduce?
| } | ||
|
|
||
| if (i - pgrid - 1 >= 0) { | ||
| Kokkos::atomic_add(&source(i - 1), sum0); |
There was a problem hiding this comment.
The i thread is updating the i-1 thread with the sum it created. Can this i loop be refactored so that thread i updates source(i) instead? I think the answer is yes and you just replace i with i+1 in the sum calculation above?
| parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, mbr, 0, mimax, [&](const int i) { | ||
| for (int j = 0; j <= i; j++) { |
There was a problem hiding this comment.
could this be a 2D par_for_inner with a screen on j <= i?
| for (int nz = 0; nz < 4; nz++) { | ||
| const int k = coag.cpod_notzero(i, j, nz); | ||
| if (k < 0) continue; | ||
| const Real kdeltajk = (j == k) ? 1.0 : 0.0; | ||
| const Real kdeltaik = (i == k) ? 1.0 : 0.0; | ||
| const Real nqs_k = | ||
| (Qp1 * coag.cpod_short(i, j, nz) + tmp1 * kdeltajk + tmp2 * kdeltaik) * Rc1; | ||
| Kokkos::atomic_add(&nQs(k), nqs_k); | ||
| } |
There was a problem hiding this comment.
| for (int nz = 0; nz < 4; nz++) { | |
| const int k = coag.cpod_notzero(i, j, nz); | |
| if (k < 0) continue; | |
| const Real kdeltajk = (j == k) ? 1.0 : 0.0; | |
| const Real kdeltaik = (i == k) ? 1.0 : 0.0; | |
| const Real nqs_k = | |
| (Qp1 * coag.cpod_short(i, j, nz) + tmp1 * kdeltajk + tmp2 * kdeltaik) * Rc1; | |
| Kokkos::atomic_add(&nQs(k), nqs_k); | |
| } | |
| Real nqs_k = 0.0; | |
| for (int nz = 0; nz < 4; nz++) { | |
| const int k = coag.cpod_notzero(i, j, nz); | |
| if (k < 0) continue; | |
| const Real kdeltajk = (j == k) ? 1.0 : 0.0; | |
| const Real kdeltaik = (i == k) ? 1.0 : 0.0; | |
| nqs_k += | |
| (Qp1 * coag.cpod_short(i, j, nz) + tmp1 * kdeltajk + tmp2 * kdeltaik) * Rc1; | |
| } | |
| Kokkos::atomic_add(&nQs(k), nqs_k); |
This one should really be a ScatterView, but I don't know if you can do that for scratchpads?
| for (int j = 0; j <= i - pgrid - 1; j++) { | ||
| const Real fett_l = CoagulationRate(i, j, kernel4, vel, stoppingTime, coag, 1); | ||
| const Real Rf1 = distri(i) * distri(j) * fett_l; | ||
| const Real dummy = coag.coagR3D(iepsfrag, i, j) * Rf1 * Q(i); | ||
| sum0 += dummy; | ||
| } | ||
|
|
||
| if (i - pgrid - 1 >= 0) { | ||
| Kokkos::atomic_add(&nQs(i - 1), sum0); | ||
| } | ||
|
|
||
| Real sum1 = -sum0; | ||
| // Full fragmentation (only negative terms) | ||
| int i1 = std::max(0, i - pgrid); | ||
| for (int j = i1; j <= i; j++) { | ||
| const Real fett_l = CoagulationRate(i, j, kernel4, vel, stoppingTime, coag, 1); | ||
| const Real Rf1 = distri(i) * distri(j) * fett_l * Q(i); | ||
| sum1 -= Rf1; | ||
| } | ||
| Kokkos::atomic_add(&nQs(i), sum1); |
There was a problem hiding this comment.
should be able to remove all atomics by incrementing the i used in the first sum calculation
|
|
||
| CoagulationSource(mbr, source, dustdens, mimax, kernel, vel, stime, coag); | ||
|
|
||
| if (coag.use_adaptive == 0 || coag.integrator == 1) { |
There was a problem hiding this comment.
Can some of these options be pulled out of this while loop done for every cell and made into a template parameter for this function?
| Real errmax; | ||
| int mimax2 = 0; | ||
| // Heun's Method | ||
| while (1) { |
There was a problem hiding this comment.
Is there a sensible maximum number of iterations that we can put here to avoid an infinite loop?
Co-authored-by: Adam M. Dempsey <adempsey@lanl.gov>
Co-authored-by: Adam M. Dempsey <adempsey@lanl.gov>
Background
Description of Changes
Checklist