No n_phi dependence in afsi_get_volume

[BCHamstra]

While trying to replicate older stellarator work in ASCOT5 (https://doi.org/10.1088/1361-6587/abd981), I stumbled on an issue in afsi.c. The afsi_get_volume method has no dependence on the number of phi bins and instead assumes the volume goes from 0 to 2pi. I believe the 2pi should be replaced with a dphi = (dist->dist_5D->max_phi - dist->dist_5D->min_phi) / dist->dist_5D->n_phi; and similarly for the thermal distribution, if I understand correctly. Though this unfortunately does not resolve the issues I am currently having.

On a somewhat related note, I do not know whether it is an issue, but I feel like I should mention it. Slightly above in the afsi_sample_thermal method I noticed the line E = -temp * (log(r1) - log(r2) * cos(CONST_PI/2)*cos(CONST_PI/2) * r3);. Aren't the later terms evaluated to 0...?

[miekkasarki]

Huge thanks for spotting that latter point! I've been wondering why AFSI gives me 337 MW of alpha particle power in one case where the power is supposed to be 330 MW. After fixing that part (the brackets were all over the place) I finally get the correct answer with AFSI.

As for the phi bins issue, what you say sounds reasonable to me. However, we are still lacking a stellarator test case and I don't have any stellarator data at hand at the moment. Would it be okay for you to make a test run with AFSI after I've fixed these issues?

Also note that there was another bug #109 in AFSI related to phi angle (the fix is still in the develop branch).

[BCHamstra]

I did see that bug, so I was already using the development branch just forgot to mention it. I will have another look at it next week, though I tried fixing it myself already and I suspect I am facing another issue as well. I will update if I find anything.

[BCHamstra]

So it would appear my other issues were mostly related to simulation settings. I changed the volume method and now seem to be getting similar results to the mentioned paper:

real afsi_get_volume(afsi_data* dist, int iR) {
    real dR, dphi, dz;

    if(dist->type == 1) {
        dR   = (dist->dist_5D->max_r - dist->dist_5D->min_r) / dist->dist_5D->n_r;
        dphi = (dist->dist_5D->max_phi - dist->dist_5D->min_phi) / dist->dist_5D->n_phi;
        dz   = (dist->dist_5D->max_z - dist->dist_5D->min_z) / dist->dist_5D->n_z;
        return dphi*dz*dR*(dist->dist_5D->min_r + dR*(iR + 0.5));
    }

    else if(dist->type == 2) {
        dR   = (dist->dist_thermal->max_r - dist->dist_thermal->min_r)
            / dist->dist_thermal->n_r;
        dphi = (dist->dist_thermal->max_phi - dist->dist_thermal->min_phi)
            / dist->dist_thermal->n_phi;
        dz   = (dist->dist_thermal->max_z - dist->dist_thermal->min_z)
            / dist->dist_thermal->n_z;
        return dphi*(dist->dist_thermal->min_r + iR * dR + 0.5*dR)*dR*dz;
    }
    return 0;
}

So I think this works.