Pr/inject species ratio

[jcd496]

I'm not sure if you have other ideas for this implementation? this method maintains the correct electron_he/electron ratio, but maybe there is a cleaner way? Specifically I don't know if we should scale according to note on line 153 in setup_particles.hxx

[jcd496]

I see that his build failed, idk though. It builds fine on summit with spack

[wrfox]

Hi John, Kai - I read the code through. I believe the following statement is not the ideal way to do it - ... if(pop == HE_population){ init_npt(non_neutralizing_population, pos, p, {jx, jy, jz}, npt); npt.n *= HE_ratio; ... Better is for all the case-specific code should be in init_npt. Therefore, init_npt should return the correct density of HE particles to add. It should not need to be modified in setupParticles. Second, there should not need to be a “non_neutralizing_population” variable Second, I propose the pseudo-code should do the loop in the following order to handle the FIXME: // positive charges For each population pop if (pop == Neutralizing population) Continue call init_npt (pop, …) If npt.kind.q < 0 continue // else positive charge n_in_cell = get_n_in_cell ... n_q_in_cell += kinds_[npt.kind].q * n_in_cell; inject.. end // negative charges For each population pop if (pop == Neutralizing population) Continue call init_npt (pop, …) If npt.kind.q >= 0 continue // else negative charge n_in_cell = get_n_in_cell // make sure don’t add non-neutralizable amount if n_q_in_cell + kinds_[npt.kind].q * n_in_cell < 0 n_in_cell = -n_q_in_cell / kinds_[npt.kind].q; n_q_in_cell += kinds_[npt.kind].q * n_in_cell; inject.. end // now add neutralizing population call init_npt (neutralizing population, …) n_in_cell = -n_q_in_cell / kinds_[npt.kind].q; assert (n_in_cell >= 0) !!!! inject... Thanks! Will

…

On Sep 8, 2020, at 9:40 PM, jcd496 ***@***.***> wrote: I see that his build failed, idk though. It builds fine on summit with spack — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[jcd496]

I don't expect this to be merged, but here it is for reference.

I tried using the following logic as an alternative

        if(kind == HE_population_){
          target_.init_npt(kind, pos, npt);
          npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]);
        }else{
          target_.init_npt(kind, pos, npt);
          npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]);
          if(kind != neutralizing_population_)  // ie kind == electron and thus scale off of HE_population
            npt.n *= (1 - HE_ratio_);
        }

But the ratio of electron_he/electron did not end up working out to the correct value. I don't know why, maybe someone can see the error?

[wrfox]

Hi John, Thanks for sending. I couldn’t quite follow this code - where does this code live?

if(kind == HE_population_){ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); }else{ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); if(kind != neutralizing_population_) // ie kind == electron and thus scale off of HE_population npt.n *= (1 - HE_ratio_); } But the ratio of electron_he/electron did not end up working out to the correct value. I don't know why, maybe someone can see the error?

Anyway, thanks for your efforts and more comments below Will

…

--- One question - how is “kind” getting set and are you sure it is always correct? Here is my (again, pseudo-code) for how to inject. (I couldn’t quickly find the population numbering / enum, so these are just my randomly chosen numbers) Also, since this code is for the particle “source term”, the code should ideally be part of the “case’ — psc_flat_foil — it is not generic code. (In Fortran, we had CASE_nVT, which returned the information for (1) the initial condition at timestep=0, then (2) the particle source term (particles to add), for t>0) // Calculate local electron density from moments ne_tot = ne + ne_he // Density to add n_add = (2.5) - ne_tot. if pop == 1 // ion n_add = n_add kind = ion // rest of parameters, Ti, etc if pop == 2: // HE electrons HE density to add = n_add * HE_factor kind = HE; // rest of parameters if pop == 3 // regular electron: n_add = n_add * (1-HE_factor) kind = regular electron // rest of parameters TESTS - I suggest the following tests: HE_ratio = 0 — no HE particles added HE_ratio = 0.01 — 1% of all electrons in the simulation are HE HE_ratio = 0.5 — 50% of all electrons in the simulation are HE — may be easier to test than 1% due to the fractional particle HE_ratio = 1 — all of the electrons in the simulation are HE

On Sep 9, 2020, at 5:19 PM, jcd496 ***@***.***> wrote: I don't expect this to be merged, but here it is for reference. I tried using the following logic as an alternative if(kind == HE_population_){ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); }else{ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); if(kind != neutralizing_population_) // ie kind == electron and thus scale off of HE_population npt.n *= (1 - HE_ratio_); } But the ratio of electron_he/electron did not end up working out to the correct value. I don't know why, maybe someone can see the error? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

[jcd496]

This part was my idea for implementing without needing to include the non_neutralizing_population. It seemed logically consistent with the commited implementation but for whatever reason didn't work correctly.

Hi John, Thanks for sending. I couldn’t quite follow this code - where does this code live?
if(kind == HE_population_){ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); }else{ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); if(kind != neutralizing_population_) // ie kind == electron and thus scale off of HE_population npt.n *= (1 - HE_ratio_); } But the ratio of electron_he/electron did not end up working out to the correct value. I don't know why, maybe someone can see the error?
Anyway, thanks for your efforts and more comments below Will

Re below: I agree that the inject implementation is case specific. Perhaps the easiest method would be to move the init_npt lambda function to psc_flatfoil_yz.cxx and pass that to the Inject constructor. in that case we would need to maintain the baseline lambda function in the master branch and keep our specific case in its own working branch. kind gets set it setup_particles.hxx, injection is wrapped in for loop iterating through each kind

I can get those test cases for tomorrows meeting

--- One question - how is “kind” getting set and are you sure it is always correct? Here is my (again, pseudo-code) for how to inject. (I couldn’t quickly find the population numbering / enum, so these are just my randomly chosen numbers) Also, since this code is for the particle “source term”, the code should ideally be part of the “case’ — psc_flat_foil — it is not generic code. (In Fortran, we had CASE_nVT, which returned the information for (1) the initial condition at timestep=0, then (2) the particle source term (particles to add), for t>0) // Calculate local electron density from moments ne_tot = ne + ne_he // Density to add n_add = (2.5) - ne_tot. if pop == 1 // ion n_add = n_add kind = ion // rest of parameters, Ti, etc if pop == 2: // HE electrons HE density to add = n_add * HE_factor kind = HE; // rest of parameters if pop == 3 // regular electron: n_add = n_add * (1-HE_factor) kind = regular electron // rest of parameters TESTS - I suggest the following tests: HE_ratio = 0 — no HE particles added HE_ratio = 0.01 — 1% of all electrons in the simulation are HE HE_ratio = 0.5 — 50% of all electrons in the simulation are HE — may be easier to test than 1% due to the fractional particle HE_ratio = 1 — all of the electrons in the simulation are HE
On Sep 9, 2020, at 5:19 PM, jcd496 @.***> wrote: I don't expect this to be merged, but here it is for reference. I tried using the following logic as an alternative if(kind == HE_population_){ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); }else{ target_.init_npt(kind, pos, npt); npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]); if(kind != neutralizing_population_) // ie kind == electron and thus scale off of HE_population npt.n *= (1 - HE_ratio_); } But the ratio of electron_he/electron did not end up working out to the correct value. I don't know why, maybe someone can see the error? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.