Output error when trying to write out averages (FD)

[rpvin]

I am getting index access issues as shown below, when I try to write out averages (Shell, Global and Azimuthal). I have traced the error to vals(self%phi_local(p),r,t), where self%phi_local has zeros and the indices of vals start at 1. Not sure why this is happening now, after I updated to the latest version of Rayleigh.

forrtl: severe (408): fort: (3): Subscript #1 of the array VALS has value 0 which is less than the lower bound of 1 Image PC Routine Line Source
libifcoremt.so.5 00001555553E0F42 for_emit_diagnost Unknown Unknown rayleigh.dbg 00000000004A4482 parallel_io_mp_ca 896 Parallel_IO.F90 rayleigh.dbg 00000000007B1B29 spherical_io_mp_s 730 Spherical_IO.F90 rayleigh.dbg 00000000007AA89F spherical_io_mp_a 582 Spherical_IO.F90 rayleigh.dbg 00000000009A8868 diagnostics_veloc 43 Diagnostics_Velocity_Field.F90 rayleigh.dbg 00000000010651DE diagnostics_inter 207 Diagnostics_Interface.F90 rayleigh.dbg 00000000011BF1B5 sphere_physical_s 164 Sphere_Physical_Space.F90 rayleigh.dbg 000000000123F4AA sphere_drivermp 181 Sphere_Driver.F90 rayleigh.dbg 00000000012CBA70 MAIN 66 Main.F90 rayleigh.dbg 0000000000402DC2 Unknown Unknown Unknown libc-2.28.so 00001555470557E5 libc_start_main Unknown Unknown rayleigh.dbg 0000000000402CCE Unknown Unknown Unknown

Also, to get to this error, I had to comment out the following in PDE_coefficients.F90, which was giving a different error, as it was a new addition that wasn't required for the anelastic FD implementation.

!Integrate dsdr to obtain a self-consistent entropy profile. Set s=0 at the upper surface. ref%entropy(1) = 0.0d0 geofac = (radius(1)3 - radius(N_R)3)/3.0d0 Do i = 2, N_R ref%entropy(i) = ref%entropy(i-1) - geofacref%dsdr(i)radial_integral_weights(i)/Radius(i)**2 Enddo ref%exp_entropy = exp(ref%entropy/pressure_specific_heat) ref%dsdr_over_cp = ref%dsdr/pressure_specific_heat Call log_deriv(ref%dsdr_over_cp(:), ref%d2s_over_cp(:), no_log=.true.)

[feathern]

Thanks. I understand the source of the second error. I can fix that later this week. As far as the indexing error goes, can you upload or copy/paste your main input file?

[rpvin]

&problemsize_namelist n_r = 2400 n_theta = 32 nprow = -1 npcol = -1 radial_grid_file = 'radius_2400.dat' / &numerical_controls_namelist chebyshev = .false. bandsolve= .false. sparsesolve= .false. deriv_cluge = .true. / &physical_controls_namelist rotation = .true. magnetism = .false. advect_reference_state = .true. viscous_heating = .false. ohmic_heating = .false. / &temporal_controls_namelist max_time_step = 1d3 min_time_step = 0.1 max_time_minutes = 600000 checkpoint_interval = 10000 quicksave_minutes = 30 num_quicksaves = 3 cflmin = 0.4d0 cflmax = 0.6d0 / &io_controls_namelist !stdout_file = "mylog" !stdout_flush_interval = 1 / &output_namelist meridional_values = 1,2,3,501 meridional_frequency = 10 meridional_nrec = 10 meridional_indices = 0, 64, 128, 192, 256, 320, 384, 448

equatorial_values = 1,2,3,501 equatorial_frequency = 10 equatorial_nrec = 10

!shellavg_values = 2,3,501,502,401 !shellavg_frequency = 10 !shellavg_nrec = 10

!globalavg_values = 401, 1101, 501 !globalavg_frequency = 10 !globalavg_nrec = 10

!shellslice_levels_nrm = 0.1, 0.2, 0.5, 0.9 !shellslice_values = 1,2,3,501,503 !shellslice_frequency = 10 !shellslice_nrec = 20

shellspectra_values = 1, 2, 3, 401, 501 shellspectra_frequency = 10 shellspectra_nrec = 10 shellspectra_levels = 2380,2300,2200,2100,2000,1970,1960,1950,1940,1930,1920,1900,1850,1800,1750,1300,1200,1100,200,100,50

full3d_values = 1,2,3,501 full3d_frequency = 25

!azavg_values = 1,2,3,501,401,1807,1808 !azavg_frequency = 10 !azavg_nrec = 10 /

&Boundary_Conditions_Namelist no_slip_boundaries = .false. strict_L_Conservation = .false. !dtdr_bottom = 0.0d0 T_Top = 0.0d0 T_Bottom = 0.0d0 fix_tvar_top = .true. fix_tvar_bottom = .true. fix_dtdr_bottom = .false. !impose_dipole_field = .true. !Br_bottom = 1000000.0d0 !C10_bottom = 1e25 !C11_bottom = 0.0d0 !C1m1_bottom = 0.0d0 / &Initial_Conditions_Namelist init_type = 7 ! change this to -1 to restart, 7 is random magnetic_init_type = -1 !mag_amp = 1.0d0 temp_amp = 1.0d1 temp_w = 0.01d4 restart_iter = 10000 ! restarts from most recent checkpoint / &Test_Namelist / &Reference_Namelist reference_type = 4 custom_reference_file='cref_from_MESA_eta_2_nu_10_kappa_10.dat' heating_type = 1 pressure_specific_heat = 3.5d8 angular_velocity = 1.8d-5

!override_constants=T !override_constant(5)= T !ra_constants(5)= 8d12 ! nu_top

!override_constant(6)= T !ra_constants(6)= 8d12 ! kappa_top

!override_constant(1)= T !ra_constants(1) = 0.5d-5 ! 2 x angular velocity

!override_constant(10)= T !ra_constants(10) = 1e6 ! luminosity boosting factor / &Transport_Namelist nu_type = 3 kappa_type = 3 eta_type = 3 /

[feathern]

Thanks. Also, I see nprow and npcol are set to -1. What configuration did you use when the indexing error occurred? That should be enough for me to have a look.

[rpvin]

Sorry, I don't understand what you mean by configuration here? Do you mean the reference state profile?

[rpvin]

config.zip

[feathern]

Nprow and npcol aren't actually -1. If you aren't setting them at the command line, the Rayleigh is doing it automatically. If that s the case, just tell me n cpu.

[rpvin]

Yes, ncpu = 8 here

[rpvin]

Right, so I just ran the code with ntheta = 256 and this problem didn't come up. We changed it to 32 to make the tests run faster but seems like there was some indexing issue with the parallel IO, when nr = 2400, ntheta = 32 and ncpu = 8.

[tukss]

Can you try either ntheta=64 or ncpu=4? Maybe the domain is just too small in this case.

[rpvin]

Sorry for the delayed reply, but the error I mentioned in the first message returned for nr = 2400, ntheta = 256 and ncpu = 600. The crash happened within 5 s after the code starts running

[feathern]

OK, so there are two issues here.

1.) The log_deriv routine, whose call you had to comment out, only refers to the Chebyshev derivative routines. I didn't get to this as planned last week, but I will work on this soon.

2.) The crash you're asking about it coming about because of how you've specfied your Meridional_Slices outputs. The indices that you are allowed to output can run from 1 up through 2*n_theta. Outputting index 0 will always cause a problem, and outputting indices as high as 448 will cause problems when n_theta is only 32.

It's probably worth having Rayleigh do a sanity check and complain, but I generally suggest using meridional_indices_nrm instead of meridional_indices. Those values can run from 0 to 1, and it makes your main_input file more resilient to resolution changes. A value of 0 corresponds to phi index=1 and a value of one corresponds to phi index=ntheta*2. Similarly so for the other 'nrm' indices that are used for the radial points.

I'll think on how to have Rayleigh deal with this after I modify the log_deriv routine.

Cheers, Nick

[feathern]

I'm going to close this issue since this has been fixed now.