Implement a dispersion relation solver

[namurphy]

When dealing with plasma waves, it would be helpful to have a dispersion relation solver.

This was recommended by John Raymond. I believe Carl Sovinec wrote a code that this several years ago, though not in python.

[tulasinandan]

I have a very simple two fluid dispersion (e.g. Stringer JNE 1963 or Rogers PRL 2002) solver. I could include that in. In the long run we could include a full kinetic solver in here too.

[namurphy]

I also talked with Kristopher Klein about a general dispersion solver written in Fortran at a conference in late 2016, though I wasn't able to find it online. We would need to port that to Python or maybe Cython, however.

[tulasinandan]

I also have Peter Gary's dispersion solver. I am currently working on a python wrapper that simply calls the executable. We could bundle 64 bit static binary along with the python wrapper until we have something that is more in line with our goals. IRFU people do similar thing with WHAMP and matlab. I would have to take Peter's permission for it though.

Kris and Daniel (Verscharen) are working on an arbitrary linear dispersion solver that they call ALPS. We could coordinate with them so that their dispersion solver could be directly compiled and imported into python.

[namurphy]

Thinking about this more...all the dispersion relation solvers I've heard of are numerical. I am wondering if we could implement something with SymPy to analytically find dispersion relationships as well, at least for situations where an analytic solution is possible. I think it would be possible, at the very least. Just a thought!

[StanczakDominik]

I've been practicing SymPy recently and I could try to implement something in that. I could use a pointer to some reference materials for the dispersion relations themselves, though.

[antoinetavant]

Hi there, My Lab is launching a project to develop a kinetic relation dispertion solver for magnetised plasma.
As it will be under free licence, it could be implemented in some way into plasmaPy !
It should be finished around mid-july !

However, we could implement natively into plasmaPy some simpler dispersion relation solver...

[lemmatum]

What is a dispersion solver and how is it different from simple dispersion relations based on the Faddeeva function?

[tulasinandan]

@antoinelpp Great! As I said earlier as well, two fluid version is almost trivial and I have a small function for that. Kinetic one is hard, especially given the computational expense in some hard parameter regimes.

@namurphy Kris Klien and Daniel Verscharen have a very general linear dispersion solver in Fortran that they plan to make open source. We should ask them to make it an affiliate package with wrappers.

@lemmatum I did not know the term Faddeeva function before your post! I believe it is the same as plasma dispersion function, right? It is different from the "dispersion relations" for different wave modes. The first one appears in the permittivity, the second one describes a relationship between wavenumber and frequency for a given wave.

[lemmatum]

@tulasinandan I see, yes, there is a difference between dispersion function and dispersion relation. I've already opened an issue for implementing dispersion functions, and I've got code to do, I just haven't gotten around to opening a PR (especially as it is for PlasmaPy v0.2).

There is this really handy wikipedia page which has a table of dispersion relations, we should use this as a reference for implementation.

[dstansby]

https://github.com/danielver02/NHDS has recently been released, and I started working on building a python wrapper for it https://github.com/dstansby/NHDSPy

but then I realised the fortran needs to be recompiled each time there's a new set of input parameters which is kinda annoying. I might still try and finish the wrapper when I get some spare time.

[namurphy]

Looks like we may run into some license issues. NHDS has been released under the GPLv3, which means that all derivative works must be released under the GPLv3 as well, and thus we won't be able to incorporate it directly into PlasmaPy unless it's released under something like an MIT or BSD license.

[tulasinandan]

I can talk to Daniel about the licensing issues. I can also talk to him about restructuring the code so that you compile it once and read parameters from an input file. Although including this or even Daniel's other code (ALPS) would be against plasmapy's goal of keeping the code c/fortran free, we could (should) include both NHDS and ALPS as affiliate packages via wrappers. So I still think it is a good and worthwhile effort.

…

On Tue, Jun 19, 2018 at 11:37 AM, Nick Murphy ***@***.***> wrote: Looks like we may run into some license issues. NHDS has been released under the GPLv3, which means that all derivative works must be released under the GPLv3 as well, and thus we won't be able to incorporate it directly into PlasmaPy unless it's released under something like an MIT or BSD license. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#11 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ACacZ7xfKq8lti7JHFzzQJ7BE_3PIUG1ks5t-Rq5gaJpZM4HBfQM> .

-- Tulasi Nandan Parashar, 217 Sharp Laboratory, 104 The Grn, University of Delaware, Newark, DE 19716 Office: +1-302-831-1498 Cell: +1-302-766-7938 http://www.flickr.com/photos/tulasinandan

[antoinetavant]

Hi there,
I'm organising a hackathon mid July for all the physicists of my lab in order to recreate WHAMP in python, improving both the user interface but also the physics modelled.
We are planing to create a standalone project, but we are kind to include it in PlasmaPy.
I'll keep you in touch with the result, but it is in pretty good shape around July the 15th.

ps: the project is here, but it is empty for now.

[tulasinandan]

NHDS makes most of standard assumptions about the state of a kinetic plasma. Daniel Told's version is more approximate compared to NHDS as it assumes fluid electrons and kinetic ions. As I said before Daniel Verscharen has another code called ALPS which can work with arbitrary distribution functions. This one is really the most general kinetic dispersion solver on this planet that I know of. If one can approximate a distribution function as a superposition of a bunch of bi-maxwellians, then we can get the linear modes from ALPS. Daniel said that he will be releasing it soon. https://arxiv.org/pdf/1803.04697.pdf We should have the two fluid dispersion relations as well. Those are trivial to write. I have them in TurbPlasma which I plan to port to PEP 8 and Python 3.6 at some point of my life :-( For now I will copy paste the trivial code here itself in case someone wants to make it PEP8 and 3.6 compliant and include it in PlasmaPy

…

On Sat, Jun 30, 2018 at 5:15 AM, David Stansby ***@***.***> wrote: I think that raises an interesting point, that there is more than one way of writing a dispersion solver depending on the assumptions that go into deriving the physics. *If* PlasmaPy chooses a dispersion solver to have as some sort of default, I think it should as a minimum have - A simple and common set of physics assumptions - Be open source (obviously) - Have a peer reviewed paper published explaining the code and verifying that it works NHDS certainly fills criteria 2 and 3, but I'm not sure about 1. I'll continue working on making a python wrapper for it, and in the meantime it might be good to have a discussion about whether to "choose" a dispersion solver for PlasmaPy, or whether to link to all open source python dispersion solvers out there. Maybe I could jump on a meeting to discuss at some point? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#11 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ACacZ0x9pVAf9pDzq6zJwbZSk2QwdP15ks5uB0GxgaJpZM4HBfQM> .

-- Tulasi Nandan Parashar, 217 Sharp Laboratory, 104 The Grn, University of Delaware, Newark, DE 19716 Office: +1-302-831-1498 Cell: +1-302-766-7938 http://www.flickr.com/photos/tulasinandan

[tulasinandan]

EDIT: Add python tag to code block for syntax highlighting and add import numpy as np
EDIT 2: FIx reference from BALE JGR 2012 to BELLAN JGR 2012

#
#
# FUNCTION TO CALCULATE PHASE SPEEDS OF THE THREE BRANCHES OF TWO FLUID
# DISPERSION RELATION (e.g. STRINGER JPP 1963, ROGERS PRL 2001, 
# and BELLAN JGR 2012)
#
#                   Tulasi Nandan Parashar
import numpy as np
def tfps(beta=0.6, ca=1., de2=0.000545, theta=0., kk=1.):
   """
      beta: Total plasma beta,
      ca: Alfven speed based on mean field, 
      de2: me/mi, 
      theta: Angle of propagation in degrees
      kk: Wavenumber of interest in units of kdi
      
      Output is frequencies of the roots and the phase speeds w/k
      The roots are w[0]:Fast/Whistler, w[1]:Alfven/KAW, w[2]: Slow/Cyclotron
   """
   tht = theta*pi/180.
   ct = np.cos(tht)
   st = np.sin(tht)
   tt = st/ct
   cs=np.sqrt(beta/2.)*ca
   di =1.
   caksq=np.cos(tht)**2*ca**2
   cmsq=ca**2 + cs**2
   
   pcs=np.zeros(4)
   D = 1 + kk**2*de2
   # Find out the root of the quadratic dispersion relation
   pcs[0] = 1.
   pcs[1] = -(ca**2/D + cs**2 + caksq*(1+kk**2*di**2/D)/D)*kk**2
   pcs[2] = caksq*kk**4*(ca**2/D + cs**2 + cs**2*(1+kk**2*di**2/D))/D
   pcs[3] = -(cs**2*kk**6*caksq**2)/D**2
   w2 = np.roots(pcs); w = np.sqrt(w2)
   speeds= w/kk
   return w,speeds
#
# Compute the dispersion relation w(k) vs k for a single theta value.
#
def tfst(beta=0.6, ca=1., de2=0.000545, theta=0., kmin=1e-2, kmax=10., npoints=200,wrt='n'):
   """
      beta: Total plasma beta,
      ca: Alfven speed based on mean field, 
      de2: me/mi, 
      theta: Angle of propagation in degrees
      kkmin: Minimum Wavenumber of interest in units of kdi
      kkmax: Maximum wavenumber of interest in units of kdi

      Output is an array with 4 columns, k, w-fast, w-alf, w-slow
   """
   import matplotlib.pyplot as plt
   kmmn=np.log10(kmin)
   kmmx=np.log10(kmax)
   kk=np.logspace(kmmn,kmmx,npoints)
   warray=np.zeros((3,npoints))
   for i in range(0,npoints):
      f,s = tfps(beta, ca, de2, theta, kk[i])
      warray[:,i]=f
   plt.loglog(kk,warray[0,:], label='Fast/Magnetosonic')
   plt.loglog(kk,warray[1,:], label='Alfven/KAW')
   plt.loglog(kk,warray[2,:], label='Slow')
   plt.xlabel('$kd_i$')
   plt.ylabel('$\omega/\omega_{ci}$')
   plt.legend(loc='best',fancybox=True,framealpha=0.2)
   plt.title('Dispersion Relation for beta='+str(beta)+' and me/mi='+str(de2))
   plt.show()
   if wrt == 'y':
      ofile=open('disp'+str(theta)+'.dat','w')
      print>> ofile,'#  k', 'Acoustic', 'Alfven', 'Fast'
      for i in range(npoints):
         print>> ofile, kk[i],warray[2,i],warray[1,i],warray[0,i]
      ofile.close()

def tfev(beta=0.6, ca=1., de2=0.000545, theta=0., k=1.,aa=0.1):
   """
      beta: Total plasma beta,
      ca: Alfven speed based on mean field, 
      de2: me/mi, 
      theta: Angle of propagation in degrees
      k: wavenumber of interest in units of kdi

      Output: Prints the two fluid eigenvector to the screen
   """
   def amp(beta,de2,k,w,theta,aa):
      th=theta*pi/180.
      bb=1-w**2*(1+de2*k**2)/(k**2*np.cos(th)**2)
      sk='sin('+str(round(k,3))+'x)'
      ck='cos('+str(round(k,3))+'x)'
      def st(a):
         return str(round(a,3))
      return 'bx=0.'+\
            '\tby = '+st(round(2*aa,3))+ck+\
            '\tbz = '+st(-np.cos(th)*bb*2*aa/w)+sk+\
            '\nux = '+st(aa*k*bb*np.sin(2*th)/(w**2-beta*k**2))+sk+\
            '\tuy = '+st(-2*aa*k*np.cos(th)/w)+ck+\
            '\tuz = '+st(2*aa*k*bb*np.cos(th)**2/w**2)+sk+\
            '\t n = '+st((k*np.cos(th)/w)*aa*k*bb*np.sin(2*th)/(w**2-beta*k**2))+sk
   f,s=tfps(beta,ca,de2,theta,k)
#     The roots are w[0]:Fast/Whistler, w[1]:Alfven/KAW, w[2]: Slow/Cyclotron
   print 'Fast/ Whistler'
   print amp(beta,de2,k,f[0],theta,aa)
   print '##################'
   print
   print 'Alfven/ KAW'
   print amp(beta,de2,k,f[1],theta,aa)
   print '##################'
   print
   print 'Slow/ Cyclotron'
   print amp(beta,de2,k,f[2],theta,aa)
   

[StanczakDominik]

@ritiek, the later part of the telecon today was kinda hectic and I didn't catch everything - did you want to work on this one, or should I?

[ritiek]

@StanczakDominik I'll be a joke for implementing physics stuff but I could help a bit with coding style and python 3.6+ compatibility.

[StanczakDominik]

Then let's just collaborate on this 😄 I'll try to start on this tomorrow.

[tulasinandan]

Riteik, You don't have to worry about physics issues. If you could make my snippet of the code 3.6 and PEP8 compliant, that should be enough to include the two fluid dispersion in the code. The code is tested for physics already. Thanks, Tulasi

…

On Tue, Jul 3, 2018, 2:50 PM Dominik Stańczak ***@***.***> wrote: Then let's just collaborate on this 😄 I'll try to start on this tomorrow. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#11 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ACacZ5FfYih1lbiExumrlZGHSXxLDJ0tks5uC7z-gaJpZM4HBfQM> .

[ritiek]

@tulasinandan Yep! You can add all your code in a file named something like dispersion.py (?) and put it in the directory PlasmaPy/plasmapy/physics/, make the changes on your fork and then make a pull request to two-fluid-dispersion branch on my fork.

I'll apply fixes to your code and then further make a pull request to main PlasmaPy repo.

(We could directly copy paste and make your code python 3.6+/pep8 compliant and save you this trouble but it would be nicer for git to display that you wrote this code :)

[ritiek]

Gentle ping @tulasinandan. :D

[tulasinandan]

Submitted a pull request to your fork.

…

On Thu, Jul 5, 2018 at 6:25 AM Ritiek Malhotra ***@***.***> wrote: Gentle ping @tulasinandan <https://github.com/tulasinandan>. :D — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#11 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ACacZ9i_mgo3THm_BrEvHMm2If2rdc_Nks5uDemZgaJpZM4HBfQM> .

-- Tulasi Nandan Parashar, 217 Sharp Laboratory, 104 The Grn, University of Delaware, Newark, DE 19716 Office: +1-302-831-1498 Cell: +1-302-766-7938 http://www.flickr.com/photos/tulasinandan

[hsxie]

Hi all. I have developed several codes relevant to plasma dispersion relation, which could be advanced than some other solvers.

1. GPDF (http://hsxie.me/codes/gpdf/): generalized plasma dispersion function, which is very fast and accurate, and can support arbitrary 1D distribution function F(v) with an one-solve-all approach, instead of the conventional only Maxwellian one.

2. PDRF (http://hsxie.me/codes/pdrf/): A general dispersion relation solver for magnetized multi-fluid plasma, which can give all the exact solutions and corresponding polarizations at once for multi-fluid plasma model very fast. The solver does not require initial guess. For example, the cold plasma wave can be produced easily for arbitrary input parameters with arbitrary species.

3. PDRK (http://hsxie.me/codes/pdrk/): A General Kinetic Dispersion Relation Solver for Magnetized Plasma, which can give all the exact solutions (except strong damped modes) and corresponding polarizations at once for drift bi-Maxwellian distribution kinetic plasma model very fast. It could be advanced than the conventional root finding solvers, such as WHAMP.

I think these solvers can be very helpful for many researchers, especially in space plasma. Could them be possible implemented to plasmapy? I did not use python too much.

Thanks.
-hsxie

[liangwang0734]

I implemented Huasheng's fluid and kinetic (electrostatic only for now) solvers in Python a while ago. I plan to issue a pull request in the near future.

[joglekara]

Hi all!

I'm sorry if I missed something but I wasn't able to figure out what the status is here.

I have a scipy-based electrostatic dispersion relation solver that I'm happy to contribute. Please let me know if this is of interest!

[namurphy]

I have a scipy-based electrostatic dispersion relation solver that I'm happy to contribute. Please let me know if this is of interest!

Yes, let's definitely talk! We're currently planning for work on the dispersion subpackage to be a project for a graduate student starting sometime this summer. We've done some planning for it during our weekly community meetings which isn't on this issue, but overall we're still just getting started with this. We'd be interested to learn more about your electrostatic dispersion relation solver.

[namurphy]

We also very recently created a chat room on Riot for discussing the dispersion relation solver.

[liangwang0734]

I made one, too, based on matrix inversion instead of searching. I'm not familiar with PlasmaPy framework, but the code itself is quite straightforward and should be easy to integrate.

https://github.com/liangwang0734/xeon

The basic algorithm was originally developed by Dr. Huasheng Xie of ENN.

[tulasinandan]

So xeon is a python translation of BO? -- Tulasi Nandan Parashar, Lecturer - Physics, 503 Laby Building, Gate 7, Kelburn Parade, Wellington 6012 Office Phone: +64-4-463-5804 http://www.flickr.com/photos/tulasinandan

…

On Fri, Apr 24, 2020 at 9:45 AM Liang Wang ***@***.***> wrote: I made one, too, based on matrix inversion instead of searching. I'm not familiar with PlasmaPy framework, but the code itself is quite straightforward and should be easy to integrate. https://github.com/liangwang0734/xeon — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#11 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AATJYZ3IQFKEG25PHSZ7GD3ROCZGXANCNFSM4BYF6QGA> .

[liangwang0734]

Yes. With a couple of minor accuracy fixes.

…

On Thu, Apr 23, 2020, 18:14 tulasinandan ***@***.***> wrote: So xeon is a python translation of BO? -- Tulasi Nandan Parashar, Lecturer - Physics, 503 Laby Building, Gate 7, Kelburn Parade, Wellington 6012 Office Phone: +64-4-463-5804 http://www.flickr.com/photos/tulasinandan On Fri, Apr 24, 2020 at 9:45 AM Liang Wang ***@***.***> wrote: > I made one, too, based on matrix inversion instead of searching. I'm not > familiar with PlasmaPy framework, but the code itself is quite > straightforward and should be easy to integrate. > > https://github.com/liangwang0734/xeon > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#11 (comment)>, > or unsubscribe > < https://github.com/notifications/unsubscribe-auth/AATJYZ3IQFKEG25PHSZ7GD3ROCZGXANCNFSM4BYF6QGA > > . > — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#11 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAHUIY7XT3TAYCZ5FNULQFLROC4TTANCNFSM4BYF6QGA> .

[tulasinandan]

That is awesome! You made our life very easy! Our plan is to first have the graduate students learn the workflow: PEP8 coding, proper docstrings, testing, going through PRs etc by porting the simple two fluid version to PlasmaPy. Once that is done, we will get in touch with you to port this into PlasmaPy. Thank you for getting in touch with us! -- Tulasi Nandan Parashar, Lecturer - Physics, 503 Laby Building, Gate 7, Kelburn Parade, Wellington 6012 Office Phone: +64-4-463-5804 http://www.flickr.com/photos/tulasinandan On Fri, Apr 24, 2020 at 11:15 AM Liang Wang <notifications@github.com> wrote:

…

Yes. With a couple of minor accuracy fixes. On Thu, Apr 23, 2020, 18:14 tulasinandan ***@***.***> wrote: > So xeon is a python translation of BO? > -- > Tulasi Nandan Parashar, > Lecturer - Physics, > 503 Laby Building, > Gate 7, Kelburn Parade, > Wellington 6012 > Office Phone: +64-4-463-5804 > http://www.flickr.com/photos/tulasinandan > > > On Fri, Apr 24, 2020 at 9:45 AM Liang Wang ***@***.***> > wrote: > > > I made one, too, based on matrix inversion instead of searching. I'm not > > familiar with PlasmaPy framework, but the code itself is quite > > straightforward and should be easy to integrate. > > > > https://github.com/liangwang0734/xeon > > > > — > > You are receiving this because you were mentioned. > > Reply to this email directly, view it on GitHub > > <#11 (comment) >, > > or unsubscribe > > < > https://github.com/notifications/unsubscribe-auth/AATJYZ3IQFKEG25PHSZ7GD3ROCZGXANCNFSM4BYF6QGA > > > > . > > > > — > You are receiving this because you commented. > Reply to this email directly, view it on GitHub > <#11 (comment)>, > or unsubscribe > < https://github.com/notifications/unsubscribe-auth/AAHUIY7XT3TAYCZ5FNULQFLROC4TTANCNFSM4BYF6QGA > > . > — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#11 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AATJYZY3HHCT6YKLGK2J6VLRODDZRANCNFSM4BYF6QGA> .

[namurphy]

There's also a dispersion relation solver from Ammar Hakim.

[namurphy]

If you want to follow up on this, there's more discussion about dispersion relation solvers in #1237 and #2471.