Skip to content

Added non-conservative Krook operator to python wrapper and JETPLUME #41

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
kgklein merged 1 commit into from
Mar 23, 2026
Merged

Added non-conservative Krook operator to python wrapper and JETPLUME #41

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7,463 changes: 7,376 additions & 87 deletions examplelinfpc.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

138 changes: 104 additions & 34 deletions linfpclib/linfpc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,21 @@ def find_nearest(array, value): #random but very useful function
idx = (np.abs(array - value)).argmin()
return idx

import math


class ParamsDict(dict):
def __init__(self, owner, *args, **kwargs):
super().__init__(*args, **kwargs)
self.owner = owner

def __setitem__(self, key, value):
old_value = self.get(key, None)
super().__setitem__(key, value)

if key == 'Kn' and old_value != value: #works like a setter
self.owner.recompute_all_nuns()

class plume_input:
#class that have dictss with all inputs for each namelist
#params
Expand All @@ -37,11 +52,28 @@ class plume_input:
def __init__(self,dataname):
self.dataname = dataname

# make these instance-level, not shared across all objects
self.params = ParamsDict(self)
self.fpc = {}
self.species = [{}]
self.maps = {}
self.scan_inputs = [{}]
self.guesses = [{}]

def read_input(): #typically load from sample
pass

def recompute_all_nuns(self):
if len(self.species) == 0:
return
if len(self.species[0]) == 0:
return

for i in range(len(self.species)):
self.species[i]['nu_ns'] = self.get_nuns(self.species[i])

def set_params(self,betap,kperp,kpar,vtp,nspec,nscan,option,\
nroot_max,use_map,writeOut):
nroot_max,use_map,writeOut,collision_type=0,Kn=None):
use_map_out = '.true.'
if(not(use_map)):
use_map_out = '.false.'
Expand All @@ -50,7 +82,8 @@ def set_params(self,betap,kperp,kpar,vtp,nspec,nscan,option,\
if(not(use_map)):
usewriteOut_out = '.false.'

self.params = {'betap':betap,
self.params = ParamsDict(self, {
'betap':betap,
'kperp':kperp,
'kpar':kpar,
'vtp':vtp,
Expand All @@ -59,8 +92,10 @@ def set_params(self,betap,kperp,kpar,vtp,nspec,nscan,option,\
'option':int(option),
'nroot_max':int(nroot_max),
'use_map':use_map_out,
'writeOut':usewriteOut_out
}
'writeOut':usewriteOut_out,
'collision_type':int(collision_type),
'Kn':Kn
})

def set_fpc(self,vperpmin=0,vperpmax=0,vparmin=0,vparmax=0,delv=0,vxmin=None,vxmax=None,vymin=None,vymax=None,vzmin=None,vzmax=None,elecdircontribution=0.):
self.fpc = {'vperpmin':vperpmin,
Expand All @@ -73,10 +108,10 @@ def set_fpc(self,vperpmin=0,vperpmax=0,vparmin=0,vparmax=0,delv=0,vxmin=None,vxm
if(vxmin!=None and vxmax!=None and vymin!=None and vymax!=None and vzmin!=None and vzmax!=None):
self.fpc['vxmin'] = vxmin
self.fpc['vxmax'] = vxmax
self.fpc['vymin'] = vymin
self.fpc['vymax'] = vymax
self.fpc['vzmin'] = vzmin
self.fpc['vzmax'] = vzmax
self.fpc['vymin'] = vymin
self.fpc['vymax'] = vymax
self.fpc['vzmin'] = vzmin
self.fpc['vzmax'] = vzmax

def set_maps(self,loggridw,omi,omf,gami,gamf,positive_roots):
loggridw_out = '.true.'
Expand All @@ -95,31 +130,61 @@ def set_maps(self,loggridw,omi,omf,gami,gamf,positive_roots):
'positive_roots':positive_roots_out
}

def get_nuns(self, spec):
muS = spec['muS']
tauS = spec['tauS']

# Fortran spec(1)%alph_s -> Python self.species[0]['alphS']
# If no species exist yet, assume this spec IS the first/reference species
if len(self.species) == 0:
alpha1 = spec['alphS']
else:
alpha1 = self.species[0]['alphS']

try:
Kn = self.params['Kn']
except:
print("Error with Kn = self.params['Kn'], returning 0 for nu_ns...")
return 0.

if(Kn is None):
return None

nu_ns = 1.0 / (math.sqrt(2.0) * Kn) * math.sqrt(muS / (tauS * alpha1))

return nu_ns

def make_species(self,tauS, muS, alphS, Qs, Ds, vvS,spec_n=-1):
tempspecies = {'tauS':tauS,
'muS':muS,
'alphS':alphS,
'Qs':Qs,
'Ds':Ds,
'vvS':vvS
'vvS':vvS,
}
if(len(self.species[0]) == 0):
if(spec_n == -1 or spec_n == 1):
print("No species found, creating first species...")
self.species = [tempspecies]
idx = 0
else:
print("Warning: spec_n ==",spec_n,"but there are no species here...")
else:
num_spec = len(self.species)
if(spec_n == -1 or num_spec+1 == spec_n):
print("Appending species to list. Total species is now ",num_spec+1)
self.species.append(tempspecies)
idx = -1
elif(num_spec+1 < spec_n):
print("Warning: there are only ",num_spec,"species but user requested we create spec number",spec_n)
print("Please input a valid spec_n...")
print("Please input a valid spec_n... Returning")
return
else:
print("Replacing species number",spec_n)
self.species[spec_n-1] = tempspecies
idx = spec_n - 1

self.species[idx]['nu_ns'] = self.get_nuns(self.species[idx])

def make_scan(self,scan_type,scan_style,swi,swf,swlog,ns,nres,heating,eigen):
swlog_out = '.true.'
Expand Down Expand Up @@ -186,8 +251,16 @@ def write_input(self,flnm,outputname,desc='',verbose=False):
line = str(key)+'='+"'"+str(outputname)+"'"+'\n'
else:
line = str(key)+'='+str(self.params[key])+'\n'
if(key == 'Kn'):
if(not(self.params['Kn']) is None):
line = str(key)+'='+str(self.params[key])+'\n'
else:
continue
f.write(line)

if(self.params['Kn'] != None and self.params['collision_type'] == 0):
print("Warning, Kn has a finite value but collision_type is 0 (i.e. collisions are off). Was this intentional?")

line = 'dataName'+"='"+str(self.dataname)+"'\n"
f.write(line)
f.write('/\n\n')
Expand All @@ -203,6 +276,8 @@ def write_input(self,flnm,outputname,desc='',verbose=False):
for spec in self.species:
f.write('&species_'+str(specidx)+'\n')
for key in spec.keys():
if(key == 'nu_ns'):
continue
line = str(key)+'='+str(spec[key])+'\n'
f.write(line)
f.write('/\n\n')
Expand Down Expand Up @@ -238,7 +313,7 @@ def write_input(self,flnm,outputname,desc='',verbose=False):
f.close()

def load_from_file(self,flnm):
paramkeys = ['betap','kperp','kpar','vtp','nspec','nscan','option','nroot_max','use_map','writeOut','outputName'] #note dataName is missing as that is set by self.dataname and defined at creation
paramkeys = ['betap','kperp','kpar','vtp','nspec','nscan','option','nroot_max','use_map','writeOut','outputName','collision_type','Kn'] #note dataName is missing as that is set by self.dataname and defined at creation
fpckeys = ['vperpmin','vperpmax','vparmin','vparmax','delv','elecdircontribution']
specieskeys = ['tauS','muS','alphS','Qs','Ds','vvS']
mapskeys = ['loggridw','omi','omf','gami','gamf','positive_roots']
Expand Down Expand Up @@ -276,7 +351,7 @@ def load_from_file(self,flnm):
tempparamdict[parse[0]] = parse[1].split('!')[0].replace('\n','')

_i += 1
self.params = tempparamdict
self.params = ParamsDict(self, tempparamdict)


if(parse[0] == '&fpc'):
Expand Down Expand Up @@ -371,15 +446,6 @@ def load_from_file(self,flnm):
_i += 1


#main dict
namelists = {}
params = {}
fpc = {}
species = [{}]
maps = {}
scan_inputs = [{}]
guesses = [{}]

dataname = 'default'

def _replace_input_aux(inputflnm,verbose=False):
Expand Down Expand Up @@ -442,19 +508,23 @@ class containing input parameters
if(verbose): print("Reading roots from ",rootflnm)

roots = []
tempfile = open(rootflnm, "r")
for line in tempfile:
parse = line.split()
try:
temp_om = float(parse[4])
except:
temp_om = float(1*10**99.)
try:
temp_gam = float(parse[5])
except:
temp_gam = float(1*10**99.)
roots.append(temp_om+temp_gam*1j)
tempfile.close()
try:
tempfile = open(rootflnm, "r")
for line in tempfile:
parse = line.split()
try:
temp_om = float(parse[4])
except:
temp_om = float(1*10**99.)
try:
temp_gam = float(parse[5])
except:
temp_gam = float(1*10**99.)
roots.append(temp_om+temp_gam*1j)
except Exception as e:
print(f"Error! Be sure to run the makefile with 'make' to compile PLUME before using this!")
raise
tempfile.close()

return np.asarray(roots)

Expand Down
23 changes: 13 additions & 10 deletions src/fpc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -374,13 +374,13 @@ subroutine compute_fpc_cart(wrootindex)
if (computemoment) then
call calc_fs1(omega, vperp, vvz(ivz), phi, ef, bf, hatV_s(is), spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, elecdircontribution, &
(1., 0.), fs0(ivx, ivy, ivz, is), fs1_SP(ivx, ivy, ivz, is), EpsilonSokhotski_Plemelj)
(1., 0.), fs0(ivx, ivy, ivz, is), fs1_SP(ivx, ivy, ivz, is), EpsilonSokhotski_Plemelj,spec(is)%nu_ns)
!We fs1_sp to correctly compute jiEi = int CorEi d3v with residual!
fs1(ivx, ivy, ivz, is) = fs1_SP(ivx, ivy, ivz, is)
else
call calc_fs1(omega, vperp, vvz(ivz), phi, ef, bf, hatV_s(is), spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, elecdircontribution, &
(1., 0.), fs0(ivx, ivy, ivz, is), fs1(ivx, ivy, ivz, is), 0.)
(1., 0.), fs0(ivx, ivy, ivz, is), fs1(ivx, ivy, ivz, is), 0.,spec(is)%nu_ns)
end if
end do
end do
Expand Down Expand Up @@ -1308,7 +1308,7 @@ subroutine compute_fpc_gyro(wrootindex)
call calc_fs1(omega, vvperp(ivperp), vvpar(ivpar), vvphi(ivphi), ef, bf, hatV_s(is),&
spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, &
elecdircontribution, exbar, fs0(ivperp, ivpar, ivphi, is), fs1(ivperp, ivpar, ivphi, is), 0.)
elecdircontribution, exbar, fs0(ivperp, ivpar, ivphi, is), fs1(ivperp, ivpar, ivphi, is), 0.,spec(is)%nu_ns)

!compute fs1 at adjacent locations in vperp1/vperp2 direction to take derivatives with later
!Note: delv may not be the best choice here when it is large.
Expand All @@ -1321,7 +1321,7 @@ subroutine compute_fpc_gyro(wrootindex)
call calc_fs1(omega, vperp_adjacent, vvpar(ivpar), phi_adjacent, ef, bf, &
hatV_s(is), spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, elecdircontribution, &
exbar, fs0_temp, fs1_plus_delvperp1(ivperp, ivpar, ivphi, is), 0.)
exbar, fs0_temp, fs1_plus_delvperp1(ivperp, ivpar, ivphi, is), 0.,spec(is)%nu_ns)
vperp1_adjacent = vvperp(ivperp)*COS(vvphi(ivphi)) - delv
vperp2_adjacent = vvperp(ivperp)*SIN(vvphi(ivphi))
phi_adjacent = ATAN2(vperp2_adjacent, vperp1_adjacent)
Expand All @@ -1330,7 +1330,7 @@ subroutine compute_fpc_gyro(wrootindex)
call calc_fs1(omega, vperp_adjacent, vvpar(ivpar), phi_adjacent, ef, bf, &
hatV_s(is), spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, elecdircontribution, &
exbar, fs0_temp, fs1_minus_delvperp1(ivperp, ivpar, ivphi, is), 0.)
exbar, fs0_temp, fs1_minus_delvperp1(ivperp, ivpar, ivphi, is), 0.,spec(is)%nu_ns)
vperp1_adjacent = vvperp(ivperp)*COS(vvphi(ivphi))
vperp2_adjacent = vvperp(ivperp)*SIN(vvphi(ivphi)) + delv
phi_adjacent = ATAN2(vperp2_adjacent, vperp1_adjacent)
Expand All @@ -1339,7 +1339,7 @@ subroutine compute_fpc_gyro(wrootindex)
call calc_fs1(omega, vperp_adjacent, vvpar(ivpar), phi_adjacent, ef, bf, &
hatV_s(is), spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, elecdircontribution, &
exbar, fs0_temp, fs1_plus_delvperp2(ivperp, ivpar, ivphi, is), 0.)
exbar, fs0_temp, fs1_plus_delvperp2(ivperp, ivpar, ivphi, is), 0.,spec(is)%nu_ns)
vperp1_adjacent = vvperp(ivperp)*COS(vvphi(ivphi))
vperp2_adjacent = vvperp(ivperp)*SIN(vvphi(ivphi)) - delv
phi_adjacent = ATAN2(vperp2_adjacent, vperp1_adjacent)
Expand All @@ -1348,7 +1348,7 @@ subroutine compute_fpc_gyro(wrootindex)
call calc_fs1(omega, vperp_adjacent, vvpar(ivpar), phi_adjacent, ef, bf, &
hatV_s(is), spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, elecdircontribution, &
exbar, fs0_temp, fs1_minus_delvperp2(ivperp, ivpar, ivphi, is), 0.)
exbar, fs0_temp, fs1_minus_delvperp2(ivperp, ivpar, ivphi, is), 0.,spec(is)%nu_ns)
end do
end do
end do
Expand Down Expand Up @@ -1632,7 +1632,7 @@ end function fs0hat
! Collin Brown and Greg Howes, 2025
!------------------------------------------------------------------------------
subroutine calc_fs1(omega, vperp, vpar, phi, ef, bf, hatV_s, q_s, aleph_s, tau_s, mu_s, &
aleph_r, elecdircontribution, exbar, fs0, fs1, epsSokhotski_Plemelj)
aleph_r, elecdircontribution, exbar, fs0, fs1, epsSokhotski_Plemelj, nu_ns)
!! Determine species perturbed VDF fs1 at given (vperp,vpar,phi)

use vars, only: betap, kperp, kpar, vtp, pi, delv
Expand Down Expand Up @@ -1689,6 +1689,9 @@ subroutine calc_fs1(omega, vperp, vpar, phi, ef, bf, hatV_s, q_s, aleph_s, tau_s
real, intent(in) :: epsSokhotski_Plemelj
!! Small value for Sokhotski_Plemelj when computing moments of fs1; zero when just computing fs1

real, intent(in) :: nu_ns
!! Neutral-species collision frequency (using Krook collision operator)

real :: epsSokhotski_Plemelj_temp
!! fixes sign of epsilon

Expand Down Expand Up @@ -1786,7 +1789,7 @@ subroutine calc_fs1(omega, vperp, vpar, phi, ef, bf, hatV_s, q_s, aleph_s, tau_s
do n = -nbesmax, nbesmax
!Calculate all parts of solution that dosn't depend on m
!epsSokhotski_Plemelj is typically 0 unless using Sokhotski-Plemelj theorem to take moment over this singularity
denom = (omega_temp - kpar_temp*vpar_temp*sqrt(mu_s/(tau_s*aleph_r)) - n*mu_s/q_s) + (0., 1.)*epsSokhotski_Plemelj_temp
denom = (omega_temp - kpar_temp*vpar_temp*sqrt(mu_s/(tau_s*aleph_r)) - n*mu_s/q_s + (0., 1.) * nu_ns) + (0., 1.)*epsSokhotski_Plemelj_temp
Wbar_s = 2.*(n*mu_s/(q_s*(omega_temp)) - 1.)*(vpar_temp - hatV_s) - 2.*(n*mu_s/(q_s*(omega_temp)*aleph_s))*vpar_temp
emult = (0., 0.)

Expand Down Expand Up @@ -1969,7 +1972,7 @@ complex function wparth_from_ratio(is,ef)
phi = ATAN2(vvy(ivy), vvx(ivx))
call calc_fs1(omega, vperp, vvz(ivz), phi, ef, ef, hatV_s, spec(is)%q_s, spec(is)%alph_s, &
spec(is)%tau_s, spec(is)%mu_s, spec(1)%alph_s, 1., &
(1., 0.), fs0(ivx, ivy, ivz, is), fs1_SP(ivx, ivy, ivz, is), EpsilonSokhotski_Plemelj)
(1., 0.), fs0(ivx, ivy, ivz, is), fs1_SP(ivx, ivy, ivz, is), EpsilonSokhotski_Plemelj, spec(is)%nu_ns)
end do
end do
end do
Expand Down
10 changes: 7 additions & 3 deletions src/functions.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ subroutine read_in_params

!default values
elecdircontribution = 0.
Kn = 1.0e30 !give default value for backwards compat (note large -> collisionless)

call get_unused_unit(input_unit_no)
call get_runname(runname)
Expand Down Expand Up @@ -189,9 +190,12 @@ subroutine spec_read(is)
spec(is)%vv_s = vvS

!calculate neutral-charged collision frequency
spec(is)%nu_ns = (sqrt(2.d0)*Kn)**(-1.d0)* &
sqrt(spec(is)%mu_s/(spec(is)%tau_s*spec(1)%alph_s))

if (Kn > 1.0e25) then
spec(is)%nu_ns = 0.0 !Treat large default Kn as effectively collisionless (streamlines fpc.f90)
else
spec(is)%nu_ns = 1.0/(sqrt(2.0)*Kn) * &
sqrt(spec(is)%mu_s/(spec(is)%tau_s*spec(1)%alph_s))
endif
end subroutine spec_read

!-=-=-=-=
Expand Down
Loading