diff --git a/demo/app/train-cloud-microphysics.F90 b/demo/app/train-cloud-microphysics.F90 index eea7565c4..fbdff9f73 100644 --- a/demo/app/train-cloud-microphysics.F90 +++ b/demo/app/train-cloud-microphysics.F90 @@ -18,7 +18,7 @@ program train_cloud_microphysics !! Internal dependencies: use phase_space_bin_m, only : phase_space_bin_t use NetCDF_file_m, only: NetCDF_file_t - use NetCDF_variable_m, only: NetCDF_variable_t + use NetCDF_variable_m, only: NetCDF_variable_t, tensors implicit none character(len=*), parameter :: usage = new_line('a') // new_line('a') // & @@ -31,7 +31,7 @@ program train_cloud_microphysics 'The presence of a file named "stop" halts execution gracefully.' type command_line_arguments_t - integer num_epochs, start_step, stride, num_bins, report_interval + integer num_epochs, start_step, stride, num_bins, report_step integer, allocatable :: end_step character(len=:), allocatable :: base_name real cost_tolerance @@ -105,7 +105,7 @@ function get_command_line_arguments() result(command_line_arguments) base_name, epochs_string, start_string, end_string, stride_string, bins_string, report_string, tolerance_string real cost_tolerance integer, allocatable :: end_step - integer num_epochs, num_bins, start_step, stride, report_interval + integer num_epochs, num_bins, start_step, stride, report_step base_name = command_line%flag_value("--base") ! gfortran 13 seg faults if this is an association epochs_string = command_line%flag_value("--epochs") @@ -122,54 +122,24 @@ function get_command_line_arguments() result(command_line_arguments) read(epochs_string,*) num_epochs - if (len(stride_string)==0) then - stride = 1 - else - read(stride_string,*) stride - end if - - if (len(start_string)==0) then - start_step = 1 - else - read(start_string,*) start_step - end if - - if (len(report_string)==0) then - report_interval = 1 - else - read(report_string,*) report_interval - end if - - if (len(bins_string)/=0) then - read(bins_string,*) num_bins - else - num_bins = 1 - end if + stride = default_integer_or_read(1, stride_string) + start_step = default_integer_or_read(1, start_string) + report_step = default_integer_or_read(1, report_string) + num_bins = default_integer_or_read(1, bins_string) + cost_tolerance = default_real_or_read(5E-8, tolerance_string) if (len(end_string)/=0) then allocate(end_step) read(end_string,*) end_step end if - if (len(start_string)==0) then - start_step = 1 - else - read(start_string,*) start_step - end if - - if (len(tolerance_string)==0) then - cost_tolerance = 5.0E-08 - else - read(tolerance_string,*) cost_tolerance - end if - if (allocated(end_step)) then command_line_arguments = command_line_arguments_t( & - num_epochs, start_step, stride, num_bins, report_interval, end_step, base_name, cost_tolerance & + num_epochs, start_step, stride, num_bins, report_step, end_step, base_name, cost_tolerance & ) else command_line_arguments = command_line_arguments_t( & - num_epochs, start_step, stride, num_bins, report_interval, null(), base_name, cost_tolerance & + num_epochs, start_step, stride, num_bins, report_step, null(), base_name, cost_tolerance & ) end if @@ -194,15 +164,17 @@ subroutine read_train_write(training_configuration, args, plot_file) enumerator :: dpotential_temperature_t=1, dqv_dt, dqc_dt, dqr_dt, dqs_dt end enum - associate(input_names => & - [string_t("pressure"), string_t("potential_temperature"), string_t("temperature"), & - string_t("qv"), string_t("qc"), string_t("qr"), string_t("qs")] & - ) + !associate(input_names => & + ! [string_t("pressure"), string_t("potential_temperature"), string_t("temperature"), & + ! string_t("qv"), string_t("qc"), string_t("qr"), string_t("qs")] & + !) + associate(input_names => [string_t("qv"), string_t("qc")]) + allocate(input_variable(size(input_names))) associate(input_file_name => args%base_name // "_input.nc") - print *,"Reading network inputs from " // input_file_name + print *,"Reading physics-based model inputs from " // input_file_name associate(input_file => netCDF_file_t(input_file_name)) @@ -222,13 +194,14 @@ subroutine read_train_write(training_configuration, args, plot_file) end associate end associate - associate(output_names => [string_t("potential_temperature"),string_t("qv"), string_t("qc"), string_t("qr"), string_t("qs")]) + !associate(output_names => [string_t("potential_temperature"),string_t("qv"), string_t("qc"), string_t("qr"), string_t("qs")]) + associate(output_names => [string_t("qv"), string_t("qc")]) allocate(output_variable(size(output_names))) associate(output_file_name => args%base_name // "_output.nc") - print *,"Reading network outputs from " // output_file_name + print *,"Reading physics-based model outputs from " // output_file_name associate(output_file => netCDF_file_t(output_file_name)) @@ -252,7 +225,7 @@ subroutine read_train_write(training_configuration, args, plot_file) block type(NetCDF_variable_t) derivative(size(output_variable)) - print *,"Calculating time derivatives" + print *,"Calculating desired neural-network model outputs" associate(dt => NetCDF_variable_t(output_time - input_time, "dt")) do v = 1, size(derivative) @@ -275,45 +248,33 @@ subroutine read_train_write(training_configuration, args, plot_file) type(input_output_pair_t), allocatable :: input_output_pairs(:) type(tensor_t), allocatable, dimension(:) :: inputs, outputs real, allocatable :: cost(:) - integer i, lon, lat, level, time, network_unit, io_status, epoch, end_step + integer i, network_unit, io_status, epoch, end_step integer(int64) start_training, finish_training associate( network_file => args%base_name // "_network.json") open(newunit=network_unit, file=network_file, form='formatted', status='old', iostat=io_status, action='read') - ! if (allocated(args%end_step)) then - ! end_step = args%end_step - ! else - ! end_step = t_end - ! end if - - ! print *,"Defining tensors from time step", args%start_step, "through", end_step, "with strides of", args%stride + if (allocated(args%end_step)) then + end_step = args%end_step + else + end_step = input_variable(1)%end_step() + end if - ! ! The following temporary copies are required by gfortran bug 100650 and possibly 49324 - ! ! See https://gcc.gnu.org/bugzilla/show_bug.cgi?id=100650 and https://gcc.gnu.org/bugzilla/show_bug.cgi?id=49324 - ! inputs = [( [( [( [( & - ! tensor_t( & - ! [ pressure_in(lon,lat,level,time), potential_temperature_in(lon,lat,level,time), temperature_in(lon,lat,level,time), & - ! qv_in(lon,lat,level,time), qc_in(lon,lat,level,time), qr_in(lon,lat,level,time), qs_in(lon,lat,level,time) & - ! ] & - ! ), lon = 1, size(qv_in,1))], lat = 1, size(qv_in,2))], level = 1, size(qv_in,3))], & - ! time = args%start_step, end_step, args%stride)] - - ! outputs = [( [( [( [( & - ! tensor_t( & - ! [dpt_dt(lon,lat,level,time), dqv_dt(lon,lat,level,time), dqc_dt(lon,lat,level,time), dqr_dt(lon,lat,level,time), & - ! dqs_dt(lon,lat,level,time) & - ! ] & - ! ), lon = 1, size(qv_in,1))], lat = 1, size(qv_in,2))], level = 1, size(qv_in,3))], & - ! time = args%start_step, end_step, args%stride)] + print *,"Defining input tensors starting from time step", args%start_step, "through", end_step, "with strides of", args%stride + inputs = tensors(input_variable, step_start = args%start_step, step_end = end_step, step_stride = args%stride) + + print *,"Defining output tensors starting from time step", args%start_step, "through", end_step, "with strides of", args%stride + outputs = tensors(output_variable, step_start = args%start_step, step_end = end_step, step_stride = args%stride) - ! print *,"Calculating output tensor component ranges." - ! output_extrema: & - ! associate( & - ! output_minima => [minval(dpt_dt), minval(dqv_dt), minval(dqc_dt), minval(dqr_dt), minval(dqs_dt)], & - ! output_maxima => [maxval(dpt_dt), maxval(dqv_dt), maxval(dqc_dt), maxval(dqr_dt), maxval(dqs_dt)] & - ! ) + print *,"Calculating output tensor component ranges." + tensor_extrema: & + associate( & + input_minima => [( input_variable(v)%minimum(), v=1,size( input_variable) )] & + ,input_maxima => [( input_variable(v)%maximum(), v=1,size( input_variable) )] & + ,output_minima => [( output_variable(v)%minimum(), v=1,size(output_variable) )] & + ,output_maxima => [( output_variable(v)%maximum(), v=1,size(output_variable) )] & + ) ! output_map: & ! associate( output_map => tensor_map_t(layer = "outputs", minima = output_minima, maxima = output_maxima)) ! read_or_initialize_network: & @@ -376,7 +337,7 @@ subroutine read_train_write(training_configuration, args, plot_file) ! " in ", count(occupied)," out of ", size(occupied, kind=int64), " bins." ! end block ! end associate output_map - ! end associate output_extrema + end associate tensor_extrema ! print *,"Normalizing the remaining input and output tensors" ! input_output_pairs = trainable_network%map_to_training_ranges(input_output_pairs) @@ -423,7 +384,7 @@ subroutine read_train_write(training_configuration, args, plot_file) ! associate(converged => average_cost <= args%cost_tolerance) ! image_1_maybe_writes: & - ! if (me==1 .and. any([converged, epoch==[first_epoch,last_epoch], mod(epoch,args%report_interval)==0])) then + ! if (me==1 .and. any([converged, epoch==[first_epoch,last_epoch], mod(epoch,args%report_step)==0])) then ! print *, epoch, average_cost ! write(plot_file%plot_unit,*) epoch, average_cost @@ -471,4 +432,30 @@ subroutine read_train_write(training_configuration, args, plot_file) end subroutine read_train_write + pure function default_integer_or_read(default, string) result(set_value) + integer, intent(in) :: default + character(len=*), intent(in) :: string + integer set_value + + if (len(string)==0) then + set_value = default + else + read(string,*) set_value + end if + + end function + + pure function default_real_or_read(default, string) result(set_value) + real, intent(in) :: default + character(len=*), intent(in) :: string + real set_value + + if (len(string)==0) then + set_value = default + else + read(string,*) set_value + end if + + end function + end program train_cloud_microphysics diff --git a/demo/src/NetCDF_variable_m.f90 b/demo/src/NetCDF_variable_m.f90 index c55b446b0..113145236 100644 --- a/demo/src/NetCDF_variable_m.f90 +++ b/demo/src/NetCDF_variable_m.f90 @@ -4,10 +4,12 @@ module NetCDF_variable_m use NetCDF_file_m, only : NetCDF_file_t use kind_parameters_m, only : default_real, double_precision use julienne_m, only : string_t + use fiats_m, only : tensor_t implicit none private public :: NetCDF_variable_t + public :: tensors type NetCDF_variable_t(k) integer, kind :: k = default_real @@ -19,16 +21,22 @@ module NetCDF_variable_m procedure, private, non_overridable :: default_real_input, double_precision_input, default_real_input_character_name, double_precision_input_character_name generic :: conformable_with => default_real_conformable_with, double_precision_conformable_with procedure, private, non_overridable :: default_real_conformable_with, double_precision_conformable_with - generic :: rank => default_real_rank, double_precision_rank - procedure, private, non_overridable :: default_real_rank, double_precision_rank - generic :: any_nan => default_real_any_nan, double_precision_any_nan - procedure, private, non_overridable :: default_real_any_nan, double_precision_any_nan - generic :: operator(-) => default_real_subtract, double_precision_subtract - procedure, private, non_overridable :: default_real_subtract, double_precision_subtract - generic :: operator(/) => default_real_divide, double_precision_divide - procedure, private, non_overridable :: default_real_divide, double_precision_divide - generic :: assignment(=) => default_real_assign, double_precision_assign - procedure, private, non_overridable :: default_real_assign, double_precision_assign + generic :: rank => default_real_rank , double_precision_rank + procedure, private, non_overridable :: default_real_rank , double_precision_rank + generic :: end_step => default_real_end_step , double_precision_end_step + procedure, private, non_overridable :: default_real_end_step , double_precision_end_step + generic :: any_nan => default_real_any_nan , double_precision_any_nan + procedure, private, non_overridable :: default_real_any_nan , double_precision_any_nan + generic :: minimum => default_real_minimum , double_precision_minimum + procedure, private, non_overridable :: default_real_minimum , double_precision_minimum + generic :: maximum => default_real_maximum , double_precision_maximum + procedure, private, non_overridable :: default_real_maximum , double_precision_maximum + generic :: operator(-) => default_real_subtract , double_precision_subtract + procedure, private, non_overridable :: default_real_subtract , double_precision_subtract + generic :: operator(/) => default_real_divide , double_precision_divide + procedure, private, non_overridable :: default_real_divide , double_precision_divide + generic :: assignment(=) => default_real_assign , double_precision_assign + procedure, private, non_overridable :: default_real_assign , double_precision_assign end type interface NetCDF_variable_t @@ -121,6 +129,18 @@ elemental module function double_precision_rank(self) result(my_rank) integer my_rank end function + elemental module function default_real_end_step(self) result(end_step) + implicit none + class(NetCDF_variable_t), intent(inout) :: self + integer end_step + end function + + elemental module function double_precision_end_step(self) result(end_step) + implicit none + class(NetCDF_variable_t(double_precision)), intent(inout) :: self + integer end_step + end function + elemental module function default_real_subtract(lhs, rhs) result(difference) implicit none class(NetCDF_variable_t), intent(in) :: lhs, rhs @@ -169,6 +189,49 @@ elemental module function double_precision_any_nan(self) result(any_nan) logical any_nan end function + elemental module function default_real_minimum(self) result(minimum) + implicit none + class(NetCDF_variable_t), intent(in) :: self + real minimum + end function + + elemental module function double_precision_minimum(self) result(minimum) + implicit none + class(NetCDF_variable_t(double_precision)), intent(in) :: self + real minimum + end function + + elemental module function default_real_maximum(self) result(maximum) + implicit none + class(NetCDF_variable_t), intent(in) :: self + real maximum + end function + + elemental module function double_precision_maximum(self) result(maximum) + implicit none + class(NetCDF_variable_t(double_precision)), intent(in) :: self + real maximum + end function + + module function tensors(NetCDF_variables, step_start, step_end, step_stride) + implicit none + type(NetCDF_variable_t), intent(in) :: NetCDF_variables(:) + type(tensor_t), allocatable :: tensors(:) + integer, optional :: step_start, step_end, step_stride + end function + + elemental module function default_real_end_time(self) result(end_time) + implicit none + class(NetCDF_variable_t), intent(inout) :: self + integer end_time + end function + + elemental module function double_precision_end_time(self) result(end_time) + implicit none + class(NetCDF_variable_t), intent(inout) :: self + integer end_time + end function + end interface end module \ No newline at end of file diff --git a/demo/src/NetCDF_variable_s.f90 b/demo/src/NetCDF_variable_s.f90 index cdd4e9fc3..0e7dec711 100644 --- a/demo/src/NetCDF_variable_s.f90 +++ b/demo/src/NetCDF_variable_s.f90 @@ -21,6 +21,11 @@ module procedure double_precision_upper_bounds end interface + interface default_if_not_present + module procedure default_integer_if_not_present + module procedure default_real_if_not_present + end interface + contains module procedure default_real_copy @@ -149,6 +154,36 @@ pure function double_precision_components_allocated(NetCDF_variable) result(allo end associate end procedure + module procedure default_real_end_step + select case(self%rank()) + case(1) + end_step = ubound(self%values_1D_,1) + case(2) + end_step = ubound(self%values_2D_,2) + case(3) + end_step = ubound(self%values_3D_,3) + case(4) + end_step = ubound(self%values_4D_,4) + case default + error stop "NetCDF_variable_s(default_real_end_step): unsupported rank" + end select + end procedure + + module procedure double_precision_end_step + select case(self%rank()) + case(1) + end_step = ubound(self%values_1D_,1) + case(2) + end_step = ubound(self%values_2D_,2) + case(3) + end_step = ubound(self%values_3D_,3) + case(4) + end_step = ubound(self%values_4D_,4) + case default + error stop "NetCDF_variable_s(double_precision_end_step): unsupported rank" + end select + end procedure + pure function default_real_lower_bounds(NetCDF_variable) result(lbounds) type(NetCDF_variable_t), intent(in) :: NetCDF_variable integer, allocatable :: lbounds(:) @@ -377,4 +412,147 @@ pure function double_precision_upper_bounds(NetCDF_variable) result(ubounds) end select end procedure + module procedure tensors + + integer t_start, t_end, t_stride + + select case(NetCDF_variables(1)%rank()) + case(4) + + t_start = default_if_not_present(1, step_start ) + t_stride = default_if_not_present(1, step_stride) + t_end = default_if_not_present(size(NetCDF_variables(1)%values_4D_,4), step_end) + + associate( longitudes => size(NetCDF_variables(1)%values_4D_,1) & + ,latitudes => size(NetCDF_variables(1)%values_4D_,2) & + ,levels => size(NetCDF_variables(1)%values_4D_,3) & + ) + block + integer v, lon, lat, lev, time + + tensors = [( [( [( [( tensor_t( [( NetCDF_variables(v)%values_4D_(lon,lat,lev,time), v=1,size(NetCDF_variables) )] ), & + lon = 1, longitudes)], lat = 1, latitudes)], lev = 1, levels)], time = t_start, t_end, t_stride)] + end block + end associate + + case default + error stop "NetCDF_variable_s(tensors): unsupported rank)" + end select + + end procedure + + module procedure default_real_end_time + select case(self%rank()) + case (1) + end_time = size(self%values_1D_,1) + case (2) + end_time = size(self%values_2D_,2) + case (3) + end_time = size(self%values_3D_,3) + case (4) + end_time = size(self%values_4D_,4) + case default + error stop 'NetCDF_variable_s(default_real_end_time): unsupported rank' + end select + end procedure + + module procedure double_precision_end_time + select case(self%rank()) + case (1) + end_time = size(self%values_1D_,1) + case (2) + end_time = size(self%values_2D_,2) + case (3) + end_time = size(self%values_3D_,3) + case (4) + end_time = size(self%values_4D_,4) + case default + error stop 'NetCDF_variable_s(double_precision_end_time): unsupported rank' + end select + end procedure + + module procedure default_real_minimum + select case(self%rank()) + case (1) + minimum = minval(self%values_1D_) + case (2) + minimum = minval(self%values_2D_) + case (3) + minimum = minval(self%values_3D_) + case (4) + minimum = minval(self%values_4D_) + case default + error stop 'NetCDF_variable_s(default_real_minimum): unsupported rank' + end select + end procedure + + module procedure double_precision_minimum + select case(self%rank()) + case (1) + minimum = minval(self%values_1D_) + case (2) + minimum = minval(self%values_2D_) + case (3) + minimum = minval(self%values_3D_) + case (4) + minimum = minval(self%values_4D_) + case default + error stop 'NetCDF_variable_s(double_precision_minimum): unsupported rank' + end select + end procedure + + module procedure default_real_maximum + select case(self%rank()) + case (1) + maximum = maxval(self%values_1D_) + case (2) + maximum = maxval(self%values_2D_) + case (3) + maximum = maxval(self%values_3D_) + case (4) + maximum = maxval(self%values_4D_) + case default + error stop 'NetCDF_variable_s(default_real_maximum): unsupported rank' + end select + end procedure + + module procedure double_precision_maximum + select case(self%rank()) + case (1) + maximum = maxval(self%values_1D_) + case (2) + maximum = maxval(self%values_2D_) + case (3) + maximum = maxval(self%values_3D_) + case (4) + maximum = maxval(self%values_4D_) + case default + error stop 'NetCDF_variable_s(double_precision_maximum): unsupported rank' + end select + end procedure + + pure function default_integer_if_not_present(default_value, optional_argument) result(set_value) + integer, intent(in) :: default_value + integer, intent(in), optional :: optional_argument + integer set_value + + if (present(optional_argument)) then + set_value = optional_argument + else + set_value = default_value + end if + end function + + pure function default_real_if_not_present(default_value, optional_argument) result(set_value) + real, intent(in) :: default_value + real, intent(in), optional :: optional_argument + real set_value + + if (present(optional_argument)) then + set_value = optional_argument + else + set_value = default_value + end if + end function + end submodule NetCDF_variable_s \ No newline at end of file