A scaling factor is needed to get correct plasma current

[leferi99]

When trying to prescribe the electric field and set the plasma current, I noticed that I have to use a scaling factor to get the plasma current I need. For example I used to do it this way:

conductivity = do.other.fluid.conductivity.getData()
jj = (Ip0 * I_p_corr) / np.trapz(2 * np.pi * radialgrid * jprof, radialgrid)
efield=jj*jprof/conductivity[-1,:]

Here Ip0 contains the correct total plasma current, jprof is a normalized current density profile and I_p_corr is the needed scaling factor. In the case of my EU-DEMO simulations with 17.75 MA total plasma current, I used ~0.68 for this factor so that I would get the 17.75 MA from the DREAMOutput.

[hoppe93]

Thanks for reporting Feri! I've seen indications of this previously as well, but have not had time to dig deeper :P Normally, the test-particle collision operator of DREAM will underestimate the ohmic current by (very roughly) 0.5. Through simulations, however, it was found that the conductivity of kinetic simulations with the test-particle collision operator are related to the full Braams-Karney conductivity via the relatively simple expression

$$\sigma{\rm kineq} = \left(1-\frac{1.406}{1.888+Z{\rm eff}}\right)\sigma_0,$$

(see section 2.2.1 of the DREAM paper), which is therefore used to compensate for the missing current in the code when ds.eqsys.j_ohm.setCorrectedConductivity(JOhm.CORRECTED_CONDUCTIVITY_ENABLED) is used.

Despite all this, it seems that the corrected conductivity is not working as it should, since the total ohmic current now seems to be underestimated (i.e. for a given electric field, we compensate by too much). Could you perhaps try to disable the corrected conductivity mode and see what correction factor you would need?

[leferi99]

I managed to find a workaround for kinetic simulations. I say it's a workaround but it's probably an intentional thing. I used the following code:

import numpy as np
from DREAM.DREAMSettings import DREAMSettings
import DREAM.Settings.Equations.OhmicCurrent as jOhm

radialgrid=np.linspace(0, a, Nr)

ds = DREAMSettings()
ds.eqsys.E_field.setType(Efield.TYPE_PRESCRIBED_OHMIC_CURRENT)
ds.eqsys.E_field.setBoundaryCondition()

ds.eqsys.j_ohm.setCurrentProfile(jprof, radius=radialgrid, Ip0=Ip0)
ds.eqsys.j_ohm.setConductivityMode(jOhm.CONDUCTIVITY_MODE_SAUTER_COLLISIONAL)
ds.eqsys.j_ohm.setCorrectedConductivity(jOhm.CORRECTED_CONDUCTIVITY_DISABLED)

Here jprof is a normalized current density profile, Ip0 is the correct total plasma current. This code correctly sets up the radial current profile and the total plasma current.

I will still check the other method with the disabled conductivity tomorrow and get back to this thread. What I can say presently regarding that issue is that WITH the corrected conductivity, I needed ~0.68 correction factor.

[leferi99]

Okay so the simulation finished with DISABLED corrected conductivity and the correction factor I needed this way was 1.15945558

[hoppe93]

Interesting! What was the $Z_{\rm eff}$ in that simulation (and by extension, the conductivity correction factor, as written above)?

[leferi99]

The effective charge was around 2.12 for most of the grid, therefore the correction factor calculated from above should be 0.6492. This number is interestingly close to the 0.68 I needed for the corrected conductivity case...

But also keep in mind that I didn't properly take into account the D-shaped geometry of the plasma (which was given with numerical profiles) when calculating the jj. Also it's entirely possible that something else in my code is affecting these numbers.

[hoppe93]

Thanks! But that must be a coincidence, right? Because the factor 0.68 is what you need when the corrected conductivity is used, while 0.6492 is the actual correction factor in that case (which gives a too high current density). Also the shaping should not matter. The desired behavior is that j = sigma*E is the same both in fluid and kinetic simulations (regardless of geometry), and since you use the "fluid" sigma to estimate E, you should get the correct current density in the kinetic simulation.

This is very strange and concerning... I will have to take a closer look (although it might take a while before I have time)

[leferi99]

I would also think it's just a coincidence if some very banal swap didn't happen in the source code regarding the corrected and the not corrected conductivity modes, and I'm not suggesting it did. As far as I can see looking at the OhmicCurrent.py file's commit history this should indeed not be the case, but I'm not that familiar with the inner C++ code.

In the meantime that workaround I mentioned should work fine right?

[hoppe93]

I am pretty sure the corrected conductivity mode has worked correctly at some point in the past. Something must have happened in the C++ code which has made it no longer work...

For the workaround you suggested, does that work in the kinetic mode (i.e. with "FULLY_KINETIC" collisions)? It looks to me as you should get the correct ohmic current in fluid/isotropic mode, but the wrong ohmic current in fully kinetic mode. For the fully kinetic mode, your original workaround should probably be better (i.e. multiplying the electric field with a correcting factor), although it will still be flawed. Since in practice, we now seem to have an incorrect conductivity, the current decay and diffusion will be different... The only way to resolve that is to dig into the C++ code unfortunately :/

[leferi99]

Sorry I got confused, my workaround works in the isotropic hottail case, not in the fully kinetic hottail case so that's not that relevant to this topic you're right.