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.



1025
1026
1027
1028
 
# File 'lib/gs2crmod/gsl_data.rb', line 1025

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 



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

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.



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
480
481
 
# File 'lib/gs2crmod/gsl_data.rb', line 451

def es_heat_flux_over_kxy_gsl_vector(options)
  raise "The write_fluxes_by_mode flag was not set." unless @write_fluxes_by_mode
  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 



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

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 



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

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 



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

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 



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

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 



861
862
863
864
865
866
 
# File 'lib/gs2crmod/gsl_data.rb', line 861

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 



868
869
870
871
872
873
 
# File 'lib/gs2crmod/gsl_data.rb', line 868

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 



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

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 



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

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 



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
648
649
 
# File 'lib/gs2crmod/gsl_data.rb', line 597

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



885
886
887
888
889
 
# File 'lib/gs2crmod/gsl_data.rb', line 885

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



877
878
879
880
881
 
# File 'lib/gs2crmod/gsl_data.rb', line 877

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.



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

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 



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

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 



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

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 



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

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 



483
484
485
486
 
# File 'lib/gs2crmod/gsl_data.rb', line 483

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 



488
489
490
491
 
# File 'lib/gs2crmod/gsl_data.rb', line 488

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 



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

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 



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
741
742
 
# File 'lib/gs2crmod/gsl_data.rb', line 715

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 



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

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 



939
940
941
 
# File 'lib/gs2crmod/gsl_data.rb', line 939

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 



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
687
688
 
# File 'lib/gs2crmod/gsl_data.rb', line 651

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 



948
949
950
951
 
# File 'lib/gs2crmod/gsl_data.rb', line 948

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 



943
944
945
946
 
# File 'lib/gs2crmod/gsl_data.rb', line 943

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



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

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 



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

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 



958
959
960
961
 
# File 'lib/gs2crmod/gsl_data.rb', line 958

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 



953
954
955
956
 
# File 'lib/gs2crmod/gsl_data.rb', line 953

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 



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
801
802
 
# File 'lib/gs2crmod/gsl_data.rb', line 744

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 



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

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 



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

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



901
902
903
904
905
 
# File 'lib/gs2crmod/gsl_data.rb', line 901

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



893
894
895
896
897
 
# File 'lib/gs2crmod/gsl_data.rb', line 893

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 



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

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 



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

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).



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
1059
1060
 
# File 'lib/gs2crmod/gsl_data.rb', line 1032

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 



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

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