Plotting example

[emaartensson]

Hello,

I am unsure if this is an issue or if, most likely, I am doing something wrong. When I try to plot the example plasma parameters (from EU-DEMO I believe), I get very strange values for the strength and neutron density in the core of the plasma. The only way I can affect it is by increasing the core ion density. The only thing that looks reasonable is the ion temperature.

I've tried changing the radius units to centimeter but same results. Unsure what I am doing wrong since I am just copying the parameters from the example.

Any help would be appreciated!

from openmc_plasma_source import TokamakSource, plotting 
import matplotlib.pyplot as plt

source = TokamakSource(
    elongation=1.557,
    ion_density_centre=1.09e20,
    ion_density_peaking_factor=1,
    ion_density_pedestal=1.09e20,
    ion_density_separatrix=3e19,
    ion_temperature_centre=45.9,
    ion_temperature_peaking_factor=8.06,
    ion_temperature_pedestal=6.09,
    ion_temperature_separatrix=0.1,
    major_radius=9.06,
    minor_radius=2.92258,
    pedestal_radius=0.8 * 2.92258,
    mode="H",
    shafranov_factor=0.44789,
    triangularity=0.270,
    ion_temperature_beta=6) 

fig = plotting.scatter_tokamak_source(source, quantity='ion_temperature')
plt.title('Ion temperature [keV]')
plt.savefig('ion_temp.png')
plt.close()

fig = plotting.scatter_tokamak_source(source, quantity='strength')
plt.title('Source strength []')
plt.savefig('strength.png')
plt.close()

fig = plotting.scatter_tokamak_source(source, quantity='neutron_source_density')
plt.title('Neutron source density [neutrons/m3/s]')
plt.savefig('neutron_density.png')
plt.close()

[RemDelaporteMathurin]

Hi @emaartensson thanks for reporting this. It is indeed a weird behaviour...

I don't have a lot of bandwith to investigate this right now, but will have a look soon.

In the meantime, could you have a look at the ion density too?

Since this is the method responsible for the neutron source density, I expect that the problem lies with the ion density: https://github.com/fusion-energy/openmc-plasma-source/blob/6be4bc815ee788a8aabad7e240054df3bb35f683/openmc_plasma_source/tokamak_source.py#L288-L303

Then, this is the method responsible for the ion density calculation: https://github.com/fusion-energy/openmc-plasma-source/blob/6be4bc815ee788a8aabad7e240054df3bb35f683/openmc_plasma_source/tokamak_source.py#L125-L159

But since you are copying the parameters from the example, this is a bit weird indeed.

[RemDelaporteMathurin]

Well, found a bit of time to at least run the example again. I confirm that I can reproduce the bug with the latest version.

Running this script with v0.2 produces the correct behaviour:

from openmc_plasma_source import TokamakSource

my_plasma = TokamakSource(
    elongation=1.557,
    ion_density_centre=1.09e20,
    ion_density_peaking_factor=1,
    ion_density_pedestal=1.09e20,
    ion_density_separatrix=3e19,
    ion_temperature_centre=45.9,
    ion_temperature_peaking_factor=8.06,
    ion_temperature_pedestal=6.09,
    ion_temperature_separatrix=0.1,
    major_radius=9.06,
    minor_radius=2.92258,
    pedestal_radius=0.8 * 2.92258,
    mode="H",
    shafranov_factor=0.44789,
    triangularity=0.270,
    ion_temperature_beta=6,
)

my_plasma.sample_sources()

import matplotlib.pyplot as plt

plt.scatter(my_plasma.RZ[0], my_plasma.RZ[1], c=my_plasma.strengths)
plt.gca().set_aspect("equal")
plt.show()

[RemDelaporteMathurin]

Right. Found the bug. Somewhere between 0.2 and the current version, an error found its way the expression for density of ions. This: https://github.com/fusion-energy/openmc-plasma-source/blob/6be4bc815ee788a8aabad7e240054df3bb35f683/openmc_plasma_source/tokamak_source.py#L146-L158

should be:

            density = np.where(
                r < self.pedestal_radius,
                (
                    (self.ion_density_centre - self.ion_density_pedestal)
                    * (1 - (r / self.pedestal_radius) ** 2)
                    ** self.ion_density_peaking_factor
                    + self.ion_density_pedestal
                ),
                (
                    (self.ion_density_pedestal - self.ion_density_separatrix)
                    * (self.major_radius - r)
                    / (self.major_radius - self.pedestal_radius)
                    + self.ion_density_separatrix
                ),
            )

I'll push a fix realy quick.

Pinging @shimwell too

[RemDelaporteMathurin]

This is fixed in the latest release 0.2.8. Thanks @emaartensson for reporting the bug!

[shimwell]

Thanks for fixing it Remi and getting a release out so quickly. super fast response rate all round.

[emaartensson]

Great stuff Remi, thank you so much for fixing it so quickly!