Skip to content

Commit

Permalink
attempted fix for negative jbar errors
Browse files Browse the repository at this point in the history
  • Loading branch information
jhmatthews committed Sep 24, 2018
1 parent 826e3b4 commit dcae954
Show file tree
Hide file tree
Showing 3 changed files with 21 additions and 13 deletions.
3 changes: 1 addition & 2 deletions source/matom.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -501,7 +501,6 @@ matom (p, nres, escape)


/* Co-moving frequency - changed to rest frequency by doppler */
/*Currently this assumed hydrogenic shape cross-section - Improve */
}
else
{ //collisional deactivation
Expand Down Expand Up @@ -1005,7 +1004,7 @@ kpkt (p, nres, escape)
/* Now (as in matom) choose a frequency for the new packet. */

//p->freq = phot_top[i].freq[0] - (log (1. - random_number(0.0,1.0)) * xplasma->t_e / H_OVER_K);
p->freq = matom_select_bf_freq(one, i);
p->freq = matom_select_bf_freq(one, i);

/* if the cross-section corresponds to a simple ion (macro_info == 0)
or if we are treating all ions as simple, then adopt the total emissivity
Expand Down
23 changes: 16 additions & 7 deletions source/resonate.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -994,7 +994,7 @@ scatter (p, nres, nnscat)
int i, n;
double p_init[3], p_final[3], dp[3], dp_cyl[3];
WindPtr one;
double prob_kpkt, kpkt_choice;
double prob_kpkt, kpkt_choice, freq_comoving;
double gamma_twiddle, gamma_twiddle_e, stim_fact;
int m, llvl, ulvl;
double v_dop;
Expand All @@ -1004,9 +1004,14 @@ scatter (p, nres, nnscat)




stuff_phot (p, &pold);
n = where_in_grid (wmain[pold.grid].ndom, pold.x); // Find out where we are

vwind_xyz (ndom, p, v); //get the local velocity at the location of the photon
v_dop = dot (p->lmn, v); //get the dot product of the photon direction with the wind, to get the doppler velocity
freq_comoving = p->freq * (1. - v_dop / C); //This is the photon frequency in the comoving frame

if (n < 0)
{
Error ("scatter: Trying to scatter a photon in grid cell %d\n", n);
Expand Down Expand Up @@ -1157,10 +1162,16 @@ scatter (p, nres, nnscat)
that goes into the electron rather than being stored as ionisation energy: this
fraction gives the selection probability for the packet. It's given by the
(photon frequency / edge frequency - 1) (SS) */

prob_kpkt = 1. - (phot_top[*nres - NLINES - 1].freq[0] / freq_comoving);


prob_kpkt =
1. - (phot_top[*nres - NLINES - 1].freq[0] / p->freq);
if (prob_kpkt < 0)
{
Error("scatter: kpkt probability (%8.4e) < 0, zeroing\n", prob_kpkt);
Log("scatter: photon comoving frequency: %8.4e, edge frequency %8.4e\n",
phot_top[*nres - NLINES - 1].freq[0], freq_comoving);
prob_kpkt = 0.0;
}


/* Now choose whether or not to make a k-packet. */
Expand Down Expand Up @@ -1226,9 +1237,7 @@ scatter (p, nres, nnscat)

if (*nres == -1) //Its an electron scatter, so we will call compton to get a direction
{
vwind_xyz (ndom, p, v); //get the local velocity at the location of the photon
v_dop = dot (p->lmn, v); //get the dot product of the photon direction with the wind, to get the doppler velocity
p->freq = p->freq * (1. - v_dop / C); //This is the photon frequency in the electron rest frame
p->freq = freq_comoving; //This is the photon frequency in the electron rest frame calculated earlier in the routine
compton_dir (p, xplasma); //Get a new direction using the KN formula
v_dop = dot (p->lmn, v); //Find the dot product of the new velocity with the wind
p->freq = p->freq / (1. - v_dop / C); //Transform back to the observers frame
Expand Down
8 changes: 4 additions & 4 deletions source/wind_updates2d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -333,7 +333,7 @@ WindPtr (w);
MPI_Pack (&plasmamain[n].kappa_ff_factor, 1, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n].nscat_es, 1, MPI_INT, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (plasmamain[n].recomb_simple, nphot_total, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (plasmamain[n].recomb_simple_upweight, nphot_total, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (plasmamain[n].recomb_simple_upweight, nphot_total, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n].kpkt_emiss, 1, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n].kpkt_abs, 1, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (plasmamain[n].kbf_use, nphot_total, MPI_INT, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
Expand Down Expand Up @@ -472,7 +472,7 @@ WindPtr (w);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].kappa_ff_factor, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].nscat_es, 1, MPI_INT, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, plasmamain[n].recomb_simple, nphot_total, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, plasmamain[n].recomb_simple_upweight, nphot_total, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, plasmamain[n].recomb_simple_upweight, nphot_total, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].kpkt_emiss, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].kpkt_abs, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, plasmamain[n].kbf_use, nphot_total, MPI_INT, MPI_COMM_WORLD);
Expand Down Expand Up @@ -1148,12 +1148,12 @@ wind_rad_init ()
if (geo.macro_simple==0 && phot_top[i].macro_info==1)
{
plasmamain[n].recomb_simple[i] = 0.0;
plasmamain[n].recomb_simple_upweight[i] = 1.0;
plasmamain[n].recomb_simple_upweight[i] = 1.0;
}
else
{ //we want a macro approach, but not for this ion so need recomb_simple
plasmamain[n].recomb_simple[i] = alpha_store = alpha_sp (&phot_top[i], &plasmamain[n], 2);
plasmamain[n].recomb_simple_upweight[i] = alpha_sp (&phot_top[i], &plasmamain[n], 1) / alpha_store;
plasmamain[n].recomb_simple_upweight[i] = alpha_sp (&phot_top[i], &plasmamain[n], 1) / alpha_store;
}
}

Expand Down

0 comments on commit dcae954

Please sign in to comment.