first of all, thank you very much for the kind support and for making Raysect rock!
Coming more specifically to my issue, I have been using CHERAB for a number of applications consisting in the evaluation of the radiation heat load over in-tokamak plasma-facing components, given a certain source of radiation in the plasma. Interest is now mainly directed in exploiting the SOLPS-ITER-computed radiation emission map as a source of radiation in CHERAB. To do so, Raysect's classes in inhomogeneous.pyx are employed.
Attached to the present, you may find an example of a standard-looking SOLPS-ITER-computed radiation emission distribution in the poloidal plane (fig. 1). As you can see, the vast very majority of the emission (log10 scale!) and of its spatial gradients (!) occur in the divertor region. Viceversa, in the core (which, in terms of volume, is the most important actor by far) there is none at all. Even if there is some non-zero core radiation (e.g. computed by ASTRA, see fig. 2 attached), the core plasma physics is such that to give smoothly varying profiles, markedly in contrast to the plasma behaviour in the divertor. Therefore, one is forced to use a small integration step along the rays, in such a way to capture the details (the maximum in particular) of the dominant divertor emission. In doing so, however, one ends up with an unnecessary high number of samples in the core.
Here below (_comput_next_point function), you may find a possible idea to select a proper integration step as a function of the local features of the radiation emission itself. The idea is to set a maximum step (e.g. 1e-3 m) which is then decreased where needed, according to the 1st and 2nd order derivatives in space of the radiation emission function. In this way, while the core could be sampled every 1e-3 m, the step in the divertor will end up being smaller (or much smaller). As I am extremely far from being able to implement this idea in inhomogeneous.pyx itself :'), I ran some homemade tests by shooting one ray at a time, to compare the baseline "uniform" and this "differential" (this bad and somewhat meaningless name is very up to discussion) approach. In a 2D case, the "differential" approach has proved to be faster than "uniform", while retaining a satisfactory accuracy. In a full 3D instance, this improvement should be even more noticeable (the weight of the core region will increase when moving from 2D to 3D).
However, I suspect that the main advantage of the "differential" strategy would lie in its higher user-friendliness: according to my experience, the combination uniform integration step and distracted user (myself in the first place) may result in an inaccurate sampling which, e.g., could miss the tiny volume where the maximum emission occurs. To avoid this, the well-educated user should perform a sampling-step-independence study, consisting in running different simulations with a gradually smaller integration step (and so gradually more computationally expensive), to hopefully show that the result is unaffected by the actual choice of the step. A "differential" (damn this name, it sounds really bad) approach would presumably take some responsibility away from the user in this sense.
#####################
if sampling == "uniform":
return p1 + step * direction
#####################
elif sampling == "differential":
p2 = p1 + step * direction
f1 = rad_function(p1.x, p1.y)
f2 = rad_function(p2.x, p2.y)
# both exactly zero
# => just skip to next sampling point
# (also because will give /0 and NaN below)
if f2 == f1:
return p2
# compute local 1st and 2nd order derivatives:
# - centred finite difference scheme (higher precision)
# - extrapolation of rad_function must be True
p0 = p1 - 0.5 * step * direction
p2 = p1 + 0.5 * step * direction
f0 = rad_function(p0.x, p0.y)
f2 = rad_function(p2.x, p2.y)
dfdx = (f2 - f0) / step
d2fdx2 = (f2 - 2 * f1 + f0) / (step ** 2)
# characteristic length of 1st and 2nd order variation
L1st = np.abs(f1 / dfdx)
L2nd = np.abs(dfdx / d2fdx2)
# most stringent requirement among the three lengths
# and no longer than step
step = np.min([L1st, L2nd, step])
return p1 + step * direction
else:
raise ValueError("Sampling strategy not recognised! Must be either uniform or differential.")
Hi all,
first of all, thank you very much for the kind support and for making Raysect rock!
Coming more specifically to my issue, I have been using CHERAB for a number of applications consisting in the evaluation of the radiation heat load over in-tokamak plasma-facing components, given a certain source of radiation in the plasma. Interest is now mainly directed in exploiting the SOLPS-ITER-computed radiation emission map as a source of radiation in CHERAB. To do so, Raysect's classes in inhomogeneous.pyx are employed.
Attached to the present, you may find an example of a standard-looking SOLPS-ITER-computed radiation emission distribution in the poloidal plane (fig. 1). As you can see, the vast very majority of the emission (log10 scale!) and of its spatial gradients (!) occur in the divertor region. Viceversa, in the core (which, in terms of volume, is the most important actor by far) there is none at all. Even if there is some non-zero core radiation (e.g. computed by ASTRA, see fig. 2 attached), the core plasma physics is such that to give smoothly varying profiles, markedly in contrast to the plasma behaviour in the divertor. Therefore, one is forced to use a small integration step along the rays, in such a way to capture the details (the maximum in particular) of the dominant divertor emission. In doing so, however, one ends up with an unnecessary high number of samples in the core.
Here below (_comput_next_point function), you may find a possible idea to select a proper integration step as a function of the local features of the radiation emission itself. The idea is to set a maximum step (e.g. 1e-3 m) which is then decreased where needed, according to the 1st and 2nd order derivatives in space of the radiation emission function. In this way, while the core could be sampled every 1e-3 m, the step in the divertor will end up being smaller (or much smaller). As I am extremely far from being able to implement this idea in inhomogeneous.pyx itself :'), I ran some homemade tests by shooting one ray at a time, to compare the baseline "uniform" and this "differential" (this bad and somewhat meaningless name is very up to discussion) approach. In a 2D case, the "differential" approach has proved to be faster than "uniform", while retaining a satisfactory accuracy. In a full 3D instance, this improvement should be even more noticeable (the weight of the core region will increase when moving from 2D to 3D).
However, I suspect that the main advantage of the "differential" strategy would lie in its higher user-friendliness: according to my experience, the combination uniform integration step and distracted user (myself in the first place) may result in an inaccurate sampling which, e.g., could miss the tiny volume where the maximum emission occurs. To avoid this, the well-educated user should perform a sampling-step-independence study, consisting in running different simulations with a gradually smaller integration step (and so gradually more computationally expensive), to hopefully show that the result is unaffected by the actual choice of the step. A "differential" (damn this name, it sounds really bad) approach would presumably take some responsibility away from the user in this sense.
What do you think about this?
Cheers, Matteo
##########################################################################################
def _compute_next_point(sampling = "uniform", step = 1.0E-02, direction = None, p1 = None, rad_function = None):