Module: CodeRunner::Gs2::GSLVectors

Included in:
CodeRunner::Gs2
Defined in:
lib/gs2crmod/gsl_data.rb

Instance Method Summary collapse

Instance Method Details

#apar2_over_time_gsl_vector(options) ⇒ Object 



224
225
226
227
228
229
230
231
232
233
234
235
 
# File 'lib/gs2crmod/gsl_data.rb', line 224

def apar2_over_time_gsl_vector(options)

  Dir.chdir(@directory) do      #Necessary options: ky
    #log 'about to open netcdf file'
    #options.setup_time_window
    phis = netcdf_file.var('apar2').get('start'=>[options[:begin_element]], 'end'=>[options[:end_element]] ).to_a
    log 'about to allocate gsl vector'
    vec = GSL::Vector.alloc(phis)
    log 'finished'
    return fix_norm(vec, 1, options)
  end
end

#drhodpsi_gsl_vector(options) ⇒ Object

This function reads in the ‘drhodpsi’ variable from the netcdf file.



1023
1024
1025
1026
 
# File 'lib/gs2crmod/gsl_data.rb', line 1023

def drhodpsi_gsl_vector(options)
  drhodpsi = netcdf_file.var('drhodpsi').get()[0]
  return drhodpsi
end

#dt_gsl_vector(options) ⇒ Object

The size of each time step, indexed by time, normalised to a/v_th1.



355
356
357
358
359
360
 
# File 'lib/gs2crmod/gsl_data.rb', line 355

def dt_gsl_vector(options)
  t = gsl_vector('t', options)
  size = t.size
  # NB t already has norm fixed
  return t.subvector(1, size - 1) - t.subvector(0, size-1)
end

#es_heat_flux_by_kx_over_time_gsl_vector(options) ⇒ Object 



400
401
402
403
 
# File 'lib/gs2crmod/gsl_data.rb', line 400

def es_heat_flux_by_kx_over_time_gsl_vector(options)
  options[:direction] = :kx
  es_heat_flux_by_kxy_over_time_gsl_vector(options)
end

#es_heat_flux_by_ky_over_time_gsl_vector(options) ⇒ Object 



405
406
407
408
 
# File 'lib/gs2crmod/gsl_data.rb', line 405

def es_heat_flux_by_ky_over_time_gsl_vector(options)
  options[:direction] = :ky
  es_heat_flux_by_kxy_over_time_gsl_vector(options)
end

#es_heat_flux_over_kx_gsl_vector(options) ⇒ Object 



438
439
440
441
 
# File 'lib/gs2crmod/gsl_data.rb', line 438

def es_heat_flux_over_kx_gsl_vector(options)
  options[:direction] = :kx
  es_heat_flux_over_kxy_gsl_vector(options)
end

#es_heat_flux_over_kxy_gsl_vector(options) ⇒ Object

This function will output the heat flux as a function of kx or ky. Default behaviour will be to average the heat flux over the time domain.



450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
 
# File 'lib/gs2crmod/gsl_data.rb', line 450

def es_heat_flux_over_kxy_gsl_vector(options)
  Dir.chdir(@directory) do
    kxy = options[:direction]
    raise "Please provide species_index " unless options[:species_index]
    if kxy==:ky
      es_heat = (netcdf_file.var('es_heat_flux_by_mode').get({'start' => [0,0,options[:species_index]-1, 0], 'end' => [-1,-1,options[:species_index]-1, -1]})) #index = [kx,ky,spec,t]
      #Need to average over time and sum over kx
      shape = es_heat.shape
      es_heat_av = []; temp = [];
      for iy in 0...shape[1]
        for ix in 0...shape[0]
          temp[ix] = es_heat[ix,iy,0,0..-1].sum / shape[3]
        end
        es_heat_av[iy] = temp.sum
      end
      return es_heat_av.to_gslv
    else
      es_heat = (netcdf_file.var('es_heat_flux_by_mode').get({'start' => [0,0,options[:species_index]-1, 0], 'end' => [-1,-1,options[:species_index]-1, -1]})) #index = [kx,ky,spec,t]
      shape = es_heat.shape
      es_heat_av = []; temp = [];
      for ix in 0...shape[0]
        for iy in 0...shape[1]
          temp[iy] = es_heat[ix,iy,0,0..-1].sum / shape[3]
        end
        es_heat_av[ix] = temp.sum
      end
      return es_heat_av.to_gslv.from_box_order
    end
  end
end

#es_heat_flux_over_ky_gsl_vector(options) ⇒ Object 



443
444
445
446
 
# File 'lib/gs2crmod/gsl_data.rb', line 443

def es_heat_flux_over_ky_gsl_vector(options)
  options[:direction] = :ky
  es_heat_flux_over_kxy_gsl_vector(options)
end

#es_heat_flux_over_time_gsl_vector(options) ⇒ Object 



825
826
827
828
829
830
831
 
# File 'lib/gs2crmod/gsl_data.rb', line 825

def es_heat_flux_over_time_gsl_vector(options)
  Dir.chdir(@directory) do

    options.setup_time_window
    return GSL::Vector.alloc(netcdf_file.var('es_heat_flux').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
  end
end

#es_heat_par_over_time_gsl_vector(options) ⇒ Object Also known as: es_heat_par_gsl_vector 



833
834
835
836
837
838
839
 
# File 'lib/gs2crmod/gsl_data.rb', line 833

def es_heat_par_over_time_gsl_vector(options)
  Dir.chdir(@directory) do

    options.setup_time_window
    return GSL::Vector.alloc(netcdf_file.var('es_heat_par').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
  end
end

#es_heat_perp_over_time_gsl_vector(options) ⇒ Object Also known as: es_heat_perp_gsl_vector 



842
843
844
845
846
847
848
 
# File 'lib/gs2crmod/gsl_data.rb', line 842

def es_heat_perp_over_time_gsl_vector(options)
  Dir.chdir(@directory) do

    options.setup_time_window
    return GSL::Vector.alloc(netcdf_file.var('es_heat_perp').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
  end
end

#es_mom_flux_over_time_gsl_vector(options) ⇒ Object 



859
860
861
862
863
864
 
# File 'lib/gs2crmod/gsl_data.rb', line 859

def es_mom_flux_over_time_gsl_vector(options)
  Dir.chdir(@directory) do
    options.setup_time_window
    return GSL::Vector.alloc(netcdf_file.var('es_mom_flux').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
  end
end

#es_part_flux_over_time_gsl_vector(options) ⇒ Object 



866
867
868
869
870
871
 
# File 'lib/gs2crmod/gsl_data.rb', line 866

def es_part_flux_over_time_gsl_vector(options)
  Dir.chdir(@directory) do
    options.setup_time_window
    return GSL::Vector.alloc(netcdf_file.var('es_part_flux').get('start' => [options[:species_index].to_i - 1, options[:begin_element]], 'end' => [options[:species_index].to_i - 1, options[:end_element]]).to_a.flatten)
  end
end

#frequency_by_kx_over_time_gsl_vector(options) ⇒ Object

The real frequency of the fluctuations, read from the .out file, indexed by time and normalised to vth_1/a. :ky_index or :kx_index must be specified in options.



307
308
309
310
 
# File 'lib/gs2crmod/gsl_data.rb', line 307

def frequency_by_kx_over_time_gsl_vector(options)
  options[:direction] = :kx
  frequency_by_kxy_over_time_gsl_vector(options)
end

#frequency_by_kxy_over_time_gsl_vector(options) ⇒ Object 



317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
 
# File 'lib/gs2crmod/gsl_data.rb', line 317

def frequency_by_kxy_over_time_gsl_vector(options)
  kxy = options[:direction]
  kxy_index = kxy + :_index
  kxys = get_list_of(kxy)
  desired_kxy = kxys[options[kxy_index]]
  raise "No k found at the desired index" if desired_kxy.nil?

  omega_reals = []
  File.open(@run_name+".out",'r') do |fileHandle|
    fileHandle.each_line do |fileLine|
      if fileLine.include?('aky=')  # Only examine the lines of the .out file that contain frequency information.

        index = fileLine.index('akx=')
        raise "akx wasn't found where it was expected in the .out file." if index.nil?
        akx = fileLine[(index+4)..-1].to_f

        index = fileLine.index('aky=')
        raise "aky wasn't found where it was expected in the .out file." if index.nil?
        aky = fileLine[(index+4)..-1].to_f

        index = fileLine.index('om=')
        raise "om wasn't found where it was expected in the .out file." if index.nil?
        omr = fileLine[(index+3)..-1].to_f
        if kxy == :kx
          # You need to be careful when testing equality of the desired k with the k in the .out file
          # since the .out file is only written to ~ 5 significant digits:
          omega_reals << omr if ((desired_kxy - akx).abs/(desired_kxy.abs + 1e-7) < 1e-4)
        else
          omega_reals << omr if ((desired_kxy - aky).abs/(desired_kxy.abs + 1e-7) < 1e-4)
        end
      end
    end
  end
  raise "No real frequencies found in the .out file for the desired k" if (omega_reals.size==0)
  GSL::Vector.alloc(omega_reals)
end

#frequency_by_ky_over_time_gsl_vector(options) ⇒ Object 



312
313
314
315
 
# File 'lib/gs2crmod/gsl_data.rb', line 312

def frequency_by_ky_over_time_gsl_vector(options)
  options[:direction] = :ky
  frequency_by_kxy_over_time_gsl_vector(options)
end

#growth_rate_by_kx_over_time_gsl_vector(options) ⇒ Object

The growth rate of the fluctuations, calculated from the potential, indexed by time and normalised to vth_1/a. :kx or :kx_index must be specified in options



280
281
282
283
 
# File 'lib/gs2crmod/gsl_data.rb', line 280

def growth_rate_by_kx_over_time_gsl_vector(options)
  options[:direction] = :kx
  growth_rate_by_kxy_over_time_gsl_vector(options)
end

#growth_rate_by_kxy_over_time_gsl_vector(options) ⇒ Object 



292
293
294
295
296
297
298
299
300
301
302
303
 
# File 'lib/gs2crmod/gsl_data.rb', line 292

def growth_rate_by_kxy_over_time_gsl_vector(options)
  # i.e. time_dependent_gr_by_ky_vs_time or phi2_by_kx_vs_time

  kxy = options[:direction]

  phi = gsl_vector("phi2_by_#{kxy}_over_time", options).log / 2.0

  size = phi.size
  dphi = phi.subvector(1, size - 1) - phi.subvector(0, size-1)
  # NB dt already has norm fixed, dphi is dimensionless
  return fix_norm(dphi/gsl_vector('dt'), 0, options)
end

#growth_rate_by_ky_over_time_gsl_vector(options) ⇒ Object

The growth rate of the fluctuations, calculated from the potential, indexed by time and normalised to vth_1/a. :ky or :ky_index must be specified in options



287
288
289
290
 
# File 'lib/gs2crmod/gsl_data.rb', line 287

def growth_rate_by_ky_over_time_gsl_vector(options)
  options[:direction] = :ky
  growth_rate_by_kxy_over_time_gsl_vector(options)
end

#growth_rate_over_kx_gsl_vector(options) ⇒ Object

The growth rate, calculated from the potential, indexed by kx. Only makes sense in linear calculations.



363
364
365
366
 
# File 'lib/gs2crmod/gsl_data.rb', line 363

def growth_rate_over_kx_gsl_vector(options)
  options[:direction] = :kx
  growth_rate_over_kxy_gsl_vector(options)
end

#growth_rate_over_kx_slice_gsl_vector(options) ⇒ Object

The growth rate, calculated from the potential, indexed by kx. Only makes sense in linear calculations.



382
383
384
385
386
387
388
389
 
# File 'lib/gs2crmod/gsl_data.rb', line 382

def growth_rate_over_kx_slice_gsl_vector(options)
  Dir.chdir(@directory) do
    slice_of_growth_rates = send(:growth_rate_at_ky_at_kx)[options[:ky]].values
    raise "Something went wrong: slice of growth rates seems empty" if slice_of_growth_rates.nil?
    return GSL::Vector.alloc(slice_of_growth_rates)
    #return GSL::Vector.alloc(send(:growth_rate_at_ky_at_kx[ky]).values)
  end
end

#growth_rate_over_ky_gsl_vector(options) ⇒ Object

The growth rate, calculated from the potential, indexed by ky. Only makes sense in linear calculations.



368
369
370
371
 
# File 'lib/gs2crmod/gsl_data.rb', line 368

def growth_rate_over_ky_gsl_vector(options)
  options[:direction] = :ky
  growth_rate_over_kxy_gsl_vector(options)
end

#growth_rate_over_ky_slice_gsl_vector(options) ⇒ Object

The growth rate, calculated from the potential, indexed by ky. Only makes sense in linear calculations.



392
393
394
395
396
397
398
 
# File 'lib/gs2crmod/gsl_data.rb', line 392

def growth_rate_over_ky_slice_gsl_vector(options)
  Dir.chdir(@directory) do
    slice_of_growth_rates = send(:growth_rate_at_ky_at_kx).values.map{|h| h[options[:kx]]}
    raise "Something went wrong: slice of growth rates seems empty" if slice_of_growth_rates.nil?
    return GSL::Vector.alloc(slice_of_growth_rates)
  end
end

#hflux_tot_over_time_gsl_vector(options) ⇒ Object Also known as: hflux_tot_gsl_vector 



813
814
815
816
817
818
819
820
821
822
 
# File 'lib/gs2crmod/gsl_data.rb', line 813

def hflux_tot_over_time_gsl_vector(options)
  Dir.chdir(@directory) do
    options.setup_time_window
    narr = netcdf_file.var('hflux_tot').get('start' => [options[:begin_element]], 'end' => [options[:end_element]])
    #eputs 'Got narr'
    #ep 'hflux_tot', hflux
    #eputs "fixing norm"
    return fix_heat_flux_norm(GSL::Vector.alloc(narr.to_a), options)
  end
end

#kpar_gsl_vector(options) ⇒ Object 



688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
 
# File 'lib/gs2crmod/gsl_data.rb', line 688

def kpar_gsl_vector(options)

  Dir.chdir(@directory) do
    if agk? or (@s_hat_input||@shat).abs  < 1.0e-5
      dk = 1
      phi = list(:theta).values
    else
      kxe = gsl_vector('linked_kx_elements', options)
      dk = 1.0/kxe.size
    phi = gsl_vector_complex('phi_along_field_line', options)
    end
    case phi.size%2
    when 0
      kpar = GSL::Vector.indgen(phi.size-1, -((phi.size-3)/2))*dk
    when 1
      kpar = GSL::Vector.indgen(phi.size-1, -((phi.size-2)/2))*dk
    end
    #ep 'kpar', kpar, 'phi.size', phi.size

    #ep 'kpar.class', kpar.class
    return kpar

  end
end

#linked_kx_elements_gsl_vector(options) ⇒ Object 



595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
 
# File 'lib/gs2crmod/gsl_data.rb', line 595

def linked_kx_elements_gsl_vector(options)
  Dir.chdir(@directory) do
    return GSL::Vector.alloc([0]) if @grid_option == "single" or agk?
    if agk? or (@s_hat_input or @shat).abs < 1.0e-5
      #p 'op1', options
      options.convert_to_index(self, :ky, :kx)
      #p 'op2', options
      #eputs "No Magnetic Shear"

#         begin
#           options.convert_to_index(:kx)
#         rescue
#           raise "Must specify kx or kx_index if no magnetics shear"
#         end
# #         theta0 = (options[:theta0] || 0)
# #         theta0 += jump(options) if @g_exb

      #theta0 = (options[:kx_index])
      #if @g_exb and @g_exb.abs > 0.0
        #theta0 += jump(options)
        #theta0 = theta0%((list(:kx).size-1)/2) if list(:kx).size > 1
      #end

      return GSL::Vector.alloc([options[:kx_index] - 1])
    end

    options.convert_to_index(self, :ky, :kx)
    nkx = netcdf_file.var('kx').dims[0].length
#       p nkx
    stride = @jtwist * (options[:ky_index] -1 )
    #stride = 3
    nlinks = [(nkx / stride).floor, 1].max
    theta0 = options[:kx_index] % @jtwist  #(options[:theta0] || 0)
    #log 'stride', stride, 'nlinks', nlinks, 'theta0', theta0
    #if @g_exb and @jtwist > 1 #and options[:t_index]
#         kx_shift = list(:ky)[options[:ky_index]]  * @g_exb
#         p list(:t)[options[:t_index]], options[:t_index], kx_shift

#         kx_shift *=  list(:t)[(options[:t_index] or list(:t).keys.max)]
#         jump = (kx_shift / list(:kx)[2]).round
      #theta0  += (@jtwist - jump(options) % @jtwist) % @jtwist

#         else
#           jump = 0
    #end
    #ep 'stride', stride, 'nlinks', nlinks, 'theta0', theta0
    #ep GSL::Vector.indgen(nlinks / 2,  nkx + theta0 - nlinks / 2 * stride, stride).connect(GSL::Vector.indgen(nlinks / 2, theta0, stride)).reverse if nlinks > 1
    #return [7,5,3,1,34].to_gslv
    return GSL::Vector.alloc([theta0 % jtwist]) if nlinks ==1
    return GSL::Vector.indgen(nlinks / 2,  nkx + theta0 - nlinks / 2 * stride, stride).connect(GSL::Vector.indgen(nlinks / 2, theta0, stride)).reverse

  end
end

#lpc_energy_gsl_vector(options) ⇒ Object

Velocity space diagnostics: fraction of dist func in higher energy harmonics



883
884
885
886
887
 
# File 'lib/gs2crmod/gsl_data.rb', line 883

def lpc_energy_gsl_vector(options)
  raise "Velocity space lpc diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.lpc"
  lpc = GSL::Vector.filescan("#@directory/#@run_name.lpc")
  return lpc[2]
end

#lpc_pitch_angle_gsl_vector(options) ⇒ Object

Velocity space diagnostics: fraction of dist func in higher pitch angle harmonics



875
876
877
878
879
 
# File 'lib/gs2crmod/gsl_data.rb', line 875

def lpc_pitch_angle_gsl_vector(options)
  raise "Velocity space lpc diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.lpc"
  lpc = GSL::Vector.filescan("#@directory/#@run_name.lpc")
  return lpc[1]
end

#mean_flow_velocity_over_x_gsl_vector(options) ⇒ Object

This function returns the mean flow velocity as a function of x (the radial coordinate). This is v_g_exb = (x - x(centre))*g_exb. The x-x(centre) ensures that the flow is zero at the middle of the box.



1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
 
# File 'lib/gs2crmod/gsl_data.rb', line 1063

def mean_flow_velocity_over_x_gsl_vector(options)
  Dir.chdir(@directory) do
    raise CRFatal.new("Need to have g_exb > 0 to have a mean flow.") unless @g_exb
    x = gsl_vector(:x)

    vec_exb_vel = GSL::Vector.alloc(x.size)
    #Take imaginary part since i k_x will lead to imaginary part being real
    vec_exb_vel = (x - x[x.size/2])*@g_exb
    return vec_exb_vel
  end
end

#par_mom_flux_over_time_gsl_vector(options) ⇒ Object 



905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
 
# File 'lib/gs2crmod/gsl_data.rb', line 905

def par_mom_flux_over_time_gsl_vector(options)
Dir.chdir(@directory) do

  options.setup_time_window
  # This is a hack... one day some one will put it in the NetCDF file (haha).
  momlines = `grep parmom #@run_name.out`
  mom = []
  momlines.scan(Regexp.new("#{LongRegexen::FLOAT.to_s}$")) do
    mom.push $~[:float].to_f
  end
  options[:end_element] = (mom.size + options[:end_element]) if options[:end_element] < 0
#       p options
  return GSL::Vector.alloc(mom).subvector(options[:begin_element], options[:end_element] - options[:begin_element] + 1)
end
end

#perp_mom_flux_over_time_gsl_vector(options) ⇒ Object 



921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
 
# File 'lib/gs2crmod/gsl_data.rb', line 921

def perp_mom_flux_over_time_gsl_vector(options)

  Dir.chdir(@directory) do
    options.setup_time_window
    # This is a hack... one day some one will put it in the NetCDF file (haha).
    momlines = `grep perpmom #@run_name.out`
    mom = []
    momlines.scan(Regexp.new("#{LongRegexen::FLOAT.to_s}$")) do
      mom.push $~[:float].to_f
    end
    options[:end_element] = (mom.size + options[:end_element]) if options[:end_element] < 0
#       p options
    return GSL::Vector.alloc(mom).subvector(options[:begin_element], options[:end_element] - options[:begin_element] + 1)
  end
end

#phi0_by_kx_by_ky_over_time_gsl_vector(options) ⇒ Object 



587
588
589
590
591
592
593
 
# File 'lib/gs2crmod/gsl_data.rb', line 587

def phi0_by_kx_by_ky_over_time_gsl_vector(options)
  Dir.chdir(@directory) do
    options.convert_to_index(self, :kx, :ky)
    phi0_array = netcdf_file.var('phi0').get.to_a.map{|arr| arr[options[:kx_index] - 1][options[:ky_index] - 1][options[:ri]]}
    return GSL::Vector.alloc(phi0_array)
  end
end

#phi2_by_kx_over_time_gsl_vector(options) ⇒ Object 



481
482
483
484
 
# File 'lib/gs2crmod/gsl_data.rb', line 481

def phi2_by_kx_over_time_gsl_vector(options)
  options[:direction] = :kx
  phi2_by_kxy_over_time_gsl_vector(options)
end

#phi2_by_ky_over_time_gsl_vector(options) ⇒ Object 



486
487
488
489
 
# File 'lib/gs2crmod/gsl_data.rb', line 486

def phi2_by_ky_over_time_gsl_vector(options)
  options[:direction] = :ky
  phi2_by_kxy_over_time_gsl_vector(options)
end

#phi2_by_mode_over_time_gsl_vector(options) ⇒ Object 



515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
 
# File 'lib/gs2crmod/gsl_data.rb', line 515

def phi2_by_mode_over_time_gsl_vector(options)
  Dir.chdir(@directory) do      #Necessary options: :ky and :kx
    #Optional options: :t_index_window
    #     eputs "got here"
    #options[:begin_element], options[:end_element] = (options[:t_index_window] ? options[:t_index_window].map{|ind| ind -1} : [0, -1])
    options.setup_time_window
    phi_t_array=nil
    if @grid_option == "single"
      phi_t_array = netcdf_file.var('phi2').get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a.flatten
    else
#         value = options[:ky]
#         eputs value
#         get_list_of(:ky)
#         index = @ky_list.find{|index,val| (val-value).abs < Float::EPSILON}[0]
      options.convert_to_index(self, :kx, :ky)
#         p options
      phi_t_array = netcdf_file.var("phi2_by_mode").get('start' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:begin_element]], 'end' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:end_element]]).to_a.flatten
#         eputs 'phi_t_array.size', phi_t_array.size
    end
    return GSL::Vector.alloc(phi_t_array)

  end
end

#phi2tot_over_time_gsl_vector(options) ⇒ Object

The square of the potential summed over all wave numbers, indexed by time, normalised to (e/T)(rho_1/a).



211
212
213
214
215
216
217
218
219
220
221
222
 
# File 'lib/gs2crmod/gsl_data.rb', line 211

def phi2tot_over_time_gsl_vector(options)

  Dir.chdir(@directory) do      #Necessary options: ky
    #log 'about to open netcdf file'
    #options.setup_time_window
    phis = netcdf_file.var('phi2').get('start'=>[options[:begin_element]], 'end'=>[options[:end_element]] ).to_a
    log 'about to allocate gsl vector'
    vec = GSL::Vector.alloc(phis)
    log 'finished'
    return fix_norm(vec, 1, options)
  end
end

#phi_along_field_line_gsl_vector(options) ⇒ Object 



713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
 
# File 'lib/gs2crmod/gsl_data.rb', line 713

def phi_along_field_line_gsl_vector(options)
  Dir.chdir(@directory) do
    complex_phi_vector= gsl_vector_complex('phi_along_field_line', options)
    case options[:imrc]
    when :im
      phi_vector = complex_phi_vector.imag
    when :mag
      _mag = true
      phi_vector = complex_phi_vector.abs2
    when :corr
      thetas = gsl_vector('theta_along_field_line', options)
      min = thetas.abs.to_a.index(thetas.abs.min)
      at_0 = complex_phi_vector[min]
#         ep at_0.class
      phi_vector = (complex_phi_vector * (at_0 / at_0.mag).conj).real
#         gsl_complex('correcting_phase', options)).real
    when :real
      phi_vector = complex_phi_vector.real
    else
      raise CRError.new("options[:imrc] was: #{options[:irmc]}")
    end
    phi_vector *= -1.0 if options[:flip]
    (phi_vector /= phi_vector.abs.max; phi_vector *= (options[:height] || 1.0)) if options[:norm]
    phi_vector = phi_vector.reverse if options[:rev]
    return phi_vector

  end
end

#phi_for_eab_movie_gsl_vector(options) ⇒ Object 



802
803
804
805
806
807
808
809
810
811
 
# File 'lib/gs2crmod/gsl_data.rb', line 802

def phi_for_eab_movie_gsl_vector(options)
  Dir.chdir(@directory) do      #options required are x_index, y_index and tm_index (Time)
    mvf_name = @run_name + '.movie.nc'
    raise CRError.new("cannot find file #{mvf_name}") unless FileTest.exist? mvf_name
    ncf = NumRu::NetCDF.open(mvf_name)
#       p ncf.var('phi_by_xmode').get.to_a[0][0][0]
    return GSL::Vector.alloc(ncf.var('phi_by_xmode').get.to_a[options[:tm_index] - 1].map{|xy_arr| xy_arr[options[:x_index] - 1][options[:y_index] - 1]})

  end
end

#scan_parameter_value_gsl_vector(options) ⇒ Object 



937
938
939
 
# File 'lib/gs2crmod/gsl_data.rb', line 937

def scan_parameter_value_gsl_vector(options)
  return GSL::Vector.alloc(netcdf_file.var('scan_parameter_value').get.to_a)
end

#spectrum_over_kpar_gsl_vector(options) ⇒ Object 



649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
 
# File 'lib/gs2crmod/gsl_data.rb', line 649

def spectrum_over_kpar_gsl_vector(options)
  Dir.chdir(@directory) do
 # , /kpar_spectrum/
    #ep 'zero?', (@s_hat_input||@shat)==0.0
    unless agk? or (@s_hat_input||@shat||0.0).abs<1.0e-5
      phi = gsl_vector_complex('phi_along_field_line', options)
      phi = phi.subvector(0,phi.size-1)
      #i = 0
      #phi = phi.collect{|re,im|
        #i+=1; GSL::Complex.alloc(Math.sin(0.1*i), Math.cos(0.1*i))+
        #GSL::Complex.alloc(Math.sin(0.4*i), Math.cos(0.4*i))

      #}
      ##GraphKit.quick_create([phi.square]).gnuplot
      phi_k = phi.forward
      phi_kr = phi_k.square
      case phi_kr.size%2
      when 0
        spec = phi_kr.subvector((phi_kr.size+2)/2, (phi_kr.size-2)/2).connect(phi_kr.subvector(0, (phi_kr.size+2)/2))
      when 1
        spec = phi_kr.subvector((phi_kr.size + 1)/2, (phi_kr.size-1)/2).connect(phi_kr.subvector(0, (phi_kr.size+1)/2))
      end
      ##spec = phi_kr
      #ep 'spec.class', spec.class
      return spec
    else

      gm = gsl_matrix('spectrum_over_ky_over_kpar', options)
      vec = GSL::Vector.alloc(gm.shape[1])
      vec.set_all(0.0)
      for ky_element in 0...gm.shape[0]
        vec+= gm.row(ky_element)
      end
      vec = vec/gm.shape[0]
      return vec
    end
  end
end

#spectrum_over_kx_avg_gsl_vector(options) ⇒ Object 



946
947
948
949
 
# File 'lib/gs2crmod/gsl_data.rb', line 946

def spectrum_over_kx_avg_gsl_vector(options)
  options[:direction] = :kx
  spectrum_over_kxy_avg_gsl_vector(options)
end

#spectrum_over_kx_gsl_vector(options) ⇒ Object 



941
942
943
944
 
# File 'lib/gs2crmod/gsl_data.rb', line 941

def spectrum_over_kx_gsl_vector(options)
  options[:direction] = :kx
  spectrum_over_kxy_gsl_vector(options)
end

#spectrum_over_kxy_avg_gsl_vector(options) ⇒ Object

spectrum averaged in time



982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
 
# File 'lib/gs2crmod/gsl_data.rb', line 982

def spectrum_over_kxy_avg_gsl_vector(options)
  Dir.chdir(@directory) do
    # i.e. spectrum_over_ky or spectrum_over_kx
    kxy = options[:direction]
    raise "Spectrum makes no sense for single modes" if @grid_option == "single"

    phi_array = netcdf_file.var("phi2_by_#{kxy}").get('start' => [0, 0], 'end' => [-1, -1]) #index = [kx or ky, t]

    shape = phi_array.shape
    phi_av = [];
    #average over time for each kx or ky individually
    for i in 0...shape[0]
      phi_av[i] = phi_array[i,0..-1].sum / shape[1]
    end

    v = GSL::Vector.alloc(phi_av)
    v = v.from_box_order if kxy == :kx
    v = v.mul(gsl_vector(kxy).square) unless options[:phi2_only]
    return v
  end
end

#spectrum_over_kxy_gsl_vector(options) ⇒ Object 



961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
 
# File 'lib/gs2crmod/gsl_data.rb', line 961

def spectrum_over_kxy_gsl_vector(options)
  Dir.chdir(@directory) do
    # i.e. spectrum_over_ky or spectrum_over_kx
    kxy = options[:direction]
#       eputs options[:t_index]
    raise "Spectrum makes no sense for single modes" if @grid_option == "single"

    options.convert_to_index(:t) if options[:t] or options[:t_element]
#       eputs options[:t_index]

    options[:t_index] ||= list(:t).keys.max
#       eputs options[:t_index]
    phi_array = netcdf_file.var("phi2_by_#{kxy}").get('start' => [0, options[:t_index] - 1], 'end' => [-1, options[:t_index] - 1]).to_a.flatten
    v = GSL::Vector.alloc(phi_array)
    v = v.from_box_order if kxy == :kx
    v = v.mul(gsl_vector(kxy).square) unless options[:phi2_only]
    return v
  end
end

#spectrum_over_ky_avg_gsl_vector(options) ⇒ Object 



956
957
958
959
 
# File 'lib/gs2crmod/gsl_data.rb', line 956

def spectrum_over_ky_avg_gsl_vector(options)
  options[:direction] = :ky
  spectrum_over_kxy_avg_gsl_vector(options)
end

#spectrum_over_ky_gsl_vector(options) ⇒ Object 



951
952
953
954
 
# File 'lib/gs2crmod/gsl_data.rb', line 951

def spectrum_over_ky_gsl_vector(options)
  options[:direction] = :ky
  spectrum_over_kxy_gsl_vector(options)
end

#theta_along_field_line_gsl_vector(options) ⇒ Object 



742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
 
# File 'lib/gs2crmod/gsl_data.rb', line 742

def theta_along_field_line_gsl_vector(options)
  Dir.chdir(@directory) do
    case @grid_option
    when "single", "range"
      theta_vector = gsl_vector(:theta)
    when "box"
      #eputs "Start theta_along_field_line"

      kx_elements = gsl_vector('linked_kx_elements', options).to_a
#         ep list(:kx).keys.max
#         ep kx_elements[0], list(:kx)[(kx_elements[0] + 1).to_i]
#         ep kx_elements[-1], list(:kx)[(kx_elements[-1] + 1).to_i]
      thetas = gsl_vector(:theta)
#         ep thetas
      #eputs "End theta_along_field_line"
      return thetas if agk? or (@s_hat_input or @shat).abs < 1.0e-5
      if gryfx?
        theta_list = ((1..kx_elements.size).to_a.map do |i|
          thetas * i
        end)
        thetas = theta_list.inject{|o,n| o.connect(n)}
        thetas -= Math::PI*(kx_elements.size-1)
        return thetas

      end
      theta_list = (kx_elements.map do |element|

        kx = list(:kx)[(element + 1).to_i]
#                      ep element
                   #ep 'kx', kx, 'shat', (@s_hat_input or @shat), 'ky',   list(:ky)[options[:ky_index]]
        thetas - 1.0 / (@s_hat_input or @shat) / list(:ky)[options[:ky_index]] * kx
      end).inject{|old, new| old.connect(new)}
#         thetas = gsl_vector(:theta) - 1.0 / @shat / list(:ky)[options[:ky_index]] * list(:kx)[(kx_elements[0] + 1).to_i] #- Math::PI*(kx_elements.size  - 1)
#         get_list_of(:ky, :t)
#         if @g_exb #and options[:t_index]

        if options[:moving]
          theta_list = theta_list  -  Math::PI * 2.0 * (jump(options) / @jtwist)
        else
#             ep 'jump % jtwist is!!', jump(options) % @jtwist
          theta_list = theta_list - Math::PI * 2.0 / @nx.to_f * ((jump(options) % @jtwist).to_f / @jtwist.to_f)
        end
#           jump = 0
#         end
#         theta_list = thetas.dup #gsl_vector(:theta) - Math::PI*kx_elements.size
#         (kx_elements.size - 1).times do
#           thetas = thetas + Math::PI * 2.0
#           theta_list = theta_list.connect(thetas)
#         end
#         pp theta_list.to_a.values_at(0, theta_list.size - 1)
#         pp theta_list.to_a.max
      theta_vector = theta_list
    end
#       theta_vector = theta_vector.reverse if options[:rev]
    theta_vector *= (@shat) if options[:z]
    return theta_vector

  end
end

#tpar2_by_mode_over_time_gsl_vector(options) ⇒ Object 



539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
 
# File 'lib/gs2crmod/gsl_data.rb', line 539

def tpar2_by_mode_over_time_gsl_vector(options)
  Dir.chdir(@directory) do      #Necessary options: :ky and :kx
    #Optional options: :t_index_window
    #     eputs "got here"
    #options[:begin_element], options[:end_element] = (options[:t_index_window] ? options[:t_index_window].map{|ind| ind -1} : [0, -1])
    options.setup_time_window
    tpar_t_array=nil
    if @grid_option == "single"
      tpar_t_array = netcdf_file.var('tpar2').get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a.flatten
    else
#         value = options[:ky]
#         eputs value
#         get_list_of(:ky)
#         index = @ky_list.find{|index,val| (val-value).abs < Float::EPSILON}[0]
      options.convert_to_index(self, :kx, :ky, :species)
#         p options
      tpar_t_array = netcdf_file.var("tpar2_by_mode").get('start' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:begin_element]], 'end' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:end_element]]).to_a.flatten
#         eputs 'tpar_t_array.size', tpar_t_array.size
    end
    return GSL::Vector.alloc(tpar_t_array)

  end
end

#tperp2_by_mode_over_time_gsl_vector(options) ⇒ Object 



563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
 
# File 'lib/gs2crmod/gsl_data.rb', line 563

def tperp2_by_mode_over_time_gsl_vector(options)
  Dir.chdir(@directory) do      #Necessary options: :ky and :kx
    #Optional options: :t_index_window
    #     eputs "got here"
    #options[:begin_element], options[:end_element] = (options[:t_index_window] ? options[:t_index_window].map{|ind| ind -1} : [0, -1])
    options.setup_time_window
    tperp_t_array=nil
    if @grid_option == "single"
      tperp_t_array = netcdf_file.var('tperp2').get('start' => [options[:begin_element]], 'end' => [options[:end_element]]).to_a.flatten
    else
#         value = options[:ky]
#         eputs value
#         get_list_of(:ky)
#         index = @ky_list.find{|index,val| (val-value).abs < Float::EPSILON}[0]
      options.convert_to_index(self, :kx, :ky, :species)
#         p options
      tperp_t_array = netcdf_file.var("tperp2_by_mode").get('start' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:begin_element]], 'end' => [options[:kx_index] - 1, options[:ky_index] - 1, options[:species_index] - 1, options[:end_element]]).to_a.flatten
#         eputs 'tperp_t_array.size', tperp_t_array.size
    end
    return GSL::Vector.alloc(tperp_t_array)

  end
end

#transient_amplification_over_kx_gsl_vector(options) ⇒ Object 



258
259
260
261
 
# File 'lib/gs2crmod/gsl_data.rb', line 258

def transient_amplification_over_kx_gsl_vector(options)
  options[:direction] = :kx
  transient_amplification_over_kxy_gsl_vector(options)
end

#transient_amplification_over_ky_gsl_vector(options) ⇒ Object 



263
264
265
266
 
# File 'lib/gs2crmod/gsl_data.rb', line 263

def transient_amplification_over_ky_gsl_vector(options)
  options[:direction] = :ky
  transient_amplification_over_kxy_gsl_vector(options)
end

#transient_es_heat_flux_amplification_over_kx_gsl_vector(options) ⇒ Object 



237
238
239
240
 
# File 'lib/gs2crmod/gsl_data.rb', line 237

def transient_es_heat_flux_amplification_over_kx_gsl_vector(options)
  options[:direction] = :kx
  transient_es_heat_flux_amplification_over_kxy_gsl_vector(options)
end

#transient_es_heat_flux_amplification_over_kxy_gsl_vector(options) ⇒ Object 



247
248
249
250
251
252
253
254
255
256
 
# File 'lib/gs2crmod/gsl_data.rb', line 247

def transient_es_heat_flux_amplification_over_kxy_gsl_vector(options)
  Dir.chdir(@directory) do      # i.e. phi2_by_ky_vs_time or phi2_by_kx_vs_time
    kxy = options[:direction].to_sym

#       ep :growth_rate_at_ + kxy
    p send(:transient_es_heat_flux_amplification_at_species_at_ + kxy)
    return GSL::Vector.alloc(send(:transient_es_heat_flux_amplification_at_species_at_ + kxy)[options[:species_index]-1].values)

  end
end

#transient_es_heat_flux_amplification_over_ky_gsl_vector(options) ⇒ Object 



242
243
244
245
 
# File 'lib/gs2crmod/gsl_data.rb', line 242

def transient_es_heat_flux_amplification_over_ky_gsl_vector(options)
  options[:direction] = :ky
  transient_es_heat_flux_amplification_over_kxy_gsl_vector(options)
end

#vres_energy_gsl_vector(options) ⇒ Object

Velocity space diagnostics: integral error due to energy resolution



899
900
901
902
903
 
# File 'lib/gs2crmod/gsl_data.rb', line 899

def vres_energy_gsl_vector(options)
  raise "Velocity space vres diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.vres"
  vres = GSL::Vector.filescan("#@directory/#@run_name.vres")
  return vres[2]
end

#vres_pitch_angle_gsl_vector(options) ⇒ Object

Velocity space diagnostics: integral error due to pitch angle resolution



891
892
893
894
895
 
# File 'lib/gs2crmod/gsl_data.rb', line 891

def vres_pitch_angle_gsl_vector(options)
  raise "Velocity space vres diagnostics not found" unless FileTest.exist? "#@directory/#@run_name.vres"
  vres = GSL::Vector.filescan("#@directory/#@run_name.vres")
  return vres[1]
end

#x_gsl_vector(options) ⇒ Object 



1004
1005
1006
1007
1008
1009
1010
1011
 
# File 'lib/gs2crmod/gsl_data.rb', line 1004

def x_gsl_vector(options)
  raise "options nakx and interpolate_x are incompatible" if options[:nakx] and options[:interpolate_x]
  kx = gsl_vector(:kx, options)
  lx = 2*Math::PI/kx.to_box_order[1]
  #ep 'lx', lx
  nx = options[:nakx]||kx.size
  GSL::Vector.indgen(nx, 0, lx/nx)
end

#y_gsl_vector(options) ⇒ Object 



1013
1014
1015
1016
1017
1018
1019
1020
 
# File 'lib/gs2crmod/gsl_data.rb', line 1013

def y_gsl_vector(options)
  raise "options naky and interpolate_y are incompatible" if options[:naky] and options[:interpolate_y]
  ky = gsl_vector(:ky, options)
  ly = 2*Math::PI/ky[1]
  ny = options[:naky]||ky.size
  ysize = ny*2-2+ny%2
  GSL::Vector.indgen(ysize, 0, ly/ysize)
end

#zf_velocity_over_x_gsl_vector(options) ⇒ Object

This function returns the zonal flow velocity as a function of x (the radial coordinate). This is v_ZF = kxfac*IFT(i k_x phi_imag), where kxfac = (qinp/rhoc)*grho(rhoc).



1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
 
# File 'lib/gs2crmod/gsl_data.rb', line 1030

def zf_velocity_over_x_gsl_vector(options)
  Dir.chdir(@directory) do
    raise CRFatal.new("Need either qinp or pk and epsl specified in order to calculate kxfac.
                      If using numerical equil use the option :kxfac to override calculation.") unless @qinp or (@pk and @eps) or options[:kxfac]

    kx = gsl_vector(:kx).to_box_order
    drhodpsi = gsl_vector('drhodpsi')

    phi = GSL::Vector.alloc(kx.size)
    for it in 0...gsl_vector(:t).size
      options[:t_index] = it
      phi += gsl_vector_complex('phi_zonal', options)
    end
    phi /= gsl_vector(:t).size

    if @qinp
      kxfac = (@qinp/@rhoc)/drhodpsi
    elsif @pk and @epsl
      kxfac = (@epsl/@pk)/drhodpsi
    elsif options[:kxfac]
      kxfac = options[:kxfac]
    end

    vec_zf_vel = GSL::Vector.alloc(kx.size)
    #Take imaginary part since i k_x will lead to imaginary part being real
    vec_zf_vel = 0.5*kxfac*(phi*kx).backward.imag
    return vec_zf_vel
  end
end

#zonal_spectrum_gsl_vector(options) ⇒ Object 



1075
1076
1077
1078
1079
1080
1081
1082
 
# File 'lib/gs2crmod/gsl_data.rb', line 1075

def zonal_spectrum_gsl_vector(options)
  Dir.chdir(@directory) do
    gmzf = gsl_matrix('spectrum_over_ky_over_kx',options)
    veczf = GSL::Vector.alloc(gmzf.shape[1])
    gmzf.shape[1].times{|i| veczf[i] = gmzf[0,i]}
    return veczf
  end
end