Class: CodeRunner::Gs2
- Inherits:
-
Run::FortranNamelist
- Object
- Run::FortranNamelist
- CodeRunner::Gs2
- Includes:
- FixNormOption, GSLComplexTensors, GSLMatrices, GSLTensors, GSLVectorComplexes, GSLVectors, ReadNetcdf, GraphKits
- Defined in:
- lib/gs2crmod/gs2.rb,
lib/gs2crmod/ingen.rb,
lib/gs2crmod/graphs.rb,
lib/gs2crmod/test_gs2.rb,
lib/gs2crmod/properties.rb,
lib/gs2crmod/gsl_data_3d.rb,
lib/gs2crmod/read_netcdf.rb,
lib/gs2crmod/calculations.rb,
lib/gs2crmod/check_convergence.rb,
lib/gs2crmod/gsl_data.rb
Overview
This module reads data from the new diagnostics output file <run_name>.cdf.
It is intended to replace a lot of the function of gsl_data.rb which reads the old netcdf file. In particular, it defines a new generic reader function which can read any variable in the new netcdf file using a standard set of index constraints
Direct Known Subclasses
Defined Under Namespace
Modules: FixNormOption, GSLComplexTensors, GSLMatrices, GSLTensors, GSLVectorComplexes, GSLVectors, ReadNetcdf, TestGs2 Classes: Astrogk, GraphKits, InputFileError, ListSubmitter, NetcdfSmartReader, OldNetcdfSmartReader, Phi, Spectrogk
Constant Summary collapse
- GS2_CRMOD_VERSION =
GS2_CRMOD_VERSION = Version.new(Gem.loaded_specs.version.to_s)
Version.new('0.5.0')
- CODE_SCRIPT_FOLDER =
MODULE_FOLDER = File.dirname(File.(__FILE__))
- NaN =
GSL::NAN
- SPECIES_DEPENDENT_NAMELISTS =
eval(File.read(folder + '/species_dependent_namelists.rb'), binding, folder + '/species_dependent_namelists.rb')
- SPECIES_DEPENDENT_VARIABLES_WITH_HELP =
SPECIES_DEPENDENT_NAMELISTS.values.inject({}) do |hash, namelist_hash| namelist_hash[:variables].each do |var, var_hash| hash[var] = var_hash[:help] end hash end
- SPECIES_DEPENDENT_VARIABLES =
SPECIES_DEPENDENT_VARIABLES_WITH_HELP.keys
- MAX_NAME_SIZE =
310
- AxisKit =
one day someone should get rid of this!
GraphKit::AxisKit
- DataKit =
GraphKit::DataKit
- GRAPHKIT_OPTIONS_HELP =
{ t_index_window: "[begin, end], window of time indices to plot (e.g. t_index_window: [0,10])", t_index: "integer, index of time at which to plot (e.g. t_index: 20)", t: "float, value of time at which to plot (e.g. t: 2.45)", ky_index: "integer, index of ky at which to plot (e.g. ky_index: 20)", ky: "float, value of ky at which to plot (e.g. ky: 0.1)", kx_index: "integer, index of kx at which to plot (e.g. kx_index: 20)", kx: "float, value of kx at which to plot (e.g. kx: 0.1)", with: "Gnuplot Option (may not apply when using other packages), e.g. with: 'lp' or with 'pm3d palette'", rgbformulae: "Gnuplot Option (may not apply when using other packages), sets colour mapping. See gnuplot help set rgbformulae", limit: "Limit the range of quantity begin plotted - any values of the quantity outside the limits will be set to the limit: eg. limit: [0,80]", flip: 'Flip the y axis, e.g. flip: true', rev: 'Reverse the x axis, e.g. rev: true', z: 'Plot quantities vs z = theta/shat rather than theta. See Beer, Cowley Hammet 1996, eg. z: true', norm: 'Normalise the graph so that its maximum is 1, e.g. norm: true', mag: 'Plot the magnitude, e.g. mag: true', species_index: "Which GS2 species to plot the graph for (1-based).", strongest_non_zonal_mode: "Plot the graph requested for the mode with the highest value of phi^2. Overrides ky, kx, ky_index, kx_index. Can be set true or false; e.g. strongest_non_zonal_mode: true", no_zonal: "Don't plot the ky=0 part (boolean, e.g. no_zonal: true)", no_kpar0: "Don't plot the kpar=0 part (boolean, e.g. no_kpar0: true)", log: "Plot the log of a given quantity (exact meaning varies). boolean", Rmaj: "The major radius in metres. This has no effect on the shape of the graph: it merely multiplies every length", n0: " The toroidal mode number of the longest y mode. In effect it is the number of periodic copies of the flux tube that will fit in the torus. Periodicity requires that n0 q is also an integer. If you specify :n0 where this is not the case, q will automatically be adjusted until it is", rho_star: " The ratio of the reference Lamour radius to the GS2 normalising length a. Cannot be specified at the same time as n0. If specified, both n0 and q will be adjusted to ensure periodicity", t_index: "The (1-based) time index", nakx: "The number of radial wave numbers to include in the plot. In effect, it is a low pass filter which reduces the resolution in the radial direction without changing the shape of the final surface. Minimum value is 4", naky: "The number of kys to include in the plot. In effect, it is a low pass filter which reduces the resolution in the y direction without changing the shape of the final surface. Minimum value is 4", gs2_coordinate_factor: "When set to 1, plot the graph in GS2 coordinates. When set to 0 plot the graph in real space. Can be set at any value between 0 and 1: the graph will smoothly distort between the two limits", xmax: "The (0-based) index of the maximum value of x to include in the plot", xmin: "The (0-based) index of the minimum value of x to include in the plot", ymax: "The (0-based) index of the maximum value of y to include in the plot", ymin: "The (0-based) index of the minimum value of y to include in the plot", thetamax: "The (0-based) index of the maximum value of theta to include in the plot", thetamin: "The (0-based) index of the minimum value of theta to include in the plot", theta_index: "integer, index of theta at which to plot (e.g. theta_index: 20)", kxfac: "float, overrides calculation of kxfac in zonal flow velocity function", ncopies: " The number of periodic copies of the flux tube to include", torphi_values: "An array of two values of the toroidal angle. The graph will be plotted in between those two values with poloidal cross sections at either end", magnify: " The magnification factor of the small section. It can take any value greater than or equal to 1", }
Constants included from GSLTensors
GSLTensors::FIELD_VALUES, GSLTensors::IRRELEVANT_INDICES, GSLTensors::TIME_VARYING_INDICES, GSLTensors::TRIVIAL_INDICES
Instance Attribute Summary collapse
-
#eigenfunctions ⇒ Object
Returns the value of attribute eigenfunctions.
-
#iphi00 ⇒ Object
Necessary for back.
-
#ky_graphs ⇒ Object
Returns the value of attribute ky_graphs.
-
#ky_list ⇒ Object
Returns the value of attribute ky_list.
-
#saturation_time ⇒ Object
Necessary for back.
-
#scan_index_window ⇒ Object
Returns the value of attribute scan_index_window.
-
#scan_parameter_value ⇒ Object
Returns the value of attribute scan_parameter_value.
-
#t_list ⇒ Object
Returns the value of attribute t_list.
-
#theta_list ⇒ Object
Returns the value of attribute theta_list.
Class Method Summary collapse
- .add_variable_to_namelist(namelist, var, value) ⇒ Object
- .cache ⇒ Object
- .check_and_update ⇒ Object
- .defaults_file_header ⇒ Object
- .generate_graphs_rdoc_file ⇒ Object
- .help_graphs ⇒ Object
- .list_code_commands ⇒ Object
- .modify_job_script(runner, runs_in, script) ⇒ Object
-
.test_gs2(*args) ⇒ Object
See TestGs2.
Instance Method Summary collapse
- #actual_number_of_processors ⇒ Object (also: #anop)
- #agk? ⇒ Boolean
- #approximate_grid_size ⇒ Object (also: #agridsze)
- #auto_axiskits(name, options) ⇒ Object
- #axiskit(name, options = {}) ⇒ Object
-
#bes_output(options = {}) ⇒ Object
This function will interpolate and output either phi or density at the outboard midplane on a 40x40 grid appropriate to analyse as experimental data.
- #box_kx_index(physical_kx_index) ⇒ Object
-
#calculate_frequencies ⇒ Object
Actually, this doesn’t calculate the frequencies but reads them from run_name.out.
- #calculate_growth_rate(vector, options = {}) ⇒ Object
- #calculate_growth_rates_and_frequencies ⇒ Object (also: #cgrf)
- #calculate_results ⇒ Object
-
#calculate_saturation_time_index(show_graph = false) ⇒ Object
(also: #csti)
I.e.
- #calculate_spectral_checks ⇒ Object (also: #csc)
- #calculate_time_averaged_fluxes ⇒ Object (also: #ctaf)
- #calculate_transient_amplification(vector, options = {}) ⇒ Object
- #calculate_transient_amplifications ⇒ Object (also: #cta)
- #calculate_transient_es_heat_flux_amplifications ⇒ Object (also: #ctehfa)
- #calculate_vspace_checks ⇒ Object (also: #cvc)
- #check_converged ⇒ Object
-
#check_parameters ⇒ Object
Eventually, this will be a full port of the ingen tool in the GS2 folder.
- #code_run_environment ⇒ Object
-
#corrected_mom_flux_stav ⇒ Object
Not needed for GS2 built after 16/06/2010.
-
#correlation_analysis(options = {}) ⇒ Object
This function will handle running the correlation analysis and writing the results to a NetCDF file.
- #ctan ⇒ Object
- #cumulative_gridpoints ⇒ Object
- #data_string ⇒ Object
-
#delete_restart_files(options = {}) ⇒ Object
Delete all the restart files (irreversible!).
- #diagnostics_namelist ⇒ Object
- #error(message) ⇒ Object
-
#estimated_nodes ⇒ Object
(also: #estnod)
Gives a guess as to the maximum number of nodes which can be can be utilized on the current system.
-
#eulerian_kx_index(options) ⇒ Object
This function is used in the presence of perpendicular flow shear.
- #generate_component_runs ⇒ Object
- #generate_input_file(&block) ⇒ Object
- #get_completed_timesteps ⇒ Object
- #get_list_of(*args) ⇒ Object (also: #list)
-
#get_run_time ⇒ Object
Try to read the runtime in minutes from the GS2 standard out.
- #get_status ⇒ Object
- #get_time ⇒ Object
- #graphkit(name, options = {}) ⇒ Object
-
#gridpoints ⇒ Object
A hash which gives the actual numbers of gridpoints indexed by their corresponding letters in the layout string.
- #gryfx? ⇒ Boolean
- #gsl_complex(name, options = {}) ⇒ Object
- #gsl_matrix(name, options = {}) ⇒ Object
- #gsl_tensor(name, options) ⇒ Object
- #gsl_vector(name, options = {}) ⇒ Object
- #gsl_vector_complex(name, options = {}) ⇒ Object
- #has_electrons? ⇒ Boolean
- #hypercoll_graphkit(options) ⇒ Object
- #hyperviscosity_graphkit(options) ⇒ Object
- #incomplete ⇒ Object
-
#ingen ⇒ Object
Run the ingen tool on the input file.
- #input_file_extension ⇒ Object
- #input_file_header ⇒ Object
- #jump(options) ⇒ Object
- #kx_indexed ⇒ Object
- #kx_shift(options) ⇒ Object
-
#latex_graphs ⇒ Object
This section defines a selection of graphs which are written to a latex file when the CR function write_report is called.
- #lenardbern_graphkit(options) ⇒ Object
-
#list_of_restart_files ⇒ Object
(also: #lorf)
Return a list of restart file paths (relative to the run directory).
- #max_es_heat_amp(species_index) ⇒ Object
-
#max_nprocs_no_x ⇒ Object
ep parallelisation.
- #max_trans_phi ⇒ Object
- #namelist_test_failed(namelist, tst) ⇒ Object
- #ncclose ⇒ Object
- #ncdump(names = nil, values = nil, extension = '.out.nc') ⇒ Object
- #netcdf_file ⇒ Object
- #netcdf_filename ⇒ Object
- #netcdf_smart_reader ⇒ Object
-
#no_restarts ⇒ Object
Returns true if this run has not been restarted, false if it has.
- #old_smart_graphkit(options) ⇒ Object
- #optimisation_flags ⇒ Object
-
#parallelizable_meshpoints ⇒ Object
Gives a guess as to the maximum number of meshpoints which can be parallelized (i.e. excluding ntheta).
- #parameter_string ⇒ Object
- #parameter_transition(run) ⇒ Object
- #percent_complete ⇒ Object
- #physical_kx_index(box_kx_index) ⇒ Object
- #plot_efit_file ⇒ Object
- #print_out_line ⇒ Object
-
#process_directory_code_specific ⇒ Object
This method, as its name suggests, is called whenever CodeRunner is asked to analyse a run directory.this happens if the run status is not :Complete, or if the user has specified recalc_all(-A on the command line) or reprocess_all (-a on the command line).
-
#recheck ⇒ Object
class ListSubmitter.
- #renew_info_file ⇒ Object
- #restart(new_run) ⇒ Object
- #restart_chain ⇒ Object
-
#run_heuristic_analysis ⇒ Object
This method overrides a method defined in heuristic_run_methods.rb in the CodeRunner source.
- #run_namelist_backwards_compatibility ⇒ Object
-
#run_namelist_tests(namelist, hash, enum = nil) ⇒ Object
Checks input parameters for inconsistencies and prints a report.
- #saturated_time_average(name, options) ⇒ Object
- #saturated_time_average_error(name, options) ⇒ Object
- #saturated_time_average_std_dev(name, options) ⇒ Object
- #sc(min) ⇒ Object
- #set_nprocs ⇒ Object
- #show_graph ⇒ Object
- #smart_graphkit(options) ⇒ Object
- #spec_chec(min, *dirns) ⇒ Object
- #species_letter ⇒ Object
- #species_type(index) ⇒ Object
- #spectrogk? ⇒ Boolean
-
#standardize_restart_files ⇒ Object
Put restart files in the conventional location, i.e.
- #stop ⇒ Object
- #test_failed(namelist, var, gs2_var, tst) ⇒ Object
- #test_variable(namelist, var, var_hash, ext, value) ⇒ Object
- #update_physics_parameters_from_miller_input_file(file) ⇒ Object
- #vim_output ⇒ Object (also: #vo)
- #vim_stdout ⇒ Object (also: #vo1)
- #visually_check_growth_rate(ky = nil) ⇒ Object
- #warning(message) ⇒ Object
- #write_input_file ⇒ Object
Methods included from GSLMatrices
#es_heat_flux_over_ky_over_kx_gsl_matrix, #growth_rate_over_ky_over_kx_gsl_matrix, #phi0_over_x_over_y_gsl_matrix, #spectrum_over_ky_over_kpar_gsl_matrix, #spectrum_over_ky_over_kx_gsl_matrix, #transient_amplification_over_ky_over_kx_gsl_matrix
Methods included from GSLVectorComplexes
#phi_along_field_line_gsl_vector_complex, #phi_zonal_gsl_vector_complex
Methods included from GSLVectors
#apar2_over_time_gsl_vector, #dt_gsl_vector, #es_heat_by_kx_over_time_gsl_vector, #es_heat_by_ky_over_time_gsl_vector, #es_heat_flux_over_time_gsl_vector, #es_heat_over_kx_gsl_vector, #es_heat_over_kxy_gsl_vector, #es_heat_over_ky_gsl_vector, #es_heat_par_over_time_gsl_vector, #es_heat_perp_over_time_gsl_vector, #es_mom_flux_over_time_gsl_vector, #frequency_by_kx_over_time_gsl_vector, #frequency_by_kxy_over_time_gsl_vector, #frequency_by_ky_over_time_gsl_vector, #frequency_over_ky_gsl_vector, #grho_gsl_vector, #growth_rate_by_kx_over_time_gsl_vector, #growth_rate_by_kxy_over_time_gsl_vector, #growth_rate_by_ky_over_time_gsl_vector, #growth_rate_over_kx_gsl_vector, #growth_rate_over_kx_slice_gsl_vector, #growth_rate_over_ky_gsl_vector, #growth_rate_over_ky_slice_gsl_vector, #hflux_tot_over_time_gsl_vector, #kpar_gsl_vector, #linked_kx_elements_gsl_vector, #lpc_energy_gsl_vector, #lpc_pitch_angle_gsl_vector, #mean_flow_velocity_over_x_gsl_vector, #par_mom_flux_over_time_gsl_vector, #perp_mom_flux_over_time_gsl_vector, #phi0_by_kx_by_ky_over_time_gsl_vector, #phi2_by_kx_over_time_gsl_vector, #phi2_by_ky_over_time_gsl_vector, #phi2_by_mode_over_time_gsl_vector, #phi2tot_over_time_gsl_vector, #phi_along_field_line_gsl_vector, #phi_for_eab_movie_gsl_vector, #scan_parameter_value_gsl_vector, #spectrum_over_kpar_gsl_vector, #spectrum_over_kx_gsl_vector, #spectrum_over_kxy_gsl_vector, #spectrum_over_ky_gsl_vector, #theta_along_field_line_gsl_vector, #tpar2_by_mode_over_time_gsl_vector, #tperp2_by_mode_over_time_gsl_vector, #transient_amplification_over_kx_gsl_vector, #transient_amplification_over_ky_gsl_vector, #transient_es_heat_flux_amplification_over_kx_gsl_vector, #transient_es_heat_flux_amplification_over_kxy_gsl_vector, #transient_es_heat_flux_amplification_over_ky_gsl_vector, #vres_energy_gsl_vector, #vres_pitch_angle_gsl_vector, #x_gsl_vector, #y_gsl_vector, #zonal_flow_velocity_over_x_gsl_vector, #zonal_spectrum_gsl_vector
Methods included from FixNormOption
#fix_heat_flux_norm, #fix_norm, #fix_norm_action
Methods included from ReadNetcdf
#new_ncclose, #new_netcdf_file, #new_netcdf_filename
Methods included from GSLComplexTensors
#field_gsl_tensor_complex, #phi_gsl_tensor_complex
Methods included from GSLTensors
#apar_gsl_tensor, #bpar_gsl_tensor, #cartesian_coordinates_gsl_tensor, #constant_torphi_surface_gsl_tensor, #correct_3d_options, #cylindrical_coordinates_gsl_tensor, #field_gsl_tensor, #field_netcdf_name, #field_real_space_gsl_tensor, #field_real_space_gsl_tensor_2, #field_species_element, #geometric_factors_gsl_tensor, #moment_gsl_tensor, #phi_real_space_gsl_tensor
Instance Attribute Details
#eigenfunctions ⇒ Object
Returns the value of attribute eigenfunctions.
410 411 412 |
# File 'lib/gs2crmod/gs2.rb', line 410 def eigenfunctions @eigenfunctions end |
#iphi00 ⇒ Object
Necessary for back. comp. due to an old bug
1056 1057 1058 |
# File 'lib/gs2crmod/gs2.rb', line 1056 def iphi00 @iphi00 end |
#ky_graphs ⇒ Object
Returns the value of attribute ky_graphs.
410 411 412 |
# File 'lib/gs2crmod/gs2.rb', line 410 def ky_graphs @ky_graphs end |
#ky_list ⇒ Object
Returns the value of attribute ky_list.
410 411 412 |
# File 'lib/gs2crmod/gs2.rb', line 410 def ky_list @ky_list end |
#saturation_time ⇒ Object
Necessary for back. comp. due to an old bug
1056 1057 1058 |
# File 'lib/gs2crmod/gs2.rb', line 1056 def saturation_time @saturation_time end |
#scan_index_window ⇒ Object
Returns the value of attribute scan_index_window.
411 412 413 |
# File 'lib/gs2crmod/gs2.rb', line 411 def scan_index_window @scan_index_window end |
#scan_parameter_value ⇒ Object
Returns the value of attribute scan_parameter_value.
411 412 413 |
# File 'lib/gs2crmod/gs2.rb', line 411 def scan_parameter_value @scan_parameter_value end |
#t_list ⇒ Object
Returns the value of attribute t_list.
410 411 412 |
# File 'lib/gs2crmod/gs2.rb', line 410 def t_list @t_list end |
#theta_list ⇒ Object
Returns the value of attribute theta_list.
410 411 412 |
# File 'lib/gs2crmod/gs2.rb', line 410 def theta_list @theta_list end |
Class Method Details
.add_variable_to_namelist(namelist, var, value) ⇒ Object
914 915 916 917 |
# File 'lib/gs2crmod/gs2.rb', line 914 def self.add_variable_to_namelist(namelist, var, value) var = :stir_ + var if namelist == :stir super(namelist, var, value) end |
.cache ⇒ Object
101 102 103 104 |
# File 'lib/gs2crmod/graphs.rb', line 101 def self.cache @cache ||= {} @cache end |
.check_and_update ⇒ Object
415 416 417 418 |
# File 'lib/gs2crmod/gs2.rb', line 415 def check_and_update old_check_and_update @readout_list = (@variables + @results - [:growth_rates_by_ky, :growth_rates, :real_frequencies, :real_frequencies_by_ky, :ky_list, :kx_list, :theta_list, :t_list]) end |
.defaults_file_header ⇒ Object
944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 |
# File 'lib/gs2crmod/gs2.rb', line 944 def self.defaults_file_header <<EOF1 ###################################################################### # Automatically generated defaults file for GS2 CodeRunner module # # # # This defaults file specifies a set of defaults for GS2 which are # # used by CodeRunner to set up and run GS2 simulations. # # # # Created #{Time.now.to_s} # # # ###################################################################### @defaults_file_description = "" EOF1 end |
.generate_graphs_rdoc_file ⇒ Object
106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 |
# File 'lib/gs2crmod/graphs.rb', line 106 def self.generate_graphs_rdoc_file File.open('graphs_rdoc.rb', 'w') do |file| graphs = self.instance_methods.find_all{|m| m.to_s =~ /_graphkit$/}.sort_by{|m| m.to_s} run = new(nil) file.puts "class #{self.to_s}::GraphKits\n" graphs.each do |graph| help = run.send(graph, command: :help) = run.send(graph, command: :options) file.puts "# #{help}" if and .size > 0 file.puts "# Options:" .each do |op| file.puts "#\n# #{op}: #{GRAPHKIT_OPTIONS_HELP[op]}" end end file.puts "def #{graph}\nend" end file.puts "end" end end |
.help_graphs ⇒ Object
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 |
# File 'lib/gs2crmod/graphs.rb', line 126 def self.help_graphs # @@runner ||= CodeRunner.fetch_runner(U: true, string = "" graphs = self.instance_methods.find_all{|m| m.to_s =~ /_graphkit$/}.sort_by{|m| m.to_s} run = new(nil) string << "-------------------------------------------\n Available Graphs For #{self.to_s}\n-------------------------------------------\n\n" graphs.each do |graph| help = run.send(graph, command: :help) = run.send(graph, command: :options) string << "\n------------------------------------\n#{graph.to_s.sub(/_graphkit/, '')}\n------------------------------------\n\n#{help}\n" if and .size > 0 string << "\n\tOptions:\n" .each do |op| string << "\t\t#{op}: #{GRAPHKIT_OPTIONS_HELP[op]}\n" end end end string.paginate end |
.list_code_commands ⇒ Object
910 911 912 |
# File 'lib/gs2crmod/gs2.rb', line 910 def self.list_code_commands puts (methods - Run.methods).sort end |
.modify_job_script(runner, runs_in, script) ⇒ Object
723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 |
# File 'lib/gs2crmod/gs2.rb', line 723 def self.modify_job_script(runner, runs_in, script) if CODE_OPTIONS[:gs2] and CODE_OPTIONS[:gs2][:list] if (list_size = CODE_OPTIONS[:gs2][:list]).kind_of? Integer raise "The total number of runs must be a multiple of the list size!" unless runs_in.size % list_size == 0 pieces = runs_in.pieces(runs_in.size/list_size) else pieces = [runs_in] end script = "" pieces.each do |runs| #ep 'there is a list' FileUtils.makedirs('job_lists') jid = "#{runs[0].id}-#{runs[-1].id}" list_file = "job_lists/gs2_list_#{jid}.list" File.open(list_file,'w') do |file| file.puts runs.size file.puts runs.map{|r| "#{r.relative_directory}/#{r.run_name}"}.join("\n") end raise "runs must all have the same nprocs" unless runs.map{|r| r.nprocs}.uniq.size == 1 runs.each do |r| # Make sure the restart file name includes the relative directory for # list runs reldir = r.relative_directory rdir = r.restart_dir #puts rdir[0...reldir.size] == reldir, rdir[0...reldir.size], reldir #raise "" if rdir r.restart_dir = reldir + '/' + rdir if not rdir[0...reldir.size] == reldir else r.restart_dir = reldir end Dir.chdir(r.directory){r.write_input_file} end np = runs[0].nprocs.split('x').map{|n| n.to_i} np[0] *= runs.size nprocs = np.map{|n| n.to_s}.join('x') @runner.nprocs = nprocs ls = ListSubmitter.new(@runner, nprocs, list_file, jid) script << ls.run_command end end return script end |
Instance Method Details
#actual_number_of_processors ⇒ Object Also known as: anop
868 869 870 871 |
# File 'lib/gs2crmod/gs2.rb', line 868 def actual_number_of_processors raise "Please specify the processor layout using the -n or (n:) option" unless @nprocs @nprocs.split('x').map{|n| n.to_i}.inject(1){|ntot, n| ntot*n} end |
#agk? ⇒ Boolean
59 60 61 |
# File 'lib/gs2crmod/gs2.rb', line 59 def agk? false end |
#approximate_grid_size ⇒ Object Also known as: agridsze
875 876 877 878 879 880 881 882 |
# File 'lib/gs2crmod/gs2.rb', line 875 def approximate_grid_size case @grid_option when "box" (2*(@nx-1)/3+1).to_i * (@naky||(@ny-1)/3+1).to_i * @ntheta * (2 * @ngauss + @ntheta/2).to_i * @negrid * 2 * @nspec else @ntheta * (2 * @ngauss + @ntheta/2).to_i * @negrid * 2 * @nspec end end |
#auto_axiskits(name, options) ⇒ Object
11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 |
# File 'lib/gs2crmod/graphs.rb', line 11 def auto_axiskits(name, ) hash = cache[:auto_axiskits] ||= {'t' => ['Time', ''], 'phi2tot_over_time' => ['Phi^2 Total', ''], 'apar2_over_time' => ['Apar^2 Total', ''], 'growth_rate_by_ky_over_time' => ['Growth Rate by ky', ''], 'growth_rate_by_kx_over_time' => ['Growth Rate by kx', ''], 'growth_rate_by_mode_over_time' => ["Growth Rate by mode", ''], # <MJL additions 2013-09-19> 'frequency_by_ky_over_time' => ['Real frequency by ky', ''], 'frequency_by_kx_over_time' => ['Real frequency by kx', ''], # </MJL> 'phi2_by_ky_over_time' => ['Phi^2 by ky', ''], 'phi2_by_kx_over_time' => ['Phi^2 by ky', ''], 'es_heat_by_ky_over_time' => ['Phi^2 by ky', ''], 'es_heat_by_kx_over_time' => ['Phi^2 by kx', ''], 'phi2_by_mode_over_time' => ["Phi^2 by mode", ''], 'tpar2_by_mode_over_time' => ["(delta T_parallel)^2 by mode", '%'], 'tperp2_by_mode_over_time' => ["(delta T_perp)^2 by mode", '%'], 'hflux_tot' => ['Total Heat Flux', ''], 'es_heat_par' => ['Parallel electrostatic heat flux', ''], 'es_heat_perp' => ['Perpendicular electrostatic heat flux', ''], 'ky' => ['ky', "1/rho_#{species_letter}"], 'kx' => ['kx', "1/rho_#{species_letter}"], 'x' => ['x', "rho_#{species_letter}", 1], 'kpar' => ['kpar', "2 pi/qR"], 'growth_rate_over_kx' => ['Growth Rate', "v_th#{species_letter}/a", 1], 'growth_rate_over_ky' => ['Growth Rate', "v_th#{species_letter}/a", 1], 'growth_rate_over_kx_slice' => ['Growth Rate', "v_th#{species_letter}/a", 1], 'growth_rate_over_ky_slice' => ['Growth Rate', "v_th#{species_letter}/a", 1], 'growth_rate_over_ky_over_kx' => ["Growth Rate", "v_th#{species_letter}/a", 2], 'frequency_over_ky' => ['Frequency', "v_th#{species_letter}/a", 1], 'transient_es_heat_flux_amplification_over_kx' => ['Transient Electrostatic Heat Amplification', "", 1], 'transient_es_heat_flux_amplification_over_ky' => ['Transient Electrostatic Heat Amplification', "", 1], 'transient_amplification_over_kx' => ['Transient Amplification', "", 1], 'transient_amplification_over_ky' => ['Transient Amplification', "", 1], 'spectrum_over_kx' => ["Spectrum at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", '', 1], 'spectrum_over_ky' => ["Spectrum at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", '', 1], 'es_heat_over_kx' => ["Heat Flux at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", 'Q_gB', 1], 'es_heat_over_ky' => ["Heat Flux at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", 'Q_gB', 1], 'es_heat_flux_over_ky_over_kx' => ["Heat flux at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", '', 2], 'spectrum_over_kpar' => ["Spectrum at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", '', 1], 'spectrum_over_ky_over_kx' => ["Spectrum at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", '', 2], 'spectrum_over_ky_over_kpar' => ["Spectrum at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", '', 2], #'phi0_over_x_over_y' => ["Phi at t = #{sprintf("%.3f" ,(options[:t] or list(:t)[options[:t_index]] or list(:t).values.max))}", '', 2], 'phi0_over_x_over_y' => ["Phi at theta = 0", '', 2], 'es_mom_flux_over_time' => ["#{species_type(([:species_index] or 1)).capitalize} Momentum Flux", '', 1], 'zonal_spectrum' => ["Zonal spectrum at t = #{sprintf("%.3f" ,([:t] or list(:t)[[:t_index]] or list(:t).values.max))}", '', 1], 'zonal_flow_velocity_over_x' => ['Zonal Flow Velocity', "", 1], 'mean_flow_velocity_over_x' => ['Mean Flow Velocity', "", 1] } return hash[name] end |
#axiskit(name, options = {}) ⇒ Object
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 |
# File 'lib/gs2crmod/graphs.rb', line 64 def axiskit(name, ={}) logf :axiskit if info = auto_axiskits(name, ) if info[2] and info[2] == 2 axis = GraphKit::AxisKit.autocreate({data: gsl_matrix(name, ), title: info[0], units: info[1]}) elsif !info[2] or info[2] == 1 axis = GraphKit::AxisKit.autocreate({data: gsl_vector(name, ), title: info[0], units: info[1]}) log 'successfully created axis' end return axis end case name when 'phi_along_field_line' title = [:imrc].to_s.capitalize + " Phi" units = "" return GraphKit::AxisKit.autocreate(data: gsl_vector(name, ), title: title, units: units) when 'theta_along_field_line' title = [:z] ? "z/l_B" : 'Theta' units = [:z] ? '' : 'radians' return GraphKit::AxisKit.autocreate(data: gsl_vector(name, ), title: title, units: units) when 'es_heat_flux' type = species_type([:species_index]).capitalize units = '' return GraphKit::AxisKit.autocreate(data: gsl_vector('es_heat_flux_over_time', ), title: "#{type} Heat Flux", units: units) # when 'spectrum_by_ky' # return AxisKit.autocreate(data: gsl_vector('spectrum_by_ky', options), title: "Phi^2 at t = #{list(:t)[options[:t_index]]}", units: '') when 'es_heat_par' puts "heat par" type = species_type([:species_index]).capitalize units = '' return GraphKit::AxisKit.autocreate(data: gsl_vector('es_heat_par_over_time', ), title: "#{type} parallel es heat flux", units: units) # when 'spectrum_by_ky' # return AxisKit.autocreate(data: gsl_vector('spectrum_by_ky', options), title: "Phi^2 at t = #{list(:t)[options[:t_index]]}", units: '') end raise CRError.new("Unknown axis kit: #{name}") end |
#bes_output(options = {}) ⇒ Object
This function will interpolate and output either phi or density at the outboard midplane on a 40x40 grid appropriate to analyse as experimental data. It called as a run_command e.g. rc ‘bes_output(options)’, j:<run number>. It will call field_real_space_poloidal_plane_graphkit for every time step, interpolate at outboard midplane, and write fields and grids out to NetCDF file.
Options: Same as field_real_space_poloidal_plane, field name must also be specified for generality. New options:
no_flux_tube_copies: ensures only one flux tube is printed out with zeroes everywhere else. amin: Minor radius (to which R,Z are normalized) so that grid is in right units output_box_size: Array of sizes of output box (in units of amin) either side of middle of fluxtube at outboard midplane in R direction and either side of outboard midplane in Z direction. output_box_points: Array of number of points in output box (R,Z). Default will be 50x50.
The interpolation routine used will only interpolate correctly inside the fluxtube and produce garbage outside. Regular points are checked for being inside or outside the fluxtube and values of the field outside the fluxtube are set to zero.
1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 |
# File 'lib/gs2crmod/gs2.rb', line 1282 def bes_output(={}) #****************** # Read in options * #****************** #In order to interpolate on constant grids, ensure constant_torphi is set to some value (default 0.0) if [:constant_torphi] == nil p 'constant_torphi not set! Setting it to 0.0.' [:constant_torphi] = 0.0 end #Check whether t_index_window is specified, if not, set to entire t range if [:t_index_window] == nil t_index_beg = 1 t_index_end = gsl_vector(:t).length else t_index_beg = [:t_index_window][0] t_index_end = [:t_index_window][1] end if [:amin] amin = [:amin] end if [:v_ref] # velocity of reference species v_ref = [:v_ref] end if [:omega] # angular velocity of plasma omega = [:omega] end if [:omega] and ([:v_ref] == nil or [:amin] == nil) raise 'Need to specify amin AND v_ref options when specifying omega to move to LAB frame!' end if [:output_box_size] and [:output_box_size].kind_of?Array _r_box_size = [:output_box_size][0] # EGH These variables are marked as unused... are they used anywhere? _z_box_size = [:output_box_size][1] #else # raise 'Option output_box_size must be specified (in units of amin) and must be an Array.' end if [:output_box_points] and [:output_box_points].kind_of?Array r_box_pts = [:output_box_points][0] z_box_pts = [:output_box_points][1] else r_box_pts = 50 z_box_pts = 50 end #Call at first time step to set up arrays and grids [:t_index] = t_index_beg kit = field_real_space_poloidal_plane_graphkit() x = kit.data[0].x.data _y = kit.data[0].y.data z = kit.data[0].z.data #Set up NetCDf file file = NumRu::NetCDF.create(@run_name + "_bes_output.nc") xdim = file.def_dim('y', x.shape[0]) zdim = file.def_dim('z', z.shape[1]) tdim = file.def_dim('t', 0) #zero means unlimited x_var = file.def_var("x", 'sfloat', [xdim, zdim]) z_var = file.def_var("z", 'sfloat', [xdim, zdim]) t_var = file.def_var("t", 'sfloat', [tdim]) field_var = file.def_var("field", 'sfloat', [xdim, zdim, tdim]) file.enddef #Write dimensions to file x_var.put(NArray.to_na(x.to_a)) z_var.put(NArray.to_na(z.to_a)) #Loop over time, load field as function of space at each time index, write to file for i in t_index_beg...t_index_end #inclusive of end Terminal.erewind(1) #go back one line in terminal eputs sprintf("Writing time index = %d of %d#{Terminal::CLEAR_LINE}", i, t_index_end-t_index_beg+1) #clear line and print time index [:t_index] = i #Need to test whether omega is specified to change torphi at each time step. If not, do nothing since torphi must be #set to a value to call the graphkit below if [:omega] [:torphi] = omega * (gsl_vector(:t)[i] - gsl_vector(:t)[0]) * (amin/v_ref) end kit = field_real_space_poloidal_plane_graphkit() t_var.put(gsl_vector(:t)[i], 'index'=>[i-t_index_beg]) #Write time to unlimited time NetCDF variable field_var.put(NArray.to_na((kit.data[0].f.data).to_a), 'start'=>[0,0,i-t_index_beg], 'end'=>[-1,-1,i-t_index_beg]) end file.close #Ignore this until interpolation issue is sorted =begin #************************** # Set up new regular grid * #************************** th_grid_size = x.shape[1] flux_tube_midpt = x[0, (th_grid_size-1)/2] + (x[-1, (th_grid_size-1)/2] - x[1, (th_grid_size-1)/2])/2 x_vec_reg = GSL::Vector.linspace(flux_tube_midpt - r_box_size, flux_tube_midpt + r_box_size, r_box_pts) z_vec_reg = GSL::Vector.linspace(-z_box_size, z_box_size, z_box_pts) x_reg = GSL::Matrix.alloc(r_box_pts, z_box_pts) z_reg = GSL::Matrix.alloc(r_box_pts, z_box_pts) field_reg = GSL::Matrix.alloc(r_box_pts, z_box_pts) for i in 0...r_box_pts for j in 0...z_box_pts x_reg[i,j] = x_vec_reg[i] z_reg[i,j] = z_vec_reg[j] end end #************************************************ # Find the field at every point on regular grid * #************************************************ #To evaluate field on a regular grid given the field on an irregular grid, need to interpolate. The rubygem #gsl_extras contains an interpolation routine called ScatterInterp which does exactly this based on a #'Radial Basis Function' method. #Have R, Z, and field on an irregular grid in the form of matrices. ScatterInterp only takes in GSL vectors #so simply convert these matrices to vectors (of size row*col) since the order of the pts don't matter. x_vec = GSL::Vector.alloc(x.shape[0]*x.shape[1]) z_vec = GSL::Vector.alloc(x.shape[0]*x.shape[1]) field_vec = GSL::Vector.alloc(x.shape[0]*x.shape[1]) for i in 0...x.shape[0] for j in 0...x.shape[1] x_vec[x.shape[1]*i + j] = x[i,j] z_vec[x.shape[1]*i + j] = z[i,j] field_vec[x.shape[1]*i + j] = field[i,j] end end #Now pass these vectors to ScatterInterp. This creates an object with instance method 'eval' which can be given an x,z coord #at which to evaluate the interpolated function. p 'Interpolating' interp = GSL::ScatterInterp.alloc(:linear, [x_vec, z_vec, field_vec], false, r0=0.1) p 'Finished interpolating' for i in 0...x_vec_reg.size for j in 0...z_vec_reg.size field_reg[i,j] = interp.eval(x_vec_reg[i], z_vec_reg[j]) end end kit = GraphKit.quick_create([x_vec_reg, z_vec_reg, field_reg]) #kit2 = GraphKit.quick_create([x_vec, z_vec, field_vec]) =end end |
#box_kx_index(physical_kx_index) ⇒ Object
1467 1468 1469 1470 |
# File 'lib/gs2crmod/gsl_data.rb', line 1467 def box_kx_index(physical_kx_index) return kx_indexed[physical_kx_index] end |
#calculate_frequencies ⇒ Object
Actually, this doesn’t calculate the frequencies but reads them from run_name.out. Requires write_line to be .true.
267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 |
# File 'lib/gs2crmod/calculations.rb', line 267 def calculate_frequencies @real_frequencies = FloatHash.new gs2_out = FileUtils.tail(@run_name + ".out", list(:ky).size*list(:kx).size) # a = gs2_out.split("\n") final_timestep_list = gs2_out #a.slice((a.size-@ky_list.size*@kx_list.size-1)..a.size-1).join("\n") log(final_timestep_list.slice(-2..-1)) # eputs final_timestep_list f = LongRegexen::FLOAT.verbatim logi(f) @frequency_at_ky_at_kx||= FloatHash.new ky_values = [] regex = Regexp.new( "^.*aky=\\s*(?<aky>#{f})\s*akx=\\s*(?<akx>#{f}).*omav=\\s*(?<re>#{f})\\s*(?<gr>#{f})") final_timestep_list.scan(regex) do aky = eval($~[:aky]) akx = eval($~[:akx]) @frequency_at_ky_at_kx[aky] = FloatHash.new unless ky_values.include? aky ky_values.push aky @frequency_at_ky_at_kx[aky][akx] = eval($~[:re]) end end |
#calculate_growth_rate(vector, options = {}) ⇒ Object
390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 |
# File 'lib/gs2crmod/calculations.rb', line 390 def calculate_growth_rate(vector, ={}) raise "This vector should be positive definite" if vector.min < 0.0 offset = 0 length = vector.length while vector[offset] == 0.0 offset+=1 return 0.0 if offset == vector.length end growth_rate = GSL::Fit::linear(gsl_vector(:t).subvector(offset, length-offset), 0.5*GSL::Sf::log(vector.subvector(offset, length - offset)))[1] divisor = 1 while (growth_rate.to_s == "NaN") #This corrects the growth rate if phi has grown all the way to NaN during the simulation divisor *= 2 length = (vector.size.to_f / divisor.to_f).floor # p length return "NaN" if length <= offset + 1 growth_rate = GSL::Fit::linear(gsl_vector(:t).subvector(offset, length-offset), 0.5*GSL::Sf::log(vector.subvector(offset, length-offset)))[1] end growth_rate end |
#calculate_growth_rates_and_frequencies ⇒ Object Also known as: cgrf
287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 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 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 |
# File 'lib/gs2crmod/calculations.rb', line 287 def calculate_growth_rates_and_frequencies return if @grid_option == "single" and @aky == 0.0 # no meaningful results Dir.chdir(@directory) do logf(:calculate_growth_rates_and_frequencies) logd calculate_frequencies # get_list_of(:ky, :kx) @growth_rates= FloatHash.new #raise CRFatal.new("Unknown value of ky read from output file: #{data[:aky].to_f}. Not in list:\n#{list(:ky).values.inspect}") # pp @ky_list # With zero magnetic shear, calculate growth rates for both kx and ky #if @shat and @shat.abs < 1.0e-5 and @nx and @nx > 1 to_calc = [:kx, :ky] @growth_rate_at_kx ||= FloatHash.new #else #to_calc = [:ky] #end @growth_rate_at_ky ||= FloatHash.new eputs # p @growth_rate_at_kx; exit to_calc.each do |kxy| growth_rates = send(:growth_rate_at_ + kxy) list(kxy).values.sort.each do |value| #p growth_rates.keys, value, growth_rates[value.to_f-0.0], #growth_rates.class, growth_rates.keys.include?(value); exit next if growth_rates.keys.include? value Terminal.erewind(1) #ep growth_rates.keys eputs sprintf("Calculating growth rate for #{kxy} = % 1.5e#{Terminal::CLEAR_LINE}", value) # Mode has 0 growth rate at ky==0 (growth_rates[value] = 0.0; next) if value == 0.0 and kxy == :ky if @g_exb_start_timestep t_index_window = [1, [(g_exb_start_timestep-1)/@nwrite, list(:t).keys.max].min] #ep "t_index_window", t_index_window else t_index_window = nil end if list(kxy).size == 1 phi2_vec = gsl_vector("phi2tot_over_time", t_index_window: t_index_window) else phi2_vec = gsl_vector("phi2_by_#{kxy}_over_time", kxy=>value, :t_index_window=> t_index_window) end (growth_rates[value] = 0.0; next) if phi2_vec.min <= 0.0 growth_rates[value] = calculate_growth_rate(phi2_vec) (eputs "\n\n----------\nIn #@run_name:\n\nphi2_by_#{kxy}_over_time is all NaN; unable to calculate growth rate\n----------\n\n"; growth_rates[value] = -1; next) if growth_rates[value] == "NaN" end end write_results # ep "growth_rate_at_ky", @growth_rate_at_ky if ENV['GS2_CALCULATE_ALL'] trap(0){eputs "Calculation of spectrum did not complete: run 'cgrf' (i.e. calculate_growth_rates_and_frequencies) for this run. E.g. from the command line \n $ coderunner rc 'cgrf' -j #{@id}"; exit} @growth_rate_at_ky_at_kx ||= FloatHash.new list(:ky).values.sort.each do |kyv| # MJL 2013-11-07: The line below originally used ||= instead of =. I'm not sure why, since ||= does not seem to work. @growth_rate_at_ky_at_kx[kyv] = FloatHash.new list(:kx).values.sort.each do |kxv| # MJL 2013-11-07: I'm not sure why this next line was originally included. It seemed to cause almost all k's to be skipped. #next if @growth_rate_at_ky_at_kx[kyv].keys.include? kxv Terminal.erewind(1) eputs sprintf("Calculating growth rate for kx = % 1.5e and ky = % 1.5e#{Terminal::CLEAR_LINE}", kxv, kyv) (@growth_rate_at_ky_at_kx[kyv][kxv] = 0.0; next) if kyv == 0.0 # Mode has 0 growth rate at ky==0 phi2_vec = gsl_vector("phi2_by_mode_over_time", {:kx=>kxv, :ky=>kyv}) (@growth_rate_at_ky_at_kx[kyv][kxv] = 0.0; next) if phi2_vec.min <= 0.0 @growth_rate_at_ky_at_kx[kyv][kxv] = calculate_growth_rate(phi2_vec) (eputs "\n\n----------\nIn #@run_name:\n\nphi2_by_#{kxy}_over_time is all NaN; unable to calculate growth rates\n----------\n\n"; @growth_rate_at_ky_at_kx[kyv][kxv] = -1; next) if @growth_rate_at_ky_at_kx[kyv][kxv] == "NaN" end write_results end trap(0){} end @growth_rates = @growth_rate_at_ky @max_growth_rate = @growth_rates.values.max @fastest_growing_mode = @growth_rates.key(@max_growth_rate) @freq_of_max_growth_rate = @real_frequencies[@fastest_growing_mode] ep @max_growth_rate, @growth_rates @decaying = (@max_growth_rate < 0) if @max_growth_rate @ky = @aky if @aky if @grid_option == "single" # ep @aky, @growth_rates @gamma_r = @growth_rates[@aky.to_f] @gamma_i = @real_frequencies[@aky.to_f] end # ep @gamma_r # eputs @growth_rates; gets end end |
#calculate_results ⇒ Object
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 |
# File 'lib/gs2crmod/gs2.rb', line 258 def calculate_results return if ENV['CODE_RUNNER_NO_ANALYSIS'] =~ /true/ eputs "Analysing run" if @nonlinear_mode == "off" calculate_growth_rates_and_frequencies calculate_transient_amplifications elsif @nonlinear_mode == "on" calculate_saturation_time_index calculate_time_averaged_fluxes begin calculate_spectral_checks calculate_vspace_checks rescue end end @growth_rates ||={} @real_frequencies ||={} end |
#calculate_saturation_time_index(show_graph = false) ⇒ Object Also known as: csti
I.e. the time at which the primary modes are saturated and the fluxes settle around a long term average.
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 |
# File 'lib/gs2crmod/calculations.rb', line 118 def calculate_saturation_time_index(show_graph = false) eprint "Checking for saturation..." #hflux = gsl_vector('hflux_tot_over_time', {}) hflux = gsl_vector('phi2tot_over_time', {}) #eputs 'got hflux' #ep 'hflux', hflux #Check if it's decayed to 0 if hflux[-1] < 1.0e-10 for i in 1..hflux.size # raise "negative heat flux: #{hflux[-i]} " if hflux[-i] < 0 (break) unless hflux[- i] < 1.0e-10 end if i > hflux.size * 1.0/10.0 #i.e if was 0 for more than a tenth of the time @saturated = true @saturation_time_index = hflux.size - i + 1 eputs "saturation time = #{list(:t)[@saturation_time_index]}" GraphKit.quick_create([gsl_vector('t',{}), hflux]).gnuplot(log_axis: 'y') if show_graph return end end # Get initial estimate for saturation time for i in 0...hflux.size rem = hflux.subvector(i, hflux.size - i) break if (hflux[i] - rem.mean).abs < rem.sd / 2.0 break if i > 3.0/4.0*hflux.size end @saturation_time_index = [i + 1, hflux.size - 2].min # fit = GSL::Fit::linear(GSL::Vector.indgen(rem.size), rem) # # slope, covar11 = fit[1], fit[4] # range = [slope + Math.sqrt(covar11), slope - Math.sqrt(covar11)] # # unless range.min < 0 and range.max > 0 # eputs "Warning: This run (#{id}) has probably not reached a saturated state: the estimated slope of the heat flux is in this range: #{range.inspect}" # @saturated = false # end # # ep fit # eputs "Saturation time estimate', @saturation_time_index = i + 1 # t_vec[@saturation_time_index - 1] max_t_index = list(:t).keys.max max_t = list(:t).values.max min_t = list(:t).values.min #hflux = gsl_vector('hflux_tot_over_time', {:t_index_window => [@saturation_time_index, max_t_index]}) hflux = gsl_vector('phi2tot_over_time', {:t_index_window => [@saturation_time_index, max_t_index]}) t_vec = gsl_vector('t', {:t_index_window => [@saturation_time_index, max_t_index]}) # p t_vec[0] i = 0 t_arr = []; conf_arr = [] loop do eprint '.' # GraphKit.autocreate(x: {data: t_vec}, y: {data: hflux}).gnuplot lomb = GSL::SpectralAnalysis::Lomb.alloc(t_vec.subvector(i, t_vec.size - i), hflux.subvector(i, hflux.size - i)) fs, periodogram = lomb.calculate_periodogram(1.0, 4.0, [0]) #(1.0) #0.1 * hflux.size / ( hflux.size - i)) # lomb.graphkit.gnuplot # eputs 'Confidence that lowest frequency is not noise is: ' # pnoise is the probability of the strength of the lowest frequency signal in the heat flux given a hypothesis of gaussian noise. If it is high there is a low likelihood that there is a signal at the lowest frequency: ie. within that window the heat flux has reached a stationary state pnoise = lomb.pnull(periodogram[0]) t_arr.push t_vec[i]; conf_arr.push pnoise (@saturated = true; break) if pnoise > 0.9 step = (hflux.size / 25.0).to_i step = 1 if step==0 i += step #(@saturated = false; i ; break) if (i >= t_vec.size or t_vec[i] > (max_t - min_t) * 2.0 / 3.0 + min_t ) (@saturated = false; break) if (i >= t_vec.size or t_vec[i] > (max_t - min_t) * 2.0 / 3.0 + min_t ) @saturation_time_index += step # ep '---i,t,size',i, t_vec[i], t_vec.size end (kit = GraphKit.autocreate({x: {data: t_vec}, y: {data: hflux / hflux.max}}, {x: {data: t_arr}, y: {data: conf_arr}}); kit.data[1].with = 'lp'; kit.gnuplot) if show_graph #(log_axis: 'y') # puts if @saturated # p i eputs "saturation time = #{list(:t)[@saturation_time_index]}" else eputs "run not saturated" end return exit # Get regularly spaced t vector # # t_delta_vec = GSL::Vector.alloc(t_vec.size - 1) # t_delta_vec.size.times.each{|i| t_delta_vec[i] = t_vec[i+1] - t_vec[i]} # # ep t_delta_vec.max, t_delta_vec.min # # even_t = GSL::Vector.linspace(t_vec.min, t_vec.max, ((t_vec.max - t_vec.min) / t_delta_vec.max).round ) # # # even_t = [] # # tm = t = t_vec[t_delta_vec.max_index] # # # loop do # # even_t.push t # # # # ep even_t.size, t_vec.size # # min_delt = t_delta_vec.min # p even_t.any?{|el| bool = (not t_vec.any?{|ele| (ele - el).abs < 1.0e-1 * min_delt}); ep el if bool; bool} # # ep t_vec.dup.delete_if{|el| not (el - 71.3).abs < 0.5} # # exit return # Calculate a series of time averaged segments pieces = hflux.pieces(20) # split into 20 pieces avgs = GSL::Vector.alloc(pieces.map{|vec| vec.sum/vec.size}) # Calculate their variance mean = (avgs.sum/avgs.size) sig = Math.sqrt((avgs.square - mean**2).sum/avgs.size) # Discount any at the start which are more than one standard deviation away from the average - they are from the linear growth phase t_index = 1 kept_avgs = avgs.dup for i in 0...pieces.size if (avgs[i] - mean).abs > sig kept_avgs.delete_at(i) t_index += pieces[i].size else break end end eputs "Warning: probably not saturated" if [kept_avgs, kept_avgs.reverse].include? kept_avgs.sort ep kept_avgs @saturation_time_index = t_index # p t_index, list(:t)[t_index] end |
#calculate_spectral_checks ⇒ Object Also known as: csc
753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 |
# File 'lib/gs2crmod/calculations.rb', line 753 def calculate_spectral_checks ky_spec = gsl_vector('spectrum_over_ky') kx_spec = gsl_vector('spectrum_over_kx') kpar_spec = gsl_vector('spectrum_over_kpar', ky_index: ky_spec.max_index + 1, kx_index: 1) @spectrum_check = [] [kx_spec, ky_spec, kpar_spec].each do |spec| begin ends_max = [spec[0], spec[-1]].max + (10.0**(-9)) p ends_max p spec.max check = (Math.log(spec.max/ends_max)/Math.log(10)).round rescue check= -10 end @spectrum_check.push check end end |
#calculate_time_averaged_fluxes ⇒ Object Also known as: ctaf
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 |
# File 'lib/gs2crmod/calculations.rb', line 19 def calculate_time_averaged_fluxes eputs 'Calculating time averaged fluxes' calculate_saturation_time_index unless @saturation_time_index return unless FileTest.exist?(netcdf_filename) @hflux_tot_stav = saturated_time_average('hflux_tot_over_time', {}) @hflux_tot_stav_error = saturated_time_average_error('hflux_tot_over_time', {}) @hflux_tot_stav_std_dev = saturated_time_average_std_dev('hflux_tot_over_time', {}) @phi2_tot_stav = saturated_time_average('phi2tot_over_time', {}) #@par_mom_flux_stav = saturated_time_average('par_mom_flux_over_time', {}) rescue nil #@perp_mom_flux_stav = saturated_time_average('perp_mom_flux_over_time', {}) rescue nil @es_mom_flux_stav = {} @es_heat_flux_stav = {} @es_mom_flux_stav_error = {} @es_heat_flux_stav_error = {} @nspec.times do |i| species_index = i + 1 @es_mom_flux_stav[species_index] = saturated_time_average('es_mom_flux_over_time', {species_index: species_index}) @es_heat_flux_stav[species_index] = saturated_time_average('es_heat_flux_over_time', {species_index: species_index}) @es_mom_flux_stav_error[species_index] = saturated_time_average_error('es_mom_flux_over_time', {species_index: species_index}) @es_heat_flux_stav_error[species_index] = saturated_time_average_error('es_heat_flux_over_time', {species_index: species_index}) end # ep @es_mom_flux_stav, @es_heat_flux_stav end |
#calculate_transient_amplification(vector, options = {}) ⇒ Object
669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 |
# File 'lib/gs2crmod/calculations.rb', line 669 def calculate_transient_amplification(vector, ={}) turning_points = {} old = vector[0] i = 0 #for i in i...vector.size #new = vector[i] #if new > old #turning_points[:first_min] = i-1 #ep "First turning point[#{i}]\n" #break #end #old = new #end #for i in i...vector.size #new = vector[i] #if new < old #turning_points[:first_max] = i-1 #ep "Second turning point[#{i}]\n" #break #end #end #unless turning_points[:first_max] # and turning_points[:first_min] #return NaN #end ##t = gsl_vector('t') ##for j in 0...vector.size ##break if t[j] > 0.2 ##end #ep "vector[0..5]: #{vector.subvector(0,5)}\n" #return Math.sqrt(vector[turning_points[:first_max]]/@phiinit) #return vector.max/@phiinit vector[0] = 0 # This ensures vector.max does not return 1st point for no transient growth return vector.max/vector[1] end |
#calculate_transient_amplifications ⇒ Object Also known as: cta
417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 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 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 |
# File 'lib/gs2crmod/calculations.rb', line 417 def calculate_transient_amplifications return if @grid_option == "single" and @aky == 0.0 # no meaningful results Dir.chdir(@directory) do # With zero magnetic shear, calculate amplifications for both kx and ky if @shat and @shat.abs < 1.0e-5 and @nx > 1 to_calc = [:kx, :ky] @transient_amplification_at_kx ||= FloatHash.new else to_calc = [:ky] end @transient_amplification_at_ky ||= FloatHash.new eputs to_calc.each do |kxy| transient_amplifications = send(:transient_amplification_at_ + kxy) list(kxy).values.sort.each do |value| #p transient_amplifications.keys, value, transient_amplifications[value.to_f-0.0], #transient_amplifications.class, transient_amplifications.keys.include?(value); exit next if transient_amplifications.keys.include? value Terminal.erewind(1) #ep transient_amplifications.keys eputs sprintf("Calculating transient amplification for #{kxy} = % 1.5e#{Terminal::CLEAR_LINE}", value) # Mode has 0 growth rate at ky==0 (transient_amplifications[value] = 0.0; next) if value == 0.0 and kxy == :ky phi2_vec = gsl_vector("phi2_by_#{kxy}_over_time", {kxy=>value}) #(transient_amplifications[value] = 0.0; next) if phi2_vec.min <= 0.0 transient_amplifications[value] = calculate_transient_amplification(phi2_vec) (eputs "\n\n----------\nIn #@run_name:\n\nphi2_by_#{kxy}_over_time is all NaN; unable to calculate growth rate\n----------\n\n"; transient_amplifications[value] = -1; next) if transient_amplifications[value].to_s == "NaN" end end write_results # ep "transient_amplification_at_ky", @transient_amplification_at_ky if ENV['GS2_CALCULATE_ALL'] trap(0){eputs "Calculation of spectrum did not complete: run 'cgrf' (i.e. calculate_transient_amplifications_and_frequencies) for this run. E.g. from the command line \n $ coderunner rc 'cgrf' -j #{@id}"; exit} @transient_amplification_at_ky_at_kx ||= FloatHash.new list(:ky).values.sort.each do |kyv| @transient_amplification_at_ky_at_kx[kyv] ||= FloatHash.new #p @transient_amplification_at_ky_at_kx[kyv] list(:kx).values.sort.each do |kxv| next if @transient_amplification_at_ky_at_kx[kyv].keys.include? kxv Terminal.erewind(1) eputs sprintf("Calculating growth rate for kx = % 1.5e and ky = % 1.5e#{Terminal::CLEAR_LINE}", kxv, kyv) (@transient_amplification_at_ky_at_kx[kyv][kxv] = 0.0; next) if kyv == 0.0 # Mode has 0 growth rate at ky==0 phi2_vec = gsl_vector("phi2_by_mode_over_time", {:kx=>kxv, :ky=>kyv}) #(@transient_amplification_at_ky_at_kx[kyv][kxv] = 0.0; next) if phi2_vec.min <= 0.0 @transient_amplification_at_ky_at_kx[kyv][kxv] = calculate_transient_amplification(phi2_vec) (eputs "\n\n----------\nIn #@run_name:\n\nphi2_by_#{kxy}_over_time is all NaN; unable to calculate growth rates\n----------\n\n"; @transient_amplification_at_ky_at_kx[kyv][kxv] = -1; next) if @transient_amplification_at_ky_at_kx[kyv][kxv].to_s == "NaN" end write_results end trap(0){} end @transient_amplifications = @transient_amplification_at_ky @max_transient_amplification = @transient_amplifications.values.max @most_amplified_mode = @transient_amplifications.key(@max_transient_amplification) #@freq_of_max_transient_amplification = @real_frequencies[@fastest_growing_mode] #ep @max_transient_amplification, @transient_amplifications #@decaying = (@max_transient_amplification < 0) if @max_transient_amplification @ky = @aky if @aky #if @grid_option == "single" ## ep @aky, @transient_amplifications #@gamma_r = @transient_amplifications[@aky.to_f] #@gamma_i = @real_frequencies[@aky.to_f] #end # ep @gamma_r # eputs @transient_amplifications; gets end end |
#calculate_transient_es_heat_flux_amplifications ⇒ Object Also known as: ctehfa
499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 |
# File 'lib/gs2crmod/calculations.rb', line 499 def calculate_transient_es_heat_flux_amplifications return if @grid_option == "single" and @aky == 0.0 # no meaningful results @transient_es_heat_flux_amplification_at_species_at_kx = [] @transient_es_heat_flux_amplification_at_species_at_ky = [] @transient_es_heat_flux_amplification_at_species_at_ky_at_kx = [] for species_index in 1..nspec Dir.chdir(@directory) do # With zero magnetic shear, calculate amplifications for both kx and ky if @shat and @shat.abs < 1.0e-5 and @nx > 1 and !@ikx_init and false to_calc = [:kx, :ky] @transient_es_heat_flux_amplification_at_species_at_kx[species_index-1] ||= FloatHash.new else to_calc = [:ky] end @transient_es_heat_flux_amplification_at_species_at_ky[species_index-1] ||= FloatHash.new eputs to_calc.each do |kxy| transient_es_heat_flux_amplifications = send(:transient_es_heat_flux_amplification_at_species_at_ + kxy)[species_index-1] list(kxy).values.sort.each do |value| #p transient_es_heat_flux_amplifications.keys, value, transient_es_heat_flux_amplifications[value.to_f-0.0], #transient_es_heat_flux_amplifications.class, transient_es_heat_flux_amplifications.keys.include?(value); exit next if transient_es_heat_flux_amplifications.keys.include? value Terminal.erewind(1) #ep transient_es_heat_flux_amplifications.keys eputs sprintf("Calculating transient amplification for #{kxy} = % 1.5e#{Terminal::CLEAR_LINE}", value) # Mode has 0 growth rate at ky==0 (transient_es_heat_flux_amplifications[value] = 0.0; next) if value == 0.0 and kxy == :ky phi2_vec = gsl_vector("es_heat_by_#{kxy}_over_time", {kxy=>value, species_index: species_index}) #(transient_es_heat_flux_amplifications[value] = 0.0; next) if phi2_vec.min <= 0.0 transient_es_heat_flux_amplifications[value] = calculate_transient_amplification(phi2_vec) (eputs "\n\n----------\nIn #@run_name:\n\nphi2_by_#{kxy}_over_time is all NaN; unable to calculate growth rate\n----------\n\n"; transient_es_heat_flux_amplifications[value] = -1; next) if transient_es_heat_flux_amplifications[value].to_s == "NaN" end end write_results # ep "transient_es_heat_flux_amplification_at_species_at_ky", @transient_es_heat_flux_amplification_at_species_at_ky if ENV['GS2_CALCULATE_ALL'] trap(0){eputs "Calculation of spectrum did not complete: run 'ctehfa' (i.e. calculate_transient_es_heat_flux_amplifications) for this run. E.g. from the command line \n $ coderunner rc 'ctehfa' -j #{@id}"; exit} @transient_es_heat_flux_amplification_at_species_at_ky_at_kx[species_index-1] ||= FloatHash.new list(:ky).values.sort.each do |kyv| @transient_es_heat_flux_amplification_at_species_at_ky_at_kx[species_index-1][kyv] ||= FloatHash.new #p @transient_es_heat_flux_amplification_at_species_at_ky_at_kx[kyv] list(:kx).values.sort.each do |kxv| next if @transient_es_heat_flux_amplification_at_species_at_ky_at_kx[species_index-1][kyv].keys.include? kxv Terminal.erewind(1) eputs sprintf("Calculating growth rate for kx = % 1.5e and ky = % 1.5e#{Terminal::CLEAR_LINE}", kxv, kyv) (@transient_es_heat_flux_amplification_at_species_at_ky_at_kx[species_index-1][kyv][kxv] = 0.0; next) if kyv == 0.0 # Mode has 0 growth rate at ky==0 phi2_vec = gsl_vector("phi2_by_mode_over_time", {:kx=>kxv, :ky=>kyv}) #(@transient_es_heat_flux_amplification_at_species_at_ky_at_kx[kyv][kxv] = 0.0; next) if phi2_vec.min <= 0.0 @transient_es_heat_flux_amplification_at_species_at_ky_at_kx[species_index-1][kyv][kxv] = calculate_transient_es_heat_flux_amplification(phi2_vec) (eputs "\n\n----------\nIn #@run_name:\n\nphi2_by_#{kxy}_over_time is all NaN; unable to calculate growth rates\n----------\n\n"; @transient_es_heat_flux_amplification_at_species_at_ky_at_kx[species_index-1][kyv][kxv] = -1; next) if @transient_es_heat_flux_amplification_at_species_at_ky_at_kx[species_index-1][kyv][kxv].to_s == "NaN" end write_results end trap(0){} end #@max_transient_es_heat_flux_amplification = @transient_es_heat_flux_amplifications.values.max #@most_amplified_mode = @transient_es_heat_flux_amplifications.key(@max_transient_es_heat_flux_amplification) #@ky = @aky if @aky end end # for species_index in 1..nspec end |
#calculate_vspace_checks ⇒ Object Also known as: cvc
772 773 774 775 776 777 |
# File 'lib/gs2crmod/calculations.rb', line 772 def calculate_vspace_checks @vspace_check = ['lpc_pitch_angle', 'vres_pitch_angle', 'lpc_energy', 'vres_energy'].map do |name| saturated_time_average(name, {}) end end |
#check_converged ⇒ Object
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 |
# File 'lib/gs2crmod/check_convergence.rb', line 5 def check_converged raise CRFatal.new("It is strongly recommended that you do not use the use_large_cache option (-U) while checking convergence. Doing so will lead to unpredictable results.") if @runner.use_large_cache Dir.chdir(@directory) do logf(:check_converged) return if @checked_converged and not @runner.recalc_all log('@runner.class:', @runner.class) unless @runner.current_request == :check_converged @runner.requests.push :check_converged log 'check_converged requested recall' logi '@runner.requests', @runner.requests logi('@runner.object_id', @runner.object_id) return end return unless @status == :Complete eputs @run_name eputs @checked_converged = true log("finding similar resolutions") @runner.generate_combined_ids(:real) case @grid_option when "box" @similar_resolutions = @runner.similar_runs([:nx, :ny, :ntheta, :negrid, :naky, :ngauss, :nperiod, :delt, :jtwist], self) when "single" @similar_resolutions = @runner.similar_runs([:ntheta, :negrid, :naky, :ngauss, :nperiod], self) else raise CRFatal.new("Unknown grid option - can't get similar runs") end logi(@similar_resolutions) unless @similar_resolutions[1] eputs @run_name @converged = Feedback.get_boolean("This is is the biggest job with these params. Do you want to say it is converged?") return end @similar_resolutions.sort! do |id1, id2| run1 = @runner.run_list[id1] run2 = @runner.run_list[id2] if @grid_option == "box" and @nonlinear_mode == "off" (run1.jtwist*run1.nx*run1.negrid*run1.ngauss*run1.ntheta*run1.delt <=> run2.jtwist*run2.nx*run2.negrid*run2.ngauss*run2.ntheta*run2.delt) elsif @grid_option == "single" and @nonlinear_mode == "off" log("using nperiod: #{run1.nperiod}; #{run2.nperiod}") run1.negrid*run1.ngauss*run1.ntheta*run1.nperiod <=> run2.negrid*run2.ngauss*run2.ntheta*run2.nperiod elsif @naky run1.nx*run1.negrid*run1.ngauss*run1.ntheta*run1.naky <=> run2.nx*run2.negrid*run2.ngauss*run2.ntheta*run2.naky else run1.nx*run1.negrid*run1.ngauss*run1.ntheta*run1.ny <=> run2.nx*run2.negrid*run2.ngauss*run2.ntheta*run2.ny end end # eputs @similar_resolutions log("finding my place") my_place = @similar_resolutions.index(@id); # eputs my_place; gets if my_place > 0 last_job = @runner.run_list[@similar_resolutions[my_place - 1]] unless last_job.status == :Complete @checked_converged = false return end else @converged = false return end log("Checking overall convergence") #graph = graphkit('phi2tot_vs_time_all_kys') + #last_job.graphkit('phi2tot_vs_time_all_kys') #graph.gnuplot eputs "\n \n Warning: there are no bigger jobs" unless @similar_resolutions[my_place + 1] #@converged = Feedback.get_boolean("Is the plot converged?") #graph.close #(@checked_converged = true; return) unless @converged log("Checking convergence by ky") orn, last_job.runner = last_job.runner, nil log('last_job', last_job.pretty_inspect) last_job.runner = orn # last_job.get_ky_graphs; last_job.get_eigenfunctions # logi(last_job.ky_graphs) catch(:quit_converge_check) do = {} list(:ky).each do |index, ky| [:ky] = ky next if index == 1 and @grid_option == "box" graph = (graphkit('phi2_by_ky_vs_time', )+last_job.graphkit('phi2_by_ky_vs_time', )) graph.gnuplot answer = Feedback.get_choice("Is the graph converged?", ["Yes", "No", "The whole run is converged, stop pestering me!"]) graph.close case answer when /No/ @converged = false throw(:quit_converge_check) when /stop/ @converged = true throw(:quit_converge_check) when /Yes/ @converged = true end cgraph = lgraph = 'efnnormmag' graph = (graphkit('efnnormmag', )+last_job.graphkit('efnnormmag', )) # graph.gnuplot loop do graph.gnuplot answer = Feedback.get_choice('Is the graph converged?', ['Yes', 'No', 'The whole run is converged, stop pestering me!', 'Show me the magnitude of the eigenfunctions', 'Show me the real part of the eigenfunctions again', 'Normalise the eigenfunctions', 'Denormalise the eigenfunctions', 'Reverse the axis of the current run', 'Flip the current run', 'Toggle xrange']) graph.close case answer when /^Yes$/ @converged = true break when /^No$/ @converged = false throw(:quit_converge_check) when /stop/ @converged = true throw(:quit_converge_check) when /magnitude/ log 'checking convergence using magnitude' lgraph += 'mag'; cgraph += 'mag' when /Normalise/ log 'normalising' lgraph += 'norm'; cgraph += 'norm' when /Denormalise/ log 'denormalising' lgraph.gsub!(/norm/, ''); cgraph.gsub!(/norm/, '') when /real/ lgraph.gsub!(/mag/, ''); cgraph.gsub!(/mag/, '') when /Reverse/ cgraph = cgraph =~ /rev/ ? cgraph.sub!(/rev/, '') : cgraph + 'rev' # graph = (@eigenfunctions[ky]+last_job.eigenfunctions[ky]) when /Flip/ cgraph = cgraph =~ /flip/ ? cgraph.sub!(/flip/, '') : cgraph + 'flip' # graph = (@eigenfunctions[ky]+last_job.eigenfunctions[ky]) when /xrange/ if [:range] [:range] = nil else [:range] = 0 end else raise CRFatal.new("couldn't match choice #{answer}") end graph = graphkit(cgraph, ) + last_job.graphkit(lgraph, ) log graph.title end end end @checked_converged =true if last_job.checked_converged last_job.ky_graphs = nil last_job.eigenfunctions = nil # last_job.t_list = nil # last_job.kx_list = nil end # finish_processing end ep self end |
#check_parameters ⇒ Object
Eventually, this will be a full port of the ingen tool in the GS2 folder. At the moment it runs a limited set of tests for common errors in the input parameters (including type checking).
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 |
# File 'lib/gs2crmod/ingen.rb', line 97 def check_parameters # Sections # Namelist Tests # Grids # Parallelisation # Initialisation # Diagnostics # Misc # Namelist Tests rcp.namelists.each do |namelist, hash| next if hash[:should_include].kind_of? String and not eval(hash[:should_include]) if en = hash[:enumerator] #ep 'en', en, namelist next unless send(en[:name]) send(en[:name]).times do |i| run_namelist_tests(namelist, hash, i+1) end else run_namelist_tests(namelist, hash) end end ############### # Grid Errors # ############### # naky warning("Setting naky when non-linear mode is on is not recommended.") if @naky and @nonlinear_mode == "on" warning("You have set both ny and naky; naky will override ny.") if @ny and @naky error("abs(shat) should not be less that 1.0e-6") if @shat and @shat.abs < 1.0e-6 and not agk? error("abs(s_hat_input) should not be less that 1.0e-6") if @s_hat_input and @s_hat_input.abs < 1.0e-6 and not agk? # delt error("Please specify delt") unless @delt error("delt <= 0") if @delt <= 0.0 warning("Nonlinear run with delt_minimum unspecified.") if @nonlinear_mode=="on" and not @delt_minimum error("delt (#@delt) < delt_minimum") if @delt and @delt_minimum and @delt < @delt_minimum # negrid warning('negrid < 8 is not a good idea!') if @negrid and @negrid < 8 # nakx warning("You have set both nx and ntheta0; ntheta0 will override nx.") if @nx and @ntheta0 warning("Do you have a reason for setting equal_arc = true (default)? If not set false.") if @equilibrium_option=="eik" and (!@equal_arc or @equal_arc.fortran_true?) warning("Recommend nperiod > 1 for linear runs.") if @nonlinear_mode == "off" and (!@nperiod or @nperiod == 1) warning("Recommend nperiod = 1 for nonlinear runs.") if @nonlinear_mode == "on" and (@nperiod > 1) warning("Consider using field_option = local and associated optimizations.") if @field_option and @field_option == "implicit" ################################# # Parallelisation/Layout Errors # ################################# # Best linear run layout is lexys warning("The best layout for linear runs is usually lexys.") if @nonlinear_mode=="off" and not @layout=="lexys" # Best nonlinear run layout is xyles warning("The best layout for nonlinear runs is usually xyles.") if @nonlinear_mode=="on" and not @layout=="xyles" # Check whether we are parallelising over x warning("Parallelising over x: suggest total number of processors should be: #{max_nprocs_no_x}") if actual_number_of_processors > max_nprocs_no_x and not @grid_option == "single" ######################### # Initialisation Errors # ######################### # Check if restart folder exists if @restart_file and @restart_file =~ /^(?<folder>[^\/]+)\// folder = $~[:folder] warning("Folder #{folder}, specified in restart_file, not present. NetCDF save may fail") unless FileTest.exist?(folder) end error("Setting @restart_file as an empty string will result in hidden restart files.") if @restart_file == "" error("ginit_option is 'many' but is_a_restart is false") if @ginit_option == "many" and not @is_a_restart error("chop_side should not be used (remove test if default changes from T to F)") if !@chop_side or @chop_side.fortran_true? ##################### # Diagnostic errors # ##################### #Check whether useful diagnostics have been omitted. not_set = [:write_verr, :save_for_restart, :write_nl_flux, :write_final_fields, :write_final_moments].find_all do |diagnostic| not (send(diagnostic) and send(diagnostic).fortran_true?) end if not_set.size > 0 str = not_set.inject("") do |s, diagnostic| s + "\n\t#{diagnostic} --- " + rcp.namelists[diagnostics_namelist][:variables][diagnostic][:description] rescue s end warning("The following useful diagnostics were not set:" + str) if str.length > 0 end warning("You are running in nonlinear mode but have not switched the nonlinear flux diagnostic.") if not (@write_nl_flux and @write_nl_flux.fortran_true?) and @nonlinear_mode == "on" #{ #write_verr: "Velocity space diagnostics will not be output for this run" #}.each do |var, warn| #warning(v"#{var} not set or .false. --- " + warn) unless send(var) and send(var).fortran_true? #end error("Please specify nwrite") unless @nwrite error("Please specify nstep") unless @nstep warning("You will write out diagnostics less than 50 times") if @nstep/@nwrite < 50 ######################## # Miscellaneous errors # ######################## error("The run name for this run is too long. Please move some of the variable settings to the local defaults file.") if @relative_directory.size + @run_name.size > MAX_NAME_SIZE warning("You are submitting a nonlinear run with no dissipation.") if @nonlinear_mode == "on" and @hyper_option=="none" and @collision_model=="none" warning("You have no spacial implicitness: (bakdif) for one of your species. Be prepared for numerical instabilities!") if (1..@nspec).to_a.find{|i| bd = send("bakdif_#{i}") and bd == 0} warning("The system will abort with rapid timestep changes...") if !@abort_rapid_time_step_change or @abort_rapid_time_step_change.fortran_true? warning("local_field_solve is an old variable that should not really be used.") if @local_field_solve and @local_field_solve.fortran_true? ############################# # Boundary Condition Errors # ############################# error("Boundary options should be linked with finite magnetic shear when running nonlinear_runs.") if (!@boundary_option or @boundary_option != "linked") and ((@s_hat_input and @s_hat_input.abs > 1.0e-6) or (@shat and @shat.abs > 1.0e-6)) and not @nonlinear_mode == "off" warning("Boundary options should be probably be linked with finite magnetic shear") if (!@boundary_option or @boundary_option != "linked") and ((@s_hat_input and @s_hat_input.abs > 1.0e-6) or (@shat and @shat.abs > 1.0e-6)) warning("You are using boundary_option = linked with grid_option = range: this is probably an error.") if @boundary_option == "linked" and @grid_option == "range" error("Set nonad_zero = true.") if @nonad_zero and not @nonad_zero.fortran_true? ################### # Spectrogk tests # ################### # if spectrogk? if @force_5d and @force_5d.fortran_true? warning("Must specify interpolation method with phi_method.") if not (@phi_method) end end ################ # Damping Rate # ################ error("Linear runs with hyperviscosity are NOT recommended!") if @nonlinear_mode=="off" and (@hyper_option and @hyper_option=="visc_only") and (@d_hypervisc and @d_hypervisc!=0) warning("Amplitude dependent part of hyperviscosity being ignored since const_amp = true") if (@hyper_option and @hyper_option=="visc_only") and (@const_amp and @const_amp.fortran_true?) ################### # Geometry Errors # ################### error("You must set bishop = 4 for Miller(local) geometry. Remember also that s_hat_input will override shat") if (@bishop!=4 and (@local_eq and @local_eq.fortran_true?)) error("Shift should be > 0 for s-alpha equilibrium.") if @equilibrium_option=="s-alpha" and (@shift and @shift < 0) error("Shift should be < 0 for Miller equilibrium.") if @equilibrium_option=="eik" and @local_eq.fortran_true? and (@shift and @shift > 0) error("irho must be 2 for Miller equilibrium.") if @equilibrium_option=="eik" and @local_eq.fortran_true? and (@irho and @irho!=2) warning("Note that shat != s_hat_input") if @shat and @s_hat_input and @shat!=@s_hat_input ################## # Species Errors # ################## error("Must set z = -1 for electron species.") if (@type_2 and @z_2 and @type_2=='electron' and @z_2 != -1) ################# # Optimisations # ################# if CODE_OPTIONS[:gs2] and CODE_OPTIONS[:gs2][:show_opt] eputs("Optimisation Summary:") optimisation_flags.each do |flag| eputs("------------------------- #{flag}: #{send(flag)}\n* #{rcp.variables_with_help[flag].gsub(/\n/, "\n\t").sub(/\A([^.]*.).*\Z/m, '\1')}") end #not_set = [:operator, :save_for_restart, :write_nl_flux, :write_final_fields, :write_final_moments].find_all do |diagnostic| #not (send(diagnostic) and send(diagnostic).fortran_true?) #end #if not_set.size > 0 #str = not_set.inject("") do |s, diagnostic| #s + "\n\t#{diagnostic} --- " + rcp.namelists[diagnostics_namelist][:variables][diagnostic][:description] rescue s #end #warning("The following useful diagnostics were not set:" + str) if str.length > 0 #end end end |
#code_run_environment ⇒ Object
90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 |
# File 'lib/gs2crmod/gs2.rb', line 90 def code_run_environment case CodeRunner::SYS when /iridis/ <<EOF module load openmpi EOF when /helios/ <<EOF module purge module load intel module load bullxmpi module load netcdf_p module load hdf5_p module load fftw/3.3.3 module load bullxde papi module load scalasca EOF #when /archer/ #<<EOF #module swap PrgEnv-cray PrgEnv-intel #module load intel/14.0.0.080 #module load fftw #module load netcdf-hdf5parallel #module load cray-hdf5-parallel #EOF else @code_run_environment end end |
#corrected_mom_flux_stav ⇒ Object
Not needed for GS2 built after 16/06/2010
413 414 415 |
# File 'lib/gs2crmod/calculations.rb', line 413 def corrected_mom_flux_stav par_mom_flux_stav - perp_mom_flux_stav end |
#correlation_analysis(options = {}) ⇒ Object
This function will handle running the correlation analysis and writing the results to a NetCDF file. Cases need to be handled differently since perp, par and full are just subsets of the full correlation function but the time correlation calculation needs to deal with each radial location separately. Time correlation uses the zonal flows in the toroidal direction to calculate the correlation time.
This function takes in the same options as field_real_space_standard_representation, along with the following new options dealing with interpolation and binning:
correlation_type: determines which subset of correlation function should be calculated (perp/par/full/time) nbins_array: array giving number of bins to use in the binning procedure. Index order (x, y, z ,t) nt_reg: Most of the time you have many more time points than you need for spatial correlations. This sets
number of new interpolation points in time.
Using this function: Since this can only be single threaded, this can be a very expensive calculation when trying to do the full correlation function, so this is not recommended for highly resolved nonlinear runs. This is why the perp/par/full splitting is implemented, allowing one dimension to be taken out essentially.
1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 |
# File 'lib/gs2crmod/gs2.rb', line 1194 def correlation_analysis(={}) #Sanity checks: #Cannot only have one bin since require difference between bins for index calculation if [:nbins_array].include?1 raise('Cannot have only one bin in nbins_array. Minuimum is two.') end #Thetamin shouldn't be equal to thetamax to avoid possibili # case [:correlation_type] when 'perp', 'par', 'full' gsl_tensor = field_correlation_gsl_tensor() shape = gsl_tensor.shape #Set up dimensions file = NumRu::NetCDF.create(@run_name + "_correlation_analysis_#{[:correlation_type]}.nc") ydim = file.def_dim('x',shape[0]) xdim = file.def_dim('y',shape[1]) zdim = file.def_dim('z',shape[2]) tdim = file.def_dim('t',shape[3]) correlation_var = file.def_var("correlation", 'sfloat', [xdim, ydim, zdim, tdim]) file.enddef #Write out array correlation_var.put(NArray.to_na(gsl_tensor.to_a)) file.close when 'time' nakx_actual = NumRu::NetCDF.open(@run_name + ".out.nc").var('kx').get kx_len = nakx_actual.length if [:nakx] == nil radial_pts = kx_len elsif [:nakx] <= kx_len radial_pts = [:nakx] else raise('nakx exceeds the total number of kx\'s in simulation') end #Check whether t_index_window is specified, if not, set to entire t range if [:t_index_window] == nil [:t_index_window] = [1, -1] end #Now loop through the radial locations and calculate the correlation function in y and t. for x in 0...radial_pts [:xmin] = x [:xmax] = x gsl_tensor = field_correlation_gsl_tensor() shape = gsl_tensor.shape if x == 0 #Write dimensions to NetCDF file file = NumRu::NetCDF.create(@run_name + "_correlation_analysis_#{[:correlation_type]}.nc") ydim = file.def_dim('x',shape[0]) xdim = file.def_dim('y',shape[1]) zdim = file.def_dim('z',shape[2]) tdim = file.def_dim('t',shape[3]) end file.redef correlation_var = file.def_var("correlation_x_#{x}", 'sfloat', [xdim, ydim, zdim, tdim]) file.enddef #Write out array correlation_var.put(NArray.to_na(gsl_tensor.to_a)) end file.close #only close after loop over radial points else raise 'Please specify correlation_type as perp/par/time/full' end end |
#ctan ⇒ Object
706 707 708 709 710 711 712 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 741 |
# File 'lib/gs2crmod/calculations.rb', line 706 def ctan list(:ky).each do |(ky_index, ky)| eputs "ky: #{ky}" phi_vec = gsl_vector("phi2_by_ky_over_time", ky_index: ky_index) t_element = 0 old = phi_vec[0] loop do t_element+=1 #print t_element, ',', phi_vec.size new = phi_vec[t_element] break if new > old or t_element == phi_vec.size - 1 old = new end if t_element == phi_vec.size - 1 @transient_amplification_at_ky[ky] = -1 eputs "No Min" next end first_min = t_element eputs "ky: #{ky}, first_min: #{first_min}" loop do t_element+=1 #print t_element, ',', phi_vec.size new = phi_vec[t_element] break if new < old or t_element == phi_vec.size - 1 end if t_element == phi_vec.size - 1 @transient_amplification_at_ky[ky] = -1 next end @transient_amplification_at_ky[ky] = phi_vec.subvector(t_element, phi_vec.size - t_element).max end end |
#cumulative_gridpoints ⇒ Object
339 340 341 342 343 |
# File 'lib/gs2crmod/ingen.rb', line 339 def cumulative_gridpoints c = 1 error("Please specify layout") unless @layout @layout.split(//).reverse.inject({}){|hash, let| c*=gridpoints[let]; hash[let] = c; hash} end |
#data_string ⇒ Object
421 422 423 424 425 426 427 428 429 |
# File 'lib/gs2crmod/gs2.rb', line 421 def data_string logf(:data_string) return "" unless @converged unless @grid_option == 'single' logi(@ky, @growth_rates, @real_frequencies) # log(:@@readout_list, @@readout_list) return rcp.readout_list.inject(""){|str,(var,_)| str+"#{(send(var) || "0")}\t"} + "\n" # @ky ? (@@variables + @@results - ).inject(""){|str,(var,type_co)| str+"#{(send(var) || "0")}\t"} + sprintf("%e\t%e\t%e\n", @ky, @growth_rates[@ky], @real_frequencies[@ky]) : (@@variables + @@results).inject(""){|str,(var,type_co)| str+"#{(send(var) || "0")}\t"} + sprintf("%e\t%e\t%e\n", @fastest_growing_mode, @max_growth_rate, @freq_of_max_growth_rate) end |
#delete_restart_files(options = {}) ⇒ Object
Delete all the restart files (irreversible!)
614 615 616 617 618 619 |
# File 'lib/gs2crmod/gs2.rb', line 614 def delete_restart_files(={}) puts 'You are about to delete the restart files for:' puts @run_name return unless Feedback.get_boolean("This action cannot be reversed. Do you wish to continue?") unless [:no_confirm] list_of_restart_files.each{|file| FileUtils.rm file} end |
#diagnostics_namelist ⇒ Object
351 352 353 |
# File 'lib/gs2crmod/ingen.rb', line 351 def diagnostics_namelist :gs2_diagnostics_knobs end |
#error(message) ⇒ Object
14 15 16 |
# File 'lib/gs2crmod/ingen.rb', line 14 def error() raise InputFileError.new("Error: " + ) end |
#estimated_nodes ⇒ Object Also known as: estnod
Gives a guess as to the maximum number of nodes which can be can be utilized on the current system
896 897 898 |
# File 'lib/gs2crmod/gs2.rb', line 896 def estimated_nodes parallelizable_meshpoints / max_ppn end |
#eulerian_kx_index(options) ⇒ Object
This function is used in the presence of perpendicular flow shear. It returns the (Eulerian) GS2 kx_index as a function of the Lagrangian kx, which is the kx_index of the mode in a shearing coordinate system, I.e. if you give it an Lagrangian kx (which is the same as the Eulerian kx at t=0) it will tell you where it has now got to. It may have left the box, in which case this function will return an error.
A given Lagrangian kx moves through the GS2 box, and thus for such a kx the response matrix varies in time. This is done because the effect of flow shear can be reduced by a shearing coordinate transformation to become merely a time varying kx.
At each timestep, phi(ikx_indexed(it)) is set equal to phi(ikx_indexed(it - jump(iky)) kx_indexed is defined in the following way.
do it=itmin(1), ntheta0
ikx_indexed (it+1-itmin(1)) = it end do
do it=1,itmin(1)-1 ikx_indexed (ntheta0 - itmin(1) + 1 + it)= it end do
In other words, what this means is that akx(ikx_indexed(0)) is the minimum kx, and that akx(ikx_indexed(ntheta0)) gives the maximum kx, kx_indexed moves the kxs out of box order.
So. remembering that jump is negative, phi(kx) is set equal phi(kx - jump * dkx) so the Lagrangian mode has moved to a lower kx. So get the Eulerian index, one starts with the Lagrangian index, and adds jump (which is negative!). This, however, must be done with indexes that are in the physical (not box) order. So this function first moves the indexes out of box order, then adds jump, then moves them back into box order so that the index returned will give the correct kx from the GS2 array.
1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 |
# File 'lib/gs2crmod/gsl_data.rb', line 1440 def eulerian_kx_index() #eputs "Start eulerian_kx_index" lagrangian_kx_index = [:kx_index] phys = physical_kx_index(lagrangian_kx_index) #ep 'jump', jump(options) index = phys + jump() raise ArgumentError.new("Lagrangian kx out of range") if index <= 0 box= box_kx_index(index) #eputs "End eulerian_kx_index" return box end |
#generate_component_runs ⇒ Object
314 315 316 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 353 354 355 356 357 358 |
# File 'lib/gs2crmod/gs2.rb', line 314 def generate_component_runs @component_runs = [] logf(:generate_component_runs) return if @grid_option == "single" and @scan_type == "none" begin list(:ky) # This will fail unless the run has output the netcdf file rescue return end return unless @status == :Complete #and @converged log(@run_name) if @grid_option == "box" and @nonlinear_mode == "off" @ky = nil # raise CRFatal.new("no @ky_list") unless @ky_list # log list(:ky) list(:ky).each do |id, ky| component_run = create_component #self.dup component_run.ky = ky component_run.gamma_r = @growth_rates[ky] component_run.gamma_i = @real_frequencies[ky] log @runner.component_ids # log('@runner.class', @runner.class) # @runner.add_component_run(component_run) end elsif (not gryfx?) and @scan_type and @scan_type != "none" t = gsl_vector('t') scan_vals = gsl_vector('scan_parameter_value') current = scan_vals[0] start = 0 for i in 0...t.size if scan_vals[i] != current component = create_component component.scan_index_window = [start+1, i] #remember indexes are elements + 1 #ep 'scan_index_window', component.scan_index_window component.scan_parameter_value = current component.growth_rate_at_ky = nil component.growth_rate_at_kx = nil component.growth_rate_at_ky_at_kx = nil component.calculate_results current = scan_vals[i] start = i end end end end |
#generate_input_file(&block) ⇒ Object
805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 |
# File 'lib/gs2crmod/gs2.rb', line 805 def generate_input_file(&block) raise CRFatal("No Input Module File Given or Module Corrupted") unless methods.include? (:input_file_text) run_namelist_backwards_compatibility if @restart_id and (not @is_a_restart or @resubmit_id) # The second test checks that the restart function has not been called manually earlier (e.g. in Trinity), but we must check that it is not in fact a resubmitted run @runner.run_list[@restart_id].restart(self) elsif @save_for_restart and @save_for_restart.fortran_true? and (not @is_a_restart or @resubmit_id) @restart_dir = "nc" #if CODE_OPTIONS[:gs2] and CODE_OPTIONS[:gs2][:list] #FileUtils.makedirs "#{@runner.root_folder}/#@restart_dir" #else FileUtils.makedirs @restart_dir #end @restart_file = "#@run_name.nc" end # Let Gs2 know how much wall clock time is available. avail_cpu_time is a GS2 input parameter. @avail_cpu_time = @wall_mins * 60 if @wall_mins # Automatically set the number of nodes to be the maximum possible without parallelising over x, if the user has left the number of nodes unspecified. set_nprocs if block ##### Allow the user to define their own pre-flight checks and changes instance_eval(&block) else ######### Check for errors and inconsistencies check_parameters ######### end write_input_file ######### Generate a report using the ingen tool if possible ingen unless block ######## end |
#get_completed_timesteps ⇒ Object
381 382 383 384 385 386 387 388 389 |
# File 'lib/gs2crmod/gs2.rb', line 381 def get_completed_timesteps #raise CRFatal.new("Couldn't find outfile #{@run_name}.out") unless FileTest.exist?(@run_name + ".out") #p 'try to get completed_timesteps', Dir.pwd, 'nwrite', @nwrite, 'delt', @delt @completed_timesteps = (list(:t).size - 1) * (@nwrite || 1) #p 'tried to get completed_timesteps' #rescue #`grep time= #@run_name.out`.split.size # File.read("#@run_name.out").scan(/^\s+time\s*=\s+/).size * @nwrite end |
#get_list_of(*args) ⇒ Object Also known as: list
462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 |
# File 'lib/gs2crmod/gs2.rb', line 462 def get_list_of(*args) #args can be any list of e.g. :ky, :kx, :theta, :t ... logf(:get_list_of) refresh = args[-1] == true ? true : false args.pop if args[-1] == true logd Dir.chdir(@directory) do args.each do |var| # eputs "Loading #{var}" list_name = var + :_list log list_name # self.class.send(:attr_accessor, list_name) next if (cache[list_name] and [:Failed, :Complete].include? status and not refresh) cache[list_name] = {} if netcdf_file.var(var.to_s) netcdf_file.var(var.to_s).get.to_a.each_with_index do |value, element| # print '.' cache[list_name][element+1]=value end else netcdf_file.dim(var.to_s).length.times.each do |element| cache[list_name][element+1]='unknown' end end # eputs send(var+:_list) end end logfc :get_list_of return cache[args[0] + :_list] if args.size == 1 end |
#get_run_time ⇒ Object
Try to read the runtime in minutes from the GS2 standard out.
285 286 287 288 289 290 291 292 293 294 295 296 |
# File 'lib/gs2crmod/gs2.rb', line 285 def get_run_time logf(:get_run_time) output = @output_file || try_to_get_output_file return nil unless output begin Regexp.new("total from timer is:\\s*#{LongRegexen::NUMBER}", Regexp::IGNORECASE).match FileUtils.tail(output, 300) logi $~ @run_time = $~[:number].to_f rescue @run_time = nil end end |
#get_status ⇒ Object
667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 |
# File 'lib/gs2crmod/gs2.rb', line 667 def get_status # eputs 'Checking Status' logf(:get_status) Dir.chdir(@directory) do if @running if FileTest.exist?(@run_name + ".out") and FileUtils.tail(@run_name + ".out", 5).split(/\n/).size > 4 and FileUtils.tail(@run_name + ".out", 200) =~ /t\=/ @status = :Incomplete else @status = :NotStarted end else if FileTest.exist?(@run_name + ".out") and FileUtils.tail(@run_name + ".out", 5).split(/\n/).size > 4 #eputs "HERE", @scan_type if @nonlinear_mode == "off" and FileUtils.tail(@run_name + ".out",200) =~ /omega converged/ eputs 'Omega converged...' @status = :Complete elsif @scan_type and @scan_type != "none" and FileUtils.tail(@run_name + ".par_scan",200) =~ /scan\s+is\s+complete/i eputs 'Scan complete...' @status = :Complete elsif @nonlinear_mode == "on" or !@omegatol or @omegatol < 0.0 or (@exit_when_converged and @exit_when_converged.fortran_false?) eputs 'No omegatol' if FileTest.exist?(@run_name + ".out.nc") #p ['pwd', Dir.pwd, netcdf_file, netcdf_file.dim('t'), netcdf_file.dims] if netcdf_file.dim('t').length > 0 get_completed_timesteps else @status = :Failed return end else eputs "Warning: no netcdf file #@run_name.out.nc" @status = :Failed return end #ep "completed_timesteps", @completed_timesteps eputs "#{percent_complete}% of Timesteps Complete" if percent_complete >= 100.0 @status = :Complete elsif percent_complete > 5 and FileUtils.tail(output_file, 200) =~ /total from timer is/ @status = :Complete else @status = :Failed end else @status = :Failed end else @status=:Failed end end end end |
#get_time ⇒ Object
362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 |
# File 'lib/gs2crmod/gs2.rb', line 362 def get_time begin lt = list(:t) return lt.values.max if lt.size>0 rescue end time = nil # eputs File.readlines(@run_name +".out").slice(-4..-1).reverse.join( "\n"); gets raise CRFatal.new("Couldn't find outfile #{@run_name}.out") unless FileTest.exist?(@run_name + ".out") tail = FileUtils.tail("#@run_name.out", 4) #File.readlines(@run_name +".out").slice(-4..-1).reverse.join( "\n") tail.sub(LongRegexen::FLOAT) do # eputs $~.inspect time = $~[:float].to_f end #if FileTest.exist? (@run_name +".out") #raise CRFatal.new("couldn't get the time from #{tail}") unless time @time = time end |
#graphkit(name, options = {}) ⇒ Object
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 |
# File 'lib/gs2crmod/graphs.rb', line 188 def graphkit(name, ={}) logf :graphkit # If an array of t, kx or ky values is provided, plot one graph for each value and then sum the graphs together [:t, :kx, :ky, :X, :Y, :e, :l, :theta].each do |var| #ep 'index', var if [var].class == Symbol and [var] == :all [var] = list(var).values elsif [var+:_index].class == Symbol and [var+:_index] == :all #ep 'Symbol' [var+:_index] = list(var).keys end if [var].class == Array return [var].map{|value| graphkit(name, .dup.absorb({var => value}))}.sum elsif [var+:_index].class == Array #ep 'Array' return [var+:_index].map{|value| graphkit(name, .dup.absorb({var+:_index => value}))}.sum end if [var].class == Symbol and [var] == :max [var] = list(var).values.max elsif [var+:_index].class == Symbol and [var+:_index] == :max ep 'Symbol' [var+:_index] = list(var).keys.max end end [:t_index] ||= [:frame_index] if [:frame_index] # Smart graphkits are defined in the file read_netcdf if name =~ /^cdf_/ return smart_graphkit( + {graphkit_name: name}) elsif name =~ /^nc_/ return old_smart_graphkit( + {graphkit_name: name}) end # If a method from the new GraphKits module can generate this graphkit use it if method = self.class.instance_methods.find{|meth| (name + '_graphkit').to_sym == meth} [:graphkit_name] = name return send(method, ) end raise "Graph #{name} not found" end |
#gridpoints ⇒ Object
A hash which gives the actual numbers of gridpoints indexed by their corresponding letters in the layout string.
329 330 331 332 333 334 335 336 337 |
# File 'lib/gs2crmod/ingen.rb', line 329 def gridpoints gridpoints = {'l' => @ngauss, 'e' => @negrid, 's' => @nspec} if @grid_option == "single" gridpoints.absorb({'x'=>1, 'y'=>1}) else gridpoints.absorb({'x' => (@ntheta0 or (2.0 * (@nx - 1.0) / 3.0 + 1.0).floor), 'y' => (@naky or ((@ny - 1.0) / 3.0 + 1.0).floor)}) end return gridpoints end |
#gryfx? ⇒ Boolean
67 68 69 |
# File 'lib/gs2crmod/gs2.rb', line 67 def gryfx? false end |
#gsl_complex(name, options = {}) ⇒ Object
1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 |
# File 'lib/gs2crmod/gsl_data.rb', line 1488 def gsl_complex(name, ={}) = eval() if .class == String # p @directory Dir.chdir(@directory) do # eputs Dir.pwd case name when /correcting_phase/ # options.convert_to_index(self, :ky) # theta0 = (options[:theta0] or 0) # # p 'options[:ky_index]', options[:ky_index] # phase_array = NumRu::NetCDF.open("#@directory/#@run_name.out.nc").var('phase').get({"start" => [0, options[:ky_index] - 1, theta0], 'end' => [1, options[:ky_index] - 1, theta0] }).to_a.flatten # p 'phase_array', phase_array # thetaelement0 = (list(:theta).key(0.0) - 1).to_i # p 'list(:theta)[thetaelement0 + 1]', list(:theta)[thetaelement0 + 1] # p 'thetaelement0', thetaelement0 # p 'theta0 - jump(options)', theta0 - jump(options) % @jtwist # p 'list(:kx)[2] * (theta0 - jump(options)%@jtwist)', list(:kx)[2] * (theta0 - jump(options)%@jtwist) # kx_element = list(:kx).key(list(:kx)[2] * (theta0 - jump(options)%@jtwist)) - 1 # at_0 = NumRu::NetCDF.open("#@directory/#@run_name.out.nc").var('phi').get({"start" => [0, thetaelement0, kx_element, options[:ky_index] - 1], 'end' => [1, thetaelement0, kx_element, options[:ky_index] - 1] }).to_a.flatten # p 'at_0', at_0 # at_0 = GSL::Complex.alloc(at_0) # p 'at_0', at_0 # return (at_0 / at_0.mag).conj # # pp 'theta0', theta0 # # pp phase_array[5][theta0] # return GSL::Complex.alloc(phase_array) # # new_options = options.dup # new_options[:imrc] = :real # thetas = gsl_vector('theta_along_field_line', new_options) # at_0 = gsl_vector_complex('phi_along_field_line', new_options)[.to_a.index(0.0)] # p at_0 exit else raise CRError.new("Unknown gsl_complex requested: #{name}") end # eputs data; gets end end |
#gsl_matrix(name, options = {}) ⇒ Object
1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 |
# File 'lib/gs2crmod/gsl_data.rb', line 1155 def gsl_matrix(name, ={}) = eval() if .class == String if [:saturated_time_average] or [:sta] raise "Not Saturated" unless @saturation_time_index tmax = list(:t).keys.max return ((@saturation_time_index..tmax).to_a.map do |t_index| gsl_matrix(name, .dup.absorb({t_index: t_index, saturated_time_average: nil, sta: nil})) end).sum / (list(:t).values.max - list(:t)[@saturation_time_index]) end if method = self.class.instance_methods.find{|meth| (name + '_gsl_matrix').to_sym == meth} [:graphkit_name] = name return send(method, ) end end |
#gsl_tensor(name, options) ⇒ Object
132 133 134 |
# File 'lib/gs2crmod/gsl_data_3d.rb', line 132 def gsl_tensor(name, ) tensor = send((name.to_s+"_gsl_tensor").to_sym , ) end |
#gsl_vector(name, options = {}) ⇒ Object
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 |
# File 'lib/gs2crmod/gsl_data.rb', line 122 def gsl_vector(name, ={}) Dir.chdir(@directory) do [:t_index_window] ||= @scan_index_window .setup_time_window if [:ky, :kx].include? name.to_sym vec = fix_norm( GSL::Vector.alloc(netcdf_file.var(name.to_s).get.to_a.sort), -1, ) # ky, ky are normalised to 1 / rho_i if i = [:interpolate_ + name.to_s.sub(/k/, '').to_sym] if name.to_sym == :ky s = (vec.size - 1)*i + 1 #return vec.connect(GSL::Vector.alloc((vec.size-1)*(i-1)) * 0.0) return (0...s).map{|k| k.to_f * vec[1]}.to_gslv else size = vec.size #vec = vec.to_box_order raise "Hmmm, kx.size should be odd" unless size%2 == 1 s = (size-1)/2 * i return (-s..s).to_a.map{|i| i.to_f * vec.to_box_order[1]}.to_gslv #new_vec = GSL::Vector.alloc((s-1)*i + 1) #new_vec *= 0.0 #for j in 0...((s-1)/2+1) #new_vec[j] = vec[j] #end #for j in 0...((s-1)/2) #new_vec[-j-1] = vec[-j-1] #end #return new_vec.from_box_order end else return vec end elsif [:theta].include? name.to_sym #ep options; gets #vec = GSL::Vector.alloc(netcdf_file.var(name.to_s).get({'start' => [options[:thetamin]||0], 'end' => [options[:thetamax]||-1]}).to_a) vec = GSL::Vector.alloc(netcdf_file.var(name.to_s).get.to_a) if gryfx? and [:periodic] #vec = vec.connect([2.0*vec[-1] - vec[-2]].to_gslv) vec = vec.connect([-vec[0]].to_gslv) end if ith = [:interpolate_theta] osize = vec.size newsize = (osize-1)*ith+1 newvec = GSL::Vector.alloc(newsize) newvec[newsize-1] = vec[osize-1]# * ith.to_f for i in 0...(newsize-1) im = i%ith frac = im.to_f/ith.to_f #iold = (i-im)/(new_shape[-1]-1)*(shape[-1]-1) iold = (i-im)/ith newvec[i] = (vec[iold] * (1.0-frac) + vec[iold+1] * frac) end vec = newvec end start = [:thetamin]||0 endv = [:thetamax]||vec.size-1 #ep ['options', options, 'vec.size', vec.size] vec = vec.subvector(start, (endv-start+1)).dup return vec elsif name.to_sym == :t #options.setup_time_window t = GSL::Vector.alloc(netcdf_file.var(name.to_s).get('start' => [[:begin_element]], 'end' => [[:end_element]]).to_a) t = t - t[0] if [:sync_time] return fix_norm(t, -1, ) # t is normalised to a/v_thi end = eval() if .class == String if [:saturated_time_average] or [:sta] raise "Not Saturated" unless @saturation_time_index tmax = list(:t).keys.max return ((@saturation_time_index..tmax).to_a.map do |t_index| gsl_vector(name, .dup.absorb({t_index: t_index, saturated_time_average: nil, sta: nil})) end).sum / (list(:t).values.max - list(:t)[@saturation_time_index]) elsif [:time_average] or [:ta] tmax = list(:t).keys.max start_t = 2 return ((start_t..tmax).to_a.map do |t_index| gsl_vector(name, .dup.absorb({t_index: t_index, time_average: nil, ta: nil})) end).sum / (list(:t).values.max - list(:t)[start_t]) end if method = self.class.instance_methods.find{|meth| (name + '_gsl_vector').to_sym == meth} [:graphkit_name] = name return send(method, ) end end raise "GSL Vector #{name} not found" end |
#gsl_vector_complex(name, options = {}) ⇒ Object
1057 1058 1059 1060 1061 1062 1063 1064 |
# File 'lib/gs2crmod/gsl_data.rb', line 1057 def gsl_vector_complex(name, ={}) = eval() if .class == String if method = self.class.instance_methods.find{|meth| (name + '_gsl_vector_complex').to_sym == meth} [:graphkit_name] = name return send(method, ) end end |
#has_electrons? ⇒ Boolean
16 17 18 |
# File 'lib/gs2crmod/properties.rb', line 16 def has_electrons? return @nspec.times.inject(false){|bool, i| bool or send(:type_ + i.to_sym) =~ /electrons/i} end |
#hypercoll_graphkit(options) ⇒ Object
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 |
# File 'lib/gs2crmod/read_netcdf.rb', line 233 def hypercoll_graphkit() case [:command] when :help "Plot of the effect of hypercollisions" when :options return [] else raise "This only works for spectrogk" unless spectrogk? [:modify_variable] = Proc.new do |varname, narray, dimhash| #dimnames = dimhash.keys p varname, dimhash if varname == "gnew2_ta" shape = narray.shape m = dimhash['m'] mmax = new_netcdf_file.var('hermite').get.to_a.size - 1 p 'shape',shape for ig in 0...shape[0] for it in 0...shape[1] for ik in 0...shape[2] for il in 0...shape[3] for ie in 0...shape[4] for is in 0...shape[5] narray[ig,it,ik,il,ie,is]*=send(:nu_h_ + (is+1).to_sym)*(m[il]/mmax)**send(:nexp_h_ + (is+1).to_sym) end end end end end end end narray end [:graphkit_name] = 'cdf_gnew2_ta' return smart_graphkit() end end |
#hyperviscosity_graphkit(options) ⇒ Object
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 |
# File 'lib/gs2crmod/read_netcdf.rb', line 196 def hyperviscosity_graphkit() case [:command] when :help "Plot of the effect of hyperviscosity" when :options return [] else raise "This only works for spectrogk" unless spectrogk? [:modify_variable] = Proc.new do |varname, narray, dimhash| #dimnames = dimhash.keys shape = narray.shape if varname == "gnew2_ta" #p dimhash #p dimhash['Y'] ky = dimhash['Y'].to_a.to_gslv kx = dimhash['X'].to_a.to_gslv shape = narray.shape for ig in 0...shape[0] for it in 0...shape[1] for ik in 0...shape[2] for il in 0...shape[3] for ie in 0...shape[4] for is in 0...shape[5] narray[ig,it,ik,il,ie,is]*=(ky[ik]**2.0 + kx[it]**2.0)**(2*@nexp)*@d_hypervisc end end end end end end end narray end [:graphkit_name] = 'cdf_gnew2_ta' return smart_graphkit() end end |
#incomplete ⇒ Object
391 392 393 |
# File 'lib/gs2crmod/gs2.rb', line 391 def incomplete return (not 100 == percent_complete) end |
#ingen ⇒ Object
Run the ingen tool on the input file
356 357 358 359 360 361 362 |
# File 'lib/gs2crmod/ingen.rb', line 356 def ingen Dir.chdir(@directory) do ing = File.dirname(File.(@executable)) + '/ingen' success = system "#{ing} #@run_name.in" warning("Could not run ingen... make sure that ingen is in the same folder as @executable and can be run on the login nodes if you want this to work") unless success end end |
#input_file_extension ⇒ Object
1263 1264 1265 |
# File 'lib/gs2crmod/gs2.rb', line 1263 def input_file_extension '.in' end |
#input_file_header ⇒ Object
919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 |
# File 'lib/gs2crmod/gs2.rb', line 919 def input_file_header run_namelist_backwards_compatibility <<EOF !============================================================================== ! GS2 INPUT FILE automatically generated by CodeRunner !============================================================================== ! ! GS2 is a gyrokinetic flux tube initial value turbulence code ! which can be used for fusion or astrophysical plasmas. ! ! See http://gyrokinetics.sourceforge.net ! ! CodeRunner is a framework for the automated running and analysis ! of large simulations. ! ! See http://coderunner.sourceforge.net ! ! Created on #{Time.now.to_s} ! by CodeRunner version #{CodeRunner::CODE_RUNNER_VERSION.to_s} ! !============================================================================== EOF end |
#jump(options) ⇒ Object
1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 |
# File 'lib/gs2crmod/gsl_data.rb', line 1393 def jump() # ep 'kx_shift', kx_shift(options) jump = ((kx_shift() / list(:kx)[2]).round) case [:t_index] when 1 return jump else if @g_exb and @g_exb.abs > 0 return jump + 1 else return 0 end end end |
#kx_indexed ⇒ Object
1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 |
# File 'lib/gs2crmod/gsl_data.rb', line 1452 def kx_indexed return cache[:kx_indexed] if cache[:kx_indexed] #kx = cache[:kx_array] ||= gsl_vector('kx').to_a #kxphys = kx.from_box_order #min_index = kx.min_index + 1 #cache[:kx_indexed] ||= kx.size.times.inject({}) do |hash, kx_element| #hash[kx_element + 1] = kxphs kx = gsl_vector('kx') size = kx.size box = GSL::Vector::Int.indgen(size) + 1 zero_element = kx.abs.min_index phys = box.subvector(zero_element, size-zero_element).connect(box.subvector(0, zero_element)) cache[:kx_indexed] = [phys.to_a, box.to_a].transpose.inject({}){|hash, (phys, box)| hash[phys] = box; hash} end |
#kx_shift(options) ⇒ Object
1386 1387 1388 1389 1390 1391 |
# File 'lib/gs2crmod/gsl_data.rb', line 1386 def kx_shift() # ep options return 0 unless @g_exb and @g_exb.abs > 0.0 #p options return - list(:ky)[[:ky_index]] * list(:t)[([:t_index] or list(:t).keys.max)] * @g_exb end |
#latex_graphs ⇒ Object
This section defines a selection of graphs which are written to a latex file when the CR function write_report is called. To add your own, simply copy one a similar looking graph and modify it to your needs. The requirements to use the latex report writing is further specified in CodeRunner.
1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 |
# File 'lib/gs2crmod/gs2.rb', line 1420 def latex_graphs #will have a different set of graphs to look at depending on whether linear or nonlinear if @nonlinear_mode == "off" #make up a list of graphs that are to be included. The order of the arguments is [code to generate graphkit, LaTeX description and ] graphs = [ #[(kit = phi2_by_mode_vs_time({kx_index:1, ky_index:1}); kit.xlabel=%[TEST]; kit.gp.term = "post eps color enhanced size 3.5in,2.33in"; kit.gp.output = "test.eps"; kit), "This is a test graph written into a \LaTeX file. \n\n \\myfigure{test.eps}"] [(kit = phi2tot_vs_time_graphkit; kit.data[0].title=""; kit.gp.logscale="y"; kit.file_name = "phi2tot.eps"; kit), "Total $\\phi^2$ versus time."], [(kit = growth_rate_vs_ky_graphkit; kit.data[0].title=""; kit.file_name = "growth_rate_vs_ky.eps"; kit), "Growth rate $\\gamma_E$ as a function of $k_y$ averaged over $k_x$ (if applicable)."], if @grid_option=="range" then [(kit = graphkit('efnmag', {norm:true, kx_index:1, ky_index: :all}); kit.data.each{|dk| dk.title=""}; kit.gp.logscale="y"; kit.file_name = "efnmag.eps"; kit.data.shift; kit), "Normalized magnitude of the eigenfunction as a function of $\\theta$ for all $k_y$'s in the simulation."] end, if @grid_option=="single" then [(kit = graphkit('efnmag', {norm:true, kx_index:1, ky_index:1}); kit.data.each{|dk| dk.title=""}; kit.gp.logscale="y"; kit.file_name = "efnmag.eps"; kit), "Normalized magnitude of the eigenfunction as a function of $\\theta$ for all $k_y$'s in the simulation."] end, ].compact else graphs = [ [(kit = ky_spectrum_graphkit; kit.gp.logscale="y"; kit.file_name = "ky_spectrum.eps"; kit), "$k_y$ spectrum at the final time step averaged over $k_x$."], [(kit = kx_spectrum_graphkit; kit.gp.logscale="y"; kit.file_name = "kx_spectrum.eps"; kit), "$k_x$ spectrum at the final time step averaged over $k_y$."], [(kit = spectrum_graphkit(no_zonal:true); kit.gp.view="map"; kit.gp.logscale="z"; kit.file_name = "spectrum.eps"; kit), "2D spectrum versus $k_x$ and $k_y$ without zonal flows."], [(kit = hflux_tot_vs_time_graphkit; kit.file_name = "hflux_tot_vs_time.eps"; kit), "Total heat flux $Q_{tot}$ as a function of time."], [(kit = es_heat_flux_vs_time_graphkit(species_index:1); kit.file_name = "es_heat_1_vs_time.eps"; kit), "Heat flux of species 1 versus time."], if @nspec > 1 then [(kit = es_heat_flux_vs_time_graphkit(species_index:2); kit.file_name = "es_heat_2_vs_time.eps"; kit), "Heat flux of species 2 versus time."] end, [(kit = es_heat_vs_ky_graphkit(species_index:1); kit.gp.logscale="y" ; kit.file_name = "es_heat_1_vs_ky.eps"; kit), "Heat flux of species 1 as a function of $k_y$."], if @nspec > 1 then [(kit = es_heat_vs_ky_graphkit(species_index:2); kit.gp.logscale="y" ; kit.file_name = "es_heat_2_vs_ky.eps"; kit), "Heat flux of species 2 as a function of $k_y$."] end, [(kit = es_heat_vs_ky_vs_kx_graphkit; kit.gp.view="map" ; kit.file_name = "es_heat_vs_ky_vs_kx.eps"; kit), "2D total heat flux spectrum as a function of $k_x$ and $k_y$."], [(kit = phi_real_space_graphkit(n0:1, thetamin:get_list_of(:theta).length/2, thetamax:get_list_of(:theta).length/2, gs2_coordinate_factor:1.0); kit.gp.view="map" ; kit.file_name = "phi_real_space.eps"; kit), "Potential fluctuations at the final time step vs GS2 $x$ and $y$ at the outboard midplane."], [(kit = density_real_space_graphkit(n0:1, species_index:1, thetamin:get_list_of(:theta).length/2, thetamax:get_list_of(:theta).length/2, gs2_coordinate_factor:1.0); kit.gp.view="map" ; kit.file_name = "density_real_space.eps"; kit), "Density fluctuations for species 1 at the final time step vs GS2 $x$ and $y$ at the outboard midplane."], if @nspec > 1 then [(kit = density_real_space_graphkit(n0:1, species_index:2, thetamin:get_list_of(:theta).length/2, thetamax:get_list_of(:theta).length/2, gs2_coordinate_factor:1.0); kit.gp.view="map" ; kit.file_name = "density_real_space.eps"; kit), "Density fluctuations for species 2 at the final time step vs GS2 $x$ and $y$ at the outboard midplane."] end, [(kit = es_mom_flux_vs_time_graphkit(species_index:1); kit.file_name = "es_mom_flux_1_vs_time.eps"; kit), "Momentum flux for species 1 as a function of time."], if @nspec > 1 then [(kit = es_mom_flux_vs_time_graphkit(species_index:2); kit.file_name = "es_mom_flux_2_vs_time.eps"; kit), "Momentum flux for species 2 as a function of time."] end, [(kit = zonal_spectrum_graphkit; kit.gp.logscale="y"; kit.file_name = "zonal_spectrum.eps"; kit), "Zonal spectrum at the final time step."], ].compact end end |
#lenardbern_graphkit(options) ⇒ Object
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 |
# File 'lib/gs2crmod/read_netcdf.rb', line 269 def lenardbern_graphkit() case [:command] when :help "Plot of the effect of Lenard Bernstein collisions" when :options return [] else raise "This only works for spectrogk" unless spectrogk? [:modify_variable] = Proc.new do |varname, narray, dimhash| #dimnames = dimhash.keys if varname == "gnew2_ta" m = dimhash['m'] shape = narray.shape for ig in 0...shape[0] for it in 0...shape[1] for ik in 0...shape[2] for il in 0...shape[3] for ie in 0...shape[4] for is in 0...shape[5] narray[ig,it,ik,il,ie,is]*=send(:nu_ + (is+1).to_sym)*m[il] end end end end end end end narray end [:graphkit_name] = 'cdf_gnew2_ta' kit = smart_graphkit() return kit end end |
#list_of_restart_files ⇒ Object Also known as: lorf
Return a list of restart file paths (relative to the run directory).
578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 |
# File 'lib/gs2crmod/gs2.rb', line 578 def list_of_restart_files Dir.chdir(@directory) do files = Dir.entries.grep(/^\.\d+$/) files = Dir.entries.grep(/\.nc(?:\.\d|_ene)/) if files.size == 0 if files.size == 0 (Dir.entries.find_all{|dir| FileTest.directory? dir} - ['.', '..']).each do |dir| files = Dir.entries(dir).grep(/\.nc(?:\.\d|_ene)/).map{|file| dir + "/" + file} break if files.size == 0 end end #if files.size == 0 # This just finds a .nc file (w/o a number) in the nc folder if using single restart file if files.size == 0 files = Dir.entries('nc').grep(/\.nc/).map{|file| 'nc' + "/" + file} end #if files.size == 0 return files end # Dir.chdir(@directory) do end |
#max_es_heat_amp(species_index) ⇒ Object
749 750 751 |
# File 'lib/gs2crmod/calculations.rb', line 749 def max_es_heat_amp(species_index) @transient_es_heat_flux_amplification_at_species_at_ky[species_index-1].values.max end |
#max_nprocs_no_x ⇒ Object
ep parallelisation
345 346 347 348 |
# File 'lib/gs2crmod/ingen.rb', line 345 def max_nprocs_no_x parallelisation = cumulative_gridpoints parallelisation[parallelisation.keys[parallelisation.keys.index('x') - 1]] end |
#max_trans_phi ⇒ Object
743 744 745 746 747 |
# File 'lib/gs2crmod/calculations.rb', line 743 def max_trans_phi phivec = gsl_vector('phi2tot_over_time') offset = 30 phivec.subvector(20, phivec.size - 20).max end |
#namelist_test_failed(namelist, tst) ⇒ Object
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 |
# File 'lib/gs2crmod/ingen.rb', line 38 def namelist_test_failed(namelist, tst) return <<EOF --------------------------- Test Failed --------------------------- Namelist: #{namelist} Test: #{tst[:test]} Explanation: #{tst[:explanation]} --------------------------- EOF end |
#ncclose ⇒ Object
54 55 56 57 |
# File 'lib/gs2crmod/gsl_data.rb', line 54 def ncclose cache[:netcdf_file].close cache.delete(:netcdf_file) end |
#ncdump(names = nil, values = nil, extension = '.out.nc') ⇒ Object
305 306 307 308 309 |
# File 'lib/gs2crmod/gs2.rb', line 305 def ncdump(names=nil, values=nil, extension = '.out.nc') names = [names] unless !names or names.class == Array names.map!{|name| name.to_s} if names pp NumRu::NetCDF.open(@run_name + extension).vars(names).to_a.sort{|var1, var2| var1.name <=> var2.name}.map{|var| values ? [var.name, var.send(values)] : var.name.to_sym} end |
#netcdf_file ⇒ Object
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 |
# File 'lib/gs2crmod/gsl_data.rb', line 32 def netcdf_file #if @runner.cache[:runs] and (open = @runner.cache[:runs].keys.find_all{|id| @runner.cache[:runs][id][:netcdf_file]}).size > 200 #ep "my id", id if (open = @runner.run_list.keys.find_all{|id| @runner.run_list[id].cache[:netcdf_file]}).size > 200 open = open.sort_by{|id| @runner.run_list[id].cache[:netcdf_file_otime]} @runner.run_list[open[0]].ncclose end if cache[:netcdf_file] and not [:Complete, :Failed].include? @status ncclose end cache[:netcdf_file_otime] = Time.now.to_i cache[:netcdf_file] ||= NumRu::NetCDF.open(netcdf_filename) cache[:netcdf_file].sync cache[:netcdf_file] end |
#netcdf_filename ⇒ Object
49 50 51 |
# File 'lib/gs2crmod/gsl_data.rb', line 49 def netcdf_filename @directory + '/' + @run_name + '.out.nc' end |
#netcdf_smart_reader ⇒ Object
171 172 173 |
# File 'lib/gs2crmod/read_netcdf.rb', line 171 def netcdf_smart_reader NetcdfSmartReader.new(new_netcdf_file) end |
#no_restarts ⇒ Object
Returns true if this run has not been restarted, false if it has. This allows one to get data from the final run of a series of restarts.
642 643 644 645 |
# File 'lib/gs2crmod/gs2.rb', line 642 def no_restarts raise NoRunnerError unless @runner !(@runner.runs.find{|run| run.restart_id == @id}) end |
#old_smart_graphkit(options) ⇒ Object
185 186 187 188 189 190 191 192 193 194 |
# File 'lib/gs2crmod/read_netcdf.rb', line 185 def old_smart_graphkit() case [:command] when :help "An old smart graphkit is a direct plot of a given variable from the old netcdf file. The name of the graphkit is the name of the variable prefixed by 'nc_'. To plot, for example, the heat flux vs time, you would give the graph name nc_hflux_tot. You can use index specifiers in the the options; for example, to plot the potential as a function of kx and ky for a given time index, you would use the graph name nc_phi2_by_mode, and the options {t_index: n}. To plot the potential as function of kx for a given ky and time would use the options {t_index, n, ky_index: m}. For each dimension you can specify the index, or a minium and/or maximum." when :options [:kx_index, :ky_index, :t_index, :e_index, :l_index, :s_index, :kxmax, :kxmin, :kx_element] else return OldNetcdfSmartReader.new(netcdf_file).graphkit([:graphkit_name].sub(/^nc_/, ''), ) end end |
#optimisation_flags ⇒ Object
307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 |
# File 'lib/gs2crmod/ingen.rb', line 307 def optimisation_flags [ :opt_redist_persist, :opt_redist_persist_overlap, :opt_redist_nbk, :opt_redist_init, :intmom_sub, :intspec_sub, #:local_field_solve, :do_smart_update, :field_subgath, :field_option, :field_local_allreduce, :field_local_allreduce_sub, :minnrow, :opt_init_bc, :opt_source ] end |
#parallelizable_meshpoints ⇒ Object
Gives a guess as to the maximum number of meshpoints which can be parallelized (i.e. excluding ntheta)
889 890 891 |
# File 'lib/gs2crmod/gs2.rb', line 889 def parallelizable_meshpoints approximate_grid_size / ntheta end |
#parameter_string ⇒ Object
905 906 907 |
# File 'lib/gs2crmod/gs2.rb', line 905 def parameter_string return "#{@run_name}.in" end |
#parameter_transition(run) ⇒ Object
395 396 |
# File 'lib/gs2crmod/gs2.rb', line 395 def parameter_transition(run) end |
#percent_complete ⇒ Object
431 432 433 |
# File 'lib/gs2crmod/gs2.rb', line 431 def percent_complete @completed_timesteps ? @completed_timesteps.to_f / @nstep.to_f * 100.0 : @percent_of_total_time end |
#physical_kx_index(box_kx_index) ⇒ Object
1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 |
# File 'lib/gs2crmod/gsl_data.rb', line 1472 def physical_kx_index(box_kx_index) return kx_indexed.key(box_kx_index) kx = cache[:kx_gslv] ||= gsl_vector('kx') return kx.from_box_order.to_a.index(kx[box_kx_index-1]) + 1 #kx = cache[:kx_gslv] ||= gsl_vector('kx') #index_of_min_kx = cache[:index_of_min_kx] ||= kx.min_index + 1 # kx.min_index returns a 0-based index #if box_kx_index < index_of_min_kx #box_kx_index + (1 + kx.size - index_of_min_kx) #else #box_kx_index - (index_of_min_kx - 1) #end end |
#plot_efit_file ⇒ Object
1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 |
# File 'lib/gs2crmod/gs2.rb', line 1118 def plot_efit_file Dir.chdir(@directory) do text = File.read(@eqfile) text_lines = text.split("\n") first_line = text_lines[0].split(/\s+/) second_line = text_lines[1].split(/\s+/) nr = first_line[-2].to_i nz = first_line[-1].to_i rwidth = second_line[1].to_f zwidth = second_line[2].to_f rmag = second_line[3].to_f nlines = (nr.to_f/5.0).ceil nlines_psi = ((nr*nz).to_f/5.0).ceil start = 5 f = text_lines[start...(start+=nlines)].join(" ").split(nil).map{|s| s.to_f}.to_gslv pres = text_lines[(start)...(start += nlines)].join(" ").split(nil).map{|s| s.to_f}.to_gslv _ = text_lines[(start)...(start += nlines)].join(" ").split(nil).map{|s| s.to_f}.to_gslv _ffprime = text_lines[(start)...(start+= nlines)].join(" ").split(nil).map{|s| s.to_f}.to_gslv psi = text_lines[(start)...(start += nlines_psi)].join(" ") q = text_lines[(start)...(start += nlines)].join(" ").split(nil).map{|s| s.to_f}.to_gslv nbound = text_lines[start...start+=1].join(" ").to_i rz = text_lines[(start)...(start += nbound*2)].join(" ").split(/\s+/) rz.shift rbound, zbound, _ = rz.inject([[], [], true]){|arr,val| arr[2] ? [arr[0].push(val), arr[1], false] : [arr[0], arr[1].push(val), true]} #rbound.shift psi = psi.split(/\s+/) psi.shift psi.map!{|v| v.to_f} psi_arr = SparseTensor.new(2) k = 0 for i in 0...nz for j in 0...nr psi_arr[j,i] = psi[k] k+=1 end end kit = GraphKit.quick_create([((0...nr).to_a.to_gslv - nr/2 - 1 )/(nr-1)*rwidth+rmag, ((0...nz).to_a.to_gslv-nz/2 + 1)/(nz-1) * zwidth, psi_arr], [rbound, zbound, rbound.map{|r| 0}]) kit.gp.contour = "" kit.gp.view = "map" #kit.gp.nosurface = "" kit.gp.cntrparam = "levels 20" kit.data[0].gp.with = 'l' kit.data[1].gp.with = 'l lw 2 nocontours' kit.gnuplot kit2 = GraphKit.quick_create([pres/pres.max],[f/f.max],[q/q.max]) kit2.data[0].title = 'Pressure/Max Pressure' kit2.data[1].title = 'Poloidal current function/Max poloidal current function' kit2.data[2].title = 'Safety factor/Max Safety factor' kit2.gnuplot #p ['f', f, 'p', pres, 'ffprime', ffprime, 'nlines', nlines, 'psi', psi, 'q', q, 'nbound', nbound, 'rbound', rbound, 'zbound', zbound] end end |
#print_out_line ⇒ Object
435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 |
# File 'lib/gs2crmod/gs2.rb', line 435 def print_out_line logf(:print_out_line) name = @run_name name += " (res: #@restart_id)" if @restart_id name += " real_id: #@real_id" if @real_id beginning = sprintf("%2d:%d %-60s %1s:%2.1f(%s) %3s%1s %1s", @id, @job_no, name, @status.to_s[0,1], @run_time.to_f / 60.0, @nprocs.to_s, percent_complete, "%", @converged.to_s) if @ky beginning += sprintf("%3s %4s %4s", @ky, @growth_rates[@ky], @real_frequencies[@ky]) elsif @nonlinear_mode == "off" beginning += sprintf("%3s %4s %4s", @fastest_growing_mode, @max_growth_rate, @freq_of_max_growth_rate) elsif @nonlinear_mode == "on" # p @hflux_tot_stav beginning += " sat:#{saturated.to_s[0]}" beginning += sprintf(" hflux:%1.2e", @hflux_tot_stav) if @hflux_tot_stav beginning += sprintf("+/-%1.2e", @hflux_tot_stav_error) if @hflux_tot_stav_error beginning += sprintf(" momflux:%1.2e", @es_mom_flux_stav.values.sum) if @es_mom_flux_stav and @es_mom_flux_stav.values[0] beginning += ' SC:' + @spectrum_check.map{|c| c.to_s}.join(',') if @spectrum_check beginning += ' VC:' + @vspace_check.map{|c| sprintf("%d", ((c*10.0).to_i rescue -1))}.join(',') if @vspace_check end beginning += " ---#{@comment}" if @comment beginning end |
#process_directory_code_specific ⇒ Object
This method, as its name suggests, is called whenever CodeRunner is asked to analyse a run directory.this happens if the run status is not :Complete, or if the user has specified recalc_all(-A on the command line) or reprocess_all (-a on the command line).
the structure of this function is very simple: first it calls get_status to determine the directory status, i.e. :Complete, :Incomplete, :NotStarted or :Failed, then it gets the time, which is the GS2 time at the end of the run, and it also gets the run_time, which is the wall clock time of the run. Finally,if non-linear mode is switched off, it calls calculate_growth_rates_and_frequencies, and if the non-linear mode is switched on, it calls calculate_time_averaged_fluxes.
231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 |
# File 'lib/gs2crmod/gs2.rb', line 231 def process_directory_code_specific run_namelist_backwards_compatibility unless @status == :Queueing get_status end eputs "Run #@status: #@run_name" if [:Complete,:Failed].include? @status try_to_get_error_file @sys = @@successful_trial_system return if @status == :NotStarted or @status == :Failed or @status == :Queueing begin percent_complete = get_completed_timesteps/@nstep @percent_of_total_time = percent_complete rescue get_time @percent_of_total_time = @time / (@delt*@nstep) * 100.0 rescue 0.0 end return if @status == :Incomplete get_run_time calculate_results end |
#recheck ⇒ Object
class ListSubmitter
788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 |
# File 'lib/gs2crmod/gs2.rb', line 788 def recheck logf(:recheck) Dir.chdir(@directory) do logi('@runner.object_id', @runner.object_id) log('@runner.class', @runner.class) #runner = @runner instance_variables.each{|var| instance_variable_set(var, nil) unless var == :@runner} begin File.delete(".code_runner_run_data") rescue Errno::ENOENT end begin File.delete("code_runner_results.rb") rescue Errno::ENOENT end logi(:@checked_converged, @checked_converged) logi('@runner.object_id after reset', @runner.object_id) log('@runner.class', @runner.class) process_directory end end |
#renew_info_file ⇒ Object
1035 1036 1037 |
# File 'lib/gs2crmod/gs2.rb', line 1035 def renew_info_file Dir.chdir(@directory){make_info_file("#@run_name.in")} end |
#restart(new_run) ⇒ Object
529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 |
# File 'lib/gs2crmod/gs2.rb', line 529 def restart(new_run) #new_run = self.dup (rcp.variables).each{|v| new_run.set(v, send(v)) if send(v)} @naming_pars.delete(:preamble) SUBMIT_OPTIONS.each{|v| new_run.set(v, self.send(v)) unless new_run.send(v)} #(rcp.results + rcp.gs2_run_info).each{|result| new_run.set(result, nil)} new_run.is_a_restart = true new_run.ginit_option = "many" new_run.delt_option = "check_restart" #if Dir.entries(@directory).include? "nc" #old_restart_run_name = (@restart_run_name or Dir.entries(@directory + '/nc').grep(/\.nc/)[0].sub(/\.nc\.\d+$/, '')) #new_run.restart_file = File.expand_path("#@directory/nc/#{old_restart_run_name}.nc") #else #new_run.restart_file = File.expand_path("#@directory/#@run_name.nc") #end new_run.restart_id = @id new_run.restart_run_name = @run_name @runner.nprocs = @nprocs if @runner.nprocs == "1" # 1 is the default so this means the user probably didn't specify nprocs raise "Restart must be on the same number of processors as the previous run: new is #{new_run.nprocs.inspect} and old is #{@nprocs.inspect}" if !new_run.nprocs or new_run.nprocs != @nprocs # @runner.parameters.each{|var, value| new_run.set(var,value)} if @runner.parameters # ep @runner.parameters new_run.run_name = nil new_run.naming_pars = @naming_pars new_run.update_submission_parameters(new_run.parameter_hash_string, false) if new_run.parameter_hash new_run.naming_pars.delete(:restart_id) new_run.generate_run_name eputs 'Copying Restart files', '' FileUtils.makedirs(new_run.directory + '/nc') #old_dir = File.dirname(@restart_file) new_run.restart_file = "#@run_name.nc" #+ File.basename(@restart_file) #.sub(/\.nc/, '') new_run.restart_dir = "nc" #files = Dir.entries(old_dir).grep(/\.nc(?:\.\d|_ene)/) #files = Dir.entries(old_dir).grep(/^\.\d+$/) if files.size == 0 files = list_of_restart_files.map do |file| @directory + "/" + file end files.each_with_index do |file , index| eputs "\033[2A" # Terminal jargon - go back one line eputs "#{index+1} out of #{files.size}" num = file.scan(/(?:\.\d+|_ene)$/)[0] #FileUtils.cp("#{old_dir}/#{file}", "nc/#@restart_file#{num}") FileUtils.cp(file, new_run.directory + "/nc/#{new_run.restart_file}#{num}") end #@runner.submit(new_run) new_run end |
#restart_chain ⇒ Object
648 649 650 651 652 653 654 655 656 657 658 659 660 |
# File 'lib/gs2crmod/gs2.rb', line 648 def restart_chain if @restart_id return @runner.run_list[@restart_id].restart_chain end chain = [] currid = @id loop do chain.push currid break unless (restrt = @runner.runs.find{|run| run.restart_id == currid}) currid = restrt.id end return chain end |
#run_heuristic_analysis ⇒ Object
This method overrides a method defined in heuristic_run_methods.rb in the CodeRunner source. It is called when CodeRunner cannot find any of its own files in the folder being analysed. It takes a GS2 input file and generates a CodeRunner info file. This means that GS2 runs which were not run using CodeRunner can nonetheless be analysed by it. In order for it to be called the -H flag must be specified on the command line.
1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 |
# File 'lib/gs2crmod/gs2.rb', line 1041 def run_heuristic_analysis ep 'run_heuristic_analysis', Dir.pwd infiles = Dir.entries.grep(/^[^\.].*\.in$/) ep infiles raise CRMild.new('No input file') unless infiles[0] raise CRError.new("More than one input file in this directory: \n\t#{infiles}") if infiles.size > 1 input_file = infiles[0] ep 'asdf' @nprocs ||= "1" @executable ||= "/dev/null" make_info_file(input_file, false) end |
#run_namelist_backwards_compatibility ⇒ Object
1098 1099 1100 1101 1102 1103 |
# File 'lib/gs2crmod/gs2.rb', line 1098 def run_namelist_backwards_compatibility SPECIES_DEPENDENT_VARIABLES.each do |var| set(var + "_1".to_sym, (send(var + "_1".to_sym) or send(var + "_i".to_sym) or send(var))) set(var + "_2".to_sym, (send(var + "_2".to_sym) or send(var + "_e".to_sym))) end end |
#run_namelist_tests(namelist, hash, enum = nil) ⇒ Object
Checks input parameters for inconsistencies and prints a report.
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 |
# File 'lib/gs2crmod/ingen.rb', line 56 def run_namelist_tests(namelist, hash, enum = nil) ext = enum ? "_#{enum}" : "" hash[:must_pass].each do |tst| error(namelist_test_failed(namelist, tst)) unless instance_eval(tst[:test]) end if hash[:must_pass] hash[:should_pass].each do |tst| warning(namelist_test_failed(namelist, tst)) unless instance_eval(tst[:test]) end if hash[:should_pass] hash[:variables].each do |var, var_hash| #gs2_var = (var_hash[:gs2_name] or var) cr_var = var+ext.to_sym value = send(cr_var) if value.kind_of? Array value.each{|v| test_variable(namelist, var, var_hash, ext, v)} else test_variable(namelist, var, var_hash, ext, value) end end end |
#saturated_time_average(name, options) ⇒ Object
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 |
# File 'lib/gs2crmod/calculations.rb', line 46 def saturated_time_average(name, ) # calculate_saturation_time_index unless @saturation_time_index # p 'sat', @saturation_time_index, 'max', list(:t).keys.max raise "saturation_time_index not calculated for #@run_name" unless @saturation_time_index [:t_index_window] = [@saturation_time_index, list(:t).keys.max - 1] #ep gsl_vector(name, {}).size #ep name, options begin vec = gsl_vector(name, ) rescue GSL::ERROR::EINVAL # IF the vector doesn't have enough values for each timestep (due to run aborting early?), this error will be thrown. [:t_index_window] = [@saturation_time_index, gsl_vector(name, {}).size] retry rescue NoMethodError eputs "Warning: could not calculate #{name} saturated time average" return nil end tvec = gsl_vector('t', ) dt = tvec.subvector(1, tvec.size - 1) - tvec.subvector(0, tvec.size - 1) trapezium = (vec.subvector(1, tvec.size - 1) + vec.subvector(0, tvec.size - 1)) / 2.0 return trapezium.mul(dt).sum / dt.sum end |
#saturated_time_average_error(name, options) ⇒ Object
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 |
# File 'lib/gs2crmod/calculations.rb', line 73 def saturated_time_average_error(name, ) # calculate_saturation_time_index unless @saturation_time_index [:t_index_window] = [@saturation_time_index, list(:t).keys.max] begin vec = gsl_vector(name, ) tavg = GSL::Vector.alloc(vec.size) vec.size.times.each{|i| tavg[i] = vec.subvector(i+1).mean} rescue NoMethodError eputs "Warning: could not calculate #{name} saturated_time_average_error" return nil end # tavg = 0.0; i = 0 # tavg_vec = vec.collect{|val| tavg += val; tavg = tavg / (i+=1); tavg} # ind = GSL::Vector.indgen(vec.size) # i = 0 # begin # fit = GSL::Fit::linear(ind.subvector(i, ind.size - i) , vec.subvector(i, ind.size - i)) # # p fit[1].abs - 100.0 * fit[4].abs # i += 1 # (eputs "Not Saturated"; break) if i > vec.size * 0.9 # end while (fit[1].abs - Math.sqrt(fit[4].abs)) > 0 # p fit # fit_vec = ind * fit[1] + fit[0] # # p tavg.size # # GraphKit.autocreate({x: {data: gsl_vector(name, {})}}) # (GraphKit.autocreate({x: {data: tavg}}) + GraphKit.autocreate({x: {data: vec}}) + GraphKit.autocreate({x: {data: fit_vec}})).gnuplot return tavg.sd end |
#saturated_time_average_std_dev(name, options) ⇒ Object
103 104 105 106 107 108 109 110 111 112 113 |
# File 'lib/gs2crmod/calculations.rb', line 103 def saturated_time_average_std_dev(name, ) # calculate_saturation_time_index unless @saturation_time_index [:t_index_window] = [@saturation_time_index, list(:t).keys.max] begin vec = gsl_vector(name, ) rescue NoMethodError eputs "Warning: could not calculate #{name} saturated_time_average_std_dev" return nil end return vec.sd end |
#sc(min) ⇒ Object
795 796 797 |
# File 'lib/gs2crmod/calculations.rb', line 795 def sc(min) return @spectrum_check.min >= min end |
#set_nprocs ⇒ Object
850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 |
# File 'lib/gs2crmod/gs2.rb', line 850 def set_nprocs if (nprocs_in = @nprocs) =~ /^x/ max = max_nprocs_no_x nodes = 0 @nprocs = "#{nodes}#{nprocs_in}" loop do nodes += 1 @nprocs = "#{nodes}#{nprocs_in}" if actual_number_of_processors > max nodes -= 1 @nprocs = "#{nodes}#{nprocs_in}" break end end end end |
#show_graph ⇒ Object
515 516 517 518 519 520 521 522 |
# File 'lib/gs2crmod/gs2.rb', line 515 def show_graph thegraph = special_graph('phi2tot_vs_time_all_kys') thegraph.title += " for g_exb = #{@g_exb.to_f.to_s}" thegraph.show sleep 1.5 # @decaying = Feedback.get_boolean("Is the graph decaying?") thegraph.kill end |
#smart_graphkit(options) ⇒ Object
175 176 177 178 179 180 181 182 183 184 |
# File 'lib/gs2crmod/read_netcdf.rb', line 175 def smart_graphkit() case [:command] when :help "A smart graphkit is a direct plot of a given variable from the new netcdf file. The name of the graphkit is the name of the variable prefixed by 'cdf_'. To plot, for example, the heat flux vs time, you would give the graph name cdf_heat_flux_tot. You can use index specifiers in the the options; for example, to plot the potential as a function of kx and ky for a given time index, you would use the graph name cdf_phi2_by_mode, and the options {t_index: n}. To plot the potential as function of kx for a given ky and time would use the options {t_index, n, Y_index: m}. For each dimension you can specify the index, or a minium and/or maximum." when :options [:X_index, :Y_index, :t_index, :e_index, :l_index, :s_index, :Xmax, :Xmin, :X_element] else netcdf_smart_reader.graphkit([:graphkit_name].sub(/^cdf_/, ''), ) end end |
#spec_chec(min, *dirns) ⇒ Object
781 782 783 784 785 786 787 788 789 790 791 792 793 |
# File 'lib/gs2crmod/calculations.rb', line 781 def spec_chec(min, *dirns) return @spectrum_check.zip([0, 1, 2]).inject(true) do |bool, (check,dirn)| unless dirns.include? dirn bool and true else unless check >= min false else bool and true end end end end |
#species_letter ⇒ Object
625 626 627 |
# File 'lib/gs2crmod/gs2.rb', line 625 def species_letter species_type(1).downcase[0,1] end |
#species_type(index) ⇒ Object
629 630 631 632 633 634 635 636 637 |
# File 'lib/gs2crmod/gs2.rb', line 629 def species_type(index) if rcp.variables.include? :type_1 type = send(:type_ + index.to_sym) else types = rcp.variables.find_all{|var| var.to_s =~ /^type/}.map{|var| send(var)} type = types[index.to_i - 1] end type end |
#spectrogk? ⇒ Boolean
63 64 65 |
# File 'lib/gs2crmod/gs2.rb', line 63 def spectrogk? false end |
#standardize_restart_files ⇒ Object
Put restart files in the conventional location, i.e. nc/run_name.proc
600 601 602 603 604 605 606 607 608 609 |
# File 'lib/gs2crmod/gs2.rb', line 600 def standardize_restart_files Dir.chdir(@directory) do FileUtils.makedirs('nc') list_of_restart_files.each do |file| proc_id = file.scan(/\.\d+$|_ene$/)[0] #p 'proc_id', proc_id FileUtils.mv(file, "nc/#@run_name.nc#{proc_id}") end end end |
#stop ⇒ Object
1106 1107 1108 |
# File 'lib/gs2crmod/gs2.rb', line 1106 def stop `touch #@directory/#@run_name.stop` end |
#test_failed(namelist, var, gs2_var, tst) ⇒ Object
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
# File 'lib/gs2crmod/ingen.rb', line 18 def test_failed(namelist, var, gs2_var, tst) return <<EOF --------------------------- Test Failed --------------------------- Namelist: #{namelist} Variable: #{var} GS2 Name: #{gs2_var} Value: #{send(var)} Test: #{tst[:test]} Explanation: #{tst[:explanation]} --------------------------- EOF end |
#test_variable(namelist, var, var_hash, ext, value) ⇒ Object
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 |
# File 'lib/gs2crmod/ingen.rb', line 76 def test_variable(namelist, var, var_hash, ext, value) gs2_var = (var_hash[:gs2_name] or var) cr_var = var+ext.to_sym if value and (not var_hash[:should_include] or eval(var_hash[:should_include])) var_hash[:must_pass].each do |tst| error(test_failed(namelist, cr_var, gs2_var, tst)) unless value.instance_eval(tst[:test]) end if var_hash[:must_pass] var_hash[:should_pass].each do |tst| warning(test_failed(namelist, cr_var, gs2_var, tst)) unless value.instance_eval(tst[:test]) end if var_hash[:should_pass] if (var_hash[:allowed_values] or var_hash[:text_options]) tst = {test: "#{(var_hash[:allowed_values] or var_hash[:text_options]).inspect}.include? self", explanation: "The variable must have one of these values"} error(test_failed(namelist, cr_var, gs2_var, tst)) unless value.instance_eval(tst[:test]) end end end |
#update_physics_parameters_from_miller_input_file(file) ⇒ Object
1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 |
# File 'lib/gs2crmod/gs2.rb', line 1005 def update_physics_parameters_from_miller_input_file(file) hash = self.class.parse_input_file(file) hash[:parameters].each do |var, val| set(var,val) end hash[:theta_grid_parameters].each do |var, val| next if [:ntheta, :nperiod].include? var set(var, val) end hash[:dist_fn_knobs].each do |var, val| next unless [:g_exb].include? var set(var, val) end hash[:theta_grid_eik_knobs].each do |var, val| next unless [:s_hat_input, :beta_prime_input].include? var set(var, val) end hash[:species_parameters_2].each do |var, val| #next unless [:s_hat_input, :beta_prime_input].include? var set((var.to_s + '_2').to_sym, val) end hash[:species_parameters_1].each do |var, val| #next unless [:s_hat_input, :beta_prime_input].include? var set((var.to_s + '_1').to_sym, val) end end |
#vim_output ⇒ Object Also known as: vo
1110 1111 1112 |
# File 'lib/gs2crmod/gs2.rb', line 1110 def vim_output system "vim -Ro #{output_file} #{error_file} #@directory/#@run_name.error #@directory/#@run_name.out " end |
#vim_stdout ⇒ Object Also known as: vo1
1114 1115 1116 |
# File 'lib/gs2crmod/gs2.rb', line 1114 def vim_stdout system "vim -Ro #{output_file} " end |
#visually_check_growth_rate(ky = nil) ⇒ Object
499 500 501 502 503 504 505 506 507 508 509 510 511 512 |
# File 'lib/gs2crmod/gs2.rb', line 499 def visually_check_growth_rate(ky=nil) logf :visually_check_growth_rate phi_vec = gsl_vector(:phi2_by_ky_over_time, {ky: ky}) t_vec = gsl_vector(:t) constant, growth_rate = GSL::Fit::linear(t_vec, 0.5*GSL::Sf::log(phi_vec)).slice(0..1) eputs growth_rate graph = @@phi2tot_vs_time_template.graph(["#{constant} * exp (2 * #{growth_rate} * x)"], [[[t_vec, phi_vec], "u 1:2 title 'phi2tot #{@run_name}' w p"]], {"set_show_commands" => "\nset log y\n", "point_size"=>'1.0'}) # eputs graph.inline_data.inspect graph.show gets graph.kill end |
#warning(message) ⇒ Object
7 8 9 |
# File 'lib/gs2crmod/ingen.rb', line 7 def warning() eputs "Warning: " + ; sleep 0.3 end |
#write_input_file ⇒ Object
846 847 848 |
# File 'lib/gs2crmod/gs2.rb', line 846 def write_input_file File.open(@run_name + ".in", 'w'){|file| file.puts input_file_text} end |